Meshing strategies for 3d geo-electromagnetic modeling in the presence of metallic infrastructure

In 3D geo-electromagnetic modeling, an adequate discretisation of the modeling domain is crucial to obtain accurate forward responses and reliable inversion results while reducing the computational cost. This paper investigates the mesh design for subsurface models, including steel-cased wells, which is relevant for many exploration settings but still remains a numerically challenging task. Applying a goal-oriented mesh refinement technique and subsequent calculations with the high-order edge finite element method, simulations of 3D controlled-source electromagnetic models in the presence of metallic infrastructure are performed. Two test models are considered, each needing a distinct version of approximation methods to incorporate the conductive steel casings of the included wells. The influence of mesh quality, goal-oriented meshing, and high-order approximations on problem sizes, computational cost, and accuracy of electromagnetic responses is investigated. The main insights of our work are: (a) the applied numerical schemes can mitigate the computational burden of geo-electromagnetic modeling in the presence of steel artifacts; (b) investigating the processes driving the meshing of models with embedded metallic infrastructures can lead to adequate strategies to deal with the inversion of such electromagnetic data sets. Based on the modeling results and analyses conducted, general recommendations for modeling strategies are proposed when performing simulations for challenging steel infrastructure scenarios.


Introduction
The controlled-source electromagnetic (CSEM) technique aims at obtaining the electrical resistivity distribution of the subsurface from simultaneous measurements of time-varying electromagnetic (EM) fields generated in the Earth by an artificial electric or magnetic source.As a geophysical imaging method, CSEM is frequently applied to hydrocarbon prospecting [1][2][3][4][5][6][7][8], mineral mining exploration [9][10][11], CO 2 storage characterization [12][13][14][15], geothermal reservoir char-Second, high resistivity contrasts between metallic infrastructure and surrounding materials lead to ill-conditioned sparse systems.As a result, efficient numerical methods, meshing strategies, and solvers are required to incorporate such steel artifacts into the modeling routines.
Four major approaches exist for solving the discrete CSEM problem: finite differences (FD; [25,26]), finite volumes (FV; [27,28]), finite elements (FE; [29,30]), and integral equations (IE; [31,32]).These numerical methods are comparable in terms of implementation effort, accuracy, grid type (e.g., structured or unstructured meshes), and computational cost.A review of these numerical schemes for CSEM problems has been presented by [33].Among them, the FE method (FEM) is one of the most popular numerical approach to simulate 3D CSEM responses [1-4, 7, 8, 34-42].The popularity of FE schemes arises due their capacity to work with fully unstructured meshes, allowing to geometrically approximate complex subsurface bodies such as topography, bathymetry, small target structures or man-made infrastructures.In addition, FE schemes also provide mathematical robustness and solid convergence rates predicted by theory [43].
Recent advances in computer technology, storage capacities, and parallelization methods enable 3D forward modeling and inversion of EM fields.However, at the last stage of inverting experimental data, problem sizes are preferred to be reduced as much as possible to save time and memory as many forward problems have to be solved during an inversion process.In discrete modeling approaches, the solution accuracy directly depends on the mesh design for which the calculations are performed.More specifically, finer meshes usually give more accurate EM responses and result in larger problem sizes (i.e., more degrees of freedom (dof) in the forward problem).Retrieving the optimal mesh design, especially for complex models, can be a time-consuming task [44].
Therefore, approaches have been developed for designing problem-specific meshes, which aim at simulating accurate EM forward responses while reducing problem sizes as much as possible: High-order discretization schemes ( prefinement [8,34,35,45,46]) and goal-oriented meshing (h-refinement [47][48][49]) are widely discussed in the EM community but rarely implemented together in one forward solver.Presenting off-shore modeling examples, the development of a forward modeler for CSEM problems using a primary/secondary field approach and an adaptive mesh refinement strategy combined with higher-order polynomial FE on tetrahedral meshes was reported by [34].The authors favor the use of goal-oriented h-refinement and secondorder polynomials.However, they also state the need to investigate the proposed approach considering complex modeling parameter distributions.Using the FEM on hexahedral elements, a modeling code for general EM problems that sup-ports goal-oriented h-as well as p-refinement was presented by [35].A detailed analysis of magnetotelluric (naturalsource EM) examples shows that goal-oriented adaptive mesh refinement in combination with second-order polynomials gives the best result (e.g., high accuracy of the forward responses with respect to the number of dof).The authors claim that the presented results are also valid for 3D CSEM modeling.A recent study on high order polynomials using the spectral element method by [46] indicates, that p-refinement is particularly efficient in improving the numerical accuracy close to sources and receivers as well as in model areas without conductivity contrasts, whereas regions with conductivity contrasts require a more focused mesh refinement approach.
p-and h-refinement can be performed in the context of inversion.For example, the open-source 3D EM modeling and inversion package custEM [8] favours second-order basis functions to combine a globally good modeling accuracy with relatively coarse meshes [50].The open-source software MARE2DEM [51] generates new adaptively refined meshes with h-refinement for each particular resistivity model of each forward operation during the inversion process using a dual-mesh approach.
Focusing on forward modeling, the core motivation for this work is to perform an analysis that provides relevant information for a better understanding of high-order FE schemes in conjunction with goal-oriented meshing to solve 3D CSEM problems in the presence of steel infrastructure.Our contribution is twofold.First, robust numerical schemes for semi-automatic mesh design are introduced, which are required to simulate forward responses for 3D CSEM models with embedded metallic structures.Second, the analyses performed provide a solid foundation to justify and guide the future development of more robust meshing schemes, such as a fully automatic hp-refinement strategy.To investigate the numerical accuracy and computational effort of the proposed meshing schemes, two modeling setups of varying complexity and relevance for academia and industry are considered.For computing synthetic EM responses, two 3D CSEM modeling routines, namely elfe3D [11] and PETGEM [7], have been utilized.Both codes are based on the edge finite element method (EFEM) developed by [29].
The paper is organized as follows: Section 2 introduces the theoretical background on the governing equations of 3D CSEM forward modeling, the 3D CSEM algorithms under consideration and the steps of our h-and p-refinement approaches.Furthermore, strategies for the numerical representation of metallic-cased wells are described.In Section 3, an analysis of synthetic EM responses for two distinct CSEM setups is conducted.For each test case, the focus is placed on the examination of their physical characteristics and the challenges they present in terms of mesh design and numerical accuracy.In Section 4, the essential aspects that govern the numerical accuracy and computational effort of the mesh-ing schemes used to solve 3D CSEM models in the presence of metallic artifacts are discussed.Based on the numerical results and analyses, meshing and modeling recommendations are formulated for addressing 3D CSEM problems involving metallic-cased wells.In Section 5, the conclusions of our study are presented.

