Search
Close
Search
Advanced Search
Search Menu
AI Discovery Assistant

SUMMARY

Since electromagnetic (EM) fields diffuse more smoothly to greater depth, it physically makes sense to apply fine discretization to model structure at near surface with an increasingly coarser grid both in horizontal and vertical directions as the depth increases for the numerical solution of EM fields. For finite-difference magnetotelluric (MT) forward modelling on regular staggered grids, this can lead to difficulties in discretizing the curl–curl equation governing the EM diffusion in the earth at regions, where the grid resolution changes. In this paper, we propose an efficient algebraic multi-resolution sampling (MRS) method for MT forward modelling. In this method, we do not require the generation of physical subgrids and merely subsample the field components. The coefficient matrix for the subsampled components can be obtained by matrix multiplication based on the initial linear system of equations and field interpolation. To guarantee divergence-free current during the iterative solution process at low frequencies, which is considered crucial for the development of efficient iterative solvers, our forward modelling is based a regularization equation obtained by augmenting the curl–curl equation with an equivalent scaled grad–div operator for electrical fields (explicitly enforcing the divergence-free condition). The correctness of our algebraic MRS algorithm is examined based on a layered model. Its stability and efficiency is exploited by comparing our results with the forward modelling on unsampled staggered grids for the Dublin Test Model 1 (DTM1) and a complex model with realistic topography, indicating a time reduction of about 42–82 per cent.

Electromagnetic theory, Magnetotellurics, Numerical modelling

1 INTRODUCTION

3-D forward modelling of electromagnetic (EM) fields is considered as an important step for EM explorations, through which we can obtain the subsurface structure information from EM survey data (Avdeev 2005; Guo et al. 2020). Due to the complex underground structure, this is generally carried out numerically, by either integral equation methods (IEM, Zhdanov et al. 2006; Kruglyakov & Kuvshinov 2018; Chen et al. 2021), finite-element methods (FEM, Ren et al. 2013; Cai et al. 2014; Grayver & Kolev 2015; Yu et al. 2022), finite-volume methods (FVM, Haber & Ruthotto 2014; Han et al. 2018; Peng et al. 2021; Hu et al. 2022) or finite-difference methods (FDM; Streich 2009; Mittet 2010; Egbert & Kelbert 2012; Kelbert et al. 2014; Li et al. 2018). Due to the simplicity and efficiency, FDM on regular grids finds wide applications in EM forward modelling and inversion. Within this method, the fine discretization of shallow structures in the horizontal directions continues to the bottom of the computational domain (Cherevatova et al. 2018), leading to an unnecessary large linear system of equations.

The well-known fact is that as EM fields diffuse more deeply, they become increasingly smoother and can be represented accurately with coarser grids (Cherevatova et al. 2018). This can be utilized to design a more practical modelling grid, which should be finely discretized in shallow area with increasingly more coarse discretization in deeper regions. The implementation for this on unstructured grids is straightforward (Ren et al. 2013; Jahandari & Farquharson 2015; Tang et al. 2022). For structured grids, previous studies (Haber & Heldmann 2007; Grayver & Bürg 2014; Cherevatova et al. 2018; Gao et al. 2020, 2021; Han et al. 2022) develop multi-resolution (MR) algorithms based on discretizing the partial differential equation with non-conforming grids. Among them, Cherevatova et al. (2018) developed a simpler multi-resolution grid (MRG) scheme, in which they discretize the model into a sequence of subgrids and grid coarsening is only carried out in horizontal directions. For each subgrid, the discretization can be implemented modularly (Cherevatova et al. 2018). However, all the MR algorithms mentioned above have to handle the ‘hanging’ edges at the interfaces of subgrids with different resolutions, posing a challenge for the numerical implementation.

To avoid the difficulty of handling the ‘hanging’ edges during the discretization of the forward modelling operator, Wang et al. (2021) proposed a multilevel forward modelling algorithm to reduce the edge unknowns. Rather than coarsening the physical grid locally, they subsample the edge components sequentially to the pre-defined maximum level, with a increase in the subsampling level resulting in half reduction of the number of edge components (for convenience, we use the name of abstract subgrids to represent the retained components for each level of subsampling). All levels of subsampling are implemented in a pre-determined region where fewer edge components can accurately describe the EM diffusion. In this case, the reduced edge components can be gradually recovered by interpolating the retrained edge components starting from the coarsest abstract grid. A relatively smaller linear system of equation can be formed for the retained edges by the multiplication of interpolation matrices with the linear system of equations on the unsampled grid. However, continuously coarsening one specific local area to the coarsest abstract grid can be problematic due to the increasingly smoother EM diffusion property as penetration depth becomes larger.

Another critical factor affecting the EM forward modelling efficiency based on iterative solvers is that whether the divergence-free current is satisfied during the solution process, particularly at low frequencies (Dong & Egbert 2019). The explicit inclusion of this as a regularization term into the initial curl–curl equation guarantees the right field jumps and can significantly improve the numerical performance of an iterative forward modelling algorithm compared with the classic iterative divergence correction technique (Smith 1996; Farquharson & Miensopust 2011; Tang et al. 2021). It is expected that the combination of MRG into the regularization technique can further improve the forward modelling efficiency.

In this paper, we propose an algebraic MRS approach combined with the regularization technique (Dong & Egbert 2019) for efficient MT forward modelling, which can be considered as an extension of the multilevel algorithm (Wang et al. 2021) and the MRG technique (Cherevatova et al. 2018). A layered geoelectric model from Cherevatova et al. (2018) is used to validate the algebraic MRS algorithm. The numerical performance of the algebraic MRS is investigated by comparing it with the modelling algorithm on regular staggered grids (SG) based on a test model from Miensopust et al. (2013) and a complex model with realistic topography modified from (Ren et al. 2017). The modelling results indicate an impressive improvement against the SG algorithm while not sacrificing accuracy.

2 BASIC THEORY

2.1 Maxwell’s equations

eiωt time-dependent EM fields are assumed. At low frequencies around less than 105 Hz, relative to the conductive current, the displacement current can be negligible. In this case, EM fields (E and H) diffuse according to the simplified Maxwell’s equations, defined by

(1)
(2)

ω represents the angular frequency, μ represents the vacuum permeability, σ represents the conductivity and J represents the external source, which is zero for MT case. By applying the curl operation to eq. (2) and substituting it into eq. (1), the Maxwell’s equations can be reduced to the second partial differential equations for electric fields, given by

(3)

Similarly, the second partial differential equations for magnetic fields can be derived, which is not considered in this paper.

By applying the divergence operator on both sides of eq. (3), we can obtain the current conservation equation, written as

(4)

This indicates that eq. (3) implicitly enforces the conservation equation, which may fail to be satisfied during the numerical solution process of eq. (3) based on iterative solvers. A more discussion can be found in Li et al. (2022). In this case, a process to enforce eq. (4) iteratively is commonly applied with the requirement of solving an additional equation, referred to as the iterative divergence correction technique (Smith 1996; Liu & Yin 2013; Zhou et al. 2021).

2.2 Regularization technique

Alternatively, Dong & Egbert (2019) weight eq. (4) properly and includes it explicitly in eq. (3). We refer to this as the regularization technique since it is similar to the regularization term used in geophysical inversions. To make the solution space consistent, we take the gradient of the weighted version of eq. (4), then include it into eq. (3) to obtain

(5)

where λ is the scaling (weighting) factor. Please note that λ can also be moved before the gradient operator as far as eq. (4) is satisfied. For more complete discussion, please see the references (Dong & Egbert 2019; Li et al. 2020, 2022). The auxiliary magnetic fields are determined by

(6)

3 NUMERICAL IMPLEMENTATION

3.1 Finite-difference discretization

With proper boundary conditions, we can discretize eq. (5) on an SG (Yee 1966), shown in Fig. 1 with electric fields defined on cell edges and magnetic fields defined on cell faces. Then we can obtain the discrete linear system of equations in matrix form, written as

(7)

where e is the vector of partitioned electric fields, and b is the vector containing applied boundary conditions. A is the coefficient matrix, given by

(8)

