Stable Calculation of Responses for a Plane Electromagnetic Wave in an Anisotropic Layered Media for the Modeling of Magnetotelluric Data

Hee June Lee; Janghwan Uhm; Junyeong Heo; Dong-Joo Min

doi:10.32390/ksmer.2024.61.4.284

Preview

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

Stable Calculation of Responses for a Plane Electromagnetic Wave in an Anisotropic Layered Media for the Modeling of Magnetotelluric Data

이방성 층서에서의 자기지전류 탐사 자료 모델링을 위한 평면 전자기파에 대한 반응 계산 안정화

Hee June Lee¹

Janghwan Uhm²^*

Junyeong Heo³

Dong-Joo Min¹⁴

이 희준¹

엄 장환²^*

허 준영³

민 동주¹⁴

¹Department of Energy Systems Engineering, Seoul National University, Seoul, Korea

²Deep Subsurface Storage and Disposal Research Center, Geology and Space Research Division, Korea Institute of Geoscience and Mineral Resources (KIGAM), Daejeon, Korea

³Mineral Exploration and Mining Research Center Mineral Resources Division, Korea Institute of Geoscience and Mineral Resources (KIGAM), Daejeon, Korea

⁴Research Institute of Energy and Resources, Seoul National University, Seoul, Korea

¹서울대학교 에너지시스템공학부

²한국지질자원연구원 심층처분환경연구센터

³한국지질자원연구원 자원탐사개발연구센터

⁴서울대학교 에너지자원신기술연구소

^{*Corresponding Author}

