Research Paper

Journal of the Korean Society of Mineral and Energy Resources Engineers. 31 August 2024. 284-295
https://doi.org/10.32390/ksmer.2024.61.4.284

ABSTRACT


MAIN

  • Introduction

  • Theory and Methods

  •   Maxwell’s Equations in 1-D Anisotropic Media

  •   Analytic Solutions using the Matrix Propagation Method

  • Stabilization strategies for numerical computation of 1-D solutions

  •   Stabilization strategy I

  •   Stabilization Strategy II

  • Numerical Experiments

  •   Verification of the two proposed strategies

  •   Numerical Experiments

  • Conclusions

Introduction

The magnetotelluric (MT) method, which is one of the passive geophysical methods utilizing natural electromagnetic fields, has extensively been employed for various purposes, such as analysis of subsurface electrical structures, and mining, geothermal and petroleum explorations. In most cases, MT data have been interpreted under the assumption of isotropic media. However, when the electrical properties of the Earth interior are remarkably anisotropic (Liu et al., 2018; Wannamaker, 2005), anisotropic properties need to be considered to accurately interpret MT data.

In recent years, a number of studies have been dedicated to MT modeling and inversion for 1-, 2- or 3-dimensional (3-D) anisotropic media (Bai et al., 2022; Cao et al., 2021; Rong et al., 2022). 1-D MT problems can be analytically solved from Maxwell’s equation through matrix propagation methods (O’Brien and Morrison, 1967; Reddy and Rankin, 1971). The EM field obtained from 1-D anisotropic modeling is employed to deal with the boundary conditions of MT modeling in 2- or 3-D anisotropic media (Li and Pek, 2008; Xiao et al., 2018).

In the case of the MT method, which in need of considering frequencies as low as 0.01~0.001 Hz, the depth of the model used for forward modeling should be of sufficient depth such that the components of EM fields become zero at the bottom boundary, i.e., approximately five to six times of skin depth. For instance, in both Dublin test model 1 (DTM1) utilized by Usui (2015) and Iheya North Knoll used by Usui et al. (2018), the depth of the computational domain was set to be 3000 km. However, the matrix propagation technique can be unstable, suffering from NaN (not a number) values, when the EM fields are computed at depths or at layers extending to several thousand kilometers.

In this study, we propose strategies to stably compute analytic solutions using the matrix propagation method for anisotropic layered models at any depth, including great depths and thick layers. The first strategy addresses instability at great depths by converting upgoing and downgoing EM fields into numerically stable exponential forms, ensuring solutions converge to zero as depth increases. The second strategy handles very thick layers by checking if the determinant of the matrix for each layer is near zero and replacing such unstable layers with an infinite half-space. By applying our strategies to 1-D anisotropic layered model, we show that our algorithm yields stable results even at great depths and thick layers. In addition, we apply the proposed algorithm to both 3-D anisotropic and isotropic media to verify the accuracy and effectiveness of the proposed strategies.

Theory and Methods

After deriving governing equations from Maxwell’s equations for an anisotropic 1-D layered model, we explain the process of computing the analytic solutions for anisotropic layers using the matrix propagation method.

Maxwell’s Equations in 1-D Anisotropic Media

For MT modeling in 1-D anisotropic media, governing equations expressed only in terms of the electric field can be derived from Maxwell’s Equations. Under the assumption that the Earth’s conductivity varies only with depth and is uniform in horizontal direction, i.e., x=y=0, the corresponding equations can be defined as follows (Pek and Santos, 2002):

(1)
2Exz2=iωμ0(AxxEx+AxyEy),
(2)
2Eyz2=iωμ0(AyxEx+AyyEy),

where

(3)
Axx=σxx+iωε-σxzσzxσzz+iωε,
(4)
Axy=σxy+iωε-σxzσzyσzz+iωε,
(5)
Ayx=σyx-σzxσyzσzz+iωε=Axy,
(6)
Ayy=σyy+iωε-σyzσzyσzz+iωε,

The electric fields for x and y directions (Ex and Ey) can be calculated using the equations (1), (2), (3), (4), (5), (6) and the electric field for the z direction is defined as Ez=-σzxEx+σzyEyσzz+iωε (Pek and Santos, 2002). The parameters, 𝜔, μ0 and 𝜀, represent the angular frequency, magnetic permeability in free space and electric permittivity, respectively (μ=μ0 for EM surveys including MT). The anisotropic media is presented by the conductivity tensor as follows:

(7)
σ~=σxxσxyσxzσyxσyyσyzσzxσzyσzz=Rz(αS)Rx(αD)Rz(αL)σ1000σ2000σ3Rz(αL)Rx(αD)Rz(αS),

and the conductivity tensor is symmetric. Thus, Axy and Ayx in equations (4) and (5) are equal. From equations (1), (2), (3), (4), (5), (6), it can be inferred that electric fields of the plane-wave are affected only by the six independent variables of the conductivity tensor, which can be simplified to three independent variables of A (Axx, Axy=Ayx and Ayy) (Pek and Santos, 2002).

Analytic Solutions using the Matrix Propagation Method

EM fields in a layered earth can be computed using a matrix propagation method which utilizes a propagation matrix to relate the electric and magnetic field at the top and bottom of each layer. This matrix is combined iteratively throughout the entire layered model to form a total propagation matrix. The matrix propagation method is required for computing EM field at any depth in any layer of an 1-D medium.

Analytic Solution for 1-D Anisotropic Media

Since the general solutions of equations (1) and (2) can be expressed in the form of exp(±kz), it is easily noted that there exist two solutions corresponding to the two different wavenumbers (k1,2):

(8)
k1,22=iωμ02(Axx+Ayy)±(Axx-Ayy)2+4AxyAyx).

In a general anisotropic layered medium, where Axy0, the x- and y-components of the electric fields (Ex and Ey), are not independent. The relationship between Ex and Ey can be expressed as follows:

(9)
Ex=exp(±kz),
(10)
Ey=qexp(±kz)=qEx,

where the coefficient q can be computed by substituting equations (9) and (10) into equations (1) and (2) as follows:

(11)
q1,2=k1,22-iωμ0Axxiωμ0Axy.

The general solution of Ex satisfying equations (1) and (2) can be expressed as follows:

(12)
Ex=d1+X(+1)+d1-X(-1)+d2+X(+2)+d2-X(-2),

where X(±j)=exp(±kjz), d1+ and d2+ are constants for two upgoing wave components, and d1- and d2- are constants for two downgoing wave components. Using the coupling effect of Ex and Ey, the system of equations for Ex, Ey, Hx and Hy can be defined as follows:

(13)
ExEyHxHy=X(+1)X(1)X(+2)X(2)q1X(+1)q1X(1)q2X(+2)q2X(2)γq1k1X(+1)γq1k1X(1)γq2k2X(+2)γq2k2X(2)γk1X(+1)γk1X(1)γk2X(+2)γk2X(2)d1+d1d2+d2,
(14)
f(z,ω)=M(z,ω)d,

where γ=-(iωμ0)-1. As the elements of M(z,ω) can be obtained from the anisotropic layered model, d of each layer remains the only variable to be solved to compute the components of the EM fields. Because of the continuity condition for the horizontal EM fields, those at zl+1 boundary computed at the lth and (l + 1)th layers must be of the same value as follows:

(15)
fl(zl+1,ω)=fl+1(zl+1,ω).

Thus, the relationship between dl for the lth layer and dl+1 for the (l + 1)th layer can be written by

(16)
dl=Ml-1(zl+1,ω)Ml+1(zl+1,ω)dl+1.

From equations (15) and (16), the EM fields at the top boundary of the air layer (0th layer) to which the boundary condition applies can be expressed by dN of the bottom layer:

(17)
f0(z0,ω)=M0(z0,ω)d0=M0(z0,ω)M01(z1,ω)M1(z1,ω)MN11(zN,ω)MN(zN,ω)dN=MTotdN.

Since only the downgoing wave components exist in the last layer (infinite half-layer), d1,N+ and d2,N+ should be zero. As the boundary condition at the top boundary of the air layer applies Esource=(Ex,source,Ey,source,Ez,source) that is usually E1,source = (1,0,0) or E2,source = (0,1,0), and dN is composed of only d1,N- and d2,N-, the 4 × 4 matrix in equation (17) can be reduced to the 2 × 2 matrix as follows:

(18)
Ex,sourceEy,source=MTot(1,2)MTot(1,4)MTot(2,2)MTot(2,4)d1,Nd2,N.

As Ex,source, Ey,source and the elements of MTot are known, dN can be computed using equation (18).

By using equation (16) and dN, it is now possible to compute d at all layers of the anisotropic model. With M(z,ω) defined from the layered model and the calculated vector d for each layer, the components of EM waves can now be computed from equation (14).

The EM field for isotropic medium can be represented in a similar manner to equation (13) as follows:

(19)
ExEyHxHy=X(+1)X(1)0000X(+2)X(2)00γkyX(+2)γkyX(2)γkxX(+1)γkxX(1)00d1+d1d2+d2+,

where X(±1)=exp(±k1z)=exp(±kxz);X(±2)=exp(±k2z)=exp(±kyz);kx2=iωμ0Axx;andky2=iωμ0Ayy.

Stabilization strategies for numerical computation of 1-D solutions

Stabilization strategy I

Stability strategy I (SS-I) is based on the open-source code ‘MARE2DEM’ (Key, 2016), which modify the exponential terms to numerically stable form to compute the electric and magnetic fields at any profound depths of layered anisotropic models. According to equations (13) and (19), the positive exponential terms of M(z,ω), X(j)=exp(kjz) approach infinity as z increases drastically. As the positive exponential terms approach infinity, numerical instabilities, i.e. NaN values, may occur when computing the inverse of M(z,ω) according to equations (16) and (17). To avoid such numerical instabilities, SS-I modifies the exponential terms X(±j)=exp(±kjz) to X'(-j)=exp(-kj(z-zi)) for downgoing wave components and X'(+j)=exp(kj(z-zi+1)) for upgoing wave components at the ith layer. By applying this strategy to equations (13) and (19), all the terms within the exponential functions are always negative and thus, become numerically stable exponential terms which converges to zero as z increases excessively. We can obtain two new matrices for anisotropic and isotropic media as follows:

(20)
ExEyHxHy=X'(+1)X'(-1)X'(+2)X'(-2)Q1X'(+1)Q1X'(-1)Q2X'(+2)Q2X'(-2)-γQ1k1X'(+1)γQ1k1X'(-1)-γQ2k2X'(+2)γQ2k2X'(-2)γk1X'(+1)-γk1X'(-1)γk2X'(+2)-γk2X'(-2)d1'+d1'-d2'+d2'-,
(21)
ExEyHxHy=X'(+1)X'(-1)0000X'(+2)X'(-2)00-γkyX'(+2)γkyX'(-2)γkxX'(+1)-γkxX'(-1)00d1'+d1'-d2'+d2'-.

Equations (20) and (21) can be expressed in a simplified manner as follows:

(22)
f(z,ω)=M'(z,ω)d',

where M(f,ω) comprises the components of EM field, M'(z,ω) is a matrix composed of the modified exponential terms and is a vector of the corresponding upgoing and downgoing wave components. By using the continuity solution, solving the components of the EM fields of equation (22) are identical to obtaining the analytic solution for 1-D anisotropic media as shown in equations (15), (16), (17). Fig. 1(a) shows how the term X'(-j) is computed for downgoing wave components; and Fig. 1(b) represents how X'(+j) for upgoing wave components is computed. As shown in Fig. 1(a) and 1(b), the terms inside the exponential functions in M'(z,ω) are always defined by negative values, which improve the stability of the computational process. By defining the exponents by negative values, the terms X'(±j) converge to zero [X'(±j)=exp(-)0 as z]. This modification allows us to stably compute the EM fields and prevents numerical overflow at large depths.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F1.jpg
Fig. 1.

Computation of general solutions of M'(z,ω) in equations (20) and (21) at an arbitrary depth (z) in the ith layer with the upper (zi) and lower (zi+1) boundaries for (a) downgoing wave components X'(-j) and (b) upgoing wave components X'(+j).

Stabilization Strategy II

As shown in equations (16) and (17), computation of the inverse matrix Mi'-1(zi+1,ω) where i varies from 0 to N-1 is necessary for solving the components of the EM fields. The components of Mi'-1(zi+1,ω) are defined as X'(j)=exp(kj(zi+1-zi+1))=exp(0)=1 and X'(-j)=exp(-kj(zi+1-zi)). When the thickness of the ith layer (zi+1-zi) is too thick, the term of X'(-j) converges to zero. The problem associated with this phenomenon is that, according to equations (20) and (21), the elements of the second and fourth columns of Mi'(zi+1,ω) also become zero. Thus, the inverse matrix of Mi'(zi+1,ω) may not be computed, for which the computation of the electric and magnetic fields for a layered model with a thick layer can be numerically unstable.

