custEM: Customizable finite-element simulation of complex controlled-source electromagnetic data

Raphael Rochlitz; Nico Skibbe; Thomas Günther

doi:10.1190/geo2018-0208.1

ABSTRACT

We have developed the open-source toolbox custEM (customizable electromagnetic modeling) for the simulation of complex 3D controlled-source electromagnetic (CSEM) problems. It is based on the open-source finite-element library FEniCS, which supports tetrahedral meshes, multiprocessing, higher order polynomials, and anisotropy. We use multiple finite-element approaches to solve the time-harmonic Maxwell equations, which are based on total or secondary electric field and gauged potential formulations. In addition, we develop a secondary magnetic field formulation, showing superior performance if only magnetic fields are required. Using Nédélec basis functions, we robustly incorporate the current density on the edges of the mesh for the total field formulations. The latter enable modeling of CSEM problems taking topography into account. We evaluate semianalytical 1D layered-earth solutions with the pyhed library, supporting arbitrary configurations of dipole or loop sources for secondary field calculations. All system matrices have been modified to be symmetric and solved in parallel with the direct solver MUMPS. Aside from the finite-element kernel, mesh generation, interpolation, and visualization modules have been implemented to simplify and automate the modeling workflow. We prove the capability of custEM, including validation against analytic-solutions, crossvalidation of all implemented approaches, and results for a model with 3D topography with four examples. The object-oriented implementation allows for customizable modifications and additions or to use only submodules designed for special tasks, such as mesh generation or matrix assembly. Therefore, the toolbox is suitable for crossvalidation with other codes and as the basis for developing 3D inversion routines.

numerical solutions, electromagnetic theory, finite-element method, 3D modeling, CSEM

INTRODUCTION

In the past few decades, marine, ground-based, and airborne controlled-source electromagnetics (CSEM) have been undergoing a steady development from simple 1D soundings to complex setups with multiple sources and receivers for large-scale multidimensional investigations. In the project Deep Electromagnetic Soundings for Mineral Exploration (DESMEX), semiairborne surveys have been conducted using multiple CSEM transmitters as well as ground-based and airborne receivers. Such CSEM setups require novel data-processing schemes as well as appropriate modeling and inversion tools taking into account the surface and subsurface geometry. Multidimensional numerical methods for Maxwell’s equations in either the time- or frequency-domain were implemented for approximately 50 years by using either the integral equation (IE), finite difference (FD), finite volume (FV), or the finite-element (FE) method. Disregarding the large number of reported 1D and 2D techniques, we restrict the review to advances in 3D CSEM modeling.

Raiche (1974) introduces the IE method for 3D CSEM problems, which is still in use and improved (Zhdanov et al., 2006; Kruglyakov et al., 2016). Although the implementation effort is nonnegligible, this method can be preferred over FD or FE methods (FEM) for modeling the scattering of a limited number of simple anomalies in a 1D layered background model. In this case, the size of the linear systems of equations to be solved is comparatively small (Avdeev, 2005).

In contrast, the FD method is most straightforward to implement, but the computational costs can become high for regular meshes with staggered grids. Within the past few decades, it was often applied in EM geophysics. Adhidjaja and Hohmann (1989) and Wang and Hohmann (1993) present the first time-domain applications. Commer and Newman (2004) and Streich (2009) report more recent usage in the frequency domain. The recent application of multiresolution techniques by Cherevatova et al. (2018) or using octree meshes, such as, for instance, used with FV by Haber and Heldmann (2007), overcomes discretization issues for structurally complex model areas without a tremendous increase in the number of degrees of freedom (dof).

The FV method (Madsen and Ziolkowski, 1990; Jahandari and Farquharson, 2014) combines the advantages of a straightforward implementation, similar to the FD method, with unstructured meshes, such as FE (Jahandari and Farquharson, 2014). Recent research shows that FV might be superior to FD in terms of computational effort, but inferior in terms of accuracy compared with FEM (Jahandari et al., 2017).

FEM is commonly accepted as being most suited to account for irregular geometries (Avdeev, 2005). Despite that FEM might be computationally most expensive, the method becomes comparatively efficient when unstructured meshes are used. In addition, the EM fields can be easily interpolated from the solution functions defined on the complete computational domain.

Pridmore et al. (1981), Gupta et al. (1989), and Livelybrooks (1993) show results of first FE implementations for CSEM problems on hexahedral meshes. For nearly two decades, FEM has rarely been applied to 3D CSEM modeling for different reasons: First, the implementation of curl-curl problems on nodal elements can lead to so-called spurious solutions that do not account for discontinuities of the electric field components (Börner, 2010; Jin, 2015). Nédélec (1980) introduces the edge elements to overcome this issue. Second, computational costs were prohibitively high for former computer architectures. Third, effective and robust mesh generators for tetrahedral elements became commonly available later on, e.g., by the development of TetGen (Si, 2015) or Gmsh (Geuzaine and Remacle, 2009).

Within the past few years, FEM has been first reconsidered in two dimensions (Key and Weiss, 2006; Kong et al., 2007; Key and Ovall, 2011; Minami and Toh, 2013) and remarkable progress has been made in adapting the FEM for modeling 3D CSEM problems. Nowadays, the separation of total fields (TFs) into primary (background) and secondary fields (SFs) is widely used to improve modeling of complex subsurface geometries and reduce computational costs. The background fields are mostly derived by applying semianalytic formulas for 1D layered-earth models. Overall, three different FE strategies have been successfully applied to solve the systems of equations for time-harmonic EM problems.

Badea et al. (2001) describe the first approach, which is based on Coulomb-gauge potentials $(A - ϕ)$ and enable the computation of SF on nodal elements by using iterative solvers. This approach, referred to as $A_{n}$ ⁠, has been taken up by Stalnaker (2005) and Puzyrev et al. (2013), who present fast and robust iterative solvers.

Schwarzbach (2009) and Schwarzbach et al. (2011) report the first successful application of the SF electric field approach, designated as $E$ ⁠, on tetrahedral elements with Nédélec basis functions. Schwarzbach (2009) accurately describes the handling of issues such as the null-space $(ω \to 0)$ and discusses direct or iterative solvers, boundary conditions, or the source implementation. The same approach was applied using hexahedral elements (da Silva et al., 2012). Um et al. (2013) and Cai et al. (2014) successfully implement iterative solution techniques for the ill-conditioned system of equations on tetrahedral and hexahedral meshes, respectively. Recently, Li et al. (2016) use a TF formulation for $E$ in the frequency domain and conversion to the time domain for efficiently simulating moving-loop setups on realistic 3D geometries. For expressing the current density, they use multiple horizontal electric dipoles (HEDs) as introduced by Chave and Cox (1982) or Streich and Becken (2011), but show only a few details about the implementation. The electric field approach is also frequently used in 3D FE magnetotelluric modeling, which becomes challenging for very low frequencies, although the implementation of plane-wave sources is less complicated compared with CSEM. Farquharson and Miensopust (2011), Ren et al. (2013), Grayver and Bürg (2014), Grayver and Kolev (2015), and Kordy et al. (2017) demonstrate important improvements regarding the null-space issue, the solution of the ill-conditioned sparse linear system of equations or adaptive mesh refinement.

