Since Yee (1966)^[1] applied the finite difference approach to solve the Maxwell's equations in isotropic media, the numerical solution of the magnetotelluric (MT) modeling has been improved significantly ^[2-16], It’s one of the most usual MT forward techniques. Many classical two-dimensional (2-D) inversion codes are based on finite difference forward subroutines, such as NLCG^[17], Rebocc ^[18], and RRI ^[19]. Currently, some effective three-dimensional (3-D) inversion algorithms such as the conjugate gradient method^[20], the nonlinear conjugate gradient method ^[21], and the data space method ^[22] are also integrated with the finite difference forward subroutines.

To apply finite difference methods to MT numerical modeling, one can solve the electric and magnetic fields simultaneously using

(1)

(2)

where μ is the air magnetic permeability, ω is the angular frequency, σ is the conductivity (the inverse of resistivity, ρ), E is the electric field, and H is the magnetic field, the time dependence for E and H is e^-iωt.

In complex domain (2-D or 3-D), to solve (1) and (2) directly requires forming a large system matrix which is time and memory consuming. In most cases, the Maxwell's equations of second order are used to solve E_p and H_p (called primary fields);

(3)

(4)

Then the auxiliary fields H_a and E_a are solved by (1) and (2), respectively.

There are several factors that lead to errors of the finite difference numerical solution in MT modeling ^[23]. Siripunvaraporn et al (2002) ^[8] indicated that the solutions obtained from formulas in terms of the electric fields are less sensitive to grid resolution than those obtained from the magnetic formulas. deGroot-Hedlin (2006) ^[24] emphasized the importance of methods used to get the auxiliary field at the surface. As 3-D MT modeling becomes popular, more and more attentions are paid to the methods used to solve linear equations; a good preconditioner together with a stable iterative solver can make the MT modeling more feasible. Han et al (2009) ^[25] presents a detailed comparison of 3-D MT modeling algorithms. How to find a better approximation to Maxwell's equations and a faster linear solver is an important challenge to current MT modeling studies ^[26].

In this paper, we make comparisons of MT finite difference numerical solutions with the analytical solutions focused on the grid resolution, ways to form system matrix and boundary conditions, aiming at what are the major factors affecting the accuracy of MT finite difference numerical solutions. We also compare the efficiency of different preconditioners combining with the linear solvers in 2-D MT cases.

The errors caused by the grid resolution is nonphysical but is inherent numerical problem in finite difference methods ^[24]. A simple 1-D model can be used to show how grid resolution affects MT numerical solutions. In a 1-D case, the MT fields to be solved can be at the top of each layer of the homogenous layered model, denoted as H₀, E₀, H₁, E₁, …, H_n, E_n in Fig. 1a. We call this mode the normal grid method. When the fields to be solved are in an alternating order as H₀, E_1/2, …, H_n-1, E_n-1/2, H_n or E₀, H_1/2, …, E_n-1, H_n-1/2, E_n in Fig. 1a, It is called the staggered grid method ^{[1, 5, 16, 27]}. If the results of 1-D forward modeling are used to prepare boundary values for the electric polarization mode of a 2-D or 3-D model, the layers have to be extended upward to the air as the 2-D or 3-D model does.

图 1

Fig. 1

Fig. 1 Relative errors of 1-D magnetotelluric fields solved by different methods of half space models in different grid resolutions. (a) A homogenous layered model with different discrete modes; (b) A 1 Ωm half space model divided with different grids respectively (Md-1 has four layers in total, the thickness of first layer is 1.6 m close to the skin depth at frequency 105 Hz; Md-2 is formed by splitting the first layer of Md-1 into two equal layers, so the first layer of Md-2 is 0.8 m; Md-3 is formed by splitting each layer of Md-1 into tow equal layers.); (c) Relative errors of electric and magnetic fields of Md-1. RE _ | Ex |; relative error of | Ex | at earth surface under Ex-Hy mode. RE_ | Hy |; relative error of | Hy | at earth surface under Ex-Hy mode. RE_ | Hy |; relative error of | Hy | at earth surface under Ex-Hy mode. RE_ | Ex |; relative error of | Ex | at the top of second layer under Ex-Hy mode. The diamonds show solutions with normal grid method, the upward triangles show solutions with staggered grid method; (d) Relative errors of electric and magnetic fields ofMd-2; (e) Relative errors of electric and magnetic fields of Md-3

