Development of a Three-Dimensional Inversion Program for Magnetotellurics Using Finite Volume and Adjoint State Method

Atsushi Suzuki

doi:10.20944/preprints202408.2236.v2

Submitted:

21 October 2024

Posted:

21 October 2024

Read the latest preprint version here

Abstract

The finite volume method (FVM) is a kind of numerical analysis method that is simple and can take into account the topography as the finite element method (FEM). The adjoint state method (ASM) can easily and efficiently calculate the gradient vector of the objective function for the inversion. We developed a program to apply these methods to 3-D Magnetotelluric (MT) inversion and verified it using test models. In the forward modelling, it was confirmed that the program could perform analysis properly and efficiently considering topography and subsurface resistivity heterogeneity. Furthermore, we conducted the inversion using a synthetic open dataset and could reconstruct the true model of the dataset. The comparison with the previous work revealed that the computation time was comparable in spite of the fact that we calculated without MPI implementation, which the previous work used.

Keywords:

Magnetotellurics

;

inversion

;

Finite volume method

;

Adjoint state method

Subject:

Environmental and Earth Sciences - Geophysics and Geology

Introduction

The magnetotelluric (MT) method is a geophysical exploration method to estimate the resistivity structure in the subsurface and is currently used in, for example, geothermal and oil and gas field exploration(Takakura 2014; Key et al., 2006). In the MT method, 3D analysis is now often conducted to estimate the resistivity in the subsurface from the acquired data, and the inversion programs were developed for this purpose (Sasaki, 2004; Siripunvaraporn et al., 2005; Kelbert et al., 2014; Usui, 2015; Usui et al., 2024). In addition, a comparison of them was also performed (Uchida and Yamaya, 2023).

In the 3D MT analysis, forward modelling is performed to analyse the propagation of electromagnetic waves into the subsurface. In general, this problem cannot be solved analytically. Hence, some numerical analysis is necessary. For this purpose, the finite difference method (FDM), FEM, and FVM are often used (e.g., Sasaki, 2004; Siripunvaraporn et al., 2005; Nam et al., 2007; Kelbert et al., 2014; Jahandari and Farquharson, 2015; Usui, 2015; Zhu et al., 2022; Usui et al., 2024). In the modelling, the assessment of the effect of the topography is necessary because it appears as a distortion of the data (Ogawa, 2002). Although FDM is a simple discretization method, it is complicated to apply the unstructured grid. Therefore, we need to express the topography with the stack of the rectangular grid in FDM. It causes the necessity of the fine grid, or the error of the calculation becomes large (Uchida, 2022; Uchida and Yamaya, 2023). On the other hand, FEM and FVM can treat the unstructured grid along the topography (Jahandari & Farquharson, 2015; Nam et al., 2007; Usui, 2015; Zhu et al., 2022; Usui et al., 2024). In these methods, we can express the topography smoothly, and the error caused by the stepped grid in FDM is expected to disappear. From these perspectives, we developed a program that adopts the hexahedral cell-centred FVM. The advantages of this method are:

(1): The unknowns are placed at the cell centre, and the spatial partial differentiation is discretized on the cell surface with a simple method similar to FDM.
(2): Unstructured meshes can be applied as FEM, allowing calculations with the topography.

Further, for the inversion, a process to update the 3D resistivity model so that the calculated data are matched to the observed one is required. MT inversion often uses the Jacobi matrix of the objective function for this purpose (e.g., Sasaki, 2004; Siripunvaraporn et al., 2005; Kelbert et al., 2014; Usui, 2015; Usui et al., 2024). However, the computational cost of acquiring the Jacobi matrix is enormous (Plessix 2006). Especially, when we use iterative solvers in the forward modelling, we need to solve equations where many different source-term vectors are applied (e.g. Rodi 1976). Uchida (2022) analysed actual data using a program which uses an iterative solver in various grid-size models. They found that some models require a large portion of the overall computation time to calculate the Jacobi matrix. We can avoid the tough solutions when we use direct solvers in which we can reuse, for example the result of LU decomposition. However, in general, direct solvers require more memory than iterative solvers and take much time in the case of insufficient computation resources.

