Magnetotelluric Forward Modeling Using a Non-Uniform Grid Finite Difference Method

Abstract

Magnetotelluric (MT) forward modeling is essential in geophysical exploration, enabling the investigation of the Earth’s subsurface electrical conductivity. Traditional finite difference methods (FDMs) typically use uniform grids, which can be computationally inefficient and fail to accurately capture complex geological structures. This study addresses these challenges by introducing a non-uniform grid-based FDM for MT forward modeling. The proposed method optimizes computational resources by varying grid resolution, offering finer grids in areas with complex geology and coarser grids in more homogeneous regions. We apply this method to both typical synthetic models and a complex fault structure case study, demonstrating its capability to accurately resolve subsurface features while reducing computational costs. The results highlight the method’s effectiveness in capturing fine-scale details that are often missed by uniform grid approaches. The conclusions drawn from this study suggest that the non-uniform grid FDM not only improves the accuracy of MT modeling but also enhances its efficiency, making it a valuable tool for geophysical exploration in challenging environments.

1. Introduction

The magnetotelluric (MT) method is a geophysical exploration technique that utilizes the Earth’s natural electromagnetic field to investigate the electromagnetic properties of subsurface materials. This method is employed to explore groundwater resources, mineral deposits, oil and gas reservoirs, environmental geology, and geological disaster prediction [1,2,3,4,5,6]. MT research is currently advancing rapidly, with a focus on multi-dimensional forward modeling and inversion techniques [7,8,9,10,11,12]. MT forward modeling calculates the response of subsurface electromagnetic fields to infer the distribution of geological structures, including mineral layers, oil and gas reservoirs, aquifers, and pollution zones. This provides theoretical support for mineral exploration, oil and gas prospecting, environmental geology, and engineering surveys. Moreover, MT forward modeling is crucial for studying the Earth’s internal physical properties, revealing subsurface structures, and investigating dynamic geological processes.

Numerical simulation methods for MT forward modeling mainly include the Finite Element Method (FEM) [13,14,15,16], Spectral Method [17,18,19], Finite Difference Method (FDM) [20,21,22,23,24,25], and Meshless Method [26,27]. The Finite Difference Method is widely utilized due to its high computational accuracy, adaptability to complex structures, simplicity of implementation, and suitability for parallel processing. However, FDMs are typically associated with structured rectangular (2D) or hexahedral (3D) grids, with only a few exceptions. These grids are simple and easy to handle but can be inefficient, especially when adapting to the geometric complexities of geological structures. Complex geological structures require high resolutions in localized regions while allowing coarser grids in background layers. This mismatch between grid resolution and geological complexity can lead to either excessive computational demands or insufficient modeling accuracy.

Non-uniform grids present a promising alternative, offering variable grid spacing that adapts to subsurface complexity. This adaptability is particularly beneficial in MT modeling, where electromagnetic fields can vary significantly over short distances due to subsurface heterogeneities. Non-uniform grids enable finer resolution in regions requiring precise modeling, such as near fault zones or highly conductive layers, while employing coarser grids in less critical regions to reduce computational demands. As a result, a non-uniform grid was developed based on the uniform grid approach. Zhang, H. and Tang, X.G. [28] (2015) introduced a novel finite difference grid subdivision method, called variable step grid subdivision, for one-dimensional layered and continuous medium models in magnetotelluric sounding. Xu, Y. et al. [29] (2017) applied non-uniform grid subdivision techniques in finite difference transient electromagnetic 3D forward modeling, conducting forward calculations on uniform half-space models to verify the algorithm’s effectiveness and accuracy. Tong, X.Z. et al. [30] (2018) derived the linear equations for one-dimensional magnetotelluric forward modeling using the non-uniform grid finite difference method and conducted numerical simulations of the magnetotelluric responses for both a uniform half-space model and a two-layer D-type geoelectric model. Singh and Dehiya [31] (2023) use a coarse mesh to calculate boundary conditions for a finer mesh to improve the condition number of the system matrix. Researchers have validated the effectiveness of the non-uniform grid algorithm through model tests, observing that errors in the finite difference results are closely related to grid step size [32]. Grid step size greatly impacts both accuracy and computational efficiency, making reasonable grid subdivision essential.

In this study, we present a detailed analysis of magnetotelluric (MT) forward modeling using the finite difference method (FDM) with non-uniform grids. The primary aim of this work is to address the challenges posed by complex geological structures, which require high spatial resolution in certain regions while allowing coarser grids elsewhere to maintain computational efficiency. We develop and validate a non-uniform grid algorithm customized for the MT method, which accurately models electromagnetic fields in heterogeneous and anisotropic subsurface environments. Through a series of typical synthetic examples and a complex fault structure case study, we demonstrate the superiority of our approach in capturing fine-scale features without compromising computational efficiency. The results highlight the crucial role of grid adaptivity in geophysical modeling, laying the foundation for more accurate and efficient MT inversion processes.

2. Control Equations and Boundary Conditions

2.1. Control Equations of the Magnetotelluric Field

The propagation of electromagnetic waves in geoelectric media is governed by Maxwell’s equations. Assuming the Earth’s medium is isotropic and non-magnetic, the frequency-domain Maxwell’s equations for magnetotelluric sounding can be expressed as follows:

$$\nabla \times \mathbf{E}=-\mu {\displaystyle \frac{\partial \mathbf{H}}{\partial t}}$$

(1)

$$\nabla \times \mathbf{H}=\sigma \mathbf{E}+\epsilon {\displaystyle \frac{\partial \mathbf{E}}{\partial t}}$$

(2)

$$\nabla ·\mathbf{E}=0$$

(3)

$$\nabla ·\mathbf{H}=0$$

(4)

