3D Airborne EM Forward Modeling Based on Time-Domain Spectral Element Method

: Airborne electromagnetic (AEM) method uses aircraft as a carrier to tow EM instruments for geophysical survey. Because of its huge amount of data, the traditional AEM data inversions take one-dimensional (1D) models. However, the underground earth is very complicated, the inversions based on 1D models can frequently deliver wrong results, so that the modeling and inversion for three-dimensional (3D) models are more practical. In order to obtain precise underground electrical structures, it is important to have a highly effective and efﬁcient 3D forward modeling algorithm as it is the basis for EM inversions. In this paper, we use time-domain spectral element (SETD) method based on Gauss-Lobatto-Chebyshev (GLC) polynomials to develop a 3D forward algorithm for modeling the time-domain AEM responses. The spectral element method combines the ﬂexibility of ﬁnite-element method in model discretization and the high accuracy of spectral method. Starting from the Maxwell's equations in time-domain, we derive the vector Helmholtz equation for the secondary electric ﬁeld. We use the high-order GLC vector interpolation functions to perform spectral expansion of the EM ﬁeld and use the Galerkin weighted residual method and the backward Euler scheme to do the space- and time-discretization to the controlling equations. After integrating the equations for all elements into a large linear equations system, we solve it by the multifrontal massively parallel solver (MUMPS) direct solver and calculate the magnetic ﬁeld responses by the Faraday's law. By comparing with 1D semi-analytical solutions for a layered earth model, we validate our SETD method and analyze the inﬂuence of the mesh size and the order of interpolation functions on the accuracy of our 3D forward modeling. The numerical experiments for typical models show that applying SETD method to 3D time-domain AEM forward modeling we can achieve high accuracy by either reﬁning the mesh or increasing the order of interpolation functions.


Introduction
Airborne electromagnetic (AEM) uses a moving platform to tow EM equipment for geophysical exploration. Because of the unnecessity of human access to the survey area, as is generally required by ground geophysics, this method has great potential working in areas of complex surface conditions, like mountains, forest-covered areas, or swamps. AEM methods can be divided into frequency-and time-domain methods based on the transmitting current passing through the transmitting coil, such as a half-sine or a trapezoid wave in time-domain or a harmonic wave in frequency-domain. Compared with the frequency-domain method, the time-domain AEM has strong capabilities in detecting the (SE) method can, however, achieve high modeling accuracy by mesh refinement and enhancement of the order of interpolation basis functions.
The spectral element (SE) method is a high-order weighted residual method for solving partial differential equations that combines the high precision characteristics of spectral method and the flexible discretization of the FE method. This method was first proposed by Patera for solving the fluid mechanics problems [29][30][31]. Two kinds of basis functions are commonly used in SE, namely the GLL (Gauss-Lobatto-Legendre) and GLC (Gauss-Lobatto-Chebyshev) polynomials. As for computational geophysics, Komatitsch and Vilotte used the SE method to simulate the propagation characteristics of the seismic wave field [32]. They formed the mass matrix by GLL basis functions and obtained the solution rapidly. Komatitsch et al. developed the open-source code SPECFEM3D_Cartesian based on SE method using the GLL basis functions [32][33][34][35]. Wang et al. used the SE method to model the seismic wave field in transversely isotropic media. They used the GLL basis functions and explicitly discretized time derivatives to improve computational efficiency [36]. In geophysical EM modeling, Yin et al. first introduced the SE method based on GLL polynomials into AEM modeling in frequency-domain, Huang et al. further applied this method for AEM anisotropic modeling [37,38], Liu et al. introduced the GLC spectral element method into 3D marine CSEM forward problem that greatly improved the modeling efficiency [39]. However, not much attention has been paid to the SE method for time-domain AEM forward modeling that is the focus of this research.
Another issue with time-domain EM modeling involves the time discretization and solution. From the time discretization, the 3D EM modeling can further be divided into three categories of frequency-time conversions, the Krylov subspace solution, and the discretized time iteration method. For the frequency-time conversion method, one needs first to calculate the EM responses for multiple frequencies, e.g., by FD, FE, or FV methods, and then uses inverse Fourier transform, Gaver-Stefest (GS) inverse Laplace transform, etc., to convert the EM responses from frequency-to time-domain [40,41]. These methods are easy to apply, but they need to calculate the EM responses at a wide frequency range, resulting in big computational cost. The Krylov subspace method solves the time-domain EM controlling equation in the Krylov subspace. It can directly solve the EM field and avoid the time discretization, so that high accuracy and efficiency can be achieved [42]. However, this method has the disadvantage that the basis functions are not easy to find and thus is difficult to be applied in EM inversions. The time iteration method starts directly from Maxwell's equations in time-domain. It uses the time difference to discretize the time derivatives in the controlling equations [43]. The most frequently used method is the backward difference that is unconditionally stable without much limitations on the time-steps, so that the numerical accuracy in the forward modeling can be assured [44,45].
In this paper, we combine the SE method with the backward time difference (called SETD method) based on GLC polynomials for 3D time-domain AEM modeling. Starting from the Maxwell's equations in time-domain, we derive the vector Helmholtz equation (i.e., the controlling equation) for the secondary electric field. We use the GLC vector interpolation basis functions to spectrally expand the electric and magnetic fields and discretize the controlling equation in space-and time-domain by the Galerkin's weighted residual method and backward difference scheme. After solving the final equations system by the multifrontal massively parallel Solver (MUMPS), we calculate the magnetic field by the Faraday's law. Further, we check the accuracy of the algorithm by comparing our results with the semi-analytical solutions for a 1D layered earth model. We also demonstrate the flexibility of our SETD algorithm in terms of the accuracy improvement by refining the mesh or increasing the order of the interpolation functions. Finally, we perform 3D forward simulations on typical models and analyze the characteristic AEM responses.

