Improved Integral Equation Method for Rapid 3-D Forward Modeling of Magnetotelluric

: Computational cost tremendously restricts the wide application of conventional integral equation (IE) method in large-scale magnetotelluric (MT) modeling. A couple of obstacles limit the developments of traditional MT modeling based on the IE method. They are: O (N2) space complexity of memory requirements for storing coefficients of dense matrix; singularity of Dyadic Green’s function; low efficiency of using digital filtering, such as Hankel transform, to calculate the Bessel function integral within the dyadic Green’s function, as well as inefficiency of accumulative calculation of 3-D discrete convolution. To solve these problems, we use an analytical formula instead of the Hankel transform to compute the integral of the Bessel function and replace a block cell by a spherical cell with the same volume to integrate through the singularity. Because the coefficient matrices are symmetric and antisymmetric three-level block-Toeplitz (BT) and Toeplitz + Hankel matrices, only non-redundant entities of the matrix are computed and stored. Afterwards, 3-D fast Fourier transform (FFT) is used to expedite matrix–vector multiplication at each successive iteration when using the contraction iterative method to solve the system of equations, which decreases memory and time consumption sharply compared with the traditional IE method.


Introduction
Magnetotelluric is playing a more and more considerable role in the interpretation of distributions of subsurface electrical conductivity, especially for large-scale surveys, such as SinoProbe [1], USArray [2], LITHOPROBE [3] and AusLamp [4]. Forward modeling as one of the important parts for interpretation of MT data has been developed for decades. The relative merits of major numerical modeling approaches (e.g., finite element (FE), finite difference (FD), finite volume (FV) and integral equation (IE) method) have been discussed in detail in previous studies [5][6][7][8][9][10][11]. From recent studies, we can see that one aim of forward modeling algorithm is to obtain electromagnetic (EM) responses of realistic geological models. Consequently, researchers developed unstructured meshes-based FE or FV methods [12][13][14][15]. By these two approaches, we are actually able to calculate EM responses of complicated models, yet as described in aforementioned publications, boundary conditions and extensional meshes are required. Compared with differential methods, as mentioned many times in previous papers, IE does not require to discrete whole modeling domain. In addition, with the development of non-uniform fast Fourier transform (NUFFT) [16] and unstructured meshes-based IE solver [17], there is a potential to develop a fast algorithm by unstructured meshes based IE method. Therefore, it is significant to establish the framework of MT forward modeling based on the uniform fast Fourier transform presented in this study.
Zhdanov et al. [18] introduced some pioneering works and made a throughout reviewed of the IE method before 2006. Afterwards, researchers continued to promote the development of IE method, which can be seen in the review work of Everett [19]. The efficiency is one of the vital criteria to assess a forward solver. Hence, some researchers managed to reduce the computational load by calculating and storing only a part of the whole matrix [20]. Meanwhile, some others used fast methods, such as the fast multipole method [21][22][23], the wavelet-based method [24,25] and the impedance matrix location method [26] to accelerate their algorithms. Another problem is that the coefficient matrix may be ill-conditioned, especially for the large conductivity contrast model. It results in slow convergence or divergence of iterative methods [27]. To overcome this obstacle, the contraction operator was deduced or used by some researchers [28][29][30][31][32].
There are some factors affecting accuracy and efficiency in conventional 3-D MT modeling based on the IE method. For example, the singularity of the dyadic Green's function [33] and the integration of Bessel function within the dyadic Green's function. In addition, the coefficient matrix is a full matrix, which is practically impossible to store and directly solve for large number of cells, while the iterative method is more suitable for this case. The dense matrix-vector multiplication also costs a lot of time and memory when using iterative method. Hursán and Zhdanov [27] have reviewed several approaches to solve the memory requirement and computation for full matrix. Hohmann [33] simply stored the upper or lower triangle matrix by neglecting the asymmetries from image currents and charges. Nevertheless, it is more reasonable for considering the inherent asymmetries. Hence, the application of the IE method is restricted to small-scale and simple structures MT modeling due to the computing burden. Nonetheless, the IE method is one of the accurate approaches for electromagnetic simulation [27,[33][34][35]. In addition, only requiring the discretization of anomalous region decreases the number of unknowns, although the coefficient matrix is a dense matrix. Meanwhile, the boundary conditions are naturally satisfied by Green's functions. Another superiority is that the coefficient matrix formed from the integral of Green's function can be a Toeplitz or Hankel matrix due to the translational invariance. Additionally, the integral of Green's function and current density is the integral of convolution type. These two factors enable us to use FFT to compute matrix-vector multiplication efficiently during each iteration when using the iterative method.
To improve the efficiency, this study uses the analytical formula to calculate the Bessel function integral, and spherical cell integration [33] to the integral Green function with current density. Combined with contraction operator [27,32], FFT is used to speed up matrix-vector multiplication [36][37][38].

IE Method Foundation for MT Modeling
The procedure of MT modeling using the IE method has reached a fairly mature state [11,18,33,34,39]. We will not expound more than what is needed here. The governing equation of electric fields becomes where E and E p are the total and incident electric fields, respectively. G E is electrical dyadic Green's function. r and r are coordinates related to observation points and sources, respectively. v is the volume of the anomalous body, and ∆σ denotes the conductivity discrepancy between the abnormal body and half-space. Note that the electrical conductivity of the subsurface media is by nature characterized as a complex second order symmetric tensor [40][41][42]. This paper primarily exhibits a fast algorithm for MT forward modeling, and hence, only the real and isotropic conductivities are considered for simplicity. It is easy to take into account the anisotropy in the IE method by replacing the scalar conductivity in Equation (2) with a complex tensor, which does not increase the computational complexity. The complex and anisotropic conductivities will be considered in the future study.
To remain as the Toeplitz and Hankel matrices, Equation (1) is divided into three parts denoting three components of electric fields: Figure 1a shows a scattered domain comprising all anomalous bodies. Block cells are used to discretize the modeling domain with N x , N y and N z cells in x-, y-and z-directions, respectively. Cells are numbered from left to right, front to back and up to down, while only anomalous bodies are discretized in Figure 1b. After discretization, Equation (3)(4)(5) can be written as: where The W in the Equations (9)-(17) can be divided into primary and secondary parts. The primary part, denoted as W ∼p E , arises from subsurface currents and charges. The secondary part, denoted as W ∼s E , arises from image currents and charges of earth-air interface, giving: where the symbol~denotes components xx to zz. Hohmann [33] has given an explicit form of the dyadic Green's function, and we rewrite it in a form that makes it easier for readers to find the Toeplitz and Hankel matrices. A more explicit form of Equation (18) is: Three-level block Toeplitz matrix arises from the first part of the right-hand side in Equation (19) which is a translational invariant form, while the second part yields threelevel block Toeplitz + Hankel matrices [43], namely, the first level is a Hankel matrix, and the second and third levels are Toeplitz matrices. In this study, the first level of these two kinds of matrices are N zr × N zs block matrices related to the number of electric field points and sources points in z-direction; analogously, the second level corresponds to x-direction, and the third level corresponds to y-direction. Taking the advantage of the block Toeplitz matrix, a matrix with constant diagonal elements, and the Hankel matrix, a matrix with constant anti-diagonal elements, we only need to compute and store non-redundant elements of the dense coefficient matrix. Specifically, the dense coefficient matrices W xxp E , W yyp E and W zzp E in Equations (9), (13) and (17)  . For symmetric and antisymmetric matrices, we just need to calculate the column or row of each level, namely, N x × N y × N z entities, and extend them to (2N x − 1) × 2N y − 1 × (2N z − 1) before using 3-D FFT to calculate matrixvector multiplication. Accordingly, the efficiency of computing coefficient matrices can be significantly improved compared with that of calculating the full matrices. For the discretization scheme in Figure 1a, more cells are required, but all EM fields at concerned positions can be obtained after solving the system of equations. It is also more suitable for inversion, since the amount of subsurface anomalous bodies is unknown. Although the case in Figure 1b calls for less cells, the structure of the coefficient matrix becomes more complicated. If one wants to use FFT for fast calculation, the matrix should be divided into Na×Na sub-matrices, where Na is the number of anomalous bodies. Sub-matrices represent interactions between each anomalous bodies, including self-interaction. Sub-matrices formed from self-interaction cases are square matrices, while others may be square or non-square matrices, which depends on the number of cells in each anomalous body. If the number of cells of two anomalous bodies are equivalent, the matrix is squared; otherwise, it is non-squared. tions between each anomalous bodies, including self-interaction. Sub-matrices formed from self-interaction cases are square matrices, while others may be square or non-square matrices, which depends on the number of cells in each anomalous body. If the number of cells of two anomalous bodies are equivalent, the matrix is squared; otherwise, it is nonsquared.

Improved Treatments
This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

Analytical Method for Computation of Bessel Function
Another time-consuming factor when calculating coefficient matrices is two Bessel function integrals [33] that are relevant to the air-earth interface. They are given by: where k 1 = −iωµ 0 σ b , ω is angular frequency, µ 0 denotes magnetic permeability of free space, σ b is the conductivity of background earth, u 1 = λ 2 − k 2 1 . To fast compute the coefficient matrix, analytical formula instead of the digital filtering algorithm [44] are used to calculate Equations (20) and (21).
According to Sommerfeld integral and the integral we can obtain the analytical formula of γ 1 and γ 2 [45,46] (see Appendix A) where In addition, where , g = ik 1 r 2 /R s . I n and K n are the first and second kind modified Bessel functions for order n with variables P and Q, respectively, which can be solved by polynomial approximation [47].

Rapid Implementation of Coefficient Matrix-Vector Multiplication
According to the aforementioned special matrices, we are able to implement fast calculation of matrix-vector multiplication using 3-D FFT. In order to use FFT, the Toeplitz + Hankel matrix needs to be transformed into a Toeplitz matrix as follows: reverse the columns order and rows order of the matrix and vector, respectively Only the first level of the Toeplitz + Hankel matrix needs to be transformed into a Toeplitz matrix in the case of this study. Subsequently, as for submatrices formed from self-interaction of anomalous bodies, we can reconstruct the coefficient matrix following the approach applied to gravity modeling [37] by changing 2-D FFT to 3-D FFT, which described the procedures to rearrange the coefficient matrix in detail. For submatrices of interaction of different anomalous bodies that method is not available anymore. If we use the discretization scheme in Figure 1a, that method still works. In the following sections, we focus on the discretization scheme shown in Figure 1b. To be adequate to all submatrices mentioned above, we presented a new method to reconstruct the non-redundant values of the full matrix as follows. Suppose we have a random 3-level asymmetric Toeplitz matrix 3 8 9 7 8 4 9 2 1  3 1 2 2 3 8 2 7 8 5 9 2  5 3 1 0 2 3 5 2 7 0 5 9  2 4 5 1 2 3 8 4 3 7 8 4  6 2 4 3 1 2 1 8 4 2 7 8  7 6 2 5 3 1 7 1 8 5 2 7  4 3 1 6 7 4 1 2 3 3 8 9  4 4 3 3 6 7 3 1 2 2 3 8  0 4 4 5 3 6 5 3 1 0 2 3  5 2 1 4 3 1 2 4 5 1 2 and a random vector x = 1 2 3 4 5 6 7 8 9 10 11 12 T .
We only need to calculate the following non-redundant values of matrix A and sort them in a 3-D array Compared with the method of Vogel [48] and Chen and Liu [37], we need no more extra zeros to pad array a. Using the discrete convolution theorem here fftn and ifftn denote 3-D FFT and inverse 3-D FFT, respectively. The sign · is the dot product operator. Extracting values from f in Equation (31), one can obtain the desired solutions. For the non-square matrix, an identical matrix-vector multiplication procedure can be used. The successive iterative method or Born series [27,32,49] is used to solve Equations (6)- (8), and the fast algorithm above can be used to speed up matrix-vector multiplication at each iteration. In consideration of the convergence of the iterative method for large conductivity contrast and large size models, the contraction operator established by Zhdanov and Fang [32] needs to be added to Equations (6)-(8) during each iteration, as shown in the following flow chart.
Neglecting displacement currents, the contraction operator α and β in Figure 2 are as follows: Consequently, magnetic fields can be obtained by taking the curl of the solved electric fields following Hohmann [33].

COMMEMI 3D-1A Model
Provide initial electric fields Consequently, magnetic fields can be obtained by taking the curl of the solved electric fields following Hohmann [33].

COMMEMI 3D-1A Model
The COMMEMI 3D-1A model ( Figure 3) with a large contrast in conductivity (ρ 1 /ρ 2 = 200) [50] is used to validate our developed IE solver. We discretize the anomalous body into 20 × 10 × 20 cells with size 100 m × 100 m × 100 m in x-, yand z-directions, respectively. Apparent resistivity and phase for 0.1 Hz after 10 iterations are shown in  Figure 4 shows that our results agree well with the results of previous papers [5,50] with a small number of discretization cells and low computational cost. Figure 5 presents the convergence of successive iteration approach when solving electric fields. Consequently, magnetic fields can be obtained by taking the curl of the solved elect fields following Hohmann [33].

COMMEMI 3D-1A Model
The COMMEMI 3D-1A model ( Figure 3) with a large contrast in conductiv ( ⁄ = 200) [50] is used to validate our developed IE solver. We discretize the anom lous body into 20 × 10 × 20 cells with size 100 m × 100 m × 100 m in x-, y-and z-dire tions, respectively. Apparent resistivity and phase for 0.1 Hz after 10 iterations are show in  Figure 4 shows that our results agr well with the results of previous papers [5,50] with a small number of discretization ce and low computational cost. Figure 5 presents the convergence of successive iteration a proach when solving electric fields.    In this part, we compare the accuracy and efficiency of calculating Equations (20) and (21) using the digital filtering method [44] as well as the formulas described in Equations (24) and (25). It is meaningless to compare the tional cost of and for one cell. Hence, we divide COMMEMI 3D-1A into 50 cells in the x-, y-and z-direction, respectively, which implies (2 × 50 − 1) × 1) × (2 × 50 − 1) = 970,299 entities to be calculated. It costs 39.0 s and 61.3 s to , , respectively, by digital filtering methods. While they are merely 0.85 s through analytical formulas. Figure 6 shows 100 random real and imaginary p and among 970,299 entities by the two methods. They are practically iden each other. Making use of the analytical formula can remarkably decrease the t culating and especially for plenty of cells. In this part, we compare the accuracy and efficiency of calculating γ 1 and γ 2 in Equations (20) and (21) using the digital filtering method [44] as well as the analytical formulas described in Equations (24) and (25). It is meaningless to compare the computational cost of γ 1 and γ 2 for one cell. Hence, we divide COMMEMI 3D-1A into 50 × 50 × 50 cells in the x-, y-and z-direction, respectively, which implies (2 × 50 − 1) × (2 × 50 − 1) × (2 × 50 − 1) = 970, 299 entities to be calculated. It costs 39.0 s and 61.3 s to calculate γ 1 , γ 2 , respectively, by digital filtering methods. While they are merely 0.85 s and 0.89 s through analytical formulas. Figure 6 shows 100 random real and imaginary parts of γ 1 and γ 2 among 970,299 entities by the two methods. They are practically identical with each other. Making use of the analytical formula can remarkably decrease the time of calculating γ 1 and γ 2 especially for plenty of cells.

