Efﬁcient goal-oriented mesh reﬁnement in 3-D ﬁnite-element modelling adapted for controlled source electromagnetic surveys

We developed a 3-D forward modelling code, which simulates controlled source electromagnetic problems in frequency domain using edge-based ﬁnite elements and a total electric ﬁeld approach. To evaluate electromagnetic data acquired across complex subsurface structures, software performing accurate 3-D modelling is required, especially for incorporation in inversion approaches. Our modelling code aims at ﬁnding a good compromise between the necessary solution accuracy at the points of interest and the general problem size by using a goal-oriented mesh reﬁnement strategy designed for models of variable electric conductivity and magnetic permeability. To formulate an improved error estimator suitable for controlled source electromagnetic problems, we developed literature approaches of mesh reﬁnement further targeting three aspects. First, to generate a roughly homogeneously ﬁne mesh discretization around all receiver sites, our new error estimator weights the adjoint source term by the approximate decay of the electric ﬁeld with increasing distance from the primal source using the expression for a homogeneous half-space. This causes almost no additional computational cost. Second, the error estimator employed in the reﬁnement approach can be optimized for models with pronounced conductivity and magnetic permeability contrasts as often encountered in, for example, mineral prospecting scenarios by optionally including terms that measure the continuity of the normal component of current ﬂow and the tangential component of the magnetic ﬁeld across interfaces of abutting elements. Third, to avoid amplitude-dependent over-reﬁning of the mesh, we formulate our element-wise error estimators relative to the local amplitude of the electromagnetic ﬁeld. In this work, we evaluate the implemented adaptive mesh reﬁnement approach and its solution accuracy comparing our solutions for simple 1-D models and a model with 3-D anomalies to semi-analytic 1-D solutions and a second-order ﬁnite-element code, respectively. Furthermore, a feasibility study for controlled-source electromagnetic measurements across ferrous mineral deposits is conducted. The numerical experiments demonstrate that our new reﬁnement procedure generates problem-speciﬁc ﬁnite-element meshes and yields accurate solutions for both simple synthetic models and realistic survey scenarios. Especially for the latter, characteristics of our code, such as the possibility of modelling extended sources as well as including arbitrary receiver distributions and detailed subsurface anomalies, are beneﬁcial.