Problem statement
The governing equations describing the 3D CSEM forward modeling problem are the frequency-domain Maxwell equations in a diffusive form, hence neglecting displacement currents.An electric dipole oscillating at a single frequency is considered in this study.By assuming a time-harmonic dependence expressed by e −iωt , the first two Maxwell equations can be stated as where E is the electric field, H is the magnetic field, i is the imaginary unit, ω is the angular frequency, and μ 0 is the free-space magnetic permeability.J s describes the source electric current, σ is the electric conductivity, which can also be expressed in terms of electric resistivity ρ = 1/σ , and σ E is the conduction current.The substitution of Eq. 1 into Eq. 2 yields to also known as the curl-curl equation in terms of the total electric field [52].This formulation can prevent numerical errors that arise when the source is located within the region of anomalous properties (e.g., models with high conductivity contrasts or with bathymetry/topography variations) [38,45].

Modeling routines
This study employs two 3D forward modeling routines: elfe3D (for goal-oriented h-refinement) and PETGEM (for p-refinement).The features of each algorithm may be summarized as follows: 1. elfe3D (modeling with the total electric field approach using finite elements in 3D) is capable of calculating forward responses for frequency-domain 3D CSEM setups using a direct forward solver PARDISO [53][54][55] or MUMPS [56].It is designed in modern Fortran using shared-memory parallelization with OpenMP.The governing equations are discretized using linear Nédélec interpolation functions.Variable model parameters are isotropic electric resistivities, magnetic permeabilities and, in case displacement currents cannot be neglected, electric permittivities.Line or loop sources are modelled along element edges in their correct extensions.elfe3D was validated in [11].The code provides the functionality of goal-oriented mesh refinement (h-refinement), which can be customized for specific model characteristics [11].Overall, it is based on error estimators obtained by calculating face jumps in the normal current density and the tangential magnetic field as well as residuals arising in the finite element formulation [35,57] combined with amplitude-dependent weights.Depending on the model, the user can choose a combination of these three error estimator components that is most expedient for the given model.Additionally, a global mesh quality improvement (q-refinement) can be applied during the refinement procedure.The factor q specifies the mesh quality, defined for each tetrahedron in a mesh by the ratio between the radius of its circumscribed sphere and the length of its shortest edge.A high q-value corresponds to a low-quality mesh and viceversa [58].