Dublin Test Model 1
In this part, the Dublin Test Model (MTNet-Home, https://www.mtnet.info/main/ (accessed on 15 March 2022)) provided at the MT 3D inversion workshop in Dublin is used to test the applicability of the presented algorithm to complex conductivity structures. The model is composed of three different blocks in a homogeneous 100 Ωm half-

Dublin Test Model 1
In this part, the Dublin Test Model (MTNet-Home, https://www.mtnet.info/main/ (accessed on 15 March 2022)) provided at the MT 3D inversion workshop in Dublin is used to test the applicability of the presented algorithm to complex conductivity structures. The model is composed of three different blocks in a homogeneous 100 Ωm half-space with the largest contrast in resistivity of 10,000/1. The sizes, locations and resistivities of three anomalous bodies are shown in Figure 7 and Table 1. We use cubic cells to discretize the three anomalous domains into 40 × 5 × 15, 15 × 25 × 5 and 15 × 25 × 30 cells in the x-, yand z-directions, respectively. With 21 periods ranging from 0.1 s to 10,000 s, we compare our apparent resistivity with those provided by Miensopust et al. [51], which includes results of 11 codes. We choose 2-11 results for comparison because the first result shows a big difference with others. In addition, the result of edge-based finite element [52] denoted as green rhombus is also used for comparison. As shown in Figure 8, all apparent resistivity and phase components, xx, xy, yx and yy, at the observational site (0, 0, 0) m after 30 iterations are closed to most of the 11 results, and the time consumption is about 25 s per period. The relative residual norms between two adjacent successive iterations for different periods are shown in Figure 9, which shows the good convergence rate of this iterative method. This example shows that our algorithm is able to solve the complex and high contrast of electrical conductivity structures efficiently.     Extend in x (km) Extend in y (km) Extend in z (km)  Figure 8. Comparison of apparent resistivity and phase components at the station located at coordinate (0, 0, 0) m of DTM1 with 21 periods ranging from 10 s to 10,000 s (4 periods per decade). Our results are compared with those provided by other authors (see [51] for details).
Minerals 2022, 12, x FOR PEER REVIEW 11 of 16 Figure 8. Comparison of apparent resistivity and phase components at the station located at coordinate (0, 0, 0) m of DTM1 with 21 periods ranging from 10 s to 10,000 s (4 periods per decade). Our results are compared with those provided by other authors (see [51] for details). In order to test the time and memory consumption of the presented algorithm, four discretization strategies shown in Table 2 are carried out for the DTM1 model. Solutions of apparent resistivities are consistent with those shown in Figure 8 and are not presented here. This experiment shows the capacity of our algorithm to solve millions of cells with acceptable time and memory use.  In order to test the time and memory consumption of the presented algorithm, four discretization strategies shown in Table 2 are carried out for the DTM1 model. Solutions of apparent resistivities are consistent with those shown in Figure 8 and are not presented here. This experiment shows the capacity of our algorithm to solve millions of cells with acceptable time and memory use.

COMMEMI 3D-2A Model
The COMMEMI 3D-2A [50] model comprising a three-layer background is used to test that our algorithm is adequate to model layered-Earth in this section. There are two anomalous prisms embedded in the top layer as shown in Figure 10. In order to use simpler half-space Green's function compared with layered Green's function, we treat the layered-Earth as anomalous bodies. We compared our results with the FE values of Farquharson and Miensopust [5]. The compared apparent resistivities for 0.1, 0.01 and 0.001 Hz are shown in Figure 11, and they agree well with each other. In this way, it indeed requires more discretized cells; however, we can calculate the responses of models with layered background using the half-space Green's function. Consequently, the aforementioned algorithms that analytical computation of Bessel function integrals and fast calculation of matrix-vector multiplication can be exploited. In theory, the layered background should be infinitely extended in the horizontal direction, but the amplitude of secondary fields induced by anomalous bodies will be exponentially reduced to zero due to the skin effect. Thus, we can choose the appropriate horizontal extension of layered-Earth as background instead of infinite extension. Moreover, the layered-Earth does not always stretch far horizontally and continuously in realistic geologic earth, therefore, it is reasonable to take layered-Earth as anomalous bodies.     [5], and circles indicate the solutions of this study. The first to third rows are for 0.1 Hz, 0.01 Hz and 0.001 Hz, respectively.

Discussion
For a layered-Earth background, such as the COMMEMI 3D-2A model, the layered media Green's function is more suitable for the IE method. However, the expressions of the layered media Green's function are more complicated compared with the half-space case. There are two significant differences in their IE implementation. Firstly, the coefficient matrices resulting from discretizing the IE equation exhibit different structures. The matrices in the layered media case can be divided into two parts. One part is related to the direct wave from the whole space, namely, when the observational points and scatter sources are at the same background layers. For this part, the structure of the matrices is identical with the primary part ( ~ in Equation (18)

Discussion
For a layered-Earth background, such as the COMMEMI 3D-2A model, the layered media Green's function is more suitable for the IE method. However, the expressions of the layered media Green's function are more complicated compared with the halfspace case. There are two significant differences in their IE implementation. Firstly, the coefficient matrices resulting from discretizing the IE equation exhibit different structures. The matrices in the layered media case can be divided into two parts. One part is related to the direct wave from the whole space, namely, when the observational points and scatter sources are at the same background layers. For this part, the structure of the matrices is identical with the primary part W ∼p E in Equation (18) of the half-space case, they are symmetric or antisymmetric three-level Toeplitz matrices, while the other part is related to the reflected wave from different layers. Different from the half-space case, these matrices are not three-level Toeplitz + Hankel matrices, but they are two-level Toeplitz matrices when the observational points and scatter sources are at the same discretized layers, and one can use 2-D FFT to conduct matrix-vector multiplication for this case. It requires N 2 z times 2-D FFT and is more time-consuming compared with the half-space case. Secondly, it is hard to find the analytical formulas for the Bessel function integral in layered media case as we do in the half-space case. As a result, the digital filtering method has to be used for calculating the Bessel function integral which significantly increases the computational time.
For the coefficient matrices with a different structure in the layered media case, we have developed corresponding efficient algorithms to conduct matrix-vector multiplication, but this is beyond the scope of this study and will be discussed in detail in future work.

Conclusions
This study presents a fast algorithm to solve MT modeling based on the IE method. To improve computational efficiency, an analytical formula is used to calculate the integral of the Bessel function in the Green's function. Coefficient matrices are symmetric and antisymmetric, three-level block-Toeplitz and Toeplitz + Hankel matrices owing to the translational invariance and convolutional nature of the free-space Green's function. Taking advantage of this special matrix structure, we simply calculate and store N x × N y × N z entities and extend them to (2N x − 1) × 2N y − 1 × (2N z − 1) entities before using 3-D FFT to calculate the matrix-vector multiplication. Afterward, we construct these nonredundant entities into a well-organized array and use 3-D FFT to compute matrix-vector multiplication efficiently at each iteration. Conventional IE method is difficult to solve a large and complex modeling region due to the preceding obstacles. This study makes the IE method a more powerful numerical approach to electromagnetic modeling.

Acknowledgments:
The authors would like to thank Farquharson for providing his 3-D MT modeling results used for validating the algorithm in this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Analytical Formula of Bessel Function Integral
Equation (22) is the Sommerfeld integral and (23) is an identical equation. Equation (20) can be written as two parts: according to Equation (23), Equation (A3) becomes Equation (A2) can be written as Based on Equation (22), we have By calculating derivatives in Equations (A4) and (A6), we have γ 1 in Equation (24). Analogously, we can obtain γ 2 in Equation (25).