3D Forward Modeling
The Maxwell's equations in time-domain can be written as where E(r, t) and H(r, t) are respectively the electric and magnetic field at position r and time t, B(r, t) is the magnetic induction vector, D(r, t) is the electric flux density, J(r, t) is the current density, while J s (r, t) is the source current density. Using the constitutive equations J(r, t) = σE(r, t), D(r, t) = εE(r, t), B(r, t) = µH(r, t) and neglecting the displacement current, we have where σ, ε, and µ represent the conductivity, dielectric constant, and magnetic permeability, respectively. In this paper, we assume that the permittivity of the earth ε = 8.85 × 10 −12 F/m and the magnetic permeability µ = 4π × 10 −7 H/m. In the following, we use the secondary electric field to formulate the AEM forward problem. For this purpose, we separate the electric field into primary and secondary parts E = E b + E s , where E s denotes the secondary electric field, while E b denotes the primary electric field produced by the background conductivity σ b . After the primary/secondary separation, we obtain from Equations (5)-(8) the following vector Helmholtz equation, i.e., Furthermore, we divide the model area into a series of non-overlapping regular hexahedral elements and apply the Galerkin weighted residual method to each element, in which we choose the weighting functions as the basis functions, so that we can write the residual for the eth element as Assuming that there are NL edges in the eth element, Nx, Ny, and Nz are the orders of the basis functions respectively in the x-, y-, and zdirection, then we have NL = Nx × (Ny + 1) × (Nz + 1) + (Nx + 1) × Ny × (Nz + 1) + (Nx + 1) × (Ny + 1) × Nz. The vector basis function Φ(r) for our SE method only contains the basis functions in the x-, y-, and z-direction, i.e., then the secondary and primary electric field at any position in the cell can be expressed as where Φ e j is the interpolation basis functions, while E e sj and E e bj are respectively the tangential secondary and primary electric field corresponding to the jth edge. Substituting Equations (12) and (13) into (10), we obtain via partial integration The interval for the integration in Equation (14) is (x, y, z) ∈ [x e , x e+1 ] × [y e , y e+1 ] × [z e , z e+1 ], while the SE basis functions composed of GLC orthogonal polynomials are defined in [−1, 1]. Thus, we need to convert the coordinates of the eth element in the physical domain (x, y, z) ∈ [x e , x e+1 ] × [y e , y e+1 ] × [z e , z e+1 ] into the coordinates of the standard element (ξ, η, ζ) ∈ −1, 1 × −1, 1 × −1, 1 in the reference domain. To distinguish between the quantities in the physical domain and in the reference domain, we add a tilde "~" to the quantities in the physical domain. Then, the mapping relationship between the two domains can be written as where J and |J| are the Jacobian matrix and the corresponding determinant, respectively. The Jacobian matrix J can be expressed as Insertion of Equations (15)- (17) into (14) yields In the matrix format, Equation (18) can be written as where A e is the stiffness matrix for the element in the reference domain, while B1 e and B2 e are the mass matrices. These matrices can be expressed as The vector interpolation basis functions in the reference coordinate system can be written as In Equation (23), i represents the ith edge in the standard element of the reference domain, r, s, and t denote the index in three directions of the standard element, e ξ , e η , and e ζ represent the unit vectors along the ξ-, η-, and ζ-directions, respectively. The reduction of one order in each direction in Equation (24) is equivalent to calculating the gradients in these directions. This can not only reveal the directionality of the electric field, but can also ensure the continuity of the tangential electric field [46]. h represent the Nξ − 1, Nη, and Nζ order of the GLC orthogonal interpolation polynomials that can be expressed as where k, t = 0, 1, 2, · · · , N, T t (ξ) = cos (t · arccos (ξ)) is the tth order of Chebyshev polynomial, ξ k = − cos(kπ/N).
Substituting Equation (24) into (20)- (22), we can obtain the specific formats for calculating the stiffness and mass matrices (see Appendix A). Finally, we assemble the residual matrix in Equation (19) from all elements and let the weighted residual R = 0, so that we can get To discretize the time derivatives in Equation (27), we choose the unconditionally stable second-order of backward Euler scheme [24] and get Equation (28) shows that the secondary electric field at the time t i+2 is determined by the secondary electric field at the time t i+1 and t i and the primary electric field at the time t i+2 , t i+1 , and t i . For simplicity, we can rewrite Equation (28) as where S is the coefficient matrix, E is the secondary electric fields at all edges, while M is the source term. To ensure the unique solution of the forward problem, we apply the following initial condition and Dirichlet boundary condition to the solution, i.e.,