the curl operator C maps cell edges to cell faces, while the adjoint of the curl operator C† maps cell faces back to cell edges, the gradient operator G maps from nodes to edges, the divergence operator D maps from edges to nodes, Λ is a diagonal matrix filled with values of scaling factor and |$\bar{\sigma }$| is the volume average conductivity of four cells surrounding each edge. For more details on the formulation of these operators, please refer to the work (Egbert & Kelbert 2012; Dong & Egbert 2019; Li et al. 2020).

The discretization of standard Yee grid.
Figure 1.

The discretization of standard Yee grid.

Open in new tabDownload slide

3.2 Implementation of algebraic multi-resolution sampling technique

For the algebraic MRS technique developed in this paper, we don’t need to generate a series of physical subgrids as in the work (Cherevatova et al. 2018). Once we generate eq. (7) on the fine grid, we use matrix multiplication to carry out the subsampling, demonstrated subsequently. If we assume that the discrete e0 in eq. (7) on the fine grid with N0 edges, can be interpolated from a subensemble of edges e1 on the first-level abstract subgrid. Then we can define e0 as the multiplication of an interpolation matrix and the subensemble vector, giving

(9)

Similarly, if we have a consecutive sequence of subsampled edge sets e1, e2,..., eN for N levels of subsampling with respective sparse interpolation matrices S1, S2,..., SN, we can form the interpolation chain as

(10)

For convenience, l can be considered as the subsampling level.

If we assume the numbers of edges included in e1, e2,..., eN are N1, N2,..., NN, respectively. Using the chain rule, we can get

(11)

S represents the total interpolation matrix with a dimension of N0 × NN. Substituting eq. (11) into eq. (7), we obtain

(12)

where AS has a dimension of N0 × NN, which is not a square matrix. In this case, if we want to solve eq. (12) for eN, we need to delete the row corresponding to the edges deleted during the subsampling procedure. This can be trivially realized by multiplying eq. (12) on both sides by an NN × N0 matrix L with only one entry of 1 at the column corresponding to one retained edge and 0 at the rest columns for each row, resulting in the following equation

(13)

Solving eq. (13) for eN, the whole original edge components e0 can be obtained using eq. (11).

So far, we have not mentioned how to implement the sparse interpolation matrices S1, S2,..., SN, which will be addressed subsequently. For convenience, we will use x-edges to demonstrate this, as shown in Fig. 2. Unlike the work (Wang et al. 2021) using a uniform subsampling at one fixed region, the subsampling at different regions is applied differently, for instance, zero-level subsampling at the near-surface, followed by increasingly higher-level subsampling as the depth increases.

An example of the algebraic MRS algorithm for electric fields ex. Each plane represents a subsampling level. The red dashed lines represent the removed edges during the first-level subsampling, the blue dashed lines represent the removed edges during the second-level subsampling and solid lines denote the retained edges.
Figure 2.

An example of the algebraic MRS algorithm for electric fields ex. Each plane represents a subsampling level. The red dashed lines represent the removed edges during the first-level subsampling, the blue dashed lines represent the removed edges during the second-level subsampling and solid lines denote the retained edges.

Open in new tabDownload slide

For easy demonstration, we use planes of single-row x-edges to show how Sl is calculated, as shown in Fig. 2 for a two-level subsampling scheme. As mentioned before, we form the abstract subgrids as in the figure just for demonstration purposes, which is not needed at all in the numerical implementation. For this purpose, we increase the level of subsampling by one as one additional layer is included. In the 0th subregion, no subsampling is applied. For the first-level subsampling, as indicated in Fig. 2, the deleted edges |$e_{x_{1,2i}}$| with i = 1,..., Ny/2 on plane 2 (assume the total x-edge number Ny along the y-direction on plane 2 is even) with the first number in the subscript for the level of sampling and the second number for the position of the electric field on the plane 2, can be interpolated from the retained edges |${e}_{x_{1,2i-1}}$| and |${e}_{x_{1,2i+1}}$| as

(14)
$$\begin{eqnarray} \text{where}\left\lbrace {\begin{array}{*{10}c}{c}_{2i-1}=\frac{{{d}_{y}}_{_{2i}}}{{{d}_{y}}_{_{2i-1}}+{{d}_{y}}_{_{2i}}} \\ {{c}_{2i+1}}=\frac{{{d}_{y}}_{_{2i-1}}}{{{d}_{y}}_{_{2i-1}}+{{d}_{y}}_{_{2i}}} \\ \end{array}} ,\right.\ \ \end{eqnarray}$$

|${{d}_{y}}_{_{2i}}$| is the edge length for the 2ith cell on the initial grid along the y-direction. Thus, we can obtain the entries of c2i − 1 and c2i + 1 at the row corresponding to the position of |${{e}_{{{x}_{1,2i}}}}$| in e0 and at the columns corresponding to the position of |${{e}_{{{x}_{1,2i-1}}}}$| and |${{e}_{{{x}_{1,2i+1}}}}$| in e1 with the rest of the columns filled with zeros for S1. For rows related to the retained edges (including edges in the air-earth region), we fill them with zeros except for setting ones for the diagonals. However, in our implementation, we use the same trick as in Wang et al. (2021), within which we first generate a N0 × N0 identity matrix S0 as shown in Fig. 3(a), then we add c2i − 1 and c2i + 1 at the row and columns corresponding to the indices of |${{e}_{{{x}_{1,2i}}}}$|⁠, |${{e}_{{{x}_{1,2i-1}}}}$| and |${{e}_{{{x}_{1,2i+1}}}}$| all in e0, respectively, with the rest entries of this row set to zeros, to generate the N0 × N0 version of S1 as shown in Fig. 3(b). Finally, we delete the zero columns to obtain N0 × N1S1.

The distribution of matrix elements based on a small-sized model (8 × 8 × 3). (a) represents of the identity matrix S0, and (b) represents the interpolation matrix S1 having the same size of S0 for the first-level subsampling. Red represents the interpolation coefficients and blue dotted lines represent ones with white for zeros.
Figure 3.

The distribution of matrix elements based on a small-sized model (8 × 8 × 3). (a) represents of the identity matrix S0, and (b) represents the interpolation matrix S1 having the same size of S0 for the first-level subsampling. Red represents the interpolation coefficients and blue dotted lines represent ones with white for zeros.

Open in new tabDownload slide

Similarly, for the second-level subsampling, the deleted edges |${{e}_{{{x}_{2,4i\text{-}1}}}}$| on plane 3 can be interpolated from the retained edges |${{e}_{{{x}_{2,4i-3}}}}$| and |${{e}_{{{x}_{2,4i+1}}}}\, (i=1,2,...,\frac{{{N}_{y}}}{{{2}^{2}}}$|⁠) as

(15)
$$\begin{eqnarray} \text{where}\left\lbrace {\begin{array}{*{10}c}{c}_{4i-3}=\frac{{{d}_{y}}_{_{4i\text{-}1}}\text{+}{{d}_{y}}_{_{4i}}}{{{d}_{y}}_{_{4i-3}}+{{d}_{y}}_{_{4i\text{-}2}}\text{+}{{d}_{y}}_{_{4i\text{-}1}}\text{+}{{d}_{y}}_{_{4i}}} \\ {{c}_{4i+1}}=\frac{{{d}_{y}}_{_{4i\text{-3}}}\text{+}{{d}_{y}}_{_{4i\text{-}2}}}{{{d}_{y}}_{_{4i-3}}+{{d}_{y}}_{_{4i\text{-}2}}\text{+}{{d}_{y}}_{_{4i\text{-}1}}\text{+}{{d}_{y}}_{_{4i}}} \\ \end{array}}.\right. \end{eqnarray}$$

We use a similar way as in the first-level subsampling to assemble S2. The only difference is that now the rows in S2 are decided by edge positions in e1, and its columns are decided by edge positions in e2. Repeat this process to form a series of Sl, until the preset highest level is reached. In implementation, we also use the similar way as used for the formulation of |${{N}_{0}}\times {{N}_{0}} \, {{\mathbf {S}}_{1}}$| to obtain |${{N}_{0}}\times {{N}_{0}} \, {{\mathbf {S}}_{2}}$|⁠. Then, we delete the zero columns and the rows corresponding to the deleted components to generate |${{N}_{1}}\times {{N}_{2}} \, {{\mathbf {S_{2}}}}$|⁠.