where $\mathbf{E}$ represents the electric field strength (V/m), $\mathbf{H}$ is the magnetic field strength (A/m), $\sigma$ denotes the electrical conductivity (S/m), $\mu$ and $\epsilon$ represent the permeability and permittivity of the medium, respectively, and $\omega$ is the angular frequency. By expanding Equations (1) and (2) of Maxwell’s equations, we derive

$$\overrightarrow{\mathit{i}}\left({\displaystyle \frac{\partial {\mathit{E}}_z}{\partial \mathit{y}}}-{\displaystyle \frac{\partial {\mathit{E}}_{\mathit{y}}}{\partial z}}\right)+\overrightarrow{\mathit{j}}\left({\displaystyle \frac{\partial {\mathit{E}}_{\mathit{x}}}{\partial z}}-{\displaystyle \frac{\partial {\mathit{E}}_z}{\partial \mathit{x}}}\right)+\overrightarrow{\mathit{k}}\left({\displaystyle \frac{\partial {\mathit{E}}_{\mathit{y}}}{\partial \mathit{x}}}-{\displaystyle \frac{\partial {\mathit{E}}_{\mathit{x}}}{\partial \mathit{y}}}\right)=\mathit{i}\mu \omega \left(\overrightarrow{\mathit{i}}{\mathit{H}}_{\mathit{x}}+\overrightarrow{\mathit{j}}{\mathit{H}}_{\mathit{y}}+\overrightarrow{\mathit{k}}{\mathit{H}}_z\right)$$

(5)

$$\overrightarrow{\mathit{i}}\left({\displaystyle \frac{\partial {\mathit{H}}_z}{\partial \mathit{y}}}-{\displaystyle \frac{\partial {\mathit{H}}_{\mathit{y}}}{\partial z}}\right)+\overrightarrow{\mathit{j}}\left({\displaystyle \frac{\partial {\mathit{H}}_{\mathit{x}}}{\partial z}}-{\displaystyle \frac{\partial {\mathit{H}}_z}{\partial \mathit{x}}}\right)+\overrightarrow{\mathit{k}}\left({\displaystyle \frac{\partial {\mathit{H}}_{\mathit{y}}}{\partial \mathit{x}}}-{\displaystyle \frac{\partial {\mathit{H}}_{\mathit{x}}}{\partial \mathit{y}}}\right)=\left(\sigma -\omega \epsilon \right)\left(\overrightarrow{\mathit{i}}{\mathit{E}}_{\mathit{x}}+\overrightarrow{\mathit{j}}{\mathit{E}}_{\mathit{y}}+\overrightarrow{\mathit{k}}{\mathit{E}}_z\right)$$

(6)

where $\overrightarrow{i}$, $\overrightarrow{j}$, and $\overrightarrow{k}$ are unit vectors and $\mathit{i}$ represents the imaginary number. In magnetotelluric sounding, the ratio of conduction current to displacement current is significantly greater than 1, which is ${\displaystyle \frac{\sigma }{\omega \epsilon }}\gg 1$. As a result, the effect of the displacement current on the field distribution in the Earth’s medium can be neglected. Therefore, $\sigma -\omega \epsilon$ can be considered as $\sigma$.

In geophysical exploration, geological structures like inclined strata, anticlines, and synclines with pronounced strike direction are typically modeled with the x-axis aligned along the strike, the y-axis perpendicular to the strike (dip direction), and the z-axis oriented vertically downward. The electrical properties of the medium vary along the y and z axes but remain constant along the strike direction, ${\displaystyle \frac{\partial }{\partial x}}=0$. Since the vector components in Equations (5) and (6) must be equal, we can derive the magnetotelluric field equations for a 2D medium

$${\displaystyle \frac{\partial {\mathit{H}}_z}{\partial y}}-{\displaystyle \frac{\partial {\mathit{H}}_y}{\partial z}}=\sigma {\mathit{E}}_x$$

(7)

$${\mathit{H}}_y={\displaystyle \frac{1}{\mathit{i}\omega \mu }}{\displaystyle \frac{\partial {\mathit{E}}_x}{\partial z}}$$

(8)

$${\mathit{H}}_z=-{\displaystyle \frac{1}{\mathit{i}\omega \mu }}{\displaystyle \frac{\partial {\mathit{E}}_x}{\partial y}}$$

(9)

$${\displaystyle \frac{\partial {\mathit{E}}_z}{\partial y}}-{\displaystyle \frac{\partial {\mathit{E}}_y}{\partial z}}=\mathit{i}\omega \mu {\mathit{H}}_x$$

(10)

$${\mathit{E}}_y={\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{\partial {\mathit{H}}_x}{\partial z}}$$

(11)

$${\mathit{E}}_z=-{\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{\partial {\mathit{H}}_x}{\partial y}}$$

(12)

Equation (7) through Equation (9) represent the partial differential equation for the magnetotelluric field in TE polarization mode, where the electric field component is perpendicular to the propagation direction of the electromagnetic wave, implying that only the magnetic field component exists along the propagation direction. Equation (10) through Equation (12) represent the partial differential equation for the magnetotelluric field in TM polarization mode, where the magnetic field component is perpendicular to the propagation direction of the electromagnetic wave, implying that the electric field component is zero in the propagation direction. Transforming Equations (7) to (12) into second-order partial differential equations that include only $E_x$ and $H_x$ yields

$${\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{1}{\mathit{i}\omega \mu }}{\displaystyle \frac{\partial {\mathit{E}}_x}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{1}{\mathit{i}\omega \mu }}{\displaystyle \frac{\partial {\mathit{E}}_x}{\partial z}}\right)+\sigma {\mathit{E}}_x=0$$

(13)

$${\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{1}{\sigma -\mathit{i}\omega \mu }}{\displaystyle \frac{\partial {\mathit{H}}_x}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{1}{\sigma -\mathit{i}\omega \mu }}{\displaystyle \frac{\partial {\mathit{H}}_x}{\partial z}}\right)+\mathit{i}\omega \mu {\mathit{H}}_x=0$$

(14)

