Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Electromagnetic methods for imaging subsurface injections Heagy, Lindsey J. 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2019_february_heagy_lindsey.pdf [ 77.27MB ]
Metadata
JSON: 24-1.0374170.json
JSON-LD: 24-1.0374170-ld.json
RDF/XML (Pretty): 24-1.0374170-rdf.xml
RDF/JSON: 24-1.0374170-rdf.json
Turtle: 24-1.0374170-turtle.txt
N-Triples: 24-1.0374170-rdf-ntriples.txt
Original Record: 24-1.0374170-source.json
Full Text
24-1.0374170-fulltext.txt
Citation
24-1.0374170.ris

Full Text

Electromagnetic methods for imaging subsurfaceinjectionsbyLindsey J. HeagyB.Sc. Geophysics, University of Alberta, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Geophysics)The University of British Columbia(Vancouver)November 2018c© Lindsey J. Heagy, 2018The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, the thesis entitled:Electromagnetic methods for imaging subsurface injectionssubmitted by Lindsey J. Heagy in partial fulfillment of the requirements for the degreeof Doctor of Philosophy in Geophysics.Examining Committee:Douglas W. Oldenburg, Earth, Ocean and Atmospheric SciencesSupervisorEldad Haber, Earth, Ocean and Atmospheric SciencesSupervisory Committee MemberMichael Bostock, Earth, Ocean and Atmospheric SciencesSupervisory Committee MemberRoger Beckie, Earth, Ocean and Atmospheric SciencesUniversity ExaminerIan Mitchell, Computing ScienceUniversity ExamineriiAbstractThe extraction of hydrocarbons from low-permeability formations is commonly achievedthrough hydraulic fracturing. In a hydraulic fracturing operation, fluid and particles,called proppant (commonly sand) are injected to create fracture pathways and to keepthose pathways open so that hydrocarbons can flow from the reservoir to the wellbore.One of the key unknowns in hydraulic fracturing operations is the distribution and extentof proppant within the reservoir. If the electrical conductivity of the injected materialsis distinct from the host rock, then electromagnetic geophysical methods can be used.In order for electromagnetics to be a viable imaging technique for this application, wemust be able to: (a) collect data that are sensitive to the injected materials, and (b) havea method for estimating a representative model of the injected materials from those datathrough an inversion process. Numerical modelling is an essential tool for assessingfeasibility of electromagnetics and for developing a suitable inversion procedure for ex-tracting meaningful information from the collected data.A complicating factor of using electromagnetics in reservoir settings is that steel-cased wells are commonly present. Steel has vastly different electrical and magneticproperties than the surrounding rock and therefore significantly alters the behavior ofelectromagnetic fields and fluxes. The success of electromagnetic methods for imagingsubsurface injections, therefore, heavily relies on our ability to understand and simulateiiithe physical behavior of fields and fluxes in these settings.Using hydraulic fracturing as a motivating application, this thesis examines aspectsof both the imaging problem for subsurface injections as well as the fundamental be-havior of electromagnetic fields and fluxes in settings with steel-cased wells. I presenta strategy for estimating the electrical conductivity of a fractured volume of rock andincorporate this into the inverse problem. I also develop a numerical approach for ac-curately simulating electromagnetic surveys in settings with steel-cased wells. Usingthis software, I examine aspects of the fundamental physics, including how the mag-netic properties of a pipe complicate the behavior of the fields and fluxes, and how thisimpacts measured data.All of the software developed during the course of this research is open-source.ivLay SummaryHydraulic fracturing is used for producing hydrocarbons from low-permeability reser-voirs. Fluid and proppant (commonly sand) are injected to create fractures and to keepthose fracture pathways open. One of the unknowns in this process is the distributionof the injected materials. If the electrical conductivity of the injected materials is dis-tinct from the reservoir rock, electromagnetic geophysical data, which are sensitive tothe distribution of those materials, can be collected. Through an inversion process, a3D model of the injected materials can be estimated from those data. The presence ofsteel-cased wells complicates this procedure because steel has vastly different electricaland magnetic properties than the reservoir rock. As a result, it significantly alters thebehavior of the electromagnetic fields. This thesis examines aspects of the fundamentalphysics of electromagnetics in settings with steel-cased wells and as well as strategiesfor estimating the distribution of injected materials from electromagnetic data.vPrefaceThe research I conducted for this dissertation was completed at the Geophysical Inver-sion Facility at the University of British Columbia under the supervision of Dr. DouglasOldenburg. During my time as a PhD student, I also completed two internships, the firstwith the Schlumberger Electromagnetics Imaging group, under the advisement of Dr.Michael Wilt, and the second at Schlumberger Doll Research, under Dr. Dzevat Omer-agic. These internships informed my understanding of hydraulic fracturing and industrypractices for electromagnetic imaging. The research presented in this thesis has resultedin three peer-reviewed publications, seven conference proceedings and several auxiliaryworks.Chapter 2 presents an approach for constructing a physical property model of a frac-tured volume of rock. Dr. Wilt had the initial idea to use effective medium theoryto estimate the coarse-scale conductivity of a fractured volume of rock. I developed thetwo-step workflow, the algorithm used to solve for effective conductivity, and performedall numerical simulations in this chapter. Dr. Oldenburg advised on all of these devel-opments. Earlier versions of this work were presented at three conferences (Heagy andOldenburg, 2013; Heagy et al., 2014a; Wilt et al., 2014) and included in the patent Wiltet al. (2015).Chapter 3 discusses a finite volume approach for simulating Maxwells equations onvicylindrical meshes. This work has been submitted to a peer-reviewed journal and isavailable on arXiv (Heagy and Oldenburg, 2018b). Preliminary versions of this workwere included in two conference abstracts (Heagy et al., 2015, 2017b). I implementedthe software within the SimPEG framework and the code-review was conducted by Dr.Rowan Cockett. Dr. Oldenburg advised on the software development, testing, and con-tent of the chapter. I performed the numerical simulations that comprise the examples inthis chapter. Much of the numerical work in this chapter is supported by course materialand instruction from Dr. Eldad Haber and Dr. Uri Ascher (Haber, 2014; Ascher, 2008).Chapter 4 discusses aspects of the fundamental physics of DC resistivity with steel-cased wells and demonstrates survey design considerations. This work was submitted toa peer-reviewed journal and is available on arXiv (Heagy and Oldenburg, 2018a). I con-ducted all of the numerical experiments in this chapter with advisement from Dr. Old-enburg. The application of DC resistivity to the casing integrity problem was promptedby conversations with Dr. Wilt.Chapters 5 explores the physics of electromagnetics in settings with steel-casedwells. I conducted all of the experiments in this chapter with advisement from Dr.Oldenburg. Earlier versions of this work were included in (Heagy et al., 2015, 2017b).Chapter 6 discusses the inverse problem for imaging a fractured volume of rock. Ideveloped the idea of using effective medium theory in the inversion, implemented thenecessary software components and performed the numerical experiments in this chap-ter, all with advisement from Dr. Oldenburg. An early version of the idea of usingeffective medium theory to invert for fracture concentration was presented in a confer-ence publication (Heagy et al., 2014b). Dr. Cockett contributed to the development ofexamples in the conference publication.Appendix B and Appendix C contain details relevant to Chapters 2 and 6. I per-viiformed both derivations with advisement from Dr. Oldenburg. The derivation in B wasinspired by a question from Dr. Frank Morrison.Appendix D describes the SimPEG electromagnetics module which is the backbonefor the numerical experiments conducted in this thesis. This work was published inComputers & Geosciences (Heagy et al., 2017a) and was conducted collaborativelywith Dr. Cockett, Dr. Seogi Kang, and Mr. Gudni Rosenkjaer with advisement fromDr. Oldenburg. I led the development and implementation of the forward-simulationframework as well as the drafting of the manuscript. Dr. Kang performed the inversionof the Bookpurnong data as described in Section D.5.2. All authors contributed edits tothe text.Appendix E provides an overview of the GeoSci.xyz project (https://geosci.xyz)which is a collaborative effort to develop educational resources for the geosciences.This effort has been led by myself, Dr. Oldenburg, Dr. Kang, and Dr. Cockett, and hasbenefitted from a community of contributors. Notably, Mr. Dominique Fournier, Mr.Devin Cowan, Mr. Thibaut Astic, Dr. Sarah Devriese, Mr. Michael Mitchell, and Dr.Dianne Mitchenson, have invested significant efforts in contributing content.Software development has been a critical part of this work. In particular, the re-search conducted in this thesis builds upon and contributes to the open-source SimPEGecosystem, described in Cockett et al. (2015), of which I am a co-author. Excerpts ofthe background information on inversion in Section 1.6 are adapted from Cockett et al.(2015). My personal contribution to SimPEG exceeds 100,000 lines of code. In addi-tion, I help maintain SimPEG and lead community development efforts by responding tocode-usage questions and conducting code-reviews. The research in this thesis benefitsfrom the efforts of all of those who are involved in the SimPEG community. Notably,Dr. Cockett and Dr. Kang have made substantial contributions to the architecture andviiifunctionality of SimPEG that have been used to conduct the research in this thesis; theyhave also reviewed and provided input on much of the code that I have written.The research in this thesis makes use of tools in the scientific Python ecosystem in-cluding Jupyter, Numpy, SciPy, and Matplotlib (Perez et al., 2015; van der Walt et al.,2011; Jones et al., 2001; Hunter, 2007) and relies upon the global community of con-tributors who develop and maintain these open-source resources.ixTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . liiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . liv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Hydraulic Fracturing . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Challenges to electromagnetics with steel-cased boreholes . . . . . . . 71.4 Research motivation and thesis structure . . . . . . . . . . . . . . . . . 121.4.1 Homogenization strategy for hydraulic fractures . . . . . . . . . 13x1.4.2 Steel casings . . . . . . . . . . . . . . . . . . . . . . . . . . . 141.4.3 Inversions for subsurface injections . . . . . . . . . . . . . . . 171.4.4 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.4.5 A note on reproducibility and dissemination . . . . . . . . . . . 171.5 Background: Electromagnetic geophysics . . . . . . . . . . . . . . . . 181.5.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . 191.5.2 Geophysical surveys . . . . . . . . . . . . . . . . . . . . . . . 221.6 Background: Geophysical Inversions . . . . . . . . . . . . . . . . . . . 251.6.1 Formulating and solving the inverse problem . . . . . . . . . . 252 A physical property model for a fractured volume of rock . . . . . . . . . 322.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322.1.1 Proppant selection . . . . . . . . . . . . . . . . . . . . . . . . 332.1.2 Numerical Modelling . . . . . . . . . . . . . . . . . . . . . . . 342.1.3 Chapter overview . . . . . . . . . . . . . . . . . . . . . . . . . 362.2 Homogenization workflow using effective medium theory . . . . . . . . 372.2.1 Summary of self-consistent effective medium theory . . . . . . 382.2.2 Step 1: Effective conductivity of a proppant-fluid mixture . . . . 432.2.3 Step 2: Effective conductivity of fractured volume of rock . . . 482.2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 552.3 Feasibility of EM for detecting fractures . . . . . . . . . . . . . . . . . 552.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623 Modelling electromagnetics on cylindrical meshes . . . . . . . . . . . . . 643.1 Numerical tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 673.1.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 69xi3.1.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 743.2 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 773.2.1 DC Resistivity Part 1: Electric fields, currents and charges in along well . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 783.2.2 DC Resistivity Part 2: Finite Length Wells . . . . . . . . . . . . 813.2.3 Frequency Domain Electromagnetics Part 1: Comparison withscale model results . . . . . . . . . . . . . . . . . . . . . . . . 843.2.4 Frequency Domain Electromagnetics Part 2: Conductivity andpermeability in the inductive response of a well . . . . . . . . . 903.3 Summary and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . 954 Direct current resistivity with steel-cased wells . . . . . . . . . . . . . . . 984.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 984.2 DC resistivity for casing integrity . . . . . . . . . . . . . . . . . . . . . 1014.2.1 Basic experiment . . . . . . . . . . . . . . . . . . . . . . . . . 1024.2.2 Survey design considerations . . . . . . . . . . . . . . . . . . . 1074.2.3 Factors influencing detectability . . . . . . . . . . . . . . . . . 1094.2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1194.3 Survey design for exciting targets at depth . . . . . . . . . . . . . . . . 1214.3.1 Source location . . . . . . . . . . . . . . . . . . . . . . . . . . 1224.3.2 Target properties . . . . . . . . . . . . . . . . . . . . . . . . . 1264.4 Coarse-scale approximations of the well . . . . . . . . . . . . . . . . . 1424.4.1 Replacing a hollow-cased well with a solid cylinder . . . . . . . 1434.4.2 Cartesian grid . . . . . . . . . . . . . . . . . . . . . . . . . . . 1474.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151xii5 Electromagnetics with steel cased wells . . . . . . . . . . . . . . . . . . . 1545.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1545.2 EM response of a conductive casing . . . . . . . . . . . . . . . . . . . 1575.2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1635.2.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1695.3 EM response of permeable casings . . . . . . . . . . . . . . . . . . . . 1705.3.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1735.3.2 Down-hole source . . . . . . . . . . . . . . . . . . . . . . . . . 1765.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1775.4 Approximations to a steel cased well . . . . . . . . . . . . . . . . . . . 1805.4.1 Electrical conductivity . . . . . . . . . . . . . . . . . . . . . . 1815.4.2 Magnetic permeability . . . . . . . . . . . . . . . . . . . . . . 1875.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1926 An inversion approach for subsurface injections . . . . . . . . . . . . . . 1946.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1946.2 Choosing an inversion model . . . . . . . . . . . . . . . . . . . . . . . 1976.3 Inversions with steel-cased wells . . . . . . . . . . . . . . . . . . . . . 2046.3.1 Voxel inversion for conductivity . . . . . . . . . . . . . . . . . 2066.3.2 Parametric Inversion . . . . . . . . . . . . . . . . . . . . . . . 2116.3.3 Inversion for fracture concentration . . . . . . . . . . . . . . . 2176.3.4 Adding a-priori volume information . . . . . . . . . . . . . . . 2246.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2257 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236xiiiA Source code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 249B Concentric spheres in a uniform electrostatic field . . . . . . . . . . . . . 251B.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251B.2 Solving for the Potential . . . . . . . . . . . . . . . . . . . . . . . . . 253B.2.1 Solving for the coefficients . . . . . . . . . . . . . . . . . . . . 254B.2.2 Sanity Checks . . . . . . . . . . . . . . . . . . . . . . . . . . . 266B.3 Effective Conductivity Expression . . . . . . . . . . . . . . . . . . . . 272C Self consistent effective medium theory derivatives . . . . . . . . . . . . . 275D A framework for simulation and inversion in electromagnetics . . . . . . 278D.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 278D.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282D.2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 285D.3 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 287D.4 Simulation Framework . . . . . . . . . . . . . . . . . . . . . . . . . . 288D.4.1 Model and Physical Properties . . . . . . . . . . . . . . . . . . 294D.4.2 Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 296D.4.3 Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 300D.4.4 Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302D.4.5 Receivers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303D.4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303D.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304D.5.1 Cylindrically Symmetric Inversions . . . . . . . . . . . . . . . 305D.5.2 Bookpurnong Field Example . . . . . . . . . . . . . . . . . . . 308D.5.3 Steel-Cased Well: Sensitivity Analysis for a Parametric Model . 315xivD.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 327E Open source practices for education . . . . . . . . . . . . . . . . . . . . . 329E.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 329E.2 GeoSci.xyz Textbooks . . . . . . . . . . . . . . . . . . . . . . . . . . . 332E.3 Interactive content . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334E.4 Open development . . . . . . . . . . . . . . . . . . . . . . . . . . . . 337E.5 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . . 339xvList of TablesTable 4.1 Integrated secondary charge over a target adjacent to the casing, asshown in Figure 4.14. . . . . . . . . . . . . . . . . . . . . . . . . . 129Table 4.2 Integrated secondary charge over a target that is not electrically con-nected to the casing, as shown in Figure 4.17. . . . . . . . . . . . . . 135Table D.1 Comparison of the maximum amplitude of the sensitivity with re-spect to each model parameter, and the approximate perturbation inthat parameter required to produce a 10−9 V/m change in the mea-sured data. The conversion from a perturbation in log-conductivityto conductivity is given by equation D.26. The perturbation in con-ductivity is also provided in terms of a percentage of the true modelconductivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 328xviList of FiguresFigure 1.1 Conventional reservoirs contain oil and gas that have migrated up-wards under pressure until they are trapped by a cap rock (left),while non-conventional tight or shale oil and gas reservoirs con-tain hydrocarbons that are trapped in low permeability formations(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 1.2 The choice of whether to solve the EM equations in the time-domain(TDEM) or the frequency domain (FDEM) depends upon the trans-mitter waveform; a step-off-type waveform (left) lends itself to asolution approach in the time domain and a harmonic signal (right)is typically addressed in the frequency domain. . . . . . . . . . . . 23Figure 1.3 Time-domain inductive source experiment. (Left) a steady-state cur-rent is passed through the transmitter loop, generating a primarymagnetic field. (Right) This current is shut-off, causing a change inmagnetic flux. The changing magnetic flux induces secondary cur-rents in conductors, which in turn create secondary magnetic fieldsthat can be measured at receivers above, on, or within the earth. . . 23xviiFigure 1.4 (Left) Direct current resistivity experiment in which a steady-statecurrent is injected into the earth, and (Right) a grounded source elec-tromagnetic experiment which uses a time-varying source current. . 24Figure 1.5 Overview of a geophysical inversion workflow. Adapted from Cock-ett et al. (2015). . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.1 Fracture complexity, from simple planar fractures on the left to com-plex fracture networks on the right. The blue dot is the injectionpoint in a horizontal well. After Cipolla et al. (2008a). . . . . . . . 35Figure 2.2 Constructing a physical property model for a fractured volume ofrock using effective medium theory. The electrical conductivity ofthe proppant fluid mixture is given by σ2, the conductivity of thebackground reservoir rock by σ1. Using effective medium theory,the coarse-scale anisotropic conductivity, Σ∗ describing the frac-tured volume of rock is computed. . . . . . . . . . . . . . . . . . . 39Figure 2.3 Effective conductivity of a proppant-fluid mixture for five differ-ent proppant conductivities, each indicated in the legend. Panel (a)shows the conductivity on a linear scale, and panel (b) uses a log-scale for the conductivity. Both panels use a linear scale for thevolume fraction of the proppant (ϕ). The conductivity of the fluid is3 S/m, similar to the conductivity of sea-water. . . . . . . . . . . . 44Figure 2.4 Impact of the conductivity of the fluid on the effective conductivityof a proppant-fluid mixture. Panels (a) and (c) show the effectiveconductivity for mixtures with a 103 S/m proppant and panels (b)and (d) show the effective conductivity for mixtures with a 104 S/mproppant. The conductivity of the fluid is indicated by the legend. . 45xviiiFigure 2.5 Effective conductivity of a 3-phase proppant-fluid mixture consist-ing of resistive proppant (10−6 S/m), conductive proppant (105 S/m)and saline fluid (3 S/m). The legend indicates the percentage ofconductive proppant in the proppant mixture and the x-axis is thevolume fraction of proppant in the proppant-fluid mixture. . . . . . 47Figure 2.6 Effective conductivity of a 3-phase proppant-fluid mixture consist-ing of resistive proppant (10−6 S/m), conductive proppant (105 S/m)and saline fluid (3 S/m). The proppant mixture contains 25% con-ductive proppant and 75% resistive proppant. The legend indicatesthe aspect ratio of the elongated conductive proppant filler (prolatespheroids). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48Figure 2.7 Oblate spheroid with normal along the x-axis. The y and z semi-axesare equal, and aspect ratio is α = b/a< 1. . . . . . . . . . . . . . 50xixFigure 2.8 Effective, anisotropic conductivity for a fractured rock with spheroidalcracks whose normal is oriented along the x-axis for five differ-ent aspect ratios, indicated by the legend. The black dashed linesshow the upper and lower Wiener Bounds, which are identical to thevolume-weighted arithmetic and harmonic averages of the conduc-tivity of the rock (0.1 S/m) and proppant-fluid mixture (2500 S/m).The black dash-dot lines are the anisotropic Hashin-Shtrikman up-per and lower bounds computed using an aspect ratio of 10−5 inequation 2.12. Note that the upper-bound coencides with σy,z forsmall aspect ratios (e.g. the blue line). Panels (a) and (c) showthe full range 0≤ ϕ ≤ 1, and panels (b) and (d) zoom in to 0≤ ϕ ≤0.005. Note that in (b), any curve that deviates from the lower boundis the σy,z component; all σx values lie on the lower bound for theconcentration range shown. . . . . . . . . . . . . . . . . . . . . . 52Figure 2.9 Effective, isotropic conductivity for a fractured rock with randomlyoriented spheroidal cracks for five different aspect ratios, indicatedby the legend. The black dashed lines show the upper and lowerWiener Bounds, which are identical the volume-weighted arithmeticand harmonic averages of the conductivity of the rock (0.1 S/m)and proppant-fluid mixture (2500 S/m).The black dash-dot lines arethe isotropic Hashin-Shtrikman upper and lower bounds. The semi-transparent dotted lines show σy,z from Figure 2.8 Panels (a) and (c)shows the full range 0 ≤ ϕ ≤ 1, and panels (b) and (d) zooms in to0≤ ϕ ≤ 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54xxFigure 2.10 Setup for the crosswell electromagnetic simulation. Ten fractures,each with 3 mm thickness and 50 m radius are spaced over a 10 minterval along the injection well. The dipole transmitter is positionedinside of the injection well and oriented along the x-axis. Receiversare located in an offset well and measure the real and imaginaryparts of the x-component of the magnetic field. . . . . . . . . . . . 57Figure 2.11 Simulated secondary magnetic field data for: (a) the anisotropicfracture model and (b) the isotropic fracture model, in a crosswellsurvey conducted at 100Hz. The top row shows the real componentand the bottom row shows the imaginary component of the mea-sured magnetic flux (nT). In all plots, the x-axis shows the trans-mitter location relative to the center of the fracture, the y-axis showsthe receiver location, and the color indicates the secondary magneticflux with respect to a 0.1 S/m whole-space. Values beneath the noisefloor of 10−4 nT have been masked and display as white. To displaythe signal as a percentage of the primary, we have included a semi-transparent overlay between 0% and 5%; low percentage values plotas darker greys. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Figure 2.12 Simulated secondary magnetic field data for: (a) the anisotropicfracture model and (b) the isotropic fracture model, in a crosswellsurvey conducted at 200Hz. The top row shows the real componentand the bottom row shows the imaginary component of the measuredmagnetic flux (nT), similar to Figure 2.11 . . . . . . . . . . . . . . 61xxiFigure 3.1 Anatomy of a finite volume cell in a (a) cartesian, regtangular mesh,(b) cylindrically symmetric mesh, and (c) a three dimensional cylin-drical mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 3.2 Construction of a 3D cylindrical mesh from a cartesian mesh. . . . . 73Figure 3.3 Depth slice (left) and cross section (right) through the 3D cylindri-cal mesh used for the comparison with Commer et al. (2015). Thesource and recievers are positioned along the θ = 90◦ line. Themesh extends 17km raidally and 30km vertically to ensure that thefields have sufficiently decayed before reaching the boundaries. . . 76Figure 3.4 Time domain EM response comparison with (Commer et al., 2015).Each of the different line colors is associated with a different loca-tion; offsets are with respect to the location of the well. . . . . . . . 77Figure 3.5 (a) Total charge density, (b) secondary charge density, (c) electricfield, and (d) current density in a section of the pipe near the sourceat z=-500 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 3.6 (a) Current along a well for 5 different wellbore lengths. The x-axisis depth normalized by the length of the well. The black dashed lineshows the short-well approximation (equation 45 in Kaufman andWightman (1993)) for a 200 m long well. The black dash-dot lineshows the long-well approximation (equation 53 in Kaufman andWightman (1993)) for a 4000 m well. (b) Charge per unit lengthalong the well for 5 different wellbore lengths. . . . . . . . . . . . 82xxiiFigure 3.7 Field strength ratio (FSR), the ratio of the measured vertical mag-netic field with the free space magnetic field, as a function of fre-quency for two different receiver locations. In (a), the receiver isin the same plane as the source, in (b), the receiver is 1.49 m offsetfrom the source. . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 3.8 Magnetic flux density at 0.1 Hz in the region of the pipe near theplane of the source for (a) the “infinite” pipe, where the source islocated at -4.5 m and the pipe extends from 0 m to -9 m, (b) a “semi-infinite” pipe, where the source is located at 0 m and the pipe extendsto -9 m. In (c), we zoom in to the top 5 cm of the “semi-infinite”pipe, and (d) shows the 5 cm at the top-end of the “infinite” pipe. . 89Figure 3.9 Field strength ratio, FSR, for a receiver positioned 3 cm beneaththe plane of the source. For comparison, we have plotted the FSRfor the permeable pipe when the source and receiver lie in the sameplane (L=0.00 m) with the semi-transparent orange lines. Note thatthe infinite-pipe solutions for L=0.03 m and L=0.00 m overlap. . . 90Figure 3.10 Normalized secondary field, NSF, as a function of frequency for twowells. The NSF is the ratio of the secondary vertical magnetic fieldwith the primary magnetic field at the receiver location (z=-500 m);the primary is defined as the whole-space primary. . . . . . . . . . 93Figure 3.11 Secondary magnetic flux density (with respect to a whole-space pri-mary) at five different frequencies for a conductive pipe (top row)and for a conductive, permeable pipe (bottom row). . . . . . . . . . 94xxiiiFigure 3.12 Normalized secondary field (NSF) through time. In the time-domain,we compute the NSF by taking the difference between the total mag-netic flux at the receiver and the whole-space response and then tak-ing the ratio with the whole-space magnetic flux prior to shutting offthe transmitter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Figure 3.13 Secondary magnetic flux density for a conductive well (top row)and a conductive, permeable well (bottom row) through time. Thesource waveform is a step-off waveform. . . . . . . . . . . . . . . 96Figure 4.1 Cross section showing: (a) electrical conductivity, (b) current den-sity, (c) charge density, and (d) electric field for a top-casing DCresistivity experiment over (top) an intact 1000 m long well and(bottom) a 1000 m long well with a 10 m flaw at 500 m depth. . . . 103Figure 4.2 (a) Charge along the length of the intact well (black), a 500 m well (“short”, grey dash-dot), and a well with a 10 m flaw at 500 m depth(blue), in a top-casing DC resistivity experiment. (b) Secondaryelectric field due on the surface of the earth due to the flaw in thecasing. The primary is defined as the electric field due to the 1000m long intact well. The return electrode is 2000 m away from thewell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105xxivFigure 4.3 Charge along the length of a 50 m long intact well (black), a 25m well (“short”, grey dash-dot), and four wells, each with a flawstarting at 25 m depth and extending the length indicated by thelegend (5× 10−1 m (blue), 5× 10−2 m (orange), and 5× 10−3 m(green)) in a top-casing DC resistivity experiment. For reference,the diameter of the casing is 10−1 m and its thickness is 10−2 m.The return electrode is 50 m away from the well and a cylindricallysymmetric mesh was used in the simulation. . . . . . . . . . . . . 106Figure 4.4 (Top row) primary electric field, (second row) secondary electricfield, and (third row) secondary electric field as a percentage of theprimary radial electric field for a return electrode that is offset (a)2000 m, (b) 750 m, (c) 500 m, and (d) 250 m from the well. Theprimary is defined as the response due to the 1000 m long, intactwell. Electrode locations are denoted by the red dots. In the thirdrow, the colorbar has been limited between 20% and 100%. Thefourth and fifth rows show radial electric field data collected alongthe θ = 90◦ azimuth (white dotted lines in the top three rows). Thefourth row shows the primary (black line), the total electric field dueto the flawed well (blue line), and the secondary radial electric field(orange line). The fifth row shows the secondary as a percentage ofthe primary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108xxvFigure 4.5 Radial electric field as the depth of the flaw along a 1km long well isvaried. The positive electrode is connected to the top of the casing,the negative electrode is positioned 500m away and data are mea-sured along a line 90◦ from the source electrodes. In (a), we showthe total electric field for four flawed wells, each with a 10m flaw atthe depth indicated on the legend. The black line shows the radialelectric field due to an intact well; we define this as the primary. In(b), the secondary radial electric field is plotted and in (c), we showthe secondary radial electric field as a percentage of the primary. . . 111Figure 4.6 Radial electric field as the conductivity of the background is variedfor a 1 km well with a 10 m flaw at 500 m depth. The positive elec-trode is connected to the top of the casing, the negative electrode ispositioned 500 m away and data are measured along a line 90◦ fromthe source electrodes. In (a), we show the total electric field forfive different background conductivities, each indicated on the leg-end. The solid lines indicate the response of the flawed well and thedashed lines indicate the response of the intact well (the primary). In(b), the secondary radial electric field is plotted and in (c), we showthe secondary radial electric field as a percentage of the primary. . . 113xxviFigure 4.7 Radial electric field as the conductivity of a 50 m thick layer posi-tioned at 400 m depth is varied. The positive electrode is connectedto the top of the casing, the negative electrode is positioned 500 maway and data are measured along a line 90◦ from the source elec-trodes. In (a), we show the total electric field five different layerconductivities. The black line shows the scenario where the layerhas the same conductivity as the background. The dashed-lines in-dicate the intact well and the solid lines indicate the flawed well. In(b), the secondary radial electric field is plotted (with respect to anintact well primary) and in (c), we show the secondary radial electricfield as a percentage of the primary. . . . . . . . . . . . . . . . . . 114Figure 4.8 Cross section showing: (a) electrical conductivity, (b) current den-sity, (c) charge density, and (d) electric field for a top-casing DCresistivity experiment over models with a conductive layer (top tworows) and a model with a resistive layer (bottom two rows). Thelayer extends from 400 m to 450 m depth. The plots in the secondand fourth rows correspond to models with a 10 m flaw at 500 mdepth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116Figure 4.9 (a) Charge along the length of wells with three different conductiv-ities (each indicated by a different color in the legend). The intactwells are denoted with dashed lines and the flawed wells are de-noted with solid lines. (b) Secondary charge along the flawed andshort wells. The primary is defined as the electric field due to the1000 m long intact well. The return electrode is 2000 m away fromthe well. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117xxviiFigure 4.10 Radial electric field as the conductivity of the casing is varied for a1 km well with a 10 m flaw at 500 m depth. The positive electrodeis connected to the top of the casing, the negative electrode is posi-tioned 500 m away and data are measured along a line 90◦ from thesource electrodes. In (a), we show the total electric field for threedifferent casing conductivities, each indicated on the legend. Thesolid lines indicate the response of the flawed well and the dashedlines indicate the response of the intact well (the primary). In (b),the secondary radial electric field is plotted and in (c), we show thesecondary radial electric field as a percentage of the primary. . . . . 118Figure 4.11 Radial electric field as the vertical extent of the flaw is varied. Thepositive electrode is connected to the top of the casing, the negativeelectrode is positioned 500 m away and data are measured alonga line 90◦ from the source electrodes. In (a), we show the totalelectric field four different flaw extents. The black line shows theresponse of the intact well. The dashed-lines indicate the partiallyflawed wells (50% of the circumference is compromised) and thesolid lines flawed wells in which the entire circumference of the wellhas been compromised. In (b), the secondary radial electric field isplotted (with respect to an intact well primary) and in (c), we showthe secondary radial electric field as a percentage of the primary. . . 120xxviiiFigure 4.12 Electrode locations to be compared. The top casing electrode (blue),centered electrode (green, 500 m depth), and downhole electrode(red, 500 m depth) are connected to the casing. The surface elec-trode (orange) is offset from the well by 0.1 m. The remaining elec-trodes are positioned along the axis of the casing. Panel (a) showsthe entire length of the casing, while (b) zooms in to the bottom ofthe casing to show the separation between the electrodes beneath thecasing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123Figure 4.13 Total current density along a vertical line offset (a) 25 m, (b) 50 mand (c) 100 m from the axis of the casing, which extends from thesurface (0 m) to 1000 m depth. The electrode locations correspondto those shown in Figure 4.12. For reference, a simulation with anelectrode 20 m below the casing when there is no casing present isshown in black. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125Figure 4.14 Cross section showing: (a) electrical conductivity, (b) current den-sity, (c) charge density, and (d) electric field for a DC resistivityexperiment with a resistive target (top) and a conductive target (bot-tom). The positive electrode is positioned in the casing at the 912.5m depth. The casing is shown by the black line that extends to 1 kmdepth in panel (a). . . . . . . . . . . . . . . . . . . . . . . . . . . 129xxixFigure 4.15 Cross section showing: (a) electrical conductivity, (b) current den-sity, (c) charge density, and (d) electric field for a DC resistivityexperiment with a conductive target (top) and a resistive target (bot-tom). The positive electrode is positioned at 912.5 m depth. Nocasing is included in this simulation. Note that the colorbars for thecharge density (c) and electric field (d) are different than those usedin Figure 4.14. For the resistive target, the colorbar is saturated, thecharge density over the resistive target is on the order of 10−13 C/m3. 130Figure 4.16 Radial electric field at the surface as the conductivity of a cylindri-cal target, in contact with the well, is varied. The target has a radiusof 25 m and extends in depth from 900 m to 925 m. The returnelectrode is on the surface, 500 m from the well and data are mea-sured along a line perpendicular to the source. The panels on theleft show (a) the total electric field, (b) the secondary electric fieldwith respect to a primary that does not include the target, and (c) thesecondary electric field as a percentage of the primary for a surveyin which the positive electrode is positioned downhole at 912.5 mdepth. The panels on the right similarly show (d) the total electricfield, (e) the secondary electric field, and (f) the secondary electricfield as a percentage of the primary for a top-casing experiment. . . 131xxxFigure 4.17 Cross section showing: (a) electrical conductivity, (b) current den-sity, (c) charge density, and (d) electric field for a DC resistivityexperiment with a conductive target (top) and a resistive target (bot-tom) which is not in contact with the well. The positive electrode ispositioned in the casing at the 912.5 m depth. The casing is shownby the black line that extends to 1 km depth in panel (a). . . . . . . 134Figure 4.18 Radial electric field at the surface as the conductivity of a cylindricaltarget, which is not contact with the well, is varied. The target hasa radius of 25 m and extends in depth from 900 m to 925 m. Thereturn electrode is on the surface, 500 m from the well and data aremeasured along a line perpendicular to the source. The panels on theleft show (a) the total electric field, (b) the secondary electric fieldwith respect to a primary that does not include the target, and (c) thesecondary electric field as a percentage of the primary for a surveyin which the positive electrode is positioned downhole at 912.5 mdepth. The panels on the right similarly show (d) the total electricfield, (e) the secondary electric field, and (f) the secondary electricfield as a percentage of the primary for a top-casing experiment. Thedata shown in Figure 4.16, for the target in contact with the well, areplotted in the dashed, semi-transparent lines for reference. . . . . . 136xxxiFigure 4.19 Depth slice showing the primary electric field due to a downholeelectrode and a return electrode located on the surface at x =-500m, y =0 m. The red line indicates the azimuth of the source. Weexamine the 3 different target azimuths shown by the white outlines.The solid line indicates the target inline with the source, the dashedis 90◦ from the source line, and the dotted line is 180◦ from thesource line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138Figure 4.20 Integrated anomalous current density (excitation), as defined in equa-tion 4.2, for a 50 m × 50 m × 25 m target at 900 m depth in aDC experiment with the positive electrode (a) downhole at 912.5 mdepth, and (b) a top-casing electrode. The return electrode is posi-tioned on the surface 500 m from the well. Each line color indicatesa different target-conductivity. The different line styles correspondto different target azimuths relative to the plane of the source elec-trodes and correspond to those show in Figure 4.19. The solid lineindicates a target inline with the source, the dashed is 90◦ from thesource, and the dotted line is 180◦ from the source. Offset is calcu-lated from the center of the well to the edge of the target closest tothe well. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139xxxiiFigure 4.21 (a) Sum of the primary and secondary radial electric field due to adipolar target with moment of 38 Am centered 50 m from the well,either calculated with the casing (blue) or simply a dipole in a half-space (orange). (b) Secondary radial electric field due to a dipolartarget in a halfspace with casing (blue) and without casing (orange).Secondary radial electric field as a percentage of the primary. Thetarget is along a line 90◦ from the source electrodes; this is the sameline along which we measure data at the surface. . . . . . . . . . . 141Figure 4.22 Currents (top row) and charges (bottom row) along the length of ahollow steel-cased well (solid lines), solid cylinder with conductiv-ity equal to that of the steel-cased well (dashed-lines), and a solidcylinder with a conductivity such that the product of the conductiv-ity and the cross sectional area of the cylinder is equal to that of thehollow-pipe (dotted lines). Each of the line-colors corresponds to adifferent casing length, as indicated in the legend. In (a), we showthe vertical current in the casing, (b) shows the difference from thetrue, hollow-cased well in the vertical current within the casing, and(c) shows that difference as a percentage of the true currents. In (d),we show the charge per unit length along the casing, (e) shows thedifference from the true, hollow-cased well and (f) shows that dif-ferences as a percentage of the true charge distribution. The x-axison all plots is depth normalized by the length of the casing. . . . . 144xxxiiiFigure 4.23 Radial electric field measured at the surface for a model of a hol-low steel-cased well (solid lines), a solid cylinder with conductivityequal to that of the steel-cased well (dashed-lines), and a solid cylin-der with a conductivity such that the product of the conductivity andthe cross sectional area of the cylinder is equal to that of the hollow-pipe (dotted lines). Each of the line-colors corresponds to a differentcasing length, as indicated in the legend. In (a), we show the totalradial electric field, (b) shows the difference in electric field fromthat due to the true, hollow-cased well, and (c) shows that differenceas a percentage of the true electric fields. The x-axis on all plots isdistance from the well normalized by the length of the casing. . . . 146Figure 4.24 Three realizations of a 2 km long casing in a layered background,where the conductivity of the layers is assigned randomly. Eachlayer is 50 m thick, and the mean conductivity of the backgroundis 0.1 S/m. The color of the title corresponds to the plots of thecurrents and charges in Figure 4.25 . . . . . . . . . . . . . . . . . 147Figure 4.25 (a) Total vertical current through the casing for the three layered-earth models shown in Figure 4.24. The solid lines indicate the re-sponse of the true, hollow steel cased-well and the dotted lines indi-cate the response of a solid cylinder having the same cross-sectionalconductance as the hollow well. (b) Difference between the currentsalong the casing in the solid well approximation and the true, hollowwell. (c) Charge per unit length for each of the models. (d) Differ-ence in charge per unit length between the true model of the casingand the approximation which preserves cross-sectional conductance. 148xxxivFigure 4.26 Currents (top row) and charges (bottom row) along the length of asteel cased well. The “true” hollow-cased well is simulated on a 3Dcylindrical mesh and has 4 cells across the width of the casing thick-ness (black line). The colored lines correspond to the currents andcharges computed along the well represented on a cartesian meshwith cell widths shown in the legend. The finest vertical discretiza-tion is 2.5 m in all simulations. To represent the hollow cased wellon the cartesian mesh, the cells intersected by the casing are as-signed a conductivity that preserves the product of the conductivityand cross-sectional area of the well. In (a), we show the verticalcurrent in the casing, (b) shows the difference from the true, hollow-cased well in the vertical current within the casing, and (c) showsthat difference as a percentage of the true currents. In (d), we showthe charge per unit length along the casing, (e) shows the differencefrom the true, hollow-cased well and (e) shows that differences as apercentage of the true charge distribution. . . . . . . . . . . . . . . 149Figure 5.1 Current density for a time domain experiment over a 10−2 S/m half-space. The positive electrode is on the surface where the well-headwill be and the return electrode is at x = 1000 m. Each row corre-sponds to different time, as indicated in the plots in panel (c). Panel(a) shows a cross section through the half-space along the same lineas the source-wire. Panels (b) and (c) show depth-slices of the cur-rents at 54 m and 801 m depth. . . . . . . . . . . . . . . . . . . . 158xxxvFigure 5.2 Current density for a top-casing time domain EM experiment witha conductive well (5×106 S/m). (a) Cross section in the immediatevicinity of the well. (b) Cross section through the formation. (c – d)depth-slices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160Figure 5.3 Difference in current density for a time domain experiment whichincludes a conductive steel-cased well (as in Figure 5.2) and an ex-periment over a halfspace (as in Figure 5.1). . . . . . . . . . . . . 162Figure 5.4 Simulated electric field at the surface of the earth at (a) 0.01 ms and(b) 5 ms after shut-off for the halfspace model. (c) Radial electricfield data measured along the survey line shown in white in (a) and(b) at 0.01 ms (blue) and 5 ms (green). (d) Radial electric field as afunction of time at 300 m along the survey line (shown in the red dotin (a)). The red dots in (c) correspond to the data observed at 300m offset and the blue and green dots in (d) correspond to the 0.01ms and 5 ms data. Similar information is shown in (e), (f), (g) and(h) for the model with the conductive casing. The difference in theradial electric field data (casing minus halfspace) is shown in (i), (j),(k) and (l). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164xxxviFigure 5.5 Sketch of the plan-view geometry of the early-time galvanic cur-rents (left) and image currents (right) in a time-domain EM experi-ment over a half-space. Through time, both current systems diffusedownwards and outwards as they decay. The wire follows a straightline between the negative and positive electrodes. The green dashedline shows where we are measuring the radial electric field data.Prior to shut-off, current in the wire flows from the negative to thepositive electrode. In the earth, the galvanic currents are dipolar innature and flow from the positive to the negative electrode. Alongthe survey-line, the radial component of the galvanic currents al-ways points outwards. Immediately after shut-off, image currentsare induced in the earth. They are oriented in the same direction asthe original current in the wire and are directed away from the nega-tive electrode towards the positive. Along the survey-line, the radialcomponent of the image currents is always pointed inwards. . . . . 165Figure 5.6 Simulated db/dt at (a) 0.01 ms and (b) 5 ms after shut-off for thehalfspace model. (c) Tangential db/dt measured along the whiteline in (a) and (b) at 0.01 ms (blue) and 5 ms (green). (d) Tangentialdb/dt as a function of time at 300 m along the survey line (shownin the red dot in (a)). Similar information is shown in (e), (f), (g)and (h) for the model with the conductive casing. The differencein the db/dt data (casing minus halfspace) is shown in (i), (j), (k)and (l). The difference is also shown as a percentage of the halfspacesolution at 0.01 ms and 5 ms in (m) and through time at 300 m offsetin (n). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168xxxviiFigure 5.7 Current density for a top-casing time domain EM experiment over a10−2 S/m half-space and a 1 km-long, conductive, permeable steel-cased well (5×106 S/m, 100µ0). . . . . . . . . . . . . . . . . . . 171Figure 5.8 Difference in current density for a time domain experiment whichincludes a conductive, permeable steel-cased well (as in Figure 5.7)and an experiment over a conductive well (as in Figure 5.2). . . . . 172Figure 5.9 Sketch demonstrating how a vertical circulation of current can ariseinside of a permeable casing in a top-casing TDEM experiment. Asource current is applied and (a) currents flow downwards throughthe pipe. (b) Currents generate rotational magnetic fields accord-ing to Ampere’s law. (c) Magnetic flux is concentrated in the per-meable pipe according to the constitutive relation between~b and~h.(d) The magnetic flux varies with time after shut-off, and the time-varying magnetic flux creates rotational electric fields according toFaraday’s law. (e) Currents associated with those electric fields areconcentrated in conductive regions of the model in accordance withOhm’s law. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174Figure 5.10 Simulated electric field at the surface of the earth, as was presentedin 5.4. The top row (a – d) are the data simulated with the conductivecasing. The second rows (e – h) are the data simulated with theconductive, permeable casing. The third row (i – l) is the difference(permeable and conductive minus conductive only), and the fourthrow (m and n) show the difference as a percentage of the conductivesolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175xxxviiiFigure 5.11 Simulated db/dt at the surface of the earth, as was presented in 5.6.The top row (a – d) are the data simulated with the conductive cas-ing. The second rows (e – h) are the data simulated with the con-ductive, permeable casing. The third row (i – l) is the difference(permeable and conductive minus conductive only), and the fourthrow (m and n) show the difference as a percentage of the conductivesolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176Figure 5.12 Current density for a down-hole time domain experiment with a con-ductive, permeable steel-cased well (5×106 S/m, 100µ0). The pos-itive electrode is connected to the casing at 950 m depth and thereturn electrode is on the surface at x= 1000 m. . . . . . . . . . . 178Figure 5.13 Difference in current density for a down-hole time domain exper-iment with a conductive, permeable steel-cased well (as in Figure5.12) and a similar experiment with a conductive well. . . . . . . . 179Figure 5.14 Simulated electric field at the surface of the earth, as was presentedin 5.10 for a TDEM experiment in which the positive electrode iscoupled to the casing at a depth of 950m. . . . . . . . . . . . . . . 180xxxixFigure 5.15 (a) Downward-directed current (a) within the casing at 0.01 ms (blue)and 5 ms (green) and (b) as a function of time at 500 m depth for thetrue, hollow-cased well. (c), (d) Current in the approximate modelwhich treats the casing as a solid cylinder. The conductivity of thecylinder approximating the well is chosen so that both models haveequal products of the conductivity and the cross sectional area. (e),(f) Difference between the current in the approximate model and thetrue model (approximate minus true). (g), (h) Difference as a per-centage of the true solution. . . . . . . . . . . . . . . . . . . . . . 183Figure 5.16 Simulated electric field at the surface of the earth, as was presentedin 5.4. The top row (a – d) are the data simulated with the true,conductive casing. The second rows (e – h) are the data simulatedwith the solid cylinder approximating the well. The third row (i –l) is the difference (approximate minus true), and the fourth row (mand n) show the difference as a percentage of the true solution. . . . 184Figure 5.17 (a) True current density (b) difference between the currents in thesimulation with the approximation to the conductive well and thetrue solution (secondary currents), (c) True charge density, (d) sec-ondary charge density, (e) true electric field, (f) secondary electricfield. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186Figure 5.18 Sketch of the geometry of the dominant geometry of the currentdensity, ~j and magnetic flux density,~b, within the casing. Most ofthe current flows vertically while most of the magnetic flux densityis rotational. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188xlFigure 5.19 (a) Downward-directed current (a) within the casing at 0.01 ms (blue)and 5 ms (green) and (b) as a function of time at 500 m depth for thetrue, permeable, hollow-cased well. (c), (d) Current in the approxi-mate model which treats the casing as a solid cylinder. The conduc-tivity of the cylinder approximating the well is chosen so that bothmodels have equal products of the conductivity and the cross sec-tional area. The permeability of the cylinder is chosen so that bothmodels have equal products of the permeability and the thickness.(e), (f) Difference between the current in the approximate modeland the true model (approximate minus true). (g), (h) Difference asa percentage of the true solution. . . . . . . . . . . . . . . . . . . 189Figure 5.20 Simulated electric field at the surface of the earth, as was presentedin 5.4. The top row (a – d) are the data simulated with the true,conductive, permeable casing. The second rows (e – h) are the datasimulated with the solid cylinder approximating the well. The thirdrow (i – l) is the difference (approximate minus true), and the fourthrow (m and n) show the difference as a percentage of the true solu-tion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190Figure 5.21 (a) True current density, (b) difference between the currents in thesimulation with the approximation to the conductive, permeable welland the true solution (secondary currents), (c) true charge density,(d) secondary charge density, (e) true electric field, (f) secondaryelectric field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191xliFigure 6.1 Approximation to a step function using an arctangent, as describedin equation 6.6. The values in the legend indicate the value of theslope, a. Smaller values have a more gradual transition between zeroand one, while larger values have a more rapid transition. . . . . . 200Figure 6.2 Model of an electrically conductive propped fracture zone (3 S/m) ina halfspace (100 Ω m) with a steel-cased well. The well is modelledas a solid cylinder with a conductivity of 1.4×104 S/m. The meshhas 4 cells across the radius of the casing. The fractured regionextends vertically from 950 to 960 m depth and has a radius of 50 m.Panel (a) has a radial extent of 80 m to show the fractured zone andpanel (b) has a radius of 0.4 m to show the casing. The a-electrodelocations are shown in panel (b). . . . . . . . . . . . . . . . . . . . 206Figure 6.3 (a) Synthetic data for the down-hole casing experiment for the back-ground, prior to the fracture (solid lines), and after the fracture (dots),(b) secondary electric field (fracture - background), and (c) sec-ondary electric field as a percentage of the primary (background).The color of the lines or dots indicates the depth of the source. . . . 207xliiFigure 6.4 (a) Observed and predicted radial electric field data, (b) differencebetween the observed and predicted data (V/m), (c) difference be-tween the observed and predicted data as a percentage of the ob-served data, and (d) conductivity model recovered in the inversion.The colors in (a), (b), and (c) indicate the source location as shownin Figure 6.3. The white outline in (d) outlines the true geometry ofthe fracture zone and the grey line shows the location of the well-bore. The data are fit to a global χ-factor < 0.1 and the inversiontook 3 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 209Figure 6.5 Integrated sensitivity for (a, b) a survey conducted in a 100Ωm half-space and (c, d) the true model with the conductive target, shown inFigure 6.2. The inversion model is log-conductivity. The casingis shown by the grey line, and the source locations are shown bythe red dots. The integrated sensitivity is computed by squaringthe elements of the sensitivity matrix, summing along the rows ofthe sensitivity matrix and then taking the square-root. This gives avector whose length is the number of model parameters. . . . . . . 210Figure 6.6 Tikhonov inversion result, similar to that shown in 6.4 that fits thedata to a global χ-factor < 0.05. The inversion took 5 iterations. . . 211Figure 6.7 Tikhonov inversion result, similar to that shown in 6.6 that uses αx=100. The inversion took 5 iterations. . . . . . . . . . . . . . . . . . 212xliiiFigure 6.8 Parametric inversion result for a starting model centered at 980 mdepth with a thickness of 5 m and a radius of 5 m. The initial back-ground conductivity is 10−2 S/m and the initial conductivity of thetarget is 3×10−2 S/m. The geometry of the starting model is shownby the white dashed-lines and the true model is shown by the solidwhite outline. The inversion reached a χ-factor < 0.05 and took 8iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213Figure 6.9 Parametric inversion result for a starting model centered at 955 mdepth with a thickness of 5 m and a radius of 5 m. The initial back-ground conductivity is 10−2 S/m and the initial conductivity of thetarget is 3×10−2 S/m. The geometry of the starting model is shownby the white dashed-lines and the true model is shown by the solidwhite outline. The inversion reached a χ-factor < 0.05 and took 8iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214Figure 6.10 Parametric inversion result for a starting model centered at 955 mdepth with a thickness of 10 m and a radius of 5 m. The initialbackground conductivity is 10−2 S/m and the initial conductivity ofthe target is 3× 10−2 S/m. The geometry of the starting model isshown by the white dashed-lines and the true model is shown by thesolid white outline. The inversion reached a χ-factor < 0.05 andtook 12 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . 215xlivFigure 6.11 Parametric inversion result where the center of the target is fixed ata depth of 955 m. The starting model has a thickness of 5 m and aradius of 75 m. The initial background conductivity is 10−2 S/m andthe initial conductivity of the target is 3×10−2 S/m. The geometryof the starting model is shown by the white dashed-lines and the truemodel is shown by the solid white outline. The inversion reached aχ-factor < 0.05 and took 6 iterations. . . . . . . . . . . . . . . . . 216Figure 6.12 Self-consistent effective medium theory estimation of electrical con-ductivity as a function of fracture concentration. The conductivityof the background is 10−2 S/m and the conductivity of the materialfilling the fractures is 2500 S/m. We use an aspect ratio of 3×10−5,and the fractures are assumed to be randomly oriented. Panels (a)and (b) show the conductivity on a linear scale while panels (c) and(d) show the conductivity on a logarithmic scale. Panels (a) and (c)show the full range of possible ϕ values from 0 to 1, and panels (b)and (d) zoom into a smaller range, 0 ≤ ϕ ≤ 0.005. . . . . . . . . . 218Figure 6.13 Integrated sensitivity, similar to that shown in 6.5, except here, theinversion model is fracture concentration. Panels (a, b) are the sensi-tivity for a ϕ = 0 half-space and (c, d) show the sensitivity calculatedusing the true model which includes the fracture. . . . . . . . . . . 219xlvFigure 6.14 (a) Observed and predicted radial electric field data, (b) differencebetween the observed and predicted data (V/m), (c) difference be-tween the observed and predicted data as a percentage of the ob-served data, (d) fracture concentration model, ϕ , recovered in theinversion, and (e) electrical conductivity model obtained by con-verting ϕ to electrical conductivity using equation 6.9 The colors in(a), (b), and (c) indicate the source location as shown in Figure 6.3.The white outline in (d) outlines the true geometry of the fracturezone and the grey line shows the location of the wellbore. The dataare fit to a global χ-factor < 0.05 and the inversion took 5 iterations. 220Figure 6.15 Inversion result using a parametric model of the fracture concentra-tion. The black dashed line outlines the geometry of the startingmodel. The initial fracture concentration is 10−10 in the backgroundand 10−4 in the target. The data are fit to a global χ-factor < 0.05and the inversion took 11 iterations. . . . . . . . . . . . . . . . . . 221Figure 6.16 Inversion result using a parametric model of the fracture concentra-tion. The concentration is converted to electrical conductivity usingthe effective medium theory mapping shown in equation 6.9. Theblack dashed line outlines the geometry of the starting model. Theinitial fracture concentration is 10−10 in the background and 10−4in the target. The data are fit to a global χ-factor < 0.05 and theinversion took 12 iterations. . . . . . . . . . . . . . . . . . . . . . 222xlviFigure 6.17 Parametric inversion result where the center of the target is fixedat a depth of 955 m. The concentration is converted to electricalconductivity using the effective medium theory mapping shown inequation 6.9. The black dashed line outlines the geometry of thestarting model. The initial fracture concentration is 10−10 in thebackground and 10−4 in the target. The data are fit to a global χ-factor < 0.05 and the inversion took 6 iterations. . . . . . . . . . . 223Figure 6.18 Parametric inversion using effective medium theory and a volumedata-misfit in the statement of the inverse problem. The uncertaintyon the volume datum was 5%. The starting model was the same asused in Figure 6.17 The inversion reached a global χ-factor < 0.05after 5 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 226Figure B.1 Problem setup. Concentric spheres in a uniform electric field. . . . . 252Figure D.1 Inversion approach using the SIMPEG framework. Adapted fromCockett et al. (2015) . . . . . . . . . . . . . . . . . . . . . . . . . 283Figure D.2 Forward simulation framework. . . . . . . . . . . . . . . . . . . . 289Figure D.3 Location of variables in the finite volume implementation for both aunit cell in (a) cartesian and (b) cylindrical coordinates (after Heagyet al. (2015)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 290Figure D.4 (a) Contributions of each module to the sensitivity. (b) process forcomputing Jv and (c) J>v; stars indicate where the source deriva-tives are incorporated. . . . . . . . . . . . . . . . . . . . . . . . . . 292xlviiFigure D.5 Mapping an inversion model, a 1D layered, log conductivity modeldefined below the surface, to electrical conductivity defined in thefull simulation domain. . . . . . . . . . . . . . . . . . . . . . . . . 295Figure D.6 Diagram showing the entire setup and organization of (a) the fre-quency domain simulation; (b) the time domain simulation; and(c) the common inversion framework used for each example. Themuted text shows the programmatic inputs to each class instance. . . 306Figure D.7 (a) True and recovered models for the FDEM and TDEM inversions;predicted and observed data for (b) the FDEM example, and (c) theTDEM example. In (b) the magnetic field data are in the negativez-direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 309Figure D.8 400 Hz In-phase RESOLVE data at (left) and High Moment SkyTEMdata at 156 µ s. The white dot at (462100m, 6196500m) on bothimages is the location of the stations chosen to demonstrate the 1Dinversions in frequency and time. . . . . . . . . . . . . . . . . . . 311Figure D.9 (a) Models recovered from from the 1D inversion of RESOLVE(back) and SkyTEM (blue) data at the location (462100m, 6196500m).(b) Observed (lines) and predicted (points) frequency domain data.(c) Observed and predicted time domain data. (d) Source waveformused in for the SkyTEM inversion, the x-axis is time (µ s) on a linearscale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313Figure D.10 (a) Conductivity model 9.9m below the surface from a stitched 1Dinversion of RESOLVE data. (b) Real component of the observedRESOLVE data at 400Hz. (c) Real component of the predicted dataat 400Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314xlviiiFigure D.11 Setup of parametric models and calculation of the sensitivity for aprimary secondary approach of simulating 3D geology and steel cas-ing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316Figure D.12 Cross sectional slice of primary (casing + background) real currentdensity. The colorbar is logarithmically scaled and shows the am-plitude of the real current density. . . . . . . . . . . . . . . . . . . 320Figure D.13 Depth slice at z=-950m showing the source current density for thesecondary problem. . . . . . . . . . . . . . . . . . . . . . . . . . . 321Figure D.14 Simulated real electric field data as measured at the surface using aprimary secondary approach for casing and a conductive target (out-lined in white). The upper panels show the total Ex (a) and Ey (b);the lower panels show the secondary (due to the conductive block)Ex (c) and Ey (d). Note that the colorbars showing the secondaryelectric fields are not on the same scale. The limits of the colorbarshave been set so that the zero-crossing is always shown in the samecolor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322Figure D.15 Sensitivity of surface real Ex (left) and Ey (right) data with respectto the physical properties, ((V/m)/(log(σ))) . . . . . . . . . . . . 323Figure D.16 Sensitivity of surface real Ex (left) and Ey (right) data with respectto the layer geometry, ((V/m)/m) . . . . . . . . . . . . . . . . . . 324Figure D.17 Sensitivity of surface real Ex (left) and Ey (right) data with respectto the block geometry, ((V/m)/m) . . . . . . . . . . . . . . . . . . 325xlixFigure E.1 Web-page on a conducting sphere in a uniform electric field (https://em.geosci.xyz/content/maxwell2 static/fields from grounded sources dcr/electrostatic sphere.html). The plot shown on this page is compiledfrom Python code and can be downloaded and run by users, or up-date by maintainers if there are improvements that can be made. . . 333Figure E.2 Text source from which the web-page shown in Figure E.1 is built.The code used to build the plot starts under the plot-directive on line134. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334Figure E.3 A sample Jupyter notebook “app” that simulates a DC resistivityover a halfspace with a cylinder. The widgets control the visualiza-tion. The user can change the model parameters (e.g. the resistivityof the cylinder and of the background), the survey geometry, andselect the visualization (e.g. electric field, currents, charges, etc.) . 336lFigure E.4 Seven step framework used for the case histories in https://em.geosci.xyz.(1) Setup: describe the geoscientific question and objectives thatare to be addressed. (2) Properties: identify the diagnostic physicalproperties (e.g. density, electrical conductivity, magnetic suscepti-bility, etc.). (3) Survey: design a survey that is suitable for detectingthe physical property contrasts relevant to the application. (4) Data:carry out the field survey and collect the data set(s). (5) Processing:plot the data and apply the analysis steps needed to interpret the data(e.g. invert the data). (6) Interpretation: interpret the results in termsof the identified physical properties and original geoscience objec-tive. (7) Synthesis: integrate the interpretation with geologic andother information relevant to the application in order to address theoriginal geoscience objective. This image and caption are adaptedfrom: https://em.geosci.xyz/content/geophysical surveys/fundamentals/seven steps.html . . . . . . . . . . . . . . . . . . . . . . . . . . . 340liAcknowledgmentsI am very fortunate to have been supported by a tremendous group of people during myacademic and personal journeys as I conducted this research.First and foremost, I am deeply grateful to my supervisor Dr. Doug Oldenburg – yourendless enthusiasm and desire to make a positive impact, both in the field of geophysicsand the daily lives of those around you, have enriched my life. You have opened doorsfor me and encouraged me to pursue projects I find meaningful. For these things andmore, thank you.Next, I thank Dr. Eldad Haber; the numerical work in this thesis has its foundationin your courses and resources on computational electromagnetics. Thank you for all youhave taught me. Thanks to Michael Bostock for serving on my committee and for yourinput which helped improve this work.To Dr. Michael Wilt, conversations with you helped shape the questions I pursuedand the final content of this thesis. I am grateful for the investment you have made inme.There are numerous educators who have inspired me and motivated my interest inpursuing science. In particular, Linda Legendre, David Westra, and Dr. David Lawrie,you have each had a meaningful impact on my career.A special thanks to the friends who have supported me through my PhD. To name aliifew: Dom Fournier, Sarah Devriese, Thibaut Astic, Gudni Rosenkjaer, Mike McMillan,Mike Mitchell, Devin Cowan, Luz Angelica Caudillo-Mata, Dikun Yang, Kris Davis,Dave Marchant, Jenn Fohring, Justin Granek, and Julie Nutini. I especially thank RowanCockett and Seogi Kang, I have learned so much from you both; I am grateful to countyou as partners in much of this work and as supportive friends. To the growing numberof contributors to SimPEG and resources in the open source ecosystem, I am inspiredby your efforts.To my family, who have been unwavering in their support including my encouragingbrother, Josh Heagy, and my grandparents.Finally, to my parents, I could not have done this without you. I am forever gratefulfor your love and support.Funding for this work was provided by Vanier Canada Graduate Scholarship Pro-gram, the University of British Columbia Four Year Fellowship.liiiDedicationTo my parents.livChapter 1IntroductionThe purpose of computing is insight, not numbers.— Richard W. Hamming (1962)Numerical Methods for Scientists and Engineers (Preface)1.1 MotivationExploration and development of unconventional resources such as shale gas formationshas increased significantly over the past decades. The combination of horizontal drillingand hydraulic fracturing are key technologies for extracting hydrocarbons from shaleand low permeability, “tight” reservoirs. As a byproduct of fracturing, wastewater, if notrecycled, must be disposed of. In many cases this is achieved by injecting the wastewaterinto the subsurface through an injection well. Similarly, with recognition of the hazardthat carbon-dioxide poses as a contributor to global warming, efforts are being madeto capture and store CO2 in the subsurface. In each of these scenarios, there are both1environmental and economic motivations for characterizing the distribution of materialsinjected into the subsurface.Electrical conductivity is a potentially diagnostic physical property in each of theseapplication; the conductivity of the injected material is often distinct from the host rockit is being injected into. Electrical and electromagnetic geophysical techniques thereforehave the potential of being viable methods for estimating the geometry of injected ma-terials. To focus the research questions addressed in this thesis, I take the application ofhydraulic fracturing as the primary motivating context, keeping in mind the connectionsand similarities with other subsurface injections. In particular, reservoir settings oftenhave steel-cased wells which complicate the electrical and electromagnetic responses.In the sections that follow, I provide background and context on hydraulic fracturingand the questions I aim to address using electrical and electromagnetic (EM) methods,and discuss the challenges introduced when attempting to use electrical and electro-magnetic methods in settings with steel-cased wells. At the end of this introductorychapter, I briefly introduce fundamental concepts of electromagnetism and geophysicalinverse problems which are the foundation for my research and are applicable to all ofthe chapters in the thesis.1.2 Hydraulic FracturingHydraulic fracturing is used to extract hydrocarbons from tight (low-permeability) andshale formations where oil and gas will not easily flow. In such settings, hydraulicfracturing is used to create pathways for the hydrocarbons to flow (Figure 1.1). Theprocess of inducing a fracture involves sealing off a section of the well and pumpingfluid into that section under high pressure until the rock fails and cracks open up in thedirection of the minimum principal stress. Typically, once the rock has fractured, sand2or ceramic particles, referred to as proppant, are pumped into the formation to keep thenewly created pathways open. Many of the wells drilled in the past two decades arehorizontal wells, and typically 15 to 30 fracture stages (in some cases up to 60) areperformed along the length of the well (Maxwell, 2014).The extent of the fractures, complexity of the fractures, and distribution of proppantand fluid within these fractures are key factors that determine the available pathways foroil or gas to flow to the well, and thus, are important for the resulting production of theresource (Brannon and Starks, 2008; Cipolla et al., 2009). One of the industry termsused to try to quantify the extent of the fractures within the reservoir is the “StimulatedReservoir Volume” (SRV). Estimates of the SRV are used to make engineering decisionson the spacing between wells and completion design such as the number and spacingbetween fracture stages along the length of a well (Palisch et al., 2016). Hoversten et al.(2015) estimate that a 5% improvement in the characterization of the SRV for a 1 billionbarrel field translates to over 0.5 billion U.S. dollars over 24 years with oil at US$50 perbarrel. Typically, an analysis to estimate the SRV is conducted using microseismic and/ or tiltmeter data. Microseismic is sensitive to acoustic events generated as the fracturepropagates through the reservoir (c.f. Mayerhofer et al. (2010); Cipolla et al. (2009);Maxwell et al. (2002); Warpinski (1996)), while tiltmeters are used to characterize thedeformation of the rock due to the presence of a fracture or a change in the stress dis-tribution (Mayerhofer et al., 2010; Wright et al., 1998). Though both techniques arevaluable for estimating fracture geometry, complexity, and orientation, they provide nomechanism for delineating the extent of the proppant within the fractures; at most, theyare a proxy for the extent of the fluid (Palisch et al., 2016). Two of the key factors infracture performance are the fracture area and the hydraulic conductivity (Cipolla andWallace, 2014), both of these are controlled by the distribution of proppant. Therefore,3Figure 1.1: Conventional reservoirs contain oil and gas that have migrated up-wards under pressure until they are trapped by a cap rock (left), while non-conventional tight or shale oil and gas reservoirs contain hydrocarbons thatare trapped in low permeability formations (right).in order to assess how effectively the fracture treatment has stimulated the reservoir andhow efficiently resources, such as water, have been used to create the fracture, we needa method to delineate the extent of the proppant and fluid within the reservoir.To accomplish this task, I propose to use EM geophysical techniques. For EM tobe a viable method for imaging the distribution of proppant and fluid within a fracturedvolume of rock, we require that: (1) the fractured volume of rock have physical proper-ties which are distinct from the background, host rock, (2) the survey must be sensitiveto this contrast, and (3) once the data have been collected, they must be interpreted orinverted in a meaningful manner. Note that this is a time-lapse problem; that is, by in-ducing a fracture, the physical properties of the reservoir have been altered. In order tocharacterize such a change, we must view the imaging problem as a time-lapse one, andcollect data to provide us with a before and an after data-view of the reservoir.Variations in subsurface electrical conductivity have been used as a diagnostic phys-4ical property in sedimentary settings for characterizing geologic formations, and theproperties and distribution of fluids within those formations. Hydrocarbons are muchmore resistive than saline formation fluids. In enhanced oil recovery projects, fluidsare injected into the formation, which may be less resistive than the hydrocarbons theyreplace. These contrasts have been the target of cross-well, surface-to-borehole andborehole-to-surface electromagnetic (EM) methods for reservoir monitoring and char-acterization applications (cf. Bevc and Morrison (1991); Wilt et al. (1995); Marsalaet al. (2008, 2011, 2014)).In the case of hydraulic fracturing, the physical properties of the fractured volumeof the reservoir depend upon the properties of the injected fluid and proppant particles.Saline water may be used, as is often the case when recycled water is used, and elec-trically conductive proppant may be manufactured and injected (Cannan et al., 2014;Vengosh et al., 2014; King, 2010). One or both of these may be used to create a phys-ical property contrast between the host reservoir rock and the fractured volume of thereservoir. This contrast is what we aim to excite using an electromagnetic survey.Using additives to make a hydraulic fracture a geophysical target is not a new idea.Byerlee and Johnston (1976) suggested using magnetic particles and a magnetic surveyto estimate fracture orientation at distances larger than can be determined by tracers andwell logs. To create a fracture with a significant magnetic susceptibility, they suggestedusing finely crushed magnetite suspended in the fracturing fluid or iron shot, sphericalparticles of iron which do not crush easily and have a high magnetic susceptibility.They treated the fracture as a plate with a known magnetic susceptibility and collectedmeasurements before and after the fracture operation; measurements of the horizontalmagnetic field inside of the injection borehole are used to indicate the orientation of thefracture and surface measurements are used to estimate the fracture geometry. Using an5analytic model for the magnetic response of a circular disc in a uniform inducing field,they demonstrated the potential of this technique for simple analog models of two fieldsites: an engineered geothermal project at Los Alamos, where fractures are induced tocirculate fluid through a dry geothermal reservoir and a hydraulic fracture operation fora tight-gas reservoir at Rio Blanco.Similarly, Bartel et al. (1976) suggest using electrically conductive fracturing fluidand measuring electric potentials at the surface in a direct current (DC) resistivity ex-periment where one electrode is connected to the injection well and a return electrode isconnected to a distant well casing. Potential electrodes are arranged in two concentricrings centered about the injection well with potential differences being measured be-tween electrodes along the same azimuth in the inner and outer rings. Variations in theamplitude of the potential difference with azimuth are used to estimate the orientationof the fracture and to detect asymmetry (e.g. if the fracture extends further one sideof the well than the other). Field tests were performed at two sites and demonstrateddetectability of fractures in an experiment in the Wattenburg field for fractures at depthsof ∼2200 m depth. The source electrode was connected to the casing and centered atthe depth of the induced fracture (see also Smith et al. (1978)).Since these initial developments in the late 70’s, there have been significant im-provements in data quality as well as our ability to model and invert electrical and elec-tromagnetic data in 3D. In conjunction, there have been advancements in fracture oper-ations and imaging techniques such as the use of microseismic, tiltmeters, and pressure-transient analysis to characterize hydraulic fracture. As a result, the questions beingasked are more detailed. How complex is the fracture? Is it a simple planar fracture oran extensive network? What is the extent of the propped volume of rock?To assess if EM can provide insights for helping answer these questions, we need6to be able to perform 3D numerical simulations of EM for fractured reservoirs. Thisprompts a question to be examined in the thesis:• How do we construct a physical property model of a fractured volume of rock thatcan be used in numerical simulations?1.3 Challenges to electromagnetics with steel-casedboreholesHydraulic fracture operations are typically conducted in steel cased wells. These wellspresents several significant challenges for modelling and inverting EM geophysical data.Most wells are cased with carbon steel, which has both a large electrical conductivity (∼5.5×106 S/m) and magnetic permeability (> 50µ0) (Wu and Habashy, 1994). This is alarge contrast compared to typical geologic settings, which typically have conductivitiesless than 1 S/m and permeabilities similar to that of free space, µ0. As a result, the casingwill have a significant impact on the behavior of the EM fields and fluxes. The propertiesof the casing, in particular, magnetic permeability, as well as casing thickness, can varyalong the length of the well. They also depend on the quality of the steel and the stateof corrosion; this adds another level of complexity to the situation. Well logging tools,such as that described by Brill et al. (2012) have been developed to characterize thesevariations.Not only does casing introduce a large, variable, physical property contrast, its ge-ometry and scale add to the difficulty. Well casing is cylindrical in shape, typically ∼1cm in thickness, 10− 20 cm in diameter, and may extend several kilometers in depth,while the geologic structures we aim to characterize with the geophysical survey havethree dimensional variations in electrical conductivity on the scale of hundreds of metersto kilometers. Approaches tailored to accurately modelling the impact of the casing in7simple 1D geologic settings may be geometrically incompatible or computationally pro-hibitive to use when reservoir-scale, three-dimensional geologic structures and targetsare included in the conductivity model.In much of the early literature, the casing was viewed as a nuisance which distortsthe EM signals of interest. Distortion of surface DC resistivity and induced polarization(IP) data, primarily in hydrocarbon settings, was examined in Wait (1983); Holladay andWest (1984); Johnston et al. (1987) and later extended to grounded source EM and IPin Wait and Williams (1985); Williams and Wait (1985); Johnston et al. (1992). Also inhydrocarbon applications, well-logging in the presence of steel cased boreholes is mo-tivation for examining the behavior of electromagnetic fields and fluxes in the vicinityof casing. Initial work focussed on DC resistivity with Kaufman (1990); Schenkel andMorrison (1990); Kaufman and Wightman (1993); Schenkel and Morrison (1994), andinductive source frequency domain experiments with Augustin et al. (1989a). Kaufman(1990) derives an analytical solution for the electric field at DC in an experiment wherean electrode is positioned along the axis of an infinite-length well. The mathematical so-lutions presented show how, and under what conditions, horizontal currents leak into theformation outside the well. Moreover, Kaufman (1990) showed, based upon asymptoticanalysis, which fields to measure inside the well so that information about the formationresistivity could be obtained. This analysis is extended to include finite-length wells inKaufman and Wightman (1993). Schenkel and Morrison (1994) show the importance ofconsidering the length of the casing in borehole resistivity measurements, and demon-strate the feasibility of cross-well DC resistivity. They also show that the presence of asteel casing can improve sensitivity to a target adjacent to the well. In frequency domainEM, Augustin et al. (1989a) consider a loop-loop experiment, where a large loop is po-sitioned on the surface of the earth and a magnetic field receiver is within the borehole.8Magnetic permeability is included in the analysis and a “casing correction”, effectively afilter due to the casing on inductive-source data, is introduced. This work was built uponfor considering cross-well frequency domain EM experiments (Uchida et al., 1991; Wiltet al., 1996).For larger scale geophysical surveys, steel cased wells have been used as “extendedelectrodes.” Rocroi and Koulikov (1985) used a pair of well casings as current elec-trodes for reservoir characterization in hydrocarbon applications. In near-surface set-tings Ramirez et al. (1996); Rucker et al. (2010); Rucker (2012) considered the useof monitoring wells as current and potential electrodes for a DC experiment aimed atimaging nuclear waste beneath a leaking storage tank. Similarly, Ronczka et al. (2015)considers the use of groundwater wells for monitoring a saltwater intrusion and investi-gates numerical strategies for simulating casings as long electrodes. Imaging hydraulicfractures has been a motivator for a number of studies at DC or EM, among them Weisset al. (2016); Hoversten et al. (2017). Some of these have suggested the use of casingsthat include resistive gaps so that currents may be injected in a segment of the well andpotentials measured across the other gaps along the well (Nekut, 1995; Zhang et al.,2018).There has also been an increase in interest in examining the use of electrical orelectromagnetic methods deployed on the surface to non-invasively look for flaws orbreaks in the casing. Wilt et al. (2018b) introduces the idea of using electrical or elec-tromagnetic methods for casing integrity which is further expanded upon in Wilt et al.(2018a). They show that low-frequency electromagnetic methods are sensitive to vari-ations in wellbore length and demonstrate that their numerical simulations agree withfield data collected over two different wellbores at the Containment and Monitoring In-stitute (CaMI) field site in southern Alberta, Canada. This work provides motivation for9further delving into the physics and assessing under which circumstances we can expectto detect a flaw along a wellbore using electrical or electromagnetic methods.As computing resources increased, our ability to forward-simulate more complexscenarios has improved. However, the large physical property contrasts and disparatelength scales introduced when a steel cased well is included in a model still presentcomputational challenges. Even the DC problem, which is relatively computationallylight, has posed challenges; those are exacerbated when solving the full Maxwell equa-tions in the frequency (FDEM) or time domain (TDEM) and can become crippling foran inversion. For models where the source and borehole are axisymmetric, cylindri-cal symmetry may be exploited to reduce the dimensionality, and thus the number ofunknowns, in the problem (e.g. Pardo and Torres-Verdin (2013); Heagy et al. (2015)).To reduce computational load in a 3D simulation, a number of authors have em-ployed simplifying assumptions. Several authors replaced the steel-cased well with asolid borehole, either with the same conductivity as the hollow-cased well (e.g. Umet 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 maybe used; Haber et al. (2016) similarly replaces the borehole with a coarser conductiv-ity structure and adopts an OcTree discretization to locally refine the mesh around thecasing. Yang et al. (2016a) uses a circuit model and introduces circuit components toaccount for the steel cased well in a 3D DC resistivity experiment. Also for DC resistiv-ity, Weiss (2017) develops a hierarchical finite element approach in which conductivityvalues can be defined on edges, facets, and volumes on the mesh. In doing so, mul-tiple thin boreholes, which may be deviated or horizontal, can be accurately modelledin a DC experiment. Another approach has been to replace the well with an “equiva-lent source”, for example, a collection of representative dipoles, inspired from Cuevas10(2014b), or with linear charge distributions for a DC problem (Weiss et al., 2016). Forthe frequency domain electromagnetic problem, a method of moments approach, whichreplaces 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 explicitlydiscretize the casing have been demonstrated, and they have been achieved at signifi-cant computational cost. Recent examples, including Commer et al. (2015); Um et al.(2015); Puzyrev et al. (2017), perform time and frequency domain simulations withfinely-discretized boreholes; they required the equivalent of days of compute-time for asingle forward simulation to complete. While these codes will undoubtedly see improve-ments in efficiency, there remains a need to accurately discretize the casing at modestcomputational cost both to serve as tool which these codes can compare against as wellas a need for researchers to be able to interact with and explore the behavior of the cur-rents, charges, electric and magnetic fields and fluxes to develop an understanding of thephysics in these high-contrast settings.There are a host of research questions to be explored in the context of highly con-ductive, permeable, steel cased wells in electromagnetic applications. This thesis is con-cerned with developing an understanding of the physics of electromagnetics in settingswith steel cased wells. In order to delve into the physics, we must have a mechanism forsimulating Maxwell’s equations which raises the question:• How do we develop a numerical approach that allows us to sufficiently discretizethe steel casing so that we can not only simulate DC and EM data, but also explorethe finer details of the currents, charges, electric and magnetic fields?Assuming a numerical strategy can be developed, we can then begin to explore ques-tions about the physics. At the DC limit,11• How do the currents, charges and electric fields behave in settings with steel casedwells? What are the implications for survey design?• Are there physical insights which, when incorporated into numerical algorithms,reduce the computational cost of simulations with steel-cased wells?Finally, electromagnetic experiments introduce the complication of time-variation inthe fields and fluxes as well as the additional concern of variable magnetic permeabilityin the model. Though the influence of the magnetic permeability of steel-casings hasbeen explored for inductive source EM (e.g. where the source is a circular loop ormagnetic dipole) (Wu and Habashy, 1994), the influence of magnetic permeability ongrounded-source EM surveys (e.g. where two electrodes are used to inject current intothe earth, similar to a DC resistivity experiment) is relatively unexplored. This promptsquestions like:• How do the fields and fluxes behave in an EM experiment which includes a steel-cased well?• What impact does magnetic permeability have on the behavior of the fields andfluxes? How important is it to consider in numerical simulations?These questions will serve as motivation for three chapters on steel-cased wells inthe thesis.1.4 Research motivation and thesis structureThe common theme throughout the thesis is the motivation to image subsurface injec-tions – specifically proppant and fluid injected during a hydraulic fracturing operation.To this end, I consider three main topics that comprise five chapters in the thesis:12• Building a physical property model for a fractured volume of rock (Chapter 2)• Developing a physical understanding of electrical and electromagnetic methodswhen steel-cased wells are present (Chapter 3, 4, and 5)• Formulating and solving the inverse problem for a fractured volume of rock (Chap-ter 6)Each of these topics is discussed in further detail below.The research chapters in this thesis make use of fundamental knowledge that needs tobe presented. Rather than insert that introductory material into the individual chapters,they are included at the end of this chapter. The background elements pertain to:• Electromagnetic geophysics: introducing Maxwell’s equations and EM geophys-ical surveys• Geophysical inversions: formulating the inverse problemThese are successively presented in Sections 1.5 and 1.6 after I have introduced theresearch material in this thesis.1.4.1 Homogenization strategy for hydraulic fracturesChapter 2 develops a strategy for estimating a bulk conductivity of a fractured volume ofrock. I break the problem into two steps, first estimating the effective conductivity of amixture of electrically conductive proppant and fluid and second, estimating the conduc-tivity of a fractured volume of rock where the fractures are filled with the proppant-fluidmixture. Having constructed a physical property model for a modest-sized fracture op-eration, I then demonstrate feasibility of detecting EM signal sensitive to the fractureusing a well-established cross-well EM survey. At this stage, I neglect the influence of13steel-cased wells. The aim of this chapter is to provide a reasonable estimate for theconductivity of a fractured volume of rock; future lab studies may provide an improvedrelationship between the conductivity of the proppant, fluid, host rock and their relativeconcentrations.The homogenization strategy taken in this chapter is based on established methodsin effective medium theory (Bruggeman, 1935) which have previously been consid-ered for fractured rocks (Berryman and Hoversten, 2013). The main contributions ofthis chapter are: algorithmic improvements to the effective conductivity calculation pre-sented by Berryman and Hoversten (2013), and the discussion of the two-stage workflowwhich allows various proppant-fluid mixtures to be considered (e.g. a fully conductiveproppant-pack, mixture of standard sand or ceramic proppant and conductive proppant).1.4.2 Steel casingsThe intention in this thesis is to advance our understanding of the physics of EM in set-tings with steel-cased wells. This understanding can then be used to develop appropriateapproximations in order to reduce computational cost, and to help calibrate our expecta-tions of the results of more advanced numerical techniques. To this end, I consider threetopics, each corresponding to one chapter in the thesis, to advancing our physical under-standing. This thesis does not contend with the challenge of efficiently simulating theDC or EM equations for settings with deviated or horizontal wells, nor do I consider set-tings with multiple wells or surface infrastructure. In a field survey, these complexitieswill need to be considered.Chapters 3, 4, and 5 contend with some of the challenges of performing DC andEM surveys in settings with steel-cased wells. In Chapter 3, I present a mimetic finitevolume implementation for simulating Maxwell’s equations at DC, in frequency and14in time, on cylindrically symmetric and 3D cylindrical meshes (which include an az-imuthal discretization). Although cylindrically symmetric meshes are common in EMmodelling, cylindrical meshes which include an azimuthal discretization are not. Acylindrical mesh conforms to the geometry of a vertical borehole and thus, can greatlyreduce the size of the mesh, and resultant cost of the computation, while still allowingfor significant refinement in the vicinity of the borehole. The purpose of this imple-mentation is to facilitate investigation into the physical behavior of EM fields and fluxesin this high-contrast setting. I consider two forms of validation: numerical validationand validation of the physical behavior. For numerical validation, I compare results ofa forward simulation with independently developed codes (Haber et al., 2007; Commeret al., 2015; Yang et al., 2016a). For validation of the physical behavior, I compare fea-tures of the simulated currents, charges, and electric fields to the behavior expected fromthe asymptotic analysis in Kaufman (1990); Kaufman and Wightman (1993) at DC, andsimilarly, compare the behavior of the magnetic fields and fluxes in the presence ofconductive and permeable pipes with the lab studies and analytical work contained inAugustin et al. (1989a) for inductive source frequency-domain electromagnetics.Chapter 4 focuses on developing an understanding of the physics of casings in a DCresistivity experiment. I start by considering the related application of casing integrity.In a casing integrity experiment, the aim is to detect a flaw or break in the casing whichmight be due to corrosion or failure under mechanical stresses. Wellbore integrity sur-veys are typically conducted with wireline tools, but more recently, Wilt et al. (2018b)introduced the idea of conducting a survey from the surface. This application providesthe opportunity to build a physical understanding about the behavior of the currents,charges, and electric fields in a DC experiment with casing. I investigate factors that in-fluence the feasibility of detecting a flaw, including the conductivity of the background15and casing, depth to the flaw, and impact on the signal if the whole circumference oronly a portion of the circumference of the casing is compromised. I also consider thelocation of the return electrode and its impact on the measured data. Many of the con-cepts developed in this section inform survey design elements in the next section, whichconsiders DC resistivity survey for a conductive or resistive target adjacent to the well.I look at survey design elements, including the impact of the source electrode locationon the resultant currents over the depth-interval of interest, and consider the impact ofparameters such as the conductivity of the target and whether the target is electricallyconnected to the casing or not. The final topic in this chapter considers approximationsto the steel-cased well in order to reduce the computational cost of the forward model.There is discrepancy among current literature as to whether the contrast between thebackground and the casing is the important quantity to conserve (e.g. Um et al. (2015))or if the product of the conductivity and the cross-sectional area of the casing shouldbe preserved (e.g. Swidinsky et al. (2013)). I compare both approaches to resolve thisquestion, demonstrate that the product of the cross-sectional area of the casing and itsconductivity is the important quantity to preserve in DC resistivity experiments, anddiscuss implications of choosing an incorrect approximation.Chapter 5 explores important aspects of electromagnetics in settings with steel casedwell. I begin by demonstrating the behavior of currents through time in time-domainEM simulation with a conductive well. The geometry of the currents through time israther complex, even in this seemingly simple experiment. Building from this, I thenconsider magnetic permeability and demonstrate how permeability alters the behaviorof the currents and in turn, discuss the impact permeability can have on data measured atthe surface. Finally, I revisit approaches for approximating steel cased wells developedat DC and explore some of the subtleties that arise in EM experiments.161.4.3 Inversions for subsurface injectionsChapter 6 returns to the application of hydraulic fracturing and examines the inverseproblem for a fractured volume of rock. Using a DC resistivity experiment, I exploreboth voxel-based and parametric strategies for extracting information about the frac-tured volume of rock. These inversions demonstrate some of the challenges caused bysteel cased wells and their influence on the sensitivity. Finally, I show how effectivemedium theory can be incorporated into the inversion and frame the inverse problemas an inversion for fracture concentration. This alters the sensitivity and provides theopportunity to incorporate a constraint on the volume of injected proppant and fluid.1.4.4 AppendicesThe thesis additionally includes 5 Appendices. Appendix A provides a listing of thesource code used to perform all of the examples shown in the thesis. Following this,there are two technical appendices, Appendix B and C, which support content in Chap-ters 2 and 6, respectively. Appendix D describes the open-source software used in thisthesis to simulate and invert Maxwell’s equations. The final appendix, Appendix E, dis-cusses how the open-source software development practices adopted in this thesis andthe broader SimPEG community are similarly used for the development of interactiveeducational resources for geophysics in the GeoSci.xyz project (https://geosci.xyz).1.4.5 A note on reproducibility and disseminationAlthough simulating and inverting data in settings with steel-cased wells has practicalapplications, within the realm of geophysics, and even within EM geophysics, it is anarrow topic. The value in focussing on specific applications is that each applicationhas its own subtleties and tangible details that must be worked through. In their own17right, the discoveries that are made by working through these details are of value tothe application-at-hand. However, this does not have to be the limit of their impact. Ifyou choose to take the perspective that the problem-at-hand is an instance of a largerclass of research questions, this allows you to view the methods and software neededto make incremental scientific discoveries as reusable components. It is then worththe investment of effort to develop and refactor these pieces so that they are modularbuilding blocks that fit into a larger framework and are of use to a broader community.This is the perspective I have sought to take as I conducted the research described in thisthesis.The simulations and inversions build on and contribute to the SimPEG framework,which is described in Cockett et al. (2015); Heagy et al. (2017a); Cockett (2017). Sim-PEG is written in python and is open-source, licensed under the permissive MIT license.Source code for all of the examples shown are also openly licensed and are available onGitHub as a combination of python scripts and Jupyter notebooks. Each chapter has acorresponding code-repository; Appendix A provides a summary of where these can beaccessed. My intention in providing all of the source-code is to be transparent aboutthe work completed; in addition, I hope these developments are of use and value to thebroader geophysics community.1.5 Background: Electromagnetic geophysicsBefore moving into the research content of the thesis, I will provide a brief overview ofthe governing equations in electromagnetics and recommend Ward and Hohmann (1988)and GeoSci.xyz (2018) for more in-depth discussion on the physics. For an overview ofhow to numerically solve Maxwell’s equations, I refer the reader to Haber and Ruthotto(2018).181.5.1 Governing equationsThe equations which govern the physics of EM are Maxwell’s equations:∇×~e=−∂~b∂ t∇×~h− ∂~d∂ t= ~j(1.1)where~e is the electric field (V/m),~b is the magnetic flux density (T),~h is the magneticfield (A/m), ~d is the electric displacement (C/m2), ~j is the current density (A/m2). Thefirst equation is Faraday’s law; it describes how a time-varying magnetic flux gener-ates rotational electric fields. The second equation is the Ampere-Maxwell equation(Maxwell’s contribution was the addition of the displacement current term, ∂ ~d/∂ t). Itdescribes how currents generate rotating magnetic fields.It is sometimes also convenient to consider Maxwell’s equations in the frequencydomain. Following Ward and Hohmann (1988), I use the eiωt Fourier transform conven-tion, namely,F(ω) =∫ ∞−∞f (t)e−iωtdtf (t) =∫ ∞−∞F(ω)e−iωtdω(1.2)where capital letters denote functions expressed in the frequency domain and lower-case letters denote functions expressed in the time-domain. In the frequency domain,Maxwell’s equations are given by∇×~E =−iω~B∇× ~H− iω~D= ~J(1.3)19These partial differential equations are coupled through the constitutive relationswhich relate the fields and fluxes through physical properties:~J = σ~E~B= µ~H~D= ε~E(1.4)where σ is electrical conductivity (S/m), µ is magnetic permeability (H/m), and ε is thedielectric permittivity (F/m). Electrical conductivity varies over many orders of magni-tude; two examples on the extreme ends are the conductivity of air,∼ 10−8 S/m, and theconductivity of steel, ∼ 106 S/m. The reciprocal of electrical conductivity is resistivity(ρ = 1/σ ), which has units of Ωm. Throughout the thesis, resistivity and conductivitywill be used interchangeably. The value of magnetic permeability is often taken to bethat of free-space, µ0 = 4pi × 10−7 H/m, in geophysical electromagnetic applications,as the permeability of most earth-materials ranges from 1-10 µ0 (Telford et al., 1990).The permeability of steel, however, is ∼ 100µ0; thus, for settings with steel-cased wellsthe simplifying assumption that µ = µ0 cannot necessarily be justified. The value ofdielectric permittivity varies between the free-space value of ε0 = 8.85×10−12 F/m andthe dielectric permittivity of water, 80ε0 (Telford et al., 1990). Variations in dielectricpermittivity are relevant at high frequencies, such as those used in ground penetratingradar (GPR) surveys (> 105 Hz), but in most geophysical electromagnetic surveys, late-time or low frequencies are employed as the high-frequency signals attenuate rapidly inconductive earth-materials. As such, the quasi-static approximation, which neglects dis-placement currents (the ∂ ~d/∂ t-term in equation 1.1 and the iω~D term in equation 1.3),is commonly adopted, and will be used throughout the thesis. Under the quasi-static20assumption, Maxwell’s equations are given by:∇×~e=−∂~b∂ t∇×~h= ~j(1.5)in time, and:∇×~E =−iω~B∇× ~H = ~J(1.6)in frequency.In the electrostatic limit, the time-derivative terms vanish. Taking the divergence ofAmpere’s law, and recognizing that the electric field~e is curl-free and can therefore beexpressed as the gradient of a scalar potential φ , gives us the governing equations forthe DC resistivity problem:∇ ·~j = I (δ (~r−~rs+)−δ (~r−~rs−))~e=−∇φ(1.7)Here ~j is the current density, I is the magnitude of the source current, and~rs+ and~rs−are the location of the positive and negative source electrodes, respectively. The electricfield and the current density are related through Ohm’s law (the first equation in equation1.4), which we can invoke to reduce the two first-order partial differential equations inequation 1.7 to a single, second order equation in φ :∇ ·σ∇φ =−I (δ (~r−~rs+)−δ (~r−~rs−)) (1.8)21This Poisson equation can be solved for the electric potential, φ , from which valuesof the currents and electric fields can be directly computed. In addition to consideringthe current density and electric fields, I will also present results in terms of chargesthroughout the thesis. The charge density is related to the electric field through∇ ·~e= ρ fε0(1.9)1.5.2 Geophysical surveysTo generate data, we require that a source be used to excite a response and that receiversmeasure the resultant fields and/or fluxes. Sources can be natural, plane-wave sources, asin the case of the magnetotelluric method, or they can be controlled, man-made sources.Broadly speaking, two categories of controlled sources exist: inductive sources andgrounded sources. Inductive sources consist of a loop of wire through which a time-varying current is passed, this in turn generates time-varying magnetic fields which actas the source fields, as shown in an airborne EM example in Figure 1.3. The time-varying magnetic field created by these systems induces currents in the conductive earthwhich, in turn, create secondary magnetic fields which can be measured at or abovethe surface. Within reservoir settings, cross-well EM is a commonly applied inductivesource survey (Wilt et al., 1995).Grounded sources require an electrical connection between the source and the earth.Two electrodes are positioned on, or in, the ground, the positive electrode injects cur-rent and a return electrode is a current-sink, as shown in Figure 1.4. In the electrostaticlimit, a grounded-source experiment is simply a DC resistivity experiment. Galvaniccurrents flow through the earth and cause charges to build up at conductivity interfaces,generating electric potentials which can be measured at the surface or within boreholes.22Figure 1.2: The choice of whether to solve the EM equations in the time-domain(TDEM) or the frequency domain (FDEM) depends upon the transmitterwaveform; a step-off-type waveform (left) lends itself to a solution approachin the time domain and a harmonic signal (right) is typically addressed in thefrequency domain.Figure 1.3: Time-domain inductive source experiment. (Left) a steady-state cur-rent is passed through the transmitter loop, generating a primary magneticfield. (Right) This current is shut-off, causing a change in magnetic flux.The changing magnetic flux induces secondary currents in conductors, whichin turn create secondary magnetic fields that can be measured at receiversabove, on, or within the earth.23Figure 1.4: (Left) Direct current resistivity experiment in which a steady-state cur-rent is injected into the earth, and (Right) a grounded source electromagneticexperiment which uses a time-varying source current.Electromagnetics introduces a time-varying component. The time-varying current flow-ing through the wire generates a time-varying primary magnetic field. This inducesvortex currents in conductive structures within the earth. Receivers can measure electricfields, magnetic flux, or time-variations in the magnetic flux or some combination ofthese.Two things are needed to generate data that are sensitive to a target or structureof interest. First, the source must be capable of getting EM energy to the target toexcite a response, and second, the response must reach the receivers with sufficientamplitude to be measurable. A measurable signal is one which is (a) above the noisefloor of the receivers, and (b) comprises a significant percentage of the primary (that is,the response with no target present). Numerical simulation of Maxwell’s equations isa critical component for assessing feasibility of detecting a target, and is an essentialelement of the inverse problem in which we aim to estimate a model from measureddata. Within the thesis, I discuss the machinery used to perform numerical simulationsin Chapter 3 and in Appendix D.241.6 Background: Geophysical InversionsOnce data have been collected, an inverse problem can be formulated. The goal of theinversion is to extract information about the subsurface from the data. Formulating, im-plementing and solving the inverse problem can be viewed as a workflow consisting ofinputs, implementation, and evaluation, as shown in Figure 1.5. The inputs are com-posed of the data, the governing equations, and prior knowledge or assumptions aboutthe setting. In the case of the fracturing problem, this may include well-log resistivitymeasurements which provide information about the background, knowledge of wherethe fracture was initiated, and the volumes of proppant and fluid pumped to create thefracture. The implementation consists of two broad categories: the forward simulationand the inversion. The forward simulation is the means by which we solve the governingequations given a model, and the inversion components evaluate and update this model.I consider a gradient based approach, which updates the model through an optimizationroutine. The output of this implementation is a model, which, prior to interpretation,must be evaluated. This requires considering, and often re-assessing, the choices andassumptions made in both the input and implementation stages (c.f. Oldenburg and Li(2005); Haber (2014); Cockett et al. (2015)).1.6.1 Formulating and solving the inverse problemIn this section, I provide a brief overview of geophysical inversions, adapted from Cock-ett et al. (2015); for more detail and examples, the reader is referred to Oldenburg andLi (2005); Cockett et al. (2015) as well as Appendix D for details specific to forwardand inverse modelling in electromagnetics.The aim of a geophysical inversion is to use the collected data to extract information25Figure 1.5: Overview of a geophysical inversion workflow. Adapted from Cockettet al. (2015).about the subsurface. In a given survey, a datum can be written asFi[m]+ εi = di (1.10)where F [·] is the forward simulation operator; for an electromagnetic problem, it sim-ulates Maxwell’s equations given a source and samples the relevant fields and fluxesat the receiver locations. The physical properties of the subsurface are captured by thevariable m, which I refer to as the inversion model. The noise is described by εi, anddi is the observed datum. A survey usually includes multiple sources and receivers, re-sulting in the observed data dobs = [d1, ...,dN ] and some estimate of their uncertainties– often assumed to be Gaussian. If the noise is Gaussian, then an appropriate measureof the data misfit is the l2-norm of the difference between the predicted data obtained26through a forward simulation and the observed data, namelyφd(m) =12‖Wd(F [m]−dobs)‖2 (1.11)Wd is a diagonal matrix whose elements are equal to Wdii = 1/εi where εi is an esti-mated standard deviation of the ith datum. Ideally, the noise model should also capturedeficiencies in the ability of the forward simulation to accurately describe the physics. Acommon choice in geophysical inversions is to assign a εi = f loor+%|di|. Percentagesare generally required when there is a large dynamic range of the data. A percentagealone can cause great difficulty for the inversion if a particular datum acquires a valueclose to zero, and therefore we include a floor.In addition to a metric that evaluates the size of the misfit, it is also required that wehave a tolerance, φ∗d ; models satisfying φd(m)≤ φ∗d are considered to adequately fit thedata (Parker, 1994). If the data errors are Gaussian and we have assigned the correctstandard deviations, then the expected value of φ∗d ∼ N/2, where N is the number ofdata. Note that the division by 2 is because the statement of the data misfit includes thefactor of 1/2 to simplify derivatives. When describing the “data-fit” of an inversion, it iscommon to quote a χ-factor, which is defined asφ invd = χφ∗d (1.12)Where φ invd is the final data-misfit in the inversion. Finding a model that has a misfitsubstantially lower than this will result in a solution that has excessive and erroneousstructure, that is, we are fitting the noise. Finding a model that has a misfit substantiallylarger than this will yield a model that is missing structure that could have been extractedfrom the data (see Oldenburg and Li (2005) for a tutorial).27The goal of an inversion is to estimate the earth-model, m from the data. In reality,the physical property distribution of the subsurface is complex and estimating this modelfrom a finite number of data is an ill-posed problem, meaning no unique model explainsthe data. Thus, in order to obtain a meaningful model from the data, assumptions andadditional information must be included. There are several mechanisms by which thiscan be achieved. One of the most common is to include consideration of a model regu-larization, φm in the inverse problem. This norm can penalize variation from a referencemodel, spatial derivatives of the model, or some combination of these. For example, theTikhonov-style regularization function (Tikhonov and Arsenin, 1977) can be expressedasφm(m) =αs2‖Ws(m−mref)‖2+ αx2 ‖Wxm‖2+αy2‖Wym‖2+ αz2 ‖Wzm‖2 (1.13)The first term is referred to as the “smallness” and penalizes difference between theinversion model and a reference model mref. The matrix Ws is a diagonal matrix; inthe simplest case it is the identity matrix. The remaining three terms are the first-ordersmoothness in the x, y, and z directions; the matrices Wx, Wy and Wz approximate thefirst order spatial derivatives in each direction. The α parameters weight the relativecontribution of each term to the regularization. Their values should consider the length-scales in the problem (c.f. Oldenburg and Li (2005)); for a typical 3D problem αx =αy = αz = 1 and αs is generally chosen to be several orders of magnitude smaller thanthe smoothness weights.To define the inverse problem, I take a deterministic approach to the inversion andtreat it as an optimization problem. Additional strong constraints on the model such asupper and lower bounds (mu, ml) are also considered. The general form of the objective28function I use combines the data misfit and regularization with a trade-off parameter, β ,between them, giving a problem of the formminimizemφ(m) = φd(m)+βφm(m)such that φd ≤ φ∗d , ml ≤m≤mu(1.14)Since the value of β is not known a priori, the above optimization problem can be solvedat many values of β to produce a trade-off, or Tikhonov, curve (cf. Parker (1994)). Anoptimum value, β ∗, can be found so that solving equation 1.14 with β ∗ produces a modelwith misfit φ∗d . One approach to finding the value of β∗ is to use cooling techniqueswhere the β is progressively reduced from some high value and the process stoppedwhen the tolerance is reached.The optimization problem posed in equation 1.14 is non-linear for DC resistivity andelectromagnetic forward simulations requiring that iterative optimization techniques beemployed (c.f. Nocedal and Wright (1999)). Gradient-based techniques are commonlyemployed. In particular, Gauss-Newton methods are effective in geophysical inversions.To ease notation, I consider a more compact description of the model regularization, andwrite our objective function asφ(m) =12‖Wd(F [m]−dobs)‖2+ β2 ‖Wmm‖2 (1.15)Note that if mref = 0 and Wm = [αsW>s ,αxW>x ,αyW>y ,αzW>z ]>, then the regularizationis equivalent to that stated in equation 1.13. The gradient is given byg(m) = J[m]>W>d Wd(F [m]−dobs)+βW>mWm(m) (1.16)29where J[m] is the sensitivity or Jacobian. The components J[m]i j specify how the ithdatum changes with respect to the jth model parameter,J[m] =dF [m]dm(1.17)I discuss the derivation of the sensitivity for time and frequency domain electromagneticproblems in depth in Appendix D.At the kth iteration, beginning with a model mk, we search for a perturbation δm thatreduces the objective function. Linearizing the forward simulation by Taylor expansion,F [mk+δm]' F [mk]+J[mk]δm+O(δm)2 (1.18)and setting the gradient equal to zero yields the standard Gauss-Newton equationsto be solved for the perturbation δm:(J[mk]>W>d WdJ[mk]+βW>mWm)δm =−g(mk) (1.19)The updated model is given bymk+1 = mk+ γδm (1.20)where γ ∈ (0,1] is a coefficient that can be found by a line search. Setting γ = 1 is thedefault and a line search is necessary if φ(mk+1)≥ φ(mk).The iterative optimization process is continued until a suitable stopping criterion isreached. Completion of this iterative process yields a minimization for particular valueof the trade-off parameter, β . If we are invoking a cooling schedule, and if the desired30misfit tolerance is not yet achieved, β is reduced and the iterative numerical optimizationprocedure is repeated.The forward simulation, computation of the sensitivities, and inversion machinerythat I use throughout this thesis are implemented in the open source software package,SimPEG (Cockett et al., 2015; Heagy et al., 2017a).31Chapter 2A physical property model for afractured volume of rock2.1 IntroductionFor electromagnetic methods to be sensitive to a propped, fractured volume of rock,the fractured volume of rock must have physical properties which are distinct from thebackground, host rock. For EM methods, this means the electrical conductivity, mag-netic permeability, or dielectric permittivity of the fractured rock must be distinct. Di-electric permittivity only plays a significant role when the frequency of the source issufficiently high, in the hundreds of kilohertz to megahertz range, as is used in groundpenetrating radar. Over the length-scales we need to consider for imaging fractures,hundreds to thousands of meters, attenuation of the EM signals due to skin depth effectsmake GPR impractical. Thus, we will work at lower frequencies, in the quasi-staticregime of Maxwell’s equations, and concern ourselves only with magnetic permeabilityand electrical conductivity. The materials traditionally used as proppant, typically sand32and ceramics, tend to have similar physical properties to the reservoir that they are beingpumped into, making it difficult to detect them on the scale of the reservoir. However, ifthe proppant were made electrically conductive or magnetically permeable, for instanceby coating it with graphite or including magnetite particles, it may create a sufficientphysical property contrast that can be imaged using EM.2.1.1 Proppant selectionTo make the proppant electromagnetically distinct from the host rock, either the mag-netic permeability or the electrical conductivity of the proppant could be targeted. Za-wadzki and Bogacki (2016) provides an overview of possible magnetic proppants andhighlights some practical considerations for the choice of material. For any proppant oradditive, the constraints include: the mechanical strength of the particles must be able towithstand the pressure of the closing fracture without crushing and clogging the fracturepathways; the material should not be toxic or reactive; and the price-point should bereasonable since high volumes are needed. Zawadzki and Bogacki (2016) also providea general classification of material types that could be considered: feedstock material,materials that are mixed with conventional proppant or replace conventional proppant,could include magnetite or steel particles. Both have a significant magnetic perme-ability, however, magnetite crushes easily, and although steel has significant mechani-cal strength, it is challenging to manufacture particles that are small enough (typically< 2mm in diameter). They also consider ferrofluids, which contain microscopic fer-romagnetic particles in suspension, and magnetic nanoparticles. Both are sufficientlysmall and remain in suspension so clogging is less of a concern. However, the cost ofeither material is quite significant and steps must be taken to reduce the environmentalhazard posed. This is particularly important for nanoparticles as they tend to be much33more reactive than particles of the same material but larger in size (Zawadzki and Bo-gacki, 2016). Continuing advances in nanotechnology has prompted some authors, e.g.(Rahmani et al., 2014), to pursue analysis of the use of magnetic markers for mappinghydraulic fractures using electromagnetic techniques.Magnetic permeability of common materials tends to vary over an order of magni-tude (Telford et al., 1990), while electrical conductivity of common materials can varyover > 8 orders of magnitude. Comparatively, there are many more materials that areelectrically conductive than there are with a significant magnetic permeability. Materialssuch as coke breeze, a hardened graphite coating applied to conventional proppants, havebeen considered by numerous authors, as have a range of manufactured electrically con-ductive proppants (Pardo and Torres-Verdin, 2013; Hoversten et al., 2015; Weiss et al.,2015; Labrecque et al., 2016; Hu et al., 2018). For these reasons, I focus on electricallyconductive proppants and treat the fractured region of the reservoir as an electricallyconductive geophysical target.2.1.2 Numerical ModellingNumerical modelling is a critical component for assessing the feasibility of detecting afractured volume of rock with an electromagnetic survey. To run a simulation, we needto discretize the modelling domain and represent the electrical conductivity of the earthon the simulation mesh. It is important to consider the large range of scales at play whenconsidering fractures. The proppant which fills the fractures is micro-to-millimeters indiameter. The fractures are millimeters thick and we are aiming to characterize a re-gion of the reservoir which extends hundreds of meters from the injection point at thewell, tens to hundreds of meters in height, and tens of meters along the length of thewell-bore. Furthermore, the numerical simulation domain must extend far enough to34Figure 2.1: Fracture complexity, from simple planar fractures on the left to com-plex fracture networks on the right. The blue dot is the injection point in ahorizontal well. After Cipolla et al. (2008a).satisfy boundary conditions (typically that the fields have sufficiently decayed). Thisproblem is exacerbated when one considers an inversion, where many forward simula-tions must be performed. We cannot expect a method to be capable of imaging bothsuch a substantial volume of the subsurface while also having the resolution on the scaleof the proppant particles. Thus, we require a characterization of the bulk impact of theconductive proppant within a fractured volume of rock. How we construct such a bulkphysical property model depends on the geometry and complexity of the induced frac-tures, for instance, different approaches should be considered if the fracture is a simpleplanar fracture versus a complex fracture network, as depicted in Figure 2.1 (c.f. Cipollaet al. (2008b)).To overcome these challenges, there are several approaches that may be taken; theappropriate choice will depend upon both the fracture complexity and the purpose of thesimulation. If the fine-scale geometry of the fractures is defined, for example in a feasi-bility study built from a synthetic fracture model, then numerical upscaling (Durlofsky,2003; Caudillo-Mata et al., 2014; Caudillo-mata et al., 2016) or multiscale techniques(Haber and Ruthotto, 2018) can be employed. In general, though, the fine-scale fracturegeometry is not known; only a handful of studies have “ground-truthed” the geometryof an induced fracture by mining the fractured region (Cipolla et al., 2008a). Further-35more, in an inversion, we cannot expect to resolve individual fractures. Rather, the aimis to characterize the bulk impact due to electrically conductive fractures in a conductivemedium. To estimate a bulk conductivity, effective medium theory is one approach thatcan be taken (c.f. Torquato (2002); Milton (2002); Berryman and Hoversten (2013)).These methods assume that the composite material is composed of a collection of ran-domly distributed spheres or ellipsoids (which may be preferentially aligned). The esti-mated conductivity depends upon the conductivity of each of the materials, their shape,and volume in the composite material. Looking forward to the inverse problem, thesemethods have the added benefit that they provide a conduit for incorporating a-priori in-formation such as expected primary-fracture orientation, and volumes of proppant andfluid.In this work, I adopt an effective medium theory approach. I take into account theconductivity of the host rock, the fracturing fluid and the proppant. The main simpli-fying assumption I employ is on the geometry of the fractures: I assume that they arecomposed of a collection of ellipsoidal cracks which may be preferentially or randomlyaligned. The derivation I present follows similar derivations presented in the literature(in particular, Torquato (2002)), here I include details for preferentially aligned fracturesresulting in a fully anisotropic conductivity.2.1.3 Chapter overviewThe purpose of this chapter is twofold: (1) I develop a workflow for estimating thephysical properties of a fractured volume of rock containing proppant and fluid thathave distinct electromagnetic properties from the host and (2) I assess if, using a com-mercially available electromagnetic survey, we can detect the anomaly introduced bythe fracture in the data. All of the computations shown in this chapter are open-source36and available as Jupyter notebooks (see Appendix A).2.2 Homogenization workflow using effective mediumtheoryIn this approach, I treat the fractures as being composed of a collection of preferentially(or randomly) oriented ellipsoidal cracks and based on the density of cracks within agiven volume of rock, construct an anisotropic description of the coarse-scale conduc-tivity.Effective medium approximations range from applying simple harmonic or arith-metic averaging of conductivity values, for example when constructing a representativevoxel conductivity model from well-log measurements, to more involved analytic orempirical relationships such as Archie’s law (Archie, 1942), which is commonly ap-plied for estimating the conductivity of a fluid-filled rock. Although discovered em-pirically, the simplest version of Archie’s law is one example of a differential effectivemedium approximation which can be derived analytically. It assumes a backgroundmatrix, and uses an incremental approach to constructing a homogenized conductivity(c.f. Torquato (2002); Milton (2002)). However, for describing a fractured volume ofrock, differential effective medium approximations are not appropriate as they assumethat the rock-matrix is always connected (Torquato, 2002); in the composite material weare considering, a single computational voxel may be intersected by a fracture, mean-ing the rock matrix is not connected. The Maxwell-approximation (Maxwell, 1873) isyet another effective-medium approximation. It again makes a distinction between thebackground and the included phases, and assumes no interaction between inclusions.I opt to consider self-consistent effective medium theory (SCEMT, also sometimesreferred to as the Coherent Potential Approximation, CPS, or Bruggeman mixing). This37is an effective medium approach which makes no distinction between background andincluded phases (e.g. see Torquato (2002)). Berryman and Hoversten (2013) similarlysuggest using self-consistent effective medium theory for fractured rocks where the frac-tures are filled with fluid.Milton (1985) demonstrated the physical realizability of SCEMT; he showed thatSCEMT is asymptotically the exact solution for the effective conductivity for fractal-like composites that are self-similar at many scales. Physical realizability means thatthe estimates given by SCEMT will always be within the Hashin-Shtrikman bounds(Torquato, 2002); at a minimum, this provides confidence that the conductivity estimatesit provides are physically within a reasonable realm.I construct the physical property model for a fractured volume of rock in two steps,as shown in Figure 2.2. Given the conductivity of the fluid and proppant particles, Iestimate the effective conductivity, σ2, of a mixture of proppant and fluid. Next, I treatthe fracture as consisting of a collection of ellipsoidal cracks filled with the proppant-fluid mixture. The ellipsoidal cracks may be preferentially aligned in a single or multipledirections or they may be randomly oriented. In both cases, I use the self-consistenteffective medium theory, originally due to Bruggeman (1935) and further developed byLandauer (1952, 1978). In the following sections, I develop the theory and demonstrateits application for computing the effective conductivity of a proppant-fluid mixture aswell as a fractured volume of rock.2.2.1 Summary of self-consistent effective medium theoryChapter 18 of Torquato (2002) provides an overview of effective medium theory ap-proaches. The discussion presented in this section follows their presentation, but ad-ditionally considers preferentially aligned cracks, resulting in an anisotropic effective38Figure 2.2: Constructing a physical property model for a fractured volume of rockusing effective medium theory. The electrical conductivity of the proppantfluid mixture is given by σ2, the conductivity of the background reservoirrock by σ1. Using effective medium theory, the coarse-scale anisotropicconductivity, Σ∗ describing the fractured volume of rock is computed.conductivity. Berryman and Hoversten (2013) present a similar overview but introducesseveral simplifying assumptions tailored for estimating the effective conductivity of anaturally fractured rock where the fractures are a relatively small concentration withrespect to the host rock. Here, I avoid such assumptions and work with the general for-mulation with the aim of being suitable for calculating both the effective conductivityof a proppant-fluid mixture as well as arbitrarily oriented fractures, each at potentiallyhigh concentrations.Each material, or phase, in the composite is assumed to be made up of spherical orellipsoidal particles with a known aspect ratio. Starting from the solution for a sphereor an ellipsoid in a uniform electric field, the effective conductivity of a heterogeneousmedium is chosen to be the conductivity for which the average perturbation to the elec-tric field – the difference between the electric field in the homogenized medium and the39true conductivity model – is zero. That is,N∑j=1ϕ j(Σ∗−σ jI)R( j,∗) = 0 (2.1)where N is the number of different phases of materials, ϕ j is the volume fraction of thej-th phase, and σ j is the electrical conductivity of the j-th phase. Σ∗ is the 3×3 effectiveconductivity tensor, and I is the 3× 3 identity matrix. The matrix R( j,∗) is the electricfield concentration tensor, and depends both on the shape of the inclusions (ie. proppantparticles or cracks composing a fracture) and conductivity of the j-th phase, as well asthe effective conductivity Σ∗.For spherical particles, the electric field concentration tensor R( j,∗) reduces to ascalar, namely,R( j,∗) =[I+13Σ∗−1(σ jI−Σ∗)]−1(2.2)The resultant effective conductivity expression in 2.1 then reduces to a scalar equation.If instead ellipsoidal inclusions are considered, the electric field concentration tensoris given byR( j,∗) =[I+AΣ∗−1(σ jI−Σ∗)]−1(2.3)Where A is the de-polarization tensor. For simplicity, I will assume that we are workingwith spheroids, either oblate (pancake-like) or prolate (needle-like) spheroids whichhave two semi-axes that are equal. The general solution for spheroids with three distinctsemi-axes is presented in chapter 17 of Torquato (2002) and additionally discussed inBerryman and Hoversten (2013). For a spheroid with semi-axes a1 = a2 = a and a3 = b40that is aligned so that a3 lies along the z-axis, the depolarization tensor is given byA =QQ1−2Q (2.4)For prolate spheroids, with aspect ratio α = b/a> 1, Q is given byQ=12(1+1α2−1[1− 12χbln(1+χb1−χb)])(2.5)and for oblate spheroids (α = b/a< 1),Q=12(1+1α2−1[1− 1χatan−1 (χa)])(2.6)withχ2a =−χ2b =1α2−1 (2.7)For preferentially aligned spheroids, the matrix A can be rotated as to align with thespheroidal axis using standard coordinate rotations. If the spheroids are randomly ori-ented, then we replace R( j,∗) with 1/3trace(R( j,∗)) in equation 2.1, and the effectiveconductivity expression reduces to a scalar equation. Note, that for spherical inclusions,Q= 1/3 and thus A reduces to a scalar equal to 1/3, showing that 2.3 is consistent with2.2.To solve for the effective conductivity, which is an implicit expression for Σ∗, werearrange equation 2.1 toΣ∗N∑j=0ϕ jR( j,∗) =N∑j=0ϕ jσ jR( j,∗) (2.8)41and solveΣ∗ =N∑j=0ϕ jσ jR( j,∗)(N∑j=0ϕ jR( j,∗))−1(2.9)using fixed-point iteration as R( j,∗) depends on Σ∗. Gradient-based optimization tech-niques could additionally be considered to speed-up convergence, but as this is a simple3×3 matrix equation and the fixed-point iteration is sufficiently fast for the problems Iconsider. Note that the matrix inverse in 2.9 is a 3×3 matrix inverse and thus is cheap toexplicitly compute. The fixed-point iteration is performed until the recovered effectiveconductivity converges within a predefined tolerance. Berryman and Hoversten (2013)similarly use a fixed-point iteration to solve for the effective conductivity, however intheir formulation, they isolate Σ∗ by pulling out the first term in the summation in 2.1,leading to an update of the formΣ∗ = σ0I− 1ϕ0 R(0,∗)−1N∑j=1ϕ j(Σ∗−σ jI)R( j,∗)This is equivalent to equation 10 in Berryman and Hoversten (2013) with the sim-plifications that R( j,∗) is replaced by R( j,0) under the assumption that the fractures com-pose a small volume fraction of the fractured rock. In practice, this approach is suitablefor low concentrations of included phases, but can cause instability in the algorithm athigher concentrations (it is possible for updates to be negative, and therefore unphysi-cal); the authors noted challenges with algorithm convergence when the concentrationof inclusions exceeded 0.2.In the following sections, I use this formulation to estimate the effective conductivityof a range of proppant-fluid mixtures as well as for a fractured volume of rock.422.2.2 Step 1: Effective conductivity of a proppant-fluid mixtureIn general, an induced fracture will be filled with two types or phases of material: prop-pant and fluid. Often this mixture is referred to as a slurry. There are two principal typesof mixtures I will consider, one where all of the proppant is conductive, and the secondwhen there is a conductive filler added to a conventional proppant.I start by considering the case of a proppant with uniform conductivity, for instanceif a conductive proppant, or coated proppant were used. I assume the conductivity ofthe proppant is known; Appendix B includes a derivation for an effective conductivityof two concentric spheres for the scenario where the conductivity of a proppant particleand its coating are known independently.For spherical proppant particles, the tensor-values in equation 2.1 reduces to a scalar-valued equation, and the resulting effective conductivity is isotropic. In Figure 2.3, Ishow the effective conductivity found using equations 2.1 and 2.2 for a proppant-fluidmixture as the concentration of proppant in the mixture (ϕ) is varied. The conductivityof the fluid is 3 S/m (similar to that of sea-water), and the proppant conductivity is variedlogarithmically from 10 S/m to 105 S/m. For example, coke-breeze, a graphite basedmaterial has conductivities ∼ 3000 S/m, and other authors have considered the use ofcontrast agents that reach conductivities of 105 S/m Weiss et al. (2016) and 106 S/mPardo and Torres-Verdin (2013), which are similar to the conductivity of steel.When the volume fraction of proppant is less than 1/3, the conductivity of the fluidis the dominant control on the resulting effective conductivity. Above a volume fractionof 1/3, the conductivity of the proppant is the primary contributor to the effective con-ductivity. The threshold between these behaviors is the percolation threshold. Belowit, the concentration of conductive material is low enough that it is quite likely discon-nected, above 1/3, the concentration is high enough to start forming connected, elec-43Figure 2.3: Effective conductivity of a proppant-fluid mixture for five differentproppant conductivities, each indicated in the legend. Panel (a) shows theconductivity on a linear scale, and panel (b) uses a log-scale for the conduc-tivity. Both panels use a linear scale for the volume fraction of the proppant(ϕ). The conductivity of the fluid is 3 S/m, similar to the conductivity ofsea-water.trically conductive pathways, causing a large jump in the effective conductivity of thesystem. Although proppant typically composes 10% to 20% of the injected slurry, someof the injected fluid leaks off into the surrounding geologic formation leaving proppantconcentration that can be 50% in the fractures Novotny (1977); Hoversten et al. (2015).The conductivity of the fluid also changes the resultant effective conductivity of themixture. In Figure 2.4, I compare the effective conductivity for a mixture of conductiveproppant (panel (a): 103 S/m, panel (b): 104 S/m) and four different fluid conductivities,ranging from 0.3 S/m to 300 S/m, as indicated in the legend. Although the conductivity44Figure 2.4: Impact of the conductivity of the fluid on the effective conductivity of aproppant-fluid mixture. Panels (a) and (c) show the effective conductivity formixtures with a 103 S/m proppant and panels (b) and (d) show the effectiveconductivity for mixtures with a 104 S/m proppant. The conductivity of thefluid is indicated by the legend.of the fluid makes a significant difference at low proppant concentrations, above thepercolation threshold of 33%, the curves start to converge, particularly when the contrastbetween the conductivity of the proppant and the fluid exceeds 3 orders of magnitude.Thus, if the proppant can be made sufficiently conductive, its conductivity will be thecontrolling factor on the effective conductivity of the slurry that remains in the fractures.The previous examples considered a 2-phase mixture in which all of the proppantwas electrically conductive, however, depending on the setting and the cost to man-ufacture conductive proppant, it may be mixed in with conventional, resistive prop-pant. To examine this, I consider a proppant-fluid mixture composed of three materials,45fluid (3 S/m), conventional, resistive proppant (10−6 S/m) and conductive proppant (105S/m). The effective conductivities of proppant-fluid mixtures for five different proppantblends, where the relative concentration of the conductive proppant is varied from 0%to 100% of the proppant phase, are shown in Figure 2.5. Again, we see the impactsof the percolation threshold; when the conductive proppant composes less than 1/3 ofthe proppant pack, the effective conductivity is dominated by the resistive proppant.When the conductive proppant composes more than 1/3 of the proppant pack, we seethat with increasing proppant concentration, the effective conductivity of the mixture isdominated by the conductive proppant. However, the percolation threshold for each ofthese mixtures is different. This is because it is the volume-fraction of the conductiveproppant in the three-phase mixture, not the ratio of proppant to fluid, that determineswhen connected, electrically conductive pathways may be formed.Another factor influencing the effective conductivity of a mixture is the shape ofthe materials. For the previous examples, the proppant was assumed to be composedof spherical particles. If elongated, conductive particles were included, we expect thatconnected, conductive pathways would form at lower concentrations. For instance, con-sider a 3 phase proppant mixture consisting of fluid (3 S/m), resistive, spherical proppant(10−6 S/m), and elongated, electrically conductive proppant (105 S/m). Assume that theratio of conductive to resistive proppant is 0.25 (below the percolation threshold forspherical particles). If the elongated particles (prolate spheroids) are randomly oriented,then the resulting effective conductivity is isotropic, meaning it is independent of thedirections of the inducing field and resulting current. The conductivity predicted by ef-fective medium theory for mixtures with five different aspect ratios is shown in figure2.6. The aspect ratio of the conductive particles influences the concentration at whichwe observe percolation. The more elongated the particles, the lower the concentration46Figure 2.5: Effective conductivity of a 3-phase proppant-fluid mixture consistingof resistive proppant (10−6 S/m), conductive proppant (105 S/m) and salinefluid (3 S/m). The legend indicates the percentage of conductive proppant inthe proppant mixture and the x-axis is the volume fraction of proppant in theproppant-fluid mixture.at which percolation occurs.In summary, there are several approaches that can be taken to create an electricallyconductive proppant-fluid mixture. Spherical proppant particles which are themselveselectrically conductive or coated with a conductive material can comprise the entireproppant pack. If conductive proppant is mixed with a conventional, resistive sand orceramic particles, then at least 1/3 of the proppant mixture needs to be comprised ofelectrically conductive particles to create a conductive mixture. This ratio can be re-duced if elongated particles, such as metallic strips, are used in the proppant mixture.47Figure 2.6: Effective conductivity of a 3-phase proppant-fluid mixture consistingof resistive proppant (10−6 S/m), conductive proppant (105 S/m) and salinefluid (3 S/m). The proppant mixture contains 25% conductive proppant and75% resistive proppant. The legend indicates the aspect ratio of the elongatedconductive proppant filler (prolate spheroids).2.2.3 Step 2: Effective conductivity of fractured volume of rockThe next step is to estimate the effective conductivity of a fractured volume of rock.I again employ self-consistent effective medium theory as described in section 2.2.1and consider the induced fractures to be composed of spheroidal cracks. Based on theanalysis in the previous section, I consider a proppant-fluid mixture that has a conduc-tivity of 2500 S/m. This could be achieved with spherical proppant particles having aconductivity of 104 S/m in a 50/50 mixture with water of 3 S/m (see Figure 2.3). Sim-ilar conductivities could be achieved with elongated particles mixed with conventional48proppant as shown in Figure 2.6. Lab measurements conducted by Zhang et al. (2016)showed that a proppant-fluid mixture composed of petroleum coke particles and seawa-ter which fills the pore-spaces reached a conductivity of ∼ 1000 S/m at 37.6% porosity.With further increase in confining pressure (thus reducing porosity and increasing theconcentration of proppant), the observed conductivities range from 3000− 5000 S/m.These results provide further confidence that conductivities > 1000 S/m for the mixturefilling the hydraulic fractures are attainable.For the following example, I will assume that the conductivity of the host-rock is0.1 S/m. There are two scenarios I will consider. In the first, I assume the cracks arepreferentially aligned, with the thin dimension of the oblate spheroid oriented along thex-axis, as shown in Figure 2.7. In this case, the recovered effective conductivity will beanisotropic, described by a diagonal matrix with entries σy = σz ≥ σx :Σ∗ =σxσyσz (2.10)Note that arbitrary, non-axes aligned, orientations can be considered; all that is requiredis that the depolarization tensor described in 2.4 is rotated to the desired orientation.To estimate the effective conductivity of a fractured volume of rock, we must alsospecify the aspect ratio of the cracks. In estimating this, assume a fractal-like approxi-mation, where the aspect ratio of the fracture is representative of the aspect ratio of thecracks that compose it. For example, if the fracture extends 50m laterally and has awidth on the order of millimeters, then the aspect ratio is on the order of 10−5. In Figure2.8, I have plotted the diagonal elements of the effective conductivity for five differentaspect ratios, indicated in the legend, as a function of the volume fraction of conductive49Figure 2.7: Oblate spheroid with normal along the x-axis. The y and z semi-axesare equal, and aspect ratio is α = b/a< 1.fractures in the rock volume sampled. Panel (a) shows the full range from 0 ≤ ϕ ≤ 1and panel (b) zooms in to lower concentrations (0 ≤ ϕ ≤ 0.01) which are more repre-sentative of a fractured rock volume, on the scale that we will consider for numericalmodelling (e.g. if 10 fractures, each with 3mm width intersect a 10m × 10m × 10mcell, then ϕ = 0.003). In each of the plots, I have also included the upper and lowerWiener bounds (see equation 21.14 in Torquato (2002); originally attributed to Wiener(1912)):σ+W =N∑j=0ϕ jσ jσ−W =(N∑j=0ϕ jσ j)−1 (2.11)in the black dashed lines; these can be understood as similar to parallel and series circuit50approximations to the conductivity. In the black dash-dot lines are the upper and lowerHashin-Shtrikman bounds for 2-phase anisotropic media in the black dash-dot (see equa-tions 21.25 and 21.26 in Torquato (2002), which is the anisotropic generalization of theisotropic bounds derived by Hashin and Shtrikman (1962)):σ+HS = σ+W I+(σ1−σ0)2A˜ ·[σ1I+σ1−σ0ϕ0A˜]−1σ−HS = σ+W I+(σ1−σ0)2A˜ ·[σ0I+σ0−σ1ϕ1A˜]−1 (2.12)For σ1 ≥ σ0 andA˜ =−ϕ0ϕ1A1 (2.13)where A1 is the depolarization tensor for the ellipsoidal cracks given by equation 2.4.For the bounds shown in the plot, the smallest aspect ratio, 10−5 was used to calculatethe depolarization tensor. For the very significant aspect ratios used here, the Hashin-Shtrikman bounds are nearly identical to the Wiener bounds. Each component of therecovered effective conductivity should fall within these bounds.For aspect ratios less than 10−3, we see very little distinction between the estimate ofσy,z and σx; the difference between the recovered effective conductivities for the aspectratios of 10−4 and 10−5 is less than 1% for all values of ϕ . This indicates that for suf-ficiently thin cracks, the exact aspect ratio is not critical for estimating a representativeconductivity of the fractured rock.At the concentrations we expect to be observing in a hydraulic fracturing scenario(Figure 2.8b), we see that the effective conductivity along the normal of the fracturescoincides with the lower bounds and remains nearly identical to the conductivity of thehost rock for all aspect ratios shown, as may be expected. For the components of theconductivity along the cracks (σy,z), the effective conductivity mimics the behavior of51Figure 2.8: Effective, anisotropic conductivity for a fractured rock with spheroidalcracks whose normal is oriented along the x-axis for five different aspect ra-tios, indicated by the legend. The black dashed lines show the upper andlower Wiener Bounds, which are identical to the volume-weighted arith-metic and harmonic averages of the conductivity of the rock (0.1 S/m)and proppant-fluid mixture (2500 S/m). The black dash-dot lines are theanisotropic Hashin-Shtrikman upper and lower bounds computed using anaspect ratio of 10−5 in equation 2.12. Note that the upper-bound coencideswith σy,z for small aspect ratios (e.g. the blue line). Panels (a) and (c) showthe full range 0 ≤ ϕ ≤ 1, and panels (b) and (d) zoom in to 0 ≤ ϕ ≤ 0.005.Note that in (b), any curve that deviates from the lower bound is the σy,zcomponent; all σx values lie on the lower bound for the concentration rangeshown.52the upper bounds. These components of the conductivity are similar to the behaviorexpected from a parallel-circuit approximation.For settings where planar fractures are expected, a single orientation of the inclu-sions produces similar behavior to what would be expected if we performed simpleseries-and-parallel circuit approximations to the components perpendicular to and alongthe fracture. However, if more complex fracture networks are expected, it may be moreappropriate to assume the cracks are randomly oriented. In this case, an isotropic con-ductivity describes the fractured volume of rock. Using the same aspect ratios shown inFigure 2.8 above, I compute the isotropic effective conductivity as a function of volumefraction of fractures in Figure 2.9. For reference, the conductivity along the fractures,σy,z from Figure 2.8, is plotted in semi-transparent lines. As in the previous figure, panel(a) shows the full range from 0≤ ϕ ≤ 1 and panel (b) zooms in to lower concentrations:0≤ ϕ ≤ 0.01.Again, we see that for aspect ratios smaller than 10−3 there is little difference be-tween the estimated effective conductivity of the fractured rock volume. The estimateof the effective conductivity largely follows the behavior of the upper bounds, while themagnitude of the effective conductivity is slightly reduced as compared to the compo-nent of the conductivity along the fractures in the anisotropic case. If we consider a10 m × 10 m × 10 m computational cell with fractures having an aspect ratio of 10−5and composing 0.3% of the total volume (e.g. 10 fractures, each with 3 mm width), weobtain σy = σz = 5.2 S/m, σx = 0.1 S/m for the case where the cracks are preferentiallyaligned while for the case where the cracks are randomly oriented, we obtain σ∗ = 3.5S/m.Due to the large contrast between the conductive fractures and the host rock, theconductivity of the background has marginal effect on σy,z in the anisotropic example53Figure 2.9: Effective, isotropic conductivity for a fractured rock with randomlyoriented spheroidal cracks for five different aspect ratios, indicated by thelegend. The black dashed lines show the upper and lower Wiener Bounds,which are identical the volume-weighted arithmetic and harmonic averagesof the conductivity of the rock (0.1 S/m) and proppant-fluid mixture (2500S/m).The black dash-dot lines are the isotropic Hashin-Shtrikman upper andlower bounds. The semi-transparent dotted lines show σy,z from Figure 2.8Panels (a) and (c) shows the full range 0 ≤ ϕ ≤ 1, and panels (b) and (d)zooms in to 0≤ ϕ ≤ 0.01.or on σ∗. If we consider the example with ϕ = 0.003, as before, and instead use abackground conductivity of 0.01 S/m, the effective anisotropic conductivity is σy =σz = 5.1 S/m, σx = 0.01 S/m and for the isotropic case σ∗ = 3.4 S/m. However, if thebackground is made more conductive and the contrast between the conductive fracturesand the host reduced, we do see some impact. Setting the background to 1 S/m for thissame example, we obtain σy = σz = 6.4 S/m, σx = 1 S/m for the anisotropic case and54σ∗ = 4.7 S/m for the isotropic case.2.2.4 SummaryTo construct an approximate conductivity model of a fractured volume of rock, I usea two step process: first, I estimate the effective conductivity of a mixture of proppantand fluid that fills the fractures and second, I estimate the effective conductivity of afractured volume of rock.In the section that follows, I will compare numerical simulations of a crosswell elec-tromagnetic survey for both anisotropic and isotropic conductivity models of a fracturedvolume of rock, representing planar and complex fracture networks, respectively.2.3 Feasibility of EM for detecting fracturesI now examine the feasibility of detecting a conductive, fractured region of a reservoirwith an electromagnetic survey. The purpose of this section is to demonstrate detectabil-ity for a simple model of a fracture with existing electromagnetic imaging techniques.The survey I consider is a frequency-domain crosswell electromagnetic survey. In acrosswell EM survey a transmitter coil, which produces a time-varying magnetic field ata single frequency is positioned inside of one well. Time-varying magnetic fields inducetime-varying currents in the earth whose distribution and magnitude depends upon theconductivity of the earth. These currents in turn generate secondary magnetic fields. In asecond well, an array of receiver induction coils, which measure the time rate-of-changeof the magnetic flux (db/dt, in units of volts), are positioned. The measured voltagescan be converted to magnetic field values. To conduct a survey, the position of the re-ceivers is fixed and the transmitter is moved at a fixed rate along the borehole (typically3-5 meters per minute). The receivers are repositioned, and the process repeated (Wilt55et al., 1995). Both the transmitter and receivers are oriented along the axis of the bore-hole; for vertical wells, this means the transmitter is a vertical magnetic dipole, and thereceivers measure the z-component of the magnetic field. Since its development, cross-well surveys have been used for a range of reservoir imaging problems for enhancedoil recovery including monitoring steam injections (Wilt et al., 1997) and water floods(Wilt et al., 2012), including a survey conducted in a pair of horizontal wells (Marsalaet al., 2015; Marsala* et al., 2015). In this section, I consider a crosswell survey andexamine the feasibility of detecting signal due to a fractured volume of rock.The setup I use is shown in Figure 2.10; a crosswell EM survey conducted in apair of horizontal wells that are spaced 500 m apart, as shown in Figure 2.10. Thetransmitter is oriented along the x-axis and positioned in the same well as where thefracture stimulation was performed; the receivers are in the offset well and measure thereal and imaginary parts of the x-component of the magnetic field. In the numericalmodelling, I assume that the wells are sufficiently deep so that the background canbe treated as a whole-space. I do not consider steel-cased boreholes in this exampleand instead assume that the boreholes are open-holes or cased with fibreglass casingas in Wilt et al. (2012). The additional complications due to steel cased wells will bediscussed in the following three chapters.Most hydraulic fracture operations involve injecting ∼800 m3 to 4000 m3 (5000 to25000 bbls) of proppant-fluid slurry into the reservoir rock, with the proppant compris-ing 10% to 20% of this mixture (Hoversten et al., 2015). In this example, I consider afracture that is on the lower-end, with an 800 m3 slurry comprised of 15% proppant byvolume. A significant portion of the fluid typically leaks off into the surrounding forma-tion (Novotny, 1977), so I will assume that the final mixture filling the cracks contains50% fluid and 50% proppant. This gives a total fracture volume of 240 m3. To distribute56Figure 2.10: Setup for the crosswell electromagnetic simulation. Ten fractures,each with 3 mm thickness and 50 m radius are spaced over a 10 m intervalalong the injection well. The dipole transmitter is positioned inside of theinjection well and oriented along the x-axis. Receivers are located in anoffset well and measure the real and imaginary parts of the x-component ofthe magnetic field.this volume, I assume that 10 fractures, each with 3 mm width are positioned within a 10m interval along the well. Conserving volume and assuming circular fractures gives us a50 m radius for each of the 10 fractures; I work with circular fractures that are centeredalong the borehole, as the resultant model is cylindrically symmetric. This significantlyreduces the computational cost of the simulation.For the physical properties, I treat the background as a whole-space with a conduc-tivity of 0.1 S/m. Within the fractures, I use a two-phase proppant-fluid mixture, whereall of the proppant is assumed to be spherical particles as in Figure 2.3. I use a prop-pant conductivity of 104 S/m and a fluid conductivity of 3 S/m; in a 50-50 mixture, the57effective conductivity computed using SCEMT, is 2500 S/m. To estimate the effectiveconductivity of the fractured volume of rock, I consider two scenarios: (1) the fracturesare preferentially aligned and I estimate an anisotropic effective conductivity as in Fig-ure 2.8; and (2) the fractures are sufficiently complex that I can approximate them as acollection of randomly oriented cracks, as in Figure 2.9. I estimate a single conductivityfor the 10 m wide cylinder with 50 m radius that encloses the fractures. Within thisregion, the fractures comprise 0.3% of the total volume. To estimate the aspect ratio, Iassume that the shape of the large fracture, with 3 mm width and 50 m radius (100 m di-ameter) is representative of the shape of the cracks that comprise it, thus the aspect ratiois 3×10−5. The estimated effective conductivity for the anisotropic case is σy = σz = 5S/m parallel to the fractures and σx = 0.1 S/m perpendicular to the fractures. For thesecond scenario, with randomly oriented cracks, the estimated isotropic conductivity is3 S/m.In terms of survey parameters, the dipole-moment of the transmitter is 5000 Am2and I use a transmitter frequency of 100 Hz, consistent with values stated in Wilt et al.(1995); Wilt (2003). For a signal to be detectable, we require that: (a) the secondarysignal, that is, the difference between the data measured with and without the fracturespresent, be above the noise floor of the receivers, and (b) that the secondary signal com-prises a significant percentage of the primary (model without the fractures). Marsala*et al. (2015) demonstrated noise levels of 0.5×10−5 nT and 0.5% repeatability betweendata profiles for a cross-well survey conducted between two open-hole horizontal wells,citing that “this was the lowest noise ever recorded with a crosswell EM system.” Usingthis as an upper-bound on our expectations, I adopt a more conservative noise level of1× 10−4 nT and select a percent-threshold of 5%. Thus, I consider signal above 10−4nT and comprising > 5% of the primary as detectable.58Figure 2.11 shows the simulated secondary magnetic field data computed for (a) theanisotropic model and (b) the isotropic model at 100 Hz. The top row shows the realcomponent and the bottom row shows the imaginary component. The x-axis shows thetransmitter location relative to the center of the fracture, the y-axis shows the receiverlocation, and the color indicates the secondary magnetic field with respect to a 0.1 S/mwhole-space. Regions of the plot which are colored, and are not covered by the semi-transparent overlay, indicate signal which is both above the noise floor and comprises asignificant percentage of the primary.As might be expected, when the transmitter and receiver are centered with the frac-ture, the secondary signal reaches its maximum amplitude in both the real and imag-inary components. The anisotropic model produces stronger signals; the geometry ofthe survey is perfectly coupled to the radial components of the conductivity (σy andσz) as the resultant currents and electric fields are purely azimuthal in orientation. Anisotropic simulation using σ = 5 S/m produces identical results to the anisotropic modelfor this survey geometry. Interestingly, the region over which data is measurable differsbetween the real and imaginary components. The real component is measurable overa larger range of transmitter positions, while the imaginary component is measurableover a larger range of receiver locations. The partition of EM energy into the real andimaginary components is a function of the conductivity of the medium as well as thefrequency. For example, if I increase the frequency to 200 Hz, the amplitude of theimaginary component is larger than the real, and it is measurable over a larger set ofsource-receiver combinations, as shown in Figure 2.12.For the simple setup shown here, there are a significant number of data that can becollected in a cross-well survey which are sensitive to both the anisotropic and isotropicfracture models. In both the 100 Hz and 200 Hz surveys, the extent over which either or59Figure 2.11: Simulated secondary magnetic field data for: (a) the anisotropic frac-ture model and (b) the isotropic fracture model, in a crosswell survey con-ducted at 100Hz. The top row shows the real component and the bottomrow shows the imaginary component of the measured magnetic flux (nT).In all plots, the x-axis shows the transmitter location relative to the center ofthe fracture, the y-axis shows the receiver location, and the color indicatesthe secondary magnetic flux with respect to a 0.1 S/m whole-space. Valuesbeneath the noise floor of 10−4 nT have been masked and display as white.To display the signal as a percentage of the primary, we have included asemi-transparent overlay between 0% and 5%; low percentage values plotas darker greys.60Figure 2.12: Simulated secondary magnetic field data for: (a) the anisotropic frac-ture model and (b) the isotropic fracture model, in a crosswell survey con-ducted at 200Hz. The top row shows the real component and the bottomrow shows the imaginary component of the measured magnetic flux (nT),similar to Figure 2.1161both of the real and imaginary components of the data are above the noise floor and com-prise a significant percentage of the primary field spanned ∼100 m for the transmitterlocation and > 800 m for the receiver locations. These results provide confidence thatEM is a potentially diagnostic tool for characterizing the propped region of a reservoir.Naturally, the success, or not, of EM in a given setting will depend upon the backgroundconductivity, volume and conductivity of the proppant and fluid, as well as noise con-ditions in the field. Thus, forward modelling tailored to the setting is required to assessfeasibility.2.4 ConclusionThe large range of scales involved in describing a propped, fractured volume of rockas a geophysical target for an electromagnetic survey poses a challenge for performingnumerical simulations. I approached this problem by using effective medium theoryto estimate the electrical conductivity of a fracture volume of rock in two steps. First,I estimated the conductivity of a slurry containing proppant and fluid. I demonstratedseveral approaches for creating an electrically conductive slurry, including using con-ductive particles (or particles coated with a conductive material) as the proppant or in-cluding conductive particles, which might be elongated to enhance the conductivity atlower concentrations. General agreement with the values published in the laboratoryexperiment by Zhang et al. (2016) gives confidence that the estimates produced usingeffective medium theory are representative of the conductivities that can be achieved; amore rigorous comparison would require control on the conductivity of each of the con-stituents. In the second step, I estimated the effective conductivity of a fractured volumeof rock by assuming that the fracture is composed of a collection of ellipsoidal cracks.These cracks may be preferentially oriented, if simple planar fractures are expected, or62randomly oriented if a more complex fracture network is anticipated.Using this approach, I developed isotropic and anisotropic electrical conductivitymodels considering proppant and fluid volumes that are representative of current hy-draulic fracture operations and simulated a crosswell EM survey. In this example, thefractured volume was detectable if I assume noise estimates that have demonstrated tobe achievable in horizontal wells for a waterflood monitoring experiment. This pro-vides evidence that EM has the potential to be a suitable method for characterizing thepropped region of a hydraulically fractured reservoir.The examples examined here neglect the effects of steel-cased well on EM signals.In the vast majority of settings where hydraulic fracturing is conducted, the wells arecased with steel, and therefore, this complication needs to be addressed. For inductionproblems, the concern is attenuation of the electromagnetic signals, as eddy currentsmay be induced in the steel cased well. However, there is increasing interest in usingthe steel casing as an “extended electrode” to help deliver current to depth. In the fol-lowing chapters, I develop a strategy for including steel cased wells in the EM numericalsimulations and examine their influence on EM signals.63Chapter 3Modelling electromagnetics oncylindrical meshesA number of geophysical electromagnetic (EM) problems lend themselves to cylindri-cal geometries. Airborne EM problems over a 1D layered earth or borehole-loggingapplications fall into this category; in these cases cylindrical modelling, which removesa degree of freedom in the azimuthal component, can be advantageous as it reducesthe computation load. This is useful when running an inversion where many forwardmodellings are required, and it is also valuable when exploring and building up an un-derstanding 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 reducesfeedback 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 footprintof the survey is primarily cylindrical, but 2D or 3D variations in the physical propertymodel may be present. For example if we consider a single sounding in an Airborne EM64survey, the primary electric fields are rotational and the magnetic fields are poloidal,but the physical property model may have lateral variations or compact targets. Moreflexibility is required from the discretization to capture these features. In this case, a 3Dcylindrical geometry, which incorporates azimuthal discretization may be advantageous.It allows finer discretization near the source where we have significant sensitivity andthe fields are changing rapidly. Far from the source, the discretization is coarser, butit still conforms to the primary behaviour of the EM fields and fluxes and captures therotational electric fields and poloidal magnetic flux.In other cases, the most significant physical property variations may conform to acylindrical geometry, for example in settings where vertical metallic well-casings arepresent, or in the emerging topic of using geophysics to “look ahead” of a tunnel boringmachine. Of particular interest to the thesis is understanding the behavior of electromag-netic fields and fluxes in the presence of steel-cased wells. In addition to the hydraulicfracturing application, understanding EM in settings with steel-cased wells is of interestacross a range of applications, from characterizing lithologic units with well-logs (Kauf-man, 1990; Kaufman and Wightman, 1993; Augustin et al., 1989a), to identifying ma-rine 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 captureand storage (Pardo and Torres-Verdin, 2013; Bo¨rner et al., 2015; Um et al., 2015; Weisset al., 2016; Hoversten et al., 2017; Zhang et al., 2018). Carbon steel, a material com-monly used for borehole casings, is highly electrically conductive (106−107 S/m) andhas a significant magnetic permeability (≥ 100 µ0) (Wu and Habashy, 1994); it there-fore can have a significant influence on electromagnetic signals. The large contrasts inphysical properties between the casing and the geologic features of interest, along withthe large range of scales that need to be considered to model both the millimeter-thick65casing walls while also capturing geologic features, provide interesting challenges andcontext for electromagnetics in cylindrical geometries.In this chapter, I introduce an approach and associated open-source software im-plementation for simulating Maxwell’s equations over conductive, permeable modelson 2D and 3D cylindrical meshes. The software is written in Python (Van Rossumand Drake Jr, 1995) and is included as an extension to the SimPEG ecosystem (Cock-ett et al., 2015); Appendix D provides a description of the SimPEG electromagneticmodule. Within the context of current research connected to steel-cased wells, my aimwith the development and distribution of this software is two-fold: (1) to facilitate theexploration of the physics of EM in these large-contrast settings, and (2) to provide asimulation tool that can be used for testing other EM codes. The large physical prop-erty contrasts in both conductivity and permeability means the physics is complicatedand often non-intuitive; as such, we prioritize the ability of the researcher to access andvisualize fields, fluxes, and charges in the simulation domain. This is particularly usefulwhen the software is used in conjunction with Jupyter notebooks which facilitate explo-ration of numerical results (Perez et al., 2015). As the mesh conforms to the geometryof a vertical borehole, a fine discretization can be used in its vicinity without resultingin a onerous computation. This provides the opportunity to build an understanding ofthe physics of EM in settings with vertical boreholes prior to moving to settings withdeviated and horizontal wells. I demonstrate the software with examples at DC, in thefrequency domain, and in the time domain. Source-code for all examples is providedas Jupyter notebooks as outlined in Appendix A; they are licensed under the permissiveMIT license with the hope of reducing the effort necessary by a researcher to compareto or build upon this work.This chapter is organized in the following manner. In section 3.1, I introduce the66governing equations, Maxwell’s equations, and describe their discretization in cylin-drical coordinates. I then compare our numerical implementation to the finite elementand finite difference results shown in (Commer et al., 2015) as well as a finite volumeOcTree simulation described in (Haber et al., 2007). Section 3.2 contains numerical ex-amples of the DC, frequency domain EM, and time domain EM implementations. Thetwo DC resistivity examples (sections 3.2.1 and 3.2.2) are built upon the foundationalwork in (Kaufman, 1990; Kaufman and Wightman, 1993) which use asymptotic analy-sis to draw conclusions about the behavior of the electric fields, currents, and chargesfor a well where an electrode has been positioned along its axis. The final two examples,in sections 3.2.3 and 3.2.4, consider a frequency domain experiment inspired by (Au-gustin et al., 1989a). These examples demonstrate the impact of magnetic permeabilityon the character of the magnetic flux within the vicinity of the borehole and discussesthe resulting magnetic field measurements made within a borehole.3.1 Numerical toolsThe governing equations under consideration are Maxwell’s equations. I provided anoverview in Section 1.5. This section is meant to provide a brief review and set thecontext for what is required of our numerical tools.Under the quasi-static approximation, Maxwell’s equations are given by:∇×~e=−∂~b∂ t∇×~h−~j =~se(3.1)where ~e is the electric field,~b is the magnetic flux density,~h is the magnetic field, ~j isthe current density and ~se is the source current density. Maxwell’s equations can also67be formulated in the frequency domain, using the eiωt Fourier Transform convention(equation 1.2), they are∇×~E =−iω~B∇× ~H− ~J =~Se(3.2)The fields and fluxes are related through the physical properties: electrical conductivity(σ , or its inverse, resistivity ρ) and magnetic permeability (µ), as described by theconstitutive relations~J = σ~E~B= µ~H(3.3)At the zero-frequency limit, we also consider the DC resistivity experiment, describedby∇ ·~j = I (δ (~r−~rs+)−δ (~r−~rs−))~e=−∇φ(3.4)where I is the magnitude of the source current density,~rs+ and~rs− are the locations ofthe current electrodes, and φ is the scalar electric potential.Of our numerical tools, we require the ability to simulate large electrical conductiv-ity contrasts, include magnetic permeability, and solve Maxwell’s equations at DC, infrequency and in time in a computationally tractable manner. Finite volume methodsare advantageous for modelling large physical property contrasts as they are conserva-tive and the operators “mimic” properties of the continuous operators, that is, the edgecurl operator is in the null space of the face divergence operator, and the nodal gradient68Figure 3.1: Anatomy of a finite volume cell in a (a) cartesian, regtangular mesh,(b) cylindrically symmetric mesh, and (c) a three dimensional cylindricalmesh.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. Horeshand Haber (2011); Haber and Ruthotto (2018); Jahandari and Farquharson (2014) andreferences within), and will be my method of choice.3.1.1 DiscretizationTo represent a set of partial differential equations on the mesh, I use a staggered-gridapproach (Yee, 1966) and discretize fields on edges, fluxes on faces, and physical prop-erties at cell centers, as shown in Figure 3.1. Scalar potentials can be discretized at cellcenters or nodes. I consider both cylindrically symmetric meshes and fully 3D cylindri-cal meshes; the anatomy of a finite volume cell for these scenarios is shown in Figure3.1 (b) and (c).To discretize Maxwell’s equations in the time domain (equation 3.1) or in the fre-quency domain (equation 3.2), I invoke the constitutive relations to formulate our system69in terms of a single field and a single flux. This gives a system in either the electric fieldand magnetic flux (E-B formulation), or the magnetic field and the current density (H-Jformulation). For example, in the frequency domain, the E-B formulation isCe+ iωb = smC>M fµ−1b−Meσe = se(3.5)and the H-J formulation isC>M fρ j+ iωMeµh = smCh− j = se(3.6)where e,b,h, j are vectors of the discrete EM fields and fluxes; sm and se are the discretemagnetic and electric source terms, respectively; C is the edge curl operator, and thematrices Me, fprop are the edge / face inner product matrices. In particular, variable electri-cal conductivity and variable magnetic permeability are captured in the discretization.The time domain equations are discretized in the same manner; for time-stepping, afirst-order backward Euler approach is used. Although the midpoint method, which issecond-order accurate, could be considered, it is susceptible to oscillations in the solu-tion, which reduce the order of accuracy, unless a sufficiently small time-step is used(Haber et al., 2004; Haber and Ruthotto, 2018). Appendix D provides more details onthe discretization of Maxwell’s equations in both the frequency and time-domains.At the zero-frequency limit, each formulation has a complementary discretizationfor the DC equations. For the E-B formulation the discretization leads to a nodal dis-70cretization of the electric potential φ , giving−G>Meσe = qe =−Gφ(3.7)where G is the nodal gradient operator, and q is the source term, defined on nodes.Note that the nodal gradient takes the discrete derivative of nodal variables, and thus theoutput is on edges. The H-J formulation leads naturally to a cell centered discretizationof the electric potentialVDj = qM fρ j = D>Vφ(3.8)Where D is the face divergence operator, V is a diagonal matrix of the cell volumes, qis the source term, which is defined at cell centers as is φ . Here, the face divergencetakes the discrete derivative from faces to cell centers, thus its transpose takes a variablefrom cell centers to faces. For a tutorial on the finite volume discretization of the DCequations, see (Cockett et al., 2016b).For the EM simulations, natural boundary conditions are employed; in the E-B for-mulation, this means ~B×~n = 0|∂Ω, and in the H-J formulation, we use ~J×~n = 0|∂Ω.Within the DC simulations, there is flexibility on the choice of boundary conditions em-ployed. In the simplest scenario, for the nodal discretization, we use Neumann bound-ary conditions, σ~E ·~n = 0|∂Ω, and for the cell centered discretization, we use Dirichletboundary conditions φ = 0|∂Ω.When employing a cylindrical mesh, the distinction between where the electricand magnetic contributions are discretized in each formulation has important implica-71tions. If we consider the cylindrically symmetric mesh (Figure 3.1b) and a magneticdipole source positioned along the axis of symmetry (sometimes referred to as the TEmode), we must use the E-B formulation of Maxwell’s equation to simulate the result-ing toroidal magnetic flux and rotational electric fields. If instead, a vertical currentdipole is positioned along the axis of symmetry (also referred to as the TM mode), thenthe H-J formulation of Maxwell’s equations must be used in order to simulate toroidalcurrents and rotational magnetic fields. The advantage of a fully 3D cylindrical mesh isthat it provides additional degrees of freedom, with the discretization in the azimuthaldirection. This allows us to simulate more complex responses. However, in order toavoid the need for very fine discretization in the azimuthal direction, we should selectthe most natural formulation of Maxwell’s equations given the source geometry beingconsidered. For a vertical steel cased well and a grounded source, we expect the major-ity of the currents to flow vertically and radially, thus the more natural discretization toemploy is the H-J formulation of Maxwell’s equations.Haber and Ruthotto (2018) provides derivations and discussion of the differentialoperators and inner product matrices; though they are described for a cartesian coor-dinate system and a rectangular grid, the extension to a three dimensional cylindricalmesh is straightforward. Effectively, a cartesian mesh is wrapped so that the x com-ponents become r components, and y components become θ components, as shown inFigure 3.2.The additional complications that are introduced are: (1) the periodic boundarycondition introduced on boundary faces and edges in the azimuthal direction, (2) theremoval of radial faces and azimuthal edges along the axis of symmetry, and (3) theelimination of the degrees of freedom of the nodes and edges at the boundary and aswell as the nodes and vertical edges along the axis of symmetry. The implementa-72Figure 3.2: Construction of a 3D cylindrical mesh from a cartesian mesh.tion of the 3D cylindrical mesh is provided as a part of the discretize package(http://discretize.simpeg.xyz), which is an open-source python package that contains fi-nite volume operators and utilities for a variety of mesh-types. All differential operatorsare tested for second order convergence and for preservation of mimetic properties (asdescribed in Haber and Ruthotto (2018)). discretize is developed in a modular,object-oriented manner and interfaces to all of the SimPEG forward modelling and in-version routines, thus, once the differential operators have been implemented, they canbe readily used to perform forward simulations (Cockett et al., 2015). One of the bene-fits of SimPEG for forward simulations is that values of the fields and fluxes are readilycomputed and visualized, which enables researchers to examine the physics as well asto simulate data. Development within the SimPEG ecosystem follows best practices formodern, opens-source software, including: peer review of code changes and additions,73versioning, automated testing, documentation, and issue tracking.3.1.2 ValidationTesting for the DC, TDEM, and FDEM implementations includes comparison with an-alytic solutions for a dipole in a whole-space. These examples are included as supple-mentary examples with the distributed notebooks. I have also compared the cylindri-cally symmetric implementation at low frequency with a DC simulation from a ResistorNetwork solution developed in MATLAB with (Figure 3 in Yang et al. (2016a)).Here, I include a comparison with the time domain electromagnetic simulation shownin Figures 13 and 14 of Commer et al. (2015). A 200m long well, with a conductivityof 106 S/m, outer diameter of 135 mm, and casing thickness of 12 mm is embeddedin a 0.0333 S/m background. For the material inside the casing, I use a conductivityequal to that of the background. The conductivity of the air is set to 3×10−4 S/m andthe permeability of the casing is ignored (µ = µ0). A 10 m long inline electric dipolesource is positioned on the surface, 50 m radially from the well. The radial electric fieldis sampled at 5 m, 10 m, 100 m, 200 m and 300 m along a line 180◦ from the source.Two simulations are included in Commer et al. (2015): a finite element (FE) anda finite difference (FD) solution. Both simulation meshes capture the thickness of thecasing with a single cell or single tetrahedral element. The finite element solution meshconsisted of over 8 million tetrahedral elements and the simulation completed in 63hours on a single core of an Intel Xeon X5550 processor (2.67 GHz). For the finitedifference solution, a conservative time-stepping was used (∆t = 3×10−10 s), resultingin a total of >120 million time steps. This simulation took 23.2 hours using 512 coreson an Intel Xeon architecture (2.33 GHz).Additionally, I include a comparison with the 3D UBC finite volume OcTree time74domain code (Haber et al., 2007). The OcTree mesh allows for adaptive refinement ofthe mesh around sources, receivers, and conductivity structures within the domain, thusreducing the number of unknowns in the domain as compared to a tensor mesh. Themesh in the UBC simulation included 5,011,924 cells, with the finest cells being equalto the width of the casing; 154 time steps were taken and 10 different step-lengths wereused (requiring 10 different matrix factorizations). This simulation took 57 minutes torun on a single Intel Xeon X5660 processor (2.80GHz).For the 3D cylindrical simulation (SimPEG), I use a mesh that has 4 cells radiallyacross the width of the casing, 2.5 m vertical discretization, and azimuthal refinementnear the source and receivers (along the θ = 90◦ line), as shown in Figure 3.3. Themesh has a total of 314,272 cells. For the time discretization, the smallest time-step Iuse is 10−6 s; the time-mesh is coarsened at later times. I used a moderately conserva-tive time-stepping scheme with 187 time-steps total. Seven different step-lengths wereemployed, requiring seven matrix factorizations. To solve the system matrix, the directsolver PARDISO was used (Petra et al., 2014; Cosmin et al., 2016). The simulation took14 minutes to run on a single Intel Xeon X5660 processor (2.80GHz).In Figure 3.4, I show the absolute value of the radial electric field sampled at fivestations; each of the different line colors is associated with a different location, and off-sets are with respect to the location of the well. Solutions were interpolated to the sameoffset using nearest neighbor interpolation.The 3D cylindrical simulation (SimPEG) isplotted with a solid line and overlaps with the UBC solution (dash-dot line) for all timesshown. The finite element (FE) solution from Commer et al. (2015) is shown with thedashed lines, and the finite difference (FD) solution is plotted with dotted lines. The3D cylindrical (SimPEG) and UBC solutions are overall in good agreement with thesolutions from Commer et al. (2015). There is a difference in amplitude and position of75Figure 3.3: Depth slice (left) and cross section (right) through the 3D cylindricalmesh used for the comparison with Commer et al. (2015). The source andrecievers are positioned along the θ = 90◦ line. The mesh extends 17kmraidally and 30km vertically to ensure that the fields have sufficiently de-cayed before reaching the boundaries.the zero-crossing (the v-shape visible in the blue and orange curves) between the Com-mer solutions and the SimPEG / UBC solutions at the shortest two offsets in the earlytimes. At such short offsets from a highly conductive target, details of the simulationand discretization, such as the construction of the physical property matrices in each ofthe various approaches become significant; this likely accounts for the discrepancies buta detailed code-comparison is beyond the scope of this chapter. My aim with this com-parison is to provide evidence that the numerical simulation is performing as expected;the overall agreement with Commer’s and UBC’s results provides confidence that it is.This example demonstrates agreement between the 3D cylindrical solution and so-lutions obtained with independently developed codes. Importantly, it also shows how,by using a cylindrical discretization which conforms to the conductivity structure of in-terest, 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. Min-76Figure 3.4: Time domain EM response comparison with (Commer et al., 2015).Each of the different line colors is associated with a different location; offsetsare with respect to the location of the well.imizing computation time was not a main focus in the development of the software andthere are still opportunities for improving efficiency. As an open-source project, contri-butions from the wider community are encouraged.3.2 Numerical ExamplesThe numerical experiments I consider in this section are motivated by the need to usethis software to delve into the physics of DC and EM in settings with steel cased wells.As such, the examples build upon examples in analytical derivations and asymptoticanalyses in Kaufman (1990); Kaufman and Wightman (1993) at DC and Augustin et al.(1989a) for frequency domain EM. I do this to provide confidence that the physicalphenomena we observe in the simulations are expected by theory as well as to build afoundation for discussion in the subsequent chapters.773.2.1 DC Resistivity Part 1: Electric fields, currents and charges ina long wellIn his two seminal papers on the topic, Kaufman uses transmission line theory to drawconclusions about the behaviour of the electric field when an electrode is positionedinside of an infinite casing. In this first example, I will revisit some of the physicalinsights discussed in (Kaufman, 1990; Kaufman and Wightman, 1993) that followedfrom an analytical derivation and compare those to my numerical results. In the secondexample, I look at the distribution of current and charges as the length of the well isvaried and compare those to the analytical results discussed in (Kaufman and Wightman,1993)I start by considering a 1 km long well (106 S/m) in a whole space (10−2 S/m), withthe conductivity of the material inside the borehole equal to that of the whole space.For modelling, I will use a cylindrically symmetric mesh. The positive electrode ispositioned on the borehole axis in the mid-point of a 1 km long well; a distant returnelectrode is positioned 1 km away at the same depth.Kaufman discusses the behavior of the electric field by dividing the response intothree zones: a near zone, an intermediate zone and a far zone (Kaufman, 1990; Kauf-man and Wightman, 1993). In the near zone, the electric field has both radial and ver-tical components, negative charges are present on the inside of the casing, and positivecharges are present on the outside of the casing. The near zone is quite localized andtypically, its vertical extent is no more than∼ 10 borehole radii away from the electrode.To examine these features in our numerical simulation, I have plotted in Figure 3.5: (a)the total charge, (b) secondary charges, (c) electric field, and (d) current density in aportion of the model near the source.78Figure 3.5: (a) Total charge density, (b) secondary charge density, (c) electric field,and (d) current density in a section of the pipe near the source at z=-500 m.79DiscussionWithin the near-zone, the total charge is dominated by the large positive charge at thecurrent electrode location and negative charges that exist along the casing wall wherecurrent is moving from a resistive region inside the borehole into a conductor. The extentof the negative charges along the inner casing wall is more evident when we look at thesecondary charge, which is obtained by subtracting the charge that would be observedin a uniform half-space from the total charge (Figure 3.5b). Inside the casing, we cansee the transition from near-zone behavior to intermediate zone behavior approximately0.5 m above and below the source; that is equal to 10 borehole radii from the sourcelocation, which agrees with Kaufman’s conclusion.In the intermediate zone, Kaufman discusses a number of interesting aspects withrespect to the behavior of the electric fields and currents which we can compare with theobserved behavior in Figure 3.5. Among them, he shows that the electric field within theborehole and casing is directed along the vertical axis; as a result no charges accumulateon the inner casing wall. Charges do, however, accumulate on the outer surface of thecasing; these generate radially-directed electric fields and currents, often referred to asleakage currents, within the formation. At each depth slice through the casing and bore-hole, 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 zonedepends on the resistivity contrast between the casing and the surrounding formationand extends beyond several hundred meters before transitioning to the far zone, wherethe 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 well-logging because they delineate theregion in which measurements can be made to acquire information about the formation80resistivity outside the well. Within the intermediate zone, fields behave like those dueto a transmission line (Kaufman, 1990), and multiple authors have adopted modellingstrategies that approximate the well and surrounding medium as a transmission line(Kong et al., 2009; Aldridge et al., 2015). I will extend this analysis in the next exampleand discuss how the length of the well impacts the behavior of the charges, fields, andfluxes.3.2.2 DC Resistivity Part 2: Finite Length WellsIn (Kaufman and Wightman, 1993), the transmission-line analysis was extended to con-sider finite-length wells. Inspired by the interest in using the casing as an “extendedelectrode” for delivering current to depth (e.g. Schenkel and Morrison (1994); Um et al.(2015); Weiss et al. (2016); Hoversten et al. (2017)), here I consider a 3D DC resistivityexperiment where one electrode is connected to the top of the well. I will examine thecurrent and charge distribution for wells ranging in length from 250 m to 4000 m andcompare those to the observations in (Kaufman and Wightman, 1993). The conductivityof the well is selected to be 106 S/m. A uniform background conductivity of 10−2 S/m isused and the return electrode is positioned 8000 m from the well; this is sufficiently farfrom the well that we do not need to examine the impact of the return electrode locationin this example. A 3D cylindrical mesh was used for the simulation.Kaufman and Wightman (1993) derives a solution for the current within a finitelength well and discusses two end-member cases: a short well and a long well. “Short”versus “long” are defined on the product of αLc, where Lc is the length of the casing andα = 1/√ST , where S is the cross-sectional conductance of the casing and has units ofS·m (S= σc2pia∆a, for a casing with radius a and thickness ∆a), and T is the transverseresistance. The transverse resistance is approximately equal to the resistivity of the sur-81rounding formation (for more discussion on where this approximation breaks down, seeSchenkel and Morrison (1994)). For short wells, αLc 1, the current decreases linearlywith distance, whereas for long wells, where αLc 1, the current decays exponentiallywith distance from the source, with the rate of decay being controlled by the param-eter α . In Figure 3.6 (a), I show current in the well for 5 different borehole lengths.The x-axis is the distance from the source normalized by the length of the well. I alsoshow the two end-member solutions (equations 45 and 53) from Kaufman and Wight-man (1993). There is significant overlap between the 250m numerical solution and theshort well approximation. As the length of the well increases, exponential decay of thecurrents becomes evident. Since α is quite small, for this example α = 2×10−3 m−1,the borehole must be very long to reach the other end member which corresponds to theexponentially decaying solution.Figure 3.6: (a) Current along a well for 5 different wellbore lengths. The x-axis isdepth normalized by the length of the well. The black dashed line shows theshort-well approximation (equation 45 in Kaufman and Wightman (1993))for a 200 m long well. The black dash-dot line shows the long-well approx-imation (equation 53 in Kaufman and Wightman (1993)) for a 4000 m well.(b) Charge per unit length along the well for 5 different wellbore lengths.In Figure 3.6 (b), I have plotted the charges along the length of the well. In the short-82well regime, the borehole is approximately an equipotential surface and the chargesare uniformly distributed; in the long well the charges decay with depth. What wassurprising to me was the noticeable increase in charge accumulation that occurs nearthe bottom of the well. This is especially evident for the short well. Initially, I wassuspicious and thought this might be due to problems with our numerical simulation;there was no obvious physical explanation that I was aware of. However, investigationinto the literature revealed that the increase in charge density at the ends of a cylinderis a real physical effect, but an exact theoretical solution does still not appear to exist(Griffiths and Li, 1997) (see figure 4, in particular).DiscussionThe results shown in Figure 3.6 have implications when testing approaches for reducingcomputational load by approximating a well with a solid tube 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, aslong as the borehole is approximated by a sufficiently conductive target, the behaviourof the fields and fluxes will be representative of the fine-scale model. However, asthe length of the well increases, the cross-sectional conductance of the well becomesrelevant as it controls the rate of decay of the currents in the well and thus the ratethat currents leak into the formation. A similar result holds when a line of chargesis used to approximate the well as a DC source; a uniform charge is suitable for asufficiently short or sufficiently conductive well, whereas a distribution of charge whichdecays exponentially with depth needs to be considered for longer wells. Thus, whenattempting to replace a fine-scale model of a well with a coarse-scale model, either witha conductivity structure or by some form of “equivalent source”, validations should be83performed on models that have the same length-scale as the experiment to ensure thatboth behaviors are being accurately modeled.3.2.3 Frequency Domain Electromagnetics Part 1: Comparisonwith scale model resultsIn the DC example, we discussed how charges are distributed along the well and cur-rents flow into the formation. From a historical perspective, practical developments inEM were pursued in the frequency domain; the mathematics is more manageable in thefrequency domain, and technological advances were being made in the development ofinduction well-logging tools (Doll, 1949; Moran and Kunz, 1962). Although the con-ductivity of the pipes generally plays the most dominant role in attenuating the signal,the magnetic permeability is non-negligible (Wait and Hill, 1977); it is the product ofthe conductivity and permeability that appears in the description of EM attenuation.Also, the fact that permeable material becomes magnetized in the presence of an exter-nal field complicates the problem. Augustin et al. (1989a) is one of the first papers oninduction logging in the presence of steel cased wells that aims to understand and isolatethe EM response of the steel cased well. Using a combination of scale modelling andanalytical mathematical modelling, they examine the impacts of conductivity and mag-netic permeability on the magnetic field observed in the pipe. In this example and theone that follows, I attempt to unravel this interplay between conductivity and magneticpermeability.The first experiment Augustin et al. (1989a) discusses is a scale model using two dif-ferent pipes, a conductive copper pipe and a conductive, permeable iron pipe; each pipeis 9 m in length. The copper pipe had an inner diameter of 0.063 m and a thickness of0.002 m, while the iron pipe had a 0.063 m inner diameter and 0.0043 m wall thickness.84A source-loop, with radius 0.6 m was co-axial with the pipe and in one experiment po-sitioned at one end of the pipe (which they refer to as the “semi-infinite pipe” scenario).In another experiment the source loop is positioned at the midpoint of the pipe (whichthey refer to as the “infinite pipe” scenario); for both experiments, magnetic field dataare measured as a function of frequency at the central axis of the pipe. Their results arepresented in terms of a Field Strength Ratio (FSR), which is the ratio of the absolutevalue of the magnetic field at the receiver with the absolute value of the magnetic fieldif no pipe is present (Figure 3 in Augustin et al. (1989a)). At low frequencies, for thedata collected within the iron pipe, static shielding (FSR < 1) was observed for the mea-surements where the receiver was in the plane of the source loop for both the “infinite”and “semi-infinite” scenarios. When the receiver was positioned within the pipe, 1.49m offset from the plane of the source loop, static enhancement effects (FSR > 1) wereobserved for both the infinite and semi-infinite scenarios. Using this experiment for con-text, I will compare the behaviour of our numerical simulation with the observations in(Augustin et al., 1989a) and examine the nature of the static shielding and enhancementeffects.For our numerical setup, the pipes are 9 m in length and have an inner diameter of0.06 m. The copper pipe has a casing-wall thickness of 0.002 m and the iron pipe has athickness of 0.004 m. Following the estimated physical property values from Augustinet al. (1989a), I use a conductivity of 3.5× 107 S/m and a relative permeability of 1for the copper pipe. For the iron pipe, a conductivity of 8.0× 106 S/m and a relativepermeability of 150 is used. A background resistivity of 104 Ωm is assumed. Thecomputed FSR values for the axial magnetic field as a function of frequency are shownin Figure 3.7.Consider the response of the conductive pipe. At low frequencies, the FSR for the85Figure 3.7: Field strength ratio (FSR), the ratio of the measured vertical magneticfield with the free space magnetic field, as a function of frequency for twodifferent receiver locations. In (a), the receiver is in the same plane as thesource, in (b), the receiver is 1.49 m offset from the source.copper pipe (blue lines) is 1 for both the infinite (solid line) and semi-infinite (dashedline) scenarios, as the field inside the copper pipe is equivalent to the free-space field.With increasing frequency, eddy currents are induced in the pipe which generate a mag-netic field that opposes the primary, causing a decrease in the observed FSR. When thesource and receiver are in the same plane (L=0.00 m), the rate of decrease is more rapidin the infinite scenario than the semi-infinite. Since there is conductive material on bothsides of the receiver in the infinite case, we would expect attenuation of the fields to86occur more rapidly than in the semi-infinite case. This observation is consistent withFigure 3a in Augustin et al. (1989a). For the offset receiver (L=1.49 m), they observeda slight separation in the infinite and semi-infinite curves which I do not; however, theyattributed this to potential errors in magnetometer position. Thus, overall, the numericalresults for the copper pipe are in good agreement with the scale model results observedby Augustin et al. (1989a).Next, we examine the response of the conductive, permeable pipe. In Figure 3.7b,we observe a static enhancement effect (FSR > 1) at low frequencies. The enhancementis larger in the infinite scenario than the semi-infinite scenario; this is in agreement withFigure 3b in Augustin et al. (1989a). There is however, a significant discrepancy be-tween our numerical simulations and the scale model for the semi-infinite pipe when thesource and receiver lie in the same plane(Figure 3.7a). Augustin et al. (1989a) observeda static shielding effect for both the infinite and semi-infinite scenarios, whereas we ob-serve a static shielding for the infinite scenario, but a significant static enhancement forthe semi-infinite case. To examine what might be the cause of this, I will examine themagnetic flux density in this region of the pipe.In Figure 3.8, I have plotted: (a) the secondary magnetic flux in the infinite-pipe sce-nario near the source (z =-4.5 m), (b) the secondary magnetic flux in the semi-infinitescenario (z=0 m for the source), and (c) top 5 cm of the semi-infinite pipe. All plotsare at 0.1 Hz. The primary magnetic field is directed upwards within the regions Ihave plotted, so upward-going magnetic flux indicates a static enhancement effect, anddownward-oriented magnetic flux indicates static shielding effects. In (a) we see a tran-sition between the static shielding in the vicinity of the source to a static enhancementapproximately 0.5 m above and below the plane of the source. Similarly in (b), we no-tice a sign-reversal in the z-component of the secondary magnetic flux at a depth of 0.687m.DiscussionThe behaviors observed in Figure 3.8 are quite comparable to Augustin et al.’s obser-vation of a transition from shielding to enhancement occurring at distances greater than0.8 m from the source. Numerical experiments show that the vertical extent of the re-gion over which static shielding is occurring increases with increasing pipe diameter,and similarly increases with increasing loop radius while the magnitude of the effectdecreases. This can be understood by considering how the pipe is magnetized; for asmall loop radius, the magnetization is largely localized near the plane of the source andrapidly falls off with distance from the plane of the source. Localized, large amplitudemagnetization causes the casing to act as a collection of dipoles around the circum-ference of the casing. As the radius of the loop increases, the magnetization spreadsout along the length of the well resulting in longer, lower-amplitude dipoles, thus bothincreasing the extent of the region over which static shielding is occurring as well asdecreasing its amplitude.This explains the nature of the static enhancement and static shielding effects, butto explain the discrepancy between the static shielding observed in the semi-infinitepipe when L=0 m by Augustin et al. (1989a), and the static enhancement we observein Figure 3.7a, I examine the magnetic flux density in the top few centimeters of thepipe. Figure 3.8c shows the top 5 cm of the secondary magnetic flux in the semi-infinitepipe; the source is in the z=0 m plane. Zooming in reveals there is yet another signreversal near the end of the pipe. This is evident even in the infinite-pipe scenario(Figure 3.7d), where the source is offset by several meters from the end of the pipe.This edge-effect perhaps bears some similarities to what we observed in Figure 3.6b,88Figure 3.8: Magnetic flux density at 0.1 Hz in the region of the pipe near the planeof the source for (a) the “infinite” pipe, where the source is located at -4.5m and the pipe extends from 0 m to -9 m, (b) a “semi-infinite” pipe, wherethe source is located at 0 m and the pipe extends to -9 m. In (c), we zoomin to the top 5 cm of the “semi-infinite” pipe, and (d) shows the 5 cm at thetop-end of the “infinite” pipe.where we saw a build up of charge near the end of the pipe in the DC scenario. At theend of the pipe, we encounter the situation where the normal component of the flux (~j,~b)from the pipe to the background needs to be continuous both in the radial and verticaldirections at the end of the pipe as does the tangential component of the fields (~e,~h). Theinterplay of these two constraints at the end of the pipe results in more complexity in theresultant fields and fluxes. Within the span of a few centimeters we transition from staticenhancement at the top of the pipe to a static shielding further down. An error as smallas a few centimeters in the position of the magnetometer causes a reversal in behavior;in Figure 3.9, I have plotted the FSR for a magnetometer positioned 3 cm beneath theplane of the source, and the static-shielding behavior observed for the semi-infinite pipeis much more aligned with that observed in Figure 3a in Augustin et al. (1989a).89Figure 3.9: Field strength ratio, FSR, for a receiver positioned 3 cm beneath theplane of the source. For comparison, we have plotted the FSR for the per-meable pipe when the source and receiver lie in the same plane (L=0.00 m)with the semi-transparent orange lines. Note that the infinite-pipe solutionsfor L=0.03 m and L=0.00 m overlap.3.2.4 Frequency Domain Electromagnetics Part 2: Conductivityand permeability in the inductive response of a wellThe experiments shown in the previous section revealed some insights into the complex-ity of the fields within the pipe and illustrated the role of permeability in the characterof the responses at low frequency. Next, we move to larger scales and examine the roleof conductivity and permeability in the responses we observe in the borehole.In this example, I consider a 2 km long well with an outer diameter of 10 cm andthickness of 1 cm in a whole-space which has a resistivity of 104 Ωm. A loop withradius 100 m is coaxial with the well and positioned at the top-end of the well. Areceiver measuring the z-component of the magnetic flux density is positioned 500 mbelow the transmitter loop, along the axis of the well. I will consider both time domainand frequency domain responses.In electromagnetics, the product of permeability and conductivity is an important90quantity in characterizing the main features of the fields and fluxes. For example, theskin depth, which describes the decay of a plane wave in a homogeneous medium, isgiven by δ =√2/(ωµσ). Most EM experiments only consider variable conductivity,however for steel-casing, both parameters are variable, which undoubtedly complicatesthe behaviour of EM fields and fluxes. To assess the contribution of each to the measuredresponses, I will investigate two scenarios. In the first, the well has a conductivity of108 S/m and a relative permeability of 1, and in the second, the well has a conductivityof 106 S/m and a relative permeability of 100; thus the product of conductivity andpermeability is equivalent for both wells.Similar to the analysis done by Augustin et al. (1989a) when looking at the role ofborehole radius in the behaviour of the magnetic response (e.g. figure 8), I will examinethe normalized secondary field (NSF) which is the ratio of the secondary field with theamplitude of the primary, where the primary is defined to be the free-space response. InFigure 3.10, I have plotted the normalized secondary field for the two pipes considered,the conductive pipe (blue) and the conductive, permeable pipe (orange). Let us start byexamining the conductivity response in Figure 3.10. Where the value of the NSF is zero,the primary dominates the response; this is the case at low frequencies where inductionis not yet contributing to the response. As frequency increases, currents are induced inthe pipe which generate a secondary magnetic field that opposes the primary, hence theNSF becomes negative. When the real part of the NSF (solid line) is -1, the secondarymagnetic field is equal in magnitude but opposite in direction to the free-space primaryand the measured real field is zero. Values less than -1 indicate a sign reversal in the realmagnetic field. Similarly, when the imaginary part of the response function goes abovezero, there is a sign reversal in the imaginary component. Note that these sign reversalsoccur even in a half-space and are a result of sampling the fields within a conductive91medium; in this case the receiver was 500 m below the surface.As compared to the conductive pipe, the frequency at which induction sets in ishigher for the conductive, permeable pipe. We also notice that the amplitude variationof both the imaginary and real parts is larger for the permeable pipe. To examine the con-tribution of conductivity and permeability to the responses, I have plotted the real part ofthe secondary magnetic flux density, b, in Figure 3.11. The top row shows the responsewithin 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 secondaryfields generated are oriented downwards. Similar to the previous example, we see thatat 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 isno measurable contribution at the receiver. At 1 Hz, we can see that induction is startingto contribute to the signal for the conductive pipe, while for the permeable pipe, it isnot until ∼10 Hz that we begin to observe the contribution of induction. At 100 Hz, thesecondary magnetic field is stronger in amplitude than the primary, and the NFS is lessthan -1 for both the conductive and permeable pipes. The amplitude of the secondarywithin the permeable pipe is stronger than that in the conductive pipe. At 1000 Hz, wehave reached the asymptote of NSF=-1 for both the conductive and permeable pipes; thesecondary 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 responsesas a function of time. For this experiment, a step-off waveform is employed and data aremeasured after shut-off, the NSF is plotted in Figure 3.12. Note here that the secondaryfield is in the same direction as the primary, so after the source has been shut off, thesecondary field is oriented upwards, as shown in Figure 3.13. Shortly after shut-off,the rate of increase in the secondary field is the same for both the conductive and the92Figure 3.10: Normalized secondary field, NSF, as a function of frequency for twowells. The NSF is the ratio of the secondary vertical magnetic field with theprimary magnetic field at the receiver location (z=-500 m); the primary isdefined as the whole-space primary.conductive, permeable wells. A maximum normalized field strength of approximately 1is reached for both cases. The responses begin to differ at 10−3 s where the conductivewell maintains a NFS∼ 1 for approximately 1 ms longer than the permeable well beforethe fields decay away.DiscussionIt is important to note that although the product of the conductivity and permeability isidentical for these wells, the geometry of the well and inducing fields results in differentcouplings for each of the parameters. For a vertical magnetic dipole source, the electricfields are purely rotational while the magnetic fields are primarily vertical. An approx-imation we can use to understand the implications of these geometric differences is toassume the inducing fields are uniform (e.g. the radius of the source loop is infinite) andto examine the conductance and permeance of the pipe. For rotational electric fields, the93Figure 3.11: Secondary magnetic flux density (with respect to a whole-space pri-mary) at five different frequencies for a conductive pipe (top row) and for aconductive, permeable pipe (bottom row).conductance isS = σtL2pir(3.9)where t is the thickness of the casing, r is the radius of the casing and L is the length-scale of the pipe segment contributing to the signal. For vertical magnetic fields, thepermeance isP = µt2pirL(3.10)As the length-scale, L, is larger than the circumference of the pipe (2pir) the geometriccontribution to the conductance is larger than that to the permeance.An important take-away from this example is that the contributions of conductivity94Figure 3.12: Normalized secondary field (NSF) through time. In the time-domain,we compute the NSF by taking the difference between the total magneticflux at the receiver and the whole-space response and then taking the ratiowith the whole-space magnetic flux prior to shutting off the transmitter.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. Thusto accurately model conductive, permeable pipes, over a range of frequencies or times,a numerical code must allow both variable conductivity and variable permeability to beconsidered.3.3 Summary and OutlookI have developed software for solving Maxwell’s equations on 2D and 3D cylindricalmeshes. The medium can have variable electrical conductivity and magnetic permeabil-ity. The 2D solution is especially computationally efficient and has a large number ofpractical applications. When cylindrical symmetry is not valid, the 3D solution can beimplemented; a judicious design of the mesh can often generate a problem with fewercells than would be required with a cartesian tensor or OcTree mesh, thus reducing thecomputational cost of a simulation. I demonstrated the versatility of the codes by mod-95Figure 3.13: Secondary magnetic flux density for a conductive well (top row) anda conductive, permeable well (bottom row) through time. The source wave-form is a step-off waveform.elling the electromagnetic fields that result when a highly conductive and permeablecasing is embedded in the earth.I presented a number of different experiments involving DC, frequency-domain, andtime-domain sources. The first two examples considered a simple DC resistivity exper-iment. In the first, I 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 currentsand charges between short and long wells. Even in this relatively simple example, thephysics was more complex than I originally anticipated; I had not intuitively expected96to see the large increase in charge density that was observed near the ends of the well.In the subsequent examples, I considered electromagnetic experiments and incor-porated magnetic permeability in the simulations. I showed that for a conductive andpermeable casing, excited by a circular current source, there is a complicated magneticfield that occurs in the top few centimeters of the pipe. Furthermore, the role of conduc-tivity and permeability in the observed responses is more complex than their product;the source geometry and coupling with the casing are important to consider.As new strategies and software are developed to handle more complex well-geometries,such as deviated or horizontal wells, it is important that we establish an understandingof the physics. Of critical importance is the ability to plot the charges, fields, and fluxesin the simulations. This is valuable for understanding the responses obtained from theexperiment and it is a solid foundation for designing a field survey. I anticipate thatthe software provided with this thesis can be a resource for building understanding andadditionally, serve as a tool for testing 3D simulations with boreholes present.The software implementation is included as a part of the SimPEG ecosystem. Sim-PEG also includes finite volume simulations on 3D tensor and OcTree meshes as wellas machinery for solving inverse problems. This means that the cylindrical codes can bereadily connected to an inversion and additionally, simulations and inversions of morecomplex 3D geologic settings can be achieved by coupling the cylindrical simulationwith a 3D tensor or OcTree mesh using a primary-secondary approach (e.g. example 3 inAppendix D). Beyond modelling steel cased wells, the 3D cylindrical mesh could proveto be useful in conducting 3D airborne EM inversions where a domain-decompositionapproach, similar to that described in Yang et al. (2014), is adopted.97Chapter 4Direct current resistivity withsteel-cased wells4.1 IntroductionSubsurface resistivity can be a valuable part of a geologic interpretation, whether that beidentifying lithologic units, characterizing changes within a reservoir, or imaging sub-surface injections associated with carbon capture and storage or hydraulic fracturing. Inmany of these settings, steel-cased wellbores are present. Steel has a significant elec-trical conductivity, which is generally six or more orders of magnitude larger that ofthe surrounding of the geologic formation. Clearly, such a large contrast is importantto consider when conducting a direct current (DC) resistivity survey. On one-hand, therole of the steel casing may be viewed as “distortion” which complicates the signalsof interest (Wait, 1983; Holladay and West, 1984; Johnston et al., 1987). In other sce-narios, a wellbore may be beneficial in that it can serve as an “extended electrode” sothat current-injection and sampling of the resultant electrical potentials can take place98beneath near-surface heterogeneities (Ramirez et al., 1996; Rucker et al., 2010; Rucker,2012; Ronczka et al., 2015) or so that currents injected at the surface can reach signif-icant depths (Schenkel and Morrison, 1994; Weiss et al., 2016; Hoversten et al., 2017).The use of casings as extended electrodes extends back several decades. Sill and Ward(1978) used the well casing as a buried electrode for their mis-a`-la-masse experimentat the Roosevelt Hot Springs geothermal field in Utah, as did Kauahikaua et al. (1980)for their mis-a`-la-masse mapping of a high temperature geothermal reservoir in Hawaii.Sill (1983) used the well as a source to monitor an injection test at Raft River, Idaho todetermine if measurable changes that might indicate the direction of fluid flow could beobserved. Rocroi and Koulikov (1985) delineated a known resistive hydrocarbon depositin the USSR by injecting current into two cased wells. More recently, applications forhydraulic fracturing, enhanced oil recovery and carbon capture and storage have beenof much interest (Commer et al., 2015; Um et al., 2015; Weiss et al., 2016; Hoverstenet al., 2017).To build a physical understanding of electrical and electromagnetic methods in set-tings where steel-cased wells are present, there are several areas to be investigated. First,the significant conductivity of the steel will impact the behavior of the charges, currents,and electric fields. This is true at the electrostatic limit, relevant to DC resistivity sur-veys, as well as when the source fields are time-varying, as in electromagnetic (EM)surveys. When considering EM surveys, induction effects also influence the responses,and magnetic fields and fluxes become relevant, meaning that the magnetic permeabilityof the steel then introduces further complexity into the signals we measure. This chapteris concerned with the first set of physical phenomena: understanding the physics of steelcasings at DC.Much of the initial theory and understanding of the behavior of electric fields, cur-99rents, and charges, was developed in the context of well-logging. Kaufman (1990) andKaufman and Wightman (1993) provide a theoretical basis for our understanding; thefirst paper derives an analytical solution for a DC experiment where an electrode is po-sitioned along the axis of an infinite length well, and discusses where charges accumu-late and how currents leak into the surrounding formation. From this, Kaufman (1990)shows that by measuring the second derivative of the electric potential, informationabout the formation resistivity can be obtained. The second paper extends the analysisfor finite length wells. Schenkel and Morrison (1990); Schenkel (1991); Schenkel andMorrison (1994) pioneered numerical work analyzing the influence of steel-cased wellson geophysical data using an integral equation approach for solving the DC resistivityproblem. They expand upon the logging-through-casing application and discuss limita-tions of the transmission line solution presented in Kaufman (1990) for this application.They also explored the feasibility of cross-hole and borehole-to-surface surveys whereone electrode is placed within, or beneath, a cased borehole. These examples demon-strated that the casing can improve detectability of a conductive target as compared tothe scenario where no cased well is present.In this chapter, I focus on three aspects of DC resistivity in the presence of steel-cased well. In section 4.2, I examine the feasibility of conducting a surface DC surveyto detect a flaw in the casing and discuss factors that influence detectability of a flaw. Insection 4.3, I examine the use of DC resistivity for geophysical imaging when a steel-cased well is present. Finally, in section 4.4, I assess strategies applied in the literaturefor approximating a steel-cased well with a coarse-scale model to reduce computationalcost.Source codes for all of the simulations shown are open source, licensed under theMIT license, and are available as Jupyter notebooks (see Appendix A). The examples100in the paper have been selected with an emphasis on examining physical principles;however, I envision that the Jupyter notebooks included with this chapter could serve asuseful survey design tools.4.2 DC resistivity for casing integrityDegraded or impaired wells can pose environmental and public-health hazards. A flawin the cement or casing can provide a conduit for methane to migrate from depth intogroundwater aquifers or into the atmosphere. This is particularly of concern for shalegas wells. Elevated levels of thermogenic methane, which are attributed to deep sources(rather than biogenic methane which can be generated closer to the surface), in ground-water wells in Pennsylvania has been positively correlated with proximity to shale gaswells in the Marcellus and Utica (Osborn et al., 2011; Jackson et al., 2013), and failurerates of unconventional wells (e.g. shale gas wells) is estimated to be 1.57 times largerthan that of a conventional well drilled in the same time-period (Ingraffea et al., 2014).Wells can fail if there is a compromise in the cement or the casing. To diagnose the in-tegrity of a well with electrical methods, we require a contrast in electrical conductivityto be associated with the flaw, thus I will focus attention to detecting flaws in the highlyconductive casing.Under what circumstances should we be able to detect a flaw in the casing using DCresistivity from the surface? To address this question, I begin by examining how a flawwhich comprises the entire circumference of the pipe along some depth interval changesthe charge distribution and thus the resultant electric fields we measure on the surface.From there, I investigate the role of parameters including the depth of the flaw and thebackground conductivity on our ability to detect it from the surface. Finally, I examinethe scenario in which only a portion of the circumference of the pipe is flawed.1014.2.1 Basic experimentThe experiment I consider is a “top-casing” DC resistivity experiment where one elec-trode is connected to the wellbore at the surface and a return electrode is positionedsome distance away. The concept and basic physics is the same as a mis-a`-la-masse sur-vey in which the positive electrode is connected to a conductive target. When the sourceis turned on, positive charges are distributed on the interface between the conductivetarget and the resistive host. Electric potentials are measured on the surface and thesedata are then used to infer information about the extent of the conductor (Telford et al.,1990). Applying the same principles to a casing integrity experiment, I connect a posi-tive electrode to the casing, and for an intact casing, positive charges will be distributedon the outer interface of the casing along its entire length. If corrosion causes a flawacross the diameter of the casing, the continuity of the conductive flow-path for chargesis interrupted, thus we expect a larger charge on the top portion of the flawed casingthan we would if it were intact. This results in a larger electric field at the surface thanwould be observed if the casing were intact. The difference in electric field (or electricpotentials) from the expected electric field that results from an intact well could then bean indicator that there is a problem with the well.To demonstrate the principles, I start by considering a simple model of a casing ina half-space. The intact well is 1 km long, has an outer diameter of 10 cm, a thicknessof 1 cm and a conductivity of 5× 106 S/m. The background is 10−1 S/m, and theconductivity of the inside of the well is taken to be equal to that of the background.The positive electrode is connected to the top of the casing and the return electrode ispositioned 2 km away. To simulate the physics, the 3D cylindrical DC code describedin Chapter 3 was employed. In Figure 4.1 I show cross-sections of the: (a) electricalconductivity model, (b) current density, (c) charge density, and (d) electric field for102Figure 4.1: Cross section showing: (a) electrical conductivity, (b) current density,(c) charge density, and (d) electric field for a top-casing DC resistivity ex-periment over (top) an intact 1000 m long well and (bottom) a 1000 m longwell with a 10 m flaw at 500 m depth.the intact well (top row) and a flawed well (bottom row) that contains a 10 m gap in thecasing at 500 m depth. As expected, the introduction of a resistive flaw prevents currentsfrom reaching the bottom portion of the well. This results in increased currents, chargedensity and thus electric fields within the top 500 m.To quantify the charge along the length of the well, I have plotted the charge as afunction of depth for the intact well (black), flawed well (blue), and also a short well of500 m length (grey dash-dot) in Figure 4.2a. In each of the wells, we observe that there103is an increase in charge density near the end of the discontinuity along the length of thewell. This was also noted in Griffiths and Li (1997) and in Chapter 3 – this behavioris attributed to edge effects. At an interface between materials with two different con-ductivities, the normal component of the current density must be conserved, as well asthe tangential component of the electric field; the discontinuity at the end of the pipe,and at the location of the flaw, means the continuity conditions must be preserved si-multaneously in the radial and vertical directions, and this complicates the behavior ofthe fields, fluxes and charges. Another observation is that the flawed and short wellshave nearly identical charge distributions in the top 500 m. In the bottom portion ofthe flawed well, where the remaining conductive material is, a small dipolar charge isintroduced, but this is nearly an order of magnitude smaller than the charge in the topportion of the pipe. The signal due to the flaw can be defined as the difference betweenthe total response due to a flawed well and the total response due to an intact well (theprimary); I will refer to this difference as the secondary response. The secondary chargeis dipolar in nature with positive charge above the flaw and negative charge beneath theflaw. Note that the charge distributions along the short well, truncated where the flawstarts at 500 m depth, and along the top portion of the flawed well are almost identical;these charges are the source of signal for a surface electric field measurement. Thissuggests that an inversion strategy, where one attempts to estimate the length of a well,may be an effective approach for characterizing the depth to a flaw.Impact of the vertical extent of the flawA 10 m flaw is quite long and it is of interest to see how the results are changed if theflaw has a smaller vertical extent. The distribution of charges shown in Figure 4.2 hintsthat the flaw may not need to be very long in order to still significantly influence the104Figure 4.2: (a) Charge along the length of the intact well (black), a 500 m well (“short”, grey dash-dot), and a well with a 10 m flaw at 500 m depth (blue), ina top-casing DC resistivity experiment. (b) Secondary electric field due onthe surface of the earth due to the flaw in the casing. The primary is definedas the electric field due to the 1000 m long intact well. The return electrodeis 2000 m away from the well.response. To confirm this, I adopt a much finer vertical discretization in order to modelsmaller flaws. Here, I use a shorter, 50 m long well in order to reduce computationalload. The flaw is positioned at 25 m depth, and the length of the impairment is varied.This simulation is conducted on a cylindrically symmetric mesh, the positive electrodeis connected to the casing, and a return electrode is positioned 50 m away.The resultant charge distributions are shown in Figure 4.3. For comparison, I againshow the charge on a well that is truncated at the location of the flaw; this is the “short”well and results are displayed using the grey dash-dot line. The charge distribution issimilar for all of the flawed-well scenarios, even for flaws smaller than the thickness ofthe casing (10−2 m). We see similar behavior to that shown in Figure 4.2, where posi-tive charge accumulates within the top portion of the well and a small dipole charge ispresent in the bottom portion of the well. There are minor differences in amplitude asthe vertical extent of the flaw is changed; as the extent of the flaw decreases, the ampli-105Figure 4.3: Charge along the length of a 50 m long intact well (black), a 25 mwell (“short”, grey dash-dot), and four wells, each with a flaw starting at25 m depth and extending the length indicated by the legend (5× 10−1 m(blue), 5× 10−2 m (orange), and 5× 10−3 m (green)) in a top-casing DCresistivity experiment. For reference, the diameter of the casing is 10−1 mand its thickness is 10−2 m. The return electrode is 50 m away from the welland a cylindrically symmetric mesh was used in the simulation.tude of the dipolar charge on the bottom portion of the well increases slightly while theamplitude of the positive charge on the top portion of the well decreases. These distinc-tions, however, are small in magnitude, and even if the background is more conductive,the casing is still orders-of-magnitude larger in conductivity than any geologic materialwe are likely to encounter. Thus, I conclude that, so long as the impairment affects theentire circumference of the casing, the extent of that flaw has little impact on the chargethat accumulates in the top portion of the well. As such, I will proceed in our analysisusing a 10 m flaw in the 1 km well so that a fine vertical discretization is not necessary.1064.2.2 Survey design considerationsWhen examining detectability of a signal, there are two aspects to consider: (1) thesignal must be larger than the noise floor of the instrument, and (2) the signal must bea significant percentage of the primary; for the casing integrity experiment, the primaryis the signal due to the intact well. Due to the cylindrical symmetry of the charge onthe well, we expect the electric field at the surface to be purely radial, thus only radialelectric field data need be collected at the surface.In Figure 4.4, I have plotted the primary field (top row), secondary field (secondrow) and secondary field as a percentage of the primary (third row) for four differentreturn electrode locations. In (a), the return electrode is 2000 m offset from the well,in (b) the offset is 750 m, in (c) the offset is 500 m, and in (d) the offset is 250 m. Inaddition to the plan view, I have plotted the primary electric field (black), total electricfield for the flawed well (blue) and secondary radial electric field (orange) along theθ = 90◦ azimuth in the fourth row of Figure 4.4. The fifth row shows the secondary asa percentage of the primary.At the furthest offset (Figure 4.4a), there is nearly complete cylindrical symmetry inthe primary field. With complete cylindrical symmetry there is no preferential directionalong which to collect data. As I move the return electrode closer, for example to 750 mfrom the well, we notice that the secondary electric field does not change substantially.However, if I examine the ratio of the secondary to the primary (second and fifth rows),we see that the ratio has increased. Although the primary field has similar, if not largeramplitude near the well, it also has considerable curvature. As a result, the proportionof the primary field that is in the radial direction has decreased in amplitude. Hence theimportant characteristic, the ratio of the secondary to primary of the radial components,has increased. The above principles are further enhanced as the return current is brought107Figure 4.4: (Top row) primary electric field, (second row) secondary electric field,and (third row) secondary electric field as a percentage of the primary radialelectric field for a return electrode that is offset (a) 2000 m, (b) 750 m, (c)500 m, and (d) 250 m from the well. The primary is defined as the responsedue to the 1000 m long, intact well. Electrode locations are denoted by thered dots. In the third row, the colorbar has been limited between 20% and100%. The fourth and fifth rows show radial electric field data collectedalong the θ = 90◦ azimuth (white dotted lines in the top three rows). Thefourth row shows the primary (black line), the total electric field due to theflawed well (blue line), and the secondary radial electric field (orange line).The fifth row shows the secondary as a percentage of the primary.108closer to the well as in panels (c) and (d), where the return electrode is brought to 500m and 250 m from the well. Again, for all of these examples the amplitude of thesecondary field at the surface is quite similar. However, the choice of azimuth for thesurvey line will greatly affect the size of the ratio. In terms of survey design, we cantake advantage of the return electrode to reduce coupling with the primary.For our following examples I will place the return electrode at 500 m from the welland collect radial data along a line that is perpendicular to the source line. I will examineseveral factors influencing detectability of a flaw, including the depth of the flaw andthe conductivity of the background in the following sections. I will also examine thescenario where only a portion of the circumference of the well has been compromised.4.2.3 Factors influencing detectabilityDepth of the flawThe introduction of a flaw in the well changes the distribution of charges along the lengthof the well and causes a secondary dipolar charge centered about the flaw. The positionand strength of this dipole will affect our ability to detect the flaw. To examine this,I have taken the same model of a 1km pipe in a 10−1 S/m background and varied thedepth of the flaw from 300 m to 900 m. In Figure 4.5, I plot radial electric field resultsalong a line perpendicular to the source electrodes; the return electrode is positioned500 m from the well. In (a), I show total radial electric field, in (b) the secondary radialelectric field (with the primary being the electric field resulting from the intact well,shown in black in panel a), and in (c) I show the secondary radial electric field as apercentage of the primary. I have also indicated where values fall below a 10−7 V/mnoise floor on Figure 4.5 (a) and (b), as well as those that fall below a 20% threshold in(c). A threshold of 20% may be conservative, however, it does depend on knowledge109of the background conductivity as well as the geometry and physical properties of thewell. In many scenarios, these may not be well-constrained, thus I select a conservativethreshold for this analysis. Any detectability analysis will be site-dependent and I havetherefore made all source-code available so that a similar workflow may be followedand adapted to include setting-specific parameters.When a well is impaired, the total radial electric field is larger than that due to thebaseline, intact well. The strength of the secondary response decreases as the depth ofthe flaw increases. For this example of a 1000 m long well in a 10−2 S/m background, aflaw at 900 m depth is not detectable; there is no overlap between the region in which thesecondary electric field (Figure 4.5b) is above the noise floor and the region in whichthe secondary comprises a significant percentage of the primary (Figure 4.5c). Thismight be expected, as the difference between the charges distributed along a 900 m longsegment versus the 1000 m long well are not drastically different. For a flaw at 700 mdepth, there is a window between 400 m offset and 800 m offset over which the radialelectric field data are sensitive to the flaw. As the depth to the impairment decreases,both the spatial extent over which data are sensitive to the flaw, and the magnitude ofthe secondary response in those data, increase.Background conductivityThe total charge on the well is controlled by the contrast in conductivity between thesteel-cased well and the surrounding geology. Increasing the conductivity of the back-ground reduces that contrast thus reducing the amount of charge on the well. The resultis a decrease in the total electric field at the surface. Similarly, the strength of the sec-ondary dipolar charge introduced with the presence of an impairment also depends uponthe available charge and will also be reduced with increasing background conductivity.110Figure 4.5: Radial electric field as the depth of the flaw along a 1km long wellis varied. The positive electrode is connected to the top of the casing, thenegative electrode is positioned 500m away and data are measured along aline 90◦ from the source electrodes. In (a), we show the total electric fieldfor four flawed wells, each with a 10m flaw at the depth indicated on thelegend. The black line shows the radial electric field due to an intact well;we define this as the primary. In (b), the secondary radial electric field isplotted and in (c), we show the secondary radial electric field as a percentageof the primary.111In Figure 4.6, I have adopted the same model of a 1 km well with a 10 m impairmentat 500 m depth, and show the radial electric field for the flawed (solid lines) and intact(dashed lines) well as the background conductivity is varied. A resistive backgroundpromotes the strongest total and secondary signals. As the conductivity increases, de-tectability becomes more challenging; at a conductivity of 3×10−1 S/m, the flaw at 500m depth is undetectable as there is no overlap in the regions where the secondary signalis above the noise floor and where it comprises a significant percentage of the primary.Variations in the background geology will also influence the distribution of chargesand thus the measured signal at the surface. To examine the challenges introduced whenvariable geology is considered, I will introduce a layer into the model and vary its con-ductivity. The layer is 50 m thick and its top is at 400 m depth. The flaw will again bepositioned at 500 m depth, and the background conductivity is 10−1 S/m. The returnelectrode is 500 m from the well, and radial electric field data are measured along aline perpendicular to the source. In Figure 4.7, I show data for a flawed well (solid)and intact well (dashed) for scenarios in which a conductive or resistive layer is posi-tioned above the flaw. The presence of a resistive layer improves detectability, while aconductive layer reduces detectability.To understand the physical phenomena governing this, I have plotted a cross sectionthrough: (a) the model, (b) the currents, (c) the charges, and (d) the electric field inFigure 4.8 for (first row) a model of an intact well with a conductive layer present,(second row) a flawed-well model including a conductive layer, (third row) a model ofan intact well with a resistive layer, and (fourth row) a flawed-well model including theresistive layer. For the comparison, there is two orders of magnitude difference betweenthe background and the layer. When a conductive layer is present, we see that it actsto “short-circuit” the system as there is significant current leak-off into that layer. This112Figure 4.6: Radial electric field as the conductivity of the background is varied fora 1 km well with a 10 m flaw at 500 m depth. The positive electrode is con-nected to the top of the casing, the negative electrode is positioned 500 maway and data are measured along a line 90◦ from the source electrodes. In(a), we show the total electric field for five different background conductiv-ities, each indicated on the legend. The solid lines indicate the response ofthe flawed well and the dashed lines indicate the response of the intact well(the primary). In (b), the secondary radial electric field is plotted and in (c),we show the secondary radial electric field as a percentage of the primary.113Figure 4.7: Radial electric field as the conductivity of a 50 m thick layer positionedat 400 m depth is varied. The positive electrode is connected to the top ofthe casing, the negative electrode is positioned 500 m away and data aremeasured along a line 90◦ from the source electrodes. In (a), we show thetotal electric field five different layer conductivities. The black line shows thescenario where the layer has the same conductivity as the background. Thedashed-lines indicate the intact well and the solid lines indicate the flawedwell. In (b), the secondary radial electric field is plotted (with respect to anintact well primary) and in (c), we show the secondary radial electric field asa percentage of the primary.114reduces the amount of current that reaches the flawed section of the well and decreasesthe total charge on the well, which is the source of our signal. Conversely, when aresistive layer is present, there is less leak-off of currents. In fact, Yang et al. (2016a)showed that rather than leaking-off, currents can enter the casing if a resistive layer ispresent. In terms of detecting a flaw beneath a resistive layer, this means that the currentdensity and charge along the well increases, thus amplifying the response due to theflaw.Conductivity of the casingThe conductivity of the casing is also relevant to how the charges are distributed alongits length. For highly conductive wells, the charge along the length of the well is ap-proximately uniform, for more resistive wells, the charges follow an exponential decay,as shown in Figure 4.9. Schenkel (1991) described the decay of currents, and thus thedistribution of charges along the length of a well, in terms of the conduction length,δL =√Scσ0=√2pirtσcσ0(4.1)Where Sc is the cross-sectional conductance of the casing (Sc= 2pirtσc for a casing withradius r, thickness t, conductivity σc and has units of [S · m]) and σ0 is the conductivityof the background. The conduction length is akin to skin depth in electromagnetics andis the depth at which the amplitude of currents have decreased by a factor of e−1. Casingconductivities of 5×105 S/m, 5×106 S/m, and 5×107 S/m correspond to conductionlengths of ∼ 180 m, 560 m, 1800 m. For the most resistive well shown, 5× 105 S/m,the vast majority of current has decayed well before it reaches the flaw; the majority ofcharges are concentrated where the currents leak off, near the top of the well. Corre-spondingly, there is greater sensitivity to a flaw in a conductive well than in a resistive115Figure 4.8: Cross section showing: (a) electrical conductivity, (b) current density,(c) charge density, and (d) electric field for a top-casing DC resistivity exper-iment over models with a conductive layer (top two rows) and a model witha resistive layer (bottom two rows). The layer extends from 400 m to 450 mdepth. The plots in the second and fourth rows correspond to models with a10 m flaw at 500 m depth.116Figure 4.9: (a) Charge along the length of wells with three different conductivities(each indicated by a different color in the legend). The intact wells are de-noted with dashed lines and the flawed wells are denoted with solid lines. (b)Secondary charge along the flawed and short wells. The primary is definedas the electric field due to the 1000 m long intact well. The return electrodeis 2000 m away from the well.well, as is reflected in the radial electric field data shown in Figure 4.10.Partial flawThe above examples considered an impairment that affects the entire circumference ofthe casing. This may be suitable in some scenarios where a particular geologic unitsubjects the well to corrosive conditions, however, flaws may also be vertical cracksalong the well (e.g. if pipe burst occurs). This is a much more challenging problem forDC resistivity because, if only a portion of the circumference is impaired, there is stilla high-conductivity pathway for currents to flow along the entire length of the well. Toexamine the feasibility of detecting a partial flaw, I have run simulations where half ofthe circumference of the casing is compromised, leaving the other-half intact.I consider four different depth extents of the flaw between 10 m and 300 m; in allscenarios the top of the flaw is at 500 m. In Figure 4.11a, I have plotted the total radialelectric field resulting from an intact well (black), wells where the entire circumference117Figure 4.10: Radial electric field as the conductivity of the casing is varied for a 1km well with a 10 m flaw at 500 m depth. The positive electrode is con-nected to the top of the casing, the negative electrode is positioned 500 maway and data are measured along a line 90◦ from the source electrodes. In(a), we show the total electric field for three different casing conductivities,each indicated on the legend. The solid lines indicate the response of theflawed well and the dashed lines indicate the response of the intact well (theprimary). In (b), the secondary radial electric field is plotted and in (c), weshow the secondary radial electric field as a percentage of the primary.118is compromised (solid) and wells in which 50% of the circumference has been compro-mised (dashed); (b) and (c) show the secondary radial electric field and the secondaryas a percentage of the primary, respectively. We see that the depth-extent of the flaw haslittle impact on the fully-compromised wells, which is consistent with the observationsin our previous examples. However, if the well is partially flawed, we do see variationin the secondary response. By compromising 50% of the circumference of the well,we have reduced the effective cross-sectional conductance over that portion of the well.Numerical experiments show that if, instead of introducing a flaw which comprises 50%of the circumference of the well, we reduce the conductivity of the intact well by 50%over the same depth extent as the flaw, we obtain similar, but not identical, responsesat the surface. Although for extensive flaws, there is a small region over which the sec-ondary signal is above the noise floor, there are no regions where this coincides withmeasurements where the secondary fields are a significant percentage of the primary.There may be a subset of circumstances, such as if the flaw is near to the surface, or ifthe background geology is sufficiently well-known so that the percent threshold can bereduced, where a partial flaw may be diagnosed, however, these results demonstrate thata partial flaw is a challenging target for a DC resistivity survey.4.2.4 SummaryIn summary, I provided an overview of the fundamental physics governing the behaviorof currents, charges, and electric fields in a top-casing DC resistivity experiment todetect an impairment in the well. If a flaw comprises the entire circumference of somedepth interval along the casing, then the charges are concentrated in the portion of thewell above the flaw, and to first approximation, the charge distribution is equal to that ofa well which has been truncated at the depth of the flaw. This excess charge is the source119Figure 4.11: Radial electric field as the vertical extent of the flaw is varied. Thepositive electrode is connected to the top of the casing, the negative elec-trode is positioned 500 m away and data are measured along a line 90◦ fromthe source electrodes. In (a), we show the total electric field four differentflaw extents. The black line shows the response of the intact well. Thedashed-lines indicate the partially flawed wells (50% of the circumferenceis compromised) and the solid lines flawed wells in which the entire cir-cumference of the well has been compromised. In (b), the secondary radialelectric field is plotted (with respect to an intact well primary) and in (c),we show the secondary radial electric field as a percentage of the primary.120of our signal. As it is cylindrically symmetric, the resultant secondary electric fields dueto the flaw are purely radial. In terms of survey-design, we can take advantage of thisknowledge and use the return electrode location to reduce coupling with the primaryelectric field in our data (as shown in Figure 4.4). Our ability to detect a flaw acrossthe entire circumference of the casing depends upon the conductivity of the backgroundand casing, as well as the depth of the flaw. Larger contrasts between the casing andthe background (e.g. a more resistive background and / or a more conductive casing)increase the secondary response, as does decreasing the depth of the flaw. If only aportion of the circumference is impaired, leaving a conductive pathway connecting thetop and bottom portions of the casing, the secondary signal is small and thus will bechallenging to detect under most circumstances.For the subset of scenarios where we do have data sensitivity to the flaw, an inverseproblem can be solved to estimate the depth of the impairment. One approach would beto use a reduced modeling procedure whereby only a few parameters are sought. For thecase presented here, we might invert for a smooth background, the length of the well,and potentially the conductivity of the casing, if it is not known a-priori.In the next section, I transition from viewing the casing as the target to working onthe scale of a geophysical imaging application in reservoir monitoring and viewing thecasing as a high-conductivity feature present in that setting.4.3 Survey design for exciting targets at depthThere are many problems in hydraulic fracturing, carbon capture and storage and en-hanced oil recovery that require targets to be illuminated and data to be acquired andinverted. Typically these experiments include steel-cased wells and the target of interestcould be resistive or conductive. The target could be immediately adjacent to a well or121offset from it, and the survey may employ electrodes on the surface or positioned down-hole. Similarly, receivers may be positioned on the surface or in adjacent boreholes.Each of these factors influences our ability to detect a target in our data.Detectability of a target requires two steps: (1) source fields must excite the target,and (2) receivers must be positioned so that the secondary response is measurable. Inthis section, I focus on the first point – exciting the target. I will examine the impact ofsource electrode locations, the physical properties of the target and the geometry of thetarget on our ability to excite a response.4.3.1 Source locationI begin by examining the impact of the source electrode location on our ability to delivercurrent to a region of interest in the model. The model I consider is 1 km long well in a10−1 S/m background. The well has a conductivity of 5×106 S/m, an outer diameter of10 cm thickness, and a 1 cm thickness; these are the same parameters used for the casingintegrity experiment described in the previous section. The conductivity of the fluidfilling the casing is identical to that of the background. I am interested in effects nearthe well and thus the modeling can be carried out using the 2D cylindrical mesh providedthat the return electrode is sufficiently far away. The return electrode is physically a discof current at a radius equal to the distance of the return electrode from the well, in thiscase 2 km. The assumption of cylindrical symmetry and the use of a distant returnelectrode has similarly been applied in Schenkel (1991).To examine the impact that the source electrode location has on our ability to excitea target, I consider the electrode locations shown in Figure 4.12. Three of the electrodesare connected to the casing (tophole - blue, centered - green, and downhole - red); theremaining electrodes are not connected to the casing; these include the surface electrode122Figure 4.12: Electrode locations to be compared. The top casing electrode (blue),centered electrode (green, 500 m depth), and downhole electrode (red, 500m depth) are connected to the casing. The surface electrode (orange) is off-set from the well by 0.1 m. The remaining electrodes are positioned alongthe axis of the casing. Panel (a) shows the entire length of the casing, while(b) zooms in to the bottom of the casing to show the separation between theelectrodes beneath the casing.(orange) as well as the five electrodes near the end of the pipe (purple - within the pipe,brown, pink, grey and yellow are beneath the end of the pipe). The surface electrode isoffset from the well by 0.1 m.To assess the ability of each electrode configuration to excite a geologic target ofinterest, I will examine the current density in the formation. In Figure 4.13, I haveplotted the amplitude of the current density along a vertical line (a) 25 m, (b) 50 m, and(c) 100 m radially offset from the well. In terms of survey design, we wish to choosea source location that maximizes the total current density within the depth region ofinterest. If the target is near the surface, we choose an electrode which is connectedto the top of the casing, or near the casing at the surface. Interestingly, at depth, thereis little distinction between these two scenarios. Thus, if one is limited to deploying123electrodes at the surface, and for practical purposes, connecting infrastructure to thewell-head presents a challenge, then grounding the electrode near the well still results ina survey that benefits from the well acting as a high-conductivity pathway to help delivercurrent to depth. If the aim however, is to excite a deeper target, we see that positioningthe electrode downhole can significantly increase the current density delivered to thatdepth. For example, if we have a target near 500 m depth, positioning the electrodenear that depth nearly doubles the current density as compared to an electrode at thesurface. If a target is near the end of the well, between 800 m and 1000 m depth, thenpositioning an electrode near the end of the well triples the current density. This effectwill be amplified if the well is lengthened, since we observe exponential decay of thecurrents carried along according to the conduction length (equation 4.1).Kaufman (1990) pointed out that the difference between an electrode positionedalong the axis of the casing and one coupled to the casing at depth is highly localizedaround the source, and thus is not an important distinction at the scales we consider fora geophysical imaging survey. I can test this numerically by comparing the currentsarising from the electrode which is connected to the casing 5 m above the bottom of thecasing (red in Figures 4.12 and 4.13), and the electrode positioned along the axis of thecasing 1.25 m above the bottom of the casing (purple in Figures 4.12 and 4.13). Indeed,we see that the red and purple lines overlap for all offsets in Figure 4.13, indicating thatboth situations result in the same distribution of currents within the formation.For electrodes beneath the casing, the distribution of currents is significantly differ-ent. For electrodes 1.25 m, 5 m, 10 m and 20 m below the pipe, we see that within∼ 100m above and below the electrode location, the currents are nearly symmetric, followingthe expected response of a point source. I have included a simulation with the electrode20 m below the pipe when there is no casing present; this is shown in black in Figure124Figure 4.13: Total current density along a vertical line offset (a) 25 m, (b) 50 mand (c) 100 m from the axis of the casing, which extends from the surface(0 m) to 1000 m depth. The electrode locations correspond to those shownin Figure 4.12. For reference, a simulation with an electrode 20 m belowthe casing when there is no casing present is shown in black.4.13. The main difference between the distribution of currents for each of these scenar-ios is the reduction in current density in the top 1000 m, with increasing electrode depth;as the electrode is moved deeper, less current is channeled into the casing. Schenkel andMorrison (1990) noted that for electrodes positioned beneath a well, if the electrode ismore than 100 casing diameters beneath the casing, then the casing has little impact onthe fields below or far from the pipe. The current is much more localized if the electrode125is beneath the casing, and thus if a target is beneath or very near the end of the well, thenit is advantageous to position the electrode beneath the well.Not surprisingly, if the source electrode can be positioned near the depth region ofinterest, the current density delivered to that region is larger. Numerical experimentsshow that the position of the return electrode makes minimal impact on the currents atdepth. However, if the return electrode is within tens of meters of the well, the nearsurface currents are significantly altered. This is consistent with our observations insection 4.2, where I showed that the return electrode location has little impact on themagnitude of the secondary signals, but its position alters the geometry of the sourcefields and this can be used to reduce coupling of receivers to the primary field.4.3.2 Target propertiesThe physical property contrast between the target and the background, the target’s ge-ometry and proximity to the well, all influence our ability to observe its impact on thedata we measure. The purpose of this section is to explore the impact of these factorson the excitation and detection of the target. In the first example, I examine the role ofthe conductivity of a cylindrical target which is in contact with the well. The second ex-ample is again a cylindrically symmetric co-axial disc target but there is a gap betweenthe casing and the target. The final example is fully 3D; the target is a block and I lookat the excitation as a function of the distance of a block from the well.Target in contact with the wellFirst, I consider a cylindrical target that is in contact with the well. Schenkel and Mor-rison (1994) examined such a scenario for a conductive target (e.g. a steam injection orwater flood) in a mis-a`-la-masse type experiment where a source electrode is connectedto the casing at the same depth as the center of the target. They considered a cross-well126experiment with potential electrodes in an offset, uncased well, and compared two sce-narios for the source well: one in which the source well is an open-hole, and the secondin which it was cased. They demonstrated that the casing enhances the response, andthus the data sensitivity to the target, as compared to an experiment where current isinjected directly into the target and no casing is present. In this example, I build uponthose findings and examine the role of the conductivity of the target on our ability toexcite it as well as the impact on the data if the target is not directly in contact with thewell.The model I use is a 1 km casing in a half-space with a target. The target extends25 m vertically and has a 25 m radius and the depth to its top is 900 m. The modelis cylindrically symmetric and thus we expect that the secondary electric field at thesurface due to the target will be purely radial. As such, I apply the lessons learned fromthe casing integrity example and use the return electrode to reduce coupling with theprimary field along a line perpendicular to the source. I position the return electrode500 m from the well-head and compare both top-casing and down-hole source electrodelocations.We begin by examining the physical behavior governing the DC response of a con-ductive and resistive target. Figure 4.14 shows the (a) conductivity model, and resultant:(b) current density, (c) charge density, and (d) electric fields for a conductive target (10S/m, top row) and a resistive target (10−3 S/m, bottom row) in a down-hole experimentwhere the source electrode is positioned at the center of the target. The extent of thesteel-cased well is noted by the vertical black line in panel (a). For the conductive tar-get, we see an accumulation of positive charges along the radial and vertical boundariesof the target. This is consistent with currents that have been channeled into the conduc-tor and exit into a more resistive background. Conversely, for the resistive target, we127see an accumulation of negative charge on the radial boundary, consistent with currentmoving from a resistive region to a more conductive material. We also notice somepositive charge accumulation on the top and bottom boundaries of the target; some ofthe currents deflected around the resistor do enter from the top and bottom, resulting inan accumulation of positive charge. This is not observed in a traditional mis-a`-la-masseexperiment, where a point source is positioned within the target. Figure 4.15 shows thecurrent density, charges and electric fields for a mis-a`-la-masse experiment in which nosteel cased well is present.In a DC experiment, the electric field response we measure is a result of the distribu-tion of charges within the domain. As a metric for quantifying excitation, I integrate thesecondary charge over this depth interval containing the target. In Table 4.1, I show thesecondary charge integrated over the depth interval containing the target; the secondarycharge on the casing within this region is included in the calculation. To examine howthe charge relates to the electric field data, I have plotted (a) total radial electric field,(b) secondary radial electric field (with respect to a primary that includes the casing in ahalfspace), and (c) the secondary radial electric field as a percentage of the primary fora down-hole source and similarly for a top-casing source (d, e, f) in Figure 4.16. I adoptthe same noise floor and percent threshold as in the casing integrity examples (10−7V/m and 20%, respectively). For time-lapse surveys where a baseline survey has beentaken and the background is well-characterized, this threshold could likely be reduced.The black line in panels (a) and (d) corresponds to the baseline model in which no tar-get is present; each of the colored lines corresponds to a different target conductivity asindicated in the legend.First, we examine the impact of the conductivity of the target and notice that there isan asymmetry between secondary charge on conductive targets and resistive targets. For128Figure 4.14: Cross section showing: (a) electrical conductivity, (b) current density,(c) charge density, and (d) electric field for a DC resistivity experimentwith a resistive target (top) and a conductive target (bottom). The positiveelectrode is positioned in the casing at the 912.5 m depth. The casing isshown by the black line that extends to 1 km depth in panel (a).Table 4.1: Integrated secondary charge over a target adjacent to the casing, asshown in Figure 4.14.integrated secondary charge (C)target conductivity (S/m) downhole source top-casing source1e-03 -4.24e-12 -1.08e-121e-02 -3.82e-12 -9.68e-131e-01 0.00e+00 0.00e+001e+00 1.75e-11 4.46e-121e+01 3.26e-11 8.28e-12129Figure 4.15: Cross section showing: (a) electrical conductivity, (b) current density,(c) charge density, and (d) electric field for a DC resistivity experimentwith a conductive target (top) and a resistive target (bottom). The positiveelectrode is positioned at 912.5 m depth. No casing is included in thissimulation. Note that the colorbars for the charge density (c) and electricfield (d) are different than those used in Figure 4.14. For the resistive target,the colorbar is saturated, the charge density over the resistive target is onthe order of 10−13 C/m3.130Figure 4.16: Radial electric field at the surface as the conductivity of a cylindricaltarget, in contact with the well, is varied. The target has a radius of 25m and extends in depth from 900 m to 925 m. The return electrode ison the surface, 500 m from the well and data are measured along a lineperpendicular to the source. The panels on the left show (a) the total electricfield, (b) the secondary electric field with respect to a primary that doesnot include the target, and (c) the secondary electric field as a percentageof the primary for a survey in which the positive electrode is positioneddownhole at 912.5 m depth. The panels on the right similarly show (d) thetotal electric field, (e) the secondary electric field, and (f) the secondaryelectric field as a percentage of the primary for a top-casing experiment.131a 1 S/m target, which is one order of magnitude more conductive than the background,the integrated secondary charge is 1.75× 10−11 C, while for a 1× 10−2 S/m target,which is one order of magnitude more resistive than the background, the integrated sec-ondary charge is −3.82× 10−12 C for the downhole casing experiment. Thus, there isa factor of 4.6 between the magnitude of the secondary charge for these targets; this isequivalent to the ratio we see between the secondary electric field measurements at thesurface observed in Figure 4.16b. When also considering the influence of the primaryelectric field on our ability to detect a target, we see that for a down-hole casing experi-ment, the conductive targets are detectable; they both have a significant region where thesecondary is above the noise floor and the secondary comprises a significant percentageof the primary. The resistive targets, however, are not. Although within 200 m of thewell, the secondary signal is above the noise floor, this also corresponds to where theprimary field is large; the percent threshold would need to be reduced to less than 5% inorder to have confidence in the signals due to the resistive targets.When comparing the downhole source to the top-casing source experiments for afixed conductivity, there is a factor of 3.9 between the integrated secondary chargeshown in 4.1; this is reflected in the secondary electric field data in Figure 4.16b &e. For the top-casing experiment, none of the targets is detectable. There are two fac-tors that make this a more challenging experiment than the downhole scenario: (1) lesscurrent is available to excite the target, as reflected in Table 4.1 and (2) the primary fieldis stronger at the receivers (200 m from the well the primary field has an amplitude of10−5 V/m, while for the down-hole source experiment, the primary has an amplitudeof 2×10−6 V/m). Addressing the excitation of the target requires that the source elec-trode be positioned downhole, closer to the target. The second point may be overcome ifreceivers can be positioned closer to the target, for example within an adjacent borehole.132In summary, the integrated secondary charge provides a metric for a survey’s abilityto excite a target, and shows that conductive targets are easier to excite than resistivetargets. As expected, if the source electrode can be positioned near the target, excitationis enhanced. This also has the added benefit of reducing the strength of the primaryelectric field at the surface, as compared to a top-casing survey; this increases the po-tential for detecting a target with surface-based receivers. In the next section, I examinethe significance of the electrical connection between the casing and the target.Target not in contact with the wellHow significant is the electrical connection between the casing and the target for ourability to excite a response? To examine this, I introduce a small gap equal to the thick-ness of the casing (1 cm) between the casing and the target. This has negligible effecton the volume of the target, but it changes the electrical characteristics of the problem.Consider a conductive target; if it is in-contact with the well, we are effectively con-ducting a mis-a`-la-masse experiment, and the conductor will have a net positive charge.When the target is isolated from the casing, the total charge on the target must be zero,and thus dipolar effects, in which negative charges build up on the inner interface of thecylinder target and positive charges build up on the outer interface of the target, will bethe source of our signal. This is demonstrated in Figure 4.17.The corresponding secondary charge integrated over the target depth and radial elec-tric field data are shown in Table 4.2 and Figure 4.18. For comparison, the data result-ing from the target in contact with the well are plotted in the dashed, semi-transparentlines. While there is little difference in the integrated secondary charge or the electricfield measurements for the resistive targets, we see that there is a factor of 1.3 difference(i.e. 30%) between the integrated secondary charges and correspondingly, the secondary133Figure 4.17: Cross section showing: (a) electrical conductivity, (b) current density,(c) charge density, and (d) electric field for a DC resistivity experimentwith a conductive target (top) and a resistive target (bottom) which is not incontact with the well. The positive electrode is positioned in the casing atthe 912.5 m depth. The casing is shown by the black line that extends to 1km depth in panel (a).electric fields, from a 10 S/m target in contact with the well versus not. Similarly, thereis a factor of 1.2 between a 1 S/m target in contact with the well versus not for both thedownhole and top-casing sources. Increasing the gap between the target and the casingdecreases the integrated charge and correspondingly reduces the secondary electric fieldat the surface. The integrated secondary charge for a 10 S/m target with a 10 cm gap be-tween the target and casing in a downhole source experiment is 1.7×10−11 C, which is134Table 4.2: Integrated secondary charge over a target that is not electrically con-nected to the casing, as shown in Figure 4.17.integrated secondary charge (C)target conductivity (S/m) downhole source top-casing source1e-03 -4.24e-12 -1.08e-121e-02 -3.80e-12 -9.64e-131e-01 0.00e+00 0.00e+001e+00 1.49e-11 3.79e-121e+01 2.51e-11 6.39e-12a factor of 2.2 smaller than the connected target; correspondingly the electric field dataat the surface are reduced by a factor of 2.2 as compared to the connected target. Thus,a direct, electrical connection between the target and the well in which we connect thesource is preferable for exciting and detecting conductive targets. In the next section, Ifurther examine the impact of the separation between the target and casing.Target offset from the wellThe examples thus far have focused on a particular geometry where the target is sym-metric about the well and is either connected or not. The more general case is wherethere is a target located anywhere in the medium and we wish to use DC or EM tocharacterize it. For example, in some scenarios, instrumenting a well for a geophysi-cal survey may not be possible if it is also actively being used for an injection. Usinganother well, offset from the injection well, may then be preferable for positioning elec-trodes. In such circumstances, the physical property model is fully 3D and there aremore factors that influence our ability to excite the target; in addition to the conductivityand geometry of the target, the distance between the well where the source electrode ispositioned is now relevant. To address these potential applications, I examine a fully 3Dscenario in which a target block is located away from the source well. Our primary goalis to compare relative excitations that arise from using different survey parameters. It135Figure 4.18: Radial electric field at the surface as the conductivity of a cylindricaltarget, which is not contact with the well, is varied. The target has a radiusof 25 m and extends in depth from 900 m to 925 m. The return electrodeis on the surface, 500 m from the well and data are measured along a lineperpendicular to the source. The panels on the left show (a) the total electricfield, (b) the secondary electric field with respect to a primary that doesnot include the target, and (c) the secondary electric field as a percentageof the primary for a survey in which the positive electrode is positioneddownhole at 912.5 m depth. The panels on the right similarly show (d) thetotal electric field, (e) the secondary electric field, and (f) the secondaryelectric field as a percentage of the primary for a top-casing experiment.The data shown in Figure 4.16, for the target in contact with the well, areplotted in the dashed, semi-transparent lines for reference.136is sufficient to evaluate an electric dipole moment of the target that is evaluated with aBorn approximation (Born, 1933).For a target offset from the pipe, we expect the secondary response due to that targetto be dipolar. Thus, a natural proxy for excitation is the dipole moment of the target. Iwill adopt a Born-approximation approach to quantify the excitation and take the normof the integrated anomalous current density over the target volume, that is,m=∣∣∣∣∣∣∣∣∫ ~ja dV ∣∣∣∣∣∣∣∣=∣∣∣∣∣∣∣∣∫ σs~ep dV ∣∣∣∣∣∣∣∣ (4.2)where σs = σ−σp is the secondary conductivity (the difference between the conductiv-ity of the target and the conductivity of the background),~ep is the primary electric field(the electric field due to the source, casing, and half-space background), and ~ja is theanomalous current density which is non-zero only over the volume where the target islocated.The target I consider is 50 m wide in both horizontal dimensions and is 25 m inheight. Its top is at 900 m depth, as in the previous examples. I will examine bothdownhole source and top-casing experiments. A depth slice showing the primary electricfield for the downhole electrode is shown in Figure 4.19; the return electrode is onthe surface at x =-500 m, y =0 m. The three different target positions relative to theborehole are outlined in white. The solid line shows the target which is inline with thereturn electrode, the dashed is 90◦ from the source electrode and the dotted line showsthe target which is 180◦ from the return electrode. I vary the distance from the wellto the target and have plotted the Born-approximated dipole moment for four differenttarget conductivities in Figure 4.20 for (a) a downhole experiment, and (b) a top-casing137Figure 4.19: Depth slice showing the primary electric field due to a downhole elec-trode and a return electrode located on the surface at x =-500 m, y =0 m.The red line indicates the azimuth of the source. We examine the 3 differ-ent target azimuths shown by the white outlines. The solid line indicatesthe target inline with the source, the dashed is 90◦ from the source line, andthe dotted line is 180◦ from the source line.experiment. The offset is calculated from the center of the well to the nearest edge ofthe target that is excited by a downhole source. As before, conductive targets are mucheasier to excite than resistive targets, and for a given conductivity, a downhole sourceprovides greater excitation than a top-casing source. Naturally, as the target is movedfurther from the well, the geometric decay of source fields reduces our ability to excitethe target. Positioning the return electrode along the same azimuth as the target acts tomitigate some of these effects for targets that are at distances greater than 200 m fromthe well, while for targets nearer to the well, the return electrode location has little effecton the excitation.The next step to consider is detection of the secondary response due to this target.Consistent with the Born-approximation approach, I simulate the target as a dipole with138Figure 4.20: Integrated anomalous current density (excitation), as defined in equa-tion 4.2, for a 50 m × 50 m × 25 m target at 900 m depth in a DC exper-iment with the positive electrode (a) downhole at 912.5 m depth, and (b) atop-casing electrode. The return electrode is positioned on the surface 500m from the well. Each line color indicates a different target-conductivity.The different line styles correspond to different target azimuths relative tothe plane of the source electrodes and correspond to those show in Figure4.19. The solid line indicates a target inline with the source, the dashed is90◦ from the source, and the dotted line is 180◦ from the source. Offset iscalculated from the center of the well to the edge of the target closest to thewell.a moment computed with equation 4.2 and compare the secondary electric field dataat the surface for models with, and without, the casing. For this example, I select themodel of a conductive target (10 S/m) with center 50 m offset from the well. The targetis along a line 90◦ from the source (e.g. along the same line as the dashed-outline inFigure 4.19). This gives a dipole moment of 38 Am for the target. The electric field data,measured along the same line as the target, are shown in Figure 4.21. The secondaryresponse with the casing is shown in blue, and the response of the same dipole in a half-space is shown in orange. The secondary response due to the dipole in a half-space fallsbelow the 10−7 V/m noise floor for all offsets, whereas, when the casing is included,the secondary response is above the noise floor until beyond offsets of 600 m from thewell. The casing not only helps excite a target, as was demonstrated in Schenkel and139Morrison (1994), it also provides a conductive pathway for the secondary currents, thusincreasing the secondary signal observed at the surface; this was similarly noted in Yanget al. (2016a). In the model with the casing, the secondary signal comprises a significantpercentage of the primary between offsets of∼ 200 m and∼ 650 m, providing a windowbetween 200 m and 600 m where the secondary electric field is above the noise floor andcomprises a significant percentage of the primary.If I move the target further from the well, positioning its center 75 m from the well,then the dipole moment is reduced to 24 Am. If I use the criterion that the secondaryelectric field must be above the noise floor and be at least 20% of the primary field,then the region over which we can expect to collect data is reduced to a 50 m windowbetween 300 m and 350 m offset from the well. For this given survey, then, we canconsider ∼ 25 Am as a threshold for detectability of a target.SummaryThe examples presented here showed that conductive targets are much easier to excitethan resistive targets. For deep targets, a downhole electrode is preferable to a top-casing source as it delivers more current at depth to excite the target and reduces thestrength of the primary at the surface; this makes the secondary field a larger percentageof the primary. For targets in close proximity to the well, if the target is in contact withthe well, that electrical connection enhances the response. Additionally, I showed thatbeyond helping excite a target, as was demonstrated by Schenkel and Morrison (1994),the casing also improves detectability of secondary signals at the surface.Designing a survey for a specific setting may require incorporation of 3D geologicstructures and may include inversions to examine a survey’s ability to recover a target.In this case, it is desirable to have a coarse-scale representation of the steel-cased well140Figure 4.21: (a) Sum of the primary and secondary radial electric field due to adipolar target with moment of 38 Am centered 50 m from the well, eithercalculated with the casing (blue) or simply a dipole in a half-space (orange).(b) Secondary radial electric field due to a dipolar target in a halfspace withcasing (blue) and without casing (orange). Secondary radial electric field asa percentage of the primary. The target is along a line 90◦ from the sourceelectrodes; this is the same line along which we measure data at the surface.141on the simulation mesh. This is the topic of the next section.4.4 Coarse-scale approximations of the wellWhen approaching the inverse problem, many forward simulations are required, andtypically, a 3D cartesian mesh, with cells that vary on the length scales of the geol-ogy, is desired. Thus, rather than performing a fine-scale simulation of the steel-casedwell, we may wish to represent the well on a coarse mesh. In the literature, two mainapproaches arise: the first approximates the well as some form of “equivalent source,”such as a charge distribution (e.g. Weiss et al. (2016)); the second approach representsthe well as a conductivity feature on the coarse-mesh (e.g. Swidinsky et al. (2013); Umet al. (2015); Yang et al. (2016a); Kohnke et al. (2017); Puzyrev et al. (2017), amongothers). Here, I will focus on the second approach, noting that a charge distributionalong the length of the well can be computed with the 2D or 3D cylindrical code de-scribed in Chapter 3. Within the literature, there is disagreement among approaches forselecting the conductivity of the coarse-scale feature approximating the well. Um et al.(2015) replaces the fluid-filled cylinder with a solid rod having the same conductivityas the casing, arguing that it is the contrast between the conductivity of the well andthe conductivity of the surrounding geology that is the most important factor; Puzyrevet al. (2017) also adopts this approach. Other authors have opted to preserve the cross-sectional conductance of the well (Swidinsky et al., 2013; Kohnke et al., 2017); this isconsistent with the transmission-line model of the well discussed in Kaufman (1990).The aim of this section is to analyze these approaches.1424.4.1 Replacing a hollow-cased well with a solid cylinderI consider a steel-cased well with a conductivity of 5× 106 S/m that is embedded ina 0.1 S/m halfspace; the conductivity of the material that fills the well is the same asthe background. The well has an outer diameter of 10 cm and a thickness of 1 cm,and I will vary its length. I will perform a top-casing experiment, where the positiveelectrode is connected to the casing at the surface. The return electrode is positioned8 km away, and a cylindrically symmetric mesh is used in the simulations. I examineapproximations that treat the casing as a solid cylinder with the same outer-diameter asthe true, hollow-cased well.The distribution of charges, or equivalently, the current in the casing, is the sourceof the electric response of the casing. Thus to judge if two models of the casing are“equivalent”, I examine the current and charges as a function of depth. In Figure 4.22,I have plotted the vertical current and charges along the casing for the true, hollowcased well (solid), solid cylinder with conductivity equal to that of the casing, 5× 106S/m (dashed), and solid cylinder with a conductivity that preserves the product of theconductivity and the cross-sectional area of the conductor, 1.8× 106 (dotted), for fourdifferent casing lengths, each indicated by a different color. Figure 4.22 shows: (a) thevertical current along the casing, (b) the difference in current between the approximatemodel and the true model, (c) that difference as a percentage of the true solution (d)the charge per unit length, (e) difference in charge per unit length and (f) difference incharge per unit length as a percentage of the true solution.For short wells, we see that the current decays linearly and that the charge distri-bution is nearly uniform above the end of the well, while for longer wells, the decayof the current is exponential in nature, as is the charge distribution. This behavior isconsistent with that predicted by the transmission line solution described in Kaufman143Figure 4.22: Currents (top row) and charges (bottom row) along the length of ahollow steel-cased well (solid lines), solid cylinder with conductivity equalto that of the steel-cased well (dashed-lines), and a solid cylinder with aconductivity such that the product of the conductivity and the cross sec-tional area of the cylinder is equal to that of the hollow-pipe (dotted lines).Each of the line-colors corresponds to a different casing length, as indi-cated in the legend. In (a), we show the vertical current in the casing, (b)shows the difference from the true, hollow-cased well in the vertical currentwithin the casing, and (c) shows that difference as a percentage of the truecurrents. In (d), we show the charge per unit length along the casing, (e)shows the difference from the true, hollow-cased well and (f) shows thatdifferences as a percentage of the true charge distribution. The x-axis on allplots is depth normalized by the length of the casing.and Wightman (1993). Kaufman and Wightman (1993) showed that the transition be-tween the linear decay of currents and the exponential decay of currents is controlled bythree factors: the cross sectional conductance of the well, the resistivity of the surround-ing formation, and the length of the well. Schenkel (1991) similarly summarized thisbehavior in the definition of the conduction length (equation 4.1), which is the lengthover which the currents in the casing have decayed by a factor of 1/e. For sufficientlyconductive and short wells (e.g. Lc/δ  1, where Lc is the length of the casing), thecurrent decay is linear and independent of the conductivity, whereas for longer wells,144(Lc/δ  1), the rate of decay of the currents is controlled by the conduction length (seeequations 45 and 53 in Kaufman and Wightman (1993)).In preserving the cross-sectional conductance, we see that the difference in cur-rents and charges along the length of the well is negligible; the maximum differencein currents for the 2000 m long well which has equivalent cross-sectional conductanceis 7× 10−7 A as compared to the difference of 0.18 A when using the conductivity ofthe casing. This difference is important as it changes how much current is available toexcite a target at depth. For a 2000 m long well, the current is overestimated by > 150%if the well is replaced by a solid cylinder with the same conductivity of the steel-casedwell. It also changes the distribution of charges and thus the electric field due to thewell. Figure 4.22e shows us that the extra conductance introduced when approximatingthe well using the conductivity equal to the casing results in a secondary dipolar chargeon the casing. This in turn reduces the electric field we observe at the surface, as shownin Figure 4.23. For a long well, the difference can be as large as 40% near the well.The numerical time-domain EM experiment used in Um et al. (2015) to justify theapproximation of the well by a solid, conductive rod having the same conductivity as thesteel-cased well used a 200 m long well with a thickness of 12.223 mm, outer diameterof 135 mm, conductivity of 106 S/m in 0.033 S/m half-space. The conduction lengthof this well is 560 m; this is more than twice the length of the well. Therefore, thebehavior of the currents falls into the linear regime, where the decay of currents is mostlyindependent of the conductivity, and thus the difference between using the conductivityof the casing or preserving cross-sectional conductance is less significant. However, iflonger wells such as those typically employed in hydrocarbon settings, are considered,the behavior of the currents and charges depends upon the conductance of the casing,and thus that is the quantity that should be conserved in an approximation of the hollow-145Figure 4.23: Radial electric field measured at the surface for a model of a hollowsteel-cased well (solid lines), a solid cylinder with conductivity equal tothat of the steel-cased well (dashed-lines), and a solid cylinder with a con-ductivity such that the product of the conductivity and the cross sectionalarea of the cylinder is equal to that of the hollow-pipe (dotted lines). Eachof the line-colors corresponds to a different casing length, as indicated inthe legend. In (a), we show the total radial electric field, (b) shows the dif-ference in electric field from that due to the true, hollow-cased well, and (c)shows that difference as a percentage of the true electric fields. The x-axison all plots is distance from the well normalized by the length of the casing.cased well by a solid rod.In order to confirm that this conclusion is valid for variable geology, I have includeda simulation with a 2 km long casing in a layered background. Each layer is 50 m thickand the conductivity was assigned randomly; three instances are included, as shown inFigure 4.24. The mean of the background conductivity is 0.1 S/m for each of the models.The currents and charges along the length of the well for the true model, and a modelapproximating the well as a solid cylinder with equal cross-sectional conductance, areshown in Figure 4.25. For all of the models shown, the difference in both the casingcurrents and the charges are 5 orders of magnitude less than the amplitude of the totalcurrents and charges; thus I conclude that approximating a hollow cylindrical steel cas-ing by a solid cylinder with a conductivity that preserves cross-sectional conductance isvalid for models with variable geology.146Figure 4.24: Three realizations of a 2 km long casing in a layered background,where the conductivity of the layers is assigned randomly. Each layer is50 m thick, and the mean conductivity of the background is 0.1 S/m. Thecolor of the title corresponds to the plots of the currents and charges inFigure 4.254.4.2 Cartesian gridIn the previous section, I showed that a hollow, cylindrical steel-cased well can be ap-proximated by a solid cylinder with equal cross-sectional conductance. In this section,I move to a coarser, cartesian mesh, such as might be employed when solving a 3D in-verse problem. I examine a simple approximation of a steel cased well on a cartesiangrid. I employ 4 tensor meshes, each with progressively larger cell widths for the finestcells that capture the casing. On each of the cartesian meshes, I approximate the casingby preserving the product of the conductivity and the cross sectional area on the mesh.For comparison, I run a fine-scale simulation on a 3D cylindrical mesh that accuratelydiscretizes the casing; it uses 4 cells across the casing-wall. The casing model is similarto that used in previous examples: it is 1 km long, has an outer diameter of 10 cm, athickness of 1 cm, and is embedded in a 0.1 S/m half-space. The positive electrode is147Figure 4.25: (a) Total vertical current through the casing for the three layered-earthmodels shown in Figure 4.24. The solid lines indicate the response of thetrue, hollow steel cased-well and the dotted lines indicate the response ofa solid cylinder having the same cross-sectional conductance as the hollowwell. (b) Difference between the currents along the casing in the solid wellapproximation and the true, hollow well. (c) Charge per unit length for eachof the models. (d) Difference in charge per unit length between the truemodel of the casing and the approximation which preserves cross-sectionalconductance.connected to the top of the casing and a return electrode is positioned 1 km from thewell-head.The resultant currents and charge per unit length are shown in Figure 4.26. In the toprow, panel (a) shows the total current in a region approximating the well, along with thetotal current in the “true” cylindrical well (black line), (b) shows the difference betweenthe current through the cartesian cells and the true model, and (c) shows the differenceas a percentage. Similarly, in the bottom row, I show (d) the charge per unit length alongthe cylindrical well (black line) and cartesian-prism approximations, (e) the difference148Figure 4.26: Currents (top row) and charges (bottom row) along the length of asteel cased well. The “true” hollow-cased well is simulated on a 3D cylin-drical mesh and has 4 cells across the width of the casing thickness (blackline). The colored lines correspond to the currents and charges computedalong the well represented on a cartesian mesh with cell widths shown inthe legend. The finest vertical discretization is 2.5 m in all simulations.To represent the hollow cased well on the cartesian mesh, the cells inter-sected by the casing are assigned a conductivity that preserves the productof the conductivity and cross-sectional area of the well. In (a), we showthe vertical current in the casing, (b) shows the difference from the true,hollow-cased well in the vertical current within the casing, and (c) showsthat difference as a percentage of the true currents. In (d), we show thecharge per unit length along the casing, (e) shows the difference from thetrue, hollow-cased well and (e) shows that differences as a percentage ofthe true charge distribution.in charge per unit length from the charge per unit length on the true cylindrical model,and (f) that difference as a percentage of the charge per unit length on the cylindricalwell.The approximation of the cylindrical well by a rectangular prism with width equalto the diameter of the casing introduces minimal error in the currents and charges com-puted using a finite volume approach, even though the casing is only captured by one cell149across its width. Comparing the current along the length of the well for the 3D cylindri-cal well and the cartesian simulation with 0.1 m cells, we see that the error introducedis < 2.5% (until the end of the well where the current approaches zero). Similarly,the difference in the charge per unit length is < ±1.25%. As successively coarser dis-cretizations are used, accuracy is gradually lost; by doubling the cell sizes to 0.2 m, theerror in the currents is 6% at its maximum and < ±3% in the charge along the casing.A factor of 8 increase in cell size (0.8 m cells) results in a maximum error of 15% inthe currents. It is important to note that the forward simulation is conducted using afinite volume approach; other approaches such as finite difference or integral equationapproaches may have worse agreement if care is not taken to handle large physical prop-erty contrasts, captured by a single cell, in the simulation. Note that the behavior of theerrors depends upon the properties of the casing (e.g. conductivity and length) as wellas the conductivity of the background. This might be expected from the description ofthe casing conduction length (equation 4.1). If the conduction length is large relativeto the length of the well, the currents decay linearly, and the geometry and conductivityof the well are less significant in the behavior of the currents. Alternatively, if the con-duction length is comparable to the length of the well, the currents decay exponentiallywith a decay rate that depends on the geometry and conductivity of the well. For exam-ple, if the background is more resistive, increasing the contrast between the casing andbackground, the errors are reduced. Using a background conductivity of 100Ωm, themaximum error introduced in the current is < 1% with 0.1 m cells and < 2% with 0.8m cells.Depending on the level of accuracy required in a 3D simulation, there are severalstrategies that one might take to reduce this error. In some cases, local refinement can beachieved with a tetrahedral mesh, as is often employed when using finite element tech-150niques (e.g. Weiss et al. (2016)), or an OcTree mesh (Haber et al., 2007). Other, moreadvanced approaches including upscaling and multiscale could also be considered. In anupscaling approach, one inverts for a conductivity model, which might be anisotropic,that replicates the physical behavior of interest (Caudillo-Mata et al., 2017a). Multiscaletechniques translate conductivity information from a fine-scale mesh to a coarse-scalemesh, on which the full simulation is to be solved, using a coarse-to-fine interpolationthat is found by solving Maxwell’s equations on the fine mesh locally for each coarsegrid cell (Haber and Ruthotto, 2018; Caudillo-Mata et al., 2017b). The 3D cylindricalforward simulation code described in Chapter 3 and used in this example can serve atool for validating and refining an approach to achieve the desired level of accuracy.4.5 ConclusionThe work in this chapter is motivated by the increasing use of steel cased wells in geo-science problems, including monitoring applications such as carbon capture and storageand hydraulic fracturing. For geophysical imaging of targets at depth, the wells arebeneficial as they can be used to channel currents to depth and enhance signals at thesurface for targets that otherwise would be undetectable from a surface-based survey.Additionally, there is interest in considering the casing itself as the geophysical targetin casing integrity experiments; here the aim is to detect flaws or breaks in the casing.These applications, coupled with advances in modeling capabilities, open up the poten-tial for advancing the utility of electrical and electromagnetic imaging in settings withmetallic-cased wells.Despite this potential, the reality is that electric fields, especially if measured at theEarth’s surface, are small. Secondary fields might only be a few percent of the primaryfield, and thus too insignificant to reliably detect the target of interest. The success of151using a DC or EM survey then depends upon many details that pertain to understandingthe basic physics, the effects of parameters of the casing, the background conductivity,location of the current electrodes, and discerning which fields should be measured. DCresistivity is the starting point, as it allows us to examine the currents, charges, andelectric fields in the electrostatic limit, prior to introducing inductive effects and theinfluence of magnetic permeability in an EM signal, and as such, was the focus of thispaper. Regarding the physics, a DC survey involves attaching a current generator toa conductive medium. This establishes a steady state current; the signal to which weare sensitive is the electric field that arises from charges that accumulate at interfacesseparating regions of different conductivity. For this reason, most of our results are firstpresented as currents and charges.The large contrasts in physical properties and significant variation in length scalesdue to long, thin, cylindrical, steel-cased wells prompt a number of questions about howthe EM fields behave. In many cases, the finer details about the physical responses haschallenged my intuition. With respect to the casing integrity application there were basicquestions: how does a flaw in the pipe affect the currents and electric fields measuredat the surface? Does the extent of the flaw change our ability to detect it (e.g. if ithas a vertical extent of several meters versus a vertical extent of centimeters)? Whathappens if the flaw only comprises a part of the well, leaving some connectedness inthe casing? When considering a geophysical experiment for imaging a target: is there asignificant difference in the currents at depth between scenarios where a source electrodeis connected to the well-head at the surface and one where the source electrode is offsetfrom the well by a few meters? The surface signals are small; are there preferentialgeometries for the source and receiver electrodes that maximize signal/noise? A majorgoal of the DC survey will be to excite and detect target bodies. For problems, such152as CO2 sequestration, enhanced oil recovery, or hydraulic fracturing, the target may ormay not be in contact with the well; how significant is an electrical contact between atarget and the well in the data we measure at the surface?Looking towards solving inverse problems in settings with steel-cased wells, it isadvantageous to reduce the computational cost of the forward simulation because aninversion requires many forward simulations. There are several approaches that can betaken to achieve this; one common approach is to approximate the finely-discretizedwell with an approximation on a coarse-scale mesh. Currently, there is disagreementwithin the literature as to how this should be done. Some authors have advocated that theconductivity contrast between the casing and the background should be preserved andthus replace a hollow steel-cased well with a solid rod that has a conductivity equal tothat of the casing. Other authors have opted to preserve the product of the conductivityand cross-sectional area of the casing, following the conclusions of the transmission-line solution shown in Kaufman (1990). What is the correct conductivity needed forsubstitution?Some of the above questions have been addressed in theoretical papers extendingback a few decades but numerical verification was often limited or carried out withsimplifying assumptions. Other questions require the ability to carry out numericalmodeling in 2D or 3D environments – these tools are just now becoming available. Mygoal with this chapter has been to examine the scientific questions above and to promoteinsight about the solution by plotting the currents, charges, and electric fields.153Chapter 5Electromagnetics with steel cased wells5.1 IntroductionIn the previous chapter, I focussed on the physics of direct current (DC) resistivity ex-periments in settings with steel cased wells. This provided a fundamental understandingin terms of the charges, currents and electric fields. In this chapter, I shift attention tothe electromagnetic (EM) problem where the fields and fluxes vary in time. EM canbe advantageous in that for a given survey geometry, the response can be sampled atmultiple frequencies or times, providing richer information than is available by onlythe electrostatic response in a DC experiment. The richness in information comes bothfrom the additional physics introduced as inductive processes become relevant, as wellas potentially multiple excitation directions as the formation currents change throughtime.In addition to time-variation of fields and fluxes in an EM experiment, the physics insettings with steel casing is further complicated by the significant magnetic permeabilityof the steel. Magnetic permeability is often neglected in EM simulations and inversionsas the variations in permeability are typically much less significant than the role of con-154ductivity in the data. However, steel has a permeability of ≥ 100 µ0 (Wu and Habashy,1994).The role of permeable casings has been explored for inductive sources and receivers,primarily in the context of cross-well EM surveys (Augustin et al., 1989a; Wu andHabashy, 1994; Wilt et al., 1996; Becker et al., 1997; Lee et al., 1998, 2005; Kim andLee, 2006). For grounded sources, the influence of permeability is much less explored.More recently, some authors have included magnetic permeability in simulations withcasings, (e.g. (Kong et al., 2009)), but most of the published EM simulations of metal-lic cased wells neglect to include magnetic permeability (e.g. (Swidinsky et al., 2013;Commer et al., 2015; Um et al., 2015; Fang et al., 2017; Puzyrev et al., 2017)). Withinthe publications that include permeability in grounded-source simulations, there is min-imal analysis on how permeability alters the behavior of the currents in the earth or howit affects data measured at the surface. A notable exception is the work in Wait andWilliams (1985); Williams and Wait (1985). They developed an analytical model fora dipole-dipole frequency domain electromagnetic experiment over a halfspace with aconductive, permeable, and polarizable casing. Their simulations showed that if the cas-ing is permeable, the apparent resistivity and phase measurements made at the surfacedata can be biased upwards as compared to only a conductive casing. Expanding on ourunderstanding of how the permeability impacts the subsurface currents and the resultantmeasured data in a grounded-source experiment with a steel cased well is a prime goalof this chapter.For numerical modelling, particularly when considering the inverse problem, it is ad-vantageous to reduce the size of the mesh to reduce the computational cost. In Chapter4, I showed that for a DC resistivity experiment, preserving the product of the conduc-tivity and the cross-sectional area of the casing is a viable strategy if the cells which155encapsulate the casing are approximately equal in size to the diameter of the casing. Itis questionable if this same approximation holds for EM experiments. Furthermore, it isunclear how to approximate the magnetic permeability on a coarse scale. Schwarzbachand Haber (2018) suggests an upscaling strategy based upon that presented in Caudillo-Mata et al. (2017a), where one inverts for the conductivity and the permeability of thecasing. By analyzing the behavior of the fields and fluxes, we may be able to calibrateour expectations of upscaling strategies for permeability.We have the choice to focus our analysis in the time-domain or the frequency-domain. Here, I choose to conduct the analysis in the time-domain. The physics ofEM are arguably more intuitive in the time-domain, as all fields are real and we can stepthrough changes with time. The frequency-domain requires that we consider the parti-tioning of electromagnetic energy into real and quadrature components; this additionalstep can be a hurdle to building intuition.The chapter is organized as follows. In section 5.2, I consider the time-domain EM(TDEM) response of a conductive well in a “top-casing” experiment where one elec-trode is connected to the top of the well and a return electrode is positioned some dis-tance away on the surface. I examine the currents within the casing and in the surround-ing geologic formation through time; this provides the foundation for understanding theEM casing response. Following from this, I introduce magnetic permeability in Section5.3 and demonstrate how it alters the current and the impact this has on data collected atthe surface. Finally, in section 5.4, I examine the the problem of approximating a con-ductive, permeable, steel-cased well and discuss some of the challenges that are uniqueto EM as compared to DC.Source code for all of the examples shown in this chapter are available as Jupyternotebooks as outlined in Appendix A.1565.2 EM response of a conductive casingAs a starting point for this discussion, I first perform a simple numerical experiment todemonstrate the behavior of the currents in a time-domain EM experiment over a half-space. For this simulation, a positive electrode is positioned at the origin (where thecasing will be) and a return electrode is positioned 1 km away (at x =1000 m). Theconductivity of the background is 10−2 S/m, and the conductivity of the air is set to1× 10−4 S/m. When we come to simulations with the casing, I have observed thatcontrasts near or larger than ∼ 1012 S/m leads to erroneous numerical solutions, andthus it is preferable to set the air to be more conductive than what is typically employedin EM simulations (often ∼ 10−8 S/m). A step-off waveform is used and the currentswithin the formation are plotted through time in Figure 5.1. Panel (a) shows a verticalcross section along the line of the wire, (b) shows a horizontal depth slice at 50 m depthand (c) shows a depth slice at 800 m depth. The images in panels (a), (b) and (c) are onthe same color scale.At time t =0 s, the currents behave according to the DC solution. In panel (a), wesee the flow of current from the source electrode, through the formation to the returnelectrode. Immediately after the source-current is shut-off, an image current developsin the formation; this image current acts to preserve the magnetic flux initially set-upby the source current. The image current flows in the same direction as the originalcurrent in the wire. This is opposite to currents in the formation, causing a circulationof current. The center of this circulation is visible as the null at x=500 m (the mid-waypoint between the two current electrodes) which propagates downwards through time.Similarly, in the depth slices, we can see the image currents diffusing downwards andoutwards with time. For example, between 1 ms and 5 ms, the image current passes 800m depth.157Figure 5.1: Current density for a time domain experiment over a 10−2 S/m half-space. The positive electrode is on the surface where the well-head will beand the return electrode is at x= 1000 m. Each row corresponds to differenttime, as indicated in the plots in panel (c). Panel (a) shows a cross sectionthrough the half-space along the same line as the source-wire. Panels (b) and(c) show depth-slices of the currents at 54 m and 801 m depth.158Having established the background response, I now consider a model with a steel-cased well. The well has a conductivity of 5×106 S/m and is 1 km long; it has an outerdiameter of 10 cm and a 1 cm thick casing wall. The mud which infills the well hasthe same conductivity as the background, 10−2 S/m. At this point, I focus on electricalconductivity only and set the permeability of the well to µ0. I will introduce variablemagnetic permeability in the following section. The source electrode is connected tothe top of the casing, and the return electrode is in the same position as the previousexample (x =1000 m). The same step-off experiment is run and the currents plotted inFigure 5.2. I have added a fourth panel (a) which zooms in to the wellbore. Panels (b),(c) and (d) again show a cross-section of the currents in the formation and depth slicesat 50 m and 800 m depth, respectively.At time, t =0 s, we can see the increased current density along the length of the well,which correspondingly increases the current density deep within the formation. At theDC initial condition, currents flow away from the well, and eventually curve back to thereturn electrode. As in the previous example, an image current develops in the forma-tion immediately after shut-off, which can be seen in Figure 5.2 (b). There is again acirculation of current, however the geometry of this circulation and the correspondingnull is more complex than in the previous example. In Figure 5.2 (a), we see the back-ground 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.01ms, the zero crossing, which distinguishes the depth between incoming and outgoingcurrent in the casing, occurs at z=90 m, at t =0.1 ms it is at 225 m and by t =1 ms, thezero crossing approaches the midway point in the casing and is at 470 m depth. At latertimes, the downward propagation of this zero-crossing slows as the image currents arechanneled into the highly conductive casing.159Figure 5.2: Current density for a top-casing time domain EM experiment with aconductive well (5× 106 S/m). (a) Cross section in the immediate vicinityof the well. (b) Cross section through the formation. (c – d) depth-slices.160Further out in the formation, we similarly see the interaction of the diffusing for-mation currents set up at DC and the image currents. Rather than being centered at themidway point between the positive and negative electrodes as in the half-space example(Figure 5.1), the center of the circulating currents shifts with time. At early times, it iscloser to the well, while at later times, it is closer to the mid-way point. By 10 ms, theimpact of the well is much less evident and the currents are very similar in character tothose obtained in the half-space experiment.On the side of the well opposite to the wire, we also see a null develop; it is visiblein the cross sections in panel (a). To help understand this, we examine the depth slicesin panel (c). Behind the well, we see that as the image currents diffuse downwards andoutwards, some of those currents are channeled back towards the well; this is visible inthe depth slice at 10−4 s. These channeled currents are opposite in direction to those theformation currents set up at t =0 s, which also are diffusing downwards and outwards;where these two processes intersect, there is a current shadow.We can isolate the “casing response” by taking the difference between the currentsgenerated when the casing is present (Figure 5.2) and those generated in the half-spacemodel (Figure 5.1). This is plotted in Figure 5.3. Within the casing and the surround-ing formation, these secondary currents are cylindrically symmetric. This might not beimmediately intuitive given the complexity of the behavior observed in Figure 5.2, how-ever, it is a highly conductive, cylindrically symmetric target, and the main impact it hason the physical behavior is that it channels currents along its length. At all times, thereis a downward-going secondary current within the casing. This results in dipolar-likecurrents within the surrounding formation.161Figure 5.3: Difference in current density for a time domain experiment which in-cludes a conductive steel-cased well (as in Figure 5.2) and an experimentover a halfspace (as in Figure 5.1).1625.2.1 DataOn the surface, we can measure electric field data and/or the time-rate of change ofthe magnetic field (db/dt) as a function of time. Due to the cylindrical symmetry ofthe currents in Figure 5.3, we can expect that the radial electric field and the tangentialdb/dt will be most impacted by the presence of the casing.Electric fieldIn Figure 5.4, I show surface electric field data for 2 times, an early time (0.01 ms)and a later time (5 ms). The top row shows data from the background model. The twopanels on the left show the electric field at the surface at (a) 0.01 ms and (b) 5 ms. On theright, (c) shows the radial electric field measured along the x=0 m line (indicated by thewhite-dashed lines in (a) and (b)), and (d) shows a time-decay for a receiver positioned300 m from the well (indicated by the red dot in (a) and (b)). Data that are positive areplotted with the solid lines and data that are negative are plotted with dashed lines. Thered dots in panel (c) correspond to the 300 m offset location where the data in (d) aretaken from. Similarly, the blue and green dots in (d) correspond to the time-channelsshown in (c). The middle row shows the same information but for the simulation withthe conductive well, and the bottom row shows the difference between the casing andbackground data (casing minus background).At early times (0.01ms), the behavior of the geometry of the electric field is fairlycomplex (panel a). The electric field we observe results from a combination of thediffusing galvanic and image currents. The galvanic currents, set up at DC, are dipolarin nature, with field lines directed radially outward from the positive electrode locatedat x =0 m. The image current is directed in the negative x-direction and causes someinward deflection of the fields (this can also be observed at 0.01 ms in Figure 5.1). A163Figure 5.4: Simulated electric field at the surface of the earth at (a) 0.01 ms and(b) 5 ms after shut-off for the halfspace model. (c) Radial electric field datameasured along the survey line shown in white in (a) and (b) at 0.01 ms(blue) and 5 ms (green). (d) Radial electric field as a function of time at300 m along the survey line (shown in the red dot in (a)). The red dots in(c) correspond to the data observed at 300 m offset and the blue and greendots in (d) correspond to the 0.01 ms and 5 ms data. Similar information isshown in (e), (f), (g) and (h) for the model with the conductive casing. Thedifference in the radial electric field data (casing minus halfspace) is shownin (i), (j), (k) and (l).164Figure 5.5: Sketch of the plan-view geometry of the early-time galvanic currents(left) and image currents (right) in a time-domain EM experiment over a half-space. Through time, both current systems diffuse downwards and outwardsas they decay. The wire follows a straight line between the negative andpositive electrodes. The green dashed line shows where we are measuringthe radial electric field data. Prior to shut-off, current in the wire flows fromthe negative to the positive electrode. In the earth, the galvanic currentsare dipolar in nature and flow from the positive to the negative electrode.Along the survey-line, the radial component of the galvanic currents alwayspoints outwards. Immediately after shut-off, image currents are induced inthe earth. They are oriented in the same direction as the original currentin the wire and are directed away from the negative electrode towards thepositive. Along the survey-line, the radial component of the image currentsis always pointed inwards.sketch of both current systems is shown in Figure 5.5The radial electric field data at offsets less than 300 m are primarily influenced by thegeometry of the image currents, while beyond 300 m, the data are primarily influencedby the dipolar galvanic currents. Diffusion of both current systems continues throughtime. By 5 ms the image current has passed and the geometry of the fields reflects thatof the diffusing galvanic currents. This means that the electric fields are a diffuse dipolarshape that is centered at the midpoint of the wire (e.g. as in Figure 5.1b at 1 and 5 ms).165Along our survey line, the electric field is oriented mainly in the negative y-direction, butwith a slight outward deflection due to the dipolar nature of the fields, and are thereforepositive. Through time, we can similarly see the impact of the passing image current.At very early times, the data are positive, as the galvanic currents and associated electricfield are the main contributor to the response. When the image current arrives, just past0.01 ms, the electric field is deflected inwards, reversing the sign of the radial electricfield data. The second sign-reversal occurs at ∼ 0.3 ms; at this point, the image currenthas diffused considerably, and we again return to a geometry that is mainly influencedby the diffusing galvanic currents.I now consider the data when a casing is present. The casing introduces furthercomplexity as it channels both galvanic and image currents. We can see the impact inthe geometry of the fields in Figure 5.4(e). The channeled currents are directed inwards,towards the well, giving the large negative radial electric field response observed atshort offsets in (g). In panel (g), we see that this additional current-channeling pushesthe sign reversal that we observed in the early time in (c) for the background model to afurther offset. At later times, we continue to see the impact of the casing as the diffusingcurrents are deflected inwards (panel f). The result is that the radial electric field dataremain negative through time as can be seen in (g) and (h).The story simplifies if we examine the difference between the data from the modelwith the casing and those simulated with the background. Essentially, by subtractingthe background data from the casing data, we remove the 3D survey geometry from thebehavior of the fields, and are left with the current-channeling behavior. The well is ahighly conductive feature and once the source current is shut off, the diffusing imagecurrent and galvanic currents are channeled towards it. The result is the cylindricallysymmetric fields shown in (g) and (h). The electric field is directed radially inwards,166and thus the secondary “casing-response” in (k) and (l) is always negative through time.Time rate-of-change of the magnetic fluxWe can examine the db/dt data at the surface. Figure 5.6, shows the db/dt data for thebackground (top row), casing model (second row), and difference between them (thirdrow). Additionally, I include plots of the db/dt data difference as a percentage of thebackground response in (m) and (n). The data plotted in the rightmost column (plots(c), (d) and below) are the azimuthal component of db/dt along the white-dashed linesshown in panels (a) and (b).The db/dt response is quite 3-dimensional in its geometry. At early times (panel(a), 0.01 ms), we have a rotational component that arises from the downward-directedcurrent density at the positive electrode, located at the origin. We also have a con-tribution from the horizontal image currents and galvanic currents, which complicatesthe behavior between x =0 m and x =1000 m. Near the surface, there is a y-directedcomponent from these current systems. Later in time (5 ms) the db/dt data in (b) areprimarily influenced by the diffusing horizontal currents and thus db/dt is oriented inthe y-direction along the wirepath. The difference between the data from the simulationwith the casing and without is much less dramatic than we observed in the electric fielddata. It is really only the late times at short offsets that are significantly impacted; be-yond 400 m offset from the well, the difference between the db/dt data with the casing(f) and without (c) is less than 10% for the 1 ms data. In contrast, the electric field datain Figure 5.4 differed by several orders of magnitude at all offsets shown.To summarize the above results, we see that each data-type is sensitive to differentaspects of the physics. The radial electric field data are influenced by radial currents,which, as seen in Figure 5.3, are significantly altered near the surface over a large range167Figure 5.6: Simulated db/dt at (a) 0.01 ms and (b) 5 ms after shut-off for thehalfspace model. (c) Tangential db/dt measured along the white line in (a)and (b) at 0.01 ms (blue) and 5 ms (green). (d) Tangential db/dt as a functionof time at 300 m along the survey line (shown in the red dot in (a)). Similarinformation is shown in (e), (f), (g) and (h) for the model with the conductivecasing. The difference in the db/dt data (casing minus halfspace) is shownin (i), (j), (k) and (l). The difference is also shown as a percentage of thehalfspace solution at 0.01 ms and 5 ms in (m) and through time at 300 moffset in (n).168of offsets when the casing is present. The tangential db/dt data, however, are primarilyinfluenced by vertical currents, which are most altered near the wellbore.5.2.2 DiscussionThere are a number of points to highlight in these examples. The first, which has beennoted 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 mech-anisms: (1) at DC, prior to shut-off, the casing acts as an “extended electrode” leakingcurrent into the formation along its length; (2) after shut-off, it channels the image cur-rents and increases the current density within the vicinity of the casing. The secondpoint to note is that there are several survey design considerations raised by examiningthe currents: targets that are positioned where there is significant current will be mostilluminated. If the target is near the surface and offset from the well, a survey wherethe source wire runs along the same line as the target will have the added benefits of theexcitation 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 amplitudeand direction of the currents with time generate different excitations of the target. Thisshould provide enhanced information in an inversion, as compared to a single excitationthat is available from a DC survey. For deeper targets in this experiment, the passingimage current has diffused significantly, and thus it appears that the wire location hasless impact on the magnitude of the current density with location. However, it is possi-ble that increasing the wire length could be beneficial. This extension is straightforwardand could be examined with the provided script. Often, little attention is paid to thewirepath in an EM survey, and only the electrode positions are considered as a part ofthe survey design. These examples show that the image currents, whose geometry is169dictated by the geometry of the wirepath, have a significant influence on the behaviorof the currents and the resultant data. Thus, these examples also provide motivation forthinking of the wirepath as a component of the EM survey design.5.3 EM response of permeable casingsThe previous section demonstrated the behavior of the currents in a TDEM experimentwith a conductive casing; in this section, I am interested in building an understandingof how the magnetic permeability of the pipe affects the currents and the EM response.Using the same survey setup and casing geometry as in Figure 5.2, I now include mag-netic permeability in the simulation. The permeability of the casing is set to µ = 100µ0and the resultant currents are shown in Figure 5.7.The early-time behavior of the currents is very similar to what we observed in thecase of a conductive well (Figure 5.7). As we move to later times, however, we cansee that the permeability of the steel introduces further complexity into the behavior ofthe currents. Notably, in Figure 5.7(a) at 5ms and 10ms, there is a vertical circulationof currents within the well-casing; near the inner wall of the casing, the currents aremoving downwards, while near the outer wall, the currents are moving upwards. This isparticularly noticeable near 750m depth in the 5ms image. The behavior of the currentswithin the formation is also more complex. By 10ms, the impact of the conductivecasing was minimal on the behavior of the formation currents, whereas here, we see thatthe permeable and conductive casing still significantly alters the formation currents (ascompared to Figure 5.1b).To isolate the influence of magnetic permeability, I take the difference between thecurrents observed in the example with the conductive, permeable well (Figure 5.7) andthose observed in the conductive well (Figure 5.2). The result is shown in Figure 5.8.170Figure 5.7: Current density for a top-casing time domain EM experiment over a10−2 S/m half-space and a 1 km-long, conductive, permeable steel-casedwell (5×106 S/m, 100µ0).171Figure 5.8: Difference in current density for a time domain experiment which in-cludes a conductive, permeable steel-cased well (as in Figure 5.7) and anexperiment over a conductive well (as in Figure 5.2).172At the DC limit, magnetic permeability has no influence on the behavior of the cur-rents. However, as we move to later times, EM effects become relevant and we beginto see the impact that permeability has on the behavior of the currents. The genera-tion of poloidal currents within the casing wall is not an effect that can be caused byconductivity alone. Often in EM problems the magnetic field,~h and the magnetic fluxdensity ~b are treated interchangeably, which is sensible when µ = µ0 everywhere inthe domain. However, when variable permeability is introduced, the distinction is im-portant. Figure 5.9 demonstrates how the vertical circulation of current can arise dueto the interplay between the currents and corresponding rotational magnetic fields andfluxes. The downward directed current density within the casing (a) generates a toroidalmagnetic field everywhere in space (b) according to Amperes law. The magnetic flux isrelated to the magnetic field through the constitutive relation, and since the pipe is muchmore permeable than the background, the magnetic flux density is concentrated withinthe pipe. Once the transmitter is shut-off, the magnetic flux decreases, and the timevariation of the magnetic flux generates rotating electric fields (d) according to Faradayslaw. Finally, the electric field is related to the current density through Ohms law, andtherefore the rotating current density is concentrated within the pipe (e). If the pipe isonly conductive and not permeable, then there is no concentration of the magnetic flux,as in (c). Instead, the magnetic flux has identical geometry to the magnetic field, and bysymmetry, the rotation of the electric field cancels, leaving only a vertical component.5.3.1 DataAs in Figure 5.3, the secondary currents due to the permeability of the well are cylin-drically symmetric. In the data, we then expect that the radial electric field will beminimally impacted at early times, but at later times, we will see the influence of per-173Figure 5.9: Sketch demonstrating how a vertical circulation of current can ariseinside of a permeable casing in a top-casing TDEM experiment. A sourcecurrent is applied and (a) currents flow downwards through the pipe. (b)Currents generate rotational magnetic fields according to Ampere’s law. (c)Magnetic flux is concentrated in the permeable pipe according to the con-stitutive relation between ~b and ~h. (d) The magnetic flux varies with timeafter shut-off, and the time-varying magnetic flux creates rotational electricfields according to Faraday’s law. (e) Currents associated with those elec-tric fields are concentrated in conductive regions of the model in accordancewith Ohm’s law.meability. Figure 5.10 shows the radial electric field data for the conductive casing (toprow), and permeable casing (second row). As well as the difference (permeable minusconductive) and difference as a percentage of the data simulated with the conductivewell in the third and fourth rows, respectively. At early times, the impact of permeabil-ity is minimal. For example, at 0.01 ms, it is only the data very near the wellbore (lessthan 100m offset) that are marginally impacted. At the later time-channel, however, thedifference is significant. At 5 ms the difference between including permeability and notis over 40% for all offsets shown (panel m). Through time, at the 300 m offset, we cansee that the response decays more slowly, particularly in the time window between 2 msand 100 ms. The maximum difference is > 4000% at 20 ms; in the conductive well-case(d), the response has decayed to ∼ 10−8 V/m, while in the permeable well (h), the value174Figure 5.10: Simulated electric field at the surface of the earth, as was presentedin 5.4. The top row (a – d) are the data simulated with the conductivecasing. The second rows (e – h) are the data simulated with the conductive,permeable casing. The third row (i – l) is the difference (permeable andconductive minus conductive only), and the fourth row (m and n) show thedifference as a percentage of the conductive solution.of the electric field is nearly two orders of magnitude larger.Although one might typically associate magnetic permeability with a large impact onthe magnetic fields and fluxes, the difference in the dbθ/dt data is much less significantand localized to short-offsets, as can be seen in Figure 5.11.175Figure 5.11: Simulated db/dt at the surface of the earth, as was presented in 5.6.The top row (a – d) are the data simulated with the conductive casing. Thesecond rows (e – h) are the data simulated with the conductive, permeablecasing. The third row (i – l) is the difference (permeable and conductiveminus conductive only), and the fourth row (m and n) show the differenceas a percentage of the conductive solution.5.3.2 Down-hole sourceThe impact of the magnetic permeability of the casing is much more significant if weconsider an experiment in which one electrode is positioned within the borehole. In thiscase, the wire path is through the borehole, and the vertical current through it generatesa azimuthal magnetic field that is essentially perfectly coupled to the permeable casing.176I demonstrate this with an experiment where the positive electrode is attached to the per-meable casing at 950m depth. The currents through time are shown in Figures 5.12 and5.13 which show the total currents and the difference between an experiment with thepermeable pipe and one with a conductive pipe whose permeability is µ0, respectively.Shortly after shut-off, currents along the entire length of the pipe contribute to thedifference in the response between the permeable pipe and a conductive pipe. This isreflected in the radial electric field data shown in Figure 5.14In the radial electric field data, shown in Figure 5.14m, this corresponds to a differ-ence of over 60% at 1000 m offset in the 5 ms data.5.3.3 DiscussionMagnetic permeability introduces yet another complication into the physics of EM insettings with steel-cased wells. In the top-casing experiment, we saw that by introduc-ing permeability, a circulation of current in which current travelled downwards nearthe inner wall of the well-casing and upwards along the outer wall of the well casing,was generated. Such an effect is not observed if the casing is only conductive. Thisindicates that it is important to model both permeability and conductivity in numericalsimulations to capture the physics; this is not commonly done in geophysical electro-magnetic simulations. Furthermore, even low-frequency or late-time data, which aresometimes treated as DC, may be impacted by magnetic permeability. Typically onemight associate magnetic permeability with most significantly impacting magnetic fielddata, however these results show that the electric field in a top-casing experiment can beoff by >40% at late times, even 1000 m from the well if permeability is not consideredin the simulation. On the positive side however, one benefit that the permeability of thesteel may have in time-domain EM imaging experiments is that it increases the current177Figure 5.12: Current density for a down-hole time domain experiment with a con-ductive, permeable steel-cased well (5× 106 S/m, 100µ0). The positiveelectrode is connected to the casing at 950 m depth and the return electrodeis on the surface at x= 1000 m.178Figure 5.13: Difference in current density for a down-hole time domain experi-ment with a conductive, permeable steel-cased well (as in Figure 5.12) anda similar experiment with a conductive well.179Figure 5.14: Simulated electric field at the surface of the earth, as was presentedin 5.10 for a TDEM experiment in which the positive electrode is coupledto the casing at a depth of 950m.density within the formation at later times. As a result, there is a long time-window overwhich a target may be excited.5.4 Approximations to a steel cased wellThe above simulations accurately discretize the casing, with 4 cells across the thicknessof the casing wall. This is suitable for simulations in which the aim is to examine thephysics. When considering an inversion, which requires many forward simulations andin which it is often preferable to use a cartesian grid to capture the geologic structures180of interest, it may not be feasible to use such a highly-refined mesh to capture the finerdetails of the physics.In the previous chapter, I showed that if we treat the well as a solid cylinder with aconductivity that preserves the product of the conductivity and cross-sectional area ofthe casing, we achieved an identical solution (within numerical error) at DC for the totalcurrent and charge-per-unit-length along the casing (see Section 4.4). This raises thequestion about whether, or not, this same approximation can be used for an electromag-netic experiment. Also, is there a straight-forward approximation that can be used toestimate the magnetic permeability for a coarse-scale representation of the well? In thissection, I examine these two questions, first by considering a conductive well, and thenintroducing permeability.5.4.1 Electrical conductivityWe consider the same model of a 1 km long, conductive well (5× 106 S/m) in a 10−2S/m half-space as used in Section 5.2 and perform a similar TDEM experiment in whichone electrode is connected to the top of the casing and a return electrode is positioned 1km away. The approximate model I use is a solid cylinder which has a conductivity of1.8×106 S/m, so that both the true and approximate models have equal cross-sectionalconductances. Both simulations are performed on identical meshes; this minimizes nu-merical differences that might otherwise occur from employing different discretizations.To compare the responses, I first consider the total current along the length of thewell. If the approximation is valid, we should expect that the total current density inthe hollow-cased well is equal to the total current density in the cylinder approximatingthe well through time. Figure 5.15 shows (a) the current at two times (0.01 ms and 5ms) along the length of the well, and (b) current at 500 m depth through time for the181true solution. Similarly, (c) and (d) show the current for the approximation to the well,(e) and (f) show the difference (approximation minus true), and (g) and (h) show thatdifference as a percentage of the true solution. At early times, there is good agreementbetween the two solutions; at 0.01 ms, the percent-error in (g) is negligible and beyond200 m depth in (e), the error is 5 orders of magnitude smaller than the true current.Above 200 m depth, we do see some coherent error that reaches a maximum of ∼ 10−3A near the surface. However, this is still several orders of magnitude smaller than thetrue current at 0.01 ms. At later times the approximate solution underestimates the truecurrent. At 5 ms, the current is underestimated by 8% along the entire length of thewell, as can be seen in (g). Through time, the error shown in (h) reaches a maximum ofnearly 10% at 8 ms. There is a large increase in percent difference beyond 30 ms, this issimply because we are dividing the difference by the very small current at the late times.In electric field data measured at the surface, similar errors are observed. Figure 5.16shows the electric field data for the true model (top row), approximate model (secondrow), difference between the approximate and true models (third row) and difference asa percentage (fourth row). At early times (0.01 ms), there is minimal difference, but at 5ms, there is an 8% difference between the data generated with the true model and thosegenerated with the approximate model. At 300 m offset, the time-period between 1 and10 ms is when we observe the largest error (panel n).Errors of 8-10% are not as egregious as neglecting magnetic permeability, but theanomalies that we will be looking for when considering subsurface injections are smalland may only comprise 10% to 20% of the background signal, so an error of 10% iscertainly not negligible.In order to provide more context for where this approximation is breaking down, Ihave plotted the current density, charge density and electric fields through time in Figure182Figure 5.15: (a) Downward-directed current (a) within the casing at 0.01 ms (blue)and 5 ms (green) and (b) as a function of time at 500 m depth for the true,hollow-cased well. (c), (d) Current in the approximate model which treatsthe casing as a solid cylinder. The conductivity of the cylinder approxi-mating the well is chosen so that both models have equal products of theconductivity and the cross sectional area. (e), (f) Difference between thecurrent in the approximate model and the true model (approximate minustrue). (g), (h) Difference as a percentage of the true solution.183Figure 5.16: Simulated electric field at the surface of the earth, as was presentedin 5.4. The top row (a – d) are the data simulated with the true, conduc-tive casing. The second rows (e – h) are the data simulated with the solidcylinder approximating the well. The third row (i – l) is the difference (ap-proximate minus true), and the fourth row (m and n) show the difference asa percentage of the true solution.5.17. Panel (a) shows the true current density, (b) shows the secondary current density(the difference between the simulation using the approximation of the casing and thatusing the true casing), (c) shows the true charge density (d) shows the secondary chargedensity, (e) shows the true electric field, and (f) shows the secondary electric field. Notethat (a) and (b) are on the same logarithmic color scales as are (e) and (f). Panels (c)and (d) are each have their own associated linear colormap. Panels (b), (c) and (e) are184effectively looking at “what physics are we introducing by making this approximation?”At the DC limit, the solutions are in agreement, with the exception that the currentdistribution between the true casing and the approximate casing are different inside ofthe well, as expected. There is also a small secondary electric field at the bottom ofthe pipe – this is because by treating the pipe as a solid cylinder, we have introduceda conductivity contrast, and thus small secondary charge at that bottom interface. Thishowever has minimal impact on the currents in the formation or data measured at thesurface.Shortly after shut-off, however, the approximation begins to break down. There is adipolar charge introduced on the casing with a concentrated negative charge density nearthe top and a positive charge density beneath that. Corresponding to the positive charges,there are similarly outward-directed currents and electric fields. As time progresses, thesecondary dipolar charge and associated currents and electric fields spread out along thelength of the well.The passing image current changes the geometry of the currents and distributionof charges. At DC, the well has a positive charge along its length and all currents areoutward from the well. The incoming image current is a source of the negative charges;the currents enter the conductive casing from the resistive background, thus negativecharges build up at the interface. This indicates that the conductivity contrast betweenthe background and the casing is important when considering EM effects. At DC, theapproximation of the well by a cylinder with equal product of the conductivity and thecross sectional area is correct. Any alteration of the conductivity would introduce errorsin the initial condition. However, the approximation does not hold through time as thegeometry of the current-systems changes. This suggests that perhaps a conductivitywhich varies with time, and possibly along the length of the casing, may be necessary185Figure 5.17: (a) True current density (b) difference between the currents in thesimulation with the approximation to the conductive well and the true so-lution (secondary currents), (c) True charge density, (d) secondary chargedensity, (e) true electric field, (f) secondary electric field.186in order to accurately approximate the casing on a coarse scale.5.4.2 Magnetic permeabilityClearly, if we do not have a suitable approach for estimating a coarse-scale conductivityof a well in an EM experiment, we cannot expect to build up a strategy for a conduc-tive, permeable pipe. If the currents are incorrect, then the magnetic fields and fluxesare incorrect, and the interplay between the currents and magnetic permeability that Idemonstrated in Section 5.3 (Figure 5.12, in particular) will not be captured. Therefore,the purpose of this section is to introduce some of the additional challenges that areraised when permeability is considered and see if, by examining the physics, there areany insights to be gleaned about this problem.Simply based on the geometry of the currents and magnetic fields, one reasonable ap-proach to estimating the permeability of a solid cylinder which approximates the hollow-cased well is to preserve the product of the permeability and the casing thickness. Themagnetic fields and fluxes are rotational within the casing, so the thickness of the casingis a primary control on the magnetic flux within the casing.To introduce a first-order approximation for a conductive, permeable well, I preservethe product of the conductivity and cross-sectional area, as in the previous section, andpreserve the product of the permeability (100 µ0 for the true well) and the thickness ofthe casing (1 cm), giving a permeability of 20µ0. The experiment is again a top-casingsimulation, as in Section 5.3. The resultant current within the pipe and data from the trueand approximate simulation are shown in Figures 5.19 and 5.20, respectively (note thatthe y-limits on this plot are different than those shown in the conductive-well example,Figures 5.15 5.16). As might be expected, the discrepancy, in terms of percentage ofthe true solution, at late times is nearly twice what we observed in the conductive well187Figure 5.18: Sketch of the geometry of the dominant geometry of the current den-sity, ~j and magnetic flux density,~b, within the casing. Most of the currentflows vertically while most of the magnetic flux density is rotational.approximation. The direction of the difference is opposite to the conductive case. Here,the currents within the casing are overestimated, creating an outward-directed secondaryelectric field which reduces the magnitude of the radial electric field data relative to thetrue solution.To see where the solution is breaking down, I plot the currents, charges, and electricfields in Figure 5.21. Note that the colorbar for the secondary charges in panel (d) has arange that is an order of magnitude larger than that shown in Figure 5.17. The behaviorof the secondary currents charges and electric fields due to the approximation is oppositeto what we observed in Figure 5.17. Here, there is a positive secondary charge at thetop of the casing and a negative secondary charge beneath it, and as the image currentpasses, we evolve to a state where there are outward-directed secondary currents andelectric fields in the top portion of the casing and inward-directed secondary currentsand electric fields in the bottom portion of the casing.The interplay between currents, magnetic flux and magnetic permeability is clearlycomplex. Before developing a strategy for upscaling conductive, permeable wells, there188Figure 5.19: (a) Downward-directed current (a) within the casing at 0.01 ms (blue)and 5 ms (green) and (b) as a function of time at 500 m depth for the true,permeable, hollow-cased well. (c), (d) Current in the approximate modelwhich treats the casing as a solid cylinder. The conductivity of the cylinderapproximating the well is chosen so that both models have equal productsof the conductivity and the cross sectional area. The permeability of thecylinder is chosen so that both models have equal products of the perme-ability and the thickness. (e), (f) Difference between the current in theapproximate model and the true model (approximate minus true). (g), (h)Difference as a percentage of the true solution.189Figure 5.20: Simulated electric field at the surface of the earth, as was presentedin 5.4. The top row (a – d) are the data simulated with the true, conductive,permeable casing. The second rows (e – h) are the data simulated with thesolid cylinder approximating the well. The third row (i – l) is the difference(approximate minus true), and the fourth row (m and n) show the differenceas a percentage of the true solution.are a host of questions that will need to be addressed: Can a conductivity which variesvertically along its length and through time be used to approximate a casing? Does thisapproximation still hold if the conductivity of the background, and thus distribution ofthe currents within the formation is variable? If an appropriate description of the con-ductivity is found, is there a simple description of the permeability that can be adopted(e.g. is capturing the complexity of the currents the main item of concern?)?190Figure 5.21: (a) True current density, (b) difference between the currents in thesimulation with the approximation to the conductive, permeable well andthe true solution (secondary currents), (c) true charge density, (d) secondarycharge density, (e) true electric field, (f) secondary electric field.1915.5 ConclusionsThis chapter poses more questions than answers. Even for the (seemingly) simple modelof a conductive casing in a half-space, there are subtleties to the behavior of the currentsin an EM experiment. The interaction of the diffusing formation currents (set up at DC,before shut-off), the image current generated when the source current is shut-off, andthe highly conductive well produce a “shadow-zone” along a line 180◦ from the wire.The geometric complexities are slightly simplified if we subtract-off the background re-sponse of a grounded source over a half-space; the resulting “casing-response” is cylin-drically symmetric through time.Steel has a significant permeability, and the simulations in this chapter demonstratethat including permeability significantly alters the currents, particularly at later times.Within the well, a vertical circulation of current develops; this cannot be explained byconductivity. Within the formation, the currents do not diffuse as rapidly as they didwhen the well was only conductive. By 10 ms, in a 10−2 S/m half-space with a conduc-tive well, there was very little distortion of the currents due to the presence of the well.However, in the simulation with a conductive, permeable well – the currents within theformation were still significantly altered at 10 ms. Correspondingly, the late-time dataare much more significantly affected by permeability than the early-time data (upwardsof 40% difference at 5 ms for a top-casing experiment and 60% for a down-hole ex-periment). Many of the publications considering EM in settings with steel cased wellsneglect permeability. These results show that it is essential to consider. Furthermore,late-time or low frequency data that are simulated and inverted as DC data could becontaminated with permeability-effects.Using sufficient spatial discretization to accurately simulate data on cartesian meshesis costly, and can render the inversions, which require many forward simulations, infea-192sible. As a result, it is advantageous to reduce the size of the mesh and approximatethe casing. The simple rule-of-thumb of preserving the product of the conductivity andthe cross-sectional area of the casing was suitable for DC resistivity problems, but thisbreaks down when we consider EM problems. The incoming image current changesthe geometry of the fields and their interaction with the pipe through time. An effectivereplacement conductivity for a simulation on a coarse mesh may need to vary with timeand with distance along the pipe. Matters are more complicated when permeability is in-troduced, as there is an interplay between the currents, magnetic field, and magnetic fluxthrough time. Using a geometrically-motivated approximation for the permeability andpreserving the product of the permeability and the thickness of the casing introduce er-rors in the opposite direction to those observed when only conductivity was considered.I showed that secondary positive charges were introduced due to the approximation ofa conductive well, but negative charges were introduced as a result of the approxima-tion to the conductive, permeable well. How to approximate the permeability is unclearat this point, and will first require that a suitable strategy be developed for electricalconductivity.There is much yet to be investigated and understood about electromagnetics in set-tings with conductive, permeable, steel cased wells. Central to this work is the ability tosimulate, visualize, and interact with aspects of the physics. This chapter provides someglimpses into the behavior of the currents, fields and fluxes, and includes source codeand Jupyter notebooks (see Appendix A) that can be used to continue the exploration.193Chapter 6An inversion approach for subsurfaceinjections6.1 IntroductionThe previous three chapters have focussed on understanding the physics of electromag-netics over conductive, permeable, steel-cased wells. In this chapter, we return to thegoal of imaging subsurface injections through geophysical inversion. We view this as atime-lapse problem, in which one data set is obtained prior to the injection and seconddata set is obtained after the injection. The aim of the inversion, then, is to characterizethe changes in the earth model due to the injection.Time-lapse direct current (DC) resistivity, and in some cases electromagnetics (EM),is commonly used in the groundwater hydrology community. In particular, DC resistiv-ity has been used, in salt-tracer experiments aimed at understanding groundwater flow(e.g. Slater et al. (2002); Kemna et al. (2002); Singha and Gorelick (2005); Doetschet al. (2012)) and for characterizing time-lapse vadose zone processes (e.g. Daily et al.194(1992); Park (1998); Binley et al. (2002)). Within the context of hydrogeophysics, typ-ically multiple time-lapse surveys are conducted. A range of interpretation techniqueshave been applied. In terms of algorithmic complexity, the most straightforward ap-proach is to invert each time-snapshot independently. The recovered models can thenbe differenced and the difference interpreted (Cassiani et al., 2006), or an intermedi-ate interpretation, for example estimating the center of mass of a plume (Singha andGorelick, 2006; Doetsch et al., 2012), can serve as an indicator that is tracked throughtime. Daily et al. (1992) presented an approach for inverting ratios between initial andsubsequent datasets, and (LaBrecque and Yang, 2000) proposed an inversion for theresistivity difference between two subsequent data sets by inverting the difference be-tween the two data sets. A common approach is to use the inversion result from aninitial timestep as a reference and starting model for the following datasets (Loke, 2001;Oldenborger et al., 2007), this is sometimes referred to as a “cascaded inversion” (Milleret al., 2008). More advanced techniques simultaneously invert all of the time-snapshotsand apply both spatial and temporal regularization (Kim et al., 2009; Loke et al., 2014).Hayley et al. (2011) provides an overview and comparison of common time-lapse inver-sion approaches and demonstrates a comparison of several approaches for a syntheticexample of the remediation of a saline plume.In the context of reservoir imaging applications, the “time-lapse” aspect of the prob-lem is somewhat simpler than that typically encountered in hydrogeophysics. Very fewtemporally dense DC or EM data sets have been collected for reservoir imaging. No-table exceptions include the cross-well DC resistivity surveys described in Carrigan et al.(2013) and Tøndel et al. (2014), for monitoring of a deep CO2 injection in a gas fieldand monitoring steam chamber growth in a steam-assisted gravity drainage operationin the Athabasca Oil sands, respectively. Much more common are time-lapse data sets195consisting of only two times: pre-injection and post-injection. Cross-well EM has beenapplied for such surveys to monitor water-floods (Wilt et al., 2005, 2012) and steaminjections (Wilt et al., 1996, 1997; Marion et al., 2011). By far the most common in-version approach involves inverting for a background model with the initial data set,typically including constraints from well-logs and potentially interfaces from seismicdata and using this as the starting model for the inversion, typically in 2D (e.g. Wiltet al. (2012)).Although the “time-lapse” aspect of injection-monitoring applications is generallysimpler than applications in hydrogeophysics, the large depths considered in reservoirapplications means that we are usually working with small signals. Furthermore, thepresence of steel cased wells complicates signals. Efforts to reduce the impact of steelcased wells on DC surveys include using a coating on the casing which the electrodesare connected to and thus insulate the casing from the survey (Tøndel et al., 2014).For cross-well EM applications, fibreglass casings may sometimes be used (Wilt et al.,2012), or when a single well is cased, a “casing-correction” is applied to the data col-lected at a fixed frequency (Augustin et al., 1989b; Becker et al., 1997).For the emerging application of grounded-source DC or EM in the monitoring ofsubsurface injections, very few inversions of synthetic or field studies have been pub-lished, leaving many open questions. How does the casing affect our sensitivities in theinversion? In particular, Rucker (2012), have shown that in near surface studies wherethe casings are used as long electrodes, depth information is lost in the inversion. Ad-ditionally, forward modelling shows that the currents spread out along the length of thewell. These prompt the question: can we expect to resolve the location of the targetalong the well? Specific to subsurface injections, there is also further a-priori infor-mation that may be included, for example the conductivity of the injected material and196the volume of material injected should be known. What is the best way to include thisinformation in the inversion?For DC data, Weiss et al. (2015) works with the scattered potential due to the fractureand suggests an approach which approximates the well as a line-charge distribution andthe fracture as a point charge. The information content of a DC survey with only a fewsources is quite limited, so reducing the number of free parameters, as in Weiss et al.(2015), is appropriate for obtaining meaningful information from the data. Here, I wishto work towards the use of electromagnetic data. In addition to providing informationat the electrostatic limit provided by DC data, EM data include the role of inductiveprocesses and are therefore richer in their information content. As such, I am interestedin using an inversion approach which enables us to interpret geometric and physicalproperty information from the result.To focus discussion, we will consider synthetic examples related to hydraulic frac-turing. The goal of the inversion in this case is to delineate the extent and geometry of thepropped region of the fractured reservoir. Changes in electrical conductivity are viewedas a proxy for estimating the concentration of propped fractures. I will consider voxel-based and parametric inversions, as well as alternate parameterizations using effectivemedium theory. Using this range of approaches, I aim to develop an understanding ofthe non-uniqueness we face for grounded-source surveys in settings with steel casedwells. The analysis conducted in this chapter builds on the SimPEG framework and allof the examples shown are available in the form of Jupyter notebooks (see Appendix A).6.2 Choosing an inversion modelBackground material on the general inversion approach I follow was provided in theintroduction (Section 1.6) with further details included in Cockett et al. (2015) and Ap-197pendix D. The overview I provided laid out a general strategy for solving for the inver-sion model m, but I have not yet specified how m is chosen. This chapter will exploreseveral different approaches for defining an inversion model including parametric mod-els and using effective medium theory to invert for fracture concentration. This sectionprovides mathematical background on how these different model parameterizations aretreated in the inversion.The physical property which we aim to characterize in an electromagnetic inver-sion is electrical conductivity σ (or equivalently, its inverse, resistivity ρ). In a DC oran EM inversion, however, it is common to invert for log-conductivity on the forwardsimulation mesh, that isσ =M (m) (6.1)where M (m) = expm. We refer to M (·) as a mapping. Mappings have two implica-tions in the inversion. One implication is in the model regularization: we have changedthe space in which we are applying the regularization, for this example, we regularizeon log-conductivity values rather than linear conductivity. As the conductivity of com-mon earth materials varies over several orders of magnitude, it is preferable to penalizejumps in orders of magnitude between voxels rather than penalizing linear values. Thesecond implication is in the computation of the sensitivity. The forward simulation inan EM or DC problem depends upon electrical conductivity, thus the mapping modifiesthe sensitivity via the chain ruleJ[m] =dF [M (m)]dM (m)dM (m)dm=dF [σ(m)]dσdσdm(6.2)Mappings can be composed, for example if inactive cells are included in the mod-elling domain, such as air cells or cells capturing known structures such as a steel-cased198well, thenσ =M2(M1(m)) (6.3)where M1 injects in the log-conductivity values of the inactive cells and M2 takes theexponential. The sensitivity is appropriately modified by adding another step to thechain rule.In addition to mappings routinely employed in electromagnetic inversions, such asthose used for working with log-conductivity values and handling inactive cells in themodelling domain, there are two mapping which we will make use of in this chapter: aneffective medium theory mapping, based on the homogenization technique for propped,fractured reservoirs discussed in Chapter 2 and parametric mappings.Parametric mappingsThe other type of mapping I will make extensive use of in this chapter are paramet-ric maps. I consider model parameterizations of blocks and ellipsoids, similar to thatdescribed in McMillan et al. (2015a); McMillan (2017). For example, using a simpleparametrization of a block, the model ism = [mback,mblock,x0,∆x,y0,∆y,z0,∆z]> (6.4)where mback is the model value of the background, mblock, (x0,y0,z0) is the center of theblock and (∆x,∆y,∆z) are the widths of the block in each dimension. The mapping isthenM (m) = mback+(mblock−mback) s(τ(m)) (6.5)199Figure 6.1: Approximation to a step function using an arctangent, as described inequation 6.6. The values in the legend indicate the value of the slope, a.Smaller values have a more gradual transition between zero and one, whilelarger values have a more rapid transition.where s(·) is a differentiable approximation to a step-function and τ is a level set func-tion of the block. To approximate a step function, I use the arctangent function,s(τ) =1pitan−1(aτ)+12(6.6)where a controls the slope of the transition between 0 and 1. Small values of a result ina gradual transition while larger values give a sharper transition, as shown in Figure 6.1.For more robust performance of the Gauss-Newton inversion, I choose a such that thetransition happens over a multiple cells in the simulation mesh (McMillan, 2017).To define a level set function for a block, I need to define a function, τ , whichevaluates to τ < 0 when outside the target and τ > 0 inside the target. To define a block,I could, for example useτ = 1−(∣∣∣∣∣∣∣∣x− x0∆x/2∣∣∣∣∣∣∣∣2∞+∣∣∣∣∣∣∣∣y− y0∆y/2∣∣∣∣∣∣∣∣2∞+∣∣∣∣∣∣∣∣z− z0∆z/2∣∣∣∣∣∣∣∣2∞)(6.7)200however, the infinity norm is not differentiable. Thus, I approximate the infinity normusing an Ekblom norm (Ekblom, 1973),τ = 1−[(x− x0∆x/2)2+ ε2]p/2+[(y− y0∆y/2)2+ ε2]p/2+[(z− z0∆z/2)2+ ε2]p/2(6.8)where ε is a small constant and p is a constant describing the approximate norm, forexample, if p = 2, then equation 6.8 describes an ellipsoid. To represent a block, Ichoose a p that is sufficiently large. For the length scales I consider, a value of p = 4is appropriate. The value of ε is chosen to be large enough, and the value of p smallenough so that the derivatives of the mapping are stable and second-order for the lengthscales of the problem. In addition, rotations can be included in the model, as describedin McMillan (2017) 1.When employing parametric mappings, there are two important implications to notein the setup of the inversion. Since the mathematical statement of the inverse problemis, in principle, overdetermined (there are more data than model parameters), I fix thevalue of β at zero and do not employ a regularization. The second point is that for astarting model, it is important to start with the block and background having differentconductivities. This was similarly discussed in McMillan (2017).Effective medium theory mappingIn Chapter 2, I introduced a two-step process for estimating the electrical conductivityof a propped, fractured volume of rock. The first step involved estimating the effectiveconductivity of a mixture of electrically conductive proppant and fluid, and in the secondstep, I estimated the effective conductivity of a volume of rock which has fractures filled1McMillan (2017) also employs a weighting scheme to scale the model parameters. I have not foundthis to be necessary for the examples I consider and thus do not use any weights or scaling201with the proppant-fluid composite.Assuming the electrical conductivity of the fluid and proppant are known, then ratherthan inverting for electrical conductivity, we can invert for the concentration of conduc-tive fractures. To avoid introducing additional non-uniqueness into the problem, I usea fixed ratio of proppant and fluid within the propped region of the reservoir and treatthe fracture concentration ϕ as the inversion model (e.g. m = ϕ). In a voxel-based in-version, ϕ is a vector with a value for the concentration in each cell. For simplicity, Iassume that the fractures are randomly oriented and work only with isotropic conduc-tivities. In this case, the mapping requires that we solve the two-pahse effective mediumtheory approximation,(1−ϕ)(σ∗−σ0)R(0,∗)+ϕ(σ∗−σ1)R(1,∗) = 0 (6.9)for the effective conductivity, σ∗. The background has conductivity σ0 and the conduc-tive, proppant-filled cracks have conductivity σ1. Note that σ0, the conductivity of thebackground, does not need to be a scalar, it can be a vector with a background conduc-tivity value for each voxel in the mesh. These values can be obtained by first invertingthe pre-fracture data. The electric field concentration tensor R(i,∗) captures the geometryof the particles that compose each phase. Note that for randomly oriented fractures R(i,∗)is a scalar (1/3traceR(i,∗), where R(i,∗) is given in equation 2.3:R( j,∗) =13trace([I+Aσ∗−1(σ j−σ∗)]−1)(6.10)For the background, I use an aspect ratio of 1, assuming a spherical geometry for theparticles that compose it, and for the fractures, I use a small aspect ratio (∼ 10−4−10−3)and treat them as ellipsoidal cracks. In Section 2.2.3, I demonstrated that for sufficiently202thin fractures, the exact aspect ratio is not significant. For a gradient based inversion, wealso require the derivative of the effective conductivity, σ∗, with respect to our model,which is the fracture concentration, e.g.∂σ∂m=∂σ∗∂ϕ(6.11)This derivation is outlined in Appendix C.The general time-lapse inversion workflow using effective medium theory is twosteps:• Construct a background conductivity model, σ0, by inverting pre-fracture data.This inversion is a voxel-inversion for log-conductivity.• Invert the post-fracture data for fracture concentration, φ . Note that using a start-ing model and reference model of m0 = 0 is equivalent to starting the post-fractureinversion with the pre-fracture model.Self-consistent effective medium theory is the method I adopt to connect the concen-tration of fractures with the effective conductivity of a fractured volume of rock, how-ever, other relationships, such as an empirical relationship estimated from a lab study,could equally be employed. One interesting implication of relating the concentration ofthe fractures to the change in conductivity is that this provides a conduit for bringingin a-priori information about the volume of proppant injected into the reservoir. Thepredicted volume isVpred =∫ϕ dV (6.12)203I can then define a volume data misfit term,φV =12∣∣∣∣∣∣∣∣ 1εV (Vpre