A workaround is that we use the gradient, not the Jacobi matrix, obtained by ASM to match between calculated data and observed data. ASM is an efficient way of finding the gradient vector (Plessix 2006). The advantage of this method is that the computationally expensive part is only one solution of the adjoint equation, so we do not need to solve the equations many times. In addition, the scale of the adjoint equation is the same as the forward modelling, regardless of the number of model parameters or data. Therefore, we can use iterative solvers without necessity of many solutions with different source-term vectors. Once the gradient is obtained, we can update the model in for example, BFGS, which is one of Quasi-Newton methods, or the gradient descent method like Adam (Kingma and Ba, 2014), which is recently often used in the field of machine learning.

We developed a program that applied the FVM and ASM together to achieve efficient analysis considering the topography. This paper reviews and discusses the methodology and results of some example models.

Method

Forward Modelling

In MT forward modelling, on the assumption that displacement current is ignored and Ohm’s law is established, two equations in Maxwell’s equations in the frequency domain are considered.

\nabla \times E = - i ω μ_{0} H

(1)

\nabla \times H = \frac{1}{ρ} E

(2)

where E is the electric field, H is the magnetic field, ω is the angular frequency, μ₀ is the magnetic permeability of the vacuum, and ρ is specific resistance.

We usually combine these equations and formulate a second-order equation. Here, there are two choices that we select the electric or magnetic field as the unknown of the equation (Siripunvaraporn et al., 2002). In this program, we choose the H-forming equation (3), namely,

\nabla \times ρ \nabla \times H = - i ω μ_{0} H

(3)

Based on (3), we need to discretise (3) for numerical analysis.

First, In FVM, (3) is integrated in volume in both sides,

\int_{V} \nabla \times ρ \nabla \times H d V = \int_{V} - i ω μ_{0} H d V,

(4)

where V is each cell volume.

Next, the left side of (4) is converted to surface integration by integral formula.

\int_{S} n \times ρ \nabla \times H d S = \int_{V} - i ω μ_{0} H d V,

(5)

where S is the surface of each cell and n is the unit vector orthogonal to the surface. We can discretise (5), namely,

\sum_{i = 1}^{6} n_{i} \times ρ_{i} {(\nabla \times H)}_{i} ∆ S_{i} = - i ω μ_{0} H_{c} ∆ V,

(6)

where, subscript i and c mean each surface of the hexahedral cell and the value at the cell centre respectively,

∆ S_{i}

is the area of each surface and

∆ V

is the volume of the cell.

In (6), because we locate the unknowns H at the cell centre because we adopt cell-centred FVM,

H_{c}

is the unknown which we should solve. Hence, we can obtain the discretised equations of each cell when discretising the term of

ρ_{i} {(\nabla \times H)}_{i}

on the surface properly.

Note that (6) is valid with arbitrary geometry; therefore, it is possible to use unstructured meshes and local mesh refinement as FEM in Usui et al. (2024).

For the air layer, (3) should be modified because of ill-condition caused by the huge resistivity of the air. We can achieve it by converting (3), assuming

ρ i s c o n s t a n t i n t h e a i r l a y e r

, and applying the condition

\nabla \cdot H = 0

as follows:

\nabla^{2} H = i ω μ_{0} σ H,

(7)

where

σ = 1 / ρ .

We can discretise it as (8) in the same manner of the subsurface,

\sum_{i = 1}^{6} n_{i} \cdot {(\nabla H)}_{i} ∆ S_{i} = i ω μ_{0} σ H_{c} ∆ V,

(8)

Equation (7) is general Poisson equation, hence using (8) instead of (6) in the air layer should contribute to the computational stability. We solve (6) and (8) in the subsurface and the air layers, respectively.

In (6),

ρ