3-D CSEM modelling with mesh refinement
1625 and complex distributions of material parameters in the subsurface. To evaluate CSEM data measured on the surface and in boreholes across complex-shaped mineral deposits, we developed an efficient and practical 3-D forward modelling algorithm optimized for mineral exploration purposes. However, the algorithm is applicable to all land-based CSEM set-ups, for example in volcanological, hydrothermal and groundwater studies.
The governing equations for CSEM modelling can be formulated in terms of total fields (Cai al. 2017;Li et al. 2017;Um et al. 2017), secondary fields (Schwarzbach et al. 2011;Grayver & Kolev 2015;Noh et al. 2016;Castillo-Reyes et al. 2018;Li et al. 2018;Soloveichik et al. 2018), coupled vector-scalar potentials using a total field formulation (Ansari & Farquharson 2014) or coupled vectorscalar potentials using a secondary field formulation (Puzyrev et al. 2013;Mukherjee & Everett 2018). In total-field FE modelling of CSEM data, the source, that generates the electromagnetic (EM) field and is typically located at the Earth's surface, is approximated by small current-carrying segments. In both total and secondary field approaches, the electric or magnetic field is modelled at cell nodes or along cell edges in the entire domain and is then interpolated to the locations of the receivers, where auxiliary quantities are calculated from the field values. Rochlitz et al. (2018) tested several of these approaches for 3-D FEM modelling and conclude, that total field approaches show the best performance with regard to computational cost and accuracy. Combining the total electric field approach with tetrahedral elements and vector edge interpolation functions, so-called Nédélec functions (Nédélec 1980;Jin 2014), ensures continuity of the tangential electric field components. Yet, it allows for discontinuity of some field components across element interfaces, which can be utilized for error estimation (Ren et al. 2013;Grayver & Kolev 2015).
Large 3-D domains and complex models result in a high number of elements, so that the number of unknowns to be solved for grows swiftly with increasing domain size and mesh density. Therefore, a balance between solution accuracy, resulting from sufficiently small model cells and a sufficiently large model domain, and problem size has to be found. This can be achieved by implementing numerically challenging higher order Nédélec elements (Schwarzbach et al. 2011;Grayver & Kolev 2015;Rochlitz et al. 2018) or by using adaptive mesh refinement, which is an acknowledged method in FEM (Rheinboldt 1981;Cendes & Shenton 1985;Larson & Bengzon 2013) and has as well been applied in geo-electromagnetics (e.g. Ovall 2006;Schwarzbach et al. 2011). Starting from a coarse mesh, global refinement approaches compute a global error estimator and iteratively refine all cells (e.g. by splitting each cell into two cells) until the global error estimator drops below a preconceived threshold that is assumed to be appropriate for the given problem.
In EM modelling, the receiver locations are the points of interest, where the field has to be determined with sufficient accuracy. In the rest of the computational domain, the fields only have to be resolved well enough not to reduce the solution accuracy at the receivers. Hence, the so-called goal-oriented refinement approach is interesting for EM modelling with FE approximations, because it aims at optimizing the solution accuracy by estimating errors in a defined measure of interest (Nitsche & Schatz 1974;Babuška & Miller 1984) such that global influences on local errors are considered (Ovall 2006). In goal-oriented refinement, element-wise errors can be estimated by expressing the error of the discrete problem solution in the quantity of a linear goal functional. This linear functional of the error can be related to the residual of the FE solution, which is itself a linear functional of the dual space to the function space of the primal boundary value problem (BVP). Thus, it becomes necessary to solve an adjoint or dual problem additionally to the primal BVP (e.g. Oden & Prudhomme 2001;Li & Key 2007;Ren et al. 2013). The solution of the adjoint problem is often called the "influence function" (Oden & Prudhomme 2001;Ren et al. 2013). The practical purpose of using the goal-oriented approach for error estimation in EM problems is to weight the element-wise error estimators with the influence functions to improve the solution accuracy at the receivers, which cannot be achieved by simply refining the mesh in the proximity of the receivers (Li & Key 2007). The source term of the adjoint problem can be defined such that adjoint sources are located at the points of interest (Ren et al. 2013;Grayver & Kolev 2015). During the mesh refinement procedure, only elements with error estimates above a given threshold are refined, with the locations of the adjoint sources having major influence on the spatial distribution of element-wise errors. Specifically for CSEM, there exist 2-D/2.5-D (Li & Key 2007;Key & Ovall 2011) and 3-D codes (Zhang 2018) using goal-oriented adaptive mesh refinement.
In the majority of published 3-D EM forward modelling algorithms, magnetic permeability μ is not reckoned a variable model parameter, presumably because contrasts in magnetic permeability are small compared with conductivity contrasts in most environments. However, ferrous ore bodies and mineralized zones often exhibit strong contrasts in magnetic permeability to their host rocks (Mooney & Bleifuss 1953;Zhdanov 2018). Neglecting magnetic permeability in inversions of data containing a strong signature of permeability can in fact lead to artefacts in the models (Farquharson et al. 2003), and it can thus be advantageous to include (Huang & Fraser 2001;Noh et al. 2016). Mukherjee & Everett (2018) state the importance of allowing for variable magnetic permeability when modelling for anthropogenic structures in the subsurface. Recently, computing EM responses of electrically conductive steel-cased wells became of interest (Puzyrev et al. 2016;Patzer et al. 2017;Couchman et al. 2018;Heagy & Oldenburg 2019), since they are often present or sometimes even utilized in EM measurements. They can also exhibit a significant contrast in magnetic permeability to the surrounding rocks. The magnetic properties of such artificial subsurface structures are easy to investigate in the laboratory and, in mineral exploration, magnetic surveys are often carried out before EM investigations, so that gross values of magnetic permeability are known prior to EM modelling. Consequently, keeping electric conductivity and magnetic permeability variable in EM forward modelling is beneficial.
Most publications associated with 3-D CSEM FE forward modelling focus mainly on verifying the accuracy of the forward responses as well as the time and memory efficiency of the software for specific and simple examples. This has the obvious advantage, that solutions from simple models can be compared to (semi-)analytic solutions or examples from the literature. However, especially when demonstrating land-based set-ups, the test models often have unrealistically high conductivities and well-shaped 3-D anomalies in combination with densely and regularly distributed receivers (Ansari & Farquharson 2014;Cai al. 2017;Mukherjee & Everett 2018) which is convenient for code validation but far from any real field set-up for onshore CSEM measurements. Rochlitz et al. (2018) present a realistic 3-D semi-airborne set-up for a conductive dike model. Applications to practical land-based measurement set-ups, realistic conductivity contrasts and anomaly shapes are often not published in articles dealing with forward modelling, even though this is of fundamental importance.
Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 Our forward modelling code is based on a total electric field approach and linear Nédélec functions on tetrahedral elements, which makes it feasible to accurately discretize topography, anomalies and extended sources. We keep both electric conductivity and magnetic permeability as variables to avoid possible errors when modelling set-ups with highly magnetic bodies. Dirichlet boundary conditions and a direct solver keep the set-up and solution of the system of equations simple. However, we aim at constraining the number of unknowns with an efficient goal-oriented automatic mesh refinement strategy. To this end, we set up an adjoint problem specifically for CSEM modelling. We furthermore adapt potent, published approaches of error estimation in magnetotelluric modelling (Ren et al. 2013;Grayver & Kolev 2015) to CSEM.
The overall aim of this paper is to report the development of a practical and efficient forward modelling algorithm for frequencydomain CSEM focusing on the design of a problem-specific goaloriented mesh refinement procedure and the application to a real field example. Section 2 describes the refinement approach for 3-D FE modelling including the formulation of the adjoint problem and our approach for calculating the error estimates. The numerical experiments in Section 3 validate our 3-D FEM CSEM modelling code and goal-oriented adaptive mesh refinement scheme by comparing the 3-D results of simple models to reference solutions. Additionally, we apply the forward modelling code to a realistic ore-body set-up with the aim of investigating the detectability of known deposits with CSEM. In Section 4, results are discussed and a general conclusion including an outlook towards the implementation of the refinement approach into the inversion framework emilia (Kalscheuer et al. 2008(Kalscheuer et al. , 2010 is provided.

3 -D F E M O D E L L I N G A N D R E F I N E M E N T A P P ROA C H
In this section, we describe the implemented FE formulation of the BVP for CSEM and introduce the methodology of our refinement approach consisting of (i) the establishment of an adjoint problem formulation for CSEM modelling based on the idea by Ren et al. (2013) which we complemented with a simple but efficient geometric weighting component crucial to CSEM set-ups, (ii) the definition of error estimators based on residuals and face jumps of normal current density and tangential magnetic field, similar to the approaches of Ren et al. (2013) and Grayver & Kolev (2015), but introducing weighting components to obtain error estimators relative to the local field amplitudes and (iii) the formulation of a general mesh refinement strategy based on goal-oriented refinement and an optional overall improvement of mesh quality during the refinement process.

The primal and adjoint problems
To describe the primal and adjoint CSEM problems, we use the weak form of the curl-curl equation of the total electric field in quasistatic formulation, first-order Nédélec vector FEs and homogeneous Dirichlet boundary conditions on a tetrahedral mesh to derive a system of linear equations suitable for solving with direct solvers. This procedure leads to the following discrete forms of the weak FE formulations of the primal problem and, by taking the symmetry of the left-hand side operator into account, the adjoint problem valid for each element (e) of volume V e and surface S e in the modelling domain using an e −iωt time dependence (e.g. Jin 2014;Cai al. 2017;Um et al. 2017). ω denotes the angular frequency, E e j the unknown coefficient for the total electric field and W e j the unknown coefficient for the adjoint solution (i.e. the influence function) along the jth edge of element e. N e i and N e j are the first-order Nédélec basis functions with the element edge counters i, j = 1, ..., 6. Considering the entire computational domain, the surface integrals cancel along the inner domain boundaries due to the continuity of the electric field. This approach accounts for spatially variable but element-wise constant isotropic electric conductivity (σ ) and magnetic permeability (μ).
The primal source current density J e s can describe galvanically and/or inductively coupled sources. J e r is the adjoint source, which we define as artificially injected vector influence sources at the receiver sites. In CSEM set-ups, the electric field component, that is approximately parallel to the orientation of the dipole source in the primal problem, decays with 1 r 3 in a homogeneous half-space (Zonge & Hughes 1991), where r denotes the distance between source midpoint and receiver location, and shows strong angular dependence. Therefore, we modify the approach of Ren et al. (2013) and Grayver & Kolev (2015) and implement artificial source currents J e r at each receiver location weighted by r 3 simulating a point source in each of the observed electric field component directions.
Details on the FE formulation including the application of the first vector Green's theorem to include boundary effects, the assembly of the global system of linear equations, the source and receiver implementations and comments on the adjoint problem formulation can be found in Appendix A. The r 3 scaling is further elaborated on in Section 2.4.

Error estimators
Identifying the cause of possible numerical errors in forward modelling is key to formulate error indicators, which guide the mesh Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 Figure 1. Two elements, labelled 1 and 2, sharing one face. The continuity of (a) the normal component of the current density J n and (b) the tangential component of the magnetic field H n× cannot be guaranteed across the common interface using an edge-based finite-element formulation for the electric field. Hence, face jumps in J n and H n× are used as optional terms in the error estimator guiding the mesh refinement. refinement. The continuity conditions of Maxwell's equations (Appendix A1) are fulfilled, when the differential equation for the primal CSEM problem is solved in its strong formulation (eq. A5). However, the discrete FE formulation using edge-based vector shape functions (eq. 1) only guarantees the continuity of the tangential electric field and the normal magnetic flux density across interior element interfaces (Ren et al. 2013). The normal component of the current density and the tangential component of the magnetic field are not ensured to be continuous across interior element interfaces in the discrete FE function space. Thus, so-called face jumps in normal current density and in tangential magnetic fields across element interfaces ( Fig. 1) cause errors, especially when the model contains contrasts in material parameters σ and μ. These face jumps and an additional residual arising in the FE approximation can be used to design an error estimator. The mathematical concept of error analysis is described in detail in, for example, Key & Ovall (2011) and Ren et al. (2013).
Our adaptive refinement approach is based on an error estimator η e L that is calculated for each element (e) by η e L = η e E η e W .
Here, η e E denotes the element error estimator derived from the primal problem. This part guides the refinement depending on the accuracy of the electric field solution in the entire modelling domain. The estimator η e W is the factor that adds the goal-orientedness to the procedure. It is obtained from the adjoint problem solution and is high for those elements that affect the solution accuracy at the receivers.
Both element error estimators are made up of three different contributions: the residual of the solution r e E (primal) or r e W (adjoint) in V e , the face jumps of the normal component of the current density j f e E (primal) or j f e W (adjoint) across element interfaces and the face jumps of the tangential magnetic fields h f e E (primal) or h f e W (adjoint) across element interfaces: Depending on the model, we can choose which (combination) of these three different components to include in the computation. For example, face jumps due to discontinuous tangential magnetic fields are supposedly small in case the model does not contain any magnetic permeability contrasts. Compared to other error estimation approaches using these three components (cf. e.g. Ren et al. 2013;Grayver & Kolev 2015), we do not use scaling factors for element residuals and face jumps consisting of maximum diameters of faces and tetrahedra. Instead, we introduce weights based on the local field amplitudes, which are defined in the following subsections and examined in Section 2.4.

Residuals
The primal residual r E corresponds to the difference between the right-hand side and the left-hand side of the primal boundary value problem (eq. A5) defined in the volume integral L 2 -norm (eq. A8): where the factor a E is a weighting factor , which we define as where c E is a positive non-zero constant that avoids division by zero. It can be chosen as any small non-zero positive number. The element residuals r e E are obtained by calculating the weighted sum of the residuals in the FE function space at four points within each element determined from Gaussian quadrature of second-order accuracy (Jin 2014, pp. 162-163). The element residuals of the adjoint problem r e W are obtained analogously.

Face jumps of normal current density
Using the constitutive relation J = σ E (cf. Appendix A1), the face jumps in current density for the primal ( j f e E ) and adjoint problems ( j f e W ) defined in the surface integral L 2 -norm (eq. A8) are determined by respectively, for a single face s l between two neighbouring elements. The indices 1 and 2 refer to the conductivities and the solutions of the two elements on either sides of the face.n denotes the normal unit vector on the face. The weights b E and b W are given by the Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 respectively, where the positive constants c E and c W are chosen the same way as they were for the residuals. The element face-jump error estimator consists of the sum of the face jumps to all four neighbouring elements. Face jumps are multiplied by 1 2 , as one face jump originates from two neighbouring elements:

Face jumps of tangential magnetic fields
Element-wise face jumps in the tangential magnetic field for the primal (h f e E ) and adjoint problem (h f e W ) are determined using the expressions respectively, for a single face s l between two neighbouring elements. The weights d E and d W are given by the expressions respectively, where the positive constants c E and c W are chosen the same way as they were for the residuals. h f e E and h f e W are obtained analogously to eq. (9).

The mesh refinement strategy
Elements with the most inaccurate solutions are identified by high element error estimators η e L , which exceed a certain threshold τ determined by The factor β remains constant during the refinement procedure and determines how large an element error estimator has to be, relative to the maximum element error estimator of the current mesh, for the corresponding element to be refined. The volumes of the elements chosen for refinement are bisected in the next meshing iteration. Once the element error estimators (eq. 3) are calculated, the average relative global error estimator η G can be obtained by where M is the total number of elements. Additionally, we define the average relative error estimator at the receivers η R as where N r is the number of receivers, and the index r denotes all elements containing receivers. An improvement in solution accuracy during the refinement can be expected if both, η G and η R , decrease during the refinement procedure. The iterative mesh refinement strategy starts with a coarse input mesh in iteration 1, where the elements are small and well-shaped only around the receiver and source locations. On the input mesh, primal and adjoint problems are solved and the element-wise error estimators are calculated (eq. 3). Based on these, the volume constraints for the new mesh are obtained and a new, finer mesh is generated, which serves as the input mesh for the next iteration. The procedure stops, when one of the below termination criteria is fulfilled. In this case, no new mesh is generated and the electric and magnetic field components at the receivers are calculated on the final mesh. Up to four termination criteria for the refinement loop can be selected: a preset maximum number of refinement steps, a maximum problem size, an unchanged problem size between successive iterations or a desired accuracy, that is if the average relative global error (eq. 13) is smaller than a predefined threshold.
In addition to the local goal-oriented refinement, a uniform improvement of mesh quality can be achieved during the refinement procedure. A low-quality mesh is characterized by elements with large ratios between the radii of the circumscribed spheres of the tetrahedra and the lengths of their shortest edges, which corresponds to a high value of the mesh quality factor q (Si 2015). Low-quality meshes ensure relatively small problem sizes and short computational times for the forward problems during the refinement procedure. However, to acquire accurate forward responses for the final mesh, globally well-shaped high-quality elements are required. Thus, one can choose whether one wants to apply timeand memory-consuming goal-oriented refinement on high-quality meshes, improve the mesh quality step-wise or improve the mesh quality in a last step during the refinement procedure. All three approaches add a global-refinement component to the goal-oriented one.
The improvement of accuracy (i.e. the decrease in relative deviation from a reference solution) with respect to the increase in degrees of freedom (dof) defines the convergence rate. During goaloriented refinement, it is highly dependent on the error estimation approach applied and the complexity of the model. In order to set the refinement procedure, that we propose, into context, an average convergence rate slightly better than −1 on a double logarithmic plot is common to the published goal-oriented refinement approaches and the number of refinement steps is usually not higher than 20 (Schwarzbach et al. 2011;Ren et al. 2013;Grayver & Kolev 2015;Zhang 2018). The convergence rates of global and non-goaloriented adaptive mesh refinement techniques have a theoretical quasi-optimal slope of − 1 Table 1. Qualitative investigation of the introduced weighting factors in our error estimator for all models tested in this work (cf. Section 3) using a decrease in average error estimators η G and η R , reasonable anomaly and receiver refinement as well as an accuracy improvement (decreased deviation from reference solution) as control indicators for a successful error estimator. Only r 3 -weighting of the adjoint source term (cf. Section A5) in combination with amplitude weighting of the error estimator components (eqs 6, 8 and 11) fulfils () all analysed control indicators.

Weighting
Decrease Decrease Receiver Anomaly Accuracy factors in η G in η R refinement refinement improvement r 3 and amplitudes r 3 Amplitudes None (minor) do not require such small elements, but refinement over a larger region. According to Key & Weiss (2006), a refined mesh can be used for two decades of frequencies. One could for example calculate the refined mesh for the central frequency of the two decades or for the highest frequency. Hereinafter, we refer to the signal frequency for which a refined mesh is generated as the refinement frequency.

Scaling/weighting factors in the error estimator
Our proposed error estimator is a combination of proven concepts (residuals and face jumps) and empirical scaling/weighting factors (r 3 -weighting of the adjoint source term and amplitude weighting of the error estimator components). The weights ensure, that the error estimator operates independently of the absolute field amplitudes in the computational domain as necessary for CSEM problems to avoid over-refining in regions of high amplitudes and under-refining in regions of low amplitudes (Key & Ovall 2011). Whereas the r 3 -weighting of the adjoint source achieves that all receiver locations are weighted equally in the adjoint problem with respect to their distance to the physical source, the purpose of the amplitude weighting is to obtain element error estimators relative to the local field amplitudes. We investigated, that we can only guarantee reasonable refinement behaviour, a drop in average error estimators and a decrease in deviation from a reference solution of the forward responses during the refinement of all models tested for this work (Section 3), when applying both introduced weighting factors (cf. Table 1). Applying only r 3 -weighting of the adjoint source results in unbalanced refinement of the central model region. In this case, the element error estimators η e E derived from the primal problem are large only exactly at the physical source and the elemental error estimators η e W derived from the adjoint problem are large only at the receivers, but they do not balance each other since their magnitudes differ by several orders. Maximum diameters of faces and tetrahedra as scaling factors in the error estimator (e.g. Ren et al. 2013;Grayver & Kolev 2015) lead, in our approach, to high element error estimators for the large elements at the domain boundary, which is not expedient for reasonable refinement.

Implementation
We developed the forward modelling algorithm in modern Fortran using automatic vectorization, shared-memory parallelization with OpenMP and linking to external software libraries. Thus, it could be conveniently implemented in our object-oriented inversion framework emilia (Kalscheuer et al. 2008(Kalscheuer et al. , 2010Kalscheuer 2020). For the mesh generation, we use the open-source software packages FacetModeller (Leliévre et al. 2018) and tetgen (Si 2015).  Table 2.
We apply the high-performance direct solver package PARDISO 6.1.0 (Schenk & Gärtner 2006;Coninck et al. 2016;Verbosio et al. 2017;Kourounis et al. 2018) to solve the system of equations. Furthermore, we adapted and included some sparse matrix routines from the SPARSEKIT package (Saad 1994) to work for complex matrices.

N U M E R I C A L E X P E R I M E N T S
To validate the refinement procedure and to evaluate the accuracy of the numerical solutions, we present the results of our algorithm for three different test models (Fig. 2): a homogeneous half-space (hhs) model, a layered (lay) model and a 3-D (bodies) model, which consists of one conductive and one resistive body embedded in a half-space. With the latter example, we want to investigate how the refinement procedure performs in the presence of local subsurface anomalies. We compare our solutions for the homogeneous halfspace and the layered model to semi-analytic reference solutions obtained with emilia (Kalscheuer et al. 2012(Kalscheuer et al. , 2015, which uses fast Hankel transforms (FHT, Christensen 1990). The reference solution used for the model with 3-D bodies was obtained with the custEM software (Rochlitz et al. 2018), which applies second-order Nédélec shape functions (p2 elements). The custEM mesh consists of ca. 2.5 million dof and the solution is thus accurate enough to serve as a reference for our code. All modelling domains have the size of 60 × 60 × 60 km 3 to avoid boundary effects. The magnetic permeability is set to vacuum permeability (4π × 10 −7 H m -1 ) everywhere. Source and receiver positions are listed in Table 2. The refinement threshold factor β is set to 0.85 for all models. Key & Weiss (2006) use a slightly different goal-oriented adaptive refinement approach, however, they suggest that β = 0.85 gives optimal results, which is confirmed by the results of our tests. It is well known that the mesh around the source and the receivers has to be Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 Table 2. Measurement set-up of the test models (Fig. 2). The dipole source discretized with 0.25 m long segments represents a long grounded cable. Receivers are distributed at a spacing of 25 m in line with the horizontal source.
x [m] y  Fig. 2 and three different combinations of error estimators, that is rJH, J and rJ. The terms r, J and H measure the residual in the weak form of the modelled partial differential equation, face jumps in the normal current density and face jumps in the tangential magnetic field, respectively. The 1-D semi-analytic reference solution is obtained using emilia. The meanings of the double logarithmic slopes − 1 3 and −1 are elaborated upon in Section 2.3. The relative deviation decreases most for the rJ approach.
sufficiently fine to achieve accurate solutions (e.g. Puzyrev et al. 2013;Soloveichik et al. 2018) and that reasonably adapted meshes result in higher accuracy than oversampled meshes (Castillo-Reyes et al. 2018). However, to ensure an accurate representation of the source even for low-quality meshes, the edges of the FE mesh coinciding with the source are just 0.25 m long and the receiver element sizes are constrained to have a maximum diameter of 1 m in the initial mesh. These settings correspond to the mesh characteristics of our 3-D reference solution.

Choice of optimal refinement strategy
The error estimator can be used with all components or combinations of the three components (eq. 4). To find the optimal approach, we refine the homogeneous half-space model using a medium quality factor of q = 1.4 for 9 refinement steps (10 mesh levels) and investigate the accuracy improvement for three cases: the threecomponent error estimator (rJH), the error estimator with only face jumps in current density (J) and with residuals and face jumps in current density (rJ). As an example, the results are shown for the real part of the E x component in Fig. 3. The reason why we predominantly show accuracy improvements for the real parts of the modelled fields is, that zero transitions in the imaginary parts can lead locally to low accuracy. All three approaches show convergence rates better than − 1 3 and −1 on a double logarithmic plot (cf. Section 2.3). At the 10th mesh level, the relative deviation from  Table 3) for all test models in Fig. 2. Reference solutions are the semi-analytic solution of emilia (hhs, lay) and the 3-D FEM solution of the custEM software (Rochlitz et al. 2018) for the bodies model. Markers are plotted for the initial mesh, the medium sized mesh after 45 refinement steps with q = 1.6 and the final high-quality mesh. During the refinement from the initial to the medium-sized mesh, the relative deviation decreases more for the models containing anomalies (lay, bodies). In all cases, the final refinement to a high-quality mesh leads to both an additional reduction in relative deviation and a significant increase in number of dof. the semi-analytic reference solution was reduced most for the rJ approach, even though the rJH approach behaves similarly. Using the 3-D (bodies) model as an example, we ascertained that face jumps in the tangential magnetic field are several decades smaller than face jumps in current density and increased neither at the resistive nor the conductive anomaly (Figs B1 and B2) for models that do not vary in magnetic permeability. Therefore, we use the rJ approach when verifying our refinement procedure on the test models (Fig. 2).
The concept of improving mesh quality in the last refinement step as proposed in Section 2.3 is investigated also for the homogeneous half-space model (Fig. B3). We found, that this procedure can result in a strong decrease of error estimators and in a small relative deviation from the semi-analytic solution. In terms of run times, this approach is favourable compared to refining on a high-quality mesh for this particular model.

Analysis of results for optimal refinement strategy
In this section, we investigate the convergence of the relative deviation from the reference solution at the refinement frequency (Fig. 4), the mesh density (Figs 5 and 6) and the element error estimator distribution (Fig. 7) during the refinement of our three validation models (Fig. 2). Modelling and refinement parameters are shown in Table 3.
We observe a decrease in the relative deviation from the reference solution from the initial mesh to the medium sized mesh after the refinement with q = 1.6 and further to the final high-quality mesh. The convergence is exemplarily shown for the E x component in  Fig. 2 at the Earth's surface around the receiver at x = 925 m (red dot). As a result of adjoint source and amplitude weighting in the error estimator, the elements around the receiver sites are homogeneously refined. Figure 6. Element volumes for the initial mesh (top), the medium sized mesh after refinement with q = 1.6 (middle) and the final mesh (bottom) of the layered test model lay in Fig. 2. The view plane is a vertical slice through the domain at y = 0 m. There is significant refinement around the layer interfaces at 500 and 1 000 m depth.
the low-quality mesh, which is more goal-oriented, than in the last refinement step, that includes the necessary improvement of global mesh quality.
In the homogeneous half-space model, the region around the receivers undergoes the strongest refinement (Fig. 5), as there are no other model regions that require special refinement. The transmitter is discretized extremely finely already in the initial mesh and Figure 7. Error estimators rJ of the initial mesh (left) and after refinement with q = 1.6 (right) for a vertical slice at y = 0 m around the conductive 3-D anomaly of the model bodies in Fig. 2. As the mesh is refined, the element error estimator values decrease around the receiver sites and around and in the block anomaly. no subsurface anomalies are present. In the layered model, successive refinement takes place additionally along the upper and lower bounds of the conductive layer between 500 and 1000 m depth (Fig. 6). How the refinement leads to a decrease in the element error indicators can best be observed in the model with 3-D bodies at the conductive block anomaly (Fig. 7). In the initial mesh, elements in and around the conductive block are large and element error estimates are high. The final mesh is refined especially in and around the anomaly. Element error estimates at the resistive block are smaller which is why this anomaly is less refined (cf. element face jumps in Figs B1 and B2).

Validation
For geometrical reasons, the electric component in x-direction (E x ) and the magnetic component in y-direction (H y ) are the two strongest field components in the measurement set-ups of our test models. Therefore, we check the E x and H y components along the receiver line against the reference solutions at 1, 10 and 100 Hz in Figs 8-10, respectively. Note that no receivers are located between x = −200 and 200 m in the vicinity of the source. The numerical results are in good agreement with the reference solutions. Deviations for all frequencies are predominantly less than 2 per cent for all three test models (Table 4, Fig. 4), except for the imaginary part of the E x component at 100 Hz. The latter results from the fact, that the imaginary part of the E x component has zero transitions at 100 Hz (Fig. 10) and is thus challenging to simulate.

Performance
Runtimes for the refinement procedure as well as memory consumption are reasonable. Compared to the examples and approaches by Schwarzbach et al. (2011) and Zhang (2018), our refinement approach appears to take longer (+18 and +90 per cent, respectively), but the memory consumption is in the same range for similar final Table 3. Modelling and refinement parameters for the test models in Fig. 2. The numbers of degrees of freedom (dof) are stated for the initial mesh and the refined (final) mesh. The refinement procedure terminated after 45 + 1 refinement steps, representing 45 iterations at a low mesh quality (q = 1.6) and one final iteration towards a high-quality mesh (q = 1.2 or q = 1.16).

Modelling parameters
Refinement termination criteria   Table 2. Grey symbols indicate relative deviations. Average deviations are listed in Table 4.      Fig. 12). The finite-element discretization was obtained using FacetModeller (Leliévre et al. 2018).
problem sizes (Table B1). Note, however, that a detailed quantitative comparison is not possible because of the differences between the example models, the start model design and the computer architectures. Since the number of refinement steps and the relative increase in dof during the refinement is intentionally a lot higher in our approach to ensure good goal-oriented behaviour during the mesh adaption and since relative deviations are low, the additional run time seems justified.

3-D test example with magnetic anomalies
To assess the behaviour of the error estimation approaches in a magnetically variable Earth, we investigate the bodies model (Fig. 2) in a scenario, where the two anomalies not only differ in resistivity but also in magnetic permeability (μ r = 5 , i.e. μ = 5 μ 0 ) from the surrounding half-space (μ = μ 0 ). We proceed similar to the test in Section 3.1, but use the error estimation approaches rJH, rJ, and JH. The relative deviation to the reference solution drops most for the JH approach (Fig. 11). All three tested error estimation approaches show a similar strong accuracy improvement with the first refinement step. By investigating the refined meshes, we see that the regions around the source and some receiver locations are refined by adding only few dof in the first refinement step. In the subsequent refinement steps, only the JH approach successfully refines the anomalies, so that the solution accuracy improves further.

Detectability study: ferrous mineral deposits
The study area, owned by Nordic Iron Ore AB, is located in the Bergslagen region in central Sweden, where metamorphosed volcano-sedimentary rocks are the dominant host rocks (Bastani et al. 2019). The deposits contain magnetite-and hematite-rich iron-oxide minerals. The general shapes, mineral contents and dip directions as well as magnetic susceptibilities of the mineral deposits are well known from borehole-logging, seismic measurements and susceptibility modelling (Maries et al. 2017;Bastani et al. 2019;Markovic et al. 2020). Deep reflection seismic investigations suggest deeper parts of mineralization and more mineralized horizons (Markovic et al. 2020) than previously known. Maries et al. (2017) derive a host rock resistivity of more than 1000 m from borehole resistivity logging, while the resistivities measured in the mineralized zones cover a wide range from 100 to 10 000 m. CSEM measurements planned for 2020 (Fig. 12) within the Smart Exploration project aim at determining more exact resistivities of the mineral deposits and investigating a possible extension to greater depth. The goals of our detectability study are to test the applicability of the forward modelling algorithm on a real field setting and to investigate the feasibility and practicability of the planned measurement set-up to detect the mineral deposits (Fig. 13).

Model characteristics
A susceptibility model by Bastani et al. (2019) suggests an average susceptibility of 0.3 SI for the mineral deposits ( Fig. 13) which is a relatively low value for magnetite-and hematite-rich deposits. More common are values of 1.5 up to 4 SI, for example, an iron ore formation in Minnesota containing 25 vol per cent magnetite (Mooney & Bleifuss 1953) and high-density skarn iron ores in northern Sweden ). Using the relationship between magnetic susceptibility κ and magnetic permeability μ (product of free-space permeability μ 0 [H/m] and relative permeability μ r ), we consider two different values of μ r for the mineral deposits in our model, that is 1.3 or 1.0, and use μ r = 1.0 within the host rock. Furthermore, we test three different scenarios for the resistivities of the mineral deposits: 10, 100 and 1 000 m. The planned measurement set-up (Fig. 12) includes two electric bipoles arranged in an L-shape and 37 distributed receivers measuring the horizontal electric and all three magnetic field components. The source transmits currents of 10 A and the source aerials are modelled as extended sources. In the initial mesh, the sources are discretized into 1 m long  14. Distribution of element volumes for the initial mesh (top) and the final refined mesh (bottom) for the model with magnetic mineral deposits (ρ = 1 000 m, μ r = 1.3; cf. Table 5). The view plane is a vertical slice through the domain at y = 2230 m. The view angle on the mineral deposits is the same as in Fig. 13. Nearby receiver locations are projected on the view plane and marked with red triangles. Enforcing volume constraints in the complex-shaped anomalies and applying refinement with a gradually increasing quality factor leads to reasonable refinement of anomalies and regions with receiver locations.
source-edge segments and each receiver is located in a tetrahedron of less than 1 m 3 volume.

Refinement
For investigating the refinement procedure, we use the long grounded cable source E12-C1 (Fig. 12), a homogeneous half-space model of 10 000 m and two models in which the mineral deposits have a resistivity of 1 000 m and a relative magnetic permeability of 1.0 or 1.3. For the homogeneous half-space model and the model with magnetic mineral deposits, the modelling and refinement parameters are listed in Table 5. The refinement frequency is set to 1 kHz, which is the logarithmic centre frequency of the frequency range planned for the measurements. Derived from a farfield skin-depth estimation (Zonge & Hughes 1991), a 1 kHz signal is supposed to penetrate to depth beneath most parts of the mineral deposits. As most of our receiver sites (Fig. 12) are however located in the transition zone, the effective penetration depth of a 1 kHz signal is expected to be less than 1000 m. The refinement of the homogeneous half-space model is performed similar to the non-magnetic test models in the previous section using the rJ error estimator. The convergence of the relative deviation from the reference solution computed using emilia and averaged over the receivers is presented in Table 6. Final relative deviations are below 3 per cent for the refined model with q = 1.2 and 883 914 dof. For comparison, we calculated forward responses on a high-quality mesh with q = 1.2 and 733 000 dof without applying refinement resulting in errors for the electric fields higher than 3 per cent. It is important to note that the mesh generator we use allows to set a target quality factor q, but not a target number of dof. Hence, we decided to base the comparison to an example without refinement on an identical quality factor.
The refinement strategy for the models containing mineral deposits has to be customized differently. To analyse the element error estimators, we compare the models with mineral deposits of 1000 m resistivity and relative magnetic permeabilities of 1.0 or 1.3. On the initial meshes, which are identical, the contrast of magnetic permeability of the model with magnetic deposits leads to element face jumps in the tangential magnetic field that are increased as compared to the model with non-magnetic mineral deposits. Given the rather small magnetic permeability contrast in the model with magnetic mineral deposits, this increase is small, but exceeds one decade in some elements at the surface of the mineral bodies (cf. Fig. B5). As illustrated in Section 3.5, the error estimator JH behaves reasonably for a model with resistivity and magnetic Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 Figure 15. Comparison of E y field components (amplitudes and phases) generated by the source E12-C1 at four different receiver locations, A3, A37, A7 and A36 (Fig. 12), for the homogeneous half-space model and three different mineral deposit scenarios (10, 100 and 1000 m and μ r = 1.3). The deviation of the FEM solution for the homogeneous half-space (hhs) model from the 1-D reference solution (obtained with emilia) is plotted as a yellow envelope.

Figure 16.
Comparison of E x field components (amplitudes and phases) generated by the source C1-E11 at four different receiver locations, A3, A37, A7 and A36 (Fig. 12), for the homogeneous half-space model and the three different mineral deposit scenarios (10, 100 and 1000 m and μ r = 1.3). The legend of Fig. 15 is valid here as well. permeability anomalies. Therefore, this error estimator is used for the model containing the magnetic mineral deposits. For the models with mineral deposits, refinement on an initial high-quality mesh (>1 × 10 6 dof) would be extremely time consuming on a desktop computer. Therefore, we use a coarse mesh of lower quality for initializing the refinement, but enforce volume constraints within the anomalies to avoid excessive refinement at the anomaly surfaces. We apply 10 refinement steps increasing the quality factor by 0.05 in each step. This procedure results in a favourable refinement behaviour: the regions around the source, the receivers and the anomaly are refined (Fig. 14), yielding approximately 1.1× 10 6 dof in the refined model with magnetic mineral deposits (cf. Table 5). Note, that refinement and the calculation of final forward responses is only performed for the model with magnetic mineral deposits.

Results
To evaluate the detectability of the mineral deposits with the planned measurement set-up (Fig. 12), we exemplarily consider the dominant field components of receivers A3 (between source and deposits), A37 (on top of the deposits), A7 (at the northern boundary of the deposits) and A36 (in the northeast of the measurement area, i.e., behind the deposits as seen from the source). Amplitudes and phases of the homogeneous half-space and the three different resistivity scenarios of magnetic mineral deposits are plotted (Figs 15-17 for the E y , E x and H z components, respectively). We also plot the absolute values of the deviations of the FEM solution for the half-space model from the semi-analytic solution (cf. Table 6) as envelopes around the amplitude values of the halfspace FEM solutions. If the responses of the anomaly models fall Figure 17. Comparison of H z field components (amplitudes and phases) generated by the source E12-C1 at four different receiver locations, A3, A37, A7 and A36 (Fig. 12), for the homogeneous half-space model and the three different mineral deposit scenarios (10, 100 and 1000 m and μ r = 1.3). The legend of Fig. 15 is valid here as well. Table 5. Modelling and refinement parameters for a homogeneous half-space model and the mineral deposits model with μ r = 1.3 for the deposits (Figs 12 and 13) with identical source-receiver set-ups. The number of dof is stated for the initial mesh and the refined mesh (final mesh). For the homogeneous half-space model, which we use as a reference for the detectability study, the refinement procedure terminated after 26 + 1 refinement steps, representing 26 iterations on a low-quality mesh (q = 1.6) and one final iteration towards a high-quality mesh (q = 1.2). For the model containing the mineral deposits, a series of 10 refinement steps increasing the mesh quality factor gradually from q = 1.6 up to q = 1.2 was identified to work best using the JH error estimator.

Model
Initial mesh Modelling frequencies 100 Hz to 10 kHz Refinement frequency 1 kHz Domain size 60 × 60 × 60 km q-factor 1.6 to 1.2 Table 6. Average deviations of E x and H y components from the 1-D reference solution (semi-analytic solution of emilia using fast Hankel transforms) in per cent for the homogeneous half-space model (cf . Table 5) at the refinement frequency of 1 kHz for three refinement levels (initial mesh, refined with q = 1.6, after the final refinement step with q = 1.2). Deviations of a discretization with q = 1.2 without refinement are given for comparison.

Initial mesh Refined mesh Final mesh
Mesh without refinement outside these envelopes, the measurements of the respective field components will potentially contribute to the detection of the mineral deposits. This is the case for all resistivity scenarios at receiver A37 and A7. Fields at receiver A3 and A36 are mainly influenced by the host rock and show only minor differences. For a thorough analysis on the detectability of the mineral deposits, measurements of the background noise need to be performed, and the deviations between the responses of the homogeneous half-space and the models including the mineral deposits need to be evaluated with respect to the noise measurements. As an example for potential borehole measurements, magnetic field responses are calculated at 400 m depth below receiver A37 (Fig. B6). Especially the H z field component at this location is influenced strongly by the mineral deposits.

Forward solutions
We have presented a forward modelling algorithm for 3-D CSEM problems with varying isotropic resistivity and magnetic permeability, which applies first-order edge-elements and a total electric field approach. It is capable of including grounded cable and loop transmitters of arbitrary geometry. For validation, we compared the forward responses of our FEM code to semi-analytic 1-D solutions and solutions from another FEM code. The algorithm produces accurate results for our three test models with a homogeneous receiver distribution and for a real measurement set-up using long grounded cables as sources. Most authors claim that higher order elements are beneficial in terms of accuracy and computational times. This is still to be investigated for our approach, as we primarily focused on reducing the problem size while increasing the accuracy using goal-oriented mesh refinement.

Refinement and performance
A goal-oriented refinement scheme was developed and tested. For the error estimation, we formulate an adjoint problem adapting the ideas from Key & Weiss (2006), Ren et al. (2013) and Grayver & Kolev (2015) to frequency domain CSEM. We introduced a new weighting of the adjoint source term based on the approximate decay of an electric field in a homogeneous half-space of 1 r 3 . As a result, our approach considers all receiver locations equally for the tested models without local anomalies. Another approach to define the adjoint source would include the electric field values at the receiver sites from the solution of the primal forward problem as weights. Whether this is a superior approach in our error estimator formulation is still to be investigated. After solving the forward problems, the primal and adjoint field amplitudes are known, which we utilize to obtain relative element-wise error estimates. Thus, model regions with anomalies and receivers are refined and the accuracy of the forward responses is sufficiently improved. Since we start with a fine mesh around source and receiver locations at the beginning of the refinement procedure, our initial deviations from the reference solutions are comparably low already. We think, this starting point leads to a more realistic expectation of the benefits that can be obtained by automatic mesh refinement than other published approaches. The strategy of increasing the mesh quality during the refinement process worked well for the presented models. We have to mention though, that refinement on a high-quality mesh is always preferable, if the computational power is available. Refining on a low-quality mesh and increasing the mesh quality only at the end of the refinement procedure requires more refinement steps and thus slightly more time than other published approaches. We found, that certain model properties such as pronounced subsurface anomalies require flexibility in the choice of parameters for the refinement procedure. Testing the error estimator with different combinations of its components revealed, that the specific model geometries and parameter contrasts determine which combination is most suitable for a particular model. For our mineral exploration example with comparatively low magnetic susceptibilities, the element face jumps in the tangential magnetic field at the magnetic anomaly contributed only with a mere fraction to the entire error estimator. However, as the error estimator that measures the strength of artefactual face jumps in normal current density and tangential magnetic field worked well for the block model example, we see significant potential in this error estimator when refining models with magnetically susceptible ore bodies or steel-cased wells. The introduced weighting approaches in our error estimator are more of a heuristic nature rather than being based on strict mathematics. This leads, for example, to a mixing of units when summing up the error estimator components. It also has to be investigated further, whether the constants introduced to avoid a division by zero in the error estimator components, need to be chosen more carefully. To this end, we implemented the refinement algorithm in a way, that it can easily be improved and customized to fulfil different modelling requirements.

Field example
Testing the forward modelling code on a realistic measurement setup with extended sources and irregular receiver spacing demonstrated the accuracy of the calculated responses also for such more complex measurement settings. We included the mineral deposits as detailed subsurface anomalies, as their shapes are well known up to a certain depth level. Precisely discretized subsurface anomalies are beneficial for forward modelling studies. Nevertheless, they can never be recovered in such detail by inversion schemes for data from diffusive EM methods collected on ground surface.
The detectability study for the ferrous mineral deposits leads to the conclusion, that the anomalies are detectable with the planned measurement set-up from the numerical point of view. Amplitude differences are, however, small in some field components, especially for the model with the lowest contrast in resistivity. For the amplitude of the E x component at receiver A37, the average relative difference to that of the half-space is only 17 per cent. Thus, to make meaningful statements about the detectability, the noise conditions in the measurement area have to be taken into account. We furthermore suggest, instead of having many receivers in the northeast of the measurement area, to place more receivers in the direct vicinity of the deposits and farther to the south than currently planned. Additional sporadic recordings of magnetic field components at depth would be beneficial and relatively easy to deploy, as several boreholes already exist in the measurement area. As our modelled results show major amplitude and phase differences only in the direct vicinity of the deposits, a dense coverage of this region will be needed to constrain the models inverted from field data well enough in the areas where anomalies are expected. A more advanced feasibility study including surface topography and considering different source locations and directions is in progress.

Towards inversion
The experiments demonstrated, that the presented refinement approach and strategy is efficient leading to accurate forward modelling. However, in an inversion, the related problem sizes can become much larger than in the presented examples. Some published CSEM inversion approaches apply refinement at every single inversion step, which can be time consuming and lead to artefacts in the models. One strategy worth to investigate is the following: The inversion is run on a coarse mesh without refinement, and the bestfitting model is then refined. Subsequently, this refined model serves as the start model for a second inversion again without refinement. This second inversion would be expected to converge relatively swiftly, as the initial model is already obtained by an inversion. The functionality of this concept will depend on the inversion approach and we intend to study it, when performing 3-D inversions Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 using the presented adaptive forward modelling algorithm within our inversion framework emilia.

A C K N O W L E D G E M E N T S
We thank Raphael Rochlitz for providing the custEM results for the bodies model. Nordic Iron Ore AB provided the mineral deposit model and Mehrdad Bastani the planned measurement set-up and the susceptibility model for the detectability study in Section 4.2. This work is funded by the Smart Exploration project. Smart Exploration has received funding from the European Union's Horizon 2020 Framework Programme under grant agreement no. 775971. The computations were partially enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPP-MAX partially funded by the Swedish Research Council through grant agreement no. 2018-05973. We highly appreciate the comments from Jochen Kamm and Rene-Edouard Plessix, which helped to substantially improve our manuscript.
Author contribution statement: The general forward modelling algorithm was implemented by PR with advice from TK. The implementation of the refinement approach was performed by PR and LMB. The experiments were set up and run by PR, who wrote as well most parts of the manuscript. LMB contributed in writing parts of Sections 1-3. The work and the paper was discussed with and proofread by TK.

DATA AVA I L A B I L I T Y
Synthetic data generated by the forward modelling algorithm can be provided upon request, except for data of the detectability study as they contain detailed geometric parameters of the deposits. We are not allowed to make the modelling code freely available due to the directives of our funding project. variable isotropic electric conductivity and magnetic permeability: where ∂ is the boundary of the modelling domain at which a homogeneous Dirichlet boundary condition holds.n is the outwardpointing normal vector on ∂ . Eq. (A5) describes a current source J s generating an electric field (E) in frequency domain inside , which decays with increasing distance to the source and is assumed to approximate zero at the domain surface ∂ . The boundary value problem is solved for the electric field with an edge-based total field formulation using the FEM.

A2 FE formulation
The edge-based vector FE formulation of eq. (A5) is described in detail in, for example, Jin (2014), Ren et al. (2013) and Key & Weiss (2006). We focus on the main points necessary to understand the relevant concepts. To retrace every step of the FE formulation including the calculation of the element matrices, we refer to chapters 5 and 8 in Jin (2014), where the 3-D vector FE analysis for tetrahedral elements is described.
To obtain the weak formulation of eq. (A5), we take the vector dot product with a vector test function v ∈ H 0 (curl, ) and integrate by parts using the first vector Green's identity (Jin 2014) for a scalar u and vectors a and b, which states: This procedure results in the following weak formulation of the curl-curl equation. Find E ∈ H 0 (curl, ) such that The Sobolev and solution spaces are defined as: H(curl, ) := {v ∈ [L 2 ( )] 3 , ∇ × v ∈ [L 2 ( )] 3 } and H 0 (curl, ) := {v ∈ H (curl, ),n × v = 0 on ∂ }, respectively. L 2 ( ) represents the Hilbert vector space of square-integrable functions with defined integral L 2 -norms for volumes and surfaces: Limiting the volume of integration in eq. (A7) to an element, the surface integrals of abutting elements are found to cancel another because the tangential electric field along the inner domain boundaries is continuous. Considering the entire computational domain as in eq. (A7), the surface integral over the outer surface ∂ is used to implement boundary conditions. To formulate the discrete equations, we use first-order Nédélec basis functions for tetrahedral elements (Nédélec 1980) as edge interpolation functions. For the jth edge with length l j of a given element e and using linear interpolation functions L jstart and L jend of the adjacent nodes, the first-order Nédélec basis functions are with where a, b, c and d are the coefficients for linear interpolation calculated from the nodal coordinates of element e. V e represents the element volume. L jend is obtained analogously. The point coordinates (x, y, z) represent any point in the modelling domain as, for example, a receiver location. The Nédélec basis functions are globally defined in a finitedimensional function space as piecewise-linear functions in 3-D for a tetrahedral element family (e.g. Nédélec 1980;Schwarzbach et al. 2011). They are vector edge functions and thus depend on the global definition of the direction of each edge. This global direction is usually determined by the mesh generator. Additionally, each global edge is assigned to the elements in the modelling domain containing this edge and numbered in a local element numbering scheme, requiring look-up tables to be generated that relate local and global edge and node indices. We use the local edge definition scheme as described in chapter 8.3.2 in Jin (2014). Each time, operations including the Nédélec functions on element level are executed, the global and local edge directions have to be included in the calculation (Anjam & Valdman 2015;Castillo-Reyes et al. 2018).
In the finite-dimensional function space using edge-based vector Nédélec functions and defining constant conductivity and magnetic permeability within each element, we achieve continuity of the electric and magnetic fields inside the elements and continuity of the electric fields tangential to the edges and faces. As a consequence of the continuity of the tangential electric field components and Faraday's law, the continuity of the normal component of the magnetic flux density across element interfaces is guaranteed. In contrast, the continuities of the normal component of current density and the tangential component of the magnetic field across element interfaces cannot be ensured (Ren et al. 2013). As the approximated electric field in one elementẼ e within the Nédélec function space is a vector field (Ẽ x ,Ẽ y ,Ẽ z ) T and defined as we obtain the explicit expression (eq. 1) for each element that has to be solved.

A3 Assembly of the global system of equations
We can define an element stiffness matrix as K e i j = V e 1 μ e (∇ × N e i ) · (∇ × N e j ) dV , an element mass matrix as M e i j = V e N e i · N e j dV and an element source term as b e i = iω V e N e i · J e s dV with i, j = 1, ..., 6, and, thus, an element system matrix is generated as The local systems of equations for all elements have to be assembled to a global one with the global system matrix A of size [dof × dof], which is complex symmetric, and the global source vector b: This system has to be solved for the array x containing the coefficients E e j of the Nédélec basis functions on the element edges (eq. A11).

A4 Source and receiver implementation
Electric dipole (ED) sources are implemented on the mesh edges similar to Rochlitz et al. (2018). If an edge with centre point (x 0 , y 0 , z 0 ) T contains an ED with a current density J s = I dl δ(x − x 0 )δ(y − y 0 )δ(z − z 0 ) (Ward & Hohmann 1988), where I is the source current, dl = (dx, dy, dz) T is the edge length and direction, and δ is the Dirac delta function indicating the source location, we insert this expression into the right-hand side of the weak formulation (eq. A7) to obtain The integral over the inner product N i · J s of the two parallel or antiparallel vectors N i and J s , which are defined along the same edge, is ±I|dl|. Assuming a source moment I|dl| = 1, each entry of the global right-hand side vector b is equal to ±iω for all source edges. The sign depends on the global direction of each source edge. For modelling extended electric dipole sources and loop sources, the source path can be placed along several connected edges ) of equal length. In this case, each source edge contributes I |dl| #source edges to the entire extended source. The final outputs of the forward modelling algorithm are the approximated electric and magnetic field components in Cartesian directions at the receiver sites. Electric field components are calculated using eq. (A11) and magnetic field components are obtained using eq. (A11) and Faraday's law (eq. A1). The receivers are located at points at the Earth surface in the modelling domain at different distances to the source. The initial mesh is designed such that the receiver locations are a few centimetres below the midpoint of the face of a small element coinciding with of the Earth's surface. Placing the receivers inside an element ensures continuity of all EM field quantities at the receiver sites. Ren et al. (2013) and Grayver & Kolev (2015) use a similar adjoint problem formulation for magnetotelluric modelling, where the primary electromagnetic field is generated by a distant source, such that the strength of the primary field is approximately the same at all receiver sites. In CSEM problems, however, the primary fields are generated, for example, by an electric dipole source at the Earth's surface. Thus, several orders of magnitude variation in field strength can be expected across a distributed array of sensors (Key & Ovall 2011). To compensate for the geometric decay of the fields in the adjoint problem and to achieve an equal influence of each adjoint source with respect to the primal source location, we implement adjoint source currents at each receiver location of amplitude r 3 , where r is the distance of the primal source to a given receiver. The adjoint sources are distributed over the edges of the corresponding receiver elements with the edge interpolation functions N e i (x, y, z) evaluated at the receiver coordinates.

A P P E N D I X B : C O M P L E M E N TA RY F I G U R E S A N D TA B L E S
Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 Table B1. Refinement and performance specifications of our modelling code for the example of the non-magnetic bodies model (Table 3 and Fig. 2) run on an hp 2240 (workstation) with an Intel R Core TM i7-6700 CPU @ 3.4 GHz quad-core processor and 64 GiB memory: number of refinement steps (ref. steps), computation time for the refinement (ref. time), number of dof of the refined mesh, relative increase in dof during the refinement, final average deviation from a reference solution of the E x component and peak memory consumption for the final forward calculation. For comparison, specifications for a published example of a layered seafloor model run with MARE3DEM by Zhang (2018) on an iMac with an Intel R Core TM i7 quad-core processor @ 3.5 GHz and a disc model by Schwarzbach et al. (2011) Figure B1. Element face jumps in (a) normal current density (jf E ) and (b) tangential magnetic field (hf E ) of the primal problem on the initial mesh (left) and on the final refined mesh (right) along a vertical slice at y = 0 m through the region around the 'resistive' but non-magnetic anomaly of the model with 3-D bodies in Fig. 2. Elements jf E are higher than hf E and the latter are not increased at the resistive anomaly. Both are reduced by more than two decades during the refinement, which only takes residuals and face jumps in current density into account. The grey line indicates the location of the anomaly.
(a) (b) Figure B2. Element face jumps in (a) normal current density (jf E ) and (b) tangential magnetic field (hf E ) of the primal problem on the initial mesh (left) and on the final refined mesh (right) along a vertical slice at y = 0 m through the region around the 'conductive' but non-magnetic anomaly of the model with 3-D bodies in Fig. 2. Elements jf E are highest in and around the conductive block. Element hf E are three decades smaller than element jf E and not increased at the anomaly. Both are reduced by more than two decades during the refinement, which only takes residuals and face jumps in current density into account. The grey line indicates the location of the anomaly. Figure B3. Behaviour of average global relative error estimator η G (top-left panel, cf. eq. 13) and average relative error estimator at the receivers η R (bottom left panel, cf. eq. 14) of the homogeneous half-space model in Fig. 2 using different refinement strategies and the error estimator rJ (cf. Fig. 3). η G and η R are normalized with the initial values of the mesh generated with q = 1.6. Tested are four different refinement scenarios (blue lines) and a high-quality mesh without refinement (orange diamond) is shown as a reference. The panels on the right associate the final relative deviations from the semi-analytic solution of the real part of the E x component with the computational run time (top right) and the memory consumption (bottom right) on an hp 2240 (workstation) with an Intel R Core TM i7-6700 CPU @ 3.4 GHz quad-core processor and 64 GiB memory for the different scenarios. Performing 45 refinement steps on a low-quality mesh with q = 1.6 and one last refinement step with q = 1.2 to generate a final high-quality mesh (blue +) is the most favourable refinement scenario for this model taking time, dof, memory, relative deviation from the semi-analytic solution and error estimates into account.  Figure B4. Behaviour of average error estimators η G and η R during the refinement procedure for the three test models (Fig. 2, Table 3). η G and η R are normalized to the initial values of the three models and drop significantly with increasing number of dof.
Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021 Figure B5. Initial distribution of element face jumps in tangential magnetic field of the primal problem hf E at the surfaces of the mineral deposits (1000 m) without a contrast in magnetic permeability (left, μ r = 1.0 for the mineral deposits) and with a contrast in magnetic permeability (middle, μ r = 1.3 for the mineral deposits) to the host rock (μ r = 1.0 in either case), and, for comparison, the initial distribution of element face jumps in normal current density of the primal problem jf E (right) for the non-magnetic mineral deposit. As the jf E distribution is the same for the magnetic and non-magnetic case, we only show the latter. The view angle on the mineral deposits is the same as in Fig. 13. Nearby receiver locations are marked with red triangles. Face jumps in tangential magnetic field increase across the surface of the anomaly, if the model contains contrasts in magnetic permeability. They thus seem in general useful to be included in the error estimation, although they are three decades smaller in magnitude than element jf E in this particular example. Figure B6. H x and H z field components (amplitudes and phases) generated by the source E12-C1 at receiver location A37 (Fig. 12) in a borehole at 400 m depth for the homogeneous half-space model and the three different mineral deposit scenarios. Large differences in amplitude and phase, especially in the H z component, suggest that sporadic recordings of the magnetic field components in some of the existing boreholes in the vicinity of the mineral deposits would be beneficial for determining more exact geometries of the deposits and their resistivities. Downloaded from https://academic.oup.com/gji/article/227/3/1624/6317613 by guest on 01 September 2021