deGroot-Hedlin et al. (1990) ^[28] pointed out that the spacing of the first layer should be approximately one-third of the skin depth (δ≈500 , here ρ is the resistivity (Ωm), and f is the frequency (Hz)) and a logarithmic depth scale is suggested. Aprea et al. (1997)^[3] used the stricter grid that the thickness of first layer should be one-tenth of the skin depth in order to model the complex topography. In the following, we take the example of a 1 Ωm haff space model subdivided by different grids, and solve the fields for two modes as E_x-H_y and H_x-E_y. The former mode is used to construct both lateral boundary values for 2-D TE (transverse electric) mode and include the fields' values upward through the air; the latter mode is used to construct both lateral boundary values for the 2-D TM (transverse magnetic) mode.

In Fig. 1b, the homogeneous half space model of 1 Ωm is divided with different grids as Md-1, Md-2 and Md-3. The numerical solutions are obtained with the normal grid method and staggered grid method at a frequency range from 10⁵~10^-5Hz. Then we compare the MT fields by the numerical modeling with the analytical solutions. Fig. 1 (c, d, e) show the relative errors of the MT fields solved by different methods for the three half space models in Fig. 1b. The fields E_x, H_y, and E_y are at the earth's surface, that is, the top of the first layer, while the filed H_x is at the top of the second layer because the top boundary in the H_x-E_y mode is at the earth's surface and H_x at the top of the first layer is constant. In the H_x-E_y mode, the relative error of the primary field E_x (as RE_ | H_x |) of Md-2 has little improvement compared with that of Md-1, the RE_ | H_x | of Md-3 has obvious improvement compared with that of Md-2. The relative errors of the auxiliary field H_y (as RE_ | E_y |) of Md-2 and Md-3 are similar and reduced greatly compared with that of Md-1. In the E_x-E_y mode, the primary field is E_x, the RE _ | E_x | is improved steadily from Md-1 to Md-2 and Md-3, while the RE_ | E_y | has little improvement from Md-1 to Md-2, but is improved much from Md-2 to Md-3. From the images of RE-| E_y |, it seems that the staggered grid method can produce a more accurate auxiliary field in the E_x-E_y mode. In sum, by only decreasing the thickness of the first layer of the model, we can only increase the accuracy of the auxiliary field in the H_x-E_y mode and primary field in the E_x-E_y mode; but if we decrease the thickness of the first layer and keep the difference between successive rows in a reasonable range (2 times around in Md-3) simultaneously, the accuracy of both the primary field and auxiliary fields will be increased.

When we use the five point difference equation to solve the second order Maxwell's equation in a 2-D case, we can define the fields in normal grids or staggered grids, the former can be divided into the normal-node method (Fig. 2a) and normal-center method (Fig. 2b). In the latter case, the fields are defined in an alternating order as Fig. 2c, we call it the staggered grid method.

图 2

Fig. 2

Fig. 2 Comparison of 2-D apparent resistivity curves of two different models with normal grid and staggered grid methods, (a) The primary and auxiliary fields are defined at the same nodes of the retangular blocks, called as normal-node method; (b) The primary and auxiliary fields are defined at the center of the blocks' top edges, called as norma-center method; (c) The primary fields are defined at the center of blocks, while the auxiliary fields are defined along the blocks' edges, called as staggered grids. (d) The apparent resistivity curves of Md-4. Solid line shows analytic solution (obtained according to Refs. [11, 12]), diamonds show solutions with normal-center method, crosses show solutions with normal-node method, the upward triangles show solutions with staggered grids. (e) The apparent resistivity curves of Md-5.

In a 2-D case, it requires some special treat-ment of the model/ is resistivity at grid nodes to obtain accurate MT solutions. From the differential form of the second order Maxwells equation in the 2-D TE mode as (5), we need to know the conductivity (σ, the inverse of resistivity ρ) at the points of the defined primary fields. The diferential form of the TM mode as (6) is a little complex; we need to know the derivatives of resistivity along both directions (y and z in Fig. 2) of the model, and the resistivity at the points of the defined primary fields.

(5)

(6)

The problem arises when we use the normal-center method solving equation (6). It can produce accurate solutions if the models are homogenous or homogenous in layers, but the numerical solutions deviate much from the accurate solutions when the resistivity is discontinuous laterally. We have a comparison of the numerical solutions in two cases on two 2-D models, which are Md-4 (a 2-D model with a 0. 2 km wide low resistivity abnormal zone (10 Ωm) in a 1000 Ωm haff space) and Md-5 (a 2D model with a 0.2 km wide high resistivity abnormal zone (1000 Ωm) in a 10 Ωm half space) as shown in Fig. 2 (d, e), both of which are discretized with the 47X36 mesh. In the TE mode, the apparent resistivity curves from the three methods coincide well with the analytic solutions for both Md-4 and Md-5. In the TM mode, the apparent resistivity curves of the normal-node method and staggeredgrid method coincide well with the analytic solutions, while the responses of the normalcenter method above the abnormal zone deviate more from the analytic solutions. The significantly different resistivity deviation can be found between response above the low resistivity abnormal zone in the model Md-4 at 10³ Hz and response above the high resistivity abnormal zone in the model Md-5 at 10^-2Hz.