Schwarzbach (2009) not only introduces $E$ ⁠, but also a TF potential approach on mixed Nédélec and nodal elements, for the vector and scalar potentials, respectively. This approach, referred to as $A_{m}$ ⁠, is more stable against the null-space issue, the implementation of sources for the TF approach, and the solution of large systems with iterative solvers. Mukherjee and Everett (2011) successfully apply this concept, also taking changes of the magnetic permeability into account. Ansari (2014) and Ansari and Farquharson (2014) report the implementation of electric and magnetic sources by integrating the contribution of the current density to related elements and a stable iterative solver for $A_{m}$ ⁠. Um et al. (2010) demonstrate the first application to 3D time-domain CSEM modeling by using the TF $E$ formulation with an explicit source term implementation and direct solution techniques on tetrahedral meshes. Börner et al. (2015) present advanced Krylov subspace methods to solve the $E$ system for transient EM data efficiently. Cai et al. (2017) report the implementation of the total electric field approach with inhomogeneous Dirichlet boundary conditions.

Most recently, Castillo-Reyes et al. (2018) introduce the development of PETGEM, the first open-source modeling toolbox for 3D (marine) CSEM problems. It is fully parallelized, uses an edge-based secondary $E$ formulation, and can handle complex geometries with unstructured meshes. Accordingly, there are lots of similarities to customizable electromagnetic modeling (custEM), presented in this work. However, currently only first-order polynomials, isotropic conductivities, and HED sources without surface topography seem to be supported in PETGEM.

Despite that significant progress has been made to model 3D CSEM and MT problems, Miensopust (2017) recently summarizes ongoing challenges in the field of geophysical EM data acquisition, modeling, and inversion. Real dipole or loop sources are often approximated with a single (infinitesimal) HED or a vertical magnetic dipole. This approximation can be incorrect for real applications, depending on the frequencies used and the subsurface conductivity distribution (Streich and Becken, 2011). Resistivity changes beneath finite transmitters are mostly not taken into account (Miensopust, 2017). Most importantly, crosschecking of forward modeling solutions for complicated 3D MT models including topography and bathymetry has revealed huge discrepancies and for marine CSEM, 3D modeling and inversion has been reported rarely (Miensopust, 2017).

The reasons might be that most implementations have been optimized for, or are restricted to, either land, marine, or airborne applications and use either SF or TF formulations. SF formulations are particularly powerful for modeling scattering of 3D anomalies embedded in a horizontally layered-earth background model. This leads to highly accurate solutions in the vicinity of transmitters, but as soon as topography or bathymetry effects occur, errors can become significant (Stalnaker, 2005). TF formulations are necessary to take topography or bathymetry into account, but they are assumed to be inferior in terms of accuracy and computational costs. The lack of open-source developments for complex 3D geophysical EM problems makes it difficult to establish crossvalidation. Miensopust (2017) clearly shows the necessity of enhancing the reliability of numerical solutions and making more codes openly available with a proper documentation, such as the multigeophysics modeling and inversion tools provided by SimPEG (Cockett et al., 2015) and pyGIMLi (Rücker et al., 2017). Important developments for geophysical EM purposes are the P223F suite (Raiche et al., 2007), Dipole1D (Key, 2009), AarhusInv (Auken et al., 2014), empymod (Werthmüller, 2017), and PETGEM. However, none of these tools is capable or suited for modeling arbitrary 3D setups.

For this reason, we developed the custEM toolbox. It is written in the programming language Python and based on the well-established and constantly evolving open-source FE library FEniCS (Logg et al., 2012). FEniCS supports nodal, vector and mixed elements, higher order basis functions, source implementation in terms of current density or primary fields, anisotropy, parallelization, and more. However, deriving equivalent real-valued formulations for complex arithmetics is required.

This paper focuses on the methodology and implementation. We modified and implemented the approaches $E$ ⁠, $A_{m}$ ⁠, and $A_{n}$ as SF and TF formulations with the quasistatic approximation. In addition, we introduce a secondary magnetic field approach on Nédélec-elements, which is designated as $H$ ⁠. The resulting symmetric linear systems of equations accelerate the solution process. For TF formulations, we incorporate the source-current density on dof of the Nédélec-function-space as an alternative to known techniques (Ansari, 2014; Li et al., 2016). For SF formulations, we implemented 1D layered-earth solutions for arbitrarily shaped grounded and ungrounded transmitters. Apart from the FEM framework, meshing tools were implemented that feature automated mesh generation for various modeling setups. Interpolation and visualization tools support presenting and validating the FEM results. First, the different FE approaches for the time-harmonic Maxwell’s equations are introduced. Afterward, implementation details of custEM are described. Four examples are presented for validation and comparison and to demonstrate the capabilities. All results can be reproduced by executing the corresponding python scripts in the attached “examples” directory in the custEM repository.

METHODOLOGY

E-field approach — E

Because CSEM methods use frequencies

< 10 kHz

⁠, the quasistatic approximation of Maxwell’s equations is valid. Combining Ampere’s and Faraday’s law yields the Helmholtz equation in terms of the electric field

\nabla \times μ^{- 1} \nabla \times E + i ω σ E = - i ω j_{e},

(1)

where

E

denotes the total electric field,

ω

is the angular frequency,

σ

is the electric conductivity,

μ

is the magnetic permeability, and

j_{e}

is the source current density. For implementation reasons, only sources in terms of electric current density

j_{e}

are considered for deriving the underlying FE equations. Because superposition is valid for

E

⁠, it can be split into a primary (background) field

E^{0}

and a secondary (anomalous) field

E^{s}

(e.g., Newman and Alumbaugh, 1995; Alumbaugh et al., 1996), corresponding to a split of the conductivity

σ = σ_{0} + Δ σ

⁠. Incorporating these terms into equation , the Helmholtz equation for the secondary electric field is derived (e.g., Grayver and Bürg, 2014)

\nabla \times μ^{- 1} \nabla \times E^{s} + i ω σ E^{s} = - i ω Δ σ E^{0} .

(2)

For a comprehensive instruction on FEM for electromagnetic applications, we refer to Monk (2003) or Jin (2015). The weak FE formulations are solved using the Galerkin method (Schwarzbach, 2009). We only present the final real-valued, symmetric system of equations with the structure

Ax = b

that needs to be solved (Grayver and Bürg, 2014). The symmetric submatrices C and D contain terms including trial and test functions

(ϕ)

on a Nédélec function space for the contributions of the coupled real and imaginary parts

[\begin{matrix} C & - D \\ - D & - C \end{matrix}] [\begin{matrix} E^{r} \\ E^{i} \end{matrix}] = [\begin{matrix} b^{r} \\ b^{i} \end{matrix}],

(3)

with

C_{i j} = \int_{Ω} (μ^{- 1} \nabla \times ϕ_{i}) \cdot (\nabla \times ϕ_{j}) d Ω,

(4)

D_{i j} = \int_{Ω} ω (σ ϕ_{i}) \cdot ϕ_{j} d Ω .

(5)

The system matrix is identical for SF and TF. Analogously, the solutions

E^{r}

and

E^{i}

correspond to the TF

E

and the SF

E^{s}

⁠, respectively. For TF, we assumed the current density to be comprised in the imaginary part of the right-side vector

b^{r} = 0,

(6)

b^{i} = \int_{Ω} ω j_{e} \cdot ϕ_{i} d Ω,

(7)

or the primary fields for SF

b^{r / i} = \int_{Ω} ω (Δ σ E^{0 i / r}) \cdot ϕ_{i} d Ω .

(8)

Nonnormalized complex-valued currents are currently not considered, but adding this feature is possible. The magnetic field

H

is related to

E

by

H = - \frac{1}{i ω μ} \nabla \times E,

(9)

where

H

is calculated on the same Nédélec-function-space as used for

