Finite Difference Algorithm on Non-Uniform Meshes for Modeling 2 D Magnetotelluric Responses

A finite-difference approach with non-uniform meshes was presented for simulating magnetotelluric responses in 2D structures. We presented the calculation formula of this scheme from the boundary value problem of electric field and magnetic field, and compared finite-difference solutions with finite-element numerical results and analytical solutions of a 1D model. First, a homogeneous half-space model was tested and the finite-difference approach can provide very good accuracy for 2D magnetotelluric modeling. Then we compared them to the analytical solutions for the two-layered geo-electric model; the relative errors of the apparent resistivity and the impedance phase were both increased when the frequency was increased. To conclude, we compare our finite-difference simulation results with COMMEMI 2D-0 model with the finite-element solutions. Both results are in close agreement to each other. These comparisons can confirm the validity and reliability of our finite-difference algorithm. Moreover, a future project will extend the 2D structures to 3D, where non-uniform meshes should perform especially well.


Introduction
The magnetotelluric method is a passive electromagnetic exploration technique that measures orthogonal components of the electric and magnetic fields on the Earth's surface [1].The field source is naturally generated by variations in Earth's magnetic field, which provide a wide and continuous spectrum of electromagnetic field waves.These fields induce currents into the Earth, which are measured at the surface and contain information about subsurface resistivity structures.With rapid advances in numerical modeling and ill-posed regularized inversion, the magnetotelluric method has become one of the most important tools for geophysical exploring [2,3].
The magnetotelluric numerical modeling aims to solve the boundary value problem derived from frequency-domain Maxwell's equations and calculate the spatial distribution of electric and magnetic fields in the subsurface for a given conductivity distribution and a range of frequencies.Numerical modeling approaches such as finite difference (FD), finite element (FE), have been developed and applied as the process of forward modeling for 2D magnetotelluric regularized inversion [4][5][6][7][8].The FD method based upon the differential form of the partial differential equations (PDEs) is to be solved.Each derivative can be replaced with an approximate difference formula and the computational domain is usually divided into rectangular cells.The efficiency and accuracy of the FD method for modeling magnetotelluric responses were proven for the 2D geo-electric model [9,10].The FE method is another numerical approach that is often used for 2D and 3D magnetotelluric modeling, which involves assumed functional forms for the model and fields in small regions of specified geometry [11][12][13][14].
The FE method is more accurate for modeling a complicated medium, especially for topography, because of the greater flexibility of mesh discretization.Meanwhile, there are some different approaches to the numerical approximation of the 2D magnetotelluric forward problem [15][16][17][18].There are some comprehensive reviews of these numerical methods and their advantages and disadvantages [19,20].The references stated here represent a general view on magnetotelluric modeling.
In this paper, a non-uniform mesh approach in the FD numerical method was presented for a general 2D magnetotelluric modeling.The main challenge in applying the method for solving the magnetotelluric boundary value problem is to calculate spatial derivatives.To verify the accuracy of the FD forward algorithm, the resulting numerical was compared to both an analytical solutions and the FE numerical solutions.

Electromagnetic Equations
Considering a right-handed coordinate system, with z-axis pointing downwards and x-axis along geo-electrical strike, and assuming time dependence as e −iωt and neglecting displacement currents, electromagnetic Maxwell's equations can be expressed in the frequency domain as [21] ∇ where E is electric field and H is magnetic field.ω is an angular frequency (ω = 2π/T, T is a period).µ is the magnetic permeability in free space and µ = µ 0 = 4π × 10 −7 H/m.σ is the electric conductivity in (S/m) and is varied only in the direction of the y-axis and z-axis; i.e., σ = σ(y, z).
For a 2D conductivity structure assuming the x-axis is the geo-electrical strike direction (i.e., ∂E/∂x ≡ 0 and ∂H/∂x ≡ 0), expanding the curl operators in Equations ( 1) and (2), the governing Maxwell's equations are given by and These two modes are commonly referred to as transverse magnetic (TM) mode and transverse magnetic electric (TE) mode.We calculate the magnetic field H x (TM-mode) or electric field E x (TE-mode) parallel to the geo-electrical strike of the conductivity model.According to Equations ( 3) and ( 4), the associated partial differential equations can be written as ∂ ∂y ∂ ∂y According to these relations, the magnetotelluric impedances for the TM and TE modes are The associated apparent resistivities and impedance phases can be calculated as