PETGEM (Parallel Edge-based Tool for Geophysical
Electromagnetic Modelling) is a parallel and high-order routine for active-source and passive-source geophysical EM modeling on realistic setups.To discretize the governing equations, the code uses a so-called mixedorder scheme proposed by Nédélec [59].The capacity and performance of the software have been assessed in shallow marine hydrocarbon exploration [7,44,45], geothermal reservoir exploration [17], 3D marine CSEM modeling in the presence of vertical transverse isotropy and complex bathymetry [60], and magnetotelluric modeling for test cases with large-resistivity contrasts [61].

Refinement schemes
The workflow followed to conduct 3D CSEM meshing and modeling experiments is depicted in Fig. 1 and consists of four overarching stages: 1. Subsurface model definition, composed of closed volumetric regions with a constant resistivity for each.The source parameters, measurement setup, and metallic artifact locations are defined at this stage.2. Input mesh generation (Initial mesh), accomplished using calls to TetGen routines [58].The resulting mesh satisfies minimum quality criteria in regions close to the source and receiver locations.3. Tailored mesh generation (refined mesh), the goaloriented refinement process of the initial mesh is performed using calls to routines from elfe3D algo-rithm [11].Error estimates based on element residuals are utilized in our approach.A refined and higher-quality mesh is obtained at this stage.4. Computation of synthetic EM responses, accomplished using the PETGEM code [7].Simulations for both input meshes (Initial and Refined mesh) using different highorder FE basis are carried out. 5. Performance assessment, which consists of analyzing the numerical accuracy, computational effort (run-time and memory needs), and practicability.Additionally, forward responses for medium-quality and high-quality initial meshes are calculated in our study.Different high-order FE basis functions are employed to assess the potential benefits of goal-oriented refinement for the presented models.The objective is to determine whether high-quality meshes and high-order elements yield comparable accuracy and computational effort, or if the goal-oriented refinement approach provides advantages.

Numerical representation of metallic-cased wells
Modeling EM responses in the presence of metallic artifacts has been a longstanding problem in EM research.The main challenge is to discretize the transition from poorly or moderately conductive cells with relatively large volumes of geological structures to highly conductive elements of small volumes representative of the steel casing.This issue results inescapably in extensive mesh refinement and, therefore, increased computational expense.Several approximations have been proposed to face these difficulties, each having its strengths and weaknesses.Table 1 describes the most relevant approaches published in the literature.Those numerical approximations fall into two categories: (a) algorithms that preserve the exact geometry and full physics (high-computational cost), and (b) modeling routines that prioritize efficiency by relaxing accuracy in some regions of the domain (low-computational cost).Among the numerical approximations presented in Table 1, the choice for modeling metallic casings is to use linewise perfect electric conductors (PEC).The metallic-cased well is represented as a linewise infinitely conductive artifact in PEC approximations.This volumeless approach positively impacts the computational effort since excessive mesh refining is avoided.In this study, the implementation of the PEC approximation is carried out as follows: 1. Target edges: determine the edges of the mesh that coincide with the metallic-cased well axis 2. Target dofs: determine the dofs at the target edges 3. Zero-electric-field condition: force the electric field to zero at the target dofs (PEC dofs).Here, the diagonal entries of the left-hand side of Eq. 3 that correspond to the PEC dofs are set to 1, and the non-diagonal entries are set to 0. Consequently, PEC dofs of the right-hand side are set to 0. 4. Solve the resulting modified version of Eq. 3 Inspired by [63], the PEC approach has been used before for approximating steel-cased wells by [64], where several combinations of solid rectangular prisms and PEC have been assessed for CSEM models with deep metallic-cased wells.One conclusion of this work is that PEC-only approximations are not necessarily accurate enough, in case a source is located close to the well.Therefore, the part of the borehole closest to the source is modeled with thin, conductive prisms to better account for induced current on the casing surface.How large this prism part of the casing compared to the PEC part has to be is not only dependent on the distance to the source but also on the source frequency and conductivity structure of the background geology [64].However, a steel-casing sufficiently distant from a source can be effectively modeled with a PEC-only approximation.In this paper, the assessment of the effectiveness of PEC approximations is expanded to include combinations with high-order edge FE methods and goal-oriented meshing strategies.