E

by solving the FE problem

\int_{Ω} ϕ_{i} \cdot ϕ_{j} d Ω \cdot H^{r / i} = \pm ω^{- 1} \int_{Ω} (μ^{- 1} \nabla \times E^{i / r}) \cdot ϕ_{j} d Ω .

(10)

$H$ -field approach — H

The Helmholtz equation in terms of magnetic fields

H

reads

\nabla \times σ^{- 1} \nabla \times H + i ω μ H = \nabla \times σ^{- 1} j_{e} .

(11)

Similar to the electric field, the magnetic field

H

can also be split into a primary (background) field

H^{0}

and a secondary (anomalous) field

H^{s}

related to the contributions of

σ_{0}

and

Δ σ

⁠. After rearranging equation , the secondary magnetic field can be calculated by solving

\nabla \times σ^{- 1} \nabla \times H^{s} + i ω μ H^{s} = \nabla \times (σ_{0}^{- 1} - σ^{- 1}) \nabla \times H^{0} .

(12)

Equation is derived in Appendix A. Similar to the

E

-field approach, the resulting linear system of equations reads

[\begin{matrix} C_{H} & - D_{H} \\ - D_{H} & - C_{H} \end{matrix}] [\begin{matrix} H^{r} \\ H^{i} \end{matrix}] = [\begin{matrix} b_{H}^{r} \\ b_{H}^{i} \end{matrix}],

(13)

where

C_{H i j} = \int_{Ω} (σ^{- 1} \nabla \times ϕ_{i}) \cdot (\nabla \times ϕ_{j}) d Ω,

(14)

D_{H i j} = \int_{Ω} ω (μ ϕ_{i}) \cdot ϕ_{j} d Ω .

(15)

Again, the submatrices

C_{H}

and

D_{H}

are equal for SF and TF formulations. The solutions

H^{r}

and

H^{i}

either correspond to

H

or

H^{s}

⁠, respectively. The right side for the TF formulation reads

b_{H}^{r} = 0,

(16)

b_{H}^{i} = \int_{Ω} \nabla \times σ^{- 1} \cdot j_{e} \cdot ϕ_{i} d Ω,

(17)

and for the SF formulation

b_{H}^{r / i} = \int_{Ω} (\nabla \times H^{0 r / i}) \cdot (\nabla \times (σ_{0}^{- 1} - σ^{- 1}) ϕ_{i}) d Ω,

(18)

where

E

can be calculated using Ampere’s law in absence of external sources by

E = \frac{1}{σ} \nabla \times H .

(19)

Note that, alternatively, equation can be expressed in terms of the primary electric field as derived in equation in Appendix A. Vice versa,

H^{0}

can be used for the

E

-field and

A - ϕ

approaches.

The $A - ϕ$ approach on mixed elements — $A_{m}$

Expressing the electric field in terms of a vector potential

A

and a scalar potential

ϕ

is a known approach (Nabighian, 1988) to modify equation into

E = i ω (A - \nabla ϕ),

(20)

μ H = \nabla \times A .

(21)

Applying this method, Ansari and Farquharson (2014) derive the following FE formulation based on potential equations:

\nabla \times \nabla \times A + i ω μ σ A + μ σ \nabla ϕ = μ j_{e},

(22)

i ω \nabla \cdot (σ A) + \nabla \cdot (σ \nabla ϕ) = - \nabla \cdot j_{e} .

(23)

Contributions to

A

are incorporated, identically to the

E

-field approach, by trial and test functions on a Nédélec space, whereas contributions of

ϕ

use functions on a nodal space. The assembled system for the TF approach is larger, but more suited for iterative solvers in contrast to the ill-conditioned system of the

E

-field approach, even for low frequencies (Ansari, 2014). The SF formulation reads (Badea et al., 2001)

\nabla \times \nabla \times A^{s} + i ω μ σ A^{s} + μ σ \nabla ϕ^{s} = i ω μ Δ σ A^{0},

(24)

i ω \nabla \cdot (σ A^{s}) + \nabla \cdot (σ \nabla ϕ^{s}) = - i ω \nabla \cdot (Δ σ A^{0}) .

(25)

Ansari and Farquharson (2014) derive the resulting real-valued system of equations, which we modified to be symmetric

[\begin{matrix} C & - D & G & 0 \\ - D & - C & 0 & - G \\ G^{T} & 0 & 0 & H \\ 0 & - G^{T} & H & 0 \end{matrix}] [\begin{matrix} A^{r} \\ A^{i} \\ ϕ^{r} \\ ϕ^{i} \end{matrix}] = [\begin{matrix} c^{r} \\ c^{i} \\ d^{i} \\ d^{r} \end{matrix}],

(26)

where C and D are identical to

E

before (equations and ). Introducing additional test-functions

η

on a nodal space, the other submatrices are defined by

G_{i k} = \int_{Ω} (\nabla ϕ_{i}) \cdot (σ η_{k}) d Ω,

(27)

H_{k l} = \int_{Ω} ω^{- 1} η_{l} \cdot (σ η_{k}) d Ω .

(28)

Even though

G

is not symmetric, the complete system matrix (equation ) for either total potentials

(A - ϕ)

or secondary potentials

(A^{s} - ϕ^{s})

⁠. The right side reads for the TF

c^{r} = \int_{Ω} j_{e} \cdot ϕ_{i} d Ω,

(29)

c^{i} = 0,

(30)

d^{i} = 0,

(31)

d^{r} = - \int_{Ω} ω^{- 1} (\nabla \cdot j_{e}) \cdot η_{i} d Ω .

(32)

For the SF formulation, there is no contribution from

ϕ^{0}

⁠, which is zero for sources used in CSEM applications (Puzyrev et al., 2013) and the right side reads

c^{r / i} = \int_{Ω} ω Δ σ A^{0 i / r} \cdot ϕ_{i} d Ω,

(33)

d^{i / r} = \mp \int_{Ω} Δ σ A^{0 i / r} \cdot \nabla η_{i} d Ω .

(34)

The $A - ϕ$ approach on nodal elements — $A_{n}$

Badea et al. (2001) derive a formulation for the potential approach based on nodal elements, thus reducing implementation effort. By using vector identities and the Coulomb gauge condition

\nabla \cdot A = 0

⁠, which is fulfilled if zero-Dirichlet boundary conditions are used, equations and can be modified to (Badea et al., 2001; Puzyrev et al., 2013)

\nabla^{2} A^{s} + i ω μ σ A^{s} + μ σ \nabla ϕ^{s}) = - i ω μ Δ σ A^{0},

(35)

i ω μ \nabla \cdot (σ A^{s}) + \nabla \cdot (σ \nabla ϕ^{s}) = - \nabla \cdot (i ω μ Δ σ A^{0}) .

(36)

The curl-curl operator is replaced by the Laplace operator. Therefore, curl-conforming Nédélec-elements are not required anymore. Similar to

A_{m}

on the mixed Nédélec-nodal space, a TF formulation can be derived for equations and :

\nabla^{2} A + i ω μ σ A + μ σ \nabla ϕ = - μ j_{e},

(37)

i ω μ \nabla \cdot (σ A) + \nabla \cdot (σ \nabla ϕ) = - \nabla \cdot j_{e} .

(38)

As before, we derived the following real-valued system of equations:

[\begin{matrix} - P & Q & 0 & 0 & 0 & 0 & 0 & X \\ Q & P & 0 & 0 & 0 & 0 & X & 0 \\ 0 & 0 & - P & Q & 0 & 0 & 0 & Y \\ 0 & 0 & Q & P & 0 & 0 & Y & 0 \\ 0 & 0 & 0 & 0 & - P & Q & 0 & Z \\ 0 & 0 & 0 & 0 & Q & P & Z & 0 \\ 0 & X^{T} & 0 & Y^{T} & 0 & Z^{T} & 0 & R \\ X^{T} & 0 & Y^{T} & 0 & Z^{T} & 0 & R & 0 \end{matrix}] [\begin{matrix} A_{x}^{r} \\ A_{x}^{i} \\ A_{y}^{r} \\ A_{y}^{i} \\ A_{z}^{r} \\ A_{z}^{i} \\ ϕ^{r} \\ ϕ^{i} \end{matrix}] = [\begin{matrix} s_{x}^{r} \\ s_{x}^{i} \\ s_{y}^{r} \\ s_{y}^{i} \\ s_{z}^{r} \\ s_{z}^{i} \\ t^{r} \\ t^{i} \end{matrix}] .

(39)

The matrices P, Q, R and X, Y, Z are

P_{i j} = \int_{Ω} μ^{- 1} (\nabla \cdot η_{i}) \cdot (\nabla \cdot η_{j}) d Ω,

(40)

Q_{i j} = \int_{Ω} ω (σ η_{i}) \cdot η_{j} d Ω,

(41)

R_{k l} = - \int_{Ω} ω (\nabla \cdot σ η_{k}) \cdot (\nabla \cdot η_{l}) d Ω,

(42)

{X, Y, Z}_{i k} = \int_{Ω} ω (σ η_{i}) \cdot \frac{\partial η_{k}}{\partial {x, y, z}} d Ω .

(43)

The right sides are for the TF formulation

s_{{x, y, z}}^{r / i} = \int_{Ω} j_{e} {\hat{e}}_{{x, y, z}} \cdot η_{i} d Ω,

(44)

t^{r / i} = - \int_{Ω} (\nabla \cdot j_{e}) \cdot η_{i} d Ω,

(45)

and for the SF formulation

s_{{x, y, z}}^{r / i} = \int_{Ω} ω Δ σ A_{{x, y, z}}^{0 i / r} \cdot η_{i} d Ω,

(46)

t^{r / i} = - \int_{Ω} ω Δ σ A^{0 i / r} \cdot (\nabla η_{i}) d Ω,

(47)

where

{\hat{e}}_{{x, y, z}}

are the unit vectors toward {x, y, z}. Similar to

A_{m}

⁠, the submatrices X, Y, and Z are not symmetric but the overall system of equations is.

Boundary conditions

We implemented the common homogeneous Dirichlet or Neumann boundary conditions. With

n

being the normal vector on the domain boundary

δ Ω

⁠, these are defined as (Schwarzbach, 2009)

Homogeneous Dirichlet : n \times E = 0,

(48)

Neumann : n \times μ^{- 1} \nabla \times E = 0 .

(49)

Alternatively, inhomogeneous Dirichlet conditions can be used for TF computations. The implementation follows directly the concept of Cai et al. (2017). Improved concepts using absorbing boundary conditions (ABCs) or perfectly matches layers (PMLs) are described by Jin (2015) and Schwarzbach (2009), respectively.

Incorporation of current density

For the total electric field approach on Nédélec spaces, relevant dof for implementing the source currents are located on the edges. If those contain segments of the source wire, the corresponding dof can be regarded as HED with the current density

j_{e}

⁠, which is (Ward and Hohmann, 1988)

j_{e} = I d l δ (x) δ (y) δ (z),

(50)

where

I

is the current,

d l

is the length of the HED or accordingly, the edge length, and

δ

is the Dirac delta function. This expression can be inserted into the source term of the weak formulation (equation )

b^{i} = \int_{Ω} ω I d l δ (x) δ (y) δ (z) \cdot ϕ_{i} d Ω .

(51)

The integral over the Dirac function equals one. Furthermore, the inner product of parallel or antiparallel vectors

(j_{e} ∥ \pm ϕ_{i})

equals

\pm 1

⁠. Hence, every element

b_{k}

of the right side vector

b

can be expressed, after assembly, as

b_{k} = \pm ω I d l .

(52)

This term is zero for nonsource dof, whereas it equals the constant

\pm ω I d l

for all dof belonging to the source. The sign needs to be evaluated with respect to the local directions of the Nédélec functions on corresponding edges. When using Nédélec spaces of higher (⁠

p

th) polynomial order, the source contribution on each edge needs to be distributed on

p

dof, and therefore, equation is to be divided by

p

⁠.

This source-incorporation technique can be also applied to $A_{m}$ ⁠. No modification is necessary for loop sources, but the contribution of the divergence term in equation 32 needs to be incorporated at the two nodes corresponding to the ends of a grounded source in an analogous way.

Incorporation of primary fields

For all SF approaches, background fields are required. In principle, custEM is independent of the algorithm used for calculating background fields. An overview on 1D layered-earth EM field calculations can be found in Key (2009) or Werthmüller (2017) along with a summary of recent developments. In custEM, the fully parallelized Python library pyhed is used for deriving primary fields based on dipole calculations according to Ward and Hohmann (1988) and optimized for irregular discretization. The semianalytical solution for the magnetic field is found using Bessel transformations for HED sources. Hankel factors according to Christensen (1990) are applied to solve the Bessel integrals. The complete source term is discretized by multiple HEDs to allow for arbitrarily shaped sources.

IMPLEMENTATION

custEM

Following the concept of reproducible research (Broggini et al., 2017), custEM is open-source and available under the GNU lesser general public license, also used by FEniCS. The complete code was documented using the Python application programmer’s interface. The documentation is available on ReadtheDocs. We also provide test, example, and tutorial files introducing the usage of custEM. Please see the “Data and materials availability” section for links.

The custEM toolbox is arranged in several submodules in the form of Python classes for different categories of tasks, as presented in Figure ¹. The meshgen submodule is designed for automated high-quality tetrahedral mesh generation using TetGen (Si, 2015) and pyGIMLi (Rücker et al., 2017). In core, the core of each FEM computation, i.e., the MOD class, is implemented as well as solvers import, export, and conversion routines. The fem submodule is used to set up function spaces, define the weak formulations with FEniCS, filling the right-side vector for TF formulations and calling a wrapper for pyhed (or potentially other codes) for primary field calculations. Note that a stable version of the independent development of pyhed is currently located in the repository until it is going to be added to pyGIMLi. Interpolation and visualization tools were implemented in the post submodule. In misc, supporting functions for checking the model parameter consistency, prints and serial/parallel computation switches are implemented. The serial submodule contains functionalities that can (currently) only be conducted in serial because not every computation step is supported in parallel by FEniCS, although all core components are designed for multiprocessing (MPI) using the mpi4py library.

Mesh generation

In most publications, meshes are assumed being generated by appropriate mesh generators. In contrast, the focus is on optimizing the boundary value problem and on solving the resulting system of equations. Nevertheless, the accuracy of the solution and convergence rates of the iterative solvers heavily depend on the underlying mesh (Schwarzbach, 2009). Therefore, it is necessary to carefully generate (tetrahedral) meshes with sufficient quality and dimensions. This is straightforward for simple half-space or layered-earth models, but it becomes challenging if arbitrary CSEM setups are taken into account.