To overcome this limitation, SS-II first computes the determinant of Mi'(zi+1,ω) for each layer of a given model before computing M'(z,ω) and . When the determinant of a certain thick layer is less than a given tolerance, the underlying layers do not have a significant effect on the EM fields, which means that the underlying layers can be neglected. Therefore, we can replace the thick layer with an infinite half-space (Fig. 2), which allows us to avoid the unstable calculation of the inverse matrix of Mi'(zi+1,ω) for the thick layer. Through this procedure, we can stably compute the EM fields for a given model with a thick layer. In particular, the stability condition of this problem changes depending on the frequency of the EM fields and the components of the conductivity tensor of each layer. However, because the elements of Mi'(zi+1,ω) take all the values for the frequency and the conductivity tensor into account, SS-II based on the determinant of the matrix is highly applicable and versatile.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F2.jpg
Fig. 2.

A four-layered model including the air layer: (a) the second layer is a thick layer; and (b) the second layer is replaced with an infinite half-space.

Fig. 2(a) shows a 4-layered model that we apply SS-II to. It is assumed that the second layer (Layer 2) is a thick layer and the last layer (Layer 4) is the infinite half-space. To compute the EM fields of the model shown in Fig. 2(a), we first compute the determinant of Mi'(zi+1,ω) for each layer of the model excluding the bottom layer according to SS-II. Because the determinant for the second layer of Fig. 2(a), which is too thick, is smaller than the given tolerance, we replace it with an infinite half-layer as shown in Fig. 2(b).

Numerical Experiments

To demonstrate the accuracy and effectiveness of our algorithm, we apply the strategies to another 4-layered anisotropic model.

Verification of the two proposed strategies

In the 4-layered anisotropic model shown in Fig. 3, the top layer and bottom half-space are assumed to be isotropic with resistivity of 1000 Ω·m; and the second and third layers are anisotropic with the principal resistivities of 400/800/400 Ω·m and 800/400/800 Ω·m with the Euler’s angles of 20°/30°/10° and 10°/20°/30°, respectively. The resistivity of the air layer is assigned to be 109 Ω·m and we use a total of 21 frequencies sampled at a uniform interval on the logarithmic scale ranging from 0.01 to 100 Hz in the computational process.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F3.jpg
Fig. 3.

A four-layered 1-D anisotropic model.

In Fig. 4, we compare the MT responses computed at the surface using our algorithm with those obtained from the code provided by Pek and Santos (2002). Fig. 4(a) and 4(b) show that the xy- and yx-components of the apparent resistivity and phase computed by our algorithm are highly accurate and similar to those obtained by the 1-D code of Pek and Santos (2002). This is further supported by the relative errors and phase differences displayed in Fig. 4(c) and 4(d). According to Fig. 4(c), the relative errors in apparent resistivity between our algorithm and Pek and Santos (2002) algorithms are almost negligible, approximately 0.004% for the xy-component and 0.003% for the yx-component. In Fig. 4(d), the phase differences are approximately 0.005° for both the xy- and yx-components. These results support for the accuracy of our algorithm.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F4.jpg
Fig. 4.

The xy- (the orange circles) and yx-components (the orange triangles) for (a) the apparent resistivity and (b) phase computed by our algorithm (denoted by “Lee”) with those obtained by Pek and Santos (2002) algorithm (the solid yellow line; denoted by “Pek”). (c) and (d) show the relative errors of the apparent resistivity and the phase differences between our and Pek and Santos’ algorithms for the xy- (the solid yellow line with points) and yx-components (the solid orange line with points), respectively.

Numerical Experiments

We compare the real x-component of the electric fields calculated using our algorithms with and without SS-I and SS-II to examine their effects. For SS-I, we compare the electric fields at receivers with different depths to investigate if SS-I stabilizes the numerical computation of the electric field, specifically at deep depths. For SS-II, we compare the electric fields at the surface varying the thickness of the first layer in the layered model, from which we can examine if SS-II resolves numerical instability at large thickness of the first layer. All the calculations are conducted with our algorithm based on 4-bytes to efficiently visualize the computational results.

Stabilization Strategy I

As mentioned in the subsection ‘Stabilizing Strategy I’, SS-I was designed to stably calculate the EM fields at any great depths of layered models. To confirm the effectiveness of SS-I, we compare the real part of x-component of the electric fields (Ex) obtained by the conventional matrix M(z,ω) (equations 13 and 19) and modified matrix M'(z,ω) (equations 20 and 21) according to the depth of receivers (0 ~ 10,000 km). For the numerical experiments, we use the layered model shown in Fig. 4 and a single frequency of 10 Hz.

