## 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.

## 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\u2212\varphi )$ and enable the computation of SF on nodal elements by using iterative solvers. This approach, referred to as $An$, 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 $(\omega \u21920)$ 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 $Am$, 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 $Am$. 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$, $Am$, and $An$ 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

### $H$-field approach — *H*

### The $A\u2212\varphi $ approach on mixed elements — $Am$

### The $A\u2212\varphi $ approach on nodal elements — $An$

### Boundary conditions

### Incorporation of current density

This source-incorporation technique can be also applied to $Am$. 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 1. 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 2. 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 2a). Afterward, a third dimension ($z$) is added to the 2D mesh coordinates (Figure 2b), and the height for the surface is adjusted by either interpolating from a digital elevation model or using synthetic topography functions (Figure 2c). 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 2c). 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 2d). Domain markers are always automatically included. Only for illustrating the raw mesh (Figure 2b and 2c), TetGen divides every polygon into triangles. In the final mesh with using quality options for the TetGen call, all faces are well-shaped (Figure 2e). 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 2f). 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 3. 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 $je$ 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 ($\xb11$) 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 4. Electric and magnetic fields were calculated for 1, 10, 100, and 1000 Hz on a half-space $(100\u2009\u2009\Omega m)$ with a $y$-directed HED source in the origin. The observation line extends $\xb15\u2009\u2009km$ perpendicular to the HED. Results from TF $E$ (the colored lines indicating the frequencies) p1 computations were added in Figure 4. The underlying mesh for this FEM solution has a comparatively fewer amount of dof $(approximately\u2009\u200937\u2009\u2009k)$ to provide a simple example that requires only a few minutes computation time and $<8\u2009\u2009GB$ 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(Ey)$ for 1000 Hz and $R(Hz)$ for 100 and 1000 Hz. The differences occur at offsets >$1\u2009\u2009km$. Here, the $R(Hz)$ 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(Hx)$, which results from boundary effects. Further discrepancies are mainly noticeable at 100 and 1000 Hz for $I(Ey)$ and $R(Hz)$, 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 $103$, $102$, and $104\u2009\u2009\Omega m$. The air-space resistivity was set to $107\u2009\u2009\Omega 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 3) and was located in the center of a $200\xd7200\xd7200\u2009\u2009km$ mesh.

In Figure 5, 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\u2009\u2009m2$ within the $4\xd74\u2009\u2009km$ observation area (Figure 3, the red lines) and a TetGen quality of 1.4 in terms of the radius-edge ratio for the elements.

The single vector-components $Ex$, $Ey$, $Hx$, $Hy$, and $Hz$ (Figure 5a–5e) of the FEM solution match well with the analytic solution. The overall magnitude errors (Figure 5f) are smaller than 1%, disregarding the vicinity of Tx. The value $Ez$ 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 3, the black lines) with a quality of 1.2 for and a maximum facet size of $50m2$ 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\xd7200\xd7200\u2009\u2009km$.

Figure 6 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 $Ey$ and imaginary part of $Hz$ 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 6c and 6f) were significantly more accurate compared with all p1 computations. Furthermore, the p1 results on the fine mesh (Figure 6d and 6e) exhibited smaller randomly distributed errors but equal systematic misfits compared with the coarse mesh (Figure 6a and 6d). 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 6). 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 7. 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\u2009\u2009m3$). The p2 mesh uses a maximum cell size of $10,000\u2009\u2009m3$, resulting in approximately $64\u2009\u2009k$ nodes. The computational parameters are listed in Table 3. The complete domain size was $2000\xd72000\xd72000\u2009\u2009km$ (only approximately 3% more nodes compared with a 1000th of this size), needed to suppress boundary effects for the $H$ and $An$ 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 8 and 9 present magnetic field results for 10 Hz. Considering the SF computations in Figure 8, all approaches exhibited equal amplitudes for the dominant $Hz$ component and only slight mismatches for the slightly weaker component $Hx$. Minor discrepancies between the four approaches could be observed for $Hy$, 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 8), especially for the weak $Hy$ component.

We compared the TF $E$ and $Am$ approaches only with the (summed) TFs from the SF $H$-p2 approach (black line) in Figure 9. The match of the obtained SFs was already shown in Figure 8. As expected, the mismatch of the TFs is most prominent for $Hy$. 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 $Hy$ component. Another mismatch between the TF and SF computations can be observed at $x=\u22123\u2009\u2009km$ for $I(Hx)$. 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\xd71\u2009\u2009km$ 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 10.

The $E$ and $H$ fields (Figure 10c–10f) were interpolated on the surface and on a level 50 m above the surface, respectively (Figure 10). 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$, $Am$, and $An$, 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\xd74\u2009\u2009km)$ 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\xd74\u2009\u2009km$ 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 $\u224820%$. 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 $E0$ 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\u2212\varphi )$ 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.