is defined on the surface of the cell, so we need an interpolation from the resistivity defined at the cell centre to the one on the cell surface is necessary. This program uses the average strategy like Brewitt-Taylor and Weaver (1976) and Mackie et al. (1993). Namely, we initially calculate the average resistivities on the edges of cells (Figure 1) using the cell-centred resistivities connected to the edges. Then, the resistivity on the centre of the cell surface is obtained with the mean of the edge resistivities weighted by the reciprocal of distance.

For computational stability, a transition layer is provided between the air and subsurface (Figure 2). This layer belongs to the equation (6) and has the resistivity of the subsurface in the cell centre. In this way, the resistivity of the surface in this layer change from huge on adjacent surfaces to the air to small on ones to the subsurface.

To discretize

\nabla \times H

in the subsurface layer and

\nabla H

in the air layer on the surface of the cell,

\partial H ⁄ \partial x_{j}

on the cell surface, where

j

is 1, 2, or 3 and

x_{j}

is a global coordinate system, must be calculated. First, we calculate the spatial partial differentiation in the local coordinate system

(e_{1}, e_{2}, e_{3})

, where

e

unit vectors, using the interpolation value

on the edges (Figure 1) in the same way of the resistivity and the scheme of FDM. Then, we convert the spatial partial differentiation from the local coordinate system

(e_{1}, e_{2}, e_{3})

to the global coordinate system

(e_{x 1}, e_{x 2}, e_{x 3})

, as follows:

(\begin{matrix} \frac{\partial H}{\partial x_{1}} \\ \frac{\partial H}{\partial x_{2}} \\ \frac{\partial H}{\partial x_{3}} \end{matrix}) = {(\begin{matrix} e_{1} \cdot e_{x_{1}} & e_{1} \cdot e_{x_{2}} & e_{1} \cdot e_{x_{3}} \\ e_{2} \cdot e_{x_{1}} & e_{2} \cdot e_{x_{2}} & e_{2} \cdot e_{x_{3}} \\ e_{3} \cdot e_{x_{1}} & e_{3} \cdot e_{x_{2}} & e_{3} \cdot e_{x_{3}} \end{matrix})}^{- 1} (\begin{matrix} \frac{\partial H}{\partial e_{1}} \\ \frac{\partial H}{\partial e_{2}} \\ \frac{\partial H}{\partial e_{3}} \end{matrix})

(9)

Here, we can obtain the terms of

\nabla \times H

and

\nabla H

from the left side of (9).

Boundary conditions on XZ and YZ planes are

\partial H / \partial n

equal to zero. On XY of the subsurface side,

H = 0

is adopted. On XY of the air layer side, the boundary conditions are source terms of the modelling. We set

H_{x}

and

H_{y}

equal to one in twice forward modelling to obtain response functions, respectively (e.g., Siripunvaraporn et al., 2002).

Through these operations, the linear equation

A x = b

(

A

: coefficient matrix,

x

: unknown magnetic field vector,

b

: source term vector) can be assembled. To solve the linear equation, we adopt BiCGSafe (Fujiwara et al, 2006) with block incomplete LU decomposition as a preconditioner as Toh et al. (2000). We set the convergence condition using the maximum change of the solution from the previous iteration

∆ x_{i}^{m a x} = \max (x_{i} - x_{i - 1})

, where i is iteration number and the residual

r_{i} = A x_{i} - b

.

\min (\frac{|r_{i}|}{|r_{0}|}, \frac{|∆ x_{i}^{m a x}|}{|x_{i}^{m a x}|}) < ϵ,

(10)

where

ϵ

is tolerance.

The divergence correction for subsurface layers is applied in the way in Dong and Egbert (2019). They add a term

k \nabla (\nabla \cdot H)

to equation (3), where

k

is the scaling factor. In this program,

k

is set to

s ρ_{m}

, where

s

and

ρ_{m}