Accuracy Verification
To verify the accuracy of our 3D AEM forward modeling algorithm, we use our SETD method to calculate the EM responses for a layered earth model and compare the results with 1D semi-analytical solutions [2]. The computational environment for our modeling is an Intel(R) Xeon(R) CPU E5-2650 v3 @ 2.30 GHz (double processors). The layered model is as follows: the first layer has a conductivity of 0.01 S/m and a thickness of 40 m; the second layer has a conductivity of 1 S/m and a thickness of 40 m; the underlying half-space has a conductivity of 0.01 S/m; the air layer's conductivity is assumed to be 10 −8 S/m. The AEM system takes a concentric configuration with a unit dipole moment. The altitudes of the transmitting and receiving coils are 30 m. Figure 1 shows the model parameters. The transmitting current is a step wave. We divide the time into ten segments, with each being further divided into 100 small ones. The time intervals for each segment are respectively 1.0 × 10 −7 s, 2.0 × 10 −7 s, 4.0 × 10 −7 s, 8.0 × 10 −7 s, 1.6 × 10 −6 s, 3.2 × 10 −6 s, 6.4 × 10 −6 s, 1.28 × 10 −5 s, 2.56 × 10 −5 s, 5.12 × 10 −5 s. To demonstrate the advantage of SETD in the improvement of the modeling accuracy by either refining the mesh or enhancing the order of the interpolation basis functions, we design two schemes: (1) 40 m × 40 m × 20 m mesh combined with the third order of interpolation basis functions, (2) 60 m × 60 m × 20 m mesh combined with the fourth order of interpolation basis functions. Figure 2(a1)-(d1),(a2)-(d2) show the dbz/dt and bz responses and the relative errors against the semi-analytical solutions. It is seen from the figure that the EM fields dbz/dt and bz for the two combined schemes decay steadily with time, the overall errors are less than 1.2%. This demonstrates that our SETD algorithm has high accuracy. In order to further illustrate the high-precision characteristics of the algorithm presented in this paper, we compare the solutions of our SE method with those of the finite-element algorithm based on structured grids [47]. Figure 2(a3)-(d3),(a4)-(d4) shows the relative errors of both methods against the semi-analytical solutions. It is clearly seen that our SE method has much higher accuracy than the finite-element method based on structured grids.

Influence Factors of Modeling Accuracy
To explore the factors that affect the accuracy of our SETD algorithm, we simulate again the 3-layer model shown in Figure 1. The order of the interpolation basis functions and the mesh size for Figure 3(a1)-(d3) are variable. Analysis on Figure 3(a1)-(d1),(a2)-(d2) shows that when the order of the interpolation functions is fixed, the modeling accuracy is improved by refining the mesh. In contrast, the analysis on Figure 3 (a2)-(d2),(a3)-(d3) shows that when the size of the mesh is fixed, increasing the order of the interpolation functions can improve the accuracy as well. This implies that using SETD method for AEM forward modeling, one can improve the accuracy by either refining the mesh or increasing the order of the interpolation functions.