Equation (13) represents the second-order partial differential equation for the magnetotelluric field in the TE polarization mode, and Equation (14) corresponds to the TM polarization mode.

2.2. Boundary Conditions

For the TE polarization mode, the top boundary AB must be positioned far enough from the ground surface, where the anomalous field is zero and $u$ is normalized to 1 unit at this boundary

$$u{|}_{z=0}=1$$

(15)

The bottom boundary is located at the interface between different subsurface media. When the area above the interface represents a geological anomaly while the area below consists of a homogeneous rock body, the propagation equation for electromagnetic waves below the interface is as follows:

$$u=u_0{e}^{kz}$$

(16)

where $u_0$ is a constant and k is the propagation coefficient. Taking the derivative of $u$ in Equation (16), we obtain ${\displaystyle \frac{\partial u}{\partial z}}=ku$. At the interface, ${\displaystyle \frac{\partial }{\partial z}}=-{\displaystyle \frac{\partial }{\partial n}}$. Here, $n$ is the normal direction of electromagnetic wave propagation. Therefore, the boundary condition at the subsurface interface is as follows:

$${\displaystyle \frac{\partial u}{\partial n}}+ku=0$$

(17)

The left and right boundaries are far enough from this study area, so

$${\displaystyle \frac{\partial u}{\partial n}}=0$$

(18)

In the TM polarization mode, the top boundary AB is placed directly on the ground surface, and u is normalized to 1 unit at this point; therefore

$$u{|}_{z=0}=1$$

(19)

The boundary conditions for the bottom, left, and right boundaries are the same as those for the TE polarization mode.

3. Finite Difference Discretization

3.1. Finite Difference Scheme Based on Non-Uniform Grids

In the 2D MT model, the non-uniform grid FDM is employed to solve the governing equations. First, the 2D MT model is meshed, as illustrated in Figure 1a, with non-uniform grid cells denoted as $\Delta y_i$ and $\Delta z_j$. The center of the grid region represents a geological anomaly. In contrast to uniform grid meshing, as shown in Figure 1b with a 50 m grid step, non-uniform grid meshing increases density around the geological anomaly and decreases it near the boundaries.

Let $u_{i,j}$ represent the magnetotelluric field at the node $\left(i,j\right)$. When the electrical conductivity of the medium is non-uniformly distributed, the function $u$ and its derivatives of all orders are single-valued, finite, and continuous functions of the variable $y$, and they can be expanded as a Taylor series

$$u\left(y+\Delta y\right)=u\left(y\right)+\Delta y{\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{1}{2}}\Delta y^2{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}+{\displaystyle \frac{1}{6}}\Delta y^3{\displaystyle \frac{{\partial }^{3}u}{\partial y^3}}+\cdots$$

(20)

$$u\left(y-\Delta y\right)=u\left(y\right)-\Delta y{\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{1}{2}}\Delta y^2{\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}-{\displaystyle \frac{1}{6}}\Delta y^3{\displaystyle \frac{{\partial }^{3}u}{\partial y^3}}+\cdots$$

(21)

In non-uniform grids, the step sizes $\Delta y_k$ and $\Delta z_j$ in the $y$ and $z$ directions are variable and may be unequal. So, Equations (20) and (21) are transformed into

$$u\left(y_i+\Delta y_{i+1}\right)=u\left(y_i\right)+\Delta y_{i+1}{\displaystyle \frac{\partial u\left(y_i\right)}{\partial y_i}}+{\displaystyle \frac{1}{2}}\Delta y_{i+1}^2{\displaystyle \frac{{\partial }^{2}u\left(y_i\right)}{\partial y_i^2}}+{\displaystyle \frac{1}{6}}\Delta y_{i+1}^3{\displaystyle \frac{{\partial }^{3}u\left(y_i\right)}{\partial y_i^3}}+\cdots$$

(22)

$$u\left(y_i-\Delta y_i\right)=u\left(y_i\right)-\Delta y_i{\displaystyle \frac{\partial u\left(y_i\right)}{\partial y_i}}+{\displaystyle \frac{1}{2}}\Delta y_i^2{\displaystyle \frac{{\partial }^{2}u\left(y_i\right)}{\partial y_i^2}}-{\displaystyle \frac{1}{6}}\Delta y_i^3{\displaystyle \frac{{\partial }^{3}u\left(y_i\right)}{\partial y_i^3}}+\cdots$$

(23)

From the manipulation (addition and subtraction) of Equations (22) and (23), the second-order derivatives of $u_{i,j}$ can be expressed as follows:

$${\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{1}{{\sigma }_{i,j}}}{\displaystyle \frac{\partial u_{i,j}}{\partial y}}\right)=2\left({\displaystyle \frac{u_{i+1,j}}{{\sigma }_{i+1,j}\Delta y_{i+1}\left(\Delta y_{i+1}+\Delta y_i\right)}}-{\displaystyle \frac{u_{i,j}}{{\sigma }_{i,j}\Delta y_i\Delta y_{i+1}}}+{\displaystyle \frac{u_{i-1,j}}{{\sigma }_{i-1,j}\Delta y_i\left(\Delta y_{i+1}+\Delta y_i\right)}}\right)$$

(24)

$${\displaystyle \frac{\partial }{\partial z}}({\displaystyle \frac{1}{{\sigma }_{i,j}}}{\displaystyle \frac{\partial u_{i,j}}{\partial z}})=2\left({\displaystyle \frac{u_{i,j+1}}{{\sigma }_{i,j+1}\Delta z_{j+1}\left(\Delta z_{j+1}+\Delta z_j\right)}}-{\displaystyle \frac{u_{i,j+1}}{{\sigma }_{i,j}\Delta z_j\Delta z_{j+1}}}+{\displaystyle \frac{u_{i,j-1}}{{\sigma }_{i,j-1}\Delta z_j\left(\Delta z_{j+1}+\Delta z_j\right)}}\right)$$

