Research Article, Geoinfor Geostat An Overview Vol: 14 Issue: 1

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

Atsushi Suzuki^*

Department of Materials Science, Yokohama National University, Hodogaya-ku, Yokohama, Japan

*Corresponding Author: Atsushi Suzuki
Department of Materials Science, Yokohama National University, Hodogaya-ku, Yokohama, Japan
E-mail: suzuki@tankai.co.jp

Received date: 02 September, 2024, Manuscript No. GIGS-24-147091; Editor assigned date: 05 September, 2024, PreQC No. GIGS-24-147091 (PQ); Reviewed date: 20 September, 2024, QC No. GIGS-24-147091; Revised date: 03 February, 2026, Manuscript No. GIGS-24-147091 (R); Published date: 10 February, 2026, DOI: 10.4172/2327-4581.1000444.

Citation: Suzuki A (2026) Development of a Three-Dimensional Inversion Program for Magnetotellurics Using Finite Volume and Adjoint State Method. Geoinfor Geostat: An Overview 14:1.

Abstract

Keywords: Magnetotellurics, Inversion, Finite volume method, Adjoint state method

Keywords

Magnetotellurics; Inversion; Finite volume method; Adjoint state method

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 exploration and oil and gas field exploration [1]. In the MT method, 3D analysis is now often conducted to estimate the resistivity structure in the subsurface from the acquired data, and the inversion programs were developed for this purpose [2]. In addition, a comparison of them was also performed [3].

In the 3D MT analysis, forward modeling is performed to analyze the propagation of electromagnetic waves into the subsurface. In general, this calculation needs to 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). Although FDM can be a simple discretization method, it is not easy to apply the unstructured grid. On the other hand, FEM and FVM can treat the unstructured grid and calculate electromagnetic response [4].

Further, inversion is also performed to estimate the 3D resistivity model in the subsurface from actual data. This process needs to update the current model based on the result obtained from the forward modeling. MT inversion often uses the Jacobi matrix of the objective function to update model parameters in inversion (e.g., Sasaki). However, the computational cost of acquiring the Jacobi matrix is enormous [5]. Uchida analyzed actual data using various mesh size models and found that some require a large portion of the overall computation time to calculate the Jacobi matrix. ASM is an efficient way of finding the gradient vector rather than calculating the Jacobi matrix from a point of computational cost [6]. The advantage of this method is that the computationally expensive part is only one solution of the adjoint equation, whose scale is the same as the forward modeling, regardless of the number of model parameters or data.

The author developed a program that applied the FVM and ASM to the 3D MT inversion to achieve efficient analysis considering the unstructured grid. This paper reviews and discusses the methodology and results of some example models.

Materials and Methods

Forward modelling

In general, the MT analysis combines Maxwell's equations in the frequency domain to form an equation of only the electric or magnetic field, which is then discretized and solved to obtain the electric or magnetic field.

Some discretization methods exist including the Finite Difference Method (FDM) FEM and FVM.

This program adopts the hexahedral cell-centered finite volume method as the forward modeling. The advantages of this method are:

• The unknowns are placed at the cell center and discretized on the cell surface using a simple method similar to FDM.
• Unstructured meshes can be applied as FEM, allowing calculations with topography and local subdivision.

In MT forward modeling, the second-order Maxwell’s equations in the form of the electric E or the magnetic H fields are often solved [7]. The H-forming Equation (1) is solved in this program (Figure 1).

Figure 1: Schematic cartoon that indicates how to calculate ∂H/∂e_i on the surfaces of cells and the locations of electric field. e₁ and e₂ are parallel to the surface S. Note that this calculation can be applied to not only the structured grid but also the unstructured and nonconforming grid.

where ω is the angular frequency, μ₀ is the magnetic permeability of the vacuum, and ρ is specific resistance. In FVM, (1) is discretized after volume integration of both sides of Equation (1), and the left hand side is converted to surface integration by integral formula [8].

Where n is the normal vector of each cell surface, ΔS is the Area of each cell surface, Δ and V is the Cell volume. Note that (2) is valid with arbitrary geometry; therefore, if we can properly approximate ∇ × H term on each cell surface, it is possible to use unstructured meshes and local mesh refinement as FEM in Usui et al. For the air layer, (2) should be modified because of poor convergence. They can be achieved by converting and discretizing (1), assuming ρ is constant in the air layer, and applying the condition ∇ . H=0 as follows:

where n=1⁄ρ. (3) is a general poisson equation, which should contribute to the computational stability in the air layer. This program solves (2) and (3) in the subsurface and the air layers, respectively. In (2), ρ on the cell surfaces are needed. This program uses the average as Taylor and Weaver and Mackie and Madden. For computational stability, a transition layer is provided between the air and subsurface. In this layer, the air resistivity is 106 Ω â?? m on the surfaces between the air and subsurface (Figure 2).