We developed a mesh generation approach, using a similar procedure as those Rücker et al. (2006) and Udphuay et al. (2011) demonstrate for ERT with 3D topography. It allows for automatically designing layered earth models of user-defined quality for land-based, airborne, or marine setups including 3D structures and arbitrary source and receiver configurations. Not only topography, but also varying depths for the subsurface layers can be incorporated. For even more complex geologic models, there are several suited (commercial) tools available, such as, e.g., used by Zehner et al. (2015).

The mesh generation procedure with custEM is illustrated in Figure ². First, surface and subsurface interfaces are meshed in 2D using the triangle algorithm (Shewchuk, 1996). Within the surface, source wires or receiver locations need to be already included (Figure ²). Afterward, a third dimension (⁠ $z$ ⁠) is added to the 2D mesh coordinates (Figure ²), and the height for the surface is adjusted by either interpolating from a digital elevation model or using synthetic topography functions (Figure ²). The surface mesh is extended to the computational domain corners by adding polygons with each two of the corners and all the nodes on the corresponding surface edge. Finally, the closing bottom and top quadrangle faces are added to build two conforming piecewise linear complexes (PLCs) for the air space and the subsurface (Figure ²). For the layered-earth case, the 3D extension procedure is conducted iteratively for each subsurface layer.

Within the subsurface layers, arbitrary anomalies (PLCs) can be added, which are allowed to reach the surface but not to intersect subsurface interfaces (for now) (Figure ²). Domain markers are always automatically included. Only for illustrating the raw mesh (Figure ² and ²), TetGen divides every polygon into triangles. In the final mesh with using quality options for the TetGen call, all faces are well-shaped (Figure ²). The ability to enclose the original mesh with a coarse tetrahedral boundary mesh is important for adjusting the overall domain size to fit the physical parameters. In this case, the tetrahedral element type showed its most significant advantage in the form of a rapid increase of the element size toward the boundary without reducing the mesh quality. Approximately 2%–10% of the prior amount of nodes is needed to increase a finely discretized domain with several layers by a factor of 10–1000, respectively (Figure ²). Code examples for generating meshes can be reviewed in the tutorials and examples directories in the custEM repository.

FEM implementation using the FEniCS framework

The implementation of the FE kernel using FEniCS is automated, robust, parallelized, and straightforward. This can be shown by a few basic commands, covering most of the effort needed for manually implementing the FE kernel.

At first, the python backend dolfin and a mesh need to be imported:

Afterward, function spaces, trial and test functions, as well as boundary conditions need to be defined:

The mixed function space is required for the coupled real-valued systems of equations introduced by equation 3. The greatest benefit of using FEniCS is the ability of readily implementing complicated weak formulations with help of the unified form language, below presented for the total $E$ -field approach

This FEM system can be solved automatically and exported afterward:

The code snippet above includes only a few basic commands, but covers the complete TF $E$ -field approach solution for a homogeneous domain, disregarding missing definitions of the physical parameters and the zero, source, and boundary functions. However, the implementation to solve arbitrary 3D CSEM problems requires the incorporation of multiple physical domains as well as real finite sources that could be achieved by using advanced functionalities of FEniCS.

Source terms and primary fields

Our technique of directly imposing the source current on Nédélec-dof requires incorporating the transmitter wire(s) as edges in the mesh, which are still allowed to differ in length. A crooked loop transmitter, incorporated in the surface-mesh, is presented in Figure ³. This source is used in example 2 subsequently. In principle, it is possible to implement borehole, marine, surface, or airborne sources anywhere in the modeling domain. Considering the transmitter location(s) as defined, the workflow for incorporating $j_{e}$ is structured as follows:

1)
Find all dof with coordinates on the source wire.
2)
For each of these dof, find one element containing the dof.
3)
Use dof and coordinates of 1, elements of 2, and local definitions of basis functions (Touma Holmberg, 1998) to identify the direction (⁠ $\pm 1$ ⁠) of the Nédélec dof(s) within the corresponding elements.
4)
Fill the right side vector according to equation 52 for all dof of step (1).

Because this procedure is completely parallelized, it takes little time compared with solving the system of equations. Primary fields for the SF formulations are calculated on the fly with pyhed for each mesh and parametrization (e.g., source and frequency) on the first run and cached for repeated use because the computational effort is not negligible.

Postprocessing

Once the main solution is calculated, conversion to either secondary or total electric or magnetic fields can be automatically conducted. Generating interpolation meshes in the form of either straight or crooked lines and interfaces (e.g., surface with topography) is automated as well. For interpolation, multithreading is applied, which leads to a rapid decrease of interpolation time. This procedure is suitable because, usually, multiple processes are in use anyway to solve the main problem. The model configuration is stored to make relevant information available for the visualization and record. The post submodule provides plot tools with utility functions for basic comparison and error computation.

EXAMPLES

Example 1: HED on half-space

First, we present a comparison with primary field solutions derived using pyhed (the dashed/dotted lines), Dipole1D (Key, 2009) (the dotted lines), and empymod (Werthmüller, 2017) (the dashed lines) in Figure ⁴. Electric and magnetic fields were calculated for 1, 10, 100, and 1000 Hz on a half-space $(100 Ω m)$ with a $y$ -directed HED source in the origin. The observation line extends $\pm 5 km$ perpendicular to the HED. Results from TF $E$ (the colored lines indicating the frequencies) p1 computations were added in Figure ⁴. The underlying mesh for this FEM solution has a comparatively fewer amount of dof $(approximately 37 k)$ to provide a simple example that requires only a few minutes computation time and $< 8 GB$ RAM with $< 4$ MPI processes. Therefore, high accuracy of the FEM solution cannot be expected.

The pyhed, Dipole1D, and empymod reference solutions match well aside from $I (E_{y})$ for 1000 Hz and $R (H_{z})$ for 100 and 1000 Hz. The differences occur at offsets > $1 km$ ⁠. Here, the $R (H_{z})$ amplitude is already six orders of magnitude smaller than in the vicinity of the HED. These variations are assumed to originate from different Hankel factors used. Even though the FEM mesh is minimalistic, the TF $E$ solutions match well to the semianalytical references. The FEM solution differs significantly from the analytic ones at 1 Hz for $I (H_{x})$ ⁠, which results from boundary effects. Further discrepancies are mainly noticeable at 100 and 1000 Hz for $I (E_{y})$ and $R (H_{z})$ ⁠, which indicates an insufficient discretization.

Example 2: Crooked-loop source on a three-layer earth

We verified the TF $E$ implementation for a three-layer model with a crooked loop (Tx) source because this is one of the most complex problems to be solved semianalytically. Reference solutions obtained with the pyhed module were used for validation.

The uppermost two layers had thicknesses of 300 and 700 m; the resistivities from top to bottom were $10^{3}$ ⁠, $10^{2}$ ⁠, and $10^{4} Ω m$ ⁠. The air-space resistivity was set to $10^{7} Ω m$ because we found that for usual CSEM setups, a conductivity contrast of four orders of magnitude between the airspace and the subsurface is completely sufficient to obtain accurate models. The randomly distorted loop source had a circumference of approximately 4 km (Figure ³) and was located in the center of a $200 \times 200 \times 200 km$ mesh.

In Figure ⁵, the TF $E$ solution at the surface based on second-order polynomials (p2), is depicted component-wise and compared with the 1D solution for a typical base frequency of 10 Hz. The underlying mesh, referred to as a coarse mesh, was parameterized by a maximum facet size of $1000 m^{2}$ within the $4 \times 4 km$ observation area (Figure ³, the red lines) and a TetGen quality of 1.4 in terms of the radius-edge ratio for the elements.