(25)

Let $\lambda =\mathit{i}\omega \mu$. Substituting Equations (24) and (25) into Equations (13) and (14) yields

$$\begin{array}{l}{\displaystyle \frac{2}{\lambda \Delta y_{i+1}\left(\Delta y_{i+1}+\Delta y_i\right)}}{E}_{xi+1,j}+{\displaystyle \frac{2}{\lambda \Delta y_i\left(\Delta y_{i+1}+\Delta y_i\right)}}{E}_{xi-1,j}+{\displaystyle \frac{2}{\lambda \Delta z_{j+1}\left(\Delta z_{j+1}+\Delta z_j\right)}}{E}_{xi,j+1}\\ +{\displaystyle \frac{2}{\lambda \Delta z_j\left(\Delta z_{j+1}+\Delta z_j\right)}}{E}_{xi,j-1}+({\sigma }_{i,j}-{\displaystyle \frac{2}{\lambda \Delta y_i\Delta y_{i+1}}}-{\displaystyle \frac{2}{\lambda \Delta z_j\Delta z_{j+1}}})E_{xi,j}=0\end{array}$$

(26)

$$\begin{array}{l}{\displaystyle \frac{2}{{\sigma }_{i+1,j}\Delta y_{i+1}\left(\Delta y_{i+1}+\Delta y_i\right)}}{H}_{xi+1,j}+{\displaystyle \frac{2}{{\sigma }_{i-1,j}\Delta y_i\left(\Delta y_{i+1}+\Delta y_i\right)}}{H}_{xi-1,j}+{\displaystyle \frac{2}{{\sigma }_{i,j+1}\Delta z_{j+1}\left(\Delta z_{j+1}+\Delta z_j\right)}}{H}_{xi,j+1}\\ +{\displaystyle \frac{2}{{\sigma }_{i,j-1}\Delta z_j\left(\Delta z_{j+1}+\Delta z_j\right)}}{H}_{xi,j-1}+(\lambda -{\displaystyle \frac{2}{{\sigma }_{i,j}\Delta y_i\Delta y_{i+1}}}-{\displaystyle \frac{2}{{\sigma }_{i,j}\Delta z_j\Delta z_{j+1}}})H_{xi,j}=0\end{array}$$

(27)

Equations (29) and (30) represent the finite difference scheme for the non-uniform grid used to calculate the wavefield value at node $\left(i,j\right)$.

3.2. Deriving the Difference Equations

To construct the discretized partial differential equation system, the nodes on the 2D grid need to be flattened into a 1D array, facilitating the use of linear algebra methods to solve the resulting equations. $N_y\times N_z$ represents the total number of nodes on the grid. This array will be used to construct the coefficient matrix $A$ and the right-hand vector $b$ of the equation system. After flattening the 2D grid into a 1D array, a new index vector $k=(i-1)\times nz+j,i=1,\cdots ny;j=1,\cdots nz$ can be used to represent each node.

Equation (27) can be written in matrix form as

$$\begin{array}{l}\left(\begin{array}{ccccccc}\lambda -{\displaystyle \frac{2}{{\sigma }_{i-1,j}}}\left({\displaystyle \frac{1}{\Delta y_{i-1}\Delta y_i}}+{\displaystyle \frac{1}{\Delta z_j\Delta z_{j+1}}}\right)& \cdots & 0& {\displaystyle \frac{2}{{\sigma }_{i,j}\Delta y_i\left(\Delta y_{i-1}+\Delta y_i\right)}}& 0& \cdots & 0\\ ⋮& \cdots & ⋮& ⋮& ⋮& \cdots & ⋮\\ 0& \cdots & \lambda -{\displaystyle \frac{2}{{\sigma }_{i,j-1}}}\left({\displaystyle \frac{1}{\Delta y_i\Delta y_{i-1}}}+{\displaystyle \frac{1}{\Delta z_j\Delta z_{j-1}}}\right)& {\displaystyle \frac{2}{{\sigma }_{i,j+1}\Delta z_j\left(\Delta z_{j-1}+\Delta z_j\right)}}& 0& \cdots & 0\\ {\displaystyle \frac{2}{{\sigma }_{i-1,j}\Delta y_i\left(\Delta y_{i+1}+\Delta y_i\right)}}& \cdots & {\displaystyle \frac{2}{{\sigma }_{i,j-1}\Delta z_j\left(\Delta z_{j+1}+\Delta z_j\right)}}& \lambda -{\displaystyle \frac{2}{{\sigma }_{i,j}}}\left({\displaystyle \frac{1}{\Delta y_i\Delta y_{i+1}}}+{\displaystyle \frac{1}{\Delta z_j\Delta z_{j+1}}}\right)& {\displaystyle \frac{2}{{\sigma }_{i,j+1}\Delta z_{j+1}\left(\Delta z_{j+1}+\Delta z_j\right)}}& \cdots & {\displaystyle \frac{2}{{\sigma }_{i+1,j}\Delta y_{i+1}\left(\Delta y_{i+1}+\Delta y_i\right)}}\\ 0& \cdots & 0& {\displaystyle \frac{2}{{\sigma }_{i,j}\Delta z_{j+1}\left(\Delta z_{j+2}+\Delta z_{j+1}\right)}}& \lambda -{\displaystyle \frac{2}{{\sigma }_{i,j+1}}}\left({\displaystyle \frac{1}{\Delta y_i\Delta y_{i+1}}}+{\displaystyle \frac{1}{\Delta z_{j+2}\Delta z_{i+1}}}\right)& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& \cdots & 0& {\displaystyle \frac{2}{{\sigma }_{i,j}\Delta y_{i+1}\left(\Delta y_{i+1}+\Delta y_{i+2}\right)}}& 0& \cdots & \lambda -{\displaystyle \frac{2}{{\sigma }_{i+1,j}}}\left({\displaystyle \frac{1}{\Delta y_{i+2}\Delta y_{i+1}}}+{\displaystyle \frac{1}{\Delta z_j\Delta z_{j+1}}}\right)\end{array}\right)\\ \left(\begin{array}{c}{H}_{x(i-1)\times nz+j}\\ ⋮\\ H_{xi\times nz+j-1}\\ H_{xi\times nz+j}\\ H_{xi\times nz+j+1}\\ ⋮\\ H_{x(i+1)\times nz+j}\end{array}\right)=\left(\begin{array}{c}0\\ ⋮\\ 0\\ 0\\ 0\\ ⋮\\ 0\end{array}\right)\end{array}$$

