Two-dimensional magnetotelluric inversion using unstructured triangular mesh implemented in Julia

This study presents a 2-D magnetotelluric inversion code tailored for unstructured triangular meshes, developed using Julia, a high-level, high-performance programming language designed for scientific and numerical computation. The forward modeling engine utilizes a node-based finite element method to solve electromagnetic fields throughout the modeling domain. The inversion process employs the Gauss-Newton optimization algorithm, iteratively updating the model via the minimization of a regularized least-squares objective function. We verified the accuracy of the forward modeling using two reference models and performed a synthetic inversion experiment to validate our inversion method and stress the necessity of accurately handling topography-influenced data. Through its application to constructing an electrical resistivity model beneath the northwestern end of Sumatra Island, Indonesia, we demonstrate the practicality of our code for inverting field datasets.


Introduction
The magnetotelluric (MT) method is a passive method in geophysical exploration for probing the electrical resistivity structure of the earth using natural electromagnetic (EM) fields [1] [2].Due to the wide spectrum of the EM source fields, the MT method can probe depths spanning hundreds of meters to hundreds of kilometers.As a result, MT has been used for various purposes, such as local-scale geothermal investigation in volcanic areas (e.g., [3][4] [5][6]) and regional-scale imaging of tectonically active areas (e.g., [7][8] [9]).
One of the challenges in MT modeling is the topographic effect [10].When the observation stations are located on uneven terrain instead of a flat surface, the observed MT data would be distorted due to anomalous EM field response to the topography [11].False treatment of such distorted data, e.g., treating them as if the stations were placed on a flat surface, would lead to misinterpretation of the resistivity structure.A direct and effective approach to tackle this issue is to include the topography directly into the computational mesh of the inversion (e.g., [12]).This strategy ensures that the model response combines contributions from both the resistivity distribution and the topographic effect.It can be accomplished by implementing inversion schemes with grids capable of precisely representing the topography, however complicated it might be.
The finite element method (FEM) is a technique used to solve differential equations and has found applications in various fields, including electromagnetics [13] [14].One of the significant strengths of FEM lies in its capability to handle intricate shapes.This is accomplished by utilizing unstructured meshes composed of irregular elements.These unstructured meshes are particularly advantageous for accurately representing curved objects, which pose challenges when employing Cartesian grids used by the finite difference method.Furthermore, unstructured meshes offer the flexibility to achieve higher local resolution, enabling the precise depiction of fine geometric details and rapid changes in the solution.In the context of MT modeling, this approach allows us to accurately represent topography and apply grid refinements in specific areas, such as around observation stations.Adjusting the element sizes based on data resolution becomes possible, thereby minimizing the degree of freedom.
In this study, we developed a two-dimensional (2-D) MT inversion code using unstructured triangular elements.The code was written in Julia, a high-level, open-source programming language designed for numerical and scientific computing [15].Notable for its performance and ease of use, Julia combines the simplicity of dynamic languages with the speed of compiled languages.It features a just-in-time (JIT) compiler that generates optimized machine code, enabling rapid execution of mathematical operations and data analysis tasks.
The following sections first introduce the forward modeling and inversion methods, succeeded by the implementation in Julia.Subsequently, we evaluated the accuracy of our forward modeling engine using two benchmark models.Following that, a numerical experiment of synthetic inversion was conducted to validate the functionality of our inversion algorithm and emphasize the importance of treating MT data affected by topography correctly.Finally, we used the code to interpret a field dataset acquired in the northwestern end of Sumatra Island, Indonesia.