Simulations
The validation of the proposed meshing methodology is conducted by computing the electromagnetic (EM) responses for two different models.This is achieved by comparing meshes at different refinement levels in the respective scenarios.Each model presents a particular modeling challenge, being suitable to investigate the mesh refinement strategy and the parameters of the modeling routines (e.g., mesh refinement grade, impact of polynomial order, computational effort to solve problems that contain multi-scale meshes).Model 1 includes a cross-line source-receiver setup and metalliccased wells.Model 2 consists of an in-line source-receiver setup over basin topography and a metallic-cased well.
For all simulations, the normalized root-mean-square difference (NRMSD) between two electric amplitude responses, denoted as Q 1 and Q 2 , is defined as follows to which we refer as the normalized difference throughout the article.The average NRMSD is then defined as the average of this quantity for all responses of a given model.According to the workflow presented in Fig. 1, the mesh generation process is carried out using TetGen [58].It is well known, that the model discretisation at and in the direct vicinity of the receivers strongly influences the accuracy of the forward responses [71].Also, the discretisation of extended sources has to be performed with caution, especially if the source is located close to resistivity contrasts in the model [72].Therefore, for each test case, an initial input mesh is constructed manually to contain small elements only at the receiver sites and at the extended source lines.Our tests have shown that this is an important prerequisite for rel- atively accurate forward solutions in combination with small problem sizes from the start.The initial input mesh is then refined using the goal-oriented mesh refinement described in Section 2.2 in combination with mesh quality improvement (q-refinement, see Section 2.2), which is imposed by decreasing the mesh quality factor q. The linear systems of equations arising from the 3D CSEM setups are solved using the multifrontal direct solver MUMPS.All PETGEM simulations have been performed on Marenostrum 1 using 384 CPUs.

Model 1: modeling in the presence of metallic-cased wells
It is considered numerically challenging to simulate EM responses in the presence of a steel-cased well, because the casing has a long, but thin hollow structure and it is million times more conductive than the surrounding geology.This extreme multi-scale nature makes conventional EM algorithms unsuitable to model casing EM problems.The test model under consideration is depicted in Fig. 2a.It consists of two solid vertical steel-cased wells embedded in a homogeneous half-space (8 m).The outer diameter of the well casing is 0.2 m, the thickness is 0.0122 m, the electrical resistivity is 2 −7 m and the relative magnetic permeability is 1.The two wells are 500 m deep and 195 m apart from each other.A 5 Hz y-directed 500 m long horizontal electric dipole source with a source moment of 500 Am is located at x = 150 m.The surface electric dipole source is not directly connected to the nearest wellhead but 10 m away from it.This configuration can provide reasonably strong EM coupling between the source and the wells.120 receivers are distributed cross-line along the x-axis between x = −900 m 1 Marenostrum supercomputer at BSC. Motivated by [64], borehole casings are modelled with a combination of long solid rectangular prisms (0 − 100 m depth) with a diameter of 0.14 m and PECs (100 − 500 m depth).The entire length of the casings cannot be modelled with a PEC-only representation in this model due to the close source.The solutions obtained using PETGEM are compared against the EM approximations presented in [64].This reference solution was generated by a 3D FE code, and has been extensively verified through comparison with FD solutions [22,73,74].

Model 1: mesh design
An initial mesh (referred to as Initial Mesh I, with a quality factor of q = 1.6) was designed, and a refined mesh (referred to as Refined Mesh) was generated using goaloriented refinement, gradually increasing the mesh quality up to a quality factor of q = 1.4.These meshes were used to perform the simulations.Our goal for both meshes was to generate fewer elements than the mesh of the reference solution contains, which consists of 1 544 931 elements.
It was observed that applying automatic refinement to this model using the goal-oriented approach of elfe3D proved to be challenging.First, high-error estimates are obtained at the upper parts of the wells, primarily attributed to the significant conductivity contrast between the prisms and the half-space.Second, the mixture of tiny elements around the wells and large elements in the rest of the half-space (see Fig. 2b) is problematic.The refinement of small elements in the well regions required the addition of more elements, increasing the computational cost of the numerical simula- For PETGEM, the polynomial basis order ( p), number of dof, modeling run-time (sec), memory (GB), and average NRMSD x (cross-line, %) and NRMSD y (in-line, %) are given.It is worth noting that the authors of the reference solution [64] did not provide information regarding the modeling run-time and the memory requirements tions.The generated meshes consist of 134 989 (Initial mesh I) and 401 614 (Refined mesh) elements.Test simulations studying the influence of the source discretisation on the accuracy of the forward responses are particularly important for the presented model, as the source is only 10 m away from the receiver line and the steel infrastructure.It was observed that a source discretization with at least 12.5 m for p = 1 and 25 m long segments for p = 2, 3 led to the desired result and that PEC approximations require similar segment lengths as the source.Additionally, forward responses for mediumquality (q = 1.4) and high-quality (q = 1.2) initial meshes were calculated (Initial mesh II and III, see Table 2).
Fig. 3 Cross-line forward responses for Model 1 (Fig. 2) obtained with the Refined mesh; NRMSD x for PETGEM ( p = 1, 2, 3) in the central region of the profile above the steel-cased well are included