The main difficulty in the normal-center method to solve equation (6) is how to estimate the resistivity at the solution points and the difference of model’s resistivity along the horizontal direction. Although we interpolate the resistivity along y(ρ_y, horizontal resistivity) and z(ρ_z, vertical resistivity) respectively for anisotropic model at solution points, t is hard to obtain the accurate difference values of the resistivity as the solution points along the horizontal direction spanning the resistivity boundaries. The accurate responses in the TM mode with the normal-center method can be finally obtained if we use the difference form of equations (1) and (2) at the beginning and form the second order difference equation;

(7)

where ρ is the model’ is resistivity at the point defined by ts subscripts, dy is the grid spacing along the horizontal direction, and dz is the grid spacing along the vertical direction. In equation (7), the resistivity at nodes are formed as;

(8)

Fig. 3 shows the apparent resistivity curves obtained with equation (7) in the TM mode in comparison with the analytic solutions. It indicates that the curves above the high resistivity abnormity coincide well with the analytic solutions.

图 3

Fig. 3

Fig. 3 Comparison of 2-D TM apparent resistivity curves of normal-center method with the difference equation as (7) and the analytic method. (a) The apparent resistivity curves of Md-4 at frequency 103 Hz. Solid line shows analytic solution, diamonds show solutions with normal-center method; (b) The apparent resistivity curves of Md-4 at frequency 10 -2 Hz. (c) The apparent resistivity curves of Md-5 at frequency 103 Hz. (d) The apparent resistivity curves of Md-5 at frequency 10 -2 Hz

We construct another complex 2-D model as Md-6 m Fig. 4 to have a further test of the normal-center method, the apparent resistivity and phase curves in both TE and TM modes are obtained, and we also calculate the responses with the forward subroutine of NLCG^[17]. From the results of three frequencies (as f=10⁴ Hz, 1 Hz, 10^-4 Hz) in Fig. 4, the responses of the normal-center method coincide well with the responses of the forward subroutine of NLCG except minor deviation of phase (φ) above the low resistivity abnormities at very high frequency as 10⁴ Hz.

图 4

Fig. 4

Fig. 4 Comparison of 2-D apparent resistivity and phase curves of normal-center method with the forward results of NLCG [17]. The mesh size of Md-6 is 46 × 45, only core parts of shallow level are shown here. Solid lines show the responses of Md-6 with the forward subroutine of NLCG, empty circles show solutions with normal-center method

For the MT problem, if the Dirichlet boundary condition at the top of the model and a homogenous half space at the bottom are used, the boundary values at the lateral sides are needed to be solved in 2-D or 3-D cases. To see how much the ateral boundary values affect the solutions, we make a simple test with a 1000 Ωm half space model. We obtain the electric fields with the 1-D numerical and analytic methods at the earth surface from frequency 10⁵~10^-5 Hz firstly, and then calculate the 2-D numerical electric fields of the TE mode with the boundary values provided by 1-D numerical solution and analytic solution, respectively. Fig. 5 shows the comparison of the primary fields (including real and imaginary parts) at the earth' s surface of both 1-D and 2-D results. Though the imaginary part of the boundary 1-D values change to some extent, the primary fields of the 2-D TE mode coincide well. It means that there is limitation to increase the accuracy of the 2-D MT finite difference numerical solution only by increasing the accuracy of the boundary values. In MT finite difference numerical modeling, the main difficulties and errors are not due to the approximations made in the boundary conditions^[6].

图 5

Fig. 5

Fig. 5 Comparison of real and imaginary electric fields of both 1-D and 2-D solutions for a 1000 Ωm half space. Solid lines denote 1-D analytic solutions at surface; Diamonds show 1-D numerical solutions at surface; Upward triangles show 2-D numerical solutions under TE mode at surface using 1-D analytic solution as the boundary values at lateral sides; Crosses show 2-D numerical solutions under TE mode at surface using 1-D numerical solution as the boundary values at lateral sides

If Maxwells equations are approx imated with the finite difference method in MT modeling, a linear equations Ax=b are constructed, here x is the vector made up of N unknown primary fields, and A is the system matrix with the leading dimension N. In 2-D or 3-D cases, the number N becomes very large, and the solutions of the equations are usually obtained with the preconditioned herative methods in order to avoid solving the reverse of A directly. Commonly, a preconditioner is a matrix M approximating the system matrix A to form anew system M^-1 Ax=M^-1b, and the new matrix M^-1 A maybe more favorable to convergence.