The single vector-components $E_{x}$ ⁠, $E_{y}$ ⁠, $H_{x}$ ⁠, $H_{y}$ ⁠, and $H_{z}$ (Figure ⁵–⁵) of the FEM solution match well with the analytic solution. The overall magnitude errors (Figure ⁵) are smaller than 1%, disregarding the vicinity of Tx. The value $E_{z}$ was not presented because the $E$ -field is discontinuous at the surface. For gaining the first estimates of relations between various frequencies, mesh discretization, polynomial orders, computation times, and memory requirements, we calculated the crooked-loop model at four frequencies (1, 10, 100, and 1000 Hz) for the three different mesh configurations listed in Table 1. The computation times are based on 32 parallel processes used on a Dell PowerEdge R940 server with four Intel Xeon Gold 6154 processors and 48 LRDIMM 64 GB, DDR4-2666, Quad Ranks. First, we used the same coarse mesh for p1, which results in a significantly smaller number of dof. Second, an inter(mediate) mesh with a quality of 1.3 and a 10th of the maximum facet size in the observation area was used. Third, we added a fine mesh (Figure ³, the black lines) with a quality of 1.2 for and a maximum facet size of $50 m^{2}$ in the observation area to evaluate the development of errors with increasing quality and discretization (Table 1). For each mesh, the overall domain size amounts to $200 \times 200 \times 200 km$ ⁠.

Figure ⁶ shows the errors of the p1 results on the coarse and the fine mesh in comparison with the p2 results on the coarse mesh at the surface. Only errors of the real part of $E_{y}$ and imaginary part of $H_{z}$ are shown as these components are representative for the overall results. A systematic misfit (the bluish or reddish areas) can be distinguished from randomly distributed errors (the red-white-blue pattern). Regarding effects and all frequencies, the results for p2 polynomials (Figure ⁶ and ⁶) were significantly more accurate compared with all p1 computations. Furthermore, the p1 results on the fine mesh (Figure ⁶ and ⁶) exhibited smaller randomly distributed errors but equal systematic misfits compared with the coarse mesh (Figure ⁶ and ⁶). Note that the computational effort for p2 on the coarse mesh is even lower than for p1 on the inter mesh (Table 1).

With respect to the different frequencies, the smallest systematic offsets occurred at 100 Hz and increased slightly at 10 Hz for all three configurations (Figure ⁶). At 1 Hz, the overall misfit became significant. At 1000 Hz, the amount of randomly distributed errors was highest. Table 2 contains mean error estimates for real and imaginary parts of $E$ and $H$ for the above configurations at 10 and 1000 Hz. Because of extreme deviations close to Tx, we calculated the mean error only considering data points with a relative error $< 100 %$ ⁠. The resulting number of outliers is listed in Table 2 as well. This corresponds to the quality of the solution in the proximity to Tx. The mean errors of p2 were far below all values of p1. For the latter, the “inter” compared with the “coarse” mesh lead to a remarkable improvement. The “fine” mesh led to a further increase in accuracy, but the computational effort was significantly higher.

Example 3: Finite dipole source on a two-layer earth with 3D anomalies

All implemented approaches were compared on a half-space model with a 1 km dipole Tx and two embedded conductive 3D anomalies, as illustrated and specified in Figure ⁷. The 1 km long brick on one side of the Tx has an anisotropic conductivity structure, and it was shifted by 200 m in the $y$ -direction. Therefore, both anomalies are expected to reveal strong 3D effects.

The computations of all SF approaches using p1 polynomials are based on a mesh (approximately 258 k nodes) with fine discretization for the observation line with 10 m node separation and the anomalies (maximum cell volume: $1000 m^{3}$ ⁠). The p2 mesh uses a maximum cell size of $10,000 m^{3}$ ⁠, resulting in approximately $64 k$ nodes. The computational parameters are listed in Table 3. The complete domain size was $2000 \times 2000 \times 2000 km$ (only approximately 3% more nodes compared with a 1000th of this size), needed to suppress boundary effects for the $H$ and $A_{n}$ approaches. The computation times and memory requirements are based on 32 parallel processes for p1 and p2. The given times refer only to the assembly and solution process of the main systems of equations on the machine already used for example 2.

Figures ⁸ and ⁹ present magnetic field results for 10 Hz. Considering the SF computations in Figure ⁸, all approaches exhibited equal amplitudes for the dominant $H_{z}$ component and only slight mismatches for the slightly weaker component $H_{x}$ ⁠. Minor discrepancies between the four approaches could be observed for $H_{y}$ ⁠, which can be attributed to numerical differences, i.e., nodal or edge functions and interpolation on the observation line. The best results could be achieved with p2 (black line in Figure ⁸), especially for the weak $H_{y}$ component.

We compared the TF $E$ and $A_{m}$ approaches only with the (summed) TFs from the SF $H$ -p2 approach (black line) in Figure ⁹. The match of the obtained SFs was already shown in Figure ⁸. As expected, the mismatch of the TFs is most prominent for $H_{y}$ ⁠. The real part revealed mismatches in the vicinity of the source and the imaginary part above the brick anomaly. As already observed in example 2, p2 for the TF approach led to increased accuracy in form of significantly smoother results for the weak $H_{y}$ component. Another mismatch between the TF and SF computations can be observed at $x = - 3 km$ for $I (H_{x})$ ⁠. In general, mismatches occur only in the nondominant horizontal components, whose amplitudes are at maximum 1% of the vector magnitude for this model setup.

Example 4: Loop Tx, sinusoidal topography, and conductive dike

The capability of custEM with regard to realistic 3D modeling is demonstrated by the fourth example. We computed $E$ and $H$ fields with the TF $E$ approach and p2 for a half-space-like model with sinusoidal-topography and a $1 \times 1 km$ large loop source, placed on the surface in the center of the model. Crossing the Tx at the surface, a crooked conductive dike with varying thickness was incorporated, as illustrated in Figure ¹⁰.

The $E$ and $H$ fields (Figure ¹⁰–¹⁰) were interpolated on the surface and on a level 50 m above the surface, respectively (Figure ¹⁰). The interpolation grids are regular in the horizontal directions and were chosen to simulate a semiairborne survey setup with a transmitter on the surface, ground-based $E$ -field, and airborne $H$ -field receiver locations. The influence of topography and in particular, the effect of the conductive dike, can be observed in all field magnitudes. The complete code required for the mesh generation, FEM solution, and interpolation can be viewed in Appendix B.

DISCUSSION

Reliability of the FEM solutions

Examples 1 and 2 provided a verification of the primary fields and the source term incorporation. The availability of different approaches allowed for crossvalidation of the solutions in example 3 up to a certain extent because the physics and the mathematical description of the original FE problem are differing. This results in unique linear systems of equations after assembly for all approaches. In this context, the overall match for example 3 is remarkable, aside from variations in the nondominant horizontal components. Although only magnetic fields were presented, example 3 implicitly showed the correct implementation of each FE approach because for $E$ ⁠, $A_{m}$ ⁠, and $A_{n}$ ⁠, the $H$ -fields are derived from the main quantities $E$ or $A$ ⁠.

Both examples clearly indicate significant advantages of p2, particularly in the vicinity of the source. Because EM fields decay exponentially, it is plausible that quadratic interpolation functions can approximate the real field behavior significantly better in the near-field zone of the transmitter or scattering anomaly. In contrast, if only the far field is of interest, i.e., for SF formulations without anomalies in the uppermost subsurface, the difference in accuracy is exiguous and using p1 is sufficient.