Model 1: cross-line analysis
To assess the accuracy of the forward responses, the amplitude of the x-component of electric field (E x ) along the complete profile as well as the resulting NRMSD x on the Initial mesh I and the Refined mesh in the region of the well at x = −45 m are depicted in Fig. 3. Performance measures and average NRMSD x of the E x field component are presented in Table 2.It is worth mentioning that the influence of the metallic well on the E y component is smaller than that on the E x component in this case.Therefore, including a figure showing the E y responses does not add relevant information.
To preserve brevity, we will focus on analyzing the cross-line responses in detail.However, the numerical results and conclusions derived below remain valid for in-line analysis (see Fig. 7).
Comparing the forward responses presented in Fig. 3, it becomes clear that it is most difficult to achieve good accuracy at the locations of the wells and near the source, since the amplitudes vary the most there.The normalized differences in the mentioned regions show that the mesh refinement does not lead to a good local improvement of the forward responses.Globally, a smaller NRMSD x can be achieved with the Refined mesh, while the improvement is more significant for elements in p = 1 (see Table 2).
On top of that, the main observation for Model 1 is that high-order elements ( p = 2, 3) produce the best and most stable results.Even on the Initial mesh I, the PETGEM solution with p = 2 elements has an average NRMSD x as low as 1.11 using only 854 898 dof.Simulations on the higher quality initial meshes (Initial mesh II and III) and on the Refined mesh do not result in better accuracy (NRMSD x of 1.35, 1.87 and 1.10, respectively), while the computational cost is increased (see Table 2).Accordingly, performing p = 2 computations on the Initial mesh I is the best option to generate an accurate forward solution for the presented model, both in terms of computational effort and accuracy.

Model 2: modeling of the Vallès basin
As a second example, the Vallès basin model, proposed by [17] and depicted in Fig. 4a, is considered.This 3D CSEM model exhibits basin floor topography which is a key difference compared to the previous example.The main setups and results of the Vallès basin model including real data evaluation are reported and discussed in [17].For the purposes of this study, the most realistic configuration presented by the authors is used as the base setup.Thus, the test case model consists of a 4 km thick air layer (10 8 m), resistive basement (1 000 m) with topography, and conductive sediments (20 m).A 2 Hz horizontal electric dipole is located at x = 1209 m, y = 2000 m, and z = −2 m.The transmitter moment is 1 Am.The metal well casing at profile distance 1.719 km can be modelled with a 200 m long PEC (0 − 200 m depth), as the source is too far away from the borehole to generate current coupling between the source and steelcasing.
The EM responses are compared against p = 3 approximations computed with PETGEM on a globally fine mesh, which serves as the reference solution.For the reference model, a vertical cylinder (0.001 m) to model the metallic well casing is used.The cylinder is centered at x = 1719 m, y = 2000 m, and z = −100 m, its length is 200 m and it has a radius of 0.2 m.The accuracy of PETGEM solutions, when using a cylinder representation for the well, has been previously verified for this model in [43,45].Again, the main goal is to generate meshes with fewer elements than the mesh of the reference solution, which consists of 946 632 elements and 17 674 324 dof, but to retain a comparable solution accuracy.

Model 2: mesh design
As in the previous cases, the initial mesh is designed with fine elements only along the receiver profile in x-direction and the source at x = 1 209 m.The metal well at x = 1 719 m is discretized using a PEC only (0 − 200 m depth with line segments of 25 m).It was observed that modeling the shallow well with a replacement prism or a cylinder is not necessary to achieve accurate solutions for this test case, presumably because the source is too far away from the borehole to cause strong coupling.The presence of basin floor topography produces a finer inner model region in the initial mesh compared to the previous test cases, but the PEC representation of the well needs substantially less elements than modeling a prism which requires small elements.The Initial mesh I is obtained with a global mesh quality factor q = 1.6 and consists of 235 676 elements.It is refined with goal-oriented refinement and an improvement in mesh quality towards q = 1.4 as in Model 1 (see Section 3.1) resulting in a Refined mesh of 618 526 elements.Its central region is depicted in Fig. 4b.Again, forward solutions on higher quality (q = 1.4 and q = 1.2) initial meshes were also obtained (see problem sizes of Initial mesh II and III in Table 3).

