Elsevier

Computers & Geosciences

Volume 159, February 2022, 105035
Computers & Geosciences

A hybrid grid-based finite-element approach for three-dimensional magnetotelluric forward modeling in general anisotropic media

https://doi.org/10.1016/j.cageo.2022.105035Get rights and content

Highlights

  • We develop a hybrid grid-based finite element (FE) approach for magnetotelluric anisotropic forward modeling.

  • The coupling between prismatic and tetrahedral elements is included in the hybrid grid A-φ FE system.

  • The prisms have a superior effect on discretizing the region around sharp electrical boundaries.

  • The grid greatly reduces the number of elements required to maintain solution accuracy.

Abstract

Magnetotelluric (MT) forward modeling often requires the consideration of the deviation generated by anisotropic structures to avoid misleading the detection results, which is implemented by the flexible unstructured finite-element (FE) method based on tetrahedrons. However, the unstructured FE method needs a large number of small elements near the sharp boundary of the electrical structure to produce sufficiently stable grids, and this leads to a sharp increase in the degrees of freedom (DoFs) of the FE system. To this end, we develop an FE approach based on a hybrid grid. It uses the stretched prismatic elements to divide the regions around the abrupt interfaces, which can significantly reduce the number of elements needed to capture the corresponding changes of the physical fields compared to using tetrahedrons. The remaining regions are still divided with tetrahedrons in this hybrid grid in order to keep the total elements needed for the entire computational region at a minimum. By showing numerical examples of the sea- and land-based electrical anisotropy models, as well as a real anisotropic inversion model, the advantages of our method are tested by evaluating the MT response curves, the number of elements, the DoFs, the computational time, and the consumed memory. The results show that the prismatic elements are suitable for discretizing the near-surface region, seawater layer, and electrical anisotropic blocks. This approach can maintain the high accuracy of the numerical solution and reduce the number of elements required for discretizing these regions and the DoFs, thus requiring less computer memory. With this hybrid grid, the computational efficiency of the FE method can be improved for MT forward modeling for a complicated model, which combined with the lower memory consumption is very suitable to implement inversion.

Introduction

In magnetotelluric (MT) exploration, underground geoelectrical media are often accompanied by electrical anisotropy, which often occurs in fault zones, fracture zones, mineral structures, submarine sedimentary rocks, and other geological structures (Jones, 2012; Li and Pek, 2008; Pek and Santos, 2002). If the anisotropy of a given electrical medium is not considered in MT forward modeling, the forward calculation result can produce non-negligible deviations. Consequently, the inversion result produces false anomalies, thereby yielding misleading MT detection results and geological interpretation results (Adetunji et al., 2015; Eisel and Haak, 1999; Martí, 2014; Yin, 2003).

In recent years, many studies have been focused on the influence of electrical anisotropy on electromagnetic induction signals. The main methods of forward modeling are the finite difference method (Kong et al., 2018; Newman and Alumbaugh, 2002; Pek and Verner, 1997) and finite-element (FE) method (Gallardo-Romero and Ruiz-Aguilar, 2022; Li, 2002; Puzyrev et al., 2016; Xiao et al., 2018, 2019; Zhang et al., 2021). The element discretization process of the FE method is the most flexible, and it is highly suitable for undulating terrain and anomalous bodies in complex geoelectric structures. However, the FE method can easily produce a large number of elements (NE), which leads to very large degrees of freedom (DoFs) for the FE system equation and reduces the computational efficiency of the algorithm (Grayver and Bürg, 2014).

To reduce the DoFs and improve the computational efficiency of forward modeling, adaptive unstructured tetrahedral FE algorithm is proposed, which can adaptively refine the elements that have the greatest impacts on the accuracy of the numerical solution (Key and Weiss, 2006; Li et al., 2013; Li and Dai, 2011; Li and Pek, 2008; Liu et al., 2018). A large number of numerical examples show that the adaptive process mainly refines the grid elements of the near-surface region and the anomalous body (Key, 2016; Key and Ovall, 2011). Particularly in the near-surface region, the enrichment of small cells makes the DoF of the FE system matrix increase sharply (Cherevatova et al., 2018; Liu et al., 2018; Ren et al., 2013).

The sharp increase of the NE is mainly due to the limitation of the element quality controlled by the aspect ratio, and hexahedral mesh can be used to solve this problem (Ainsworth and Coggins, 2000; Cherevatova et al., 2018; Javidinejad, 2012; Wang et al., 2013). However, the hexahedron used in these former studies is regular and loses the ability of a free tetrahedron to flexibly discretize the complex underground electrical structure (Grayver and Kolev, 2015). To retain the flexibility of the used grid as well as reduce the number of the total required elements, this study uses prismatic elements to discretize the near-surface region around the sharp air-Earth interface and free tetrahedrons to discretize the others. On the condition that the solution accuracy is guaranteed, this hybrid grid requires fewer elements both around the sharp boundary than the purely tetrahedral grid and in the remaining homogeneous region than the purely hexahedral grid, resulting in fewer total elements and smaller DoFs. Thus, the computational cost is supposed to be lower and will be verified by the numerical examples.

In the subsequent part of this article, Section 2 describes the governing equations and response calculations for MT forward modeling. Section 3 describes how to build hybrid grids and the corresponding FE system equations. In Sections 4 Numerical example of a layered model, 5 Numerical example of an outcrop block model, 6 Numerical example of a buried block model, 7 Numerical example of a real anisotropic inversion model, we use the layered model, the simple block anomaly model, the complex ocean anomaly model, and a real anisotropic inversion model to illustrate the effectiveness of a set of hybrid grids constructed by prismatic elements and free tetrahedral elements. The conclusion is given in Section 8.