For p1 and p2, there is no possibility of designing a pair of meshes with identical dof positions. Furthermore, the matrix sparsity patterns, and thus the bandwidth, are different for an equal number of dof. However, the most reasonable way for comparing accuracy and performance might be to use either the same mesh or a similar number of dof as well as equal mesh quality. According to Tables 1 and 2, it can be inferred that p2 is not only superior in terms of accuracy. The computational effort can be significantly lower compared with p1, if certain accuracy for a larger observation area $(4 \times 4 km)$ is required. At 1 Hz, the systematic errors became significant even for p2, which is an indication for an insufficient domain size. At high frequencies, the systematic misfit increased strongly toward the observation area boundaries, which is probably related to the rapid decrease of the field amplitudes within the $4 \times 4 km$ domain.

Using polynomials of third (p3) or even higher order would further increase the accuracy on even coarser meshes. However, the solution of the system matrices requires significantly more resources and time as was tested for a HED in fullspace with p3. We assume that p2 gives the optimum balance between a good approximation of the exponential EM-field decay behavior with quadratic elements and the computational effort. Grayver and Kolev (2015) conclude an identical statement from their investigations.

Further performance analysis and validation of the FEM results calculated with custEM are already in progress. The focus is especially on crosschecking with other codes for CSEM data, e.g., PETGEM.

Performance

The direct solver MUMPS had shown to be robust and efficient for all approaches with a maximum memory requirement of 1000 GB. Alternative direct solvers provided by FEniCS (i.e., UMFPACK and PETSc LU solver) are not appropriate for MPI and do not support symmetric matrices. Solution of the $E$ or $H$ approaches with direct solvers showed the highest performance due to the reduced size and bandwidth of the matrices. A machine with 256 GB RAM is recommended for realistic simulations of models with multiple layers and anomalies.

Regarding computation times and memory requirements, the provided values vary with the number of parallel processes. Even though the computation times decrease with more processes, the scaling becomes worse and the memory requirements increase. For instance, the RAM overhead using 40 instead of 20 processes can become very high (>100 GB) for a secondary $E$ computation with 10 M dof and p1, whereas the computation time is reduced only by $\approx 20 %$ ⁠. Hence, simultaneously running multiple jobs (e.g., for different frequencies) with a smaller number (12–32) of processes performs better than consecutive computations with more parallel processes.

Besides the main problem, there might be a significant overhead in time when initially computing primary fields for SF formulations. Even though this procedure is parallelized, computation times can become even higher than for solving the main system of equations, if the source discretization requires a large number of HEDs (>100) and the subsurface consists of many layers. In addition, postprocessing in the form of converting $E$ to $H$ or vice versa requires 20%–50% of the time needed to solve the main problem. For the potential approaches, conversion to $H$ is fast but deriving $E$ takes even more time. Therefore, the SF $H$ approach can be used to significantly reduce computation times, if no electric fields are required, because neither the calculation of $E^{0}$ nor the postprocessing conversion is necessary anymore. Suitable applications are inversion procedures for airborne EM data or multidimensional nuclear magnetic resonance data.

Potential of custEM and FEniCS

Various combinations of iterative solvers and preconditioners available in FEniCS were tested, but no method showed sufficient convergence rates yet. Robust iterative solution techniques have been applied successfully to solve the $(A - ϕ)$ systems (Puzyrev et al., 2013; Ansari, 2014) and the SF $E$ system of equations (Grayver and Bürg, 2014). We assume that iterative solvers can be applied for custEM as well, if effort is spent on analyzing the assembled system matrices and finding appropriate preconditioners. The use of iterative solvers can result in a reduction of computation time and memory requirements, especially for the potential approaches, which would make them not only valuable for crosschecking solutions obtained with the $E$ or $H$ approaches. Advanced boundary conditions such as ABC or PML might reduce (systematic) errors or can decrease the modeling domain size (Schwarzbach, 2009). Another important improvement in custEM would be the application of adaptive mesh refinement, which is in principle supported by FEniCS.

Aside from frequency-domain CSEM, it stands to reason to adapt custEM for transient electromagnetic modeling, either in the time domain or by using effective transformations on multiple frequency-dependent results. Moreover, the support of magnetotelluric modeling would require a manageable implementation effort on modifying the FE modules. In general, FEniCS is well-suited for multigeophysics modeling. The submodules of custEM can be further used for special tasks, e.g., high-performance assembly of the systems of equations and export of the matrix in the CRS format for custom solvers. Due to the modular design, additional features can be easily implemented for user-specific requirements.

The number of available forward modeling codes for 3D CSEM data is low, but the number of inversion codes is even smaller. The custEM toolbox can be the basis for 3D inversion, particularly because anisotropy or changes in the magnetic permeability can be easily implemented with FEniCS. Following the concept of Günther et al. (2006), the TF formulation could be applied to calculate sufficiently accurate primary fields on a fine-discretized 3D mesh, including topography or bathymetry. Afterward, the computational effort for inversion might be reduced drastically with using SF approaches on a mesh with a lesser number of dof, when only changes of resistivity in the subsurface are of interest.

CONCLUSION

We verified the implementation and provided insights into the capability of custEM. It was shown that four different SF and two TF approaches, based on varying physical and mathematical concepts, led to remarkably consistent results, even though smaller misfits, probably due to interpolation and boundary effects, can be observed. For TF formulations, we successfully incorporated the source current density on the edges of the mesh as an alternative technique to existing concepts. In any event, corresponding advantages and disadvantages should be further investigated. The use of second-order polynomials clearly outperformed p1 in terms of accuracy and computational effort, especially for TF approaches. With p2, sufficient accuracy at the observation points could be achieved by using a comparatively coarse discretization and intermediate mesh quality.

The $E$ - and $H$ -field approaches have revealed highest performance in terms of accuracy and computational costs. The potential approaches were valuable for validating the implementation. We introduced the $H$ -field approach for modeling 3D CSEM problems with FE, which saved 20%–50% of the overall calculation time without the transformation from $E$ to $H$ ⁠, if only magnetic fields are of interest as, e.g., in airborne EM applications.

Aside from the FE modules, mesh generation, interpolation, and visualization tools enable automated workflows. The developed toolbox custEM is freely available and provides not only the ability for customizable CSEM modeling and crossvalidation, but also the potential to be used for inversion or other geophysical methods.

ACKNOWLEDGMENTS

We thank the developers of FEniCS, TetGen, and pyGIMLi for their effort on developing software tools for the community over the years. The project DESMEX was funded by the German Federal Ministry of Education and Research (BMBF) in the framework of the research and development program Fona-r4 under grant 033R130D. We highly appreciate the thorough comments from R.-U. Börner, D. Werthmüller, and an anonymous reviewer that helped to significantly improve the code and manuscript.

DATA AND MATERIALS AVAILABILITY

The custEM toolbox is open-source and available at https://gitlab.com/Rochlitz.R/custEM under the GNU Lesser General Public License (LGPL). The complete code is documented utilizing the Python Application Programmer’s Interface (API). The documentation is available on ReadtheDocs at https://custem.readthedocs.io/en/latest/. Instructions, notes and tutorials are available on this page as well. The results of all presented examples can be reproduced using the provided scripts in the “examples” directory of the custEM repository.

DERIVATION OF THE SECONDARY MAGNETIC FIELD APPROACH

Faraday’s law and Ampere’s law, disregarding magnetic source terms, read

\nabla \times E = - i ω μ H,

(A-1)

