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

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 − ΦÞ 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 ðω → 0Þ and discusses direct or iterative solvers, boun-dary 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 edgebased 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élecelements, 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 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 1, the Helmholtz equation for the secondary electric field is derived (e.g., Grayver and Bürg, 2014) 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 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; or the primary fields for SF Nonnormalized complex-valued currents are currently not considered, but adding this feature is possible. The magnetic field H is related to E by where H is calculated on the same Nédélec-function-space as used for E by solving the FE problem

H-field approach -H
The Helmholtz equation in terms of magnetic fields H reads 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 11, the secondary magnetic field can be calculated by solving Equation 12 is derived in Appendix A. Similar to the E-field approach, the resulting linear system of equations reads and for the SF formulation where E can be calculated using Ampere's law in absence of external sources by Note that, alternatively, equation 18 can be expressed in terms of the primary electric field as derived in equation A-9 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 1 into Applying this method, Ansari and Farquharson (2014) derive the following FE formulation based on potential equations: 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) Ansari and Farquharson (2014) derive the resulting real-valued system of equations, which we modified to be symmetric where C and D are identical to E before (equations 14 and 15). Introducing additional test-functions η on a nodal space, the other submatrices are defined by Even though G is not symmetric, the complete system matrix (equation 26) for either total potentials ðA − ΦÞ or secondary potentials ðA s − Φ s Þ. The right side reads for the TF 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 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 ∇ · A ¼ 0, which is fulfilled if zero-Dirichlet boundary conditions are used, equations 24 and 25 can be modified to (Badea et al., 2001;Puzyrev et al., 2013) 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 35 and 36: As before, we derived the following real-valued system of equations: 2 6 6 6 6 6 6 6 6 6 6 4 The matrices P, Q, R and X, Y, Z are The right sides are for the TF formulation and for the SF formulation whereê fx;y;zg 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) 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 ¼ IdlδðxÞδðyÞδðzÞ; where I is the current, dl 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 (equa- The integral over the Dirac function equals one. Furthermore, the inner product of parallel or antiparallel vectors ðj e k AE ϕ i Þ equals AE1. Hence, every element b k of the right side vector b can be expressed, after assembly, as b k ¼ AEωIdl: evaluated with respect to the local directions of the Nédélec functions on corresponding edges. When using Nédélec spaces of higher (pth) polynomial order, the source contribution on each edge needs to be distributed on p dof, and therefore, equation 52 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 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 landbased, 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 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 (AE1) 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).

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 ΩmÞ with a ydirected HED source in the origin. The observation line extends AE5 km 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 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 airspace 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 3) and was located in the center of a 200 × 200 × 200 km mesh.
In Figure 5, the TF E solution at the surface based on secondorder 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 × 4 km 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 E x , E y , H x , H y , and H z (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 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 3, the black lines) with a quality of 1.2 for and a maximum facet size of 50m 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 × 200 × 200 km. 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 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 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 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 × 2000 × 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 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 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 8), 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 9. The match of the obtained SFs was already shown in Figure 8. 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 Table 2. Example 2: kd FEM − d ana ∕d ana k (%)considering all relative errors <100%for the real and imaginary parts of E and H, n specifies the amount of data points with a relative error <100% out of 40,401 values in total; meshes are introduced in Table 1   with the TF E approach and p2 for a half-space-like model with sinusoidal-topography and a 1 × 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 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.

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 Figure 8. Example 3: Secondary magnetic field H s at 10 Hz on 10 km x-directed observation line, calculated with the secondary E, H, A m , and A n approaches using p1 as well as E and H with p2, model: half-space, 1 km y-directed dipole transmitter, dipping plate, and either the isotropic or anisotropic brick anomaly ( Figure 6). Figure 9. Example 3: Total magnetic field H at 10 Hz on 10 km x-directed observation line, calculated with the total E and A m approaches using p1 as well as total E with p2: The results are crosschecked against SF H with p2, model: half-space, 1 km y-directed dipole transmitter, dipping plate, and either the isotropic or anisotropic brick anomaly ( Figure 6).

F28
area ð4 × 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 × 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 ≈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 cus-tEM repository.

DERIVATION OF THE SECONDARY MAGNETIC FIELD APPROACH
Faraday's law and Ampere's law, disregarding magnetic source terms, read ∇ × E ¼ −iωμH; (A-1) for the TF and