Model 2: in-line analysis
To present our results for Model 2, the amplitude of the in-line component of electric field (E x ) as well as the resulting NRMSD x on the Initial mesh I and Refined mesh in the region of the well are plotted in Fig. 5.The modeling statistics are depicted in Table 3.The results of Model 2 show that goal-oriented mesh refinement in conjunction with Fig. 5 In-line forward responses for Model 2 (see Fig. 4) calculated on the Refined mesh.NRMSD x for PETGEM ( p = 1, 2, 3) in the central region of the profile above the steel-cased well at x = 1 719 m are included For elfe3D, the refinement run-time (sec) is provided.For PETGEM, the polynomial basis order ( p), number of dof, modeling run-time (sec), memory (GB), and average NRMSD x (in-line, %) and NRMSD y (cross-line, %) are given p = 2 discretizations is a good option to obtain accurate EM responses in combination with acceptable problem sizes for complex, but appropriately discretized, models.More specifically, on the Initial mesh I, the achieved accuracy was at a non-acceptable level for all element-orders (average NRMSD x above 7.67%).The applied h-and q-refinement turned out to work best to lower the average NRMSD x values (average NRMSD x for p = 1 elements of the Refined mesh 3.92%), which could not be achieved by improving the global mesh quality only (average NRMSD x of the Initial meshes II and III for p = 1 elements are 8.24% and 4.43%, respectively).Furthermore, the accuracy is substantially improved by using p = 2 elements on the Refined mesh (average NRMSD x 0.83%).Presumably due to the PEConly representation of the well, the goal-oriented refinement resulted in nicely improved responses near the well for this Fig. 6 Cross-line forward responses for Model 2 (see Fig. 4) calculated on the Refined mesh.NRMSD y for PETGEM ( p = 1, 2, 3) in the central region of the profile above the steel-cased well at x = 1 719 m are included test case (inlays in Fig. 5), without adding extensively many elements (see Fig. 4c).The problem size (3 918 118 dof) of our preferred modeling option (e.g., Refined mesh in combination with p = 2 elements) is also acceptable for a model with topography and small scale features and is much smaller than the problem sizes of the reference solution (17 674 342 dof) and the original solution (5 344 746 dof) for Vallès basin model presented in [17].Comparable accuracy without goaloriented refinement to our preferred modeling option can also be reached on the high-quality initial mesh (Initial mesh III) with p = 3 calculations (average NRMSD x 0.98%).However, memory consumption and problem sizes for Initial mesh III and p = 3 calculations are more than three times larger than for the Refined mesh.Additionally, taking the run-time for the refinement into account, run-times are also longer for the Initial Mesh II (2784.27sec on Initial mesh III with p = 3 elements and 2114.55 sec on the Refined mesh with p = 2 elements).

Model 2: cross-line analysis
The amplitude of the cross-line component of electric field (E y ) and the resulting NRMSD y on the Initial mesh I and Refined mesh in the region of the well as well as the modeling statistics are summarized in Fig. 6 and Table 3, respectively.Be aware that the E y component is the weak component in this modeling setup (e.g., E y amplitude is overall more than two decades smaller than the E x amplitude).It is thus expected to be more difficult to obtain good accuracy for the E y component.The average start NRMSD y is in the same range as the average start NRMSD x (NRMSD x = 10.95 and NRMSD y = 10.36 on Initial mesh I for p = 1 approximations).The overall observations, i.e. that the accuracy of the responses can be improved most significantly with goal-oriented refinement in combination with p = 2 calculations, holds in the same way as for the E x component.However, the average NRMSD y for this meshing option is Fig.

Discussion
Accuracy and performance metrics were investigated using the presented 3D CSEM setups.In Model 1, EM responses are simulated in the presence of two steel-cased wells.This model is considered numerically challenging due to the large scale variations and high-conductivity contrasts.Model 2 also exhibits large scale variations due to a metallic casing and additional basin topography.Figure 7 shows a performance metrics summary (condensed results of Tables 2  and 3).Here, the average NRMSD versus run-time and memory needs of the simulations is presented for each test model, each mesh, and each basis order p = 1, 2, 3.

