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

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 deep conductors at low cost [1]. The time-domain AEM can further be divided into fixed-wing and helicopter EM methods [2]. In each method, a big loop is used to transmit the primary field into the earth and several receiving coils are used to measure the secondary EM signal reflected from the inhomogeneous conductive underground, bringing back the information on the conductivity distribution of the underground media. Right now, the time-domain AEM has been widely used in mineral exploration, oil and gas, hydrogeology, and environment and engineering [3,4]. Because of the high working efficiency and sampling rate in AEM, large amount of data are generally collected. This brings big challenges to EM geophysicists when interpreting the data. Thus, traditionally the AEM data are interpreted by one-dimensional (1D) model [5,6,7]. With the development of computational equipment and technology, three-dimensional (3D) AEM forward modeling has attracted much attention from the geophysical community. Different kinds of 3D modeling methods are developed for AEM simulation. Among them, the finite-difference (FD) method proposed by Yee was first used by Madden and Mackie in magnetotelluric modeling and inversion [8,9]. Newman and Alumbaugh first adopted FD technique in AEM modeling, Liu and Yin used FD method as engine in 3D AEM inversions [10,11]. The advantage of FD method is that it is easy to solve the discontinuous problem of EM field inside the medium, and there are many mature open-source codes that are easy to implement [12]. The integral equation (IE) method is first applied in geophysical EM modeling by Hohmann for 2D anomalous models [13]. It does not need to divide the entire calculation area, but instead only the abnormal bodies, and thereby greatly reduces the number of unknowns. Peltoniemi et al. first applied this method to AEM modeling [14]. The conventional IE method is not suitable for modeling complex geometries as the matrix system is dense. To solve this problem, Cox et al. used quasi-linear approximation to accelerate the IE method in AEM modeling and successfully implemented 3D inversion based on the “moving footprint” concept [15]. Yin et al. used the multi-grid quasi-linear approximation technology to study 3D frequency-domain forward modeling based on IE method and solved the problems of storing the large matrix for Green's function [16]. The finite-element (FE) method emerged in the 1950s was first used to solve fluid mechanics problems. Coggon adopted it to EM modeling, while Wannamaker and Unsworth used this method for MT and controlled-source EM modeling [17,18,19]. Since FE is more flexible and accurate than FD and IE, it quickly became the mainstream technique in geophysical EM modeling [20,21,22]. Especially, the adoption of unstructured tetrahedral grids and the adaptive discretization technology further stimulate researchers to do EM modeling using the FE method. Zhang et al. combined the unstructured vector FE method with adaptive technology to study the 3D AEM forward modeling in the frequency- and time-domain, respectively [23]. Yin et al. showed that the accuracy of AEM numerical simulation of undulating terrain and complex underground structures can be improved by using unstructured FE method [24]. The finite-volume (FV) method is a numerical simulation method mainly used to solve the problems of fluid mechanics and thermodynamics. Haber and Ascher first used it in geophysical EM modeling and solved the problem with magnetic discontinuity. They proved that for electrical transmitting sources, the speed of FV for potential modeling is faster than the staggered grid method [25]. To obtain more accurate AEM modeling, Haber and Schwarzbach proposed a new FV method with multiple OcTree meshes, in which the meshes near the source are refined, so that the modeling accuracy is guaranteed [26]. Haber and Ruthotto further proposed a multi-scale FV method for the Maxwell’s equation in frequency-domain. By solving the local Maxwell’s equations on a coarse grids and then interpolating onto fine grids, they largely reduced the number of unknowns and improved the efficiency [27]. To discretize the earth model more flexibly, Jahandari and Farquharson adopted the unstructured tetrahedral grids into FV for EM modeling, however, the multi-scale and unstructured FV method has not been used in AEM forward modeling yet [28]. It must be pointed out that among these algorithms, the modeling accuracy can only be achieved either by refining the mesh (h-type) or by enhancing the order of the interpolation basis functions (p-type). The spectral element (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.

The Maxwell’s equations in time-domain can be written as where ${\rm E}\left(r,t\right)$ and $H\left(r,t\right)$ are respectively the electric and magnetic field at position r and time t, $B\left(r,t\right)$ is the magnetic induction vector, $D\left(r,t\right)$ is the electric flux density, $J\left(r,t\right)$ is the current density, while $J_s\left(r,t\right)$ is the source current density. Using the constitutive equations $J\left(r,t\right)=\sigma E\left(r,t\right)$, $D\left(r,t\right)=\epsilon E\left(r,t\right)$, $B\left(r,t\right)=\mu H\left(r,t\right)$ and neglecting the displacement current, we have where $\sigma$, $\epsilon$, and $\mu$ represent the conductivity, dielectric constant, and magnetic permeability, respectively. In this paper, we assume that the permittivity of the earth $\epsilon$ = 8.85 × 10⁻¹²F/m and the magnetic permeability μ = 4π × 10⁻⁷ 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_{\mathrm{b}}$ denotes the primary electric field produced by the background conductivity ${\sigma }_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 $z$- direction, then we have $NL=Nx\times \left(Ny+1\right)\times \left(Nz+1\right)+\left(Nx+1\right)\times Ny\times \left(Nz+1\right)+\left(Nx+1\right)\times \left(Ny+1\right)\times Nz$. The vector basis function $\tilde{\Phi }\left(r\right)$ 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 ${\tilde{\Phi }}_j^e$ is the interpolation basis functions, while $E_{sj}^e$ and $E_{bj}^e$ 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 $\left(x,y,z\right)\in \left[x_e, x_{e+1}\right]\times \left[y_e, y_{e+1}\right]\times \left[z_e, z_{e+1}\right]$, while the SE basis functions composed of GLC orthogonal polynomials are defined in $\left[-1, 1\right]$. Thus, we need to convert the coordinates of the eth element in the physical domain $\left(x,y,z\right)\in \left[x_e, x_{e+1}\right]\times \left[y_e, y_{e+1}\right]\times \left[z_e, z_{e+1}\right]$ into the coordinates of the standard element $\left(\xi ,\eta ,\zeta \right)\in \lfloor -1,1\rfloor \times \lfloor -1,1\rfloor \times \lfloor -1,1\rfloor$ 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 $\left|J\right|$ 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 $B{1}^e$ and $B{2}^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 where ${\Phi }_{rst}^{\xi }\left(\xi ,\eta ,\zeta \right)$, ${\Phi }_{rst}^{\eta }\left(\xi ,\eta ,\zeta \right)$, and ${\Phi }_{rst}^{\zeta }\left(\xi ,\eta ,\zeta \right)$ are

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^{\xi }$, $e^{\eta }$, and $e^{\zeta }$ represent the unit vectors along the $\xi$-, $\eta$-, and $\zeta$-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_r^{\left(N\xi -1\right)}\left(\xi \right)$, $h_s^{\left(N\eta \right)}\left(\eta \right)$, and $h_t^{\left(N\zeta \right)}\left(\zeta \right)$ represent the $N\xi -1$, $N\eta$, and $N\zeta$ order of the GLC orthogonal interpolation polynomials that can be expressed as where $k,t=0, 1, 2, \cdots , N$, $T_t\left(\xi \right)=\cos \left(t\cdot \arccos \left(\xi \right)\right)$ is the tth order of Chebyshev polynomial, ${\xi }_k=-\cos \left(k\pi /N\right)$.

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.,

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.