Fig. 5 shows the responses of Ex computed using the conventional method and SS-I when the plane wave sources, E1,source = (1,0,0) (Fig. 5(a)) and E2,source = (0,1,0) (Fig. 5(b)), are applied. In Fig. 5(a) and 5(b), the electric fields obtained by the conventional method and SS-I are nearly identical up to a certain depth (about 447.6 km). However, the conventional method produces NaN values for Ex for both E1,source and E2,source below a depth of 447.6 km, which is denoted by the orange curly bracket. In contrast, the algorithm with SS-I stably compute the electric fields, so that they converge to zero even below a depth of 447.6 km.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F5.jpg
Fig. 5.

The real x-component of the electric field (Ex) for (a) E1,source= (1,0,0) and (b) E2,source= (0,1,0) with respect to the depth of receivers (0 ~ 10,000 km). The results computed by the conventional method and SS-I are displayed with the orange markers and yellow solid lines, respectively. The part displayed by the orange curly brackets indicates the depths with NaN values for the conventional method.

Among the terms X(±j) and X'(±j) which are the general solutions required to compute Ex with respect to the depth of the receiver for the conventional method and SS-I, Figs. 6 and 7 shows the real and imaginary components of the term X(±j) up to a depth of 10 km. It is observed that the terms X(±j) and X'(±j) show continuous results within each layer, but discontinuous results at the boundaries of layers. Because layer 1 (0~2 km) and layer 4 (7.5 km ~) are isotropic, the terms X(±2) and X'(±2) are zero (see Figs. 6(c) and 7(c)). The terms X(±j) and X'(±j) computed with both the conventional method and SS-I do not suffer from numerical instabilities at shallow depths. As z increases (z), the terms X(-1), X(-2) and X'(±j) converge to zero as exp(-)=0. In contrast, because X(+1), X(+2)=exp(-), numerical instabilities may occur.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F6.jpg
Fig. 6.

(a)-(d): The real components and (e)-(h): imaginary components of X(±j) computed using the conventional method up to a receiver depth of 10 km. The dim orange dotted vertical lines represent the depths of layer interfaces.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F7.jpg
Fig. 7.

(a)-(d): The real components and (e)-(h): imaginary components of X'(±j) computed using the conventional method up to a receiver depth of 10 km. The dim orange dotted vertical lines represent the depths of layer interfaces.

Tables 1 and 2 show the results of d and at each layer whose components are the coefficients required to compute Ex. In Tables 1 and 2, there are no numerical instabilities of d and at each layer. Because upgoing wave components at the last layer do not exist, the values for components d1+, d1'+, d2+ and d2'+ are all zero. Even if X(+1) has large values at depths shallower than 447.6 km, Ex is not affected because d1+ is zero, and thus d1+X(+1)=0. However, below a depth of 447.6 km, X(+1) becomes infinity, resulting in numerical instability in the computation of Ex. Table 2 shows that SS-I does not suffer from numerical instability in computing .

Table 1.

Vector d calculated with the conventional method

Layer Vector d for the conventional method
d1+d1-d2+d2-
1 ‒1.945×10-5
‒2.968×10-4 i
2.462×10-3
2.241×10-3 i
‒4.099×10-5
‒8.213×10-5 i
‒4.118×10-5
‒8.151×10-5 i
2 ‒1.226×10-5
‒9.726×10-6 i
2.029×10-3
1.125×10-3 i
1.215×10-5
‒4.931×10-6 i
7.278×10-4
5.818×10-4 i
3 7.633×10-7
‒3.738×10-6 i
1.574×10-3
6.570×10-4 i
‒9.380×10-7
‒2.344×10-6 i
6.858×10-4
1.067×10-3 i
4 0
0 i
8.090×10-4
1.911×10-3 i
0
0 i
‒4.096×10-5
‒6.756×10-5 i
Table 2.

Vector calculated with SS-I

Layer Vector for SS-I
d1'+d1'-d2'+d2'-
1 ‒1.990×10-4
3.989×10-4 i
2.480×10-3
‒2.258×10-3 i
‒1.043×10-4
8.974×10-5 i
‒4.148×10-5
8.212×10-5 i
2 ‒4.714×10-5
‒5.939×10-5 i
5.267×10-4
‒1.130×10-3 i
2.365×10-8
4.321×10-5 i
2.396×10-4
‒5.330×10-4 i
3 ‒3.252×10-5
‒1.722×10-5 i
‒1.238×10-4
‒3.603×10-4 i
‒1.196×10-5
‒6.165×10-6 i
‒2.162×10-4
‒3.610×10-4 i
4 0
0 i
‒4.177×10-4
‒2.180×10-4 i
0
0 i
1.454×10-5
1.050×10-5