\nabla \times H - σ E = j_{e},

(A-2)

for the TF and

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

(A-3)

\nabla \times H^{0} - σ_{0} E^{0} = j_{e}, (A 4)

(A-4)

for the primary field. With

E^{s} = E - E^{0}

⁠,

H^{s} = H - H^{0}

⁠, and

Δ σ = σ - σ_{0}

⁠, subtracting equation from equations and from equation gives

\nabla \times E^{s} = - i ω μ H^{s}, (A 5)

(A-5)

\nabla \times H^{s} - σ E^{s} = Δ σ E^{0} .

(A-6)

Next, multiply

σ^{- 1}

and equation

σ^{- 1} \nabla \times H^{s} - E^{s} = σ^{- 1} Δ σ E^{0}, (A 7)

(A-7)

and take the curl of equation

\nabla \times σ^{- 1} \nabla \times H^{s} - \nabla \times E^{s} = \nabla \times σ^{- 1} Δ σ E^{0} .

(A-8)

Inserting equation into equation leads to the governing equation for the secondary magnetic field approach in terms of primary electric fields

\nabla \times σ^{- 1} \nabla \times H^{s} + i ω μ H^{s} = \nabla \times σ^{- 1} Δ σ E^{0} .

(A-9)

In terms of primary magnetic fields, the following equation can be derived:

\nabla \times σ^{- 1} \nabla \times H^{s} + i ω μ H^{s} = \nabla \times (σ_{0}^{- 1} - σ^{- 1}) \nabla \times H^{0} .

(A-10)

Mesh	$p$	$A_{\max} (m^{2})$	$q$	Nodes (k)	dof (k)	$t$ (min)	RAM (GB)
Coarse	1	1000	1.4	60	858	$< 1$	$\approx 32$
Inter	1	100	1.3	502	7545	$\approx 6$	$\approx 307$
Fine	1	50	1.2	1283	18635	$\approx 29$	$\approx 972$
Coarse	2	1000	1.4	60	4675	$\approx 5$	$\approx 284$

Mesh	$p$	$A_{\max} (m^{2})$	$q$	Nodes (k)	dof (k)	$t$ (min)	RAM (GB)
Coarse	1	1000	1.4	60	858	$< 1$	$\approx 32$
Inter	1	100	1.3	502	7545	$\approx 6$	$\approx 307$
Fine	1	50	1.2	1283	18635	$\approx 29$	$\approx 972$
Coarse	2	1000	1.4	60	4675	$\approx 5$	$\approx 284$

Mesh	$p$	$R (\| E \|)$	$I (\| E \|)$	$n$	$R (\| H \|)$	$n$	$I (\| H \|)$	$f$
Coarse	1	3.53	2.01	39	4.28	669	3.28	at 10 Hz
Inter	1	2.40	1.70	24	1.92	207	1.84	–
Fine	1	1.53	0.76	23	1.34	142	1.21	–
Coarse	2	0.12	0.16	25	0.99	222	0.10	–
Coarse	1	2.70	5.68	43	5.18	657	6.00	at 1000 Hz
Inter	1	1.41	3.04	24	2.08	206	1.74	–
Fine	1	0.68	1.75	25	1.34	143	1.18	–
Coarse	2	0.22	0.30	25	0.90	222	0.38	–

Mesh	$p$	$R (\| E \|)$	$I (\| E \|)$	$n$	$R (\| H \|)$	$n$	$I (\| H \|)$	$f$
Coarse	1	3.53	2.01	39	4.28	669	3.28	at 10 Hz
Inter	1	2.40	1.70	24	1.92	207	1.84	–
Fine	1	1.53	0.76	23	1.34	142	1.21	–
Coarse	2	0.12	0.16	25	0.99	222	0.10	–
Coarse	1	2.70	5.68	43	5.18	657	6.00	at 1000 Hz
Inter	1	1.41	3.04	24	2.08	206	1.74	–
Fine	1	0.68	1.75	25	1.34	143	1.18	–
Coarse	2	0.22	0.30	25	0.90	222	0.38	–

Approach	$P$	dof (k)	$t$ (min)	RAM (GB)
TF $E$	1	3887	$\approx 1.9$	108
$TF A_{m}$	1	4421	$\approx 3.9$	207
SF $E$	1	3887	$\approx 1.9$	124
SF $H$	1	3887	$\approx 1.9$	133
$SF A_{m}$	1	4421	$\approx 4.0$	223
$SF A_{n}$	1	2137	$\approx 1.9$	228
TF $E$	2	5228	$\approx 3.3$	240
SF $E$	2	5228	$\approx 3.4$	250
SF $H$	2	5228	$\approx 3.4$	260

custEM: Customizable finite-element simulation of complex controlled-source electromagnetic data

ABSTRACT

INTRODUCTION

METHODOLOGY

E-field approach — E

$H$ -field approach — H

The $A - ϕ$ approach on mixed elements — $A_{m}$

The $A - ϕ$ approach on nodal elements — $A_{n}$

Boundary conditions

Incorporation of current density

Incorporation of primary fields

IMPLEMENTATION

custEM

Mesh generation

FEM implementation using the FEniCS framework

Source terms and primary fields

Postprocessing

EXAMPLES

Example 1: HED on half-space

Example 2: Crooked-loop source on a three-layer earth

Example 3: Finite dipole source on a two-layer earth with 3D anomalies

Example 4: Loop Tx, sinusoidal topography, and conductive dike

DISCUSSION

Reliability of the FEM solutions

Performance

Potential of custEM and FEniCS

CONCLUSION

ACKNOWLEDGMENTS

DATA AND MATERIALS AVAILABILITY

DERIVATION OF THE SECONDARY MAGNETIC FIELD APPROACH

COMPLETE CODE FOR GENERATING THE EXAMPLE 3 MESH: MESH_EXAMPLE_3.PY

COMPLETE CODE FOR CALCULATING THE EXAMPLE 3: RUN_EXAMPLE_3.PY

Data & Figures

About GSW

Resources

Explore

Connect

custEM: Customizable finite-element simulation of complex controlled-source electromagnetic data

ABSTRACT

INTRODUCTION

METHODOLOGY

E-field approach — E

H-field approach — H

The A−ϕ approach on mixed elements — Am

The A−ϕ approach on nodal elements — An

Boundary conditions

Incorporation of current density

Incorporation of primary fields

IMPLEMENTATION

custEM

Mesh generation

FEM implementation using the FEniCS framework

Source terms and primary fields

Postprocessing

EXAMPLES

Example 1: HED on half-space

Example 2: Crooked-loop source on a three-layer earth

Example 3: Finite dipole source on a two-layer earth with 3D anomalies

Example 4: Loop Tx, sinusoidal topography, and conductive dike

DISCUSSION

Reliability of the FEM solutions

Performance

Potential of custEM and FEniCS

CONCLUSION

ACKNOWLEDGMENTS

DATA AND MATERIALS AVAILABILITY

DERIVATION OF THE SECONDARY MAGNETIC FIELD APPROACH

COMPLETE CODE FOR GENERATING THE EXAMPLE 3 MESH: MESH_EXAMPLE_3.PY

COMPLETE CODE FOR CALCULATING THE EXAMPLE 3: RUN_EXAMPLE_3.PY

Data & Figures

About GSW

Resources

Explore

Connect

This Feature Is Available To Subscribers Only

$H$ -field approach — H

The $A - ϕ$ approach on mixed elements — $A_{m}$

The $A - ϕ$ approach on nodal elements — $A_{n}$