The preconditioned Krylov subspace methods^[29] are mostly used in solving MT problems. A good preconditioner is an art when it combines well with the iterative solver^[30], which is very essential in 3D MT finite difference forward modeling. For example, Mackie et al. (1994)^[6] used incomplete Cholesky decomposition as preconditioner for Minimum Residual (MR) algorithm and Biconju-gate Gradient (BiCG). Smith (1996)^[10] found that the incomplete Cholesky preconditioner can accelerate convergence of BICG iterative method. Siripunvaraporn et al. (2002)^[8] combined incomplete LU decomposition with the Quasi-Minimal Residual (QMR) method. We make a comparison of the validity of four iterative methods BiCG, QMR, Conjugate Gradient Squared (CGS), and Biconjugate Gradient Stabilized (BiCGSTAB) using a 2-D model as same as Md-4 with a 47 X 36 mesh, the desired residual is 10 ^-12. A Jacobi precondi-tioner (M is a diagonal matrix, its elements are the diagonal of matrix A.) is used for the four linear solvers, the convergence of both Jacobi-BICGSTAB and Jacobi-CGS are not improved at low frequencies, while Jacobi-BICG and Jacobi-QMR can reach the desired residual. Fig. 6 shows the convergence and consuming time of Jacobi-BICG and Jacobi-QMR in the TE and TM modes at frequency 1 Hz. Though the Jacobi-QMR shows a more smooth convergence than Jacobi-BICG, the former consumes more time per iteration than the latter.

图 6

Fig. 6

Fig. 6 Comparison of the convergence (left) and time (right) of Jacobi-BICG and Jacobi-QMR on a 47×36 2-D MT mesh at frequency 1 Hz. Additional 15 air layers for TE mode

The incomplete LU factorization is widely used in constructing preconditioner for the nonsymmetric system matrix, e. g. the DILU is one of the simplest incomplete LU factorization^[29]. By using the DILU preconditioner, we make a further comparison of the solvers BICG, CGS, BICGSTAB and QMR on the same 2-D model Md-4 in the TE and TM modes at frequency 1 Hz. As shown in Fig. 7, all the four solvers can reach the desired residual, and the DILU-BICGSTAB takes the least iterations to reach the desired residual. Taking the example of BICG solver, the iterative number of DILU BICG decreases notably to less than one-tenth of that of Jacobi-BICG, and the total consuming time is about one-fifth of that of Jacobi-BICG. Of course, the inconsistent decrease of iterative number and the corresponding consu-ming time implies that it takes more time to construct and apply DILU preconditioner than Jacobi preconditioner.

图 7

Fig. 7

Fig. 7 Comparison of the convergence (left) and time (right) of BICG, CGS, BICGSTAB, and QMR with a DILU preconditioner on a 47×36 2-D MT mesh at frequency 1 Hz. Additional 15 air layers for TE mode

The finite difference MT numerical modeling is quickly developed, but enhancement of the accuracy of solutions is still an important task. It is necessary to study the ways to find the resistivity at nodes and their difference values, and the ways to form the system equations with an appropriate mode, for example, to divide the model with normal grids or staggered grids. Unlike the usual method to solve E first in the 2-D TE mode, Weaver et al. (1996)^[13] calculated E at surface from the H solutions obtained on both types of grid. If the higher accuracy of primary fields is reached, the accuracy of auxiliary fields at the surface can also be obtained with a suitable method as that of Smith (1996)^[9]. In our example of solving 2-D TM fields with the normal-center method, the difference form of the first order Maxwell’s equations (1) and (2) is used to form the secondorder difference equation (7), in which themodels’ resistivity variance and the resistivity at nodes are carefully considered and accuracy of solutions is improved.

The convergence of four solvers BICG, CGS, BICGSTAB and QMR with a Jacobi preconditioner and DILU preconditioner is discussed respectively, and it is found that the linear solvers (Jacobi-CGS and Jacobi-BICGSTAB) become stagnant at low frequencies in the 2-D case. All the four solvers BICG, CGS, BICGSTAB and QMR with DILU preconditioner can reach the desired residual in the 2-D case at frequency 1 Hz, and the efficiency of DILU preconditioner is super over Jacobi precondi-toner. Of course, DILU is one of the simplest incomplete factorization methods, and some more suitable and effective preconditioners should be tested and adopted in future research.

This study has been supported by the Director Foundation of Institute of Geology, China Earthquake Administration (Grant No. DF-IGCEA-0608-2-16) and the National Natural Science Foundation of China (Grant No. 40534023 and Grant No. 40974042). Mr. Rodi and Mackie are acknowledged for providing their MT inversion codes NLCG. Mr. Weaver, Le Quang and Fischer are acknowledged for providing their codes to calculate 2-D analytic solutions. We also thank the anonymous reviewers and editors for their comments that lead to improvements in this paper.