Modular finite volume approach for 3D magnetotelluric modeling of the Earth medium with general anisotropy

Highlights

• A finite volume approach is elaborately formulated for 3D Earth models with general anisotropy.

• This approach compromises advantages of finite element method and finite difference method on structured meshes.

• Modular implementation ensures flexible further developments towards inversion.

• Physical fields of Earth models with electrical anisotropy could be directly explored.

Abstract

Electrical property of the Earth should be best depicted by anisotropic 3D models, while many numerical issues remain to be solved before it becomes practical, such as implementation complexity and computational efficiency. Magnetotellurics (MT) is the primary method to explore electrical property of the deep Earth, and MT forward modeling with 3D anisotropy has been investigated primarily by finite element methods (FEM) and finite difference methods (FDM). Here, we present a new implementation of finite volume method (FVM) accounting for general anisotropy. Theoretically, FVM can combine the advantages of FEM in mesh flexibility and superior accuracy, and the advantages of FDM in straightforward implementation and economic computation load. For demonstration, a modular framework to approximate Maxwell's equations on a structured staggered-grid (mesh flexibility ignored), which solves electric fields tangential to cell edges, is implemented in an object-oriented paradigm. Synthetic examples are presented to verify its accuracy and feasibility, which includes (quasi-) analytic solution of simplified models and numerical solution of a dyke model. A variant of COMMEMI 3D-2 model with general anisotropy is simulated to probe the characteristic behaviors of EM fields under such circumstance. This approach makes the dependence of model operators on anisotropic conductivity parameters explicit and can simplify further sensitivity derivation, and takes a step forward for the feasible development of MT inversion algorithms concerning 3D anisotropy.

Introduction

Electrical anisotropy is a direction-dependent property of the Earth, which should be distinguished from its location-dependent inhomogeneity, while both features are of great importance to learn deep geological processes and tectonic activities. Over the decades, many studies have revealed its prevalence in the Earth (e.g., Kellett et al., 1992; Mareschal et al., 1995; Bahr, 1997; Everett and Constable, 1999; Bahr and Simpson, 2002; Gatzemeier and Moorkamp, 2005; Hamilton et al., 2006; Newman et al., 2010; Roux et al., 2011; Zhang et al., 2014; Adetunji et al., 2015). The state-of-the-art research on electrical anisotropy has been reviewed elaborately by related practitioners (cf., Wannamaker, 2005; Baba, 2005; Martí, 2014; Liu et al., 2018b). From the perspective of physics, 3D models with anisotropy in different scales would be the most exhaustive and appropriate to delineate Earth's conductivity structures.

With regard to numerical aspects of 3D anisotropic modeling for geophysical electrical or electromagnetic (EM) methods, Xiong et al. (1986) and Xiong (1989) presented integral equation approaches for dipole source responses of 3D structures embedded in layered anisotropic media, which, to our knowledge, are the earliest reports on numerical calculation of 3D anisotropic geo-electromagnetic responses. A lot of efforts have been dedicated to related disciplines, such as EM induction logging for gas hydrates, DC resistivity method, and marine controlled-source EM method, and publications are continuously emerging (e.g., Davydycheva and Druskin, 1999; Wang and Fang, 2001; Weiss and Newman, 2002; Newman and Alumbaugh, 2002; Li and Spitzer, 2005; Wang et al., 2013; Yin et al., 2014; Liu and Yin, 2014; Cai et al., 2014; Yin et al., 2016; Qiu et al., 2018). Pain et al. (2003) even demonstrated a general framework of 3D anisotropic inversion for DC resistivity data, which shed lights on the anisotropic inversion of other electrical/EM methods.

Focusing on MT specifically, Martinelli and Osella (1997) used a Rayleigh-Fourier method to approximate the 3D MT responses of a multi-layered Earth with irregular interfaces, which is similar to integral equation method in better adapting to localized anomaly. Weidelt (1999) derived a 3D staggered-grid finite difference/volume algorithm solving electric fields perpendicular to cell faces, which accounts for general anisotropy but needs extra efforts to deal with the discontinuity of electric fields. And recently, Löwer and Junge (2017) investigated the influence of 3D anisotropy on the MT transfer functions; Cao et al. (2018a) presented a finite-difference method to model and invert 3D axial anisotropy; a similar implementation by Yu et al. (2018) and Han et al. (2018) accounts for general anisotropy. Xiao et al. (2018), Liu et al. (2018a) and Cao et al. (2018b) respectively reported their finite element algorithms for 3D anisotropy featured by edge-based analysis and various adaptive refinement techniques. Nevertheless, a concise, memory and computationally efficient approach is still in demand for practical use and further development towards inversion.

It's broadly acknowledged that FDM is straightforward to implement, while its accuracy is usually inferior to FEM and FVM; FEM can better adapt to topography and abrupt contrast, but could also be computationally expensive compared to FDM and FVM. In addition, tensor conductivity representing anisotropy may also substantially increase the number of non-zero elements when assembling the FEM coefficient matrix. Thus, theoretically, FVM can combine the advantages of FEM in mesh flexibility and superior accuracy compared to FDM, and the advantages of FDM in straightforward implementation and economic computation load compared to FEM (Jahandari and Farquharson, 2014; Jahandari et al., 2017; Rochlitz et al., 2019).

In light of above-mentioned justifications, along with the facts that physical meaning intuitively originates from the integral representations of differential operators, and straightforward implementation could be partly achieved by using regular meshes, we present a new staggered-grid finite volume approximation of MT responses with general anisotropy for the 3D Earth, which solves electric fields tangential to cell edges and discretizes integral form of frequency-domain Maxwell's equations. The method is implemented in an object-oriented paradigm, for which the simulation procedure can be abstracted into modules that is readily available to the efficient modular derivation of sensitivity calculation and inversion development (Egbert and Kelbert, 2012; Kelbert et al., 2014).

In Section 2, we briefly introduce the governing equations and establish the problem to be numerically solved. Then we elaborate the details about finite volume approach with respect to 3D electrical anisotropy in Section 3, with formal expressions to approximate the Curl-Curl operator and current density. In Section 4, (quasi-) analytic solutions of a half-space and a layered Earth model with general anisotropy are given to verify this approach, and additionally, a dyke model tested by us in 2D case (Guo et al., 2018) is constructed here in 3D situation for comparison. Finally presented in this Section is a variant of COMMEMI 3D-2 model with general anisotropy to probe the characteristic behaviors of EM fields under such situation.