Accuracy and performance evaluation
The numerical results indicate that a good match between the reference solutions and our synthetic EM responses can be achieved for both of the presented models.However, it was observed that it is complicated to achieve low NRMSD values near wells, especially for the weaker component of the electric field.From the analysis of the NRMSD values, it can be inferred that the goal-oriented strategy does not perform satisfactorily when there is a mixture of elements with significantly different sizes and high-resistivity contrasts (Model 1).In this case, it can be seen that highorder approximations ( p = 2) can offer a good balance between numerical accuracy and computational cost.For a more balanced mesh in terms of element sizes and resistivity distribution (Model 2), goal-oriented refinement in combination with p = 2 approximations results in the best trade-off between accuracy and computational cost.
In general, code performance depends on the input model, mesh quality, solver-type, and computational architecture.Run-time for modeling and memory consumption scale with problem size.Also, additional run-time for the mesh refinement procedure must be considered.In addition to the aforementioned insights specific to the models studied, the experiments performed here have yielded the following outcomes (illustrated in Fig. 7): 1. High-quality meshes (q-refinement only) produce less accurate solutions than goal-oriented refined meshes for the corresponding polynomial basis in case the goaloriented refinement is working adequately (Model 2), even when the refined mesh has a lower quality: for exam-ple for p = 2 elements for Model 2, the Refined mesh (q = 1.4) produces an average NRMSD x of 0.83% and the Initial mesh III (q = 1.2) results in an NRMSD x of 2.15%.This behavior is also seen in the weaker component (E y ).The memory consumption is also lower for the Refined Mesh, as problem sizes are smaller compared to the high-quality initial mesh and the same polynomial basis.Run-times for the same polynomial basis are increased due to the additional run-time consumed for the refinement procedure and the increased number of elements.They are, however, shorter than for initial meshes with a higher polynomial basis order, which also does not necessarily lead to comparably low average NRMSD values.2. A high-order polynomial basis is favorable from a practical perspective.However, its performance compromise depends on the input model and the grid quality.Based on our numerical results, order p = 2 exhibits the best accuracy/performance trade-off (see Fig. 7).This conclusion is consistent with those observed in other applications [8,34,35,45].3. The goal-oriented refinement strategy in conjunction with high-order discretizations and q-refinement is highly practical and beneficial for Model 2).The most pronounced improvement occurs from low-quality meshes with q = 1.6 to refined meshes with q = 1.4 and from p = 1 to p = 2 solutions.
When attempting to invert field data, accurate forward simulations are required.Typical instrument characteristics and noise levels can give us an idea of how accurate the forward solutions actually have to be in order to avoid inverting for errors resulting from numerical inaccuracy of the forward computation.Often, a minimum error floor on any measured field data is set to 2 %, thus there is a desire to model forward responses with a smaller numerical error.The relative noise level for CSEM measurements near the well in the Vallès basin [17] was about 5 %, which implies that our numerical approximations are sufficiently accurate for inverting the data (e.g., NRMSDx, y for refined mesh with p = 2, 3 depicted in the inlays of Figs. 5 and 6 and the average NRMSDx, y in Table 3).The local (close to the well location) and average NRMSDx, y for p = 1 approximations is also below, but quite close, to 5 % for the Refined mesh and for the E x field component of the high-quality Initial mesh II.The conclusion is that inversions with p = 1 elements on those meshes can also work, but it is important to carefully study data errors and modeling accuracy before conducting any inversion tests.The used electric field sensors have a sensitivity, i.e. absolute intrinsic noise of the receiver, of 1 × 10 −9 V/m, thus, smaller changes than this value in the electric field cannot be detected.For Model 2, the E x field component is in the range of 10 −7 V/m and the E y field component is in the range of 10 −9 V/m in the well region (see Figs. 5 and 6).Thus, for the strong field component E x , the observed NRMSD differences between Initial and Refined meshes (about 7 % for Model 2) and lower and higher-order solutions (about 2 % between p-orders) would be in an amplitude range detected by the instruments.

