Modeling electromagnetics on cylindrical meshes with applications to steelcased wells
Abstract
Simulating direct current resistivity, frequency domain electromagnetics and time domain electromagnetics in settings where steel cased boreholes are present is of interest across a range of applications including welllogging, monitoring subsurface injections such as hydraulic fracturing or carbon capture and storage. In some surveys, wellcasings have been used as “extended electrodes” for near surface environmental or geotechnical applications. Wells are often cased with steel, which has both a high conductivity and a significant magnetic permeability. The large physical property contrasts as well as the large disparity in lengthscales, which are introduced when a steelcased well is in a modeling domain, makes performing an electromagnetic forward simulation challenging. Using this setting as motivation, we present a finite volume approach for modeling electromagnetic problems on cylindrically symmetric and 3D cylindrical meshes which include an azimuthal discretization. The associated software implementation includes modeling capabilities for direct current resistivity, time domain electromagnetics, and frequency domain electromagnetics for models that include variable electrical conductivity and magnetic permeability. Electric and magnetic fields, fluxes, and charges are readily accessible in any simulation so that they can be visualized and interrogated. We demonstrate the value of being able to explore the behaviour of electromagnetic fields and fluxes through examples which revisit a number of foundational papers on direct current resistivity and electromagnetics in steelcased wells. The software implementation is open source and included as a part of the SimPEG software ecosystem for simulation and parameter estimation in geophysics.
keywords:
Finite volume, Partial differential equations, Time domain electromagnetics, Frequency domain electromagnetics, Direct current resistivity, Boreholesnoitemsep,topsep=0pt,parsep=0pt,partopsep=0pt
1 Introduction
A number of geophysical electromagnetic (EM) problems lend themselves to cylindrical geometries. Airborne EM problems over a 1D layered earth or boreholelogging applications fall into this category; in these cases cylindrical modeling, which removes a degree of freedom in the azimuthal component, can be advantageous as it reduces the computation load. This is useful when running an inversion where many forward modelings are required, and it is also valuable when exploring and building up an understanding of the behaviour of electromagnetic fields and fluxes in a variety of settings, such as the canonical model of an airborne EM sounding over a sphere, as it reduces feedback time between asking a question and visualizing results (e.g. Oldenburg et al. (2017)).
Beyond these simple settings, there are also a range of scenarios where the footprint of the survey is primarily cylindrical, but 2D or 3D variations in the physical property model may be present. For example if we consider a single sounding in an Airborne EM survey, the primary electric fields are rotational and the magnetic fields are poloidal, but the physical property model may have lateral variations or compact targets. More flexibility is required from the discretization to capture these features. In this case, a 3D cylindrical geometry, which incorporates azimuthal discretization may be advantageous. It allows finer discretization near the source where we have the most sensitivity and the fields are changing rapidly. Far from the source, the discretization is coarser, but it still conforms to the primary behaviour of the EM fields and fluxes and captures the rotational electric fields and poloidal magnetic flux.
In other cases, the most significant physical property variations may conform to a cylindrical geometry, for example in settings where vertical metallic wellcasings are present, or in the emerging topic of using geophysics to “look ahead” of a tunnel boring machine. In particular, understanding the behavior of electromagnetic fields and fluxes in the presence of steelcased wells is of interest across a range of applications, from characterizing lithologic units with welllogs (Kaufman, 1990; Kaufman and Wightman, 1993; Augustin et al., 1989), to identifying marine hydrocarbon targets (Kong et al., 2009; Swidinsky et al., 2013; Tietze et al., 2015), to mapping changes in a reservoir induced by hydraulic fracturing or carbon capture and storage (Pardo and TorresVerdin, 2013; Börner et al., 2015; Um et al., 2015; Weiss et al., 2016; Hoversten et al., 2017; Zhang et al., 2018). Carbon steel, a material commonly used for borehole casings, is highly electrically conductive ( S/m) and has a significant magnetic permeability ( ) (Wu and Habashy, 1994); it therefore can have a significant influence on electromagnetic signals. The large contrasts in physical properties between the casing and the geologic features of interest, along with the large range of scales that need to be considered to model both the millimeterthick casing walls while also capturing geologic features, provide interesting challenges and context for electromagnetics in cylindrical geometries. As such, we will use EM simulations of conductive, permeable boreholes as motivation throughout this paper.
In much of the early literature, the casing was viewed as a nuisance which distorts the EM signals of interest. Distortion of surface direct current (DC) resistivity and induced polarization (IP) data, primarily in hydrocarbon settings, was examined in (Wait, 1983; Holladay and West, 1984; Johnston et al., 1987) and later extended to grounded source EM and IP in (Wait and Williams, 1985; Williams and Wait, 1985; Johnston et al., 1992). Also in hydrocarbon applications, welllogging in the presence of steel cased boreholes is motivation for examining the behavior of electromagnetic fields and fluxes in the vicinity of casing. Initial work focussed on DC resistivity with (Kaufman, 1990; Schenkel and Morrison, 1990; Kaufman and Wightman, 1993; Schenkel and Morrison, 1994), and inductive source frequency domain experiments with (Augustin et al., 1989). Kaufman (1990) derives an analytical solution for the electric field at DC in an experiment where an electrode is positioned along the axis of an infinitelength well. The mathematical solutions presented shows how, and under what conditions, horizontal currents leak into the formation outside the well. Moreover, Kaufman (1990) showed, based upon asymptotic analysis, which fields to measure inside the well so that information about the formation resistivity could be obtained. This analysis is extended to include finitelength wells in Kaufman and Wightman (1993). Schenkel and Morrison (1994) show the importance of considering the length of the casing in borehole resistivity measurements, and demonstrate the feasibility of crosswell DC resistivity. They also show that the presence of a steel casing can improve sensitivity to a target adjacent to the well. In frequency domain EM, Augustin et al. (1989) consider a looploop experiment, where a large loop is positioned on the surface of the earth and a magnetic field receiver is within the borehole. Magnetic permeability is included in the analysis and a “casing correction”, effectively a filter due to the casing on inductivesource data, is introduced. This work was built upon for considering crosswell frequency domain EM experiments (Uchida et al., 1991; Wilt et al., 1996).
For larger scale geophysical surveys, steel cased wells have been used as “extended electrodes.” Rocroi and Koulikov (1985) used a pair of well casings as current electrodes for reservoir characterization in hydrocarbon applications. In nearsurface settings (Ramirez et al., 1996; Rucker et al., 2010; Rucker, 2012) considered the use of monitoring wells as current and potential electrodes for a DC experiment aimed at imaging nuclear waste beneath a leaking storage tank. Similarly, Ronczka et al. (2015) considers the use of groundwater wells for monitoring a saltwater intrusion and investigates numerical strategies for simulating casings as long electrodes. Imaging hydraulic fractures has been a motivator for a number of studies at DC or EM, among them Weiss et al. (2016); Hoversten et al. (2017). Some of these have suggested the use of casings that include resistive gaps so that currents may be injected in a segment of the well and potentials measured across the other gaps along the well (Nekut, 1995; Zhang et al., 2018). There is also interest in modeling casings for casing integrity applications where the aim of the DC or EM survey is to diagnose if a well is flawed or intact based on data collected on the surface (Wilt et al., 2018).
As computing resources increased, our ability to forwardsimulate more complex scenarios has improved. However, the large physical property contrasts and disperate length scales introduced when a steel cased well is included in a model still present computational challenges. Even the DC problem, which is relatively computationally light, has posed challenges; those are exacerbated when solving the full Maxwell equations in the frequency (FDEM) or time domain (TDEM) and can become crippling for an inversion. For models where the source and borehole are axisymmetric, cylindrical symmetry may be exploited to reduce the dimensionality, and thus the number of unknowns, in the problem (e.g. Pardo and TorresVerdin (2013); Heagy et al. (2015)).
To reduce computational load in a 3D simulation, a number of authors have employed simplifying assumptions. Several authors replaced the steelcased well with a solid borehole, either with the same conductivity as the hollowcased well (e.g. Um et al. (2015); Puzyrev et al. (2017)) or preserving the cross sectional conductance (e.g. Swidinsky et al. (2013); Kohnke et al. (2017)), so that a coarser discretization may be used; Haber et al. (2016) similarly replaces the borehole with a coarser conductivity structure and adopts an OcTree discretization locally refine the mesh around the casing. Yang et al. (2016) uses a circuit model and introduces circuit components to account for the steel cased well in a 3D DC resistivity experiment; Weiss (2017) adopts a similar strategy in a Finite Element scheme. Another approach has been to replace the well with an “equivalent source”, for example, a collection of representative dipoles, inspired from Cuevas (2014), or with a linecharge distribution for a DC problem (Weiss et al., 2016). For the frequency domain electromagnetic problem, a method of moments approach, which replaces the casing with a series of current dipoles, has been taken in Kohnke et al. (2017).
For 3D survey geometries, only a handful of forward simulations which accurately discretize the casing have been demonstrated, and they have been achieved at significant computational cost. Recent examples, including Commer et al. (2015); Um et al. (2015); Puzyrev et al. (2017), perform time and frequency domain simulations with finelydiscretized boreholes; they required the equivalent of days of computetime for a single forward simulation to complete. While these codes will undoubtedly see improvements in efficiency, what we present here is an alternative approach to the discretization which capitalizes on the cylindrical geometry of a borehole. Thus far, the majority of the literature has focussed on electrical conductivity, with little attention being paid to magnetic permeability, leaving open many questions about how it alters the behavior of the fields and fluxes and the impact it has on data measured at the surface. In our approach, we ensure that variable magnetic permeability can be included in order to facilitate exploration of these questions.
We introduce an approach and associated opensource software implementation for simulating Maxwell’s equations over conductive, permeable models on 2D and 3D cylindrical meshes. The software is written in Python (Van Rossum and Drake Jr, 1995) and is included as an extension to the SimPEG ecosystem (Cockett et al., 2015; Heagy et al., 2017). Within the context of current research connected to steelcased wells, our aim with the development and distribution of this software is twofold: (1) to facilitate the exploration of the physics of EM in these largecontrast settings, and (2) to provide a simulation tool that can be used for testing other EM codes. The large physical property contrasts in both conductivity and permeability means the physics is complicated and often nonintuitive; as such, we ensure that the researcher can readily access and visualize fields, fluxes, and charges in the simulation domain. This is particularly useful when the software is used in conjunction with Jupyter notebooks which facilitate exploration of numerical results (Perez et al., 2015). As the mesh conforms to the geometry of a vertical borehole, a fine discretization can be used in its vicinity without resulting in a onerous computation. This provides the opportunity to build an understanding of the physics of EM in settings with vertical boreholes prior to moving to settings with deviated and horizontal wells. We demonstrate the software with examples at DC, in the frequency domain, and in the time domain. Sourcecode for all examples is provided as Jupyter notebooks at https://github.com/simpegresearch/heagy2018emcyl (Heagy, 2018); they are licensed under the permissive MIT license with the hope of reducing the effort necessary by a researcher to compare to or build upon this work.
Our paper is organized in the following manner. In section 2, we introduce the governing equations, Maxwell’s equations, and describe their discretization in cylindrical coordinates. We then compare our numerical implementation to the finite element and finite difference results shown in (Commer et al., 2015) as well as a finite volume OcTree simulation described in (Haber et al., 2007). Section 3 contains numerical examples of the DC, frequency domain EM, and time domain EM implementations. The two DC resistivity examples (sections 3.1 and 3.2) are built upon the foundational work in (Kaufman, 1990; Kaufman and Wightman, 1993) which use asymptotic analysis to draw conclusions about the behavior of the electric fields, currents, and charges for a well where an electrode has been positioned along its axis. The next example, in section 3.3, is motivated by the interest in using a steelcased well as an “extended electrode” in a time domain EM experiment. We perform a “topcasing” experiment, with one electrode connected to the top of the well and examine the currents in the surrounding geologic formation through time. Our final two examples, in sections 3.4 and 3.5, consider a frequency domain experiment inspired by (Augustin et al., 1989). These examples demonstrate the impact of magnetic permeability on the character of the magnetic flux within the vicinity of the borehole and discusses the resulting magnetic field measurements made within a borehole.
2 Methodology
2.1 Governing Equations
The governing equations under consideration are Maxwell’s equations. Here we provide a brief overview and recommend Ward and Hohmann (1988) for more detail. Under the quasistatic approximation, Maxwell’s equations are given by:
(1) 
where is the electric field, is the magnetic flux density, is the magnetic field, is the current density and is the source current density. Maxwell’s equations can also be formulated in the frequency domain, using the Fourier Transform convention, they are
(2) 
The fields and fluxes are related through the physical properties: electrical conductivity (, or its inverse, resistivity ) and magnetic permeability (), as described by the constitutive relations
(3) 
At the zerofrequency limit, we also consider the DC resistivity experiment, described by
(4) 
where is the magnitude of the source current density, and are the locations of the current electrodes, and is the scalar electric potential.
Of our numerical tools, we require the ability to simulate large electrical conductivity contrasts, include magnetic permeability, and solve Maxwell’s equations at DC, in frequency and in time in a computationally tractable manner. Finite volume methods are advantageous for modeling large physical property contrasts as they are conservative and the operators “mimic” properties of the continuous operators, that is, the edge curl operator is in the null space of the face divergence operator, and the nodal gradient operator is in the null space of the edge curl operator (Hyman and Shashkov, 1999). As such, they are common practice for many electromagnetic simulations (e.g. Horesh and Haber (2011); Haber and Ruthotto (2018); Jahandari and Farquharson (2014) and references within), and will be our method of choice.
2.2 Finite Volume Discretization
To represent a set of partial differential equations on the mesh, we use a staggeredgrid approach (Yee, 1966) and discretize fields on edges, fluxes on faces, and physical properties at cell centers, as shown in Figure 1. Scalar potentials can be discretized at cell centers or nodes. We consider both cylindrically symmetric meshes and fully 3D cylindrical meshes; the anatomy of a finite volume cell for these scenarios is shown in Figure 1 (b) and (c).
To discretize Maxwell’s equations in the time domain (equation 1) or in the frequency domain (equation 2), we invoke the constitutive relations to formulate our system in terms of a single field and a single flux. This gives a system in either the electric field and magnetic flux (EB formulation), or the magnetic field and the current density (HJ formulation). For example, in the frequency domain, the EB formulation is
(5) 
and the HJ formulation is
(6) 
where are vectors of the discrete EM fields and fluxes; and are the discrete magnetic and electric source terms, respectively; is the edge curl operator, and the matrices are the edge / face inner product matrices. In particular, variable electrical conductivity and variable magnetic permeability are captured in the discretization. The time domain equations are discretized in the same manner as is discussed in (Heagy et al., 2017); for timestepping, a firstorder backward Euler approach is used. Although the midpoint method, which is secondorder accurate, could be considered, it is susceptible to oscillations in the solution, which reduce the order of accuracy, unless a sufficiently small timestep is used (Haber et al., 2004; Haber and Ruthotto, 2018).
At the zerofrequency limit, each formulation has a complementary discretization for the DC equations, for the EB formulation the discretization leads to a nodal discretization of the electric potential , giving
(7) 
where is the nodal gradient operator, and is the source term, defined on nodes. Note that the nodal gradient takes the discrete derivative of nodal variables, and thus the output is on edges. The HJ formulation leads naturally to a cell centered discretization of the electric potential
(8) 
Where is the face divergence operator, is a diagonal matrix of the cell volumes, is the source term, which is defined at cell centers as is . Here, the face divergence takes the discrete derivative from faces to cell centers, thus its transpose takes a variable from cell centers to faces. For a tutorial on the finite volume discretization of the DC equations, see (Cockett et al., 2016).
For the EM simulations, natural boundary conditions are employed; in the EB formulation, this means , and in the HJ formulation, we use . Within the DC simulations, there is flexibility on the choice of boundary conditions employed. In the simplest scenario, for the nodal discretization, we use Neumann boundary conditions, , and for the cell centered discretization, we use Dirichlet boundary conditions .
When employing a cylindrical mesh, the distinction between where the electric and magnetic contributions are discretized in each formulation has important implications. If we consider the cylindrically symmetric mesh (Figure 1b) and a magnetic dipole source positioned along the axis of symmetry (sometimes referred to as the TE mode), we must use the EB formulation of Maxwell’s equation to simulate the resulting toroidal magnetic flux and rotational electric fields. If instead, a vertical current dipole is positioned along the axis of symmetry (also referred to as the TM mode), then the HJ formulation of Maxwell’s equations must be used in order to simulate toroidal currents and rotational magnetic fields. The advantage of a fully 3D cylindrical mesh provides additional degrees of freedom, with the discretization in the azimuthal direction, allowing us to simulate more complex responses. However, in order to avoid the need for very fine discretization in the azimuthal direction, we should select the most natural formulation of Maxwell’s equations given the source geometry being considered. For a vertical steel cased well and a grounded source, we expect the majority of the currents to flow vertically and radially, thus the more natural discretization to employ is the HJ formulation of Maxwell’s equations.
Haber and Ruthotto (2018) provides derivations and discussion of the differential operators and inner product matrices; though they are described for a cartesian coordinate system and a rectangular grid, the extension to a three dimensional cylindrical mesh is straightforward. Effectively, a cartesian mesh is wrapped so that the components become components, and components become components, as shown in Figure 2.
The additional complications that are introduced are: (1) the periodic boundary condition introduced on boundary faces and edges in the azimuthal direction, (2) the removal of radial faces and azimuthal edges along the axis of symmetry, and (3) the elimination of the degrees of freedom of the nodes and edges at the boundary and as well as the nodes and vertical edges along the axis of symmetry. The implementation of the 3D cylindrical mesh is provided as a part of the discretize package (http://discretize.simpeg.xyz), which is an opensource python package that contains finite volume operators and utilities for a variety of meshtypes. All differential operators are tested for second order convergence and for preservation of mimetic properties (as described in Haber and Ruthotto (2018)). discretize is developed in a modular, objectoriented manner and interfaces to all of the SimPEG forward modeling and inversion routines, thus, once the differential operators have been implemented, they can be readily used to perform forward simulations (Cockett et al., 2015; Heagy et al., 2017). One of the benefits of SimPEG for forward simulations is that values of the fields and fluxes are readily computed and visualized, which enables researchers to examine the physics as well as to simulate data. Development within the SimPEG ecosystem follows best practices for modern, openssource software, including: peer review of code changes and additions, versioning, automated testing, documentation, and issue tracking.
2.3 Validation
Testing for the DC, TDEM, and FDEM implementations includes comparison with analytic solutions for a dipole in a wholespace. These examples are included as supplementary examples with the distributed notebooks. We have also compared the cylindrically symmetric implementation at low frequency with a DC simulations from a Resistor Network solution developed in MATLAB with (Figure 3 in Yang et al. (2016)).
Here, we include a comparison with the time domain electromagnetic simulation shown in Figures 13 and 14 of Commer et al. (2015). A 200m long well, with a conductivity of S/m, outer diameter of 135 mm, and casing thickness of 12 mm is embedded in a 0.0333 S/m background. For the material inside the casing, we use a conductivity equal to that of the background. The conductivity of the air is set to S/m and the permeability of the casing is ignored (). A 10 m long inline electric dipole source is positioned on the surface, 50 m radially from the well. The radial electric field is sampled at 5 m, 10 m, 100 m, 200 m and 300 m along a line from the source.
Two simulations are included in Commer et al. (2015): a finite element (FE) and a finite difference (FD) solution. Both simulation meshes capture the thickness of the casing with a single cell or single tetrahedral element. Additionally, we include a comparison with the 3D UBC finite volume OcTree time domain code (Haber et al., 2007). The OcTree mesh allows for adaptive refinement of the mesh around sources, receivers, and conductivity structures within the domain, thus reducing the number of unknowns in the domain as compared to a tensor mesh.
For the 3D cylindrical simulation (SimPEG), we use a mesh that has 4 cells radially across the width of the casing, 2.5m vertical discretization, and azimuthal refinement near the source and receivers (along the line), as shown in Figure 3. To solve the systemmatrix, the direct solver PARDISO was used (Petra et al., 2014; Cosmin et al., 2016). The simulation took 14 minutes to run on a single Intel Xeon X5660 processor (2.80GHz). The details of each simulation are shown in Table 1.
Mesh  Timestepping  Compute Resources  Compute Time  

Commer FE  8 421 559 tetrahedral elements  893 time steps  
9 factorizations  single core  
Intel Xeon X5550 (2.67 GHz)  63 hours  
Commer FD  2 182 528 cells  s  
120 598 277 timesteps  512 cores  
Intel Xeon (2.33 GHz)  23.2 hours  
UBC OcTree  5 011 924 cell  154 time steps  
10 factorizations  single core  
Intel Xeon X5660 (2.80 GHz)  57 minutes  
SimPEG  314 272 cells  187 timesteps  
7 factorizations  single core  
Intel Xeon X5660 (2.80GHz)  14 minutes 
In Figure 4, we show the absolute value of the radial electric field sampled at fives stations; each of the different line colors is associated with a different location, and offsets are with respect to the location of the well. Solutions were interpolated to the same offset using nearest neighbor interpolation.The 3D cylindrical simulation (SimPEG) is plotted with a solid line and overlaps with the UBC solution (dashdot line) for all times shown. The finite element (FE) solution from Commer et al. (2015) is shown with the dashed lines, and the finite difference (FD) solution is plotted with dotted lines. The 3D cylindrical (SimPEG) and UBC solutions are overall in good agreement with the solutions from Commer et al. (2015). There is a difference in amplitude and position of the zerocrossing (the vshape visible in the blue and orange curves) between the Commer solutions and the SimPEG / UBC solutions at the shortest two offsets in the early times. At such short offsets from a highly conductive target, details of the simulation and discretization, such as the construction of the physical property matrices in each of the various approaches become significant; this likely accounts for the discrepancies but a detailed codecomparison is beyond the scope of this paper. Our aim with this comparison is to provide evidence that our numerical simulation is performing as expected; the overall agreement with Commer’s and UBC’s results is provides confidence that it is.
This example demonstrates agreement between the 3D cylindrical solution and solutions obtained with independently developed codes. Importantly, it also shows how, by using a cylindrical discretization which conforms to the conductivity structure of interest, the size of the mesh and resultant cost of the computation can be greatly reduced. This is true even with relatively conservative spatial and temporal discretizations. Minimizing computation time was not a main focus in the development of the software and there are still opportunities for improving efficiency. As an opensource project, contributions from the wider community are encouraged.
3 Numerical Examples
We demonstrate the implementation through examples using the DC, time domain EM and frequency domain EM codes. To focus discussion, each of the examples explores an aspect of the physical behavior of electromagnetic fields and fluxes in the presence of a steelcased well.
3.1 DC Resistivity Part 1: Electric fields, currents and charges in a long well
In his two seminal papers on the topic, Kaufman uses transmission line theory to draw conclusions about the behaviour of the electric field when an electrode is positioned inside of an infinite casing. In this first example, we will revisit some of the physical insights discussed in (Kaufman, 1990; Kaufman and Wightman, 1993) that followed from an analytical derivation and compare those to our numerical results. In the second example, we look at the distribution of current and charges as the length of the well is varied and compare those to the analytical results discussed in (Kaufman and Wightman, 1993)
We start by considering a 1km long well ( S/m) in a whole space ( S/m), with the conductivity of the material inside the borehole equal to that of the whole space. For modeling, we will use a cylindrically symmetric mesh. The positive electrode is positioned on the borehole axis in the midpoint of a 1km long well; a distant return electrode is positioned 1km away at the same depth.
Kaufman discusses the behavior of the electric field by dividing the response into three zones: a near zone, an intermediate zone and a far zone (Kaufman, 1990; Kaufman and Wightman, 1993). In the near zone, the electric field has both radial and vertical components, negative charges are present on the inside of the casing, and positive charges are present on the outside of the casing. The near zone is quite localized and typically, its vertical extent is no more than borehole radii away from the electrode. To examine these features in our numerical simulation, we have plotted in Figure 5: (a) the total charge, (b) secondary charges, (c) electric field, and (d) current density in a portion of the model near the source. The behaviours expected by Kaufman are consistent with our numerical results.
Within the nearzone, the total charge is dominated by the large positive charge at the current electrode location and negative charges that exist along the casing wall where current is moving from a resistive region inside the borehole into a conductor. The extent of the negative charges along the inner casing wall is more evident when we look at the secondary charge, which is obtained by subtracting the charge that would be observed in a uniform halfspace from the total charge (Figure 5b). Inside the casing, we can see the transition from nearzone behavior to intermediate zone behavior approximately 0.5 m above and below the source; that is equal to 10 borehole radii from the source location, which agrees with Kaufman’s conclusion.
In the intermediate zone, Kaufman discusses a number of interesting aspects with respect to the behavior of the electric fields and currents which we can compare with the observed behavior in Figure 5. Among them, he shows that the electric field within the borehole and casing is directed along the vertical axis; as a result no charges accumulate on the inner casing wall. Charges do, however, accumulate on the outer surface of the casing; these generate radiallydirected electric fields and currents, often referred to as âleakage currentsâ, within the formation. At each depth slice through the casing and borehole, the electric field is uniform, however, due to the high conductivity of the casing, most of the current flows within the casing. The vertical extent of the intermediate zone depends on the resistivity contrast between the casing and the surrounding formation and extends beyond several hundred meters before transitioning to the far zone, where the influence of the casing disappears (Kaufman, 1990).
The radially directed fields from the casing, and the length of the intermediate zone, have practical implications in the context of welllogging because they delineate the region in which measurements can be made to acquire information about the formation resistivity outside the well. Within the intermediate zone, fields behave like those due to a transmission line (Kaufman, 1990), and multiple authors have adopted modeling strategies that approximate the well and surrounding medium as a transmission line (Kong et al., 2009; Aldridge et al., 2015). We will extend this analysis in the next example and discuss how the length of the well impacts the behavior of the charges, fields, and fluxes.
3.2 DC Resistivity Part 2: Finite Length Wells
In (Kaufman and Wightman, 1993), the transmissionline analysis was extended to consider finitelength wells. Inspired by the interest in using the casing as an “extended electrode” for delivering current to depth (e.g. Schenkel and Morrison (1994); Um et al. (2015); Weiss et al. (2016); Hoversten et al. (2017)), here we consider a 3D DC resistivity experiment where one electrode is connected to the top of the well. We will examine the current and charge distribution for wells ranging in length from 250 m to 4000 m and compare those to the observations in (Kaufman and Wightman, 1993). The conductivity of the well is selected to be S/m. A uniform background conductivity of S/m is used and the return electrode is positioned 8000m from the well; this is sufficiently far from the well that we do not need to examine the impact of the return electrode location in this example. A 3D cylindrical mesh was used for the simulation.
Kaufman and Wightman (1993) derives a solution for the current within a finite length well and discusses two endmember cases: a short well and a long well. “Short” versus “long” are defined on the product of , where is the length of the casing and , where is the crosssectional conductance of the casing and has units of Sm (, for a casing with radius and thickness ), and is the transverse resistance. The transverse resistance is approximately equal to the resistivity of the surrounding formation (for more discussion on where this approximation breaks down, see Schenkel and Morrison (1994)). For short wells, , the current decreases linearly with distance, whereas for long wells, where , the current decays exponentially with distance from the source, with the rate of decay being controlled by the parameter . In Figure 6 (a), we show current in the well for 5 different borehole lengths. The xaxis is the distance from the source normalized by the length of the well. We also show the two endmember solutions (equations 45 and 53) from Kaufman and Wightman (1993). There is significant overlap between the 250m numerical solution and the short well approximation. As the length of the well increases, exponential decay of the currents becomes evident. Since is quite small, for this example , the borehole must be very long to reach the other end member which corresponds to the exponentially decaying solution.
In Figure 6 (b), we have plotted the charges along the length of the well. In the shortwell regime, the borehole is approximately an equipotential surface and the charges are uniformly distributed; in the long well the charges decay with depth. What was surprising to us was the noticeable increase in charge accumulation that occurs near the bottom of the well. This is especially evident for the short well. Initially, we were suspicious and thought this might be due to problems with our numerical simulation; there was no obvious physical explanation that we were aware of. However, investigation into the literature revealed that the increase in charge density at the ends of a cylinder is a real physical effect, but an exact theoretical solution does still not appear to exist (Griffiths and Li, 1997) (see figure 4, in particular).
The results shown in Figure 6 have implications when testing approaches for reducing computational load by approximating a well with a solid cylinder or prism, as in Um et al. (2015), or replacing the well with a distribution of charges, as in Weiss et al. (2016). For a short well, the behaviour of the currents is independent of conductivity, so, as long as the borehole is approximated by a sufficiently conductive target, the behaviour of the fields and fluxes will be representative of the finescale model. However, as the length of the well increases, the crosssectional conductance of the well becomes relevant as it controls the rate of decay of the currents in the well and thus the rate that currents leak into the formation. A similar result holds when a line of charges is used to approximate the well as a DC source; a uniform charge is suitable for a sufficiently short or sufficiently conductive well, whereas a distribution of charge which decays exponentially with depth needs to be considered for longer wells. Thus, when attempting to replace a finescale model of a well with a coarsescale model, either with a conductivity structure or by some form of “equivalent source”, validations should be performed on models that have the same lengthscale as the experiment to ensure that both behaviors are being accurately modeled.
3.3 Time Domain Electromagnetics
In this example, we examine the behaviour of electric currents in an experiment where the casing is used as an “extended electrode”. Although the initial investigations with casings centered around using a DC source, greater information about the subsurface can be had by employing a frequency or time domain source. A particular application is the monitoring of hydraulic fracturing proppant and fluids, or CO; this is active research carried out by many groups worldwide (e.g. Hoversten et al. (2015); Um et al. (2015); Puzyrev et al. (2017); Zhang et al. (2018) among others). The challenge is to have efficient and accurate forward modeling; solving the full Maxwell equations is much more demanding than solving the DC problem. For our simulation, a positive electrode is connected to the top of the casing and a return electrode is positioned 1 km away. The well has a conductivity of S/m and is 1 km long; it has an outer diameter of 10 cm and a 1 cm thick casing wall. The mud which infills the well has the same conductivity as the background, S/m. The conductivity of the air is set to S/m; in numerical experiments, we have observed that contrasts near or larger than S/m leads to erroneous numerical solutions. For this example, we will focus on electrical conductivity only and set the permeability of the well to . A stepoff waveform is used, and the currents within the formation are plotted through time in Figure 7. Panel (a) shows a zoomedin cross section of the casing, (b) shows a vertical cross section along the line of the wire (c) shows a horizontal depth slice at 50 m depth and (d) shows a depth slice at 800 m depth. The images in panels (b), (c) and (d) are on the same color scale.
We begin by examining Figure 7 (b), which shows the currents in the formation. At time , we have the DC solution. Currents flow away from the well, and eventually curve back to the return electrode. Immediately after shutoff, we see an image current develop in the formation. The image current flows in the same direction as the original current in the wire; this is opposite to currents in the formation, causing a circulation of current. The center of this circulation is visible as the null propagating downwards and to the right in Figure 7 (b). In Figure 7 (a), we see the background circulating currents being channeled into the well and propagating downwards. The depth range over which currents enter the casing depends upon time. At t=0.01 ms, the zero crossing, which distinguishes the depth between incoming and outgoing current in the casing, occurs at z=90 m, at t=0.1 ms it is at 225 m and by t=1 ms, the zero crossing approaches the midway point in the casing and is at 470 m depth. At later times, the downward propagation of this null slows as the image currents are channeled into the highly conductive casing; at 5 ms it is at 520 m depth, at 10 ms, 560 m depth and by 100 ms (not shown), it is at 800 m depth.
On the side of the well opposite to the wire, we also see a null develop; it is visible in the cross sections in panel (a). To help understand this, we examine the depth slices in panel (c). Behind the well, we see that as the image currents diffuse downwards and outwards, some of those currents are channeled back towards the well; this is visible in the depth slice at . These channeled currents are opposite in direction to those the formation currents set up at t=0, which also are diffusing downwards and outwards; where these two processes intersect, there is a current shadow.
There are a number of points to highlight in this example. The first, which has been noted by several authors (e.g. Schenkel and Morrison (1994); Hoversten et al. (2015)), is that the casing helps increase sensitivity to targets at depth. This occurs by two mechanisms: (1) at DC, prior to shutoff, the casing acts as an “extended electrode” leaking current into the formation along its length; (2) after shutoff, it channels the image currents and increases the current density within the vicinity of the casing. The second point to note is that there are several survey design considerations raised by examining the currents: targets that are positioned where there is significant current will be most illuminated. If the target is near the surface and offset from the well, a survey where the source wire runs along the same line as the target will have the added benefits of the excitation due to the image currents. These benefits are twofold: (1) the passing imagecurrent increases the current density for a period of time, and (2) the changing amplitude and direction of the currents with time generate different excitations of the target. this should provide enhanced information in an inversion, as compared to a single excitation that is available from a DC survey. For deeper targets in this experiment, the passing image current has diffused significantly, and thus it appears that the wire location has less impact on the magnitude of the current density with location. However, it is possible that increasing the wirelength could be beneficial. This extension is straightforward and could be examined with the provided script. There may also be added benefit by having the target positioned along the same line as the source wire, as at later times, the direction of current reverses, changing the excitation of the target. The final point to note from this example is that although this is a simple model, the behavior of the currents is not intuitive; visualizations of the currents, fields and fluxes, allow researchers to explore the basic physics and prompts new questions.
3.4 Frequency Domain Electromagnetics Part 1: Comparison with scale model results
In the DC example, we discussed how charges are distributed along the well and currents flow into the formation. The time domain example extended the analysis of grounded sources, showed the potential importance of EM induction effects and illuminated the underlying physics. From a historical perspective, however, practical developments in EM were pursued in the frequency domain; the mathematics is more manageable in the frequency domain, and technological advances were being made in the development of induction welllogging tools (Doll, 1949; Moran and Kunz, 1962). Although conductivity of the pipes is generally plays the most dominant role in attenuating the signal, the magnetic permeability is nonnegligible (Wait and Hill, 1977); it is the product of the conductivity and permeability that appears in the description of EM attenuation. Also, the fact that permeable material becomes magnetized in the presence of an external field complicates the problem. Augustin et al. (1989) is one of the first papers on induction logging in the presence of steel cased wells that aims to understand and isolate the EM response of the steel cased well. Using a combination of scale modeling and analytical mathematical modeling, they examined the impacts of conductivity and magnetic permeability on the magnetic field observed in the pipe. In this example and the one that follows, we attempt to unravel this interplay between conductivity and magnetic permeability.
The first experiment Augustin et al. (1989) discuss is a scale model using two different pipes, a conductive copper pipe and a conductive, permeable iron pipe; each pipe is 9 m in length. The copper pipe had an inner diameter of 0.063 m and a thickness of 0.002 m, while the iron pipe had a 0.063 m inner diameter and 0.0043 m wall thickness. A sourceloop, with radius 0.6 m was coaxial with the pipe and in one experiment positioned at one end of the pipe (which they refer to as the “semiinfinite pipe” scenario). In another experiment the source loop is positioned at the midpoint of the pipe (which they refer to as the “infinite pipe” scenario); for both experiments, magnetic field data are measured as a function of frequency at the central axis of the pipe. Their results are presented in terms of a Field Strength Ratio (FSR), which is the ratio of the absolute value of the magnetic field at the receiver with the absolute value of the magnetic field if no pipe is present (Figure 3 in Augustin et al. (1989)). At low frequencies, for the data collected within the iron pipe, static shielding (FSR 1) was observed for the measurements where the receiver was in the plane of the source loop for both the “infinite” and “semiinfinite” scenarios. When the receiver was positioned within the pipe, 1.49 m offset from the plane of the source loop, static enhancement effects (FSR 1) were observed for both the infinite and semiinfinite scenarios. Using this experiment for context, we will compare the behaviour of our numerical simulation with the observations in (Augustin et al., 1989) and examine the nature of the static shielding and enhancement effects.
For our numerical setup, the pipes are 9 m in length and have an inner diameter of 0.06 m. The copper pipe has a casingwall thickness of 0.002 m and the iron pipe has a thickness of 0.004 m. Following the estimated physical property values from Augustin et al. (1989), we use a conductivity of S/m and a relative permeability of 1 for the copper pipe. For the iron pipe, a conductivity of S/m and a relative permeability of 150 is used. A background conductivity of m is assumed. The computed FSR values for the axial magnetic field as a function of frequency are shown in Figure 8.
Consider the response of the conductive pipe. At low frequencies, the FSR for the copper pipe (blue lines) is 1 for both the infinite (solid line) and semiinfinite (dashed line) scenarios, as the field inside the copper pipe is equivalent to the freespace field. With increasing frequency, eddy currents are induced in the pipe which generate a magnetic field that opposes the primary, causing a decrease in the observed FSR. When the source and receiver are in the same plane (L=0.00 m), the rate of decrease is more rapid in the infinite scenario than the semiinfinite. Since there is conductive material on both sides of the receiver in the infinite case, we would expect attenuation of the fields to occur more rapidly than in the semiinfinite case. This observation is consistent with Figure 3a in Augustin et al. (1989). For the offset receiver (L=1.49 m), they observed a slight separation in the infinite and semiinfinite curves which we do not; however, they attributed this to potential errors in magnetometer position. Thus, overall, the numerical results for the copper pipe are in good agreement with the scale model results observed by Augustin et al. (1989).
Next, we examine the response of the conductive, permeable pipe. In Figure 8b, we observe a static enhancement effect (FSR 1) at low frequencies. The enhancement is larger in the infinite scenario than the semiinfinite scenario; this is in agreement with Figure 3b in Augustin et al. (1989). There is however, a significant discrepancy between our numerical simulations and the scale model for the semiinfinite pipe when the source and receiver lie in the same plane(Figure 8a). Augustin et al. (1989) observed a static shielding effect for both the infinite and semiinfinite scenarios, whereas we observe a static shielding for the infinite scenario, but a significant static enhancement for the semiinfinite case. To diagnose what might be the cause of this, we will examine the magnetic flux density in this region of the pipe.
In Figure 9, we have plotted: (a) the secondary magnetic flux in the infinitepipe scenario near the source (z=4.5 m), (b) the secondary magnetic flux in the semiinfinite scenario (z=0 m for the source), and (c) top 5 cm of the semiinfinite pipe. All plots are at 0.1 Hz. The primary magnetic field is directed upwards within the regions we have plotted, so upwardgoing magnetic flux indicates a static enhancement effect, and downwardoriented magnetic flux indicates static shielding effects. In (a) we see a transition between the static shielding in the vicinity of the source to a static enhancement approximately 0.5 m above and below the plane of the source. Similarly in (b), we notice a signreversal in the zcomponent of the secondary magnetic flux at a depth of 0.6 m.
The behaviors observed in Figure 9 are quite comparable to Augustin et al.’s observation of a transition from shielding to enhancement occurring at distances greater than 0.8 m from the source. Numerical experiments show that the vertical extent of the region over which static shielding is occuring increases with increasing pipe diameter, and similarly increases with increasing loop radius while the magnitude of the effect decreases. This can be understood by considering how the pipe is magnetized; for a small loop radius, the magnetization is largely localized near the plane of the source and rapidly falls off with distance from the plane of the source. Localized, large amplitude magnetization causes the casing to act as a collection of dipoles around the circumference of the casing. As the radius of the loop increases, the magnetization spreads out along the length of the well resulting in longer, loweramplitude dipoles, thus both increasing the extent of the region over which static shielding is occuring as well as decreasing its amplitude.
This explains the nature of the static enhancement and static shielding effects, but to explain the discrepancy between the static shielding observed in the semiinfinite pipe when L=0 m by Augustin et al., and the static enhancement we observe in Figure 8a, we examine the magnetic flux density in the top few centimeters of the pipe. Figure 9c shows the top 5 cm of the secondary magnetic flux in the semiinfinite pipe; the source is in the z=0 m plane. Zooming in reveals there is yet another sign reversal near the end of the pipe. This is evident even in the infinitepipe scenario (Figure 8d), where the source is offset by several meters from the end of the pipe. This edgeeffect perhaps bears some similarities to what we observed in Figure 6b, where we saw a build up of charge near the end of the pipe in the DC scenario. At the end of the pipe, we encounter the situation where the normal component of the flux () from the pipe to the background needs to be continuous both in the radial and vertical directions at the end of the pipe as does the tangential component of the fields (). The interplay of these two constraints at the end of the pipe results in more complexity in the resultant fields and fluxes. Within the span of a few centimeters we transition from static enhancement at the top of the pipe to a static shielding further down. An error as small as a few centimeters in the position of the magnetometer causes a reversal in behavior; in Figure 10, we have plotted the FSR for a magnetometer positioned 3cm beneath the plane of the source, and the staticshielding behavior observed for the semiinfinite pipe is much more aligned with that observed in Figure 3a in Augustin et al. (1989).
3.5 Frequency Domain Electromagnetics Part 2: Conductivity and permeability in the inductive response of a well
The experiments shown in the previous section revealed some insights into the complexity of the fields within the pipe and illustrated the role of permeability in the character of the responses at low frequency. Next, we move to larger scales and examine the role of conductivity and permeability in the responses we observe in the borehole.
In this example, we consider a 2 km long well with an outer diameter of 10 cm and thickness of 1 cm in a wholespace which has a resistivity of m. A loop with radius 100 m is coaxial with the well and positioned at the topend of the well. A receiver measuring the zcomponent of the magnetic flux density is positioned 500 m below the transmitter loop, along the axis of the well. We will consider both time domain and frequency domain responses.
In electromagnetics, it is often the product of permeability and conductivity that we consider to be the main controlling factor on the EM responses. To assess the contribution of each to the measured responses, we will investigate two scenarios. In the first, the well has a conductivity of S/m and a relative permeability of 1, and in the second, the well has a conductivity of S/m and a relative permeability of 100; thus the product of conductivity and permeability is equivalent for both wells.
Similar to the analysis done by Augustin et al. (1989) when looking at the role of borehole radius in the behaviour of the magnetic response (e.g. figure 8), we will examine the normalized secondary field (NSF) which is the ratio of the secondary field with the amplitude of the primary, where the primary is defined to be the freespace response. In Figure 11, we have plotted the normalized secondary field for the two pipes considered, the conductive pipe (blue) and the conductive, permeable pipe (orange). Let’s start by examining the conductivity response in Figure 11. Where the value of the NSF is zero, the primary dominates the response; this is the case at low frequencies where induction is not yet contributing to the response. As frequency increases, currents are induced in the pipe which generate a secondary magnetic field that opposes the primary, hence the NSF becomes negative. When the real part of the NSF (solid line) is 1, the secondary magnetic field is equal in magnitude but opposite in direction to the freespace primary and the measured real field is zero. Values less than 1 indicate a sign reversal in the real magnetic field. Similarly, when the imaginary part of the response function goes above zero, there is a sign reversal in the imaginary component. Note that these sign reversals occur even in a halfspace and are a result of sampling the fields within a conductive medium; in this case the receiver was 500 m below the surface.
As compared to the conductive pipe, the frequency at which induction sets in is higher for the conductive, permeable pipe. We also notice that the amplitude variation of both the imaginary and real parts is larger for the permeable pipe. To examine the contribution of conductivity and permeability to the responses, we have plotted the real part of the secondary magnetic flux density, , in Figure 12. The top row shows the response within the conductive pipe and the bottom row shows the conductive, permeable pipe. The primary magnetic flux is oriented upwards and we can see that all of the secondary fields generated are oriented downwards. Similar to the previous example, we see that at low frequencies, there is magnetostatic response due to the permeable pipe. However, due to the larger length scales of the source loop and the casing in this example, there is no measurable contribution at the receiver. At 1 Hz, we can see that induction is starting to contribute to the signal for the conductive pipe, while for the permeable pipe, it is not until 10 Hz that we begin to observe the contribution of induction. At 100 Hz, the secondary magnetic field is stronger in amplitude than the primary, and the NFS is less than 1 for both the conductive and permeable pipes. The amplitude of the secondary within the permeable pipe is stronger than that in the conductive pipe. At 1000 Hz, we have reached the asymptote of NSF=1 for both the conductive and permeable pipes; the secondary magnetic flux is equal in magnitude but opposite in direction to the primary.
Conducting a similar experiment in the time domain, we can compare the responses as a function of time. For this experiment, a stepoff waveform is employed and data are measured after shutoff, the NSF is plotted in Figure 13. Note here that the secondary field is in the same direction as the primary, so after the source has been shut off, the secondary field is oriented upwards, as shown in Figure 14. Shortly after shutoff, the rate of increase in the secondary field is the same for both the conductive and the conductive, permeable wells. A maximum normalized field strength of approximately 1 is reached for both cases. The responses begin to differ at s where the conductive well maintains a NFS for approximately 1 ms longer than the permeable well before the fields decay away.
It is important to note that although the product of the conductivity and permeability is identical for these wells, the geometry of the well and inducing fields results in different couplings for each of the parameters. For a vertical magnetic dipole source, the electric fields are purely rotational while the magnetic fields are primarily vertical. An approximation we can use to understand the implications of these geometric difference is to assume the inducing fields are uniform (e.g. the radius of the source loop is infinite) and to examine the conductance and permeance of the pipe. For rotational electric fields, the conductance is
(9) 
where is the thickness of the casing, is the radius of the casing and is the lengthscale of the pipe segment contributing to the signal. For vertical magnetic fields, the permeance is
(10) 
As the lengthscale, , is larger than the circumference of the pipe () the geometric contribution to the conductance is larger than that to the permeance.
An important takeaway from this example is that the contributions of conductivity and permeability to the observed EM signals are not simply governed by their product. The geometry of the source fields plays an important role in how each contributes. Thus to accurately model conductive, permeable pipes, over a range of frequencies or times, a numerical code must allow both variable conductivity and variable permeability to be considered.
4 Discussion
The behavior of EM fields and fluxes in the presence of highly conductive, permeable, steelcased wells is often nonintuitive. In DC resistivity experiments, we demonstrated that there is a charge buildup near the end of the well. This has not previously been discussed in the geophysics literature, but it was recognized by Griffiths and Li (1997). In practice, such a charge buildup might be consequential in an inversion as it will alter the sensitivities; this is an avenue for future research.
Moving to EM experiments complicates matters in two ways: (1) the fields and fluxes vary through time introducing inductive phenomena, and (2) variable magnetic permeability alters the fields and fluxes. With respect to magnetic permeability, Augustin et al. (1989) noticed “static shielding” and “static enhancement” effects in a scalemodel experiment with an iron pipe subject to an inductive source. They observed that “the truncated pipe is more effective at shielding static and lowfrequency fields than the infinite pipe,” however, they offered no explanation as to why this is the case. In the numerical experiment in section 3.4, we were able to replicate the nature of their data, and furthermore, we demonstrated that near the ends of the pipe, the magnetic fields change rapidly over short lengthscales. Although this may seem to be a detail that is unimportant unless measurements are being made near the end of a pipe, it demonstrates some of the complexity that is introduced when both conductivity and permeability are significant in a model. To date, the geophysics literature that considers the magnetic permeability of wellcasings does so primarily in the context of inductive sources; very little research examines the impact of permeability on groundedsource experiments. Thus, there are open questions about how magnetic permeability alters the currents and impacts the data in such experiments. Building an understanding of these impacts is critical both for assessing feasibility of using EM methods to delineate targets of interest as well as for developing strategies to reduce the computational cost of 3D simulations which include steelcased wells.
In many settings where DC or EM experiments are being considered, wells are deviated or horizontal, and several wells may be present. The cylindrical discretization strategy presented here does not accommodate such geometries. Recent advancements such as the hierarchical finite element approach presented in Weiss (2017) make modeling these complex scenarios feasible for DC resistivity. However, it is not clear that approaches that work at DC will be suitable for EM. For example, Weiss (2017) points out several challenges that arise when considering the inherently more complex PDE governing EM, and in particular, points out that it is unclear how to include magnetic permeability in a hierarchical approach. For EM, other avenues for performing simulations include the upscaling approach suggested by Schwarzbach and Haber (2018); CaudilloMata et al. (2017) in which one inverts for the conductivity and permeability of a coarsescale model of the casing, as well as the method of moments approach taken by Kohnke et al. (2017). Irrespective of the strategy that is taken, it is important to have numerical tools that yield accurate, computationally efficient simulations and that are easy to use. Although tools such as COMSOL are versatile, the cognitive overhead for a researcher to set up a simple simulation to test their understanding of the physics is significant. The software presented here aims to bridge that gap and serve as a resource for researchers to calibrate their understanding of the physics, as well as to assess the assumptions that new approaches are making and to benchmark their accuracy.
5 Summary and Outlook
We developed software for solving Maxwell’s equations on 2D and 3D cylindrical meshes. Variable electrical conductivity and magnetic permeability are considered. The 2D solution is especially computationally efficient and has a large number of practical applications. When cylindrical symmetry is not valid, the 3D solution can be implemented; a judicious design of the mesh can often generate a problem with fewer cells than would be required with a tensor or OcTree mesh, thus reducing the computational cost of a simulation. We demonstrated the versatility of the codes by modeling the electromagnetic fields that result when a highly conductive and permeable casing is embedded in the earth.
We presented a number of different experiments involving DC, frequencydomain, and timedomain sources. The first two examples considered a simple DC resistivity experiment. In the first, we demonstrated that the numerically obtained currents, electric fields, and charges emulated those predicted by the asymptotic analysis in Kaufman (1990) for long wells. The second example looked at the transition in behavior of currents and charges between short and long wells. Even in this relatively simple example, the physics was more complex than we originally anticipated.
In the subsequent examples, we considered electromagnetic experiments. The third example presented a groundedsource timedomain experiment and showed the distribution of currents in the formation through time. It showed that the steelcasing can help excite a target at depth by two mechanisms: (1) the casing provides a highconductivity pathway for bringing DC currents to depth, and (2) the casing channels the image current that is created after shutoff of the source. With respect to survey design, one consequence of the second point is that there may be an advantage to positioning the wire and return electrode along the same line as where the target is expected to be located. In this way, the current direction is reversed as the image current passes; the target is thus excited from multiple excitation directions and the resultant data can be beneficial in an inversion.
The final two examples incorporated magnetic permeability in the simulations. We showed that for a conductive and permeable casing, excited by a circular current source, there is a complicated magnetic field that occurs in the top few centimeters of the pipe. Furthermore, the role of conductivity and permeability in the observed responses is more complex than their product; the source geometry and coupling with the casing are important to consider. In each of the examples, the ability to plot the charges, fields, and fluxes was of critical importance; these ground our understanding of the physics and provide a foundation for designing a field survey.
The software implementation is included as a part of the SimPEG ecosystem. SimPEG also includes finite volume simulations on 3D tensor and OcTree meshes as well as machinery for solving inverse problems. Thus, the cylindrical codes can be readily connected to an inversion and additionally, simulations and inversions of more complex 3D geologic settings can be achieved by coupling the cylindrical simulation with a 3D tensor or OcTree mesh using a primarysecondary approach (e.g. example 3 in Heagy et al. (2017)). Beyond modeling steel cased wells, we envision that the 3D cylindrical mesh could prove to be useful in conducting 3D airborne EM inversions where a domaindecomposition approach, similar to that described in Yang et al. (2014), is adopted.
SimPEG and all of the further developments described in this paper are open source and freely available. The examples have been provided as Jupyter notebooks. This not only allows all of the figures in the paper to be reproduced, but provides an avenue by which the reader can ask questions, change parameters, and use resultant images to confirm (or not) their presumed outcome. We hope that our efforts to make the software and examples accessible promotes the utility of this work for the wider community.
6 Acknowledgements
The authors thank Michael Commer and Christoph Schwarzbach for providing the simulation results shown in Figure 4 and for permission to distribute them. We also thank Thibaut Astic and Dikun Yang for their suggestions and input on the early draft of this paper. Finally, we are grateful to Rowan Cockett, Seogi Kang and the rest of the SimPEG community for their discussions and efforts on the development of the SimPEG.
We are also grateful to Raphael Rochlitz and the two other anonymous reviewers for their critiques and suggestions which improved the quality of the manuscript.
The funding for this work is provided through the Vanier Canada Graduate Scholarships Program.
7 Computer Code Availability
All of the software used in this paper is open source and was made available in 2018. The examples are provided as Jupyter notebooks and the code is written in Python. The code and installation instructions are available at https://github.com/simpegresearch/heagy2018emcyl and have been archived at https://doi.org/10.5281/zenodo.1220427. The main developer is Lindsey Heagy (email: lheagy@eos.ubc.ca, phone: (604) 8362715).
References
References

Aldridge et al. (2015)
Aldridge, D. F., Weiss, C. J., Knox, H. A., Schramm, K. A., Bartel, L. C.,
2015. Is a SteelCased Borehole an Electrical Transmission Line? SEG
Technical Program Expanded Abstracts 2015, 736–741.
URL http://library.seg.org/doi/10.1190/segam20155927874.1 
Augustin et al. (1989)
Augustin, A. M., Kennedy, W. D., Morrison, H. F., Lee, K. H., jan 1989. A
theoretical study of surfaceâtoâborehole electromagnetic logging in cased
holes. Geophysics 54 (1), 90–99.
URL http://library.seg.org/doi/10.1190/1.1442581http://library.seg.org/doi/abs/10.1190/1.1442581  Börner et al. (2015) Börner, J. H., Wang, F., Weißflog, J., Bär, M., Görz, I., Spitzer, K., 2015. Multimethod virtual electromagnetic experiments for developing suitable monitoring designs: A fictitious CO2 sequestration scenario in Northern Germany. Geophysical Prospecting 63 (6), 1430–1449.
 CaudilloMata et al. (2017) CaudilloMata, L., Haber, E., Heagy, L., Schwarzbach, C., 2017. A framework for the upscaling of the electrical conductivity in the quasistatic Maxwell’s equations. Journal of Computational and Applied Mathematics 317.

Cockett et al. (2016)
Cockett, R., Heagy, L. J., Oldenburg, D. W., aug 2016. Pixels and their
neighbors: Finite volume. The Leading Edge 35 (8), 703–706.
URL http://library.seg.org/doi/10.1190/tle35080703.1 
Cockett et al. (2015)
Cockett, R., Kang, S., Heagy, L. J., Pidlisecky, A., Oldenburg, D. W., dec
2015. SimPEG: An open source framework for simulation and gradient based
parameter estimation in geophysical applications. Computers & Geosciences
85, 142–154.
URL http://linkinghub.elsevier.com/retrieve/pii/S009830041530056X 
Commer et al. (2015)
Commer, M., Hoversten, G. M., Um, E. S., mar 2015. Transientelectromagnetic
finitedifference timedomain earth modeling over steel infrastructure.
Geophysics 80 (2), E147–E162.
URL http://library.seg.org/doi/pdf/10.1190/geo20140324.1http://library.seg.org/doi/abs/10.1190/geo20140324.1  Cosmin et al. (2016) Cosmin, G., Schenk, O., Lubin, M., Gäertner, K., 2016. An Augmented Incomplete Factorization Approach for Computing the Schur Complement in Stochastic Optimization.

Cuevas (2014)
Cuevas, N. H., sep 2014. Analytical solutions of EM fields due to a dipolar
source inside an infinite casing. Geophysics 79 (5), E231–E241.
URL http://library.seg.org/doi/pdf/10.1190/geo20130223.1http://library.seg.org/doi/abs/10.1190/geo20130223.1  Doll (1949) Doll, H., 1949. Introduction to Induction Logging and Application to Logging of Wells Drilled With Oil Base Mud. Journal of Petroleum Technology 1 (6).

Griffiths and Li (1997)
Griffiths, D. J., Li, Y., 1997. Charge density on a conducting needle.
URL http://aapt.scitation.org/doi/10.1119/1.18671 
Haber et al. (2004)
Haber, E., Ascher, U. M., Oldenburg, D. W., 2004. Inversion of 3D
electromagnetic data in frequency and time domain using an inexact
allâatâonce approach. Geophysics 69 (5), 1216–1228.
URL http://library.seg.org/doi/10.1190/1.1801938  Haber et al. (2007) Haber, E., Heldmann, S., Ascher, U., 2007. Adaptive finite volume method for distributed nonsmooth parameter identification. Inverse Problems 23 (4), 1659–1676.

Haber and Ruthotto (2018)
Haber, E., Ruthotto, L., nov 2018. A multiscale finite volume method for
Maxwell ’ s equations at low frequencies. Geophysical Journal International
m (May), 1268–1277.
URL http://academic.oup.com/gji/article/199/2/1268/616603/AmultiscalefinitevolumemethodforMaxwells 
Haber et al. (2016)
Haber, E., Schwarzbach, C., Shekhtman, R., 2016. Modeling electromagnetic
fields in the presence of casing. SEG Technical Program Expanded Abstracts
2016 (1988), 959–964.
URL http://library.seg.org/doi/10.1190/segam201613965568.1 
Heagy (2018)
Heagy, L., apr 2018. simpegresearch/heagy_2018_emcyl. Zenodo.
URL https://doi.org/10.5281/zenodo.1220427 
Heagy et al. (2017)
Heagy, L. J., Cockett, R., Kang, S., Rosenkjaer, G. K., Oldenburg, D. W., 2017.
A framework for simulation and inversion in electromagnetics. Computers and
Geosciences 107 (July), 1–19.
URL http://dx.doi.org/10.1016/j.cageo.2017.06.018 
Heagy et al. (2015)
Heagy, L. J., Cockett, R., Oldenburg, D. W., Wilt, M., aug 2015. Modelling
electromagnetic problems in the presence of cased wells. In: SEG Technical
Program Expanded Abstracts 2015. Society of Exploration Geophysicists, pp.
699–703.
URL http://library.seg.org/doi/10.1190/segam20155931035.1 
Holladay and West (1984)
Holladay, J. S., West, G. F., 1984. Effect of well casings on surface
electrical surveys. Geophysics 49 (2), 177–188.
URL http://library.seg.org/doi/10.1190/1.1441649 
Horesh and Haber (2011)
Horesh, L., Haber, E., jan 2011. A Second Order Discretization of Maxwell’s
Equations in the QuasiStatic Regime on OcTree Grids. SIAM Journal on
Scientific Computing 33 (5), 2805–2822.
URL http://epubs.siam.org/doi/10.1137/100798508  Hoversten et al. (2015) Hoversten, G. M., Commer, M., Haber, E., Schwarzbach, C., 2015. Hydrofrac monitoring using ground timedomain electromagnetics. Geophysical Prospecting 63 (6), 1508–1526.
 Hoversten et al. (2017) Hoversten, G. M., Schwarzbach, C., Belliveau, P., Haber, E., Shekhtman, R., 2017. Borehole to Surface Electromagnetic Monitoring of Hydraulic Fractures. In: 79th EAGE Conference and Exhibition 2017.
 Hyman and Shashkov (1999) Hyman, J. M., Shashkov, M., 1999. Mimetic Discretizations for Maxwell ’ s Equations 909, 881–909.

Jahandari and Farquharson (2014)
Jahandari, H., Farquharson, C. G., 2014. A finitevolume solution to the
geophysical electromagnetic forward problem using unstructured grids.
Geophysics 79 (6), E287–E302.
URL http://library.seg.org/doi/10.1190/geo20130312.1 
Johnston et al. (1987)
Johnston, R., Trofimenkoff, F., Haslett, J., jul 1987. Resistivity Response of
a Homogeneous Earth with a FiniteLength Contained Vertical Conductor. IEEE
Transactions on Geoscience and Remote Sensing GE25 (4), 414–421.
URL http://ieeexplore.ieee.org/document/4072660/  Johnston et al. (1992) Johnston, R. H., Trofimenkoff, F. N., Haslett, J. W., 1992. The Complex Resistivity Response of a Homogeneous Earth with a FiniteLength Contained Vertical Conductor. IEEE Transactions on Geoscience and Remote Sensing 30 (1), 46–54.

Kaufman (1990)
Kaufman, A. A., jan 1990. The electrical field in a borehole with a casing.
Geophysics 55 (1), 29–38.
URL http://link.aip.org/link/?GPY/55/29/1{&}Agg=doihttp://library.seg.org/doi/abs/10.1190/1.1442769  Kaufman and Wightman (1993) Kaufman, A. A., Wightman, E. W., 1993. A transmissionline model for electrical logging through casing. Geophysics 58 (12), 1739.

Kohnke et al. (2017)
Kohnke, C., Liu, L., Streich, R., Swidinsky, A., 2017. A method of moments
approach to modeling the electromagnetic response of multiple steel casings
in a layered earth. Geophysics, 1–72.
URL https://library.seg.org/doi/10.1190/geo20170303.1 
Kong et al. (2009)
Kong, F. N., Roth, F., Olsen, P. A., Stalheim, S. O., 2009. Casing effects in
the seatoborehole electromagnetic method. Geophysics 74 (5), F77–F87.
URL http://library.seg.org/doi/10.1190/1.3173807 
Moran and Kunz (1962)
Moran, J. H., Kunz, K. S., dec 1962. Basic theory of induction logging and
application to study of twocoil sondes. GEOPHYSICS 27 (6), 829–858.
URL http://library.seg.org/doi/10.1190/1.1439108 
Nekut (1995)
Nekut, A. G., 1995. Crosswell electromagnetic tomography in steelâcased
wells. Geophysics 60 (3), 912–920.
URL http://library.seg.org/doi/10.1190/1.1443826 
Oldenburg et al. (2017)
Oldenburg, D. W., Heagy, L. J., Kang, S., 2017. Geophysical Electromagnetics:
Fundamentals and Applications. SEG Distinguished Instructor Short Course.
URL https://seg.org/Education/Courses/DISC/2017DISCDoug%****␣casingSoftware_processed.bbl␣Line␣200␣****Oldenburg  Pardo and TorresVerdin (2013) Pardo, D., TorresVerdin, C., 2013. Sensitivity analysis for the appraisal of hydrofractures in horizontal wells with borehole resistivity measurements. Geophysics 78 (4), D209–D222.
 Perez et al. (2015) Perez, F., Berkeley, L., Granger, B. E., 2015. Project Jupyter : Computational Narratives as the Engine of Collaborative Data Science (April), 1–24.
 Petra et al. (2014) Petra, C. G., Schenk, O., Anitescu, M., 2014. Realtime stochastic optimization of complex energy systems on highperformance computers. Computing in Science and Engineering 16 (5), 32–42.
 Puzyrev et al. (2017) Puzyrev, V., Vilamajo, E., Queralt, P., Ledo, J., Marcuello, A., 2017. ThreeDimensional Modeling of the Casing Effect in Onshore ControlledSource Electromagnetic Surveys. Surveys in Geophysics 38 (2), 527–545.

Ramirez et al. (1996)
Ramirez, A., Daily, W., Binley, A., LaBrecque, D., Roelant, D., 1996.
Detection of Leaks in Underground Storage Tanks Using Electrical Resistance
Methods. Journal of Environmental and Engineering Geophysics 1 (3),
189–203.
URL http://library.seg.org/doi/abs/10.4133/JEEG1.3.189 
Rocroi and Koulikov (1985)
Rocroi, J. P., Koulikov, A. V., feb 1985. The use of vertical line sources in
electrical prospecting for hydrocarbon. Geophysical Prospecting 33 (1),
138–152.
URL http://doi.wiley.com/10.1111/j.13652478.1985.tb00426.x 
Ronczka et al. (2015)
Ronczka, M., Rücker, C., Günther, T., 2015. Numerical study of
longelectrode electric resistivity tomography â Accuracy, sensitivity, and
resolution. Geophysics 80 (6), E317–E328.
URL http://library.seg.org/doi/10.1190/geo20140551.1  Rucker (2012) Rucker, D. F., 2012. Enhanced resolution for long electrode ERT. Geophysical Journal International 191 (1), 101–111.
 Rucker et al. (2010) Rucker, D. F., Loke, M. H., Levitt, M. T., Noonan, G. E., 2010. Electricalresistivity characterization of an industrial site using long electrodes. Geophysics 75 (4), WA95.

Schenkel and Morrison (1990)
Schenkel, C. J., Morrison, H. F., aug 1990. Effects of well casing on
potential field measurements using downhole current sources. Geophysical
Prospecting 38 (6), 663–686.
URL http://onlinelibrary.wiley.com/doi/10.1111/j.13652478.1990.tb01868.x/abstracthttp://doi.wiley.com/10.1111/j.13652478.1990.tb01868.x  Schenkel and Morrison (1994) Schenkel, C. J., Morrison, H. F., 1994. Electrical resistivity measurement through metal casing. Geophysics 59 (7), 1072–1082.
 Schwarzbach and Haber (2018) Schwarzbach, C., Haber, E., 2018. Improved upscaling of steelcased wells through inversion. SEG Technical Program Expanded Abstracts 2018, 888–892.
 Swidinsky et al. (2013) Swidinsky, A., Edwards, R. N., Jegen, M., 2013. The marine controlled source electromagnetic response of a steel borehole casing: Applications for the NEPTUNE Canada gas hydrate observatory. Geophysical Prospecting 61 (4), 842–856.
 Tietze et al. (2015) Tietze, K., Ritter, O., Veeken, P., 2015. Controlledsource electromagnetic monitoring of reservoir oil saturation using a novel boreholetosurface configuration. Geophysical Prospecting 63 (6), 1468–1490.

Uchida et al. (1991)
Uchida, T., Lee, K. H., Wilt, M. J., 1991. Effect of a steel casing on
crosshole EM measurement. SEG Technical Program Expanded Abstracts 1991,
442–445.
URL http://library.seg.org/doi/abs/10.1190/1.1889151 
Um et al. (2015)
Um, E. S., Commer, M., Newman, G. A., Hoversten, G. M., jun 2015. Finite
element modelling of transient electromagnetic fields near steelcased
wells. Geophysical Journal International 202 (2), 901–913.
URL http://gji.oxfordjournals.org/cgi/doi/10.1093/gji/ggv193  Van Rossum and Drake Jr (1995) Van Rossum, G., Drake Jr, F. L., 1995. Python reference manual. Centrum voor Wiskunde en Informatica Amsterdam.

Wait and Hill (1977)
Wait, J., Hill, D., apr 1977. Electromagnetic Shielding of Sources within a
MetalCased Bore Hole. IEEE Transactions on Geoscience Electronics 15 (2),
108–112.
URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?%****␣casingSoftware_processed.bbl␣Line␣300␣****arnumber=4071844  Wait (1983) Wait, J. R., 1983. Resistivity Response of a Homogeneous Earth with a Contained Vertical Conductor (1).

Wait and Williams (1985)
Wait, J. R., Williams, J. T., aug 1985. EM AND IP RESPONSE OF A STEEL WELL
CASING FOR A FOURâELECTRODE SURFACE ARRAY. PART I: THEORY. Geophysical
Prospecting 33 (5), 723–735.
URL http://doi.wiley.com/10.1111/j.13652478.1985.tb00775.x 
Ward and Hohmann (1988)
Ward, S. H., Hohmann, G. W., 1988. Electromagnetic Theory for Geophysical
Applications. In: Electromagnetic Methods in Applied Geophysics, 1st
Edition. Society of Exploration Geophysicists, Ch. 4, pp. 130–311.
URL http://library.seg.org/doi/abs/10.1190/1.9781560802631.ch4 
Weiss (2017)
Weiss, C. J., 2017. Finiteelement analysis for model parameters distributed
on a hierarchy of geometric simplices. Geophysics 82 (4), E155–E167.
URL http://library.seg.org/doi/10.1190/geo20170058.1 
Weiss et al. (2016)
Weiss, C. J., Aldridge, D. F., Knox, H. A., Schramm, K. A., Bartel, L. C.,
2016. The directcurrent response of electrically conducting fractures
excited by a grounded current source. Geophysics 81 (3), E201–E210.
URL http://library.seg.org/doi/10.1190/geo20150262.1 
Williams and Wait (1985)
Williams, J. T., Wait, J. R., aug 1985. EM and IP response of a steel well
casing for a fourelectrode surface array. Part II: Numerical results.
Geophysical Prospecting 33 (5), 736–745.
URL http://doi.wiley.com/10.1111/j.13652478.1985.tb00776.x 
Wilt et al. (1996)
Wilt, M., Lee, K. H., Becker, A., Spies, B., Wang, S., jan 1996. Crosshole EM
in steelâcased boreholes. In: SEG Technical Program Expanded Abstracts
1996. Society of Exploration Geophysicists, pp. 230–233.
URL http://library.seg.org/doi/abs/10.1190/1.1826605  Wilt et al. (2018) Wilt, M., Um, E., Weiss, C., Vasco, D., Petrov, P., Newman, G., Wu, Y., 2018. Wellbore Integrity Assessment with CasingBased Advanced Sensing (Figure 1), 1–7.

Wu and Habashy (1994)
Wu, X., Habashy, T. M., mar 1994. Influence of steel casings on
electromagnetic signals. Geophysics 59 (3), 378–390.
URL http://library.seg.org/doi/pdf/10.1190/1.1443600http://library.seg.org/doi/abs/10.1190/1.1443600  Yang et al. (2014) Yang, D., Oldenburg, D. W., Haber, E., 2014. 3D inversion of airborne electromagnetic data parallelized and accelerated by local mesh and adaptive soundings. Geophysical Journal International 196 (3), 1492–1507.

Yang et al. (2016)
Yang, D., Oldenburg, D. W., Heagy, L. J., sep 2016. 3D DC resistivity modeling
of steel casing for reservoir monitoring using equivalent resistor network.
SEG Technical Program Expanded Abstracts 2016, 932–936.
URL http://library.seg.org/doi/10.1190/segam201613868475.1 
Yee (1966)
Yee, K., may 1966. Numerical solution of initial boundary value problems
involving maxwell’s equations in isotropic media. IEEE Transactions on
Antennas and Propagation 14 (3), 302–307.
URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1138693{%}5Cnhttp://ieeexplore.ieee.org/xpls/abs{_}all.jsp?arnumber=1138693http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1138693  Zhang et al. (2018) Zhang, P., Brick, Y., Sharma, M. M., 2018. Numerical study of an electrodebased resistivity tool for fracture diagnostics in steelcased wellbores 83 (2).