License (open-access, https://creativecommons.org/licenses/by-nc/4.0/):

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

Research in multi-dimensional magnetotelluric (MT) modeling and inversion, particularly within anisotropic media, is gaining substantial interest. One-dimensional (1-D) modeling is utilized as boundary conditions necessary for accurate 2-D or 3-D anisotropic modeling. Based on Maxwell’s equations and interface boundary conditions, we developed a new algorithm to calculate electromagfor both 1-D isotropic and anisotropic layered models. However, this algorithm encounters stability issues at considerable depths or within exceptionally thick layers. To address these challenges, we introduced two stabilization strategies. The first strategy alters the exponential terms in the plane-wave solutions to mitigate instabilities, while the second evaluates each layer’s contribution to numerical instability, replacing problematic thick layers with an infinite half-space. Through numerical tests on a four-layer anisotropic model, these strategies proved effective, ensuring stable computations of electromagnetic fields at any depth and in configurations with thick layers, thereby enhancing the model’s applicability for in-depth anisotropic analyses.

Keywords

magnetotelluric

1-D anisotropic layered model

stabilization strategies

electromagnetic fields

modeling

본 연구에서는 자기지전류(Magnetotelluric, MT) 탐사 자료 모델링을 위한 1차원 이방성 매질에서의 평면파 반응 계산의 안정화 전략을 제안한다. 층서구조에 대한 반응, 2차원이나 3차원 모델링에서 경계조건으로 이용되기 때문에, 이방성 매질에 대한 고차원 모델링의 정확한 해를 구하기 위해서는 이방성 층서 구조에 대한 정확하고 안정적인 해의 계산이 필수적이다. 고심도 혹은 매우 두꺼운 층들이 존재할 때 이방성 층서 구조 해의 불안정성 문제를 극복하기 위해 두 가지 안정화 전략을 제안한다. 즉, 평면파 해의 지수 항을 수식적으로 깊이가 증가함에 따른 발산 현상을 방지하거나, 층서 모델에서 불안정성을 일으키는 두꺼운 층이 존재할 때 이 층을 반무한 매질으로 치환해 수치적 안정성을 향상시킨다. 두 안정화 전략은 임의의 깊이와 두꺼운 층이 존재하는 경우에도 안정적으로 전자기장 성분들을 구할 수 있어, 제안된 알고리즘은 이방성 층서 구조의 임의의 깊이에서 1차원 해를 계산할 때 활용될 수 있을 것으로 기대된다.

키워드

자기지전류

1-D 이방성 층서 구조

안정화 전략

전자기장

모델링

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., $\frac{\partial}{\partial x} = \frac{\partial}{\partial y} = 0$ , the corresponding equations can be defined as follows (Pek and Santos, 2002):

(1)

\frac{\partial^{2} E_{x}}{\partial z^{2}} = i ω μ_{0} (A_{x x} E_{x} + A_{x y} E_{y}),

(2)

\frac{\partial^{2} E_{y}}{\partial z^{2}} = i ω μ_{0} (A_{y x} E_{x} + A_{y y} E_{y}),

where

(3)

A_{x x} = σ_{x x} + i ω ε - \frac{σ_{x z} σ_{z x}}{σ_{z z} + i ω ε},

(4)

A_{x y} = σ_{x y} + i ω ε - \frac{σ_{x z} σ_{z y}}{σ_{z z} + i ω ε},

(5)

A_{y x} = σ_{y x} - \frac{σ_{z x} σ_{y z}}{σ_{z z} + i ω ε} = A_{x y},

(6)

A_{y y} = σ_{y y} + i ω ε - \frac{σ_{y z} σ_{z y}}{σ_{z z} + i ω ε},

The electric fields for x and y directions ( $E_{x}$ and $E_{y}$ ) can be calculated using the equations (1), (2), (3), (4), (5), (6) and the electric field for the z direction is defined as $E_{z} = - \frac{σ_{z x} E_{x} + σ_{z y} E_{y}}{σ_{z z} + 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)

\begin{aligned} \tilde{σ} = (\begin{array}{ccc} σ_{x x} & σ_{x y} & σ_{x z} \\ σ_{y x} & σ_{y y} & σ_{y z} \\ σ_{z x} & σ_{z y} & σ_{z z} \end{array}) = \\ R_{z} (- α_{S}) R_{x} (- α_{D}) R_{z} (- α_{L}) (\begin{array}{cccc} σ_{1} & 0 & 0 \\ 0 & σ_{2} & 0 \\ 0 & 0 & σ_{3} \end{array}) R_{z} (α_{L}) R_{x} (α_{D}) R_{z} (α_{S}), \end{aligned}

and the conductivity tensor is symmetric. Thus, $A_{x y}$ and $A_{y x}$ 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 ( $A_{x x}$ , $A_{x y} = A_{y x}$ and $A_{y y}$ ) (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 $e x p (\pm k z)$ , it is easily noted that there exist two solutions corresponding to the two different wavenumbers ( $k_{1, 2}$ ):

(8)

k_{1, 2}^{2} = \frac{i ω μ_{0}}{2} \{(A_{x x} + A_{y y}) \pm \sqrt{(A_{x x - A_{y y}})^{2} + 4 A_{x y} A_{y x})}\} .

In a general anisotropic layered medium, where $A_{x y} \neq 0$ , the x- and y-components of the electric fields ( $E_{x}$ and $E_{y}$ ), are not independent. The relationship between $E_{x}$ and $E_{y}$ can be expressed as follows:

(9)

E_{x} = e x p (\pm k z),

(10)

E_{y} = q e x p (\pm k z) = q E_{x},

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

(11)

q_{1, 2} = \frac{k_{1, 2}^{2} - i ω μ_{0} A_{x x}}{i ω μ_{0} A_{x y}} .

The general solution of $E_{x}$ satisfying equations (1) and (2) can be expressed as follows:

(12)

E_{x} = d_{1}^{+} X (+ 1) + d_{1}^{-} X (- 1) + d_{2}^{+} X (+ 2) + d_{2}^{-} X (- 2),

where $X (\pm j) = e x p (\pm k_{j} z)$ , $d_{1}^{+}$ and $d_{2}^{+}$ are constants for two upgoing wave components, and $d_{1}^{-}$ and $d_{2}^{-}$ are constants for two downgoing wave components. Using the coupling effect of $E_{x}$ and $E_{y}$ , the system of equations for $E_{x}$ , $E_{y}$ , $H_{x}$ and $H_{y}$ can be defined as follows:

(13)

\begin{aligned} (\begin{array}{c} E_{x} \\ E_{y} \\ H_{x} \\ H_{y} \end{array}) = \\ (\begin{array}{cccc} X (+ 1) & X (- 1) & X (+ 2) & X (- 2) \\ q_{1} X (+ 1) & q_{1} X (- 1) & q_{2} X (+ 2) & q_{2} X (- 2) \\ - γ q_{1} k_{1} X (+ 1) & γ q_{1} k_{1} X (- 1) & - γ q_{2} k_{2} X (+ 2) & γ q_{2} k_{2} X (- 2) \\ γ k_{1} X (+ 1) & - γ k_{1} X (- 1) & γ k_{2} X (+ 2) & - γ k_{2} X (- 2) \end{array}) (\begin{matrix} d_{1}^{+} \\ d_{1}^{-} \\ d_{2}^{+} \\ d_{2}^{-} \end{matrix}), \end{aligned}

(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 $z_{l + 1}$ boundary computed at the l^th and (l + 1)^th layers must be of the same value as follows:

(15)

f_{l} (z_{l + 1}, ω) = f_{l + 1} (z_{l + 1}, ω) .

Thus, the relationship between d_l for the l^th layer and $d_{l + 1}$ for the (l + 1)^thlayer can be written by

(16)

d_{l} = M_{l}^{- 1} (z_{l + 1}, ω) M_{l + 1} (z_{l + 1}, ω) d_{l + 1} .

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

(17)

\begin{aligned} f_{0} (z_{0}, ω) = M_{0} (z_{0}, ω) d_{0} = M_{0} (z_{0}, ω) M_{0}^{- 1} (z_{1}, ω) \\ M_{1} (z_{1}, ω) \dots M_{N - 1}^{- 1} (z_{N}, ω) M_{N} (z_{N}, ω) d_{N} = M_{Tot} d_{N} . \end{aligned}

Since only the downgoing wave components exist in the last layer (infinite half-layer), $d_{1, N}^{+}$ and $d_{2, N}^{+}$ should be zero. As the boundary condition at the top boundary of the air layer applies $E_{source} = (E_{x, source}, E_{y, source}, E_{z, source})$ that is usually $E_{1, source}$ = (1,0,0) or $E_{2, source}$ = (0,1,0), and d_N is composed of only $d_{1, N}^{-}$ and $d_{2, N}^{-}$ , the 4 × 4 matrix in equation (17) can be reduced to the 2 × 2 matrix as follows:

(18)

(\binom{E_{x, source}}{E_{y, source}}) = (\begin{array}{ll} M_{T o t} (1, 2) & M_{T o t} (1, 4) \\ M_{T o t} (2, 2) & M_{T o t} (2, 4) \end{array}) (\binom{d_{1, N}^{-}}{d_{2, N}^{-}}) .

As $E_{x, s o u r c e}$ , $E_{y, s o u r c e}$ and the elements of $M_{Tot}$ are known, d_N can be computed using equation (18).

By using equation (16) and d_N, 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)

\begin{aligned} (\begin{matrix} E_{x} \\ E_{y} \\ H_{x} \\ H_{y} \end{matrix}) = \\ (\begin{array}{cccc} X (+ 1) & X (- 1) & 0 & 0 \\ 0 & 0 & X (+ 2) & X (- 2) \\ 0 & 0 & - γ k_{y} X (+ 2) & γ k_{y} X (- 2) \\ γ k_{x} X (+ 1) & - γ k_{x} X (- 1) & 0 & 0 \end{array}) (\begin{matrix} d_{1}^{+} \\ d_{1}^{-} \\ d_{2}^{+} \\ d_{2}^{+} \end{matrix}), \end{aligned}

where $X (\pm 1) = e x p (\pm k_{1} z) = e x p (\pm k_{x} z); X (\pm 2) = e x p (\pm k_{2} z) = e x p (\pm k_{y} z); k_{x}^{2} = i ω μ_{0} A_{x x}; and k_{y}^{2} = i ω μ_{0} A_{y y}$ .

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) = e x p (k_{j} z)$ 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 (\pm j) = e x p (\pm k_{j} z)$ to $X' (- j) = e x p (- k_{j} (z - z_{i}))$ for downgoing wave components and $X' (+ j) = e x p (k_{j} (z - z_{i + 1}))$ for upgoing wave components at the i^th 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)