(28)

Let $u={\left(u_0,u_1,\cdots ,u_{i\times nz+j},\cdots ,u_{i\times nx+nz}\right)}^{\prime }$ and $b={\left(0,0,\cdots ,0,\cdots ,0\right)}^{\prime }$, and $A$ denotes the sparse matrix. Extending the matrix of Equations (26) and (27) to all nodes, we obtain

$$Au=b$$

(29)

$$A=\left(\begin{array}{cccccc}{\displaystyle \frac{2}{{\sigma }_{0,0}\Delta z_1\left(\Delta z_2+\Delta z_1\right)}}& \lambda -{\displaystyle \frac{2}{{\sigma }_{1,0}}}\left({\displaystyle \frac{1}{\Delta y_1\Delta y_2}}+{\displaystyle \frac{1}{\Delta z_1\Delta z_2}}\right)& \cdots & 0& \cdots & 0\\ 0& {\displaystyle \frac{2}{{\sigma }_{0,1}\Delta z_2\left(\Delta z_3+\Delta z_2\right)}}& \cdots & 0& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮& \cdots & ⋮\\ 0& 0& \cdots & \lambda -{\displaystyle \frac{2}{{\sigma }_{i,j}}}\left({\displaystyle \frac{1}{\Delta y_i\Delta y_{i+1}}}+{\displaystyle \frac{1}{\Delta z_j\Delta z_{j+1}}}\right)& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& \cdots & {\displaystyle \frac{2}{{\sigma }_{nx,nz}\Delta y_{nz}\left(\Delta y_{nx-1}+\Delta y_{nx}\right)}}\end{array}\right),$$

This system consists of $(N_y-1)\times (N_z-1)$ equations and $N_y\times N_z$ unknowns.

From Equation (19), the solution for the top boundary condition is derived as follows: $u_{z=0}=1$.

From Equations (17) and (20), the solution for the bottom boundary condition is derived as follows:

$${\displaystyle \frac{\partial u_{i,j}}{\partial z}}+k{u}_{i,j}{|}_{i=ny,j=nz}=0$$

(30)

$${\displaystyle \frac{u_{i,j}-u_{i-1,j}}{\Delta z_j}}+k{u}_{i,j}=0$$

(31)

$$({\displaystyle \frac{1}{\Delta z_j}}+k)u_{i,j}{|}_{i=ny,j=nz}-{\displaystyle \frac{1}{\Delta z_j}}{u}_{i-1,j}{|}_{i=ny,j=nz}=0$$

(32)

From Equations (18) and (20), the conditions for the left and right boundaries are derived as follows:

$${\displaystyle \frac{\partial u_{i,j}}{\partial y}}{|}_{i=y\min }={\displaystyle \frac{\partial u_{i,j}}{\partial y}}{|}_{i=y\max }=0$$

(33)

Incorporating boundary conditions into Equations (19), (32) and (33), the sparse matrix $A$ and $b$ becomes

$$\begin{array}{}A=\left(\begin{array}{cccccc}1& 0& \cdots & 0& \cdots & 0\\ {\displaystyle \frac{2}{{\sigma }_{0,0}\Delta z_1\left(\Delta z_2+\Delta z_1\right)}}& \lambda -{\displaystyle \frac{2}{{\sigma }_{1,0}}}\left({\displaystyle \frac{1}{\Delta y_1\Delta y_2}}+{\displaystyle \frac{1}{\Delta z_1\Delta z_2}}\right)& \cdots & 0& \cdots & 0\\ 0& {\displaystyle \frac{2}{{\sigma }_{0,1}\Delta z_2\left(\Delta z_3+\Delta z_2\right)}}& \cdots & 0& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮& \cdots & ⋮\\ 0& 0& \cdots & \lambda -{\displaystyle \frac{2}{{\sigma }_{i,j}}}\left({\displaystyle \frac{1}{\Delta y_i\Delta y_{i+1}}}+{\displaystyle \frac{1}{\Delta z_j\Delta z_{j+1}}}\right)& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& -{\displaystyle \frac{1}{\Delta z_{nz}}}& {\displaystyle \frac{1}{\Delta z_{nz}+k}}\end{array}0\right),\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}b={\left(1,0,0,\cdots ,0,\cdots ,0\right)}^{\prime }\end{array}$$

At this stage, the system of equations contains $(N_y+1)\times (N_z+1)$ equations and $(N_y+1)\times (N_z+1)$ unknowns. Solving Equation (29) yields the magnetic field values at the nodes. Similarly, by deriving the corresponding sparse coefficient matrix from Equation (26) and solving it, the electric field values at the nodes can be obtained.

3.3. Solving the Linear System

3.3.1. The BICGSTAB Algorithm

In the 2D MT forward simulation using the non-uniform grid finite difference method, the approach to solving the linear system significantly affects simulation accuracy and computational efficiency. As the grid density increases, the number of non-zero elements in matrix A rises, resulting in greater computational complexity and longer solving times, which reduce overall computational efficiency. Consequently, choosing an appropriate solution method is crucial.

In 2D magnetotelluric (MT) forward modeling, commonly used algorithms for solving large sparse matrices include the Gaussian elimination method, the Jacobi iteration method, the Gauss–Seidel iteration method, the Conjugate Gradient (CG) method, the Bi-Conjugate Gradient (BICG) method, and the Bi-Conjugate Gradient Stabilized (BICGSTAB) method [33,34,35,36,37].