For the general Nth-level subsampling, we can interpolate the deleted edges |${{e}_{{{x}_{N,{{2}^{N-1}}\cdot (2i-1)+1}}}}$| from the retained edges |${{e}_{{{x}_{N,{{2}^{N}}\cdot (i-1)+1}}}}$| and |${{e}_{{{x}_{N,{{2}^{N}}\cdot i+1}}}}\, (i=1,2,...,\frac{{{N}_{y}}}{{{2}^{N}}})$| as

(16)
$$\begin{eqnarray} \text{where}\left\lbrace {\begin{array}{*{10}c}{c}_{{{2}^{N}}\cdot (i-1)+1}=\frac{\sum \limits _{j=({{2}^{N}}-{{2}^{N-1}})\cdot i}^{{{2}^{N}}\cdot i}{{{d}_{{{y}_{j}}}}}}{\sum \limits _{j=1}^{{{2}^{N}}\cdot i}{{{d}_{{{y}_{j}}}}}}\\ {{c}_{{{2}^{N}}\cdot i+1}}=\frac{\sum \limits _{j=1}^{({{2}^{N}}-{{2}^{N-1}})\cdot i}{{{d}_{{{y}_{j}}}}}}{\sum \limits _{j=1}^{{{2}^{N}}\cdot i}{{{d}_{{{y}_{j}}}}}}\\ \end{array}}.\right. \end{eqnarray}$$

Similarly, we can repeat this interpolation process for y-edges.

For the z-edges ez, subsampling is implemented in both the x- and y-directions (Wang et al. 2021) to guarantee the accuracy. In this case, we use a general grid with 8 × 8 cells in the horizontal plane for demonstration purpose. For simplicity, we use one plane to show the first-level subsampling as shown in Fig. 4. As indicated in figure, one component deleted is surrounded and interpolated by four retained components. For instance, |${{e}_{{{z}_{1,11}}}}$| can be interpolated from |${{e}_{{{z}_{1,2}}}}\text{, }{{e}_{{{z}_{1,10}}}},\text{ }{{e}_{{{z}_{1,12}}}}\, \mathrm{ and}\,\text{ }{{e}_{{{z}_{1,20}}}}$| as

(17)
$$\begin{eqnarray} \text{where}\left\lbrace {\begin{array}{*{10}c}{c}_{10}=\frac{{{d}_{{{y}_{2}}}}}{{{d}_{{{y}_{1}}}}+{{d}_{{{y}_{2}}}}} \\ {{c}_{12}}=\frac{{{d}_{{{y}_{1}}}}}{{{d}_{{{y}_{1}}}}+{{d}_{{{y}_{2}}}}} \\ {{c}_{2}}=\frac{{{d}_{{{x}_{2}}}}}{{{d}_{{{x}_{1}}}}+{{d}_{{{x}_{2}}}}} \\ {{c}_{20}}=\frac{{{d}_{{{x}_{1}}}}}{{{d}_{{{x}_{1}}}}+{{d}_{{{x}_{2}}}}} \\ \end{array}} \right.,{\begin{array}{*{10}c}{w}_{1}=\frac{{{d}_{{{x}_{1}}}}+{{d}_{{{x}_{2}}}}}{{{d}_{{{x}_{1}}}}+{{d}_{{{x}_{2}}}}+{{d}_{{{y}_{1}}}}+{{d}_{{{y}_{2}}}}} \\ {{w}_{2}}=\frac{{{d}_{{{y}_{1}}}}+{{d}_{{{y}_{2}}}}}{{{d}_{{{x}_{1}}}}+{{d}_{{{x}_{2}}}}+{{d}_{{{y}_{1}}}}+{{d}_{{{y}_{2}}}}} \\ \end{array}}. \end{eqnarray}$$

|${{d}_{x}}_{_{i}}$| is the edge length on the initial grid for the ith element along the x-direction. We can trivially generalize this procedure to the Nth-level subsampling for ez. In the Nth subregion, the deleted edge |${{e}_{{{z}_{N,{{2}^{N-1}}\left( 2i+{{N}_{y}} \right)+1}}}}$| during the Nth-level subsampling can be interpolated from |${{e}_{{{z}_{N,{{2}^{N-1}}\left( 2i-1 \right)+1}}}}$|⁠, |${{e}_{{{z}_{N,{{2}^{N-1}}\left( 2i+{{N}_{y}}-1 \right)+1}}}}$|⁠, |${{e}_{{{z}_{N,{{2}^{N-1}}\left( 2i+{{N}_{y}}+1 \right)+1}}}}\text{ and }{{e}_{{{z}_{N,{{2}^{N-1}}\left( 2i+2{{N}_{y}}+1 \right)+1}}}}$| as

(18)
$$\begin{eqnarray} \text{where}\left\lbrace {\begin{array}{*{10}c}{c}_{{{2}^{N-1}}\left( 2i+{{N}_{y}}-1 \right)+1}=\frac{\sum \limits _{j={{2}^{(N-1)}+1}}^{{{2}^{N}}}{{{d}_{{{y}_{j}}}}}}{\sum \limits _{j=1}^{{{2}^{N}}}{{{d}_{{{y}_{j}}}}}} \\ {{c}_{{{2}^{N-1}}\left( 2i+{{N}_{y}}+1 \right)+1}}=\frac{\sum \limits _{j=1}^{{{2}^{(N-1)}}}{{{d}_{{{y}_{j}}}}}}{\sum \limits _{j=1}^{{{2}^{N}}}{{{d}_{{{y}_{j}}}}}} \\ {{c}_{{{2}^{N-1}}\left( 2i-1 \right)+1}}=\frac{\sum \limits _{k={{2}^{(N-1)}+1}}^{{{2}^{N}}}{{{d}_{{{x}_{k}}}}}}{\sum \limits _{k=1}^{{{2}^{N}}}{{{d}_{{{x}_{k}}}}}} \\ {{c}_{{{2}^{N-1}}\left( 2i+2{{N}_{y}}+1 \right)+1}}=\frac{\sum \limits _{k=1}^{{{2}^{(N-1)}}}{{{d}_{{{x}_{k}}}}}}{\sum \limits _{k=1}^{{{2}^{N}}}{{{d}_{{{x}_{k}}}}}} \\ \end{array}} \right.,{\begin{array}{*{10}c}{w}_{1}=\frac{\sum \limits _{k=1}^{{{2}^{N}}}{{{d}_{{{x}_{k}}}}}}{\sum \limits _{k=1}^{{{2}^{N}}}{{{d}_{{{x}_{k}}}}}+\sum \limits _{j=1}^{{{2}^{N}}}{{{d}_{{{y}_{j}}}}}} \\ {{w}_{2}}=\frac{\sum \limits _{j=1}^{{{2}^{N}}}{{{d}_{{{y}_{j}}}}}}{\sum \limits _{k=1}^{{{2}^{N}}}{{{d}_{{{x}_{k}}}}}+\sum \limits _{j=1}^{{{2}^{N}}}{{{d}_{{{y}_{j}}}}}} \\ \end{array}}. \end{eqnarray}$$

The positions of the interpolation coefficients in Sl are determined similarly as for ex, mentioned previously.

An example of the algebraic MRS algorithm for electric fields ez. The blue triangles and grey polka dots are the reduced and retained vertical electric fields after the first-level subsampling on the plane, respectively.
Figure 4.

An example of the algebraic MRS algorithm for electric fields ez. The blue triangles and grey polka dots are the reduced and retained vertical electric fields after the first-level subsampling on the plane, respectively.

Open in new tabDownload slide

