Abstarct: Magnetotelluric (MT) is an electromagnetic (EM) method that uses natural EM fields for mapping the variations of electrical resistivity within the earth. High penetration depth of EM fields in MT method has made its wide application in deep target exploration such as geothermal and hydrocarbon resources. For this reason MT method has special status among the geophysical methods. Similar to other geophysical methods, modeling of the MT data is necessary for interpreting of subsurface geological structures. The MT inverse modeling is known as a non-linear and highly ill-posed problem, so to prevent the dilemma of none unique solution and to obtain meaningful results, it is normally solved using Tikhonov’s regularization method. In solving this regularized inverse problem, a regularization parameter is required to control the weight given to minimization of the model norm relative to minimization of the data misfit. Clearly, a large amount of regularization parameter favors a small solution model norm at the cost of a large data misfit, while a small amount of regularization parameter has the opposite effect. In the case of small regularization value the non-real structures may be produced in the model. Thus, the regularization parameter is an important quantity which controls the properties of the regularized solution, hence it must be chosen carefully. In the available MT data inversion algorithms, this parameter is estimated using the trial and error or by the discrepancy principle method. In this study, it is attempted to develop one or several techniques to choose the optimal value of the regularization parameter automatically for two and three-dimensional (2-D & 3-D) inversion of the MT data so that the speed, accuracy and the computational efficiency of inversion are enhanced. To achieve the goals, a comprehensive study was done to find out how this regularization parameter is selected in the general ill-posed problems and 12 suitable schemes have been selected. Once the necessary improvement has been done on some of them, they were programmed within the Matlab software to be used as a part of the main 2-D inversion program of MT data. These methods are including discrepancy principle, Generalized Cross Validation (GCV), modified GCV, L-curve with maximum curvature algorithm, L-curve with adaptive pruning algorithm, error estimation from Brezineski et al. and Reichel et al., fixed point, quasi-optimality, Chi-square, a scheme baxsed on Meju approach and ratio of misfit to model norm methods. These methods have been coded in the Matlab software and then inserted in an available magnetotelluric (MT2DInvMatlab) 2-D inversion source code. Subsequently inversion of two synthetic and two real data sets have been implemented by the modified 2-D inversion algorithm and the necessary analysis have been done. The results indicate that, the modified GCV, L-curve with adaptive pruning algorithm, ratio of misfit to model norm and a scheme baxsed on Meju approach operate better than to the other methods for selecting the regularization parameter, considering to accuracy, automaticity, and computational efficiency of 2-D inversion. In the next step, ratio of misfit to model norm method was modified to the ratio of misfit to sum of misfit and model norm and coded in the Matlab as a proposed method for choosing regularization parameter automatically in an available source code (WSINV3DMT) for efficient 3-D inversion of MT data. Once the required modification of the main 3-D inversion program was made by implementation of the proposed method for optimally selecting of regularization parameter, then the necessary investigations have been done by performing 3-D inversions on two synthetic and two real MT data sets of Sabalan geothermal field. The obtained results of the 3-D inversion using the proposed method were compared for the accuracy and computational efficiency to the inversion results of the original 3-D inversion program in which the discrepancy principle method is used to select the regularization parameter of the inversion. The results show that the proposed method considerably increases the accuracy of the modeling and in addition reduces computational time of 3-D inversion at least 30 percent compared to the original available method of inversion that both of them and particularly the latter one is very important in 3-D inversion. Using the results of the 3-D inversion obtained by the proposed method and integrating them with other exploration data, the locations of the geothermal sources (i.e. the hot zones) were determined in the north and west part of mount Sabalan. Moreover the locations of relevant alteration zones were identified by low resistivity anomalies wherein they are spread out along the faults and main structures in the area.