are the safety factor equal to 1.0 and the average of the resistivity surrounding cells respectively in the transition zone, and 0.01 and minimum of the resistivity surrounding cells respectively in the subsurface. Although Dong and Egbert (2019) does not adopt the safety factor, we do because the large value of

k

in the subsurface resulted in poor accuracy solutions in our test case.

Once

H

is obtained,

E

is necessary in the cell centre to calculate the impedance tensor. Here, (2) is adopted to obtain

E

. Then, we calculate two components of

ρ \nabla \times H

parallel to each surface are calculated (Figure 1).

E^{c e n t e r} \cdot p_{i j} = E_{i j}^{p a r a l l e l},

(11)

where

E^{c e n t e r}

is the required electric field at the cell centre,

E^{p a r a l l e l}

is the electric field on the cell surface calculated from equation (2),

p

is the unit vector parallel to the cell surface,

i

is unit vector directions

(1 o r 2)

, and

j i s c e l l s u r f a c e s n u m b e r (f r o m 1 t o 4)

. The component of the electric field orthogonal to cell surfaces is excluded because it is not included in (6) and found to cause the numerical error. From equation (11), the eight equations associated with

E^{c e n t e r}

in each cell can be obtained. Therefore,

E^{c e n t e r}

can be calculated using the least squares method Because

E^{c e n t e r}

has three components. Using

H

and

E^{c e n t e r},

The impedance tensor and magnetic transfer function can be obtained in the same way as in Usui (2015).

Since this program uses the cell-centred FVM, it is necessary to use the surrounding cell-centred magnetic field to obtain the electric field at the cell surface. This program calculates the impedance tensor and magnetic transfer function on the earth’s surface using the electromagnetic fields in the transition zone and near-surface cells, as shown in Figure 2, which is a 2D case.

For the calculation of an impedance tensor with distortion, the method adopted by Ogawa (2002) is used.

Z_{c a l c} = D Z,

(12)

where

Z_{c a l c}

is the calculated impedance tensor affected distortion,

D

is the distortion tensor, and

Z

is the originally calculated impedance tensor.

Note that C++ language is used to implement this program, and the library Eigen (Jacobi et al., 2024) is used for the internal matrix operations. Parallelization is implemented by OpenMP. The linear equations

A x = b

at each frequency are solved on each thread.

Inversion

We implemented ASM for MT inversion. The details are described in Plessix (2006). In MT inversion, an objective function

F = F (m)

, where

m \in R^{N_{m o d e l}}

is the vector of model parameter, and

N_{m o d e l}

is the number of model parameters, should be minimized. In ASM, the objective function

F

is considered as the independent variables are

m

and

x

∈ ℂ^3Ncell, where

x

is the solution of linear equation

A x = b

in forward modelling and N_cell is the total number of computation cells, that is,

F = F (m, x)

(13)

Next, the adjoint equation is solved about

λ

.

A^{*} λ = {[\frac{\partial F {(m, x)}^{*}}{\partial x}]}^{T}

(14)

Note that (14) is the same scale equation as

A x = b

. Hence, the computational cost is roughly the same as solving

A x = b

.

The gradient vector of the objective function can be obtained using

λ

as follows:

\frac{\partial F (m)}{\partial m_{i}} = \frac{\partial F (m, x)}{\partial m_{i}} - R e (λ^{*} \frac{\partial A}{\partial m_{i}} x),

(15)

where

i

is from

1

to

N_{m o d e l}

. Note that

\partial F (m, x) / \partial m_{i}

, and

\partial F {(m, x)}^{*} / \partial x

are calculated by automatic differentiation using library kv (Kashiwagi 2024).

The objective function

F

is set as below in this program.

F (m) = {|W d_{o b s} - W d_{c a l c}|}^{2} + α^{2} {|R \log (m)|}^{2} + β^{2} \sum_{n = 1}^{N_{s i t e}} \sum_{i = 1}^{2} \sum_{j = 1}^{2} (D_{i j}^{n} - I_{i j}),

(16)

where

d_{o b s}