All of these matrices of eq. (13) are combined into a single matrix LAS during the solution process. Since it is not mathematically obvious that LAS is also a sparse matrix, we examined this by plotting the sparsity pattern of LAS. Considering that the resulting matrix pattern is independent of the complexity of a model, based on a three-layered earth model (see Section 4.1) discretized with a 16 × 16 × 42 grid (plus 12 air layers), we plot the resulting matrix for the sampling scheme of [S,N] = [1; 8; 0; 16; 1; 10; 2; 20] in Fig. 5, indicating the sparsity of the resulting matrix.

Matrix sparsity of LAS in the sampling scheme of [S,N] = [1; 8; 0; 16; 1; 10; 2; 20] for the three-layered earth model.
Figure 5.

Matrix sparsity of LAS in the sampling scheme of [S,N] = [1; 8; 0; 16; 1; 10; 2; 20] for the three-layered earth model.

Open in new tabDownload slide

For convenience, we use a similar way as in work (Cherevatova et al. 2018) to represent the abstract subsampling scheme. For example, if we only consider a two-level multi-resolution subsampling in Fig. 2, then the subsampling scheme can be represented by the vector [S, Nz] = ⌊0, 1; 1, 1; 2, 1⌋, with S indicating the sampling level and Nz for the number of layers in the z-direction for this subsampling level. Each pair of numbers separated by a semicolon represent a particular choice of [S, Nz] for one abstract subregion, organized from the top to the bottom of the computational domain. As the subsampling level increases by one, the size of the resulting components in one plane is reduced by half. For this example, the computational domain in Fig. 2 is divided into three abstract subregions. The first abstract subregion [S, Nz] = [0, 1] includes one layer of 1 × 8 horizontal cells. The second abstract subregion [S, Nz] = [1, 1], indicating the first-level subsampling, includes one layer of 1 × 4 horizontal cells and the third abstract subregion [S, Nz] = [2, 1] sampled twice contains one layer of 1 × 2 horizontal cells. No subsampling is applied in the z-direction.

In general, any subsampling schemes can be used as far as the subsampling does not bring obvious accuracy deterioration. However, there are no rigorous ways to determine an optimal subsampling scheme. The physical fact is that the EM fields are generally not very sensitive to geoelectric structures deeper than one skin depth and coarser grids can be used to describe the diffusion of EM fields accurately, as the fields go deeper. In this paper, this has been utilized to design the subsampling scheme. The subsampling starts from one skin depth which physically makes sense, then increases the spacing of two neighbouring levels as the subsampling moves to deeper regions. This particular choice leads to different subsampling schemes for different periods and it is a natural choice to gradually increase the subsampling level as the subsampling goes deeper.

4 NUMERICAL EXAMPLES

The accuracy of the algebraic MRS technique is examined based on a layered geoelectric model. We then illustrate the effectiveness of the algebraic MRS technique based on DTM1 and a realistic terrain model. In this paper, BiCGstab with SSOR pre-conditioner is used to solve the linear system of equations for both SG and MRS. All tests coded in MATLAB are performed on an AMD Ryzen Threadripper 2.9 GHz computer in serial, with the preset relative residual tolerance of 10−10 and maximum iterations of 3000.

4.1 Three-layered earth model

In this section, the accuracy of the algebraic MRS technique developed in this paper is validated by a layered geoelectric model from Cherevatova et al. (2018), shown in Fig. 6. The computational domain is discretized into 16 × 16 × 42 cells (plus 12 air layers), with a uniform discretization of 500 m in the horizontal direction and non-uniform discretization in the vertical direction. This leads to a total of 45 526 discrete edge components. Periods 1, 10, 100, 1000 and 10000 s are considered for this model.

Vertical cross-section for the layered model from Cherevatova et al. (2018).
Figure 6.

Vertical cross-section for the layered model from Cherevatova et al. (2018).

Open in new tabDownload slide

For this model, we use the same mutiresolution scheme in Cherevatova et al. (2018). In this scheme, the model is subsampled with [S, N] = [1, 8; 0, 16; 1, 10; 2, 20], which divides the computational domain into four subregions. No subsampling is applied in the region close to the earth surface, containing 4 air layers and 12 earth layers. We apply the one-level subsampling in the immediate neighbouring abstract subregions on both sides of the unsampled region. In the upper abstract subregion consisting of eight air layers, the subsampling leads to 8 × 8 cells in the horizontal planes. Similarly, below the unsampled region, 10 layers are subsampled horizontally. The two-level subsampling is applied in the bottom abstract subregions of 20 layers, resulting in 4 × 4 cells in the horizontal planes. Due to the use of the algebraic MRS scheme, a total of 14 784 edges are removed.

We also have calculated the condition number for A and LAS in eq. (13) for the three-layered earth model at the considered periods. As shown in Table 1, the resulting condition numbers for A and LAS are comparable for all considered periods.

Table 1.
Open in new tab

The condition number of A and LAS of eq. (13) for the three-layered earth model.

Period
1 s10 s100 s1000 s10000 s
A1.65 × 1051.13 × 1051.14 × 1051.14 × 1051.14 × 105
LAS1.61 × 1051.01 × 1051.02 × 1051.02 × 1051.02 × 105
Period
1 s10 s100 s1000 s10000 s
A1.65 × 1051.13 × 1051.14 × 1051.14 × 1051.14 × 105
LAS1.61 × 1051.01 × 1051.02 × 1051.02 × 1051.02 × 105
Table 1.
Open in new tab

The condition number of A and LAS of eq. (13) for the three-layered earth model.

Period
1 s10 s100 s1000 s10000 s
A1.65 × 1051.13 × 1051.14 × 1051.14 × 1051.14 × 105
LAS1.61 × 1051.01 × 1051.02 × 1051.02 × 1051.02 × 105
Period
1 s10 s100 s1000 s10000 s
A1.65 × 1051.13 × 1051.14 × 1051.14 × 1051.14 × 105
LAS1.61 × 1051.01 × 1051.02 × 1051.02 × 1051.02 × 105

The reference solution is calculated using the analytical method at logarithmically spaced periods from 1 to 10000 s, based on which the accuracy of our algorithm is examined. As shown in Fig. 7, the results from our algebraic MRS method and SG are both consistent with the analytical solution. MRS and SG algorithms have almost identical relative errors, indicating the subsampling does not deteriorate the modelling accuracy obviously for this case.

Comparison of (a) apparent resistivity and (b) phase based a layered model from the algebraic MRS technique and the SG scheme with the analytical solution over logarithmically distributed periods from 1 to 10000 s.
Figure 7.

Comparison of (a) apparent resistivity and (b) phase based a layered model from the algebraic MRS technique and the SG scheme with the analytical solution over logarithmically distributed periods from 1 to 10000 s.

Open in new tabDownload slide

4.2 Dublin Test Model 1

DTM1 (Miensopust et al. 2013), is designed to investigate the performance of our algorithm. As indicated in Fig. 8, the test model includes three resistivity anomalies buried in a background of 100 Ω · m, whose details are listed in Table 2. To provide accurate EM modelling, we utilize a model-dependent discretization as in Miensopust et al. (2013). The computational domain is discretized into 160 × 160 × 80 cells (excluding 12 the air layer) with the smallest cell sizes of |$250\times 250\times 250\, {{\text{m}}^{3}}$| (≤40s) and |$2500\times 2500\times 2500\, {{\text{m}}^{3}}$| (>40s). This leads to a total of 7.18 million discrete edge components. Periods 0.1, 1, 10,100 and 1000 s are considered for this model.

Sketch of DTM1 from Miensopust et al. (2013).
Figure 8.

Sketch of DTM1 from Miensopust et al. (2013).

Open in new tabDownload slide
Table 2.
Open in new tab

The details of anomalous bodies for DTM1.

x (km)y (km)z (km)Resistivity (⁠|$\Omega \, \cdot$|m)
Body1−20 to 20−2.5 to 2.55 to 2010
Body2−15 to 0−2.5 to 22.520 to 251
Body30 to 15−22.5 to 2.520 to 5010 000
x (km)y (km)z (km)Resistivity (⁠|$\Omega \, \cdot$|m)
Body1−20 to 20−2.5 to 2.55 to 2010
Body2−15 to 0−2.5 to 22.520 to 251
Body30 to 15−22.5 to 2.520 to 5010 000
Table 2.
Open in new tab