The Gaussian elimination method transforms the matrix into echelon form through row operations, followed by back-substitution. This method is suitable for small-scale problems but becomes computationally intensive for larger ones. The Jacobi iteration method decomposes the matrix into the sum of a diagonal matrix, an upper triangular matrix, and a lower triangular matrix, and then solves it iteratively. It is effective for diagonally dominant matrices but has a slow convergence rate. The Gauss–Seidel iteration method is similar to the Jacobi method but uses updated variable values during iteration. It generally converges faster than the Jacobi method. The Conjugate Gradient (CG) method is an iterative technique based on Krylov subspaces and is suitable for symmetric positive-definite matrices. It converges quickly, particularly for large sparse problems. When using the CG method, it is essential to ensure matrix symmetry and positive definiteness. The Bi-Conjugate Gradient (BICG) method is similar to the CG method but is applicable to non-symmetric matrices. It typically converges faster than the Jacobi and Gauss–Seidel methods. When using the BICG method, it is important to consider matrix non-symmetry and convergence properties. The Bi-Conjugate Gradient Stabilized (BICGSTAB) method is an enhanced version of the BICG method, offering improved stability and convergence, and it is suitable for non-symmetric and ill-conditioned problems. When using the BICGSTAB algorithm, its stability and convergence must be considered, and an appropriate preconditioner should be selected to enhance convergence and improve solution efficiency.

To compare the computational efficiency and accuracy of the aforementioned methods, a sparse linear matrix system with a dimension of 1000 was created using a maximum of 100 iterations and a tolerance of 1 × 10⁻¹². The Gauss–Seidel, CG, BICG, and BICGSTAB methods were employed to solve the system. The relationship between the relative residual and the number of iterations is illustrated in Figure 2.

The relative error curves in Figure 2 demonstrate that the BICGSTAB algorithm converges to the maximum tolerance with significantly fewer iterations compared to other methods, thereby greatly enhancing computational efficiency. In solving the 2D MT forward modeling problem with non-uniform grid finite difference equations, the partitioning creates an asymmetric and ill-conditioned coefficient matrix, which makes direct methods prone to errors. Therefore, iterative methods are more appropriate. The BICGSTAB algorithm is particularly suited for asymmetric and ill-conditioned problems while also meeting computational efficiency requirements. Thus, this study utilizes the BICGSTAB algorithm to solve the large sparse linear matrix systems arising from non-uniform grid finite difference equations.

The basic steps of the BICGSTAB algorithm [38] are outlined in Figure 3.

First, given the linear system $Ax=b$, select an initial guess $x_0$ and compute the residual $r_0=b-A{x}_0$. Set $p_0=r_0$, $p_0=\alpha =\omega =1$, and $p=v=0$. Compute $p_k=\left(r_0,r_k\right)$, where $p_k$ is an index vector. If $p_k=0$, the algorithm fails and needs to adjust $r_0$ and continue.

If $p_k\ne 0$, the algorithm converges. Compute the coefficient $\beta ={\displaystyle \frac{p_k\alpha }{p_k{\omega }_{k-1}}}$ and update the index vector $p_k=r_{k-1}+{\beta }_k\left(p_{k-1}-{\omega }_{k-1}{v}_{k-1}\right)$. Continue by computing the intermediate vector $v_k={\mathbf{A}}^{\ast }{p}_k$ and the coefficient $\alpha ={\displaystyle \frac{p_k}{\left(r_0,v_k\right)}}$ to update the solution vector $x_k=x_{k-1}+{\alpha }^{\ast }{p}_k$.

Update the residual $s_k=r_k-{\alpha }^{\ast }{v}_k$ and check for convergence. If $‖s_k‖<\epsilon$, where $\epsilon$ is the given tolerance, the algorithm has converged. Next, compute the intermediate vector $t_k={\mathbf{A}}^{\ast }{s}_k$ and the coefficient $\omega ={\displaystyle \frac{\left(t_k,s_k\right)}{{‖t_k‖}^2}}$ to update the solution vector $x_k=x_k+\omega s_k$.

Update the residual $r_k=s_k-{\omega }^{\ast }{t}_k$. When $‖r_k‖<\epsilon$, the algorithm converges, and $x_k$ is returned as the solution. Upon completion of the iterations, $x_k$ is the approximate solution to the linear system $Ax=b$.

For asymmetric and ill-conditioned matrix problems, BICGSTAB generally provides superior convergence speed and stability compared to other iterative methods because it does not require computing a preconditioner matrix, thereby reducing computational complexity. However, the convergence speed of the BICGSTAB algorithm depends on the characteristics of the linear system. When addressing ill-conditioned problems and unfavorable eigenvalues, convergence may be slow, necessitating the use of preconditioning techniques to optimize the algorithm’s performance.

3.3.2. Incomplete LU Preconditioning Method

Incomplete LU preconditioning (ILU) is a practical technique that effectively improves the condition number of linear systems, thereby enhancing the convergence speed of iterative solvers [39]. ILU preconditioning constructs a factorization that approximates the original matrix A, i.e., $A\approx LU$, where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. The core of ILU preconditioning is that it allows $L$ and $U$ to maintain a certain level of sparsity, which helps reduce computational complexity and alleviate memory burden.

The basic steps of the ILU preconditioning method are as follows.

First, given a linear system $Ax=b$, set a fill-in threshold $\tau$ and create an initial copy of matrix $A$, denoted as $A^{’}$.

Next, traverse the non-zero elements of matrix $A^{’}$ one by one, performing LU factorization based on Gaussian elimination principles while retaining only the positions of non-zero elements in the original matrix $A$. During the factorization process, update only the positions of existing non-zero elements in the original matrix, ignoring values at other positions. If a newly computed value is below the threshold $\tau$, it is considered zero.