is the observed data vector,

d_{c a l c}

is calculated data vector, W is weight matrix, R is smoothness matrix for constraint term which is necessary in ill-posed MT inversion,

N_{s i t e}

is the total number of observation sites where impedance tensor is obtained,

D^{n} \in C^{2 \times 2}

is distortion tensor in each observation site n as described in Ogawa (2002), I is identity matrix,

α^{2}

and

β^{2}

are trade-off parameters.

The model parameters are constraint on the limitations of upper and lower boundary in the way of (Kim and Kim 2011).

To minimize the objective function F with the gradient vector, Adam (Kingma and Ba, 2014) is implemented, which can be switched easily because the library OptimLib (O’Hara 2024) is used internally.

Numerical Examples and Discussion

Forward Modelling Considering Topography

To evaluate if the developed program can consider the topography, the trapezoidal hill model was adopted. This model is often used to validate the effect of topography (e.g., Nam et al., 2007; Usui, 2015; Zhu et al., 2022, Figure 3).

The minimum size of cells was 25 m ×25 m ×10 m and the total number of cells was 519,456. The resistivity in the air and the subsurface was

10^{6}

and 100

Ω m

respectively. The calculated frequency was 2 Hz. These conditions are the same as previous works (Nam et al., 2007; Usui, 2015; Zhu et al., 2022). The computer we used was Precision 3630 Tower (CPU: Intel Core i7-9700 CPU 3.00GHz, RAM:16GB).

The results are shown in Figure 4. The results are compared to Zhu et al. (2022). Except for some points, these results are mostly consistent with each other. The gaps between them would be caused by the difference in the discretization method, as which Zhu et al. (2022) adopted FEM, and the size of cells. However, since the differences are quite small, it would not influence when the inversion is performed. Furthermore, the tendency of the change in the apparent resistivity and phase along the hill are reproduced well. Hence, we can conclude that this program can consider the topographic effect.

Forward Modelling Considering Resistivity Heterogeneity

As a test case for heterogeneity of subsurface resistivity, the model in Mackie et al. (1993) was adopted (Figure 5). The total number of cells was 370,712, with minimum size of cells 500 m ×500 m ×200 m. The resistivity in the air was

10^{6}

Ω m

. The compared frequency was 10 and 1000 seconds. The computer is the same one used in the Trapezoidal hill model.

The results of Apparent resistivity and phase in cross sections along the Y direction at X=0 m for periods of 10 and 1000 s are compared with Mackie et al. (1993) in Figure 6 and Figure 7. While there are some small differences, the results are in good agreement with those of Mackie et al. (1993). Therefore, we also conclude that the heterogeneity of the resistivity can be treated properly.

Inversion for Dublin Secret Model 2

We performed inversion using Dublin Secret Model 2 (DSM2) in Miensopust et al. (2013). This model was used to compare different programs in their work. The synthetic observed data of impedance tensor are openly available in Miensopust et al. (2013). Note that the synthetic data is not produced using our program. Hence, we can confirm not only the function of inversion, but also the forward modelling because the result of the inversion should not be consistent with the true model if the result of the forward modelling is quite different from the one which was used to produce the synthetic data in Miensopust et al. (2013).

As the settings of the inversion, the total number of cells was 268,960. The number of degrees of freedom was 806,880. The number of cells that were inverted as model parameters was 209,884. The minimum size of cells was 1000 m ×1000 m ×100 m. The initial resistivity was 50

Ω m

. The total number of calculated periods was 11, from 0.025 to 2500 s. We

initially set

α^{2}

to 100 and lowered it after the convergence criteria were achieved. The criteria were that the maximum change of model parameters or the change of objective function from the previous iteration was below 0.3 %. For simplicity,

β^{2}

was set to the same value as

α^{2}

. The computer for this calculation was Dell XPS 8950 (CPU: Intel Core i7-12700K 3.60 GHz, RAM:64GB).

The results of the inverted resistivity structure compared to the true model are shown in Figure 8. We adopted the model when