The details of anomalous bodies for DTM1.

x (km)y (km)z (km)Resistivity (⁠|$\Omega \, \cdot$|m)
Body1−20 to 20−2.5 to 2.55 to 2010
Body2−15 to 0−2.5 to 22.520 to 251
Body30 to 15−22.5 to 2.520 to 5010 000
x (km)y (km)z (km)Resistivity (⁠|$\Omega \, \cdot$|m)
Body1−20 to 20−2.5 to 2.55 to 2010
Body2−15 to 0−2.5 to 22.520 to 251
Body30 to 15−22.5 to 2.520 to 5010 000

The abstract subsampling for this model at 100 s is defined by |$[{S,N}]=[2,3;1,3;0,\, 19;1,6;\, 2,6;3,5;4,9;5,41]$| leading to eight abstract subregions with the highest subsampling level of 5. The first and second abstract subregions apply the two and one-level subsampling, respectively, each containing three air layers. Their sizes are 40 × 40 × 3 cells and 80 × 80 × 3 cells, respectively. The third subregion contains the earth surface with 6 air layers and 13 earth layers. No subsampling is applied in this region discretized with 160 × 160 × 19 cells. The sampling level for subsequent subregions is increased by one, as they move to each deeper region, which leads to respective abstract grid sizes of 80 × 80 × 6, 40 × 40 × 6, 20 × 20 × 5, 10 × 10 × 9 and 5 × 5 × 41 cells, respectively. Relative to SG, the algebraic MRS scheme reduces the edges by 3.82 million (almost half the original size) for this period. Similarly, we use [S, N] = [2, 3; 1, 3; 0, 13; 1, 3; 2, 3; 3, 2; 4, 5; 5, 60], [S, N] = [2, 3; 1, 3; 0, 21; 1, 6; 2, 4; 3, 4; 4, 5; 5, 46], [S, N] = [2, 3; 1, 3; 0, 35; 1, 7; 2, 8; 3, 6; 4, 10; 5, 20] and [S, N] = [2, 3; 1, 3; 0, 37; 1, 10; 2, 9; 3, 7; 4, 13; 5, 10] for 0.1, 1, 10 and 1000 s, respectively. Their respective abstract grid size for each subregion can be trivially obtained by repeating the procedure as used for 100 s.

Then, the numerical performance of the algebraic MRS technique is compared with that of the forward modelling on SG and MRG (Cherevatova et al. 2018) in terms of the total number of discrete edges and computation time. As shown in Table 3, the algebraic MRS technique significantly reduces the total number of discrete edges relative to SG, resulting in much less non-zero elements as indicated in Table 4. Consequently, this leads to less computation time for our method. For instance, at 1 s, the computation time for our algebraic MRS method is nearly 25 per cent of that for SG (156 s vs 613 s). Interestingly, for this period, the reduction in run time is much larger than that in edge components (75 per cent vs 52 per cent), possibly due to the nonlinear relationship between the degrees of freedom of the resulting linear system and the computation time (Carriere & Jeandel 1991). This nonlinearity accounts for the larger reduction in time for the MRS method than in DoFs. For MRG and our algebraic MRS technique, their accuracy is comparable and our algorithm performs better in terms of computation time at most periods, possibly due to the fact that the hanging edges in MRG result in worse conditioning for the resulting system of equations. In Fig. 9, we compare the convergence history for SG and MRS for the Dublin Test Model 1 (DTM1). As we noted, MRS generally takes fewer iterations than SG for all considered periods.

Convergence process of SG and MRS at periods from 0.1 to 1000 s for DTM1.
Figure 9.

Convergence process of SG and MRS at periods from 0.1 to 1000 s for DTM1.

Open in new tabDownload slide
Table 3.
Open in new tab

Numerical performance of SG scheme, MRG and algebraic MRS technique in terms of run time and the number of sampled edges for DTM1.

Period
0.1 s1 s10 s100 s1000 s
MethodTime (s) and no. sampled edges (million)
SG1737.186137.184787.182007.182547.18
MRS652.881563.432764.34883.351444.55
MRG1181.221771.902183.021351.761463.24
MRS/SG (per cent)38402548586044475763
MRG/SG (per cent)68172926464268255745
Period
0.1 s1 s10 s100 s1000 s
MethodTime (s) and no. sampled edges (million)
SG1737.186137.184787.182007.182547.18
MRS652.881563.432764.34883.351444.55
MRG1181.221771.902183.021351.761463.24
MRS/SG (per cent)38402548586044475763
MRG/SG (per cent)68172926464268255745
Table 3.
Open in new tab

Numerical performance of SG scheme, MRG and algebraic MRS technique in terms of run time and the number of sampled edges for DTM1.

Period
0.1 s1 s10 s100 s1000 s
MethodTime (s) and no. sampled edges (million)
SG1737.186137.184787.182007.182547.18
MRS652.881563.432764.34883.351444.55
MRG1181.221771.902183.021351.761463.24
MRS/SG (per cent)38402548586044475763
MRG/SG (per cent)68172926464268255745
Period
0.1 s1 s10 s100 s1000 s
MethodTime (s) and no. sampled edges (million)
SG1737.186137.184787.182007.182547.18
MRS652.881563.432764.34883.351444.55
MRG1181.221771.902183.021351.761463.24
MRS/SG (per cent)38402548586044475763
MRG/SG (per cent)68172926464268255745
Table 4.
Open in new tab

The non-zero elements in the coefficient matrix of SG and MRS for DTM1.

Period
0.1 s1 s10 s100 s1000 s
A8.47 × 1078.47 × 1078.47 × 1078.47 × 1078.47 × 107
LAS4.73 × 1075.27 × 1076.12 × 1075.19 × 1076.31 × 107
Period
0.1 s1 s10 s100 s1000 s
A8.47 × 1078.47 × 1078.47 × 1078.47 × 1078.47 × 107
LAS4.73 × 1075.27 × 1076.12 × 1075.19 × 1076.31 × 107
Table 4.
Open in new tab

The non-zero elements in the coefficient matrix of SG and MRS for DTM1.

Period
0.1 s1 s10 s100 s1000 s
A8.47 × 1078.47 × 1078.47 × 1078.47 × 1078.47 × 107
LAS4.73 × 1075.27 × 1076.12 × 1075.19 × 1076.31 × 107
Period
0.1 s1 s10 s100 s1000 s
A8.47 × 1078.47 × 1078.47 × 1078.47 × 1078.47 × 107
LAS4.73 × 1075.27 × 1076.12 × 1075.19 × 1076.31 × 107

For the point that the reduction in time is more significant at high frequencies, this is mainly caused by the use of the subsampling scheme. As mentioned before, the subsampling starts at one skin depth, which leads to different start points for different frequencies. If the same initial discretization grids are used for different frequencies (the same initial grids for a range of frequencies is used in this paper), this results in different sampled edges for different frequencies with smaller sampled edges for higher frequencies, accounting for a more significant reduction in time at high frequencies.

Finally, we examine the accuracy of our algorithm again by comparing our solution with the public result from Nam et al. (2007) at 100 s along one survey line (x=0 km) in Fig. 10. For this purpose, we also calculate the solution for SG. As indicated in the figure, both the algebraic MRS method and SG produce results in good agreement with the reference solution. The relative difference for apparent resistivity is below 3.9 per cent, and the absolute phase difference is below 0.8°.

Comparison of solution from the algebraic MRS technique, MRG technique and the SG scheme with the reference solution from Nam et al. (2007) for DTM1. (a) Apparent resistivities for XY polarization, (b) apparent resistivities for YX polarization, (d) phases for XY polarization and (e) phases for YX polarization. (c) and (f) are the relative difference for apparent resistivity and the phase difference against the reference solution, respectively.
Figure 10.

Comparison of solution from the algebraic MRS technique, MRG technique and the SG scheme with the reference solution from Nam et al. (2007) for DTM1. (a) Apparent resistivities for XY polarization, (b) apparent resistivities for YX polarization, (d) phases for XY polarization and (e) phases for YX polarization. (c) and (f) are the relative difference for apparent resistivity and the phase difference against the reference solution, respectively.