Figure 2: Schematic cartoon of a 2D case that indicates where electromagnetic fields are calculated. H denotes magnetic field.

To discretize ∇H in the air layer and ∇ × H in the subsurface layer, ∂H⁄∂xi on the cell surfaces, where i is 1, 2, or 3 and xi is a global coordinate system, should be calculated. This can be done by converting the spatial partial differentiation from the local coordinate system (e₁, e₂, e₃) to the global coordinate system, as follows:

∂H/∂e_i (i=1,2,3) were obtained by interpolation from the surrounding cells.

Boundary conditions on XZ and YZ planes are ∂H⁄∂n, where n is the normal vector. 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 modeling. The author set H_X=1 and H_Y=1 in twice forward modeling to obtain response functions, respectively.

Through these operations, (2) and (3) are discretized, and the linear equation Ax=b (A: Coefficient matrix, x: Unknown magnetic field vector, b: source term) can be assembled. We can obtain H by solving this equation.

To solve the linear equation, this program adopts BiCGSafe [9] with block incomplete LU decomposition as a preconditioner asToh et al. We set the convergence condition using the maximum change of the solution from the previous iteration

ΔΔx_i ^max=max (x_i−x_i−1), where I is iteration number and the residual r_i=Ax_i-b.

Where ε is tolerance.

The divergence correction for subsurface layers is applied as Dong and Egbert, adding a term λ∇ (∇ . H), where λ is the scaling factor, to equation (1). In this program, λ is set to ρ × min (1, kω), where k is the safety factor equal to 0.1 because the large value of λ resulted in poor accuracy solutions, especially in long periods.

Once magnetic field H is obtained, the electric field E is necessary at the cell center to calculate the impedance tensor. The following equation is used for this calculation.

Based on equation (6), the electric fields of two directions parallel to cell surfaces are calculated on each cell surface (Figure 1).

where Ecenter is the required electric field at the cell center, Eparallel is the electric field on the cell surface calculated from equation (6), p is the unit vector parallel to the cell surface, i is unit vector directions (1 or 2), and JJ is cell surfaces number (from 1 to 4). The component of the electric field orthogonal to cell surfaces is excluded because it is not continuous on cell surfaces. From equation (7), eight equations associated with Ecenter in each cell can be obtained. Therefore, Ecenter can be calculated using the least squares method Because Ecenter has three components. Using H and Ecenter, the impedance tensor and magnetic transfer function can be obtained in the same way as in Usui.

Since this program uses the cell-centered FVM, it is necessary to use the surrounding cell-centered 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 is used.

where Z_calc 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, Guennebaud, and other contributor is used for the internal matrix operations. Parallelization is implemented by OpenMP. The linear equations Ax=b at each frequency are solved on each thread.

Inversion: In this program, ASM is implemented for MT inversion. The details are described in Plessix.

In MT inversion, an objective function F=F(m), where m ∈ℜ^N model is the vector of model parameter, and Nmodel is the number of model parameters, should be minimized. In ASM, the objective function F is considered as whose independent variables are m and x ∈ C³Ncell, where x is the solution of linear equation Ax=b in forward modeling and Ncell is the total number of computation cells, that is,

Next, the adjoint equation is solved about λ.

Note that (10) is the same scale equation as Ax=b. Hence, the computational cost is roughly the same as solving Ax=b. The gradient vector of the objective function can be obtained using λ as follows:

where i is from 1 to N_model. Note that ∂F (m, x)/∂mi and ∂F (m, x)^∗/ ∂x are calculated by automatic differentiation using library kv. The objective function F is set as below in this program.

where d_obs is the observed data vector, d_calc is calculated data vector, W is weight matrix, R is smoothness matrix for constraint term, which is necessary in ill-posed MT inversion, N_site is the total number of observation sites where impedance tensor is obtained, Dⁿ ∈ C^{2 × 2} is distortion tensor in each observation site n as described in Ogawa, I is identity matrix, α² and β² are trade-off parameters.

 

To minimize the objective function F with the gradient vector, Adam is implemented, which can be switched easily because the library OptimLib is used internally.

Results and Discussions

Numerical examples

Forward modelling considering topography: To evaluate if the developed program can consider topography, the trapezoidal hill model, which is often used to validate the effect of topography Figure 3 was adopted.

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

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 106 and 100 Ω m respectively. The calculated frequency was 2 Hz. These conditions are the same as previous works. The computer we used was Precision 3630 Tower (CPU: Intel® Core™ i7-9700 CPU 3.00 GHz, RAM:16 GB).

The results are shown in Figure 4. The results are compared to Zhu et al. 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. adopted FEM, and the size of cells. However, since the differences are quite small, it would not influence when the inversion is performed.

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.

Forward modelling considering resistivity heterogeneity: As a test case for heterogeneity of subsurface resistivity, the model in Mackie et al. 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 and the subsurface was 106 Ωm. The compared frequency was 10 and 1000 s. The computer used is the same one used in the Trapezoidal hill model.

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

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. in Figures 6 and 7. While there are some small differences, the results are in good agreement with those of Mackie et al.

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 diamond denote our results. Black dashed lines represent the result in Mackie et al.

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 diamond denote our results. Black dashed lines represent the result in Mackie et al.

Inversion for Dublin secret model 2: The author performed inversion using Dublin Secret Model 2 (DSM2) in Miensopust et al. This model was used as a comparison between different programs. The synthetic observed data of impedance tensor are openly available in Miensopust et al. The author used them as observed data.

The total number of cells was 223,784. The number of cells that were inverted as model parameters is 187,172. The number of degrees of freedom was 671,352. The minimum size of cells was 2000 m × 2000 m ×100 m. The initial resistivity was 100 Ωm. The total number of calculated periods was 10, from 0.025 to 1000 s. The author 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, β² was set to the same value as α². The computer for this calculation was Dell XPS 8950 (CPU: Intel® Core™ i7- 12700 K 3.60 GHz, RAM:64 GB).

The results of the inverted resistivity structure compared to the true model are shown in Figure 8. The author adopted the model when α² was 0.1, judging from the L-curved figure between RMS and model roughness (Figure 9) as in Usui et al. The resistive and conductive anomalies shallower than 10 km are inverted. While the conductive anomaly is connected to the one deeper than 30 km, this would be caused by low sensitivity under the conductive anomaly, as mentioned in Usui. In the conductive layer deeper than 30 km, the value of our result is a little high. This may be improved when data in longer periods are used.

Figure 8: Each cross-section of the true DSM2 model in Miensopust et al. 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 α², respectively. We chose α²=0.1 (red one) as our final result.

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

The total calculation time was about 20 hours. Usui 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 without MPI parallelization, while one in Usui did with it. Nevertheless, we could obtain the results at a comparable time as Usui, although the conditions of computation, such as the computers and the total number of calculated frequencies, were different. Further, our calculation was conducted on an ordinary personal computer. They would suggest the efficient performance of our program.

Conclusion

The author developed a program for 3D MT inversion. For forward modeling and inversion, FVM, in which topography can be considered smoothly, and ASM, in which the gradient vector of the objective function can be efficiently obtained, were adopted, respectively.

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

For the test model of inversion, the author 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 efficient calculations.

Acknowledgement

The authors report there are no competing interests to declare. The data supporting the findings of this study is available from the corresponding author, upon reasonable request.

References

1. Dong H, Egbert GD (2019) Divergence-free solutions to electromagnetic forward and adjoint problems: A regularization approach. Geophys J Int 216: 906–918.

   [Crossref] [Google Scholar]

2. Fujiwara M, Yoshida M, Fujino S (2006) BiCGSafe method with crout version of ILU decomposition giving triple safe keys for convergence. IPSJ Trans Adv Comput Sys 47: 52–60.
3. Groom RW, Bailey RC (1989) Decomposition of magnetotelluric impedance tensors in the presence of local three-dimensional Galvanic distortion. J Geophys Res 94: 1913–1925.

   [Crossref] [Google Scholar]

4. Jahandari H, Farquharson CG (2015) Finite-volume modelling of geophysical electromagnetic data on unstructured grids using potentials. Geophys J Int 202: 1859–1876.

   [Crossref] [Google Scholar]

5. Kelbert A, Meqbel N, Egbert GD, Tandon K (2014) ModEM: A modular system for inversion of electromagnetic geophysical data. Comput Geosci 66: 40–53.

   [Crossref] [Google Scholar]

6. Key KW, Constable SC, Weiss CJ (2006) Mapping 3D salt using the 2D marine magnetotelluric method: Case study from Gemini prospect, Gulf of Mexico. Geophys 71.

   [Crossref] [Google Scholar]

7. Miensopust MP, Queralt P, Jones AG (2013) Magnetotelluric 3-D inversion-a review of two successful workshops on forward and inversion code testing and comparison. Geophys J Int 193: 1216–1238.

   [Crossref] [Google Scholar]

8. Mitsuhata Y, Uchida T (2004) 3D Magnetotelluric modeling using the T-Ω finite-element method. Geophys 69: 108–119.

   [Crossref] [Google Scholar]

9. Nam MJ, Kim HJ, Song Y, Lee TJ, Son JS, et al. (2007) 3D Magnetotelluric Modelling including surface topography. Geophys Prospect 55: 277–287.

   [Crossref] [Google Scholar]