α^{2}

was 0.1, judging from the L-curved figure between RMS and model roughness (Figure 9) as in Usui et al. (2024). The resistive and conductive anomalies shallower than 10 km are inverted. While the conductive anomaly is connected to one deeper than 30 km, this would be caused by low sensitivity under the conductive anomaly, as mentioned in Usui (2015) and Kelbert et al. (2014). The other structures are reconstructed well in view of the size and value of the anomalies.

Considering that the synthetic data is distorted according to Groom and Bailey (1989) and has 5% random noise of maximum impedance tensor value (Miensopust et al. 2013), the results of our program seem in agreement with the true model.

The total calculation time was about 23 hours. Usui (2015) calculated the same model in about 15 hours, with a number of degrees of freedom of 687,599, roughly the same as our calculation. Our calculation was performed on an ordinary personal computer without MPI parallelization, while Usui (2015) did with it because they used a direct solver for the forward modelling. Although the conditions of computation such as the computers and the total number of calculated frequencies were different, we could obtain the results in a comparable time as Usui (2015). This would suggest the efficient performance of our program.

From the view of the algorithm, our program accommodates to be implemented in MPI parallelization. Especially, the MPI parallelization of matrix-vector multiplication would be efficient for the calculation speed because the multiplication in the iterative solver of

A x = b

costs large computation resources and appears twice every iteration in BiCGSafe (Fujiwara et al, 2006). Our program uses the library Eigen (Jacobi et al., 2024) and its multiplication has been implemented in multi-thread parallelization. Hence, when we solve

A x = b

in each frequency in each node using MPI parallelization, the calculation speed would improve. This is a future work.

Conclusions

We developed a program for 3D MT inversion. For the forward modelling and inversion, FVM and ASM were adopted. FVM can consider the unstructured grid where the smooth topography expression can be treated. ASM is an efficient way to obtain the gradient of the objective function. Once we obtain the gradient, we can minimize the objective function in the Quasi-Newton method like BFGS or the gradient descent method like Adam.

As test models for the forward modelling, the trapezoidal hill model in Nam et al. (2007) and a model that contains resistivity heterogeneity in Mackie et al. (1993) were calculated. The differences between our results and the previous works were slight. Therefore, the program can correctly calculate electromagnetic wave propagation.

For the test model of inversion, we adopted DSM2. The results were in good agreement with the true model and the ones in previous works. The calculation time was also comparable to the previous work with an ordinary personal computer. As a future work, in addition to OpenMP parallelization, MPI parallelization would be necessary for our program to perform more fast and efficient calculations.

Figure 1. The schematic cartoon that indicates how to calculate

\partial H / \partial e_{i}

on the surface of each cell and the locations of the electric field.

e_{1}

and

e_{2}

are parallel to the surface

S

. Note that this calculation can be applied to not only the structured grid but also the unstructured and non-conforming grid.

Figure 2. Schematic cartoon of a 2D case that indicate where electromagnetic fields are calculated. H and E denotes magnetic and electric fields, respectively.

Figure 3. (a) XY and (b) YZ plane with used mesh of trapezoidal model in Nam et al. (2007).

Figure 4. (a) Apparent resistivity of YX mode and (b) XY mode. (c) The phase of YX mode and (d) ZY mode of 2 Hz. Blue lines and diamonds denote our results. Black dashed lines represent the result in Zhu et al. (2022).

Figure 5. (a) XY and (b) YZ plane of the resistivity structure model in our forward calculation after Mackie et al. (1993). Transparent black lines denote cell edges.

Figure 6. (a) Apparent resistivity of YX mode and (b) XY mode. (c) The phase of YX mode and (d) ZY mode of 10 s. Blue diamonds denote our results. Black dashed lines represent the result in Mackie et al. (1993).

Figure 7. (a) Apparent resistivity of YX mode and (b) XY mode. (c) The phase of YX mode and (d) ZY mode of 1000 s. Blue diamonds denote our results. Black dashed lines represent the result in Mackie et al. (1993).