Obtain the approximate factorization $A\approx LU$, where the sparsity of $L$ and $U$ is similar to that of the original matrix $A$.

After completing ILU preconditioning, it can be applied to iterative solvers. In each iteration, the solver needs to solve a linear system of the form $L\ast U_y=r$, where $r$ is the residual vector. Given that $L$ and $U$ are both triangular matrices, this linear system can be efficiently solved using forward and backward substitution.

Figure 4 shows the relative residual curves of the BICGSTAB algorithm with and without ILU preconditioning. Comparative analysis of the figure reveals that, under the same tolerance of 1 × 10⁻¹² and a maximum of 100 iterations, the BICGSTAB algorithm with ILU preconditioning converges faster and achieves a smaller relative residual. Therefore, the computational efficiency and accuracy of the BICGSTAB algorithm are significantly improved with ILU preconditioning.

3.4. Calculation of Electromagnetic Response

After solving Equation (29) using the ILU preconditioned BICGSTAB method to obtain the $u$ values for each node, the partial derivative of the field value in the vertical direction ${\displaystyle \frac{\partial u}{\partial z}}$ is calculated numerically. To improve accuracy, the calculations typically involve the four nodes closest to the surface.

In the TE mode, ${\displaystyle \frac{\partial u}{\partial z}}$ corresponds to ${\displaystyle \frac{\partial E_x}{\partial z}}$. By substituting into the Cagniard apparent resistivity formulas (Equations (34) through (36)), the apparent resistivity and impedance phase for the TE mode can be determined.

$$Z_{xy}=E_x/\left(\lambda {\displaystyle \frac{\partial E_x}{\partial z}}\right)$$

(34)

$${\rho }_{xy}={\displaystyle \frac{1}{\omega \mu }}|Z_{xy}{|}^2$$

(35)

$${\phi }_{xy}=\arctan {\displaystyle \frac{\mathrm{Im}\left[Z_{xy}\right]}{\mathrm{Re}[Z_{xy}]}}$$

(36)

where $Z_{xy}$ represents the impedance, ${\rho }_{xy}$ represents the apparent resistivity, and ${\varphi }_{xy}$ denotes the impedance phase in the TE polarization mode. In the TM mode, ${\displaystyle \frac{\partial u}{\partial z}}$ corresponds to ${\displaystyle \frac{\partial H_x}{\partial z}}$. By substituting into the Cagniard apparent resistivity formulas (Equations (37) through (39)), the apparent resistivity and impedance phase for the TM mode can be determined.

$$Z_{yx}=-{\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{\partial H_x}{\partial z}}/H_x$$

(37)

$${\rho }_{yx}={\displaystyle \frac{1}{\omega \mu }}|Z_{yx}{|}^2$$

(38)

$${\phi }_{yx}=\arctan {\displaystyle \frac{\mathrm{Im}\left[Z_{yx}\right]}{\mathrm{Re}[Z_{yx}]}}$$

(39)

where $Z_{yx}$ represents the impedance, ${\rho }_{yx}$ represents the apparent resistivity, and ${\varphi }_{yx}$ denotes the impedance phase in the TM polarization mode.

4. Results and Discussion

4.1. Homogeneous Half-Space Model Modeling and Analysis

To verify the accuracy and efficiency of the non-uniform grid finite difference method, a homogeneous half-space model with dimensions of 1 km by 1 km and a resistivity of 100 Ω⋅m was constructed. Both uniform and non-uniform grid partitioning methods were employed. The uniform grid partitioning had a cell size of ∆y = ∆z = 10 m, resulting in a total of 100 × 100 grids. The non-uniform grid partitioning used cell sizes of 16 m and 8 m, with a total of 75 × 75 grids. Since the model represents a homogeneous half-space without any subsurface anomalies, the same non-uniform grid was applied in both horizontal and vertical directions. Here, the non-uniform grid is employed to verify the algorithm’s capability to accurately resolve subsurface features.

By comparing the magnetic field intensity curves of the central grid cells from both partitioning methods (Figure 5a), it is evident that the non-uniform grid finite difference method achieves computational accuracy comparable to that of the uniform grid finite difference method. The numerical solutions are close, with minimal errors, and the non-uniform grid finite difference method does not exhibit dispersion due to changes in grid size.

Additionally, as shown in Figure 5b, during the 2D magnetotelluric forward modeling of the uniform half-space model, the total computation time for the non-uniform grid finite difference method is significantly less than that of the uniform grid finite difference method, with computational efficiency nearly doubled. To further enhance the comparative analysis, we introduced an alternative non-uniform grid partitioning scheme. In this scheme, the non-uniform grid employed cell sizes of 20 m and 10 m, distributed over a total of 60 × 60 grid points. Its computational efficiency is further improved, with a total runtime of only 1.57 s. This represents a twofold reduction in computation time compared to the 75 × 75 non-uniform grid, which required 3.45 s.

4.2. Layered Model Modeling and Analysis

To further validate the effectiveness of the non-uniform grid finite difference method, four distinct layered models A-type, Q-type, H-type, and K-type were constructed and analyzed using this approach. Figure 6 illustrates these models, and Figure 7 shows curves of the apparent resistivity and impedance phase of different models.

Figure 7a–d show the result of the four models. In the A-type model (Figure 7a), the apparent resistivity values increase with increasing period, corresponding to higher resistivity at greater depths. The apparent resistivity curve shows a consistent upward trend as the frequency decreases, indicating the typical characteristics of a stratified formation with increasing resistivity at greater depths. Similarly, the impedance phase curve follows a comparable pattern, gradually rising as the frequency decreases. This behavior reflects the electrical properties of a layered geological structure, where resistivity increases with depth. The excellent agreement between the numerical results and the analytical solution demonstrates the accuracy of this method in capturing the electromagnetic response of such geological settings.