Open in new tabDownload slide

Since the choice of subsampling scheme used in this paper results in the dependence of the subsampling level, the number of abstract subregions and the number of layers contained in each abstract subregion, we examined the effect of different choices of the number of abstract subregions (the maximum of subregions is limited by the initial grid of the model) on numerical performance for the DTM1 at 1 and 100 s. For this purpose, we design two different subsamplings, the subsampling defined by [S, N] = [2, 3; 1, 3; 0, 21; 1, 19; 2, 46] (referred as MRS-1) for 1 s and [S, N] = [2, 3; 1, 3; 0, 19; 1, 26; 2, 41] for 100 s and the five-level subsampling (referred as MRS-2) as used in the previous tests. MRS-2 can be considered as a result of increasing the number of abstract subregions from MRS-1. As listed in Table 5, the result shows that while adding more subregions improves the modelling efficiency in proportion to the reduction of edges, it does not obviously affect the accuracy as shown in Fig. 11.

Comparison of solution from the different MRS schemes for DTM1. (a) Apparent resistivities for XY polarization, (b) apparent resistivities for YX polarization, (d) phases for XY polarization and (e) phases for YX polarization. (c) and (f) are the relative difference for apparent resistivity and the phase difference against the reference solution, respectively.
Figure 11.

Comparison of solution from the different MRS schemes for DTM1. (a) Apparent resistivities for XY polarization, (b) apparent resistivities for YX polarization, (d) phases for XY polarization and (e) phases for YX polarization. (c) and (f) are the relative difference for apparent resistivity and the phase difference against the reference solution, respectively.

Open in new tabDownload slide
Table 5.
Open in new tab

Numerical performance of different sampling schemes in terms of run time and the number of sampled edges for DTM1.

Period
1 s100 s
Time (s) and no. sampled edges (million)
MRS-12044.231184.24
MRS-21563.43883.35
MRS-2/MRS-1 (per cent)76817579
Period
1 s100 s
Time (s) and no. sampled edges (million)
MRS-12044.231184.24
MRS-21563.43883.35
MRS-2/MRS-1 (per cent)76817579
Table 5.
Open in new tab

Numerical performance of different sampling schemes in terms of run time and the number of sampled edges for DTM1.

Period
1 s100 s
Time (s) and no. sampled edges (million)
MRS-12044.231184.24
MRS-21563.43883.35
MRS-2/MRS-1 (per cent)76817579
Period
1 s100 s
Time (s) and no. sampled edges (million)
MRS-12044.231184.24
MRS-21563.43883.35
MRS-2/MRS-1 (per cent)76817579

4.3 Realistic terrain model

In this part, we consider a realistic terrain model modified from Ren et al. (2017), called the realistic model for convenience. This model covers an area from −24 to 24 km in the horizontal directions with the highest altitude of 1163 m as shown in Fig. 12. The realistic model has an arbitrary resistivity variation below the surface, as shown in the Fig. 13. To guarantee the modelling accuracy, we extend the computation domain with an additional distance of 80 km in each direction. The computation domain is discretized into 224 × 224 × 112 cells (excluding the 12 air layers) with the smallest cell size of |$360\times 360\times 80\, {{\text{m}}^{3}}$| (≤1s) and |$360\times 360\times 200\, {{\text{m}}^{3}}\, (\gt 1\text{ s})$|⁠. This discretization results in a total of 18.88 million discrete edge components.

Illustration of the terrain for the realistic model modified from Ren et al. (2017).
Figure 12.

Illustration of the terrain for the realistic model modified from Ren et al. (2017).

Open in new tabDownload slide
3-D resistivity distribution below the surface for the realistic model. (a) y–z slice view and (b) x–z slice view.
Figure 13.

3-D resistivity distribution below the surface for the realistic model. (a) y–z slice view and (b) x–z slice view.

Open in new tabDownload slide

For the test at 2 s, the subsampling scheme is defined by the subsampling vector [S, N] = [2, 3; 1, 3;0, 30; 1, 6; 2, 6; 3, 5; 4, 7; 5, 64] with a total of eight abstract subregions and the highest subsampling level of five. The third subregion, a unsampled region, includes 6 air layers and 24 earth layers of 224 × 224 horizontal cells. The first and second abstract subregions above the third subregion have sizes of 56 × 56 × 3 and 112 × 112 × 3 cells, respectively. For subregions below the unsampled region, the sampling level increases by one at each time when the subregion goes deeper. The resulting cell sizes from the fourth abstract subregion to the bottom abstract subregion are 112 × 112 × 6, 56 × 56 × 6, 28 × 28 × 5, 14 × 14 × 7 and 7 × 7 × 64 cells, respectively. This results in a total of 9.90 million edge reduction. For 0.02, 0.2 and 20 s, the abstract subsampling schemes of [S, N] = [2, 3; 1, 3; 0, 25; 1, 2; 2, 2; 3, 2; 4, 4; 5, 83], [S, N] = [2, 3; 1, 3; 0, 31; 1, 5; 2, 5; 3, 5; 4, 12; 5, 60] and [S, N] = [2, 3; 1, 3; 0, 46; 1, 6; 2, 8; 3, 8; 4, 10; 5, 40] are used.

Again, the numerical performance is examined by comparing our algebraic MRS technique with the SG scheme for this realistic model at periods of 0.02, 0.2, 2 and 20 s, listed in Table 6. The table shows that the algebraic MRS technique exclusively outperforms the SG scheme. Especially at 0.02 s, the computation time of our algebraic MRS method is merely 18 per cent of that (1346 s vs 7334 s) for SG. Similar as in the previous case, the decrease in computation time is larger than the reduction in the edge number.

Table 6.
Open in new tab

Numerical performance of SG scheme and algebraic MRS technique in terms of run time and the number of sampled edges for the terrain model.

Period
0.02 s0.2 s2 s20 s
MethodTime (s) and no. sampled edges (million)
SG733418.886,10018.886,18018.884,05018.88
MRS13468.072,0709.031,8338.982,02410.82
MRS/SG (per cent)1843344830485057
Period
0.02 s0.2 s2 s20 s
MethodTime (s) and no. sampled edges (million)
SG733418.886,10018.886,18018.884,05018.88
MRS13468.072,0709.031,8338.982,02410.82
MRS/SG (per cent)1843344830485057
Table 6.
Open in new tab

Numerical performance of SG scheme and algebraic MRS technique in terms of run time and the number of sampled edges for the terrain model.

Period
0.02 s0.2 s2 s20 s
MethodTime (s) and no. sampled edges (million)
SG733418.886,10018.886,18018.884,05018.88
MRS13468.072,0709.031,8338.982,02410.82
MRS/SG (per cent)1843344830485057
Period
0.02 s0.2 s2 s20 s
MethodTime (s) and no. sampled edges (million)
SG733418.886,10018.886,18018.884,05018.88
MRS13468.072,0709.031,8338.982,02410.82
MRS/SG (per cent)1843344830485057

For this realistic model, the accuracy of our algorithm is once again checked against SG. Fig. 14 compares the numerical solutions from the algebraic MRS algorithm and the SG scheme at 2 s along one line (x=2 km). They are in good agreement, with the relative difference in apparent resistivity below 3.3 per cent and the absolute phase differences below 2.2°.

Comparison of the algebraic MRS technique and SG scheme for the realistic model. (a) Apparent resistivities for XY polarization, (b) apparent resistivities for YX polarization, (d) phases for XY polarization and (e) phases for YX polarization. (c) and (f) are the relative difference for apparent resistivity and the phase difference against the SG solution, respectively.
Figure 14.

Comparison of the algebraic MRS technique and SG scheme for the realistic model. (a) Apparent resistivities for XY polarization, (b) apparent resistivities for YX polarization, (d) phases for XY polarization and (e) phases for YX polarization. (c) and (f) are the relative difference for apparent resistivity and the phase difference against the SG solution, respectively.

Open in new tabDownload slide

5 CONCLUSIONS