EM Responses for Typical Models
To analyze the characteristics of AEM responses in time-domain, we design the models of a horizontal plate, a vertical dyke, two vertical dykes in a homogeneous half-space. Figure 4 shows the horizontal dyke model with a dimension of 90 m × 90 m × 30 m, buried at a depth of 30 m and centered at (0,0,45 m). The conductivity of the dyke is 0.2 S/m, the background conductivity is 0.002 S/m. We model a flight system with a concentric configuration at an altitude of 30 m. The transmitting current is a step wave. The mesh size is 30 m × 30 m × 30 m and the third order of interpolation basis functions are used. Figure 5a,b show dbz/dt and bz responses for the horizontal plate. It is seen that (1) the existence of the conductive horizontal plate produces a symmetric and broad single-peak anomaly; (2) at early time, the EM field diffuses mainly in the upper resistive half-space, so that the EM signal is weak and broad; (3) with time passing, the EM field diffuses into the location of the conductive plate, so that strong EM signal is observed; (4) at later time, the EM field diffuses again into deep resistive earth, so that the EM signal becomes weak and broad again; (5) from the locations of half-peaks we can roughly determine the boundaries of the horizontal plate. For the computer environment mentioned in Section 3.1, the time for calculating a decay curve of 1000 time channels at a single survey station takes about 41.4 min.  To further demonstrate the high accuracy of our SETD algorithm for 3D AEM modeling, we also use the FE method based on unstructured grids proposed by Qi et al. to calculate the EM responses for the model in Figure 4 at a single survey point [45]. The comparison of these two results is shown in Figure 6. It can obviously be seen that the results of the two methods agree well with each other, which again illustrates that our SETD algorithm has high modeling accuracy.  Figure 8a,b shows the dbz/dt and bz responses for the vertical dyke model. It is seen that (1) the EM responses at the early time is similar to those in Figure 5. Since the EM wave at early time has high-frequency content, the vertical dyke appears as a thick dyke, so that we see a single-peak anomaly; (2) with time passing, the EM field diffuses into the location of the dyke, the influence of the dyke becomes more obvious. In opposite to a horizontal plate, the EM responses for a vertical dyke shows two peaks. There exists a minimum over the top of the dyke due to zero coupling between the transmitter-receiver pair and the dyke; (3) at later time, the EM field passes through the dyke, the EM signal diffuses into the underlying half-space and becomes broader and weaker; (4) the center of the dyke can be determined by the minimum location of the EM responses, while the boundaries of the dyke can be roughly estimated from the middle points between the maximum and minimum peaks. This helps identify the underground targets from the AEM anomalies. For the computer environment mentioned in Section 3.1, the time for calculating a decay curve of 1000 time channels at a single survey station takes about 36.5 min.
Next, we study the AEM responses for a model with two vertical dykes embedded in a homogeneous half-space (Figure 9). Each dyke has a dimension of 30 m × 120 m × 120 m and is buried at the depth of 30 m. The distance between the two dykes is 40 m. The conductivities of the two dykes are 0.2 S/m, while the background conductivity is 0.002 S/m. We discretize the model domain into the elements of 30 m × 30 m × 30 m and take the fourth order of interpolation basis functions. Figure 10a,b show the dbz/dt and bz responses. It is seen from the figures that the EM responses demonstrate as the superposition of two dykes. At the early time, the EM field shows a weak and broad anomaly, meaning that the two dykes are coupled and not distinguishable. With time passing the anomalies from two dykes become stronger and more distinguishable. Three peaks can be seen clearly, representing the EM field of two dykes. The anomaly in the middle is larger than the neighboring two anomalies as it is the superposed anomaly from the two dykes. Similarly, the centers and the boundaries of the dykes can be distinguished respectively from the minimum anomaly locations and the middle points between the maximum and minimum peaks. This further implies that from the anomalies of a time-domain AEM system, we can distinguish the anomalous bodies in the underground. For the computer environment mentioned in Section 3.1, the time for calculating a decay curve of 1000 time channels at a single survey station takes about 43.7 min.

Conclusions
In this paper, we have successfully developed and applied the SETD method to 3D AEM forward modeling in time-domain. Using the high-order of GLC interpolation basis functions, we can calculate analytically the stiffness and mass matrices for AEM forward problem and achieve high numerical accuracy. The backward Euler scheme for the time discretization, resulting in an unconditionally stable equation for EM fields, further improves the stability and accuracy of our 3D time-domain AEM forward modeling. Comparisons with the semi-analytical solutions for a layered earth model and with the finite-element solutions for 3D models demonstrated that our SETD method has high accuracy. In addition, the modeling accuracy of our SE method can be improved by either refining the mesh or increasing the order of the interpolation basis functions. This offers more possibilities for us to achieve accurate results in time-domain AEM modeling.
Numerical experiments on typical models show that the geometries and locations of the abnormal bodies in the underground can be identified from the AEM signal.
It must be pointed out that we have in this paper only used the regular hexahedral meshes in our SE modeling, so that we cannot yet model the complex underground structures and topography. Moreover, the modeling efficiency of our algorithm is not yet satisfying that will surely restrict its application in AEM inversions. All these will be our future research focus.

Acknowledgments:
We want to thank the three reviewers and the RS editors for their constructive suggestions that help clarify this paper. We want to thank Wei Huang and Yanfu Qi for offering their model results to validate our SETD algorithm.

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

Appendix A
The stiffness and mass matrices in Equations (20)-(22) contain 9 block matrices, i.e., where the elements of the block matrices can be written as t2 (ζ) · l z l z 4 8 l x l y l z dξdηdζ · l x l x 4 8 l x l y l z dξdηdζ 8 l x l y l z dξdηdζ 8 l x l y l z dξdηdζ , t2 (ζ) · l z l z ∂ξ l y l y 4 8 l x l y l z dξdηdζ