Meshing recommendations
Based on our results, preferable meshing strategies and modeling options can be formulated for models that contain similar features as the ones presented: 1.In Model 1, our approaches are studied to simulate EM responses in the presence of two steel-cased wells discretized with a combination of highly conductive prisms and PECs and no additional complicated subsurface geometries.For this type of well approximation and a simple model structure, it can be seen those p = 2 calculations in combination with fine source and receiver elements simplify the mesh design and offer an excellent balance between numerical accuracy and computational cost.2. Model 2 exhibits large scale variations due to a metallic casing and basin topography.A simplified representation of the steel infrastructure (PEC-only approximation) was implemented.For models of this type, it is recommended to use fine source and receiver elements, apply goal-oriented refinement techniques, and utilize p = 2 elements.3. It is worth it to investigate, whether the model needs a combined prism and PEC representation of the well or whether a PEC-only presentation is sufficient.A PEConly representation is straightforward to implement in a forward modeling code, it reduces problem sizes and does not hamper goal-oriented automatic refinement.Furthermore, it might simplify regularisation approaches in an inversion: As the PEC is implemented as a surface-less representation of the well, smoothing can automatically be imposed across the well.However, it has to be carefully assessed if the PEC approximation influences the inversion result.4. Be aware that for models including steel-cased wells, the most conventional approach using p = 1 elements only might not be sufficient to obtain responses of acceptable accuracy to be used for inversion, even when using a mesh of very good quality.5.It can be recommended to include p-and h-refinement into inversion routines.To avoid an excessive increase in model parameters due to h-refinement, a combined p-and h-refinement approach would require a dual mesh concept, i.e. the separation of forward modeling and inversion parameter meshes.For 3D inversion, hrefinement might be too time consuming to conduct in every iteration.A strategy worth to investigate is to run a first inversion with an initial forward modeling mesh using p = 2 elements and then perform subsequent h-refinement of the forward modeling mesh for the bestfitting model.With this refined mesh and the preferred model from the first inversion run as a start model, a second inversion run with p = 2 elements can be started, accounting for the updated resistivity distribution of the subsurface with the improved forward modeling mesh.

Conclusions
Our results provide evidence that goal-oriented refinement and high-order elements can work together well.However, the h-refinement applied here is based on p = 1 approximation, which can lead to over-refined meshes for p = 2 and p = 3 simulations.Our results can be a motivation for the implementation of goal-oriented refinement adapted for higher order elements, which can probably lead to even better adapted meshes [34,43,45].The focus of the examples presented in this work is to investigate goal-oriented refinement and high-order elements for meshing models with steel-cased wells using variations of PEC approximations.Future studies should focus on investigating the interplay between the implemented methods for models with additional subsurface anomalies, such as anisotropy, as well as more complex source and receiver distributions.Furthermore, steel casings are modelled as nonmagnetic in this study.However, production wells often have a non-negligible magnetic permeability and measured electromagnetic fields can be substantially impacted by magnetic permeabilities of the wells for frequencies > 10 Hz [23].Based on the low frequencies (2 Hz and 5 Hz) used in the presented models to image the deep subsurface, we anticipate that our meshing recommendations will remain unchanged for magnetic steel casings.For higher frequencies in combination with magnetic casings, it has to be investigated, if a PEC approximation is still applicable.

Fig. 2
Fig. 2 Model 1: 3D model in the presence of metallic-cased wells; Borehole casings are modelled with a combination of 100 m long solid rectangular prisms (0 − 100 m depth) with a diameter of 0.2 m and 400 m long perfect electric conductors (PECs, 100 − 500 m depth) and x = 900 m with a spacing of 5 m in the central region and 20 m in the outer region.The entire modeling domain is 60 × 60 × 60 km.

Fig. 4
Fig. 4 Model 2: Vallès basin model; The borehole casing at profile distance 1.719 km is modelled with a 200 m long PEC (0 − 200 m depth).A horizontal x-directed dipole source is located at x = 1.209 km provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX partially funded by the Swedish Research Council through grant agreement No. SNIC 2021/22-883.This work benefited from valuable suggestions and comments of Dra.Pilar Queralt (University of Barcelona).In addition, O.C-R thanks the support of Dra.Pilar who shared the experimental data of Model 2 (obtained in the scope of the GEO-URBAN project under Grant PCI2018-092943).FundingOpen Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.The work of O.C-R.has received funding from the Ministerio de Educación y Ciencia (Spain) under Project TED2021-131882B-C42.The code development of P.R. has been financed by the Smart Exploration project.Smart Exploration has received funding from the European Union's Horizon 2020 Framework Programme under grant agreement N • 775971.The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX partially funded by the Swedish Research Council through grant agreement N • SNIC 2021/22-883.

Table 1
Summary of the most relevant approaches for numerical representation of metallic-cased wells solid cylindrical prism (electric sources can be placed above or below the well)

Table 2
Statistics for Model 1 (see Fig.2) at 5 Hz.For elfe3D, the refinement run-time (sec) is provided

Table 3
Statistics for Model 2

7
Performance behaviour that summarizes the improvement in average NRMSD versus run-time (modeling and h-refinement) and memory consumption for our meshing refinement strategies.Black circles indicate the preferred meshing option for each model still 2.31.It can be decreased with p = 3 approximations to NRMSD y = 1.43, but this comes with a huge increase in computational cost.Calculations on Initial mesh II and III do not result in average NRMS values lower than 3.26.