We have introduced an efficient algebraic MRS method for MT forward modelling on staggered FD grids based on the regularization equation, in which the traditional curl–curl equation is augmented with a scaled grad–div operator to enforce the divergence free current explicitly. The algebraic MRS scheme does not require the formation of subgrids as in Cherevatova et al. (2018). Once the abstract subsampling is defined, we can easily formulate the interpolation matrices based on the spatial relationship between retained and removed edges for each subsampling level. Then, we can subsample the grid based on matrix multiplication to generate linear systems of equations for different levels of subsampling. This avoids the discretization of the curl–curl equation on inconsistent subgrids as in Cherevatova et al. (2018).

The correctness of the algebraic MRS algorithm is verified by comparing our result with the reference solution based on a layered model. The accuracy and efficiency of the proposed algorithm are examined against the classic SG algorithm based on DTM1 and a realistic terrain model. In general, the algebraic MRS algorithm can reduce the degrees of freedom of the resulting linear system of equations by around half at an accuracy comparable to the SG solution, resulting in an efficiency gain of more than the reduction of the degrees of freedom.

For this algorithm, there are several obvious advantages against the MRG scheme proposed by Cherevatova et al. (2018). The most attractive one is that we do not need to discretize the curl–curl equation on different subgrids, which can be a challenge at interfaces of different subgrids. The bypass of the formation of subgrids can be considered as another benefit, simplifying the implementation of this algorithm for inversion. For the numerical performance in terms of efficiency and accuracy, no conclusion can be reached without a close examination on the two methods.

ACKNOWLEDGEMENTS

We would like to thank the associated editor Ute Weckmann, and two anonymous reviewers for their valuable comments and suggestions, which significantly improved the manuscript. We also would like to thank M. Cherevatova, G.D. Egbert and M. Yu. Smirnov for providing their code for comparison and M. J. Nam, H. J. Kim, Y. Song, T. J. Lee, J. S. Son and J. H. Suh for making valuable data publicly available. We are also grateful to the High Performance Computing Center of Central South University for partial support of this work. This work was financially supported by the National Natural Science Foundation of China (42130810, 72088101, 42074165 and 42174171), the Hunan Provincial Innovation Foundation for Postgraduate (CX20210271), the Department of Emergency Management of Hunan Province (2020-06), and the Open Research Fund Program of Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring (2021YSJS21).

DATA AVAILABILITY STATEMENT