Section snippets

MT modeling

Assuming that the time-harmonic factor is eiwt, the quasi-static electromagnetic field governing equations of MT modeling are:×E+iωμH=0×HσE=0andσ=(σxxσxyσxzσyxσyyσyzσzxσzyσzz)is the conductivity tensor in general anisotropic media. ω is the angular frequency, and μ is the permeability in free space. The conductivity tensor can be obtained from the conductivity (σx, σy, σz) in the directions of the three principal electrical axes and the rotations of the principal axes (Bai et al., 2021;

FE theory

The consideration of a hybrid grid can still follow the traditional FE theory. In FE, solving the governing equation of potential is a problem that involves solving partial differential equations (PDEs), and the residual weights of the PDEs can be minimized by using the weighted residual method to solve the system of equations. Using the FE method, Equation (7) and Equation (8) can be discretized as:Ω(×W)(×A˜)dΩ+iωμΩσWA˜dΩ+μΩσWφ˜dΩ+ΩWλdΩ=0iωΩvσA˜dΩ+Ωvσφ˜dΩ=0where A˜, φ

Numerical example of a layered model

A three-layer electrical anisotropic ocean model (Fig. 2) is used to verify the effectiveness of the proposed algorithm and illustrate the advantages of prismatic elements for the near-surface subdivision. The first layer of the layered model is a seawater layer with a thickness of 500 m and a resistivity of 0.3 Ω m. The second layer is an electrically isotropic medium with a thickness of 10 km and a resistivity of 100 Ω m. The third layer is an infinitely downwardly extending electrical

Numerical example of an outcrop block model

A three-dimensional electrical inhomogeneity model containing an electrical anisotropic block (Fig. 4) is used to illustrate the applicability of prismatic elements to discretize electrical anisotropy. The three-dimensional size of the model is 60 × 60 × 60 km, and the surrounding rock resistivity is 100 Ω m. The three-dimensional size of the anomaly is 20 × 10 × 10 km, the axial resistivities are 10/100/10 Ω m, and the anisotropy strike angle around the z-axis is 20°. The anomalous body has no

Numerical example of a buried block model

In a complex three-dimensional ocean model, an electrical anisotropic block anomaly is buried in the electrical isotropic medium below the seafloor (Fig. 7). The three-dimensional size of the model is 50 × 50 × 50 km, and the resistivity of the surrounding seabed rock is 100 Ω m. The three-dimensional size of the anomaly is 5 × 5 × 4 km, the top is 500 m away from the seafloor, and the axial resistivities are 5/5/2 Ω m. The thickness of the seawater layer is 50 m, and its resistivity is

Numerical example of a real anisotropic inversion model

A real anisotropic inversion model is obtained from Kong et al. (2021) to validate the effectiveness of the proposed hybrid grid FE method. Fig. 9 shows some horizontal and vertical slices of the axial resistivities from this anisotropic inversion model. Five frequencies 50, 10, 5, 1, and 0.1 Hz were considered in the anisotropic inversion. The finite-difference hexahedral grid size of Kong et al. (2021) is 43 × 43 × 43 (including 10 air layers). For the hybrid grid discretization, the

Conclusion

We develop and report a set of hybrid grid-based FE forward modeling algorithms for electrical anisotropic MT forward modeling. Prismatic elements and tetrahedral elements are used in this set of hybrid grids. Prismatic elements can be triangular prisms or hexahedral prisms. The bottom face of a triangular prism element can be coupled with a tetrahedron through a common triangular face, and the bottom face of a hexahedral prism element is coupled with the tetrahedron through a pyramid element.

Code availability section

All the data source codes for the models are freely available from the GitHub repository at the address https://github.com/yunianCQU/source-code, or by direct request to the corresponding author ([email protected]). In its current implementation, the programs use MATLAB and COMSOL.

CRediT authorship contribution statement

Nian Yu: Writing – review & editing, Supervision, Project administration. Ruiheng Li: Writing – original draft, Methodology. Wenxin Kong: Writing – review & editing. Lei Gao: Software. Xialan Wu: Data curation. Enci Wang: Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Many thanks are given to the two anonymous reviewers for their helpful and constructive comments. This research was jointly supported in part by the National Natural Science Foundation of China (42074081, 91755215, 42104072), the Science Foundation of Chongqing (cstc2021jcyj-jqX0018, cstc2021ycjh-bgzxm0102).

References (32)

  • A.Q. Adetunji et al.

    Reexamination of magnetotelluric responses and electrical anisotropy of the lithospheric mantle in the Grenville Province, Canada

    J. Geophys. Res. Solid Earth

    (2015)
  • M. Ainsworth et al.

    The stability of mixed hp-finite element methods for Stokes flow on high aspect ratio elements

    SIAM J. Numer. Anal.

    (2000)
  • M. Cherevatova et al.

    A multi-resolution approach to electromagnetic modelling

    Geophys. J. Int.

    (2018)
  • M. Eisel et al.

    Macro-anisotropy of the electrical conductivity of the crust: a magnetotelluric study of the German Continental Deep Drilling site (KTB)

    Geophys. J. Int.

    (1999)
  • A.V. Grayver et al.

    Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method

    Geophys. J. Int.

    (2014)
  • A.V. Grayver et al.

    Large-scale 3D geoelectromagnetic modeling using parallel adaptive high-order finite element method

    Geophysics

    (2015)
  • Cited by (0)

    View full text