In the Q-type model (Figure 7b), the apparent resistivity decreases with increasing period, reflecting a reduction in resistivity at greater depths. As the frequency decreases, the apparent resistivity curve exhibits a clear downward trend, indicating the presence of more conductive layers at depth. Similarly, the impedance phase curve follows a downward trend, further confirming the decreasing resistivity with depth.

In the H-type model (Figure 7c), the apparent resistivity curve displays a noticeable minimum in the mid-frequency range, indicating the presence of a low-resistivity layer between two high-resistivity layers. This behavior is also reflected in the impedance phase curve, which shows a corresponding minimum, further highlighting the conductive nature of the intermediate layer. Such a pattern suggests a conductive zone, possibly representing a groundwater layer or a region with high fluid content, situated beneath resistive rock formations.

In the K-type model (Figure 7d), the apparent resistivity curve exhibits a pronounced peak in the mid-frequency range, indicating a high-resistivity layer at intermediate depths. This peak reflects the presence of a resistive anomaly, such as a mineral deposit, surrounded by more conductive layers. The impedance phase curve remains relatively smooth but shows corresponding variations that align with the resistivity changes. The behavior of both curves confirms the model’s ability to capture the electromagnetic signature of a high-resistivity layer embedded in a conductive background.

The numerical simulation results (red and blue curves) closely match the analytical solutions (black curves), demonstrating the correctness and reliability of the non-uniform grid finite difference algorithm. This method effectively handles various types of layered formations, capturing the electromagnetic responses with high fidelity and validating their applicability to complex geophysical problems.

4.3. Fault Structure Model Modeling and Analysis

To further assess the adaptability of the non-uniform grid finite difference algorithm and to explore the electromagnetic response characteristics within fault systems, a geoelectrical model was constructed as depicted in Figure 8. The model dimensions are 2400 m × 1000 m, with a fault width of 100 m. The resistivity values are set uniformly at 6 Ω·m for the fault, while the resistivity of the surrounding strata from top to bottom is 200 Ω·m, 1000 Ω·m, and 2500 Ω·m (representing the basement), respectively.

Grid discretization plays a critical role in the forward modeling of fault systems. Utilizing large-scale grids similar to the surrounding rock can lead to a “stair-step” structure along the fault plane, which may distort the results. Conversely, applying uniformly fine grids across the entire model would significantly increase computational load and storage requirements, thereby reducing efficiency. To overcome this challenge, a strategy of sparse boundary and dense center grid discretization was employed. For the background strata far from the fault, a larger grid size was used, while a finer grid was applied in and near the fault structure. The lateral grid size transitions from 200 m, 100 m, 50 m, and 40 m down to 10 m, with the 10 m grid specifically refined for the fault structure. Similarly, the vertical grid size transitions from 100 m, 50 m, and 20 m down to 10 m. Specifically, the grid for the model consists of 87 nodes in the horizontal direction and 64 nodes in the vertical direction, resulting in a total of 5568 nodes. Figure 8 illustrates the grid discretization scheme used in the forward modeling.

Figure 9 presents the apparent resistivity and impedance phase pseudo-sections for the fault model. The TM mode pseudo-section reveals significant resistivity variations with depth, characterized by a prominent low-resistivity zone near the center. This anomaly suggests the presence of a fault or a conductive geological structure that deepens towards the center, potentially indicating fault-related fluid saturation or a highly conductive fault zone. In contrast, the TE mode pseudo-section displays more consistent high resistivity with fewer pronounced anomalies, indicating less sensitivity to the fault structure compared to the TM mode.

The impedance phase pseudo-section for the TM mode exhibits a more gradual variation across the frequency range, with smoother transitions and fewer localized anomalies. This pattern suggests a relatively uniform subsurface structure in the TM mode, with a consistent increase in phase at lower frequencies, indicative of deeper layers with higher resistivity. Conversely, the impedance phase pseudo-section for the TE mode reveals distinct localized high-phase anomalies around the central region. These anomalies suggest the presence of conductive features or structures, such as fault zones or highly conductive layers, which are more pronounced in the TE mode due to its sensitivity to lateral variations in resistivity. This comparison highlights the differing sensitivities of the TM and TE modes in detecting subsurface resistivity variations, with the TE mode being more responsive to lateral conductivity contrasts.

5. Conclusions

In this paper, we developed and implemented a finite difference method (FDM) based on non-uniform grids for magnetotelluric (MT) forward modeling. By introducing of non-uniform grids, we effectively addressed the limitations of traditional uniform grids, particularly when modeling complex geological structures that require variable spatial resolution. This innovative approach significantly enhances computational efficiency and result accuracy while reducing computational costs without compromising quality.

Our study validated the robustness and flexibility of the method through applications to synthetic MT data, demonstrating its capability to handle a variety of geological scenarios. The non-uniform grid approach provides a crucial advantage by balancing model resolution with computational efficiency, making it an important tool for large-scale MT surveys. Additionally, the non-uniform grid approach provides a crucial advantage by balancing model resolution with computational efficiency, making it a valuable tool for large-scale MT surveys.

However, this method faces challenges in resolving highly complex domains, particularly those with sharp variations or extreme heterogeneity. In such cases, adaptive grid refinement or alternative meshing strategies may be necessary to accurately capture detailed electromagnetic responses.

Future research will focus on several key directions: first, optimizing grid generation algorithms to improve the efficiency and adaptability of non-uniform grids for more complex geological environments; second, extending this approach to 3D modeling to further enhance its applicability and accuracy in complex geological settings. Additionally, conducting a theoretical sensitivity analysis of the parameters will help deepen our understanding of their impact on the results, thereby optimizing parameter selection. Finally, exploring the integration of adaptive grid refinement methods will improve modeling precision in highly heterogeneous areas.

Overall, this work advances MT forward modeling techniques by offering a reliable and computationally efficient method that can be readily adopted in geophysical research and exploration.