The data and plot code underlying this paper is available in the GitHub repository (https://github.com/Li205001046/algebraic-multi-resolution).

References

Avdeev
D. B.
,
2005
.
Three-dimensional electromagnetic modelling and inversion from theory to application
,
Surv. Geophys.
,
26
(
6
),
767
–
799
.
10.1007/s10712-005-1836-x

Google Scholar

Crossref
Search ADS

 

Cai
H. Z.
,
Xiong
B.
,
Han
M.
,
Zhdanov
M.
,
2014
.
3D controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method
,
Comput.Geosci.
,
73
,
164
–
176
.
10.1016/j.cageo.2014.09.008

Google Scholar

Crossref
Search ADS

 

Carriere
P.
,
Jeandel
D.
,
1991
.
A 3D finite element method for the simulation of thermoconvective flows and its performances on a vector-parallel computer
,
Int. J. Numer. Methods Fluids
,
12
(
10
),
929
–
946
.
10.1002/fld.1650121003

Google Scholar

Crossref
Search ADS

 

Chen
C. J
,
Kruglyakov
M.
,
Kuvshinov
A.
,
2021
.
Advanced three-dimensional electromagnetic modelling using a nested integral equation approach
,
Geophys. J. Int.
,
226
(
1
),
114
–
130
.
10.1093/gji/ggab072

Google Scholar

Crossref
Search ADS

 

Cherevatova
M.
,
Egbert
G. D.
,
Smirnov
M. Y.
,
2018
.
A multi-resolution approach to electromagnetic modelling
,
Geophys. J. Int.
,
214
(
1
),
656
–
671
.
10.1093/gji/ggy153

Google Scholar

Crossref
Search ADS

 

Dong
H.
,
Egbert
G. D.
,
2019
.
Divergence-free solutions to electromagnetic forward and adjoint problems: a regularization approach
,
Geophys. J. Int.
,
216
(
2
),
906
–
918
.
10.1093/gji/ggy462

Google Scholar

Crossref
Search ADS

 

Egbert
G. D.
,
Kelbert
A.
,
2012
.
Computational recipes for electromagnetic inverse problems
,
Geophys. J. Int.
,
189
(
1
),
251
–
267
.
10.1111/j.1365-246X.2011.05347.x

Google Scholar

Crossref
Search ADS

 

Farquharson
C. G.
,
Miensopust
M. P.
,
2011
.
Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction
,
J. appl. Geophys.
,
75
(
4
),
699
–
710
.
10.1016/j.jappgeo.2011.09.025

Google Scholar

Crossref
Search ADS

 

Gao
J.
,
Smirnov
M.
,
Smirnova
M. Y.
,
Egbert
G. D.
,
2020
.
3-D DC Resistivity forward modeling using the multi-resolution grid
,
Pure Appl. Geophys.
,
177
,
2803
–
2819
.
10.1007/s00024-019-02365-3

Google Scholar

Crossref
Search ADS

 

Gao
J.
,
Smirnov
M.
,
Smirnova
M. Y.
,
Egbert
G. D.
,
2021
.
3-D time-domain electromagnetic modeling based on multi-resolution grid with application to geomagnetically induced currents
,
Phys. Earth planet. Inter.
,
312
(
12
),
106651
.
DOI: 10.1016/j.pepi.2021.106651
.
10.1016/j.pepi.2021.106651

Google Scholar

Crossref
Search ADS

 

Guo
Z.
et al. ,
2020
.
Geosciences in Central South University: A state-of-the-art review
,
J. Cent. South Univ.
,
27
(
4
),
975
–
996
.
10.1007/s11771-020-4347-5

Google Scholar

Crossref
Search ADS

 

Grayver
A. V.
,
Bürg
M.
,
2014
.
Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method
,
Geophys. J. Int.
,
198
(
1
),
110
–
125
.
10.1093/gji/ggu119

Google Scholar

Crossref
Search ADS

 

Grayver
A. V.
,
Kolev
T. V.
,
2015
.
Large-scale 3D geoelectromagnetic modeling using parallel adaptive high-order finite element method
,
Geophysics
,
80
(
6
),
E277
–
E291
.
10.1190/geo2015-0013.1

Google Scholar

Crossref
Search ADS

 

Haber
E.
,
Heldmann
S.
,
2007
.
An octree multigridmethod for quasi-static Maxwell’s equations with highly discontinuous coefficients
,
J. Comput. Phys.
,
223
,
783
–
796
.
10.1016/j.jcp.2006.10.012

Google Scholar

Crossref
Search ADS

 

Haber
E.
,
Ruthotto
L.
,
2014
.
A multiscale finite volume method for Maxwell’s equations at low frequencies
,
Geophys. J. Int.
,
199
(
2
),
1268
–
1277
.
10.1093/gji/ggu268

Google Scholar

Crossref
Search ADS

 

Han
B.
,
Li
Y.
,
Li
G.
,
2018
.
3D forward modeling of magnetotelluric fields in general anisotropic media and its numerical implementation in Julia
,
Geophysics
,
83
(
4
),
F29
–
F40
.
10.1190/geo2017-0515.1

Google Scholar

Crossref
Search ADS

 

Han
X.
,
Yin
C.
,
Su
Y.
,
Zhang
B.
,
Liu
Y.
,
Ren
X.
,
Ni
J.
,
Farquharson
C.G.
,
2022
.
3D finite-element forward modeling of airborne em systems in frequency-domain using octree meshes
,
IEEE Trans. Geosci. Remote Sens.
,
60
,
5912813
.
DOI: 10.1109/TGRS.2022.3153026
.
10.1109/TGRS.2022.3153026

Google Scholar

Crossref
Search ADS

 

Hu
Y.
,
Yang
D.
,
Li
Y.
,
Wang
Z.
,
Lu
Y.
,
2022
.
3-D numerical study on controlled source electromagnetic monitoring of hydraulic fracturing fluid with the effect of steel-cased wells
,
IEEE Trans. Geosci. Remote Sens.
,
60
,
4504210
.

Google Scholar

OpenURL Placeholder Text

 

Jahandari
H.
,
Farquharson
C.
,
2015
.
Finite-volume modelling of geophysical electromagnetic data on unstructured grids using potentials
,
Geophys. J. Int.
,
202
(
3
),
1859
–
1876
.
10.1093/gji/ggv257

Google Scholar

Crossref
Search ADS

 

Kelbert
A.
,
Meqbel
N.
,
Egbert
G. D.
,
Tandon
K.
,
2014
.
ModEM: a modular system for inversion of electromagnetic geophysical data
,
Comput. Geosci.
,
66
,
40
–
53
.
10.1016/j.cageo.2014.01.010

Google Scholar

Crossref
Search ADS

 

Kruglyakov
M.
,
Kuvshinov
A.
,
2018
.
Using high-order polynomial basis in 3-D EM forward modeling based on volume integral equation method
,
Geophys. J. Int.
,
213
(
2
),
1387
–
1401
.
10.1093/gji/ggy059

Google Scholar

Crossref
Search ADS

 

Li
G.
,
Li
Y.
,
Han
B.
,
Liu
Z.
,
2018
.
Application of the perfectly matched layer in 3-D marine controlled-source electromagnetic modelling
,
Geophys. J. Int.
,
212
(
1
),
333
–
344
.
10.1093/gji/ggx382

Google Scholar

Crossref
Search ADS

 

Li
J.
,
Liu
R.
,
Guo
R.
,
Wang
Y.
,
Wang
X.
,
2020
.
3D finite difference modeling of controlled-source electromagnetic response in frequency domain based on a modified curl-curl equation
,
J. appl. Geophys.
,
183
,
104202
.
DOI: 10.1016/j.jappgeo.2020.104202
.
10.1016/j.jappgeo.2020.104202

Google Scholar

Crossref
Search ADS

 

Li
J.
,
Liu
J.
,
Guo
R.
,
Liu
R.
,
Wang
Y.
,
Chen
H.
,
2022
.
Extension of the regularization technique to controlled-source electromagnetic modeling in general anisotropic conductivity media
,
Geophysics
,
87
(
4
),
E243
–
E251
.
10.1190/geo2021-0319.1

Google Scholar

Crossref
Search ADS

 

Liu
Y.
,
Yin
C.
,
2013
.
Electromagnetic divergence correction for 3D anisotropic EM modeling
,
J. appl. Geophys.
,
96
,
19
–
27
.
10.1016/j.jappgeo.2013.06.014

Google Scholar

Crossref
Search ADS

 

Miensopust
M. P.
,
Queralt
P.
,
Jones
A. G.
,
2013
.
Magnetotelluric 3-D inversion—a review of two successful workshops on forward and inversion code testing and comparison
,
Geophys. J. Int.
,
193
(
3
),
1216
–
1238
.
10.1093/gji/ggt066

Google Scholar

Crossref
Search ADS

 

Mittet
R.
,
2010
.
High-order finite-difference simulations of marine CSEM surveys using a correspondence principle for wave and diffusion fields
,
Geophysics
,
75
(
1
),
F33
–
F50
.
10.1190/1.3278525

Google Scholar

Crossref
Search ADS

 

Nam
M. J.
,
Kim
H. J.
,
Song
Y.
,
Lee
T. J.
,
Son
J. S.
,
Suh
J. H.
,
2007
.
3D magnetotelluric modelling including surface topography
,
Geophys. Prospect.
,
55
,
277
–
287
.
10.1111/j.1365-2478.2007.00614.x

Google Scholar

Crossref
Search ADS

 

Peng
R.
,
Han
B.
,
Liu
Y.
,
Hu
X.
,
2021
.
EM3DANI: A Julia package for fully anisotropic 3D forward modeling of electromagnetic data
,
Geophysics
,
86
(
5
),
WB1
–
WB8
.
10.1190/geo2020-0621.1

Google Scholar

Crossref
Search ADS

 

Ren
Z.
,
Kalscheuer
T.
,
Greenhalgh
S.
,
Maurer
H.
,
2013
.
A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling
,
Geophys. J. Int.
,
194
(
2
),
700
–
718
.
10.1093/gji/ggt154

Google Scholar

Crossref
Search ADS

 

Ren
Z.
,
Tang
J.
,
Kalscheuer
T.
,
Maurer
H.
,
2017
.
Fast 3-D large-scale gravity and magnetic modeling using unstructured grids and an adaptive multilevel fast multipole method
,
J. geophys. Res.: Solid Earth
,
122
(
1
),
79
–
109
.
10.1002/2016JB012987

Google Scholar

Crossref
Search ADS

 

Smith
J. T.
,
1996
.
Conservative modeling of 3D electromagnetic fields, Part II: biconjugate gradient solution and an accelerator
,
Geophysics
,
61
(
5
),
1319
–
1324
.
10.1190/1.1444055

Google Scholar

Crossref
Search ADS

 

Streich
R.
,
2009
.
3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy
,
Geophysics
,
74
(
5
),
F95
–
F105
.
10.1190/1.3196241

Google Scholar

Crossref
Search ADS

 

Tang
W.
,
Li
Y.
,
Liu
J.
,
Deng
J.
,
2021
.
Three-dimensional controlled-source electromagnetic forward modeling by edge-based finite element with a divergence correction
,
Geophysics
,
86
(
6
),
E367
–
E382
.
10.1190/geo2020-0520.1

Google Scholar

Crossref
Search ADS

 

Tang
W.
,
Huang
Q.
,
Deng
J.
,
Liu
J.
,
Zhou
F.
,
2022
.
Joint application of secondary field and coupled potential formulations to unstructured meshes for 3-D CSEM forward modeling
,
IEEE Trans. Geosci. Remote Sens.
,
60
,
5921409
.
DOI: 10.1109/TGRS.2022.3212884
.
10.1109/TGRS.2022.3212884

Google Scholar

Crossref
Search ADS

 

Wang
Y.
,
Jin
S.
,
Dong
H.
,
2021
.
Multi-level down-sampling scheme for accelerated solution in magnetotelluric forward modelling
,
J. appl. Geophys.
,
192
,
104384
.
DOI: 10.1016/j.jappgeo.2021.104384
.
10.1016/j.jappgeo.2021.104384

Google Scholar

Crossref
Search ADS

 

Yee
K. S.
,
1966
.
Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media
,
IEEE Trans. Antennas Propag.
,
14
(
3
),
302
–
307
.
10.1109/TAP.1966.1138693

Google Scholar

Crossref
Search ADS

 

Yu
N.
,
Li
R.
,
Kong
W.
,
Gao
L.
,
Wu
X.
,
Wang
E.
,
2022
.
A hybrid grid-based finite-element approach for three-dimensional magnetotelluric forward modeling in general anisotropic media
,
Comput.Geosci.
,
159
,
105035
.
DOI: 10.1016/j.cageo.2022.105035
.
10.1016/j.cageo.2022.105035

Google Scholar

Crossref
Search ADS

 

Zhdanov
M. S.
,
Lee
S. K.
,
Yoshioka
K.
,
2006
.
Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity
,
Geophysics
,
71
(
6
),
G333
–
G345
.
10.1190/1.2358403

Google Scholar

Crossref
Search ADS

 

Zhou
J.
,
Hu
X.
,
Xiao
T.
,
Cai
H.
,
Li
J.
,
Peng
R.
,
Long
Z.
,
2021
.
Three-dimensional edge-based finite element modeling of magnetotelluric data in anisotropic media with a divergence correction
,
J. appl. Geophys.
,
189
,
104324
.
DOI: 10.1016/j.jappgeo.2021.104324
.
10.1016/j.jappgeo.2021.104324

Google Scholar

Crossref
Search ADS

 
© The Author(s) 2023. Published by Oxford University Press on behalf of The Royal Astronomical Society.
This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
Subject
Subject
Issue Section:
Geomagnetism and Electromagnitism
Download all slides