Mainstream numerical methods for 3D time-domain airborne electromagnetic (AEM) modeling, such as the finite-difference (FDTD) or finite-element (FETD) methods, are quite mature. However, these methods have limitations in terms of their ability to handle complex geologic structures and their dependence on quality meshing of the earth model. We have developed a time-domain spectral-element (SETD) method based on the mixed-order spectral-element (SE) approach for space discretization and the backward Euler (BE) approach for time discretization. The mixed-order SE approach can contribute an accurate result by increasing the order of polynomials and suppress spurious solutions. The BE method is an unconditionally stable technique without limitations on time steps. To deal with the rapid variation of the fields close to the AEM transmitting loop, we separate a secondary field from the primary field and simulate the secondary field only, for which the primary field is calculated in advance. To obtain a block diagonal mass matrix and hence minimize the number of nonzero elements in the system of equations to be solved, we apply Gauss-Lobatto-Legendre integral techniques of reduced order. A direct solver is then adopted for the system of equations, which allows for efficient treatment of the multiple AEM sources. To check the accuracy of our SETD algorithm, we compare our results with the semianalytical solution for a layered earth model. Then, we analyze the modeling accuracy and efficiency for different 3D models using deformed physical meshes and compare them against results from 3D FETD codes, to further show the flexibility of SETD for AEM forward modeling.

Three dimensional gravity inversion is an effective way to extract subsurface density distribution from gravity data. Different from the conventional geophysics-based inversions, machine-learning-based inversion is a data-driven method mapping the observed data to a 3D model. We have developed a new machine-learning-based inversion method by establishing a decomposition network (DecNet). Unlike existing machine-learning-based inversion methods, the proposed DecNet method is a mapping from 2D to 2D, which requires less training time and memory space. Instead of learning the density information of each grid point, this network learns the boundary position, vertical center, thickness, and density distribution by 2D-to-2D mapping and reconstructs the 3D model by using these predicted parameters. Furthermore, by using the highly accurate boundary information learned from this network as supplement information, the DecNet method is optimized into a DecNetB method. By comparing the least-squares inversion and U-Net inversion on synthetic and real survey data, the DecNet and DecNetB methods have shown the advantage in dealing with inverse problems for targets with boundaries.

The spectral-element (SE) method, which is based on the Galerkin technique, has been gradually implemented in geophysical electromagnetic (EM) 3D simulation. The accuracy and efficiency of this approach, implemented for deformed hexahedral and regular meshes, has been verified for airborne EM forward modeling. One advantage of the SE method over the conventional finite-element method is that it provides accurate results for earth models that can be adequately parameterized using a coarse mesh. However, realistic models can contain important small-scale conductivity variations or larger features with complicated boundaries. To overcome the limitations imposed by using the same mesh to parameterize the model and for implementing the forward-modeling approach, we have developed an adaptation of the conventional SE method. This is inspired by the ideas behind the element-free Galerkin (EFG) method, in which the conductivity is no longer assumed to be constant within a cell; instead, it is handled via the same kind of numerical integration as in the EFG method. This allows a coarse, regular hexahedral mesh to be used for the forward modeling for complex earth models. After presenting the theory for this new SE approach, we test it for airborne EM modeling of 1D and 3D models to verify its flexibility and accuracy. Finally, we model the Ovoid Zone massive sulfide ore body located at Voisey’s Bay, Labrador, Canada, to illustrate the flexibility and practicality of our approach.

The transient electromagnetic method (TEM) is widely used in the exploration of mineral, petroleum, and geothermal resources due to its sensitivity to low-resistivity bodies, limited site constraints, and strong resistance to interference. In practical applications, the TEM often uses a long wire source instead of an idealized horizontal electric dipole (HED) source as the excitation source. This is due to the complex external conditions and the relatively large distance between the receiving zone and the transmitter source. Compared to the HED, the long wire source can provide a larger excitation current, generating stronger signals to meet the requirements of a higher signal-to-noise ratio or deeper exploration. It also produces longer-duration signals, thereby providing better resolution. Additionally, for the interpretation of TEM data, three-dimensional forward modeling plays a crucial role. However, the mature traditional TEM forward method is based on a simple, sometimes inappropriate model, as it is well established that the induced polarization (IP) effect is widely present in the deep earth, especially in oil and gas reservoirs. The presence of the IP effect results in negative responses in field data that do not conform to the traditional theoretical decay law of TEM, which can significantly impact data processing and inversion results. To address this issue, a TEM forward modeling method considering the IP effect based on the spectral element method (SEM) has been developed in this study. Firstly, starting from the Helmholtz equation satisfied by the time domain electric field, we introduce the Debye model with polarization information into the forward modeling by utilizing the differential form of Ohm’s law. As a result, we derive the boundary value problem for the time domain electric field that considers the induced polarization effect. Using Gauss–Lobatto–Legendre (GLL) polynomials as the basis functions, the SEM is employed to discretize the governing equations at each time step and obtain spectral element discretization equations. Then, temporal discretization equations are derived using the second-order backward Euler formula, and the linear system of equations is solved using the Pardiso direct solver. Finally, the electromagnetic responses at any time channel are obtained via SEM interpolation and numerical integration, thereby achieving three-dimensional TEM forward modeling considering the IP effect. The results indicate that this method can effectively reflect the spatial distribution of polarizable subsurface media. It provides valuable references for studying the polarization parameters of subsurface media and performing a three-dimensional inversion of TEM data considering the induced polarization effect.

Previous works have demonstrated that inverting MT data jointly with gravity data can synergize the high lateral resolution of gravity and the vertical resolution of MT. However, the existing joint stabilizers usually work for structured grids instead of unstructured ones that are more powerful for characterizing complex geology. Here, we utilize Local Pearson Correlation Coefficient (LPCC) for the joint inversion of gravity and MT data based on unstructured grids. We first establish a background mesh by discretizing the research area into virtual rectangular grids and then enhance the structured similarity between density and resistivity via the LPCC. Compared to existing joint constraints, our method is more flexible in solving multiscale joint inversions thanks to the adjustable subdomain size. The synthetic experiments show that the joint stabilizer can recover the subsurface targets at a higher resolution, especially for gravity data, than the standalone inversions. This method is further applied to the joint inversion of gravity and MT data from the Yellowstone area and the inverted density and resistivity models are structurally consistent. Then, based on the inverted subsurface structures and incorporating the existing research, we infer that the two inverted zones with low density and low resistivity correspond to the partially molten rhyolitic and basaltic, respectively. The proposed joint constraint can be further extended to the inversion of other geophysical data.