Additionally, Tables 1 and 2 show that the components of d and are not affected by the depth of the receiver in a given 1-D model. In contrast, the components of X(±j) and X'(±j) are affected by the depth of the receiver, as depicted in Figs. 6 and 7. From these numerical experiments, we can conclude that the components of X(+j) degrade the stability in the computation of Ex using the conventional method at large depths, whereas the components of X'(±j) are computed in a stabilized manner through SS-I, allowing stable computation of the EM field at large depths, as shown in Fig. 5.

Stabilization Strategy II

As mentioned in subsection ‘Stabilization Strategy II’, the presence of a thick layer in a layered model can lead to unstable or inaccurate results in the computation of the electric and magnetic fields. To verify the effectiveness of the proposed algorithm with SS-II, a numerical experiment was conducted. The thickness of the first layer of the model shown in Fig. 3 varies from 2 to 1002 km in a uniform increment of 0.1 km. A receiver is located at the surface and a single frequency of 10 Hz is used. The tolerance was empirically determined, and from several numerical experiments, we recommend the tolerances of 10-27 and 10-100 for the algorithms with 4-bytes and 8-bytes, respectively.

Fig. 8 shows the real x-component of the electric field (Ex) computed without and with SS-II for both E1,source= (1,0,0) (Fig. 8(a)) and E2,source = (0,1,0) (Fig. 8(b)), when the thickness of the first layer in Fig. 3 increases from 2 to 1002 km. In Fig. 8(a) for E1,source, Ex computed with SS-II converges an approximate value of 2.5 × 10-3 V/m whereas Ex computed without SS-II produces NaN value beginning at a thickness of 187.5 km which is denoted with the orange curly brackets. In Fig. 8(b), when the thickness exceeds 187.5 km, Ex computed with SS-II for E2,source converges to zero, whereas Ex computed with the conventional method yields NaN values. As mentioned in the subsection of ‘Stabilizing Strategy II’, the corresponding strategy determines whether the determinant of the first layer (det M1'(z2,ω)) is smaller than the given tolerance. If the thickness of the first layer exceeds 166.4 km for E1,source, which means that the determinant of M1'(z2,ω) is less than the tolerance, the first layer is replaced with a new infinite half-space, which enables numerically stable computation of the EM field.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F8.jpg
Fig. 8.

The real x-component of the electric field (Ex) computed for (a) E1,source= (1,0,0) and (b) E2,source = (0,1,0) according to the thickness of the first layer in Figure 4 (from 2 to 1002 km) obtained by the algorithms without (the orange points) and with (the solid yellow lines) SS-II. The thicknesses with NaN values in (a) and (b) are denoted with the orange curly brackets.

We now discuss the factors that may cause instability for a thick layer in the computation process of Ex and how SS-II stabilizes the computational process. Because the first layer shown in Fig. 3 is isotropic, Ex at the surface can be defined as Ex=X1'+d1'++X1'-d1'-+0*d2'++0*d2'- according to equation (21). X'(+1) can be defined as exp(k1(-z2)) at the surface (z=0) and as z2 increases, X'(+1) converges to zero, and the term X'(-1) at z=0 has a constant value of 1 as exp(-k1(0-0))=exp(0)=1.

The components of for each layer can be computed using the continuity condition, which require calculating the inverse of M1'(z2,ω). M1'(z2,ω) takes the form as shown in equation (21) because the first layer is isotropic. Fig. 9 shows the real and imaginary components of X'(+1) and X'(-1) of M1'(z2,ω). Fig. 10 shows the determinant of M1'(z2,ω). Because the first layer is isotropic (i.e., k1=kx=ky), the components X'(+1)=exp(k1(z2-z2))=exp(0)=1 as shown in Fig. 9(a) and 9(c). However, the components X'(-1)=exp(-k1(z2-z1)) converges to zero, when z2 continues to increase, as shown in Fig. 9(b) and 9(d). Similarly, X'(-2) in equation (21) converges to zero as z2 increases. Consequently, as the thickness of the first layer increases, the elements of the second and fourth columns of M1'(z2,ω) become zero according to equation (21) and thus the determinant converges to zero as shown in Fig. 10. Thus, the computation of the inverse of M1'(z2,ω) becomes numerically unstable, which also result in unstable computation of of each layer.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F9.jpg
Fig. 9.

The real components of (a) X1'+ and (b) X1'+and imaginary components of (c) X1'+ and (d) X1'+ for M1'(z2,ω) computed without SS-II denoted with the orange markers and with SS-II denoted with the solid yellow lines.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F10.jpg
Fig. 10.

The determinant of M1'(z2,ω) with respect to the thickness of the first layer ranging from 2 to 1002 km.

Fig. 11 shows the real components of computed without and with SS-II at the surface, when the thickness of the first layer varies from 2 to 1002 km. As mentioned before, the components d1'+, d1'-, d2'+ and d2'- produce NaN values at thicknesses larger than 187.4 km when computed without SS-II (refer to the curly orange brackets in Fig. 11). As Ex=X1'+d1'++X1'-d1'-+0*d2'++0*d2'- in the first layer, it results in unstable computation of Ex starting at a thickness of 187.5 km as shown in Fig. 8(a) and 8(b). With the implementation of SS-II, when the thickness of the first layer exceeds 166.4 km, the corresponding layer can now be regarded as an infinite half-space. As a result, based on equation (24), the computation of does not require the inverse matrix of M1'(z2,ω), preventing numerical instabilities as the components d1'+, d1'-, d2'+ and d2'- converge to ‒0.831, 1.371, 0 and 0, respectively. Thus, SS-II allows numerically stable computation of the EM field at presence of a thick layer.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F11.jpg
Fig. 11.

The real components of (a) d1'+, (b) d1'-, (c) d2'+ and (d) d2'- of vector computed without SS-II denoted with the orange markers and with SS-II denoted with the solid yellow lines. The section where NaN values appear is denoted with the orange curly brackets.

3-D Model II (3D homogeneous model for initial model of inversion)

In 3-D MT inversion, where both low and high frequencies must be considered, accurate computation of EM fields at low frequencies need increasing size of the mesh with depth, which may cause numerical errors in forward modeling response at high frequencies. To prevent arising of such numerical instabilities, our suggested strategies are of great importance. An extremely large 3-D homogeneous isotropic model with a background resistivity of 100 Ω·m and a computational domain of 6,000 km × 6,000 km × 6,000 km as shown in Fig. 12, was used to illustrate the effectiveness and accuracy of SS-I and SS-II

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F12.jpg
Fig. 12.

(a) 3-D homogeneous isotropic model with a computational domain of 6,000 km × 6,000 km × 6,000 km, excluding the air layer and (b) its y-z cross-section.

Fig. 13(a) and 13(b) show the computed xy-components of forward modeling responses, i.e., apparent resistivities and phase, and Fig. 13(c) and 13(d) show the relative error in computed forward modeling responses between the open-source code (Bai et al., 2022) and the proposed algorithm for a total of 17 frequencies ranging from 0.0001 Hz to 1 Hz. As shown in Fig. 13, the forward modeling responses computed with the proposed algorithm and the open-source code show stable results at low frequencies. However, as shown in Fig. 13, the forward modeling responses computed using the open-source code start to diverge and show numerical instabilities at higher frequencies, which is denoted by the orange curly brackets. In contrast, the forward modeling responses computed with the proposed algorithm, denoted by the solid yellow line, show stable results even at higher frequencies. From Fig. 13(a) and 13(b), it is evident that the forward modeling responses computed with the open-source code begin to show numerical instabilities starting from a frequency of 0.1 Hz. In contrast, the proposed algorithm shows stable computation of the forward modeling responses even beyond a frequency of 0.1 Hz. The two forward modeling responses match well, with a maximum apparent resistivity error of 0.71% and a maximum phase error of 0.13° for both xy- and yx-components, respectively. From numerical experiments, it can be seen that SS-I and SS-II allow stable computation and prevent numerical instabilities even in an extremely large 3-D model, further verifying the effectiveness of the proposed algorithm.

https://cdn.apub.kr/journalsite/sites/ksmer/2024-061-04/N0330610402/images/ksmer_61_04_02_F13.jpg
Fig. 13.

Components of (a) ρxy, (b) ϕxy computed using the proposed algorithm denoted with the yellow solid line and open-source code (Bai et al., 2022) denoted with the orange markers along the dotted line. (c) and (d) represent the relative error and difference of xy-component of apparent resistivity and phase between the proposed algorithm and open-source code, respectively. The section where NaN values appear is denoted with the orange curly brackets.

Conclusions

For a stable computation of various MT responses and EM fields at an arbitrary depth of 1D layered models including both isotropic and anisotropic layers, we proposed the two strategies. The first strategy involves converting the general solutions of the second-order differential equations into computationally stable forms to prevent the components of EM field from diverging at large depths. The second strategy assesses whether each layer of the layered model is thick enough to cause instability. Through the numerical experiments on the synthetic model, we confirmed that our algorithm stably computes the EM waves even at considerable depths and yields stable and accurate results even in the presence of an extremely thick layer. The developed algorithm designed for 1-D anisotropic media can be adapted to satisfy the boundary conditions in multi-dimensional (2-D and 3-D) MT or EM modeling.

Acknowledgements

This work was supported by the Institute for Korea Spent Nuclear Fuel (iKSNF) (No. 2021040101003A) and Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korea government (Ministry of Trade, Industry and Energy (MOTIE)). (No. 20226A10100030, Development of high-performance monitoring technology for marine CO2 storage). The authors would like to thank Bai et al. (2022) for the open-source code which is used in this paper for comparison of forward modeling responses. The authors would also like to thank Pek and Santos (2002) for providing the 1-D anisotropic modeling code.

References

1

Bai, N., Zhou, J., Hu, X., and Han, B., 2022. 3D edge-based and nodal finite element modeling of magnetotelluric in general anisotropic media, Computers & Geosciences, 158, 104975.

10.1016/j.cageo.2021.104975
2

Cao, X., Huang, X., Yin, C., Yan, L., and Zhang, B., 2021. 3D MT anisotropic inversion based on unstructured finite-element method, Journal of Environmental and Engineering Geophysics, 26(1), p.49-60.

10.32389/JEEG20-006
3

Key, K., 2016. MARE2DEM: a 2-D inversion code for controlled-source electromagnetic and magnetotelluric data, Geophysical Journal International, 207(1), p.571-588.

10.1093/gji/ggw290
4

Li, Y. and Pek, J., 2008. Adaptive finite element modeling of two-dimensional magnetotelluric fields in general anisotropic media, Geophysical Journal International, 175(3), p.942-954.

10.1111/j.1365-246X.2008.03955.x
5

Liu, Y., Xu, Z., and Li, Y., 2018. Adaptive finite element modeling of three-dimensional magnetotelluric fields in general anisotropic media, Journal of Applied Geophysics, 151, p.113-124.

10.1016/j.jappgeo.2018.01.012
6

O'Brien, D.P. and Morrison, H.F., 1967. Electromagnetic fields in an n-layer anisotropic half-space, Geophysics, 32(4), p.668-677.

10.1190/1.1439882
7

Pek, J. and Santos, F.A., 2002. Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media, Computers & Geosciences, 28(8), p.939-950.

10.1016/S0098-3004(02)00014-6
8

Reddy, I.K. and Rankin, D., 1971. Magnetotelluric effect of dipping anisotropies, Geophysical Prospecting, 19(1), p.84-97.

10.1111/j.1365-2478.1971.tb00586.x
9

Rong, Z., Liu, Y., Yin, C., Wang, L., Ma, X., Qiu, C., Zhang, B., Ren, X., Su, Y., and Weng, A., 2022. Three-dimensional magnetotelluric inversion for arbitrarily anisotropic earth using unstructured tetrahedral discretization, Journal of Geophysical Research: Solid Earth, 127(8), e2021JB023778.

10.1029/2021JB023778
10

Usui, Y., 2015. 3-D inversion of magnetotelluric data using unstructured tetrahedral elements: applicability to data affected by topography, Geophysical Journal International, 202(2), p.828-849.

10.1093/gji/ggv186
11

Usui, Y., Kasaya, T., Ogawa, Y., and Iwamoto, H., 2018. Marine magnetotelluric inversion with an unstructured tetrahedral mesh, Geophysical Journal International, 214(2), p.952-974.

10.1093/gji/ggy171
12

Wannamaker, P.E., 2005. Anisotropy versus heterogeneity in continental solid earth electromagnetic studies: fundamental response characteristics and implications for physicochemical state, Surveys in Geophysics, 26, p.733-765.

10.1007/s10712-005-1832-1
13

Xiao, T., Liu, Y., Wang, Y., and Fu, L.Y., 2018. Three-dimensional magnetotelluric modeling in anisotropic media using edge-based finite element method, Journal of Applied Geophysics, 149, p.1-9.

10.1016/j.jappgeo.2017.12.009
페이지 상단으로 이동하기