Figure 8. Each cross-section of the true DSM2 model in Miensopust et al. (2013) and our results. (a) True Model and (b) Our inversion results.

Figure 9. The L-curve picture of our inversion results. The vertical axis and horizontal axis are RMS and model roughness in each

α^{2}

, respectively. We chose

α^{2} = 0.1

(red one) as our final result.

Acknowledgement

The authors report there are no competing interests to declare.

References

Brewitt-Taylor, C R, and J T Weaver. 1976. “On the Finite Difference Solution of Two-Dimensional Induction Problems.” Geophysical Journal of the Royal Astronomical Society 47:375–96. [CrossRef]
Dong, Hao, and Gary D. Egbert. 2019. “Divergence-Free Solutions to Electromagnetic Forward and Adjoint Problems: A Regularization Approach.” Geophysical Journal International 216 (2): 906–18. [CrossRef]
Fujiwara, Maki, Masahiro Yoshida, and Seiji Fujino. 2006. “BiCGSafe Method with Crout Version of ILU Decomposition Giving Triple Safe Keys for Convergence.” IPSJ Transactions on Advanced Computing Systems 47 (7): 52–60.
Groom, Ross W, and Richard C Bailey. 1989. “Decomposition of Magnetotelluric Impedance Tensors in the Presence of Local Three-Dimensional Galvanic Distortion.” JOURNAL OF GEOPHYSICAL RESEARCH 94 (B2): 1913–25. [CrossRef]
Jacobi, Benoît, Gaël Guennebaud, and other contributor. 2024. “Eigen.”.
Jahandari, H., and C. G. Farquharson. 2015. “Finite-Volume Modelling of Geophysical Electromagnetic Data on Unstructured Grids Using Potentials.” Geophysical Journal International 202 (3): 1859–76. [CrossRef]
Kashiwagi, Masahide. 2024. “Kv - a C++ Library for Verified Numerical Computation.”.
Kelbert, Anna, Naser Meqbel, Gary D. Egbert, and Kush Tandon. 2014. “ModEM: A Modular System for Inversion of Electromagnetic Geophysical Data.” Computers and Geosciences 66 (May):40–53. [CrossRef]
Key, Kerry W., Steven C. Constable, and Chester J. Weiss. 2006. “Mapping 3D Salt Using the 2D Marine Magnetotelluric Method: Case Study from Gemini Prospect, Gulf of Mexico.” Geophysics 71 (1). [CrossRef]
Kim, Hee Joon, and Younghee Kim. 2011. “A Unified Transformation Function for Lower and Upper Bounding Constraints on Model Parameters in Electrical and Electromagnetic Inversion.” Journal of Geophysics and Engineering 8 (1): 21–26. [CrossRef]
Kingma, Diederik P., and Jimmy Ba. 2014. “Adam: A Method for Stochastic Optimization,” December. http://arxiv.org/abs/1412.6980.
Mackie, Randall L, Theodore R Madden, and Wannamaker Philip E. 1993. “Three-Dimensional Magnetotelluric Modeling Using Difference Equations - Theory and Comparisons to Integral Equation Solutions.” GEOPHYSICS. Vol. 58. http://segdl.org/.
Miensopust, Marion P., Pilar Queralt, Alan G. Jones, Dmitry Avdeev, Anna Avdeeva, Ralph Uwe Börner, David Bosch, et al. 2013. “Magnetotelluric 3-D Inversion-a Review of Two Successful Workshops on Forward and Inversion Code Testing and Comparison.” Geophysical Journal International 193 (3): 1216–38. [CrossRef]
Nam, Myung Jin, Hee Joon Kim, Yoonho Song, Tae Jong Lee, Jeong-Sul Son, and Jung Hee Suh. 2007. “3D Magnetotelluric Modelling Including Surface Topography.” Geophysical Prospect 55 (2): 277–87. [CrossRef]
Ogawa, Yasuo. 2002. “ON TWO-DIMENSIONAL MODELING OF MAGNETOTELLURIC FIELD DATA.” Surveys in Geophysics 23:251–73. [CrossRef]
O’Hara, Keith. 2024. “OptimLib.”.
Plessix, Rene Edouard. 2006. “A Review of the Adjoint-State Method for Computing the Gradient of a Functional with Geophysical Applications.” Geophysical Journal International. [CrossRef]
Rodi, William L. 1976. “A Technique for Improving the Accuracy of Finite Element Solutions for Magnetotelluric Data.” Vol. 44. https://academic.oup.com/gji/article/44/2/483/681963. [CrossRef]
Sasaki, Yutaka. 2004. “Three-Dimensional Inversion of Static-Shifted Magnetotelluric Data.” Earth Planets Space. Vol. 56. [CrossRef]
Siripunvaraporn, Weerachai, Gary Egbert, and Yongwimon Lenbury. 2002. “Numerical Accuracy of Magnetotelluric Modeling: A Comparison of Finite Difference Approximations.” LETTER Earth Planets Space. Vol. 54. [CrossRef]
Siripunvaraporn, Weerachai, Gary Egbert, Yongwimon Lenbury, and Makoto Uyeshima. 2005. “Three-Dimensional Magnetotelluric Inversion: Data-Space Method.” Physics of the Earth and Planetary Interiors 150 (1-3 SPEC. ISS.): 3–14. [CrossRef]
Takakura, Shinichi. 2014. “Estimation of a Regional Geothermal System by the Electromagnetic Exploration.” BUTSURI-TANSA 67 (3): 195–203. [CrossRef]
Toh, Hiroaki, Adam SCHULTZ, and Makoto Uyeshima. 2000. “Electromagnetic Induction in Fully Heterogeneous Spheres Using the Staggered Grid Finite Difference Method.”.
Uchida, Toshihiro. 2022. “Consideration on Three-Dimensional Inversion of Magnetotelluric Data by Finite-Difference Modeling Including Topography: Data Interpretation at Northern Hakkoda Area.” Journal of the Geothermal Research Society of Japan 44 (1): 7–21.
Uchida Toshihiro, and Yusuke Yamaya. 2023. “Three-Dimensional Finite-Element and Finite Difference Inversions of Magnetotelluric Data: Application at Northern Hakkoda Geothermal Area.” Journal of the Geothermal Research Society of Japan 45 (3): 175–94. [CrossRef]
Usui, Yoshiya. 2015. “3-D Inversion of Magnetotelluric Data Using Unstructured Tetrahedral Elements: Applicability to Data Affected by Topography.” Geophysical Journal International 202 (2): 828–49. [CrossRef]
Usui, Yoshiya, Makoto Uyeshima, Hideaki Hase, Hiroshi Ichihara, Koki Aizawa, Takao Koyama, Shin’ya Sakanaka, et al. 2024. “Three-Dimensional Electrical Resistivity Structure Beneath a Strain Concentration Area in the Back-Arc Side of the Northeastern Japan Arc.” Journal of Geophysical Research: Solid Earth 129 (5). [CrossRef]
Zhu, Xiaoxiong, Jie Liu, Yian Cui, and Chunye Gong. 2022. “A Scalable Parallel Algorithm for 3-D Magnetotelluric Finite Element Modeling in Anisotropic Media.” IEEE Transactions on Geoscience and Remote Sensing 60. [CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

Development of a Three-Dimensional Inversion Program for Magnetotellurics Using Finite Volume and Adjoint State Method

Abstract

Keywords:

Subject:

Introduction

Method

Forward Modelling

Inversion

Numerical Examples and Discussion

Forward Modelling Considering Topography

Forward Modelling Considering Resistivity Heterogeneity

Inversion for Dublin Secret Model 2

Conclusions

Acknowledgement

References

MDPI Initiatives

Important Links

Subscribe