Boundary Conditions
For the TE-mode, the numerical domain is composed of both the air space and the subsurface earth.For the TM-mode, the magnetic field H x is nearly unchanged in the air space, so the air space can be eliminated from the computational domain [22].To complete the boundary problem of the H-polarization, we must supply the boundary conditions for the magnetic field component on the outer boundaries.We restrict the computational domain for Equations ( 5) and ( 6) to 2D bounded domain Ω = [y min , y max ] × [z min , z max ], as shown in Figure 1.The computational domain is reduced to an isolated rectangular block with suitable boundary conditions, Then, we can divide the boundary of the computational domain Ω can be divided into four parts expressed as follows: Γ 1 = {(y, z)|y min < y < y max , z = z min }, Γ 2 = {(y, z)|z min < z < z max , y = y min }, Γ 3 = {(y, z)|z min < z < z max , y = y max }, and Γ 4 = {(y, z)|y min < y < y max , z = z max }.We impose Dirichlet Γ 1 , Neumann Γ 2 , Neumann Γ 3 , and Robin Γ 4 boundary conditions on the computational domain.
According to these relations, the magnetotelluric impedances for the TM and TE modes are The associated apparent resistivities and impedance phases can be calculated as

Boundary Conditions
For the TE-mode, the numerical domain is composed of both the air space and the subsurface earth.For the TM-mode, the magnetic field x H is nearly unchanged in the air space, so the air space can be eliminated from the computational domain [22].To complete the boundary problem of the H-polarization, we must supply the boundary conditions for the magnetic field component on the outer boundaries.We restrict the computational domain for Equations ( 5) and ( 6

FD Representation for 2D Magnetotelluric Problem
Before we can solve Equations ( 5) and ( 6), it is necessary to split the computational area into a mesh with (N zb + N za ) × N y cells, where N y is the number of cells in the horizontal direction and N zb and N za denote the number of cells in the vertical direction for the subsurface and the air, respectively.A resistivity value is assigned to each cell or each node.
The FD approximation is the simplest form of electromagnetic modeling.In order to construct FD scheme for solving Equations ( 5) and ( 6), non-uniform meshes shown in Figure 2 are constructed.In the TM-mode, approximation of second-order derivatives for Equation ( 5) can be expressed as and respectively.A resistivity value is assigned to each cell or each node.The FD approximation is the simplest form of electromagnetic modeling.In order to construct FD scheme for solving Equations ( 5) and ( 6), non-uniform meshes shown in Figure 2 are constructed.In the TM-mode, approximation of second-order derivatives for Equation ( 5) can be expressed as Figure 2. Discretization for 2D geo-electric model with non-uniform meshes.
Substituting Equations ( 10) and (11) into Equation ( 5), the difference equation for TM-mode can be written as Substituting Equations (10) and (11) into Equation ( 5), the difference equation for TM-mode can be written as 2(σ i,j +σ i+1,j )(ui+1,j−ui,j) In the TE-mode, the difference equation resulting from approximation of Equation ( 6) is: and Algorithms 2018, 11, 203 5 of 14 Therefore, substituting Equations ( 13) and ( 14) into Equation ( 6), the difference equation for TE-mode can be written as The solution vector u, which represents the components H x or E x for the different mesh nodes, is rewritten by . . .
Considering boundary conditions of the 2D magnetotelluric forward problem, the resulting linear system of equations in the form gives a numerical solution of Equations ( 5) and ( 6).Where K is the system matrix containing electrical parameters σ, and the right-side column vector s contains information related to boundary conditions.Equation ( 17) can be solved by the direct solver for a sparse matrix in MATLAB.Finally, the magnetotelluric responses including impedance, apparent resistivity and phase at each site for each frequency are calculated by Equations ( 7)-( 9).The subroutines described here for calculating magnetotelluric responses, given in Appendices A and B, were developed in MATLAB.

Benchmark with Homogeneous Half-Space
To test the accuracy of our FD algorithm, a homogeneous half-space model is illustrated for TM-mode.The size of computational domain was designed as 20 km × 5 km.The model was homogeneous with a conductivity of 0.1 S/m.We supposed that there is one magnetotelluric site located on the surface (z = 0 km).The magnetotelluric responses of the magnetic field H x are computed in the frequency range 0.001 to 1000 Hz.The given domain was discretized into many meshes and nodes with non-uniform approach.According to the information of mesh generation, the number of nodes in y-axis and z-axis are set as 51 and 41, respectively.The total numbers of nodes generated were 2091.
The numerical results of the homogeneous half-space are shown in Figure 3.It is obvious that the magnetic field H x including its real part and imaginary part computed by the FD method with non-uniform meshes are good agree with the analytical solution.The numerical results indicate that our FD approach can provide very good accuracy for 2D magnetotelluric modeling.
meshes and nodes with non-uniform approach.According to the information of mesh generation, the number of nodes in y-axis and z-axis are set as 51 and 41, respectively.The total numbers of nodes generated were 2091.
The numerical results of the homogeneous half-space are shown in Figure 3.It is obvious that the magnetic field Hx including its real part and imaginary part computed by the FD method with non-uniform meshes are good agree with the analytical solution.The numerical results indicate that our FD approach can provide very good accuracy for 2D magnetotelluric modeling.

Comparison of FD Results and Analytical Solutions
The two-layered model is used as a test example for the comparison of FD numerical solution and analytical solution.It is supposed that the 1D geo-electric model includes two resistivity layers.The thickness of the top layer is 1 km, and its resistivity value is assumed to be 10 ohm-m.For the bottom layer, its thickness and resistivity value are assumed to be 1 km and 100 ohm-m.The computational size of the model is set as 20 km × 2 km.
The results of the two-layered model are shown in Figure 4.It is obvious that relative errors of the apparent resistivities and the impedance phases were both increased when the frequency was increased, with the maximum error at the longest frequency 1000 Hz and the minimum at the shortest frequency 0.001 Hz.At the highest frequency, the vertical grid spacing size at earth's surface must be small enough so that the linear approximation of the electric field is reasonably close to the exponential decay of the electromagnetic fields.Through the numerical examples, the first vertical

Comparison of FD Results and Analytical Solutions
The two-layered model is used as a test example for the comparison of FD numerical solution and analytical solution.It is supposed that the 1D geo-electric model includes two resistivity layers.The thickness of the top layer is 1 km, and its resistivity value is assumed to be 10 ohm-m.For the bottom layer, its thickness and resistivity value are assumed to be 1 km and 100 ohm-m.The computational size of the model is set as 20 km × 2 km.
The results of the two-layered model are shown in Figure 4.It is obvious that relative errors of the apparent resistivities and the impedance phases were both increased when the frequency was increased, with the maximum error at the longest frequency 1000 Hz and the minimum at the shortest frequency 0.001 Hz.At the highest frequency, the vertical grid spacing size at earth's surface must be small enough so that the linear approximation of the electric field is reasonably close to the exponential decay of the electromagnetic fields.Through the numerical examples, the first vertical mesh size should be suggested approximately 1/3 of the shortest skin depth of the top layer.
The two-layered model is used as a test example for the comparison of FD numerical solution and analytical solution.It is supposed that the 1D geo-electric model includes two resistivity layers.The thickness of the top layer is 1 km, and its resistivity value is assumed to be 10 ohm-m.For the bottom layer, its thickness and resistivity value are assumed to be 1 km and 100 ohm-m.The computational size of the model is set as 20 km × 2 km.
The results of the two-layered model are shown in Figure 4.It is obvious that relative errors of the apparent resistivities and the impedance phases were both increased when the frequency was increased, with the maximum error at the longest frequency 1000 Hz and the minimum at the shortest frequency 0.001 Hz.At the highest frequency, the vertical grid spacing size at earth's surface must be small enough so that the linear approximation of the electric field is reasonably close to the exponential decay of the electromagnetic fields.Through the numerical examples, the first vertical mesh size should be suggested approximately 1/3 of the shortest skin depth of the top layer.

Comparison of FD Results and FE Solutions
In this section, the reliability of the non-uniform meshes FD algorithm is confirmed by the COMMEMI 2D-0 testing model.The COMMEMI 2D-0 model, proposed by Zhdanov et al. for comparing of numerical modeling methods for electromagnetic induction [19], is illustrated in Figure 5.With a testing period of 30s, the simulated apparent resistivities by the presented FD approach are compared with the FE solution [7] and the averaged volume integral results from the COMMEMI project [19].The computational size of magnetotelluric domain was designed as 100 km × 100 km.Using the technique of the non-uniform meshes, we set discrete elements in given domain as 100 × 80 (i.e., Ny = 100 and Nz = 80), and extended 50 km to air space for TE-mode.Figure 6 shows the results of the FD computations along a profile.With a testing period of 30 s, the simulated apparent resistivities by the presented FD approach are compared with the FE solution [7] and the averaged volume integral results from the COMMEMI project [19].The computational size of magnetotelluric domain was designed as 100 km × 100 km.Using the technique of the non-uniform meshes, we set discrete elements in given domain as 100 × 80 (i.e., N y = 100 and N z = 80), and extended 50 km to air space for TE-mode.Figure 6 shows the results of the FD computations along a profile.
One can see (Figure 6) that the results calculated by the FD method and the FE method are in close agreement to each other.It shows that the apparent resistivities have no significant difference.Compared to the FE solution, the average relative error for the apparent resistivity is 0.05, which are in the acceptance limits.Meanwhile, our FD numerical results were quite comparable those of the COMMEMI 2D-0 projects.
With a testing period of 30s, the simulated apparent resistivities by the presented FD approach are compared with the FE solution [7] and the averaged volume integral results from the COMMEMI project [19].The computational size of magnetotelluric domain was designed as 100 km × 100 km.Using the technique of the non-uniform meshes, we set discrete elements in given domain as 100 × 80 (i.e., Ny = 100 and Nz = 80), and extended 50 km to air space for TE-mode.Figure 6 shows the results of the FD computations along a profile.

Conclusions
The finite difference method with non-uniform meshes has been adapted to simulate the 2D magnetotelluric responses.We presented the calculation formula of this approach from the boundary value problem of electromagnetic fields.In the first investigation, the benchmark of a finite-difference algorithm was tested by a homogeneous half-space model and the approach of non-uniform meshes can provide very good accuracy for 2D magnetotelluric modeling.Compared to the analytical solutions for a 1D geo-electric model, the relative errors of the apparent resistivity and the impedance phase were both increased when the frequency was increased.Meanwhile, the relative errors were reduced when the mesh size was reduced.The reason can be attributed to the accuracy of the FD calculation for the apparent resistivity and the impedance phase depends on the vertical grid spacing size of the near-surface.Furthermore, we presented the COMMEMI 2D-0 model in comparison to results derived by the FE calculation.Both results are in close agreement to each other.These comparisons demonstrate that our FD algorithm can be valid and reliable.

3 
in Figure1.The computational domain is reduced to an isolated rectangular block with suitable boundary conditions, Then, we can divide the boundary of the computational domain  can be divided into four parts expressed as follows: , and Robin 4  boundary conditions on the computational domain.

Figure 1 .
Figure 1.A sketch of computational domain for modeling 2D magnetotelluric responses.Figure 1.A sketch of computational domain for modeling 2D magnetotelluric responses.

Figure 1 .
Figure 1.A sketch of computational domain for modeling 2D magnetotelluric responses.Figure 1.A sketch of computational domain for modeling 2D magnetotelluric responses.

Figure 3 .
Figure 3. FD solution of magnetic field in homogeneous half-space model.(a) Numerical solution of real component; (b) Numerical solution of imaginary component; (c) Analytical solution of real component; (d) Analytical solution of imaginary component; (e) Relative error of real component; (f) Relative error of imaginary component.

Figure 3 .
Figure 3. FD solution of magnetic field in homogeneous half-space model.(a) Numerical solution of real component; (b) Numerical solution of imaginary component; (c) Analytical solution of real component; (d) Analytical solution of imaginary component; (e) Relative error of real component; (f) Relative error of imaginary component.

Figure 4 .
Figure 4. Comparison of FD solution and analytical solution for the two-layered model.(a) Apparent resistivities of TM-mode; (b) Apparent resistivities of TE-mode; (c) Phases of TM-mode; (d) Phases of TE-mode.

Figure 4 .
Figure 4. Comparison of FD solution and analytical solution for the two-layered model.(a) Apparent resistivities of TM-mode; (b) Apparent resistivities of TE-mode; (c) Phases of TM-mode; (d) Phases of TE-mode.

4. 2 .
Comparison of FD Results and FE SolutionsIn this section, the reliability of the non-uniform meshes FD algorithm is confirmed by the COMMEMI 2D-0 testing model.The COMMEMI 2D-0 model, proposed by Zhdanov et al. for comparing of numerical modeling methods for electromagnetic induction[19], is illustrated in Figure5.It consists of three segments of different conductivities with horizontal contrasts of 1:10 and 2:1, lying on a perfectly conducting basement.Algorithms 2018, 11, x FOR PEER REVIEW 9 of 15 It consists of three segments of different conductivities with horizontal contrasts of 1:10 and 2:1, lying on a perfectly conducting basement.

Figure 6 .
Figure 6.Comparison of FD solution and FE solution for the COMMEMI 2D-0 model.The upper figure displays the apparent resistivity of TM-mode, the lower figure represents the TE-mode.

Figure 6 .
Figure 6.Comparison of FD solution and FE solution for the COMMEMI 2D-0 model.The upper figure displays the apparent resistivity of TM-mode, the lower figure represents the TE-mode.