\begin{aligned} (\begin{matrix} E_{x} \\ E_{y} \\ H_{x} \\ H_{y} \end{matrix}) = \\ (\begin{array}{cccc} X^{'} (+ 1) & X^{'} (- 1) & X^{'} (+ 2) & X^{'} (- 2) \\ Q_{1} X^{'} (+ 1) & Q_{1} X^{'} (- 1) & Q_{2} X^{'} (+ 2) & Q_{2} X^{'} (- 2) \\ - γ Q_{1} k_{1} X^{'} (+ 1) & γ Q_{1} k_{1} X^{'} (- 1) & - γ Q_{2} k_{2} X^{'} (+ 2) & γ Q_{2} k_{2} X^{'} (- 2) \\ γ k_{1} X^{'} (+ 1) & - γ k_{1} X^{'} (- 1) & γ k_{2} X^{'} (+ 2) & - γ k_{2} X^{'} (- 2) \end{array}) (\begin{matrix} d_{1}^{' +} \\ d_{1}^{' -} \\ d_{2}^{' +} \\ d_{2}^{' -} \end{matrix}), \end{aligned}

(21)

\begin{aligned} (\begin{matrix} E_{x} \\ E_{y} \\ H_{x} \\ H_{y} \end{matrix}) = \\ (\begin{array}{cccc} X^{'} (+ 1) & X^{'} (- 1) & 0 & 0 \\ 0 & 0 & X^{'} (+ 2) & X^{'} (- 2) \\ 0 & 0 & - γ k_{y} X^{'} (+ 2) & γ k_{y} X^{'} (- 2) \\ γ k_{x} X^{'} (+ 1) & - γ k_{x} X^{'} (- 1) & 0 & 0 \end{array}) (\begin{matrix} d_{1}^{' +} \\ d_{1}^{' -} \\ d_{2}^{' +} \\ d_{2}^{' -} \end{matrix}) . \end{aligned}

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 d´ 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' (\pm j)$ converge to zero [ $X' (\pm j) = e x p (- \infty) \approx 0$ as $z \to \infty$ ]. 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 i^th layer with the upper ( $z_{i}$ ) and lower ( $z_{i + 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 $M_{i}'^{- 1} (z_{i + 1}, ω)$ where i varies from 0 to $N - 1$ is necessary for solving the components of the EM fields. The components of $M_{i}'^{- 1} (z_{i + 1}, ω)$ are defined as $X' (j) = e x p (k_{j} (z_{i + 1} - z_{i + 1})) = e x p (0) = 1$ and $X' (- j) = e x p (- k_{j} (z_{i + 1} - z_{i}))$ . When the thickness of the i^th layer $(z_{i + 1} - z_{i})$ 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 $M_{i}' (z_{i + 1}, ω)$ also become zero. Thus, the inverse matrix of $M_{i}' (z_{i + 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 $M_{i}' (z_{i + 1}, ω)$ for each layer of a given model before computing $M' (z, ω)$ and d´. 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 $M_{i}' (z_{i + 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 $M_{i}' (z_{i + 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 $M_{i}' (z_{i + 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 10⁹ Ω·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 $(E_{x})$ 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 $E_{x}$ computed using the conventional method and SS-I when the plane wave sources, $E_{1, source}$ = (1,0,0) (Fig. 5(a)) and $E_{2, 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 $E_{x}$ for both $E_{1, source}$ and $E_{2, 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 $(E_{x})$ for (a) $E_{1, source}$ = (1,0,0) and (b) $E_{2, 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 (\pm j)$ and $X' (\pm j)$ which are the general solutions required to compute $E_{x}$ 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 (\pm j)$ up to a depth of 10 km. It is observed that the terms $X (\pm j)$ and $X' (\pm 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 (\pm 2)$ and $X' (\pm 2)$ are zero (see Figs. 6(c) and 7(c)). The terms $X (\pm j)$ and $X' (\pm j)$ computed with both the conventional method and SS-I do not suffer from numerical instabilities at shallow depths. As $z$ increases $(z \to \infty)$ , the terms $X (- 1)$ , $X (- 2)$ and $X' (\pm j)$ converge to zero as $e x p (- \infty) = 0$ . In contrast, because $X (+ 1)$ , $X (+ 2) = e x p (- \infty)$ , 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 (\pm 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' (\pm 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 d´ at each layer whose components are the coefficients required to compute $E_{x}$ . In Tables 1 and 2, there are no numerical instabilities of d and d´ at each layer. Because upgoing wave components at the last layer do not exist, the values for components $d_{1}^{+}$ , $d_{1}'^{+}$ , $d_{2}^{+}$ and $d_{2}'^{+}$ are all zero. Even if $X (+ 1)$ has large values at depths shallower than 447.6 km, $E_{x}$ is not affected because $d_{1}^{+}$ is zero, and thus $d_{1}^{+} X (+ 1) = 0$ . However, below a depth of 447.6 km, $X (+ 1)$ becomes infinity, resulting in numerical instability in the computation of $E_{x}$ . Table 2 shows that SS-I does not suffer from numerical instability in computing d´.

Table 1.

Vector d calculated with the conventional method

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

Table 2.

Vector d´ calculated with SS-I

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

Additionally, Tables 1 and 2 show that the components of d and d´ are not affected by the depth of the receiver in a given 1-D model. In contrast, the components of $X (\pm j)$ and $X' (\pm 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 $E_{x}$ using the conventional method at large depths, whereas the components of $X' (\pm 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 $(E_{x})$ computed without and with SS-II for both $E_{1, source}$ = (1,0,0) (Fig. 8(a)) and $E_{2, 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 $E_{1, source}$ , $E_{x}$ computed with SS-II converges an approximate value of 2.5 × 10^-3 V/m whereas $E_{x}$ 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, $E_{x}$ computed with SS-II for $E_{2, source}$ converges to zero, whereas $E_{x}$ 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 $M_{1}' (z_{2}, ω)$ ) is smaller than the given tolerance. If the thickness of the first layer exceeds 166.4 km for $E_{1, source}$ , which means that the determinant of $M_{1}' (z_{2}, ω)$ 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 $(E_{x})$ computed for (a) $E_{1, source}$ = (1,0,0) and (b) $E_{2, 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 $E_{x}$ and how SS-II stabilizes the computational process. Because the first layer shown in Fig. 3 is isotropic, $E_{x}$ at the surface can be defined as $E_{x} = X_{1}'^{+} d_{1}'^{+} + X_{1}'^{-} d_{1}'^{-} + 0 * d_{2}'^{+} + 0 * d_{2}'^{-}$ according to equation (21). $X' (+ 1)$ can be defined as $e x p (k_{1} (- z_{2}))$ at the surface $(z = 0)$ and as $z_{2}$ increases, $X' (+ 1)$ converges to zero, and the term $X' (- 1)$ at $z = 0$ has a constant value of 1 as $e x p (- k_{1} (0 - 0)) = e x p (0) = 1$ .

The components of d´ for each layer can be computed using the continuity condition, which require calculating the inverse of $M_{1}' (z_{2}, ω)$ . $M_{1}' (z_{2}, ω)$ 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 $M_{1}' (z_{2}, ω)$ . Fig. 10 shows the determinant of $M_{1}' (z_{2}, ω)$ . Because the first layer is isotropic (i.e., $k_{1} = k_{x} = k_{y}$ ), the components $X' (+ 1) = e x p (k_{1} (z_{2} - z_{2})) = e x p (0) = 1$ as shown in Fig. 9(a) and 9(c). However, the components $X' (- 1) = e x p (- k_{1} (z_{2} - z_{1}))$ converges to zero, when $z_{2}$ continues to increase, as shown in Fig. 9(b) and 9(d). Similarly, $X' (- 2)$ in equation (21) converges to zero as $z_{2}$ increases. Consequently, as the thickness of the first layer increases, the elements of the second and fourth columns of $M_{1}' (z_{2}, ω)$ 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 $M_{1}' (z_{2}, ω)$ becomes numerically unstable, which also result in unstable computation of d´ 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) $X_{1}'^{+}$ and (b) $X_{1}'^{+}$ and imaginary components of (c) $X_{1}'^{+}$ and (d) $X_{1}'^{+}$ for $M_{1}' (z_{2}, ω)$ 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 $M_{1}' (z_{2}, ω)$ with respect to the thickness of the first layer ranging from 2 to 1002 km.

Fig. 11 shows the real components of d´ 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 $d_{1}'^{+}$ , $d_{1}'^{-}$ , $d_{2}'^{+}$ and $d_{2}'^{-}$ 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 $E_{x} = X_{1}'^{+} d_{1}'^{+} + X_{1}'^{-} d_{1}'^{-} + 0 * d_{2}'^{+} + 0 * d_{2}'^{-}$ in the first layer, it results in unstable computation of $E_{x}$ 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 d´ does not require the inverse matrix of $M_{1}' (z_{2}, ω)$ , preventing numerical instabilities as the components $d_{1}'^{+}$ , $d_{1}'^{-}$ , $d_{2}'^{+}$ and $d_{2}'^{-}$ 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) $d_{1}'^{+}$ , (b) $d_{1}'^{-}$ , (c) $d_{2}'^{+}$ and (d) $d_{2}'^{-}$ of vector d´ 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) $ρ_{x y}$ , (b) $ϕ_{x y}$ 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

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

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

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

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

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

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

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

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

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

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

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

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

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

Journal of the Korean Society of Mineral and Energy Resources Engineers ISSN:2288-0291(Print) 2288-2790(Online) 한국자원공학회지

Preview

Stable Calculation of Responses for a Plane Electromagnetic Wave in an Anisotropic Layered Media for the Modeling of Magnetotelluric Data

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

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

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

A four-layered 1-D anisotropic model.

Fig. 4.

Fig. 5.

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.

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.

Table 1.

Vector d calculated with the conventional method

Table 2.

Vector d´ calculated with SS-I

Fig. 8.

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.

Fig. 10.

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

Fig. 11.

The real components of (a) d1'+, (b) d1'-, (c) d2'+ and (d) d2'- of vector d´ 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.

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.

Acknowledgements

References

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

(a)-(d): The real components and (e)-(h): imaginary components of $X (\pm 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.

(a)-(d): The real components and (e)-(h): imaginary components of $X' (\pm 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.

The real components of (a) $X_{1}'^{+}$ and (b) $X_{1}'^{+}$ and imaginary components of (c) $X_{1}'^{+}$ and (d) $X_{1}'^{+}$ for $M_{1}' (z_{2}, ω)$ computed without SS-II denoted with the orange markers and with SS-II denoted with the solid yellow lines.

The determinant of $M_{1}' (z_{2}, ω)$ with respect to the thickness of the first layer ranging from 2 to 1002 km.

The real components of (a) $d_{1}'^{+}$ , (b) $d_{1}'^{-}$ , (c) $d_{2}'^{+}$ and (d) $d_{2}'^{-}$ of vector d´ 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.