Forward modeling
In 2-D earth, the governing electromagnetic induction equations can be decoupled into the TE and TM modes [1] [2].The TE (transverse electric) mode consists of an electric field component parallel to the strike, a magnetic field component perpendicular to the strike, and a vertical magnetic field component.The TM (transverse magnetic) mode, on the other hand, consists of an electric field component perpendicular to the strike, a magnetic field component parallel to the strike, and a vertical electric field.The strike direction is the direction where electrical resistivity is constant.
Generally, the 2-D magnetotelluric forward modeling problem seeks a solution to the following differential equation [2], where , , and  are scalar-valued functions whose explicit form depends on the modes.The goal is to recover the field function  throughout the modeling domain.In TE mode,  is the inverse of vacuum permeability ( ),  is the time dependency of the EM field () times the electrical conductivity (), and  is the strike-aligned electric field.Suppose the -axis is the strike direction, then  is  .So, the governing equation to be solved in TE mode has the form of In TM mode,  is the electrical resistivity (),  is  , and  is the strike-aligned magnetic field ( ).
The governing equation to be solved in the TM mode has the form of We used the node-based FEM to solve the differential equations in eqs (2) and (3) numerically.The modeling domain is discretized using unstructured 2-D triangular elements where a unique resistivity value is assigned for each triangular element.The following briefly describes how the differential equation is solved, but for details, interested readers can refer to textbooks on FEM for electromagnetic problems (e.g., [13] [14]).
Here, for simplicity, we derive the formulation with respect to the governing form in eq (1).Multiplying eq (1) by a scalar weighting function  and integrating it over the modeling domain Ω gives Using the identity ∇ •  = ∇ •  − ∇, where  is a scalar field, and  is a vector field, With respect to the second term on the left-hand side of eq (5), Gauss' theorem for surface integral states that where the right-hand side is a line integral over the boundaries of Ω.Along the boundaries of Ω, the solution is completely known.Dirichlet boundary condition is imposed for the whole boundaries.At the top, we apply a uniform value representing a uniform source field.For the side boundaries, we analytically solve the fields as in the case of one-dimensional earth.The bottom boundary is forced to be zero.Because of this, the weighting function vanishes at the boundaries, and the boundary integral in eq (6) equals zero.Therefore, Eq (7) is called the weak form of eq (1).We seek the approximate solutions of  that are expanded in the basis functions , where  is the number of nodes.We used the piecewise linear, or nodal, basis functions in 2-D space.Then, we adopted the Galerkin method, that is, setting the weighting function in eq (7) to be the basis function.After the substitution of eq (8) to eq (7), an  -dimensional linear equation of the form can be obtained, where  is the vector of unknown fields ,  is a null vector, and the element of matrix  at the th row and th column is: After solving eq (9) for , or the field value at all interior nodes, response functions at the observation stations can be calculated.The code can calculate the MT impedance, vertical magnetic field transfer function (or the Tipper), and inter-station horizontal magnetic field transfer function.Generally, a response function at the th observation station ( ) can be calculated as The vectors  and  depend on the response function.As an example, for MT impedance in TE mode, the numerator and denominator on the right-hand side of eq (11) are  and  at the th station, respectively.In this case,  consists of coefficients to perform interpolations of , whereas  consists of coefficients to perform numerical derivatives and interpolations of .

Inversion scheme
In the inversion, a regularized least-squares objective function as the following is minimized: The right-hand side of eq (12) consists of a data misfit term and a model roughness term. is an dimensional vector of observed data where  is the number of data. is a diagonal weighting matrix consisting of the inverse of error variance.The model, , is an  -dimensional vector of the logresistivities where  is the number of triangular elements with variable resistivity value. is a forward modeling functional so that () is an  -dimensional vector of calculated data based on the model .The trade-off parameter  controls the balance between the data misfit and model roughness terms.
The operator  stabilizes the inversion by measuring the model variations, thereby minimizing it will deter spurious structures in the model.Unlike in the cartesian grid, defining the form of  in unstructured triangular grid is not straightforward.After numerous trials, we found that the following works robustly: is the number of neighboring elements sharing an edge with the th-element.It equals three except for elements along the boundaries, along the air-earth interface, or along the surface of fixed-resistivity bodies.The weight  allows unbalanced roughness in horizontal and vertical directions.Setting  smaller (greater) than one implies higher (lower) penalty in the horizontal than the vertical direction.
Adopting the Gauss-Newton algorithm for minimizing the objective function, the model at thiteration can be expressed as The model increment  , which is based on a linearization of eq (12) using the first-order Taylor expansion, is expressed as Eq (15) involves the inversion of  -dimensional matrix, meaning that the computational cost depends heavily on the mesh.In practice,  is often smaller than  by one to two orders of magnitude.Thus, for efficiency, we transformed the dimension of the inversion from  to  using the Woodbury matrix identity, as implemented by previous studies (e.g., [16]).We confirmed that this technique reduces memory usage and computation time, all while maintaining accuracy.The calculation of the Jacobian matrix  is also efficient when  <  by exploiting the reciprocity property of the EM field [17] [18].

Implementation in Julia
Julia is an open-source programming language designed for numerical and scientific computing [15].It was developed by computer scientists from Massachusetts Institute of Technology (MIT) and first appeared in 2012.Notable for its performance and ease of use, Julia combines the simplicity of dynamic languages with the speed of compiled languages.It features a just-in-time (JIT) compiler that enables rapid execution of mathematical operations and data analysis tasks.Thus, Julia is suitable for the MT forward and inverse problems.Julia's syntax is now stable, and for this study, we used version 1.7.3.The primary computational burden in forward modeling and inversion is solving the system of linear equations.Performing a direct calculation for the solution was unfeasible in the past due to computational limitations.As a result, iterative solvers were frequently employed to derive an approximate solution.Yet, with the increasing computational capabilities now available, the utilization of direct solvers has become viable (e.g., [19] [20]).Direct solvers assure solution accuracy and are more dependable when dealing with ill-conditioned systems.For this reason, we preferred to use direct solvers over iterative solvers and used those supported in the LinearAlgebra library of Julia.
For each Gauss-Newton iteration, we have to solve the linear equation using the same coefficient matrix, ,  +1 times (once for forward modeling and  times for constructing the Jacobian matrix).The direct solver involves two main phases: the factorization phase and the solving phase.When dealing with a large coefficient matrix, as in 2-D MT inversion, the factorization phase incurs significant computational cost.Once the matrix is factorized, however, the solving phase becomes highly efficient.Thus, it was desirable to factorize matrix  at first.For this purpose, we used a built-in LU factorization function designed for sparse matrices.This function utilizes the UMFPACK library [21].The syntax for the factorization was simply  = (()), where  is the factorization of .Then, with  established, we proceeded with the solving phase using various right-hand side vectors, achieving efficiency in the process.Using the left division operator (\), the syntax appeared as follows:  = \.
Additionally, we must address a linear equation involving the model increment, as in eq (15), once per Gauss-Newton iteration.Unlike the sparse nature of matrix , the inverted matrix in this context, often referred to as the Hessian, tends to be dense.Consequently, we applied the LU factorization function designed for conventional dense matrices to perform the matrix factorization.Following this, we solved the system by utilizing the left division operator.
Thanks to Julia's inherent ability to handle multithreading, the code's execution would be faster on a computer system equipped with multiple CPU threads.Within the code, the execution of 'for' loops is parallelized through the utilization of the @threads macro.This macro facilitates the concurrent execution of 'for' loops in an unspecified order.Furthermore, the allocation of asynchronous tasks into numerous accessible threads is achievable through the application of @spawn.Once these tasks are finished, the results are returned to the main thread using the fetch command.

COMMEMI 2D-3B Model
The first model used to examine the accuracy of the forward computation was Model 2D-3B from the COMMEMI project (Figure 1) [22].The anomalies in the model are surficial with four distinct regions.The difficulty comes from the differentiation of the field at points where lateral conductivity boundaries reach the surface.We designed the mesh using the Delaunay triangulation method in Gmsh, which is a tool for generating unstructured mesh representations of geometric shapes commonly applied in finite element analysis and numerical simulations [23].The total number of triangular elements in both the conducting half-space and the insulating air half-space was 15754.Given that the imaging resolution decreases with depth and lateral distance from the station array, the size of the triangles varies accordingly to save the degree of freedom.The resistivity value of the air layer was 10 Ωm.The placement of the observation stations followed the reference paper [22].The result of the forward modeling of the COMMEMI 2D-3B Model is shown in Figure 2. Apparent resistivity, Tipper ( ), and inter-station magnetic transfer function ( ) at a period of 100 seconds were evaluated.No reference for phase is available in [22] for unknown reasons. analysis referred to the surface field in the left boundary of the modeling domain ( at  = 0 and  = []).The standard deviations are relatively higher for  and  , and so are the normalized RMS, probably because  and  are ratios of auxiliary fields in TE mode.Overall, the RMS misfits between our results and the reference results are less than one, indicating a good accuracy of our modeling scheme.

2-D Trapezoidal Hill Model
To evaluate the accuracy of the code in representing the topographic effect, we referred to the 2-D trapezoidal hill model proposed by [24] for the second testing model.As in Figure 3a, the hypothetical hill is 0.45 km in height.The width of the hill is 0.45 km at the top and 2 km at the base.The resistivity of the conducting earth is 100 Ωm, whereas the resistivity of the air is 10 Ωm.We calculated the response functions at 77 observation stations with 0.04 km spacing along the hill (Figure 3b).The total number of triangular elements in both the conducting and the insulating air half-space was 11856.
Different from our approach, in [24], the sloping resistivity body is accommodated by subdivision of rectangular elements into triangular elements.Yet, we are pleased to observe that the apparent resistivities calculated by our code fit their results (Figure 3c).Deviation of the apparent resistivities from 100 Ωm reflects the topographic effect.The spike in the TM mode apparent resistivity occurs at the points where the slope changes because of the discontinuity of  [24] [25].In the TE mode, the apparent resistivity anomaly is positive over the hill due to the lack of lateral current flow [24].

Synthetic inversion test
Following the validation of the reliability of the forward modeling approach, we conducted a synthetic inversion test to evaluate the inversion algorithm implemented in the code.For the synthetic dataset, we employed a test model featuring a 2-D trapezoidal hill structure with a rectangular conductive anomaly.
The hill is 3 km high.Its width spans 4 km at its highest point and extends to 20 km at its base.The anomaly is 5 km by 10 km, and its upper boundary lies at 3 km depth.The anomaly is 10 Ωm, embedded in a 100 Ωm background.On the surface of the earth, we placed seven hypothetical observation stations with consistent lateral intervals.At each station, apparent resistivity, phase, and Tipper values were computed across 15 representative periods ranging from 0.1 to 316 s, equally spaced logarithmically, amounting to a total data of 630.A 2% Gaussian noise was added to the data.We first inverted the data using the correct topography in the inversion.The starting model replicated the trapezoidal hill configuration of the true model yet featured a uniform 100 Ωm resistivity of the subsurface.The inversion mesh, however, was different from the forward modeling mesh.In the inversion mesh, the total number of model parameters, or triangular elements with variable resistivity value, was 15190.The resistivity of the air layer, which was 10 Ωm, was kept unchanged during the inversion process.
The inversion result depends heavily on the selection of the trade-off parameter  in the objective function.We run the inversion six times with a different  on each run (300, 100, 30, 10, 3, and 1), as in previous studies (e.g., [16]).It is important to note that  must remain constant during the inversion process to ensure the convergence of the Gauss-Newton algorithm.Then, to determine the best result, or the optimum  value that balances the two terms in the objective function, we adopted the L-curve method [26].We found that the result obtained using  = 30 was at the knee point of the L-curve plot, so the corresponding model was chosen as the best model (see Figure 6b).Figure 4 shows the evolution of the objective function, data misfit, model roughness, and RMS misfit of the inversion with  = 30.As is characteristic of the Gauss-Newton algorithm, the objective function saw a significant decrease at initial iterations, followed by slight reduction in the subsequent iterations (e.g., [18][27]).After five iterations, the inversion was terminated as the change of objective function was less than 1%.The final model explains the data with an RMS of 0.97.The observed and calculated data are in good agreement (Figure 5), indicating that our inversion scheme works correctly.
To understand the potential outcome of misinterpreting the topographic effect, we also inverted the data with no topography in the inversion.The starting model of this inversion was a uniform 100 Ωm model with a flat surface.For the mesh, the total number of model parameters was 10447, less than that of the inversion with the hill in the model.As before, we ran the inversion several times with several different .However, during these runs, the RMS misfit never came close to one.So, the model with the smallest RMS possible, which was 1.62, was chosen as the best model (Figure 6c).The final models of the inversions with and without the topography are shown in Figures 6b and 6c, respectively.The model in Figure 6b is close to the true model.The resistivity and position of the conductive anomaly were well resolved.The top and side boundaries of the anomaly exactly fit the true model.The bottom of the conductive anomaly is slightly deeper than the true model, but it is typical of the EM induction imaging methods.In contrast, despite the absence of the hill, the model in Figure 6c is incomparable to the true model.Due to the misinterpretation of the topographic effect in the data, the inversion overestimated the depth and the width of the anomaly.Also, the inversion resulted in artifacts, such as the surficial small-scale structures and high-resistivity bodies on the left and right sides of the conductive anomaly.This simple experiment indicates that the topographic effect in the data must be interpreted by representing the topography accurately in the inversion.

Field data application
Lastly, we confirmed the practicality of the inversion code through the field data inversion.The data were collected at the northwestern end of Sumatra Island, Indonesia [28].It comprises TE and TM modes apparent resistivities and phases at 12 MT stations and 18 periods per station, ranging from 0.01 to 100 s.As in [28], however, we only used TM mode data in the inversion to avoid three-dimensional signatures due to the surrounding ocean mainly affecting TE mode data (see Figure 7a).Also, a 5% error floor was applied to the data.To deal with the static shift in TM apparent resistivity, we loosened the roughness penalty of cells just beneath MT stations, a method proposed by [29].We set  in eq (13) to be 0.1 to allow laterally elongated structures in the model.The starting model of the inversion was a uniform 100 Ωm model with 0.25 Ωm seawater, and the total number of model parameters was 13,648.
We ran the inversion several times with a different trade-off parameter  on each run and plotted the data misfit against the model roughness of the results (the L-curve).The model shown in Figure 7b, which has an RMS of 1.22, was at the knee point of the L-curve plot, so it was chosen as the preferred result.This model was obtained after 12 Gauss-Newton iterations.The starting RMS was as high as 20.4 due to the static shift in apparent resistivity.But, as we allowed the static shift to be represented by the cells under MT stations, the RMS after one and two model updates were reduced greatly to 6.43 and 2.65, respectively.Running the computation using four Intel ® Xeon ® E5-2687W v4 threads, the computation time needed for one Gauss-Newton iteration was about 9 seconds, so the inversion finished in less than two minutes.The memory usage was about 7.0 GB.Agreements of the observed and calculated data are shown in Figure 7c-d.We are pleased to see that the model correctly explains the data.Also, we can clearly identify the correlation between the model and the data, especially the phases.For example, the resistor R2 explains the significantly low phases in the longer period data at northeastern stations.The resistor R1 and the underlying conductor C1 explain the undulating period dependence of phases at southeastern stations.The conductor C2 explains high phases and very low apparent resistivities at stations in the middle of the profile.
This experiment demonstrates the practicality of our inversion code in handling field datasets.But it is also important to acknowledge some challenges we encountered when using unstructured triangular meshes.Firstly, recovering thin surface layers like the sedimentary basins proves to be more challenging with unstructured meshes compared to regular cartesian meshes.This, however, may depend on the form of the regularization term, and we attempted to include  in eq (13) to address this problem.Second, abrupt changes in slopes around MT stations can lead to inaccurate response functions.This is often a result of incorporating extremely high-resolution topography data in the meshing process.Thus, it is important to incorporate the topography as accurately and as smoothly as possible.

Conclusion
We present a 2-D MT inversion scheme tailored for unstructured triangular meshes and its implementation in Julia programming language.The forward modeling engine utilizes a node-based FEM approach to solve EM induction equations throughout the modeling domain.The Gauss-Newton optimization algorithm is used to minimize a regularized least-squares objective function.We verified the accuracy of the modeling against reference models and validated the inversion scheme by a synthetic inversion experiment.We also used our code to invert a 2-D MT dataset in Sumatra, Indonesia.The inversion efficiently generated a 2-D resistivity structure that accurately explains the observations.

Figure 1 .
Figure 1. a Sketch of COMMEMI 2D-3B [22].Y and Z are in km.Inverted triangles denote the MT stations.bMesh representation of the model for calculating the response functions.

Figure 2 .
Figure 2. A comparison of model responses between our forward modeling results and COMMEMI reference results for the test model 2D-3B.Black dots with error bars denote the average value and standard deviation of COMMEMI reference results [22].Red circles denote the responses from our code.

10thFigure 3 .a
Figure 3. a Sketch of the trapezoidal hill model from Wannamaker et al. (1986) [24].b Mesh representation of the model for calculating the response functions.There are 77 MT stations (inverted triangles) with about 0.04 km intervals.c The calculated apparent resistivities at 50 Hz and 2 Hz are compared with those of [24].Thin vertical lines indicate the position of slope changes.

Figure 4 .
Figure 4. Evolution of the objective function, data misfit, model roughness, and RMS misfit for the inversion with topography.As the subsurface structure in the starting model was uniform, the model roughness was zero initially.The final model of the resistivity structure after five iterations is shown in Figure 6b.

Figure 5 .
Figure 5. Pseudo-section comparison between the observed response functions and the response functions calculated from the model in Figure 6b.Apparent resistivity (  ), phase (), and Tipper (  ) data were used in the inversion.The observed and calculated data are in good agreement.

10th 9 Figure 6 .
Figure 6.Numerical experiment of synthetic inversion using a trapezoidal hill model with a conductive anomaly.The synthetic data of apparent resistivity, phase, and Tipper were calculated from the model shown in a.In b and c, we present inversion results with and without the topography in the inversion, respectively.

Figure 7 .
Figure 7. a MT observation stations (inverted triangles) in the northwestern end of Sumatra [28].A green star in the index map denotes the epicenter of the 2004 M9.1 Sumatra-Andaman Earthquake.The cross-section of resistivity structure obtained from the inversion of TM mode data is shown in b.Pseudo-sections of apparent resistivities and phases are shown in c and d, respectively.