Forward Modelling and Inversion of Time-DomainElectromagnetic Geophysical Surveys in the Presence ofChargeable MaterialsbyPatrick BelliveauB.Sc., Simon Fraser University, 2010M.Sc., Memorial University of Newfoundland, 2014A 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)April 2019c© Patrick Belliveau, 2019The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled:Forward Modelling and Inversion of Time-Domain Electromagnetic Geophys-ical Surveys in the Presence of Chargeable Materialssubmitted by Patrick Belliveau in partial fulfillment of the requirements for thedegree of Doctor of Philosophy in GeophysicsExamining Committee:Eldad Haber, Earth, Ocean and Atmospheric SciencesSupervisorDouglas Oldenburg, Earth, Ocean and Atmospheric SciencesSupervisory Committee MemberRoger Beckie, Earth, Ocean and Atmospheric SciencesUniversity ExaminerChen Greif, Computer ScienceUniversity ExaminerAdditional Supervisory Committee Members:Christian Schoof, Earth, Ocean, and Atmospheric SciencesSupervisory Committee MemberiiAbstractThe work described in this thesis examines transient geophysical electromagneticforward modelling and inversion in the presence of induced polarization (IP) ef-fects. The thesis introduces a new method of modelling IP using stretched expo-nential relaxation. A three-dimensional (3D) forward modelling algorithm takingfull account of the coupling of IP and electromagnetic induction is developed. Thestretched exponential modelling algorithm has been implemented using efficientnumerical methods that allow it to tackle large-scale problems and are amenableto use in inversion. In particular, a parallel time-stepping technique has been de-veloped that allows transient electric fields at multiple time steps to be computedsimultaneously. The behavior of the stretched exponential model is demonstratedby applying it to synthetic numerical examples that simulate a grounded sourceIP survey with significant electromagnetic induction effects and a concentric-loopairborne electromagnetic sounding over a polarizable body.An inversion algorithm using the stretched exponential model was developedthat is able to recover the 3D structure of physical properties of the earth relatedto IP from transient geophysical electromagnetic data. The method is tested on asimple synthetic example problem. The thesis finishes with the development ofa novel stochastic parametric level-set inversion algorithm, which could be usefulin applying stretched exponential inversion to real world problems in the future.The algorithm addresses some of the shortcomings of the simple inversion ap-proach used for stretched exponential inversion earlier in the thesis. The stochasticparametric inversion algorithm is used to solve shape reconstruction inverse prob-lems in which the object of interest is embedded in a heterogeneous backgroundmedium that is known only approximately. Shape reconstruction is posed as aiiistochastic programming problem, in which the background medium is treated asa random field with a known probability distribution. It is demonstrated that byusing accelerated stochastic gradient descent the method can be applied to large-scale problems. The capabilities of the method are demonstrated on a simple 2Dmodel problem and in a more demanding application to a 3D inverse conductivityproblem in geophysical imaging.ivLay SummaryThis thesis describes research in the field of electromagnetic geophysical imaging.The work led to the development of novel algorithms and software to facilitatethree-dimensional imaging of the electromagnetic properties of the earth’s subsur-face. Electromagnetic imaging results are useful to geologists and other geoscienceprofessionals in their attempts to assess subsurface geology for applications such asnatural resource exploration and environmental assessment. The particular focus ofthis thesis was on improving the quality and efficiency of imaging algorithms thatconsider an electromagnetic phenomenon known as induced polarization. Inducedpolarization effects can be important indicators of geological features of interestto practitioners but they are difficult to model and standard methods either ignorethem or treat them in some simplified form. This thesis has developed the compu-tational methods necessary to fully consider the effects of induced polarization inelectromagnetic imaging.vPrefaceThis thesis is my original work. It has resulted in one peer reviewed journal publi-cation, one peer reviewed journal submission currently under review, one expandedabstract publication, and seven conference presentations.The parallel time-stepping algorithm in chapter 2 originated in discussions withDr. Haber and is based on his original idea. I developed the algorithm and imple-mented it. It was published in [11]. Chapter 2 provides more detailed analysis andmore extensive background material than was included in the published paper.The stretched exponential forward modelling algorithm in chapter 3 also origi-nated in discussions with Dr. Haber and is based on his original idea. I developed,implemented and tested the algorithm based on the discussions with Dr. Haber. Aversion of chapter 3 was published as [11]. Chapter 3 contains a more extensivecomparison of the stretched exponential model of induced polarization with theCole-Cole model than was included in the published paper.The stretched exponential inversion algorithm and workflow described in chap-ter 4 is my original work. It was carried out independently but benefited fromdiscussions with Dr. Haber and Roman Shektman.The stochastic inversion algorithm described in chapter 5 originated in discus-sions with Dr. Haber based on my original idea. I implemented and analyzed thealgorithms. A version of this chapter has been submitted for publication and iscurrently under review.The work in this thesis entailed significant software development. All devel-opment was done within the open source jInv software framework, which wascreated in Dr. Haber’s research group at UBC and is available online at https://github.com/juliaInv. Four of the conference presentations I gave over the courseviof my Ph.D. focused on the jInv software package.viiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Geophysical background . . . . . . . . . . . . . . . . . . . . . . 31.2.1 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . 51.2.2 IP models . . . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Efficient modelling of non-dispersive TEM data . . . . . . . . . . . 122.1 Discretization of the time-domain quasi-static Maxwell equations 13viii2.1.1 Spatial discretization . . . . . . . . . . . . . . . . . . . . 142.1.2 Discretization in time . . . . . . . . . . . . . . . . . . . . 172.2 Parallel time-stepping . . . . . . . . . . . . . . . . . . . . . . . . 262.2.1 Evaluating accuracy . . . . . . . . . . . . . . . . . . . . 282.2.2 Evaluating performance . . . . . . . . . . . . . . . . . . 302.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353 Forward modelling of the dispersive time-domain Maxwell Equa-tions using the stretched exponential function . . . . . . . . . . . . . 363.1 The stretched exponential formulation . . . . . . . . . . . . . . . 363.1.1 Comparison of stretched exponential and Cole-Cole models 413.2 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.2.1 Backward Euler Algorithm . . . . . . . . . . . . . . . . . 453.2.2 Solving the linear system . . . . . . . . . . . . . . . . . . 473.2.3 BDF2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.2.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . 533.3 Synthetic modelling examples . . . . . . . . . . . . . . . . . . . 543.3.1 Grounded source example . . . . . . . . . . . . . . . . . 543.3.2 Inductive source example . . . . . . . . . . . . . . . . . . 583.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614 Stretched exponential inversion . . . . . . . . . . . . . . . . . . . . . 634.1 The gradient array example . . . . . . . . . . . . . . . . . . . . . 644.2 Inversion methodology . . . . . . . . . . . . . . . . . . . . . . . 654.3 Inversion results . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735 Randomized background stochastic inversion . . . . . . . . . . . . . 755.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 755.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 775.2.1 Computing Background Samples . . . . . . . . . . . . . 815.2.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . 825.3 Application to a Simple Model Problem . . . . . . . . . . . . . . 835.3.1 SAA inversion . . . . . . . . . . . . . . . . . . . . . . . 86ix5.3.2 Optimization by stochastic gradient descent . . . . . . . . 915.4 Application to non-linear 3D problem . . . . . . . . . . . . . . . 935.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 996 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105xList of TablesTable 1.2 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . 6Table 2.1 Parallel time-stepping scheme. Six parallel forward modellingsimulations with constant step-sizes shown in this table wererun in parallel to simulate the receiver voltages at the same ob-servation times used in the reference backward Euler simula-tion. The number of observation times per parallel process var-ied but the total number of time-steps was the same in all thesimulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Table 2.2 Backward Euler solution run times. . . . . . . . . . . . . . . . 33Table 2.3 Parallel time-stepping run times using 6 parallel time-steppingprocesses run across 4 compute nodes of a computer cluster. . . 34Table 5.1 Starting guesses and parameter bounds for straight ray tomog-raphy inversion parameters. φ is the ellipse orientation angle in, x and y are the horizontal and vertical coordinates of the centreposition, and s the slowness inside the ellipse. . . . . . . . . . 87Table 5.2 Marginal statistics of SAA inversion results. The sample meansof the recovered values of each parameter are shown for the 1,10, and 100 background inversion experiments in a), with thecorresponding sample standard deviations shown in b). Recallthat the true parameter values are φ = 30◦, x= 0.5, y= 0.5, ands = 4. See table 5.1 for parameter bounds and initial guesses. . 89xiTable 5.3 Marginal statistics of SAA inversion results using the ADAMoptimization algorithm. Recall that the true parameter valuesare φ = 30◦, x = 0.5, y = 0.5, and s = 4. See table 5.1 forparameter bounds and initial guesses. . . . . . . . . . . . . . . 93xiiList of FiguresFigure 1.1 Illustration of typical grounded and inductive source surveys.Left image shows current electrodes connected to ground andenergy source, and receiver measuring the potential differencebetween separate measurement electrodes. Right image showsa concentric loop airborne electromagnetic (EM) system, wherea current run through the outer loop creates a time-varyingmagnetic field that induces a response from the earth that ismeasured by the inner coil. . . . . . . . . . . . . . . . . . . . 4Figure 1.2 DC experimental setup. From Tomoquest.com . . . . . . . . 5Figure 2.1 Yee grid. Components of the electric field are discretized onmesh edges in corresponding Cartesian coordinate directions.Magnetic fields are discretized on face centres. Figure adaptedfrom [51]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16Figure 2.2 Typical example of realistic time-domain electromagnetic (TDEM)transmitter current time-dependence plotted in black, with dotsshowing times at which solution is computed. The solid blueline shows the variation in the step-size over the full time rangeof the simulation. There are six unique step-sizes. . . . . . . . 20Figure 2.3 Transmitter current step-off waveform. Assumes transmittercurrent has been on long enough that fields are steady-state att = 0, where forward modelling begins. Shutoff takes place inone time-step. . . . . . . . . . . . . . . . . . . . . . . . . . . 22xiiiFigure 2.4 Comparison of time-stepping methods in concentric loop ex-ample. In the legend, BDF2 1 step init refers to the solutioncomputed using backward Euler for the first step only, usingsecond order backward differentiation formula (BDF2) there-after with no special treatment of the discontinuity. BDF2 3step init refers to the solution computed using backward Eu-ler time-stepping for the first three steps and BDF2 steppingthereafter. Dashed lines indicate negative values. a) shows theoverall behaviour of each method over the full 200 time-stepsimulation. c) illustrates the large oscillation observed fromthe second time-step when using BDF2 with one step initial-ization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25Figure 2.5 Parallel time-stepping accuracy on grounded source test problem. 29Figure 2.6 Time-stepping scheme for reference backward Euler solution. There are 10 steps each of sizes 2× 10−5 s, 4× 10−5 s, 8×10−5 s, 1.6×10−4 s, and 14 steps of size 3.2×10−4 s. . . . . 31Figure 2.7 Pardiso parallel scaling results when using a single sequentialbackward Euler time-stepping process. Parallel efficiency forN threads is T1/(N ·TN) where T1 is the sequential run time andTN is the run time using N threads. . . . . . . . . . . . . . . . 34Figure 3.1 a) Stretched exponential impulse response for various valuesof β . b) Cole-Cole impulse responses for various values of c.For both plots η = 0.1, τ = 0.1 and σ∞ = 0.1. . . . . . . . . . 42Figure 3.2 Cole-Cole impulse response with c = 0.5, τ = 0.1, η = 0.1,σ∞ = 0.1 plotted in blue alongside stretched exponential im-pulse responses with parameters chosen by non-linear leastsquares fitting. The green line attempted to fit the 0.001 <t < 100s time range and the magenta line the 0.1 < t < 100srange. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Figure 3.3 Error in PCG computed example problem electric fields forloose and tight PCG relative residual stopping tolerances. . . . 51xivFigure 3.4 Gradient array example survey layout. a) Cartoon showingtransmitter, receiver array bounding box, and sample cartoonreceiver. b) Plan view layout of transmitter, receiver array andtwo block synthetic model. . . . . . . . . . . . . . . . . . . . 55Figure 3.5 Visualization of OcTree mesh discretization of the conductiv-ity model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56Figure 3.6 DC inline electric field at the earth’s surface for two block ex-ample model. The locations of the buried blocks are shown byblack squares. . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 3.7 In-line component of surface electric field for gradient arrayexample. a) Plan-view in-line field at 0.01 s after current shut-off. Left and right black squares show the boundaries of thechargeable and conductive blocks, respectively. b) Plan viewfield at 0.14 s. c) Electric field decays above the centre ofthe chargeable block, with stretched exponential field in blueand corresponding non-dispersive field (η = 0) in red. Dashedlines indicate negative values. d) Electric fields above the cen-tre of the conductive block, showing that the stretched expo-nential and non-dispersive responses are almost identical. . . . 58Figure 3.8 Section view of block in a half-space with chargeable over-burden model. The overburden depth is 10 m and the blockextends from 40-90 m depth with a plan view cross-sectionof 100×100 m. The background halfspace conductivity was0.005 S/m and the conductivity of the block was 0.1 S/m. Theoverburden had σ∞ = 0.01 S/m η = 0.3, τ = 1× 10−3 s andβ = 0.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59xvFigure 3.9 Effect of chargeable overburden on ATEM response from buriedconductor. All plots show the vertical component of db/dt.The background halfspace response is shown in green in eachplot. Dashed lines indicate negative values. a) Block withno overburden. The response with the conductive block at adepth of 40 m is shown in blue. The magenta line shows theresponse obtained when the top of the block is moved to adepth of 70 m. b) Stretched exponential response with 10 mof overburden (σ∞ = 0.005, η = 0.3, τ = 1× 10−3, β = 0.8)is shown in black, along with the overburden free and back-ground responses. c) Stretched exponential response with thetop of the block at 70 m in black, along with the overburdenfree and background responses. d) Stretched exponential re-sponse with chargeable overburden and no conductive block inblack, along with background EM response. . . . . . . . . . . 60Figure 4.1 Gradient array example survey layout. a) Cartoon showingtransmitter, receiver array bounding box, and sample cartoonreceiver. b) Plan view layout of transmitter, receiver array andtwo block synthetic model. . . . . . . . . . . . . . . . . . . . 64Figure 4.2 Selected slices of true conductivity model. Left hand blockhad σ∞ = 0.018 S/m and right block σ∞ = 0.1 S/m. Left blockhad η = 0.25, β = 0.5 and τ = 0.5 s. Right block was notchargeable. a) z =−85 m b) y = 75 m c) y =−75 m. . . . . . 65Figure 4.3 Selected slices of recovered conductivity model from DC andearly off-time inversion. Thick black lines show true blockboundaries. True conductivities are 0.01 S/m for the back-ground, 0.018 S/m for the left block and 0.1 S/m for the rightblock. a) z =−85 m b) y = 75 m c) y =−75 m. . . . . . . . 68Figure 4.4 Misfit convergence curve for σ∞ inversion. . . . . . . . . . . 69xviFigure 4.5 Selected slices of recovered chargeability model from DC andmid to late off-time data inversion. Thick black lines showtrue block boundaries. Left hand block has η = 0.25. Righthand block is not chargeable. a) z = −85 m b) y = 75 m c)y =−75 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 4.6 Selected slices of recovered time constant model from late off-time data inversion. Thick black lines show true block bound-aries. Left hand block has true τ = 0.5. Right hand block isnot chargeable. a) z =−85 m b) y = 75 m c) y =−75 m. . . . 70Figure 4.7 Selected slices of recovered exponent model from late off-timedata inversion. Thick black lines show true block boundaries.Left hand block has true β = 0.5. Right hand block is notchargeable. a) z =−85 m b) y = 75 m c) y =−75 m. . . . . . 70Figure 4.8 Selected slices of recovered chargeability model from stage 4mid to late off-time data inversion. Thick black lines show trueblock boundaries. Left hand block has true η = 0.25. Righthand block is not chargeable. a) z = −85 m b) y = 75 m c)y =−75 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 4.9 Selected slices of recovered time constant model from stage4 mid to late off-time data inversion. Thick black lines showtrue block boundaries. Left hand block has true τ = 0.5. Righthand block is not chargeable. a) z = −85 m b) y = 75 m c)y =−75 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 4.10 Selected slices of recovered exponent model from stage 4 midto late off-time data inversion. Thick black lines show trueblock boundaries. Left hand block has true β = 0.5. Righthand block is not chargeable. a) z = −85 m b) y = 75 m c)y =−75 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 4.11 Voltage data for the x-directed receivers due to the x-directedtransmitter at 0.02 s after source shutoff. . . . . . . . . . . . . 73Figure 4.12 Voltage data for the x-directed receivers due to the x-directedtransmitter at 0.02 s after source shutoff. . . . . . . . . . . . . 73xviiFigure 4.13 true and predicted voltage data above the centre of the charge-able block from 0.02 s after source shutoff onward plotted onlog scale. Blue line shows observed data and red line showspredicted data. Dashed lines indicate negative values. . . . . . 74Figure 5.1 Cartoon schematic of straight ray tomography experiment show-ing straight ray paths from the transmitters to the receivers. . . 83Figure 5.2 True model used for straight ray tomography inversions. Themean model a) consists of homogeneous horizontal layers. Therandom perturbation shown in b) is added to the mean modelto form c), used as the background in the true model shown in d). 84Figure 5.3 Results of 100 independent single background inversions. Therecovered values for each inversion parameter are plotted ashistograms with accompanying kernel density estimates. Thered vertical lines show the true parameter values. Note thatthe estimated densities may have support outside the parame-ter bounds. They are included in order to give general visualimpressions of the spread of inversion results and these effectsare not of concern. . . . . . . . . . . . . . . . . . . . . . . . 88Figure 5.4 Kernel density plots showing inversion results for single back-ground, 10 background SAA, and 100 background SAA in-versions. 100 background angle results are not shown be-cause they took the upper bound value of 90◦ for all but twoof these inversions rendering the corresponding kernel densityplot meaningless. The vertical axes of the plots were clippedin order to show the single and 10 background results at rea-sonable scales. . . . . . . . . . . . . . . . . . . . . . . . . . 90Figure 5.5 Kernel density plots of inversion results for 100 backgroundand 500 background ADAM inversions with mini-batch size1. Horizontal axis limits are the same as in fig. 5.4 . . . . . . 94Figure 5.6 Mean background conductivity in the x-z plane at y = 0 in thecore region of the mesh. . . . . . . . . . . . . . . . . . . . . 97Figure 5.7 True model slices . . . . . . . . . . . . . . . . . . . . . . . . 98xviiiFigure 5.8 Results from single background and SAA ADAM inversionsof DC data. Blue dots show the recovered values of each pa-rameter from the 20 single background inversions and the starsshow the parameter values recovered by the SAA inversion. Ineach plot, the vertical axis limits are the parameter bounds, thedashed black line shows the initial guess and the red line thetrue value of the parameter. . . . . . . . . . . . . . . . . . . . 100xixGlossary2D two-dimensional3D three-dimensionalATEM airborne time-domain electromagneticBDF backward differentiation formulaBDF2 second order backward differentiation formulaDC direct currentDCIP direct current induced polarization, an IP experiment conducted alongsidea DC survey. EM induction effects are not considered.EM electromagneticEMIP coupled electromagnetic induction and induced polarizationFLCBDF2 fixed leading coefficient second order backward differentiation for-mulaFDTD finite difference time domainIP induced polarizationODE ordinary differential equationPDE partial differential equationRMS root mean squaredxxTDEM time-domain electromagneticxxiAcknowledgmentsFirst and foremost I would like to acknowledge the patient and thoughtful guid-ance of my supervisor, Eldad Haber. He helped shape my research direction whilealways encouraging me to follow my curiosity and believe in my abilities. I wasalso very lucky to have the chance to learn from Doug Oldenburg, who taught memuch about electromagnetics and about teaching.Seogi Kang and Dave Marchant were invaluable sources of sound advice on IPand EM as well as being warm, friendly, and collegial folks to chat with aroundthe office. Christoph Schwarzbach was also an invaluable source of insight onEM. These folks, along with fellow office mates and geophysics students SannaTyrva¨inen, Xiaojin Tan, Michelle Liu, Luz-Angelica Caudillo Mata, Wenke Wil-helms, Devin Cowan, Daniel Bild-Enkin, Mike MacMillan, Mike Mitchell, ThibautAstic, Dom Fournier, Lindsey Heagy, and Klara Steklova, made UBC a collegialand downright fun place to work.I would be remiss not to mention the two pre-UBC mentors in academia that in-troduced me to research, gave me much guidance and shared their love of learning:my undergraduate research supervisor Gwenn Flowers and my M.Sc. supervisorColin Farquharson.Most importantly, I could not have made it through my Ph.D. without the con-stant love and support of my parents, my sister Rachael, and my partner Caroline.xxiiDedicationTo my familyxxiiiChapter 1Introduction1.1 Research MotivationAt a broad level, this thesis is concerned with numerical modelling of transientelectromagnetic (EM) geophysical surveys and associated parameter estimationproblems. The particular emphasis is on understanding how to do this efficientlywhen the surveys are affected by the phenomenon of induced polarization (IP).Geophysical methods designed to detect IP have a long history in mineral explo-ration. They have been used in many exploration contexts but are particularlyimportant in the search for disseminated mineralization. Some important types ofdisseminated mineralization exhibit strong IP anomalies but are not sensitive toother geophysical methods—see e.g. Ward 96, Zonge et al. 106. IP is the mostimportant geophysical method in exploration for porphyry copper deposits, for ex-ample [106]. It is also increasingly important in environmental applications [54].IP is a dynamic EM phenomenon. However, in traditional processing and inter-pretation of geophysical surveys designed to detect it, it is often treated as a staticelectrical phenomenon—e.g. as in [67]. Such treatments neglect two key factorsthat affect IP data. They ignore the nature of the time or frequency dependenceof IP relaxations and the coupling between IP and EM induction. A full treatmentof the problem requires the ability to solve the quasi-static Maxwell equations—ineither the frequency or time-domain—using a constitutive law that models the EMresponse of polarizable materials.1Traditional IP surveys excite the earth using energy sources connected to theearth by grounded electrodes. IP effects can also be generated by closed loop sys-tems that excite the earth by EM induction. The increasingly common observationof IP effects in airborne transient inductive source EM data has generated renewedinterest in understanding and modelling IP in the context of transient electromag-netics. In this thesis I develop tools that are general and efficient for modellinggrounded and inductive source transient electromagnetic surveys in the presenceof polarizable media.Simulation tools that can accurately model realistic three-dimensional (3D) ge-ometries are very useful both for developing physical understanding and in design-ing surveys. Modelling the coupling of EM and IP effects can help a geophysicistunderstand in detail how EM coupling might affect a grounded source IP surveyand allow them to use inductive responses as useful signal rather than just noise.Some work has been done on modelling the quasi-static time-domain Maxwellequations in the presence of polarizable media (e.g. [25, 62]) but it remains anunder-researched area.The first two chapters of the thesis are focused on forward modelling and thelast two on geophysical inversion. Inversion, the process of estimating a subsurfacephysical property model from data, is a key step in geophysical data processing/in-terpretation. We want to be able to invert time-domain EM data for both conduc-tivity and IP parameters in 3D, whether the data comes from traditional groundedsource surveys designed to detect IP or from inductive source surveys. After de-veloping an algorithm for simulating IP effects in transient electromagnetic data,described in chapter 3 of this thesis, I describe the corresponding inversion algo-rithm in chapter 4.In addition to the need for appropriate forward modelling tools, a major diffi-culty in inverting EM data for IP parameters is underdeterminedness in the inverseproblem. Not having enough data to uniquely define a model of the subsurfaceis standard in geophysical inversion but the problem is exacerbated in the case ofmulti-parameter problems like IP inversion where it may be difficult to distinguishwhether some features in the data are due to IP or to variations in conductivity.Also, in multi-parameter inverse problems, a geophysicist is often only interestedin a subset of the parameters but is forced to estimate them all in order to create a2model that matches observed geophysical data. I address these issues in the finalchapter of this thesis by developing a stochastic inversion methodology that im-proves the precision of estimates of parameters of interest due to uncertainty in thenuisance parameter estimates.1.2 Geophysical backgroundThis section will place IP in the wider context of electromagnetic geophysics, de-scribe some general principles of EM surveying, and how IP fits into the pictureby way of constitutive relations. By electromagnetic geophysics, here we meanlow frequency (or long time) methods that are sensitive primarily to the electri-cal conductivity of the subsurface. This excludes high frequency techniques suchas ground penetrating radar which are, of course, electromagnetic in nature butfundamentally different from the range of low-frequency methods normally en-compassed by the term electromagnetic geophysics.These methods are part of the field of applied geophysics, whose major goal isthe mapping of variations in the material properties of the earth’s subsurface usingremote measurements. Such information can be useful in characterizing subsurfacegeology or locating buried objects without labour intensive and expensive directmethods such as drilling [92].The family of geophysical EM methods may be divided into natural sourcemethods, which sense the earth response to natural variations in atmospheric EMfields, and controlled source methods that measure the response due to artificialexcitations of the earth controlled by a geophysicist. This thesis will focus exclu-sively on controlled source methods. These may be further divided into groundedand inductive source methods. In grounded source surveys the energy source isdirectly connected to the earth by means of electrodes, allowing electric currentsto flow from the source into the earth. By contrast, inductive source surveys useclosed coils of wire to create time-varying magnetic fields that excite responsesin the earth through EM induction. Typical source and receiver arrangements forthese survey types are shown in fig. 1.1.These methods have been designed primarily to sense the electrical conductiv-ity of the subsurface but they can all be influenced by the IP phenomenon. When3VFigure 1.1: Illustration of typical grounded and inductive source surveys.Left image shows current electrodes connected to ground and energysource, and receiver measuring the potential difference between sepa-rate measurement electrodes. Right image shows a concentric loop air-borne EM system, where a current run through the outer loop creates atime-varying magnetic field that induces a response from the earth thatis measured by the inner coil.an electric field is applied to a conductive material, ions will move freely throughit. This is the flow of electric current. When the movement is impeded, chargeswill build up, whether that be at the interface between two homogeneous regions ofdiffering conductivities or due to more complex microscopic effects. This createselectrical polarization in the subsurface. Normally, when the applied field is re-moved, the charges will return to an electrically neutral state extremely quickly—far too quickly to be detected in an ordinary geophysical experiment. However,in some geological materials large polarization effects are observed that build upslowly under an applied current and then relax slowly enough to be measured ina geophysical EM survey. This is the IP phenomenon. The effect has complexorigins in the microscopic properties of the rocks. At the macroscopic level po-larization acts to impede current flow due to an externally applied field. Then,when the applied field is removed, the relaxation of the polarized charges to equi-librium will result in a small but measurable current flow. This effect is modelledas reduced steady-state conductivity in the static case and as frequency dependentconductivity for harmonic sources. In the time-domain the current flow at a giventime is modelled as depending on the history of the electric field driving the currentflow and not just its instantaneous value.In the canonical IP experiment, a direct current is driven into the ground using apair of electrodes connected to an energy source. The current remains switched on4Figure 1.2: DC experimental setup. From Tomoquest.comuntil the polarization has reached a steady state, at which point potential differencesbetween electrodes at various locations on the earth’s surface are measured. Thenthe current is abruptly switched off. The potential differences resulting from thecurrent created by the polarized charges returning to equilibrium is then measured.This is illustrated in fig. 1.21.2.1 Maxwell’s equationsTo describe this process quantitatively, we turn to Maxwell’s equations. Macro-scopic EM fields are governed by the macroscopic Maxwell equations, the funda-mental laws of classical electromagnetism. These may be stated as a set of partialdifferential equations (PDEs) describing the behaviour of EM fields, augmentedby a set of empirical rules called constitutive relations that define how EM fieldsare modified by materials. Maxwell’s equations are collected in table 1.2. Theequations may be expressed in the time-domain, where the fields are functions ofspace and time. What are known as the frequency-domain equations are obtainedby taking the Fourier transforms with respect to time of the time-domain equations.The macroscopic Maxwell equations relate several vector fields: the electricfield e, the magnetic induction b, the magnetic field h, the electric displacement dand the electric current j. There are three constitutive laws. The first relates the5Name Time FrequencyFaraday’s law ∇× e =−∂b∂ t(1.1.1) ∇×E =−iωB (1.1.2)Ampe`re-Maxwell ∇×h− j = ∂d∂ t(1.1.3) ∇×h− j = iωD (1.1.4)Gauss’s law ∇ ·d = ρ f (1.1.5) ∇ ·D = ρ f (1.1.6)Unnamed ∇ ·b = 0 (1.1.7) ∇ ·B = 0 (1.1.8)Table 1.2: Maxwell’s equationselectric field to the electric displacement, the second, the magnetic induction to themagnetic field, and the third the electric field to the electric current. In generalthese may be non-linear functionsd = d(e) (1.1)h = h(b) (1.2)j = j(e) (1.3)but the EM properties of most earth materials may be described by linear constitu-tive relations. These are typically writtend = εe (1.4)h = µ−1b (1.5)j = σe, (1.6)where ε is called the electrical permittivity, µ the magnetic permeability, and σ theelectrical conductivity. The constants of proportionality may be scalars or symmet-ric rank-2 tensors but depend only on position, not on time or frequency.In this thesis I will consider Maxwell’s equations in the quasi-static approxi-mation, in which the electric displacement term in Ampe`re-Maxwell equation is6ignored. This is a standard assumption in geophysical EM [97]. Most materi-als have ε on the order of the permittivity of free space, ε0 = 8.85× 10−12 F/m.The least conductive common earth materials have conductivities on the order of10-6 S/m. Larger values in the range 10-3-1 S/m are most often encountered andsome important geological materials such as nickel sulphides can have conductivi-ties as high as 105 S/m [70]. Frequency-domain geophysical experiments typicallyuse frequencies under 105 Hz and time-domain experiments operate on time-scaleslonger than 10-6 s. Consider the Ampe`re-Maxwell equation in frequency usinglinear constitutive relations∇×µ−1b− (σ + iωε)e = 0. (1.7)It is clear from the above paragraph that in typical geophysical parameter regimesσ >> ωε , making it safe to ignore the electric displacement.1.2.2 IP modelsThe constitutive law that really concerns us here is that between electric field andelectric current, which we use to model the IP effect. For common materials thatdon’t exhibit IP, Ohm’s law is typically used to describe how currents flow inresponse to applied electric fieldsj = σe (1.8)where the electrical conductivity σ may be a scalar or symmetric rank-2 tensor thatvaries in space but not in time or frequency. When IP effects are present, the flowof electric current will depend on the frequency of the applied electric field. Theseeffects can be modelled by replacing the standard Ohm’s law with a dispersiveconstitutive relation.Early IP modelling treated the phenomenon statically, as a perturbation of theDC resistivity problem. This was characterized mathematically by Seigel [85]. Hedescribed the polarization of a subsurface body as a perturbation of the conductiv-ity, writing Ohm’s law asj = σ(1−η)e (1.9)7where η is a unit-less number on the interval [0,1) called chargeability. The charge-ability quantifies the magnitude of the polarization effect, given sufficient time forcharge accumulation. This provided a powerful framework for the interpretationof IP data from the standard time-domain IP experiment described above. How-ever, Seigel himself noted that chargeability says nothing about how polarizationpotentials decay in time. In the 1970s, several authors showed how mineral typeswith similar conductivities and chargeabilities could be distinguished based on thespectral characteristics of their IP responses—see e.g. [71, 94, 107]. Pelton et al.[71] described this spectral IP effect using an Ohm’s law with a frequency de-pendent complex conductivity, where the frequency dependence is given by thephenomenological Cole-Cole modelσ(ω) = σ∞(1− η1+(1−η)(iωτ)c)(1.10)where ω is angular frequency, η is Seigel’s chargeability, τ is a time constant thatgoverns the overall time-scale of the relaxation, the exponent c controls the shapeof the spectrum, and σ∞ is the conductivity at infinite frequency or in the absenceof IP effects. This phenomenological model was introduced by Cole and Cole [22]to describe dispersive dielectric permittivity. Pelton et al. [71] had considerablesuccess in discriminating between minerals based on their Cole-Cole parametersin laboratory experiments. The Cole-Cole model has become the de facto standardmodel for characterizing the IP spectrum of geological materials. Zhdanov [105]proposed a more general model of induced polarization based on effective mediumtheory, known as GEMTIP, that includes the Cole-Cole model as a special case butit does not seem to be widely used. This may be because the Cole-Cole modelhas been seen as sufficient to fit IP data seen in the field or because it has notimproved the ability to correlate IP with geology. Given one of these models, or anyother frequency dependent conductivity model, frequency domain IP effects canbe modelled using standard methods for solving the frequency domain Maxwellequations designed without IP in mind. There are several examples of this approachin the geophysical literature, including [24, 100, 102].It should be noted that in the analysis of data collected in frequency domainIP surveys, it is not standard to model inductive effects by solving the full quasi-8static Maxwell equations. Typical experiments have traditionally been designed toreduce or eliminate the effects of EM induction. What induction did occur was con-sidered noise and had to be removed from field data using imperfect pre-processingprocedures. Classical examples of this approach include [71] and [79].Grounded source frequency domain surveys designed to detect IP effects arestill in use, especially for near-surface environmental work (see e.g. [54, 76]) buttime-domain surveys remain the more common way to collect data in the mineralexploration setting. Being able to extract information on the ground based on theshape and time-scale of the IP decay in a time-domain IP experiment, as well as be-ing able to fully model EM coupling effects motivates us to study how to efficientlymodel IP effects in the time-domain.A larger motivation for the EM geophysics community at large to study IP inthe time domain (and to study coupled IP and EM induction effects) is the increas-ing awareness of IP effects in airborne time-domain electromagnetic (ATEM) data.The well known phenomenon of negative transients in concentric loop ATEM data(see e.g. [56, 99]) cannot be modelled using non-dispersive ATEM modelling algo-rithms. More subtly, IP signals may mask conductive features at depth without thesmoking gun of negative transients. Some groups (e.g. [56, 87]) have focused onseparating induction and IP effects in ATEM data through data processing ratherthan capturing both phenomena in earth model based simulations. Uniquely, Kangand Oldenburg [53] have taken a compromise approach between data driven decou-pling and full time domain EM-IP simulation. They simulated IP affected ATEMdata as a linear perturbation of non-dispersive EM data, then used those perturba-tions as data to invert for chargeability.Unfortunately, modelling IP in the time domain is more complex than in fre-quency. The simple notion of Ohm’s law with a frequency dependent conductivityis lost and the constitutive law becomes a convolution. Most work in modellingtime-domain dispersive EM in time has attempted to avoid treating the convo-lutional Ohm’s law directly by performing the modelling in frequency domainat many frequencies and then transforming the results to the time-domain. Thiswas done in the suite of software released through AMIRA project P223F [100].More recently, Fiandaca et al. [32] used this approach in a scheme to invert ATEMdata for Cole-Cole parameters. Their inversion algorithm was computationally9tractable because they used one-dimensional forward modelling. The approach ofmodelling in frequency then transforming to time becomes very computationallyexpensive for 3D problems. Zaslavsky and Druskin [104] were able to to more effi-ciently compute time-domain responses from frequency domain information usinga Krylov subspace reduction technique.By contrast, there have been very few methods described that model time-domain IP by direct time-stepping. Recently, Commer et al. [25] modelled dis-persive TEM fields using an explicit time-stepping finite difference time-domain(FDTD) algorithm, approximating the Cole-Cole convolutional Ohm’s law usinga spectrum of Debye decays. Marchant [62] used an implicit time-stepping ap-proach, directly treating the convolutional Ohm’s law using numerical quadrature,with the convolution kernel transformed to the time-domain by the digital filteringtechnique. The only other approach we are aware of in the literature was developedby Marchant et al. [63]. In this approach the transformation of Ohm’s law to time ishandled using rational function approximations, which allow it to be described byan ordinary differential equation in time. This approach was recently extended byCai et al. [20] to make it more computationally efficient. In chapter 3 I will discussthe details of time domain EM and IP modelling in more detail and describe ournew approach to the problem.1.3 Thesis outlineThe first part of this thesis (chapters 2 and 3) considers forward simulation oftransient EM survey data, while the second part (chapters 4 and 5) discusses theinversion of such data. Chapter 2 focuses on the development of efficient timediscretization methods for the quasi-static Maxwell equations. IP is ignored inthis chapter. Modelling IP effects adds additional computational complexity to thestandard EM forward simulation problem. I develop efficient methods for the non-dispersive case before expanding to consider IP effects in chapter 3. The mainresult is a novel parallel time-stepping algorithm. I give a brief overview of thespatial discretization approach, followed by a discussion of discretization in timeand the development of the parallel time-stepping algorithm.Chapter 3 presents a novel approach to modelling IP effects in time-domain10electromagnetic data developed by my supervisor Dr. Haber and myself, based onstretched exponential relaxation. Stretched exponential relaxation has been used tomodel relaxations in disordered systems in other fields such as the study of glassysystems in solid-state physics. It has been suggested by Everett [31] and by Haber[37] as a way to model IP effects but has not been studied in detail in the geophys-ical community. I compare the stretched exponential model of IP to the widelyused Cole-Cole model and develop an efficient forward modelling algorithm usingthe stretched exponential approach. The capabilities of the algorithm are demon-strated using synthetic examples based on typical scenarios in grounded source IPsurveying and airborne EM surveying.Chapter 4 of the thesis describes preliminary investigations into inverting tran-sient EM data to recover the parameters that describe polarizable materials in thestretched exponential model. This is a difficult, multi-parameter, highly under-determined inverse problem. The work in chapter 4 presents a proof of conceptdemonstration that recovery of stretched exponential IP related parameters fromtransient EM data is possible but further work is needed in order to develop robustinversion algorithms and workflows suitable for use by practising geophysicists.Chapter 5 begins to address these issues. One approach to addressing the non-uniqueness in the stretched exponential inverse problem is to use parametric inver-sion. In a parametric inversion one assumes that the target of the inversion is a ho-mogeneous body embedded in a heterogeneous background medium and attemptsto recover its shape, position and physical properties—rather than trying to recoverphysical property values for each cell in a computational mesh. In chapter 5 I de-velop a stochastic parametric inversion algorithm that accounts for uncertainty inestimation of the heterogeneous background medium. I test the algorithm on a sim-ple linear two-dimensional (2D) inverse problem and on a 3D direct current (DC)resistivity inverse problem. Applying this approach to stretched exponential in-version is beyond the scope of this thesis but presents an excellent opportunity forfuture work.11Chapter 2Efficient modelling ofnon-dispersive TEM dataModelling the quasi-static time-domain Maxwell equations is a computationallyexpensive problem, even in the non-dispersive case. Adding induced polariza-tion (IP) effects further increases the complexity. Therefore efficient numericalmethods are essential to the successful application of coupled modelling of elec-tromagnetic (EM) and IP effects to large scale real world problems. A significantamount of research has investigated spatial discretization of Maxwell’s equationsover large domains and at high spatial resolutions—e.g domain decomposition([29, 75]) and multigrid methods ([38]). Advancing solutions in time, however,remains a bottleneck. This chapter addresses that bottleneck by developing a novelparallel in time forward modelling algorithm for the quasi-static Maxwell equa-tions. I consider only the non-dispersive case in this chapter before considering IPeffects in chapter 3.I start this chapter by describing the formulation of Maxwell’s equations Iwill solve and giving a brief overview of how the equations are discretized inspace. I use established methods for spatial discretization. I then discuss somegeneral considerations for time-discretization of Maxwell’s equations before de-scribing the parallel time-stepping algorithm for geophysical time-domain electro-magnetic (TDEM) simulations.122.1 Discretization of the time-domain quasi-staticMaxwell equationsThis section will describe how I’ve discretized the quasi-static Maxwell equationsto simulate time-domain electromagnetic geophysical surveys over non-dispersiveearth models. These surveys are governed by Faraday’s law and Ampe`re’s lawsupplemented by appropriate initial and boundary conditions. Recall Faraday’slaw and Ampe`re’s law∇× e =−∂b∂ t(2.1)∇×µ−1b− j = js, (2.2)where e is the electric field, b the magnetic flux density, j the electric current in theearth and js the source current. In this chapter I will neglect IP effects and assumethat currents in the earth are governed by Ohm’s law,j = σe, (2.3)where σ is the isotropic conductivity, a scalar function of space. SubstitutingOhm’s law into Ampe`re’s law eliminates explicit dependence on j. I eliminateexplicit dependence on b by combining eq. (2.1) and eq. (2.2) to give the singlevector-valued parabolic partial differential equation (PDE)∇× (µ−1∇× e)+σ ∂e∂ t=−∂ js∂ t. (2.4)This electric field equation is the preferred formulation of Maxwell’s equations formodelling grounded source surveys as it provides the most natural setting for dis-cretizing the sources and computing observable quantities from the electric fields.I consider this equation on a rectangular spatial domain Ω ∈ R3, on the closedtime interval [0, t∗]. I assume steady-state initial conditions. For grounded sourcesurveys the initial electric field is computed by solving the steady-state problem∇ ·σ∇φ = ∇ · js, (2.5)13where e(t = 0) = −∇φ . The initial electric field is zero in inductive source prob-lems. I impose natural boundary conditions on the electric fieldnˆ× (µ−1∇× e) = 0, (2.6)where nˆ is the unit normal vector to the boundary. The diffusive electric fields inquasi-static electromagnetics decay asymptotically toward zero away from energysources. In practice the boundary conditions are implemented approximately bychoosing a domain large enough such that the boundaries are far enough from anysources that 2.6 is approximately satisfied.Collecting the information from the last paragraphs, I can now state the initialboundary value problem that is to be solved∇× (µ−1∇× e)+σ ∂e∂ t=−∂ js∂ ton Ω× (0, t∗] (2.7a)∇ ·σ∇φ = ∇ · js on Ω× (−∞,0] (2.7b)nˆ× (µ−1∇× e) = 0 on ∂Ω. (2.7c)I take a method of lines approach to discretizing the problem. The main focusof this chapter is on modelling and utilizing the temporal information in TDEMgeophysical surveys. I have used existing numerical methods and software forspatial discretization. I will give a brief overview of the spatial discretization Ihave employed before moving on to describe time-discretization methods and theapplication to grounded source TDEM inversion.2.1.1 Spatial discretizationI perform spatial discretization using the mimetic finite volume method on locallyrefined rectilinear meshes called OcTrees. The general approach I use is describedin [37]. It is based on the mimetic finite difference algorithms developed by Hymanand Shashkov [49]. It was extended for use in modelling the quasi-static Maxwellequations on OcTree meshes in work by Haber and Heldmann [38] and by Horeshand Haber [45].14I discretize the weak form of eq. (2.7a),(µ−1∇× e,∇×v)+(σ ∂e∂ t,v) =−(∂ js∂ t,v), (2.8)where the parentheses indicate the L2(Ω) inner product, v is an arbitrary test func-tion from the Sobolev space H0(curl,Ω), and the curl operator is to be understoodin the weak sense. The definition of H0(curl,Ω) isH0(curl,Ω) = {u ∈ L2(Ω)3 : ∇×u ∈ L2(Ω)3,n×u = 0 ∈ ∂Ω}, (2.9)where L2(Ω)3 is the set of 3-component vector valued functions whose compo-nents are square integrable on Ω. A solution of eq. (2.8), called a weak solution,is a function in H0(curl,Ω) that satisfies eq. (2.8) for arbitrary test functions inH0(curl,Ω).The weak formulation of the electric field curl-curl equation admits a largerclass of solutions than the strong form. In particular it provides a natural setting tostudy non-smooth solutions. It also provides a convenient setting in which to for-mulate mimetic finite volume discretization. For a detailed and rigorous treatmentof the weak formulation of Maxwell’s equations the reader is referred to [65].The discretization of eq. (2.8) follows three steps. First the curl operator isdiscretized, followed by the inner products, giving a system of ordinary differentialequations (ODEs). Finally, the time derivatives are discretized to give an implicittime-stepping scheme to solve for the time evolution of the electric field.The electric field is approximated on mesh cell edges and its curl on face cen-tres, following the staggered grid approach introduced by Yee [103]. This is illus-trated in fig. 2.1. The use of staggered grids is standard in the numerical solutionof Maxwell’s equations. See e.g. [4, 23, 30] for examples from the geophysicalliterature using staggered grids to approximate Maxwell’s equations using spectral,finite difference, and finite element methods, respectively.I will briefly sketch the discretization of the curl operator. The curl of a vectorfield u may be defined as(∇×u) · nˆ = limA→0(1|A|∮∂Au ·dr)(2.10)15(i, j, k) (i+1, j, k)(i, j, k+1)(i+1, j+1, k)(i+1, j+1, k+1)ezexeybybxbzFigure 2.1: Yee grid. Components of the electric field are discretized on meshedges in corresponding Cartesian coordinate directions. Magnetic fieldsare discretized on face centres. Figure adapted from [51].where A is an area and the line integral is around the boundary of A. Let a f be thearea of a face fi of the finite volume mesh. The component of ∇×u normal to fiis approximated by(∇×u) · nˆf ≈ 1a f mid(u) (2.11)where mid(u) is the midpoint rule approximation of the line integral of u aroundthe boundary of the face fi. Applying this over the entire mesh gives the discretecurl operator C which maps an edge grid vector to a face grid vector.The inner product integrals are approximated using the three-dimensional (3D)trapezoidal rule. For a detailed derivation see [39]. This leads to the implicit systemof ODEsCT Mµ−1Ce+Mσ∂e∂ t=−∂ s∂ t(2.12)where e is now understood to be an edge grid vector and s is the discretization ofthe source current. The M matrices are symmetric positive definite mass matricesarising from discretization of the inner products in eq. (2.8). Their subscripts referto the PDE coefficients involved in the respective inner products.I have used the julia software package jInv ([80]) to implement the spatial dis-16cretization. Julia is a relatively new dynamic language designed for high perfor-mance technical computing [15]. It allows rapid prototyping, with a syntax similarto Matlab, but gives high performance and has built in distributed memory paral-lelism that allows for code that scales to large distributed memory systems.jInv is a software package and framework for PDE constrained parameter es-timation problems. It provides finite volume discretizations of differential opera-tors, other building blocks for developing software for the numerical solution ofpartial differential equations, and optimization routines designed specifically forPDE constrained parameter estimation problems. It uses Julia’s builtin distributedmemory parallelism and is designed to scale to large scale industrial problems. Iwas not involved in the initial development of jInv but have become a contribu-tor to the core jInv package, a primary maintainer of the JOcTree package whichprovides support for finite volume discretizations on OcTree meshes, and the pri-mary developer of the open source transient EM module used to run transient EMforward modelling and inversions in the jInv framework.2.1.2 Discretization in timeHaving approximated eq. (2.8) by the system of ODEs eq. (2.12), I now study itsdiscretization in time using backward differentiation formulas (BDFs). The quasi-static Maxwell equations in time are extremely stiff [40]. It is difficult to give aprecise and rigorous definition of stiffness of ODEs. Ascher and Greif [8] give thefollowing intuitive definition:The initial-value ODE problem is stiff if the step size needed to main-tain absolute stability of the forward Euler method is much smallerthan the step size needed to represent the solution accurately.Stiffness is a common property of systems of ODEs that result from the discretiza-tion of parabolic PDEs such as the curl-curl equations derived from the quasi-staticMaxwell equations. This seems to occur because it is often the case that the so-lutions of these equations have secondary features that decay extremely quicklywhile their primary features, that are typically of interest, evolve more slowly.Because of stiffness the choice of time-stepping method is dictated as muchby stability considerations as accuracy considerations I limit my attention to time-17stepping methods that maintain absolute stability for all step-sizes. Such methodsare called A-stable. All A-stable methods are implicit, meaning that a system oflinear or non-linear equations must be solved at each time-step. This makes themmuch more expensive per step than explicit methods, in which the solution at agiven time can be computed explicitly in terms of known quantities. However, forstiff problems, the much larger step-sizes that can be taken with implicit methodsmore than compensate for the increased cost of each step.In addition to A-stability, time-stepping methods for some very stiff problems,including the quasi-static Maxwell equations, must possess a property known asL-stability. Roughly speaking, L-stable methods are methods that when applied toODEs with solutions that decay in time, will produce approximate solutions thatdecay arbitrarily close to zero in a single step as the step-size grows toward infinity[6, 7]. A-stable methods that are not L-stable will produce decaying solutionsfor all step-sizes but may introduce spurious solution components that decay veryslowly. I have seen in practice that such methods (e.g. trapezoidal rule, TR-BDF2)fail for transient electromagnetic geophysics problems.I use the class of A and L-stable methods known as BDFs. I have used thefirst and second order BDF methods. Higher order methods could be used but inthe transient EM geophysics use case the benefit is not likely to be large enoughto outweigh the additional code complexity incurred, especially when consideringsensitivity computations in inverse problems. The simplest BDF method is thebackward Euler algorithm. It approximates the solution of the initial value problemy˙(t) = f(t,y) (2.13)y(t0) = y0 (2.14)at discrete times {tn+1} byyn+1 = yn+∆tnf(tn+1,yn+1), (2.15)where ∆tn = tn+1− tn. Applying backward Euler to eq. (2.12) gives a system of18linear equations to be solved to update the electric field at each time-step(CT Mµ−1C+1∆tnMσ )en+1 =1∆tn(Mσen+ sn− sn+1). (2.16)For second order time-stepping I used the variable step-size and fixed leadingcoefficient forms of the second order backward differentiation formula (BDF2).Defining ρn = ∆tn/∆tn−1, the variable step-size BDF2 method is1+2ρn1+ρnyn+1 = (1+ρn)yn− ρ2n1+ρnyn−1+∆t f (yn+1, tn+1). (2.17)The method maintains stability and second order accuracy for ρ < 1+√2 [43].Applying this method to eq. (2.12) gives the linear time-stepping system of equa-tions(CT Mµ−1C+1+2ρ∆tn(1+ρ)Mσ )en+1 =1∆tn(Mσ[(1+ρ)en− ρ21+ρen−1]−1+2ρ1+ρsn+1+(1+ρ)sn− ρ21+ρsn−1). (2.18)Two main factors complicate the use of BDF2 in TDEM simulations. The firstdifficulty is that the source current in TDEM surveys is usually discontinuous intime. Secondly, one is usually interested in simulating TDEM fields on logarith-mic time scales. Previous work [68] has shown that simulating TDEM geophysicalsurveys over such time-scales is most efficiently accomplished by using sparse di-rect matrix factorization methods to solve eq. (2.16). The time-range of interestis divided into regions, with the step-size kept constant in each region. For exam-ple, fig. 2.2 shows a synthetic transmitter waveform—the current in the transmitteras a function of time—similar to that employed in the commercial SkyTEM TDEMsystem [88]. A single step-size is used to discretize the on-time, during which mag-netic fields in the earth are developed. The response to the shutoff of the transmitteris simulated over a time range from 0.01 ms to 6 ms after the end of the transmittershutoff, using six distinct step-sizes. The backward Euler system matrix must befactored once for each unique step-size.Like the backward Euler system, the variable step-size BDF2 system matrix19−0.006 −0.004 −0.002 0.000 0.002 0.004 0.006time (s)0.00.20.40.60.81.0transmitter current (A)10−510−4step-size (s)Figure 2.2: Typical example of realistic TDEM transmitter current time-dependence plotted in black, with dots showing times at which solutionis computed. The solid blue line shows the variation in the step-size overthe full time range of the simulation. There are six unique step-sizes.eq. (2.19) remains constant while the step-size is constant. However, the systemmatrix at a transitional step between a time region where a constant step size ∆t1was and ∆t2 is different from the matrix for subsequent steps of size ∆t2. Matrixfactorization is computationally expensive and for TDEM simulation methods tobe efficient, it is essential to minimize the number of factorizations. There aretwo main approaches in the ODE literature to avoid factoring an additional matrixwhen the step-size changes (see e.g. [6, 43]). In the first approach, the fields atprevious times may be interpolated to compute the solution at the time that willallow the constant step-size BDF2 method for the new step size to be used at thestep-size boundary. In the second approach, which I have used, the BDF2 method isreformulated in so-called fixed leading coefficient form. This is more robust thanthe interpolation approach. In the fixed leading coefficient form, the coefficientof the mass term at the leading time is replaced by one that depends only on thecurrent step-size, at the cost of an additional term in the right hand side of the time-stepping linear system. Using the fixed leading coefficient second order backwarddifferentiation formula (FLCBDF2) in the discretization of eq. (2.16) gives the20time-stepping system of equations(CT Mµ−1C+32∆tnMσ )en+1 =1∆tnMσ[32(1+ρ2/3)en− ρ22en−1]−−dsn+1+ ρ−12(CT Mµ−1Cen+dsn). (2.19)where dsn+1 denotes the second order backward approximation to the time deriva-tive of the transmitter current at time tn+1.For efficiency, it is crucial that BDF2 time-stepping is initialized without per-forming an additional factorization. BDF2 is a multi-step method, with compu-tation of each step requiring the approximate solution at the two previous steps.Therefore the first step must be computed using a different algorithm. Standardpractice is to initialize BDF2 with a single step of backward Euler. Unfortunately,this requires an additional factorization, since the system matrix for backward Eu-ler time-stepping is different from the BDF2 matrix. However, the constant step-size and fixed leading coefficient BDF2 system matrices for step-size ∆t are equalto the backward Euler system matrix for step-size ∆tbe = 2∆t/3. I have initializedBDF2 time-stepping by taking two backward Euler steps of size ∆tbe and linearlyinterpolating between them to compute the fields at time t0 +∆t, where t0 is thetime at which initial conditions are specified. The additional backward Euler stepneeded for interpolation is much cheaper to compute than an additional factor-ization and the interpolation does not cause a degradation in the overall order ofaccuracy of the BDF2 method.Now I address the treatment of discontinuities in the electric field. BDF2,and indeed all higher order time-stepping methods, lose their convergence prop-erties in the presence of discontinuities in the solution [89]. In TDEM surveying,electrical currents in the earth are induced by rapid changes in the magnetic fieldof the transmitter. The change in magnetic field is achieved by shutting off thetransmitter current as quickly as possible. In cases where the transmitter on-time isshort, and data will be collected in early times after the transmitter shutoff, accuratemodelling of the full transmitter current time-dependence, such as that depictedin fig. 2.2 is necessary. In other cases, the transmitter current may be modelled asI(t) = 1−H(t), where H(t) is the Heaviside step function. Such a waveform is21−1.0 −0.5 0.0 0.5 1.0 1.5time (s)1e−40.00.20.40.60.81.0transmitter current (A)Figure 2.3: Transmitter current step-off waveform. Assumes transmitter cur-rent has been on long enough that fields are steady-state at t = 0, whereforward modelling begins. Shutoff takes place in one time-step.illustrated in fig. 2.3. Regardless of the exact waveform, the time-derivative of thesource current, which appears in the forward modelling equation eq. (2.16), is typ-ically discontinuous at the point where the current reaches 0. Fortunately the pointof discontinuity depends only on the transmitter waveform and so is known a pri-ori. In simulation of actual geophysical surveys the waveform will be known onlyapproximately but the approximate waveform will be taken as a fixed input to theforward modelling algorithm—thereby making the point of discontinuity a fixedinput to the algorithm. This makes it straightforward to take steps to minimize itseffect. I have found empirically that when no effort is made to address discontinu-ities, much error is introduced into the solution near the time of the discontinuitybut that the accuracy at late times is not severely affected.Because BDF2 is a backward looking two step method, it will use informationfrom both sides of a discontinuity when taking the first step immediately after thediscontinuity, giving spurious results. For example, in the approximate step-offwaveform in fig. 2.3 the true time derivative of the current at 1×10−4 s is clearly220 A/s. However, the BDF2 approximation of the derivative isdIdt≈ 32I(1×10−4)−2I(5×10−5)+ 12I(0) =12. (2.20)I address this problem by reverting to backward Euler time-stepping for twosteps after crossing a discontinuity so that the solution at the time of the disconti-nuity is not used in any BDF2 steps. In the case of a step-off transmitter I simplytake the first three steps using backward Euler.I use a simple concentric loop survey configuration above a homogeneousearth—for which an analytic solution is available—to empirically study the effec-tiveness of this approach. This model problem consists of a circular loop transmit-ter of radius r, placed at the surface of a homogeneous half space of conductivityσh. Assuming that a steady current runs through the transmitter for times t < 0,and that the current is instantaneously terminated at t = 0, the time derivative ofthe vertical component of the magnetic field for t > 0 is∂hz∂ t=1µ0σhr3(3erf(τ)− 2√piτ(3+2τ2)e−τ2), (2.21)where τ = r√σhµ04t , erf is the error function, and µ0 is the magnetic permeability offree space [97]. This idealized problem is based on the typical airborne transientEM survey in which a horizontal induction coil at the centre of a horizontal looptransmitter measures the ∂hz∂ t response when the current in the loop is abruptly shutoff.I simulated the electric fields due to a loop with radius 13.5 m at the surface ofa half-space with a conductivity of 0.01 S/m using backward Euler time-stepping,BDF2 initialized with a single step of backward Euler, and BDF2 initialized withthree steps of backward Euler. The spatial discretization was the same in all threecases. It used an OcTree mesh consisting of a core region of 5 m3 fine cells cov-ering a rectangular volume of side length 200 m centred at the loop centre. Thefields were simulated from 0.005-1 ms using 200 steps of length 0.005 ms. Thetransmitter current was set to 1 A at the initial time t = 0 (immediately prior to theinstantaneous shutoff) and 0 A for all subsequent times. The z-component of themagnetic field dhz/dt was computed from the electric fields at each time-step by23numerically computing the z-component of ∇× e at the centre of the transmitterloop. Because spatial discretization adds additional error to the numerical solution,the analytic solution eq. (2.21) cannot be used to perform a rigorous convergencestudy of my implementations of the backward Euler and BDF2 time-stepping algo-rithms. However, by using a relatively fine spatial mesh much can be learned fromcomparisons with the analytic solution.The true solution and computed dhz/dt from all three time-stepping methodsare shown in fig. 2.5a. All three methods are wildly inaccurate in the first few time-steps. This is the expected behaviour. In these early steps the solution is dominatedby the approximation of the transmitter step-off. Fortunately these artefacts arequickly damped and at later times the solution is less affected by the accuracy ofthe current waveform approximation. This is consistent with the general behaviourof BDF methods applied to parabolic problems. In these problems it is not neces-sary to resolve the details of early time behaviour in order to capture the late timebehaviour [6]. In practice, when approximating a real transmitter waveform by astep-off, one must simply choose a small enough step-size such that any artefactsfrom inaccurate approximation of the step-off are sufficiently damped before theearliest time of interest.After these early steps dominated by artefacts, both BDF2 solutions are muchmore accurate than the backward Euler solution. This can be seen by eye in fig. 2.5aand in the relative errors in the approximate solutions plotted in fig. 2.4b. Fig-ure 2.4b also shows that initializing the time-stepping with three backward Eulersteps in order to avoid using BDF2 across the discontinuity in the solution does im-prove the accuracy at later times, as well as eliminating the large oscillation seen inthe naive solution where BDF2 is applied across the discontinuity. This oscillationcan be more clearly seen in fig. 2.4c, in which the solutions are plotted on a linear,rather than a logarithmic scale. At very late times the two BDF2 solutions do con-verge but this is likely because other sources of error, such as those due to spatialdiscretization, are starting to dominate the discrepancy with the true solution.2410−510−410−3time (s)10−910−810−710−610−510−410−310−2dhz/dt(nT/m) rue solu ionbackward EulerBDF2 1 s ep ini BDF2 3 s ep ini (a) dhz/dt, plotted on log scale10−510−410−3time (s)10−210−1100101102relative error in dh_z/dt(b) Relative error0.5 1.0 1.5 2.0 2.5time (s)1e−5−0.004−0.0020.0000.0020.0040.0060.0080.010dhz/dt(nT/m)(c) dhz/dt for first 5 time-steps, plottedon linear scaleFigure 2.4: Comparison of time-stepping methods in concentric loop exam-ple. In the legend, BDF2 1 step init refers to the solution computedusing backward Euler for the first step only, using BDF2 thereafter withno special treatment of the discontinuity. BDF2 3 step init refers tothe solution computed using backward Euler time-stepping for the firstthree steps and BDF2 stepping thereafter. Dashed lines indicate negativevalues. a) shows the overall behaviour of each method over the full 200time-step simulation. c) illustrates the large oscillation observed fromthe second time-step when using BDF2 with one step initialization.252.2 Parallel time-steppingThe computational cost of transient geophysical EM forward modelling is dom-inated by the solution of the linear system of equations needed to advance theelectric fields at each time-step. The process is inherently sequential. As discussedabove, sparse direct methods have proven to be more efficient than iterative meth-ods for geophysical transient EM modelling. It is possible to use parallel sparsefactorization software but the factorization computations have limited parallel scal-ing.With these considerations in mind I seek alternative methods for the parallel nu-merical solution of the stretched exponential forward modelling problem and otherparabolic initial-boundary value problems. Parallelism across time has receivedmuch less attention from the numerical analysis community than parallelism inspace but a body of literature on the topic does exist. The greatest effort has beenin the acceleration of high order time-stepping methods. The prototypical algo-rithm of this class is the parareal method [59]. These methods are likely of littleuse here since they require serial execution of a low order time-stepping method tofacilitate parallelization of the high-order method. Here I am seeking to parallelizea first or second order method. Other methods such as space-time multigrid [46]and multigrid waveform relaxation [47] treat the entire space-time problem all atonce, which is difficult for large scale problems. Additionally, it is not clear howto reuse information from forward modelling in derivative computations with thesemethods. Relatively recently, parallel exponential integrator methods, such as themethod presented by Bo¨rner et al. [16], have shown promising performance andpotential for reuse of information during derivative computations but they are dif-ficult to differentiate. It may be that their parallelism can scale to large numbers ofprocessors but this is not clear to us.Here I propose a parallel time-stepping method for time-domain electromag-netic geophysics that is a small adjustment of the methods already discussed in thischapter. The key insight leading to the method is that due to the diffusive nature ofquasi-static transient EM fields it is possible to use independent time stepping pro-cesses to simulate electric fields at different time-scales of interest. For parabolicPDEs, with sources that step off or ramp off in time, it is possible to model late26time behaviour in a stable fashion without resolving early-time behaviour [6]. Fea-tures of the fields decay in magnitude and become spatially smoother as the timeafter source shutoff increases. For example, to simulate a survey collecting dataover the time interval 0.001-1 s, a minimum step-size of 1×10−4 or smaller wouldbe required but a step-size of 0.05 s might be sufficient to roughly approximate theelectric field 1 s after current shutoff.In my parallel time-stepping algorithm I consider modelling the electric fieldat each time of interest as an independent forward problem. I then choose a sin-gle time step size and number of steps adequate to accurately model the field atthat time. Then each time-stepping process has a single step-size. Recall that thematrix of the time-stepping linear system remains constant while the step-size re-mains constant, so each independent forward problem in the parallel time-steppingscheme requires only a single factorization, and the time-stepping processes aretrivially parallel.I have tested this parallel time stepping approach only for ideal step-off sources.As just described it should also work for more complex transmitter current time-dependence (waveform), provided that the source waveform can be adequately dis-cretized in time using the same step size required for simulating late time fields. Incases where the source waveform varies on a shorter time-scale than the step-sizerequired for simulating late time fields, a more sophisticated scheme would needto be used in order to achieve good parallel scaling. A possible approach would bethat described by Gander and Gu¨ttel [34]. The reader is directed to that work fora detailed explanation of the method. In their method, an inhomogeneous linearinitial value problem is split into two independent sub-problems that can be solvedin parallel, one homogeneous and one inhomogeneous. In this manner a small steptime-stepping process incorporating the source behaviour may be run in parallelwith a large step process run until the time of interest. The solution of the trueproblem is then given by the sum of the independent sub-problems. I have not im-plemented it but I believe my parallel time-stepping approach could be augmentedby such a scheme in order to extend it to more complex transmitter waveforms.The main determining factor of the parallel efficiency of the method is thenumber of time steps, which I call spin-up steps, needed to simulate the electricfield at a given time. If a very small step-size—and therefore an excessive number27of steps—are needed to ensure accuracy at late times then the benefit of performingmultiple factorizations in parallel will be greatly reduced. Although knowledge ofearly time-behaviour is not required to resolve late-time behaviour, the accuracy ofthe electric field at any time t will be proportional to the step-size or square of thestep-size for backward Euler and BDF2, respectively. I have found by numericalexperiment that the number of backward Euler time steps required to resolve latetime fields with sufficient accuracy was large enough to significantly dissipate thespeedup acquired through parallelization. The benefit of parallel time-stepping isgreatly improved by the use of second order time-stepping methods.2.2.1 Evaluating accuracyThe parallel BDF2 time-stepping algorithm was tested with a range of exampleproblems to get an empirical sense of the typical number of spin-up iterations re-quired to achieve sufficient accuracy over a wide range of time-scales. For each testproblem I computed a reference backward Euler solution using a very small step-size and then computed the relative error between the reference solution and BDF-2solutions with various constant step-sizes. Results for a typical test problem usinga reference solution with BE step-size of 1× 10−5 s are shown in figure fig. 2.5.The test problem used to generate those figures was that of computing the electricfield at the surface of the earth due to a 2 km long straight wire grounded trans-mitter over a homogeneous half-space. The transmitter waveform was a simplestep-off. The mesh consisted of a core region of 50× 50× 50 m cells coveringthe transmitter and a 1.2×1.2 km area centred about the centre of the transmitter.fig. 2.5a shows the component of the electric field parallel to the transmitter wire,at the point half way between the transmitter electrodes. It shows that for the firsttwo to three time-steps for each BDF2 step size, the fields differ significantly fromthe reference backward Euler solution, as expected. A step-size of 5× 10−4 s forexample is not small enough to accurately resolve the electric field at 10−3 s. How-ever, the fields for all BDF2 step sizes start to agree with the backward Euler fieldafter four steps. Figure 2.5b shows the relative difference between the reference2810−410−310−2time (s)10−410−310−210−1100voltage (V)BE dt=1×10−5BDF2 dt=1×10−4BDF2 dt=5×10−4BDF2 dt=1×10−3BDF2 dt=2.5×10−3BDF2 dt=5×10−3(a) Electric field in-line component over the time in-terval [10−4,10−1] s for reference backward Euler solu-tion computed using step-size of 1× 10−5 s and BDF-2solutions at various step-sizes shown in the legend.5 10 15 20number of time steps10−510−410−310−210−1100101Relative difference(b) Relative error between reference backward Eulersolution and BDF-2 solutions for each BDF-2 step-size.Figure 2.5: Parallel time-stepping accuracy on grounded source test problem.29backward Euler solution and the BDF2 solutions, which is computed point-wise asrel. dif. =d−drefdref. (2.22)For all step-sizes the BDF2 solutions have a relative difference with the referencesolutions of less than 1× 10−2 after 6 time steps. In the context of inverse mod-elling such accuracy is normally sufficient. Errors in the physical property modelswill overwhelm discrepancies from this parallel time-stepping scheme. Testing onmore complex conductivity models did not show a degradation in accuracy.2.2.2 Evaluating performanceI tested the performance of the algorithm on a synthetic problem with a more re-alistic time-stepping scheme. This problem is based on the gradient array styleof survey used in induced polarization experiments. It includes two perpendicularlong grounded wire transmitters that intersect at their centres. An array of electricdipole receivers is laid out on a regular grid centred at the point of intersection ofthe transmitter wire paths. The volume enclosing the survey area was meshed witha core volume of 50×50×25 m cells, with the cell size gradually expanding fromthere to fill the domain. With the domain being made large enough to approxi-mately fulfill the zero-flux boundary conditions. In total the mesh had 41120 cellsand 135233 edges. Constraining the so-called hanging edges of the mesh (a detailof spatial discretization not discussed in this chapter) then gave a linear systemwith 112849 unknowns to be solved at each time-step.This single mesh was used to simulate the data from both transmitters. Thelinear systems of equations at each time-step of the forward modelling processwere solved using the Pardiso sparse-direct linear solver software package [84].As discussed earlier in this chapter, factorization of the system matrix is the mostexpensive element of transient EM simulations when using direct methods. Thesystem matrix depends on the mesh, earth model and time-step size but not on thetransmitter or the time itself. The transmitter enters the system only through theright hand side. Therefore multiple transmitters and multiple time-steps of equalsize can be handled very efficiently by a direct solver—each additional transmitterand equally sized time-step adds a forward and backward substitution step but not300 10 20 30 40 50step number0.000050.000100.000150.000200.000250.00030step-size (s)Figure 2.6: Time-stepping scheme for reference backward Euler solution .There are 10 steps each of sizes 2×10−5 s, 4×10−5 s, 8×10−5 s, 1.6×10−4 s, and 14 steps of size 3.2×10−4 s.an additional factorization. Previous surveys of sparse matrix factorization soft-ware such as those by Gould et al. [36] and by Davis et al. [26] have found Pardisoto be among the best sparse factorization software in terms of performance and us-ability. My own previous work ([9]) and the experience of my research group havefound it to be the fastest such software for factorizing matrices and solving linearsystems arising in geophysical EM simulations.The current waveforms of both transmitters were modelled as step-offs. Re-ceiver voltages were computed at 22 time steps in the interval [2.8× 10−4,7.5×10−3] s. In the reference backward Euler solution this was modelled using 54 stepsof 5 unique sizes, as illustrated in fig. 2.6. Due to its lower accuracy, a larger num-ber of steps before the first observation time are needed to eliminate initializationartefacts in backward Euler time-stepping than with BDF2. In this example 12steps were taken to reach the first observation time.I performed parallel BDF2 time-stepping on the same mesh to compute the re-ceiver voltages at the same 22 observation times. I did not have sufficient computa-tional resources to test using a separate time-stepping process for each observationtime but I was able to divide the 22 observation times into 6 parallel time-steppingprocesses that each used could be performed in parallel, with each using a constant31step-size first obs. time last obs. time # of obs. times total # of steps3.5×10−5 2.8×10−4 4.4×10−4 3 136.5×10−5 5.2×10−4 7.6×10−4 3 131.15×10−4 9.2×10−4 1.4×10−3 4 132.15×10−4 1.72×10−3 2.68×10−3 4 133.75×10−4 3×10−3 4.28×10−3 3 136.15×10−4 4.92×10−3 7.48×10−3 5 13Table 2.1: Parallel time-stepping scheme. Six parallel forward modellingsimulations with constant step-sizes shown in this table were run in par-allel to simulate the receiver voltages at the same observation times usedin the reference backward Euler simulation. The number of observationtimes per parallel process varied but the total number of time-steps wasthe same in all the simulations.step-size and thus only a single factorization. For each parallel time-stepping pro-cess, 8 steps were taken to reach the first observation time. The step sizes used andobservation times covered by each parallel time-stepping process are summarizedin table 2.1. These tests were performed on a small computer cluster containingfour twelve core compute nodes, with two nodes running one time-stepping pro-cess each and the other two nodes each running two time-stepping processes. At-tempting to run more than two parallel time-stepping processes per node resultedin significantly reduced performance. Each time-stepping process is both computeand memory intensive and this poor performance was presumably due to low mem-ory bandwidth on the cluster nodes causing the computations to become memorybound rather than compute bound. The specific time-stepping scheme detailed intable 2.1 was chosen to minimize the number of time-steps per parallel processwhile balancing their computational loads.The reference backward was computed in approximately 89 s on a single coreof a single node of the cluster. This time, as with all timings in this performancetest, did not include initial setup computations such as the construction of the meshand differential operator matrices. This method of measuring performance waschosen because it is most representative of the typical cost of forward modellingduring an inversion. The setup computations need only be done once at the begin-ning of an inversion and subsequent simulations with updated earth models need32# of threads run time (s)1 89.2 62.4 46.6 40.12 36.Table 2.2: Backward Euler solution run times.only update the mass matrices that depend on the earth model.As previously mentioned in this chapter, aside from running parallel time-stepping processes, the linear solver itself can be parallelized. The Pardiso solveruses openMP shared memory parallelism. Using shared memory limits the parallelscaling to a single node but provides efficient parallelism within a node. Distributedmemory parallel sparse matrix factorization software that can scale across multi-ple computers, such as the MuMPS package [3], are available. Previous experiencewith MuMPS seems to indicate that good parallel scaling across multiple machineswill only occur for the largest of geophysical EM simulation problems.The parallel time-stepping algorithm provides a way to expand the parallelismof transient EM forward modelling beyond what a parallel sparse linear solver canprovide. In this example the best performance was achieved when combining thetwo modes of parallelism. To achieve a full comparison, I ran both the referencebackward Euler solution and the parallel time-stepping BDF2 solutions with Par-diso in both single and multi-threaded modes.I ran the reference backward Euler solution with 1, 2, 4, 6, and 12 threads.Recall that this simulation involved 54 time-steps of 5 unique sizes, thus requir-ing 5 factorizations. The run times for these tests are shown in table 2.2 and thecorresponding parallel efficiency plot is shown in fig. 2.7. Using Pardiso gives anon-trivial speedup but the parallel efficiency drops off precipitously as the numberof threads increases.The parallel time-stepping scheme achieved lower run times than was possi-ble by parallelizing the linear solver. This was the case despite the fact that theparallel time-stepping solution actually performed more floating point operationsthan the backward Euler solution, using a total of 78 time-steps and six factoriza-332 4 6 8 10 12# of threads020406080100parallel efficiency (%)Figure 2.7: Pardiso parallel scaling results when using a single sequentialbackward Euler time-stepping process. Parallel efficiency for N threadsis T1/(N ·TN) where T1 is the sequential run time and TN is the run timeusing N threads.threads per process run time (s)1 202 154 116 11Table 2.3: Parallel time-stepping run times using 6 parallel time-steppingprocesses run across 4 compute nodes of a computer cluster.tions across all six processes in comparison with the 54 steps and 5 factorizationsused by the backward Euler solution. The parallel time-stepping algorithm is ableto take advantage of distributed memory high performance computing hardware,while the sequential time-stepping algorithm is not. With each process runningsingle threaded, the parallel time-stepping run time was 20 s—already a substan-tial improvement over single time-stepping process running with 12 threads. Thiswas improved by running Pardiso with multiple threads in each time-stepping pro-cess. The results are summarized in table 2.3. These results show that the paralleltime-stepping method is capable of significantly improving upon the parallelismpresent in matrix factorization software.342.3 ConclusionsThe work presented in this chapter has shown that the speed and accuracy of geo-physical transient EM forward modelling can be significantly improved by the useof more sophisticated numerical methods than are typically used in the temporaldiscretization of the problem. I have investigated the use of second order time-stepping methods in initial boundary value problems involving the time-domainquasi-static Maxwell equations. I found that as long as care is taken in initializingthe time-stepping and handling discontinuities, then these methods provide supe-rior performance, relative to first order methods. The use of second order meth-ods allowed the development of a parallel time-stepping algorithm that removed amajor bottleneck in the parallel scaling of quasi-static transient EM forward mod-elling. I demonstrated that the parallel time-stepping method could achieve signifi-cantly greater speedup than parallelization of the linear algebra of the time-steppingalone could achieve.35Chapter 3Forward modelling of thedispersive time-domain MaxwellEquations using the stretchedexponential functionThis chapter will describe the algorithm I have developed for modelling time-domain quasi-static electromagnetic (EM) fields in the presence of polarizableearth materials. I start by stating the forward modelling problem precisely, dis-cussing the difficulties introduced by the consideration of induced polarization (IP),and motivating the choice to use stretched exponential relaxation to model IP de-cays. I then discuss the similarities and differences between the stretched exponen-tial model of IP and the Cole-Cole model, the de facto standard frequency domaingeophysical IP model, before describing the implementation of the stretched ex-ponential forward modelling algorithm. The chapter finishes with synthetic mod-elling examples illustrating the capabilities of the stretched exponential approach.3.1 The stretched exponential formulationI am interested in solving the quasi-static time-domain Maxwell equations to de-termine the electric field due to either grounded or inductive sources, in the pres-36ence of polarizable media. The polarization effect is contained in the constitutiverelation between the electric current j and the electric field e. I formulate the quasi-static Maxwell equations in the time-domain in terms of the electric field e, mag-netic flux density b, current in the earth j and source current js. I use zero-fluxboundary conditions. Other boundary conditions could be used but zero-flux weresufficient for the work considered in this thesis. I will assume steady state initialconditions. Formally, the forward modelling problem is described by the follow-ing initial boundary value problem on a cuboidal domain Ω ∈R3 and time intervalt ∈ [0, t∗].∇× e =−∂b∂ ton Ω× (0, t∗) (3.1a)∇×µ−1b− j = js on Ω× (0, t∗) (3.1b)∂e∂ t= 0 on (−∞,0] (3.1c)nˆ× (µ−1∇× e) = 0 on ∂Ω× (0, t∗). (3.1d)For inductive sources, I assume the electric field is zero initially and for groundedsources, I compute the initial electric field by solving the DC problem for the givenelectrode configuration. The current in the earth j is the primary concern here. Inthe non-dispersive case, explicit dependence on j may be eliminated by the useof the standard constitutive relation j = σe, where σ might be a scalar or tensorquantity but depends only on position.As briefly discussed in the introduction to this thesis, the relationship betweenelectric field and current is less straightforward when IP effects are considered. TheIP phenomenon is most often described in the frequency domain, as a frequencydependent conductivity. Here Ohm’s law remains a simple linear relationshipJ = σ(ω)E. (3.2)When transformed to the time-domain, this becomes a convolutionj(t) =∫ ∞−∞σˆ(t− τ˜)e(τ˜)dτ˜ (3.3)37where σˆ(t) is the impulse response of the conductivity. Due to causality and theinitial condition on the electric field for t ≤ 0, this may be reduced to [62]j(t) = σ∞e(t)−∫ t0σˆ(t− τ˜)e(τ˜)dτ˜. (3.4)Note that I have not assumed any particular frequency dependence here. The Cole-Cole model is the most popular choice but eq. (3.4) applies generally to any fre-quency dependence. Making a different choice will simply lead to a different im-pulse response function.The question then becomes, how to deal with this convolution. This questionhas seen some recent attention in the geophysical geophysical EM community andin the electrical engineering finite difference time domain (FDTD) modelling com-munity. As noted by, e.g Commer et al. [25], Takayama and Klaus [91], the effortsfall into three main categories: direct numerical evaluation of the impulse responseand convolution integral (e.g. [62, 101]), auxiliary differential equation methods(e.g. [63, 74, 91]), and methods based on the Z-transform (e.g. [90, 98]).These methods are all based on translating frequency domain characteriza-tions of dispersive electromagnetic media to the time-domain and they all haveadvantages and disadvantages whose importance depends on the application. Forexample, Z-transform methods can be made highly efficient only when constanttime-stepping is employed, making them inappropriate for geophysics when thevariation in EM fields over logarithmic time scales is typically of interest. Auxil-iary differential equation methods work well when the frequency dependence has asimple enough form but require the use of fractional derivatives for more complexscenarios. Direct convolution methods are quite general but can be computationallyexpensive and difficult to differentiate with respect to model parameters.Another approach is to consider a parametrization of the IP effect defined natu-rally in the time-domain. This is the approach I take in this chapter. I use stretchedexponential relaxation to describe the decay of IP currents. The stretched exponen-tial functionf (t) = e−tβ /τ , (3.5)where β ∈ (0,1] and τ is a time-constant, describes the relaxation of complex sys-38tems whose multiple components exhibit exponential decays of different time con-stants [2]. It has been used widely in solid state physics to describe relaxations inglassy systems [2], and also in other areas of physics such as time-resolved lumi-nescence spectroscopy [13] and the diffusion of water in biological tissues [12]. Inthe application to glassy systems it is often thought of as a time-domain analogueto frequency domain representations similar to the Cole-Cole model [2, 44]. Thisexperience in the physics literature motivates the choice to explore the use of thestretched exponential to describe geophysical IP effects. Everett [31] and Haber[37] have noted how the stretched exponential may be used to describe the de-cay of IP currents but to my knowledge no electromagnetic simulation algorithmsbased on the idea have been developed.The algorithm developed in this chapter incorporates stretched exponentialrelaxation into Maxwell’s equations using an auxiliary differential equation ap-proach. The implementation originated from an auxiliary differential equationalgorithm for the special case of Debye dispersion (Cole-Cole with c = 1) byMarchant et al. [63]. Here I generalize the approach and make key changes toimprove the stability of the resulting numerical algorithm. I will show that formu-lating the constitutive law as an ordinary differential equation (ODE) is mathemat-ically equivalent to the convolution representation.Before describing the implementation of the stretched exponential approach, Iwill discuss the general approach to defining the j(e) constitutive law via an ODEand how it relates to convolutions. I start with two definitions. First, I introducethe effective DC conductivity σ0 = σ∞(1−η). As noted by Seigel [85], this is theeffective conductivity observed in a DC survey of a material with conductivity σ∞in the absence of IP effects and chargeability η . Next I write the total earth currentj(t) asj(t) = σ0e(t)− r(t). (3.6)I call r(t) the residual current. It contains any deviations from standard non-dispersive current flow. This decomposition of the current will be important forachieving a stable and robust discretization of the time-domain dispersive Maxwellequations.Before specializing to the stretched exponential case, consider defining the con-39stitutive law relating e and r by the following ordinary differential equation (ODE)in r:ασ0∂e∂ t=∂r∂ t+λ (α,θ ,β ,x, t)r, (3.7)where λ is an arbitrary smooth function of position and possibly time, and α ,θ , and β are spatially varying parameters. They will be related to the Cole-Coleparameters below. I will also associate the initial condition r(0) = 0 with the ODE.This condition gives the correct initial current for both grounded and inductivesource time-domain EM experiments.The ODE (3.7) has the solutionr(t) = ασ0∫ ∞0G(λ , t− τ˜)∂e∂ τ˜dτ˜ (3.8)where G is a Green’s function. The r(0) = 0 conditions implies G = 0 for τ˜ > t,which is consistent with the physical notion that r is caused by an applied electricfield e. Thus the integration may be truncated at t, givingr(t) = ασ0∫ t0G(λ , t− τ˜)∂e∂ τ˜dτ˜. (3.9)Considering the derivative of the Green’s function in the distributional sense,one may use a generalized integration by parts to giver(t) =−ασ0∫ t0∂G(λ ,x, t− τ˜)∂ τ˜edτ˜. (3.10)Comparing this result with eq. (3.4) it is seen that −ασ0 ∂G(x,t−τ˜)∂ τ˜ is equivalent toa conductivity impulse response that arises when transforming a frequency depen-dent Ohm’s law to the time domain.In principle one could choose the function λ in eq. (3.7) to generate the Cole-Cole impulse response but except for the special cases of c = 1,0.5, this wouldlikely still lead to an integral that would need to be evaluated numerically. Stretchedexponential relaxation offers a similar range of behaviour to the Cole-Cole modeland has a simple functional form in the time-domain. Alternatively, I will choose λsuch that the residual currents decay with a stretched exponential time dependencein the absence of induction effects. One may ignore induction when ∂e∂ t is negligi-40bly small. Thus we choose lambda such that r=Ae−tβ /θ satisfies the homogeneousversion of eq. (3.7)∂r∂ t+λ (α,θ ,x, t)r = 0. (3.11)This implies λ = θ−1β tβ−1. In the general case when induction and IP effectsoccur simultaneously, r will depend on the electric field and the constitutive lawrelating current and electric field will be defined by eq. (3.6) and the inhomoge-neous ODEασ0∂e∂ t=∂r∂ t+θ−1β tβ−1r. (3.12)The forward modelling initial boundary value problem is then fully specified byMaxwell’s equations ((3.1)), the constitutive relation defined by eqs. (3.6) and (3.12),and the initial condition r(0) = 0.3.1.1 Comparison of stretched exponential and Cole-Cole modelsIt is instructive to compare the properties of the stretched exponential to the Cole-Cole model. In the case of Debye dispersion the two models are entirely equiv-alent, which may be shown with the appropriate choice of parameters. Follow-ing Marchant et al. [63] and Haber [37], I write Ohm’s law in frequency for theDebye model as1σ0(1− iωτη1+ iωτ)J = E. (3.13)In this case, Ohm’s law may be analytically transformed to the time domain, whereit becomes the differential equation(1−η)∂ j∂ t+ τ−1j = σ0(∂e∂ t+ τ−1e), (3.14)This is equivalent to eq. (3.12) if we substitute j = σ0e− r, α = 1− (1−η)−1,θ = τ(1−η), and set β = 1. This shows that in the stretched exponential model,α plays the role of chargeability and θ that of a time constant. The parameter βplays a role analogous to that of c in the Cole-Cole model.For c,β 6= 1 the models are not fully equivalent but they share many char-acteristics. These connections can be explored by comparing the Cole-Cole and4110−510−310−1101103time (s)10−1010−810−610−410−2100102impulse response (Sms)c, β=0.2c, β=0.5c, β=0.8c, β=1.0(a)10−510−310−1101103time (s)10−1010−810−610−410−2100102impulse response (Sms)(b)Figure 3.1: a) Stretched exponential impulse response for various values ofβ . b) Cole-Cole impulse responses for various values of c. For bothplots η = 0.1, τ = 0.1 and σ∞ = 0.1.stretched exponential impulse responses. The stretched exponential impulse re-sponse isσˆse(t) =−ασ0 ∂G∂ t =−ασ0θ−1β tβ−1e−(tβ )/θ . (3.15)The α and σ0 parameters control the amplitude of the response, τ controls the over-all time scale of the response, and β the shape. The impulse response is very smallfor times much greater than τ . Recall that the residual current r is the convolu-tion of the impulse response with the electric field. This means that when gainingintuition about r one should think of early times in the impulse response plots asaffecting the behaviour of r at late times. For example, a sharp drop in σˆse(t) atlate times corresponds to a negligible IP current at early times after source shutoffand σˆse(t) increasing at early times as β decreases corresponds to slower decayof r with decreasing β . This is illustrated in fig. 3.1a, which shows the stretchedexponential impulse response σˆse(t) when η = 0.1, τ = 0.1 and σ0 = 0.1(1−η)plotted for several values of β . The Cole-Cole impulse responses for the same pa-rameter values are shown in fig. 3.1b. The responses are very similar at early timesbut the stretched exponential curves decay much more quickly at late times thantheir Cole-Cole equivalents.The Cole-Cole impulse responses were computed based on analytic series rep-resentations derived by Hilfer [44]. To my knowledge his paper is the only source42to give a general analytic expression for the Cole-Cole impulse response. He de-rives frequency and time domain representations of several models of relaxation indisordered systems (including the stretched exponential and Cole-Cole models) interms of contour integrals of a specific form known as H-functions. Infinite seriesrepresentations of these relaxation functions may be derived from the H-functionrepresentations. The series form of the Cole-Cole impulse response isddtEc [−(t/τ)c] , (3.16)where Ec is the Mittag-Leffler functionEc(x) =∞∑k=0xkΓ(ck+1). (3.17)While this series representation is exact, it is very ill-behaved in certain parameterregimes and generally difficult to work with numerically. To my knowledge, morerobust analytic formulae for the Cole-Cole conductivity impulse response only ex-ist for the special cases c = 1 and c = 0.5.Although the general form of the Cole-Cole impulse response is rather un-wieldy, Marchant [62] was able to show that for t τ and any c ∈ (0,1), it has theapproximate formσˆcc ≈ mtc−1+d (3.18)for some constants m and d. Inspection of eq. (3.15) shows thatσˆse(t)≈−ασ0θ−1β tβ−1 (3.19)for t τ , the same time dependence as the Cole-Cole response. This implies that atlate times the Cole-Cole and stretched exponential models generate approximatelyequivalent current flow, again recalling that early time behaviour of the impulseresponse corresponds to late time current flow.The relationship between the Cole-Cole and stretched exponential impulse re-sponses can be further illustrated through the exercise of trying to fit the stretchedexponential impulse response function to a given Cole-Cole response. I used non-linear least squares to find the stretched exponential parameters that best fit the4310−410−310−210−1100101102time (s)10−1010−810−610−410−2100102impulse response (Sms)Figure 3.2: Cole-Cole impulse response with c= 0.5, τ = 0.1, η = 0.1, σ∞ =0.1 plotted in blue alongside stretched exponential impulse responseswith parameters chosen by non-linear least squares fitting. The greenline attempted to fit the 0.001 < t < 100s time range and the magentaline the 0.1 < t < 100s range.Cole-Cole impulse response with parameters c = 0.5, σ∞ = 0.1 S/m, η = 0.1, andτ = 0.1 s. Figure 3.2 shows the Cole-Cole response along with two fitted stretchedexponential models, one to the time range 0.001 < t < 100 s and one to the range0.1 < t < 100 s. As expected, when attempting to fit the wide time range it is easyto fit the early-time Cole-Cole response at the expense of capturing its late timepersistence. The best fitting stretched exponential model had α = 0.088, β = 0.47,θ = 0.23. Using 1−α = 1/(1−η) and θ = τ(1−η), α corresponds to η = 0.081and θ corresponds to τ = 0.25. Again, recall that this excellent fit at early timesin the impulse response corresponds to equivalent residual current decays at latetimes.It is not possible to fit the entire Cole-Cole model with a single stretched ex-ponential. Trying to fit only the late time Cole-Cole impulse response gives avery rough approximation, as shown in fig. 3.2. This analysis also shows thatthe stretched exponential parameters do not correspond exactly to Cole-Cole pa-rameters, even when the responses coincide. However, if one is interested in thetime-scale where the models are equivalent it would be straightforward to find the44Cole-Cole parameters that match a certain stretched exponential electric or mag-netic field decay. If one is concerned with matching the Cole-Cole behaviour veryclosely over a wide range of time-scales then a residual current model based on asum of stretched exponential decays with different time constants could be used.I have not explored that possibility in this thesis but it would be a straightforwardextension.3.2 DiscretizationI present two discretizations of the coupled system of differential equations com-posed of eqs. (3.1a), (3.1b) and (3.12). In both cases the spatial operators arediscretized using the mimetic finite volume method on locally refined rectilinearmeshes known as OcTrees, as discussed in chapter 2. In previous work of my Ph.D.supervisor’s research group (e.g. [41]), which this thesis builds upon, the back-ward Euler method has been used for discretization in time of the non-dispersivequasi-static Maxwell equations. Backward Euler was chosen for its simplicity andexcellent stability properties. The quasi-static Maxwell equations in time are verystiff and implicit time-stepping was used to avoid having to take excessively smallsteps to preserve stability.The backward Euler discretization described here extends that approach to in-clude the stretched exponential auxiliary ODE formulation of Ohm’s law. A similaralgorithm for the special case of Debye dispersion and for the Cole-Cole model de-scribed by Pade´ approximation was described by [63]. I will show how the systemsof equations arising from this discretization can be solved efficiently using a combi-nation of direct and iterative methods. In the second approach time-discretization isaccomplished by the second order backward differentiation formula (BDF2). Thismore accurate time-stepping method may be used to parallelize the time-steppingprocess by the method developed in chapter 2.3.2.1 Backward Euler AlgorithmI begin the first scheme by discretizing the time derivatives in Maxwell’s equationsand Ohm’s law. I will then manipulate this semi-discrete system of equations toeliminate explicit dependence on the magnetic flux density b and form an equation45for the electric field e at time t in terms of quantities at earlier times. This equationwill then be discretized in space to form a system of linear-algebraic equations tobe solved for e at each time-step.Discretizing eqs. (3.1a), (3.1b) and (3.12) in time using the backward Eulermethod but keeping spatial variables continuous yields a system of semi-discreteequationsbn+1−bn∆tn=−∇× en+1 (3.20a)∇×µ−1bn+1−σ0en+1 = jsn+1− rn+1 (3.20b)rn+1− rn+∆tntβ−1n+1 θ−1rn+1 = ασ0(en+1− en), (3.20c)where ∆tn is the size of the nth time-step and the n+1 subscript on field variablesrefers to the field at time tn+1 =∑nj=1∆t j. These three equations can be combined toform a single equation expressing en+1 in terms of the residual current and electricfield at earlier times, and the source current:∇×µ−1∇× en+1+∆t−1n σ0(1− γα)en+1 =−∆t−1n[jsn+1− jsn+(1− γ)rn− (1− γα)σ0en], (3.21)where I define γ = (1+∆tntβ−1θ−1)−1 for brevity.Next, following the methods described in section 2.1.1, I convert eq. (3.21) toweak form and discretize it. This results in a system of linear algebraic equationsthat can be solved for the approximate electric field discretized on cell edges(CT Mµ−1C+∆t−1n Mσ0(1−γα))en+1 =−∆t−1n[sn+1− sn+M(1−γ)rn−Mσ0(1−γα)en]. (3.22)The matrix C is the discrete curl operator. All M matrices are mass matrices. Theirsubscripts indicate the PDE coefficients they discretize. The discrete source currents is computed by discretizing the transmitter wire path onto mesh edges.After solving for the electric field at a given time-step the residual current mustalso be updated, for use in the right hand side of the next electric field update.46To derive the residual current update I solve eq. (3.20c) for rn+1 and discretize inspace, giving the linear systemMrn+1 = Mγrn+Mσ0γα(en+1− en), (3.23)where M is a unit coefficient mass matrix. M is not diagonal in general but it isvery sparse and has a structure amenable to efficient factorization by sparse directmethods. It may be factored once at the beginning of a simulation and eq. (3.23)can be solved quickly by triangular substitution thereafter.Initial conditions are the last ingredient in the discretization. I assume steadystate initial conditions, which have the convenient implication that the initial resid-ual current is always zero. The initial electric field may either be zero or given bythe solution to a DC resistivity problem, as occurs e.g. when modelling step-offgalvanic sources. In this case I use a mimetic finite volume method to solve for theDC potential on mesh nodes. The initial electric field is then given by the nega-tive discrete gradient of the potential. In this chapter I present numerical examplesonly for step off sources but the implementation is capable of handling arbitrarytransmitter waveforms.3.2.2 Solving the linear systemI now move to solving the electric field-update system eq. (3.22). For realistic pa-rameter values and step-sizes the stretched exponential time stepping system canbe considered a perturbation from the corresponding system arising from back-ward Euler discretization of the non-dispersive Maxwell equations. As such, ex-perience with the non-dispersive equations has guided the approach here and Iwill frequently refer to the non-dispersive problem in the following discussion. Innon-dispersive TEM simulation and inversion algorithms based on implicit time-stepping it is often preferable to solve the electric field update system using directfactorization methods rather than iterative techniques such as the conjugate gradi-ent method. Stretched exponential modelling creates additional challenges but Ican still make efficient use of direct methods and argue that they are preferable toiterative methods.Over the last fifteen years, several excellent sparse matrix factorization soft-47ware packages have been developed [26] that make direct methods practical forvery large matrices such as those encountered in 3D electromagnetic simulations.In this work I have chosen to use the MUMPS [3] software package because ofits performance and numerical stability. Matrix factorization is expensive both interms of memory and CPU time but once it has been completed, the factors canbe reused to solve systems with the same matrix and different right hand sides ex-tremely efficiently. In backward Euler discretization of the non-dispersive Maxwellequations (and in the dispersive TEM algorithm of [62]), the electric field updatesystem depends on the step-size but not on the current time. Therefore, as longas the step size remains constant, one may factor the system matrix once and thenadvance the solution many time steps for relatively little computational cost. Inan inversion using gradient based optimization, the factors may also be reused toaccelerate the computation of cost function gradients.In transient electromagnetic simulation, one is generally interested in a widerange of time-scales. Airborne systems typically collect data over several orders ofmagnitude in time. To maximize information about both inductive and IP effectsin a grounded source survey one may need to simulate electric fields over four orfive decades in time. In such scenarios, very small step-sizes must be used to accu-rately model early time behaviour and it is not feasible to use a single step-size forthe entire simulation. However, using a small number of step-sizes, factoring thesystem matrix each time the step-size changes and advancing the solution severalsteps for each step size has still been more efficient than using iterative methods.This strategy was analyzed in detail by Oldenburg et al. [68].Unfortunately, time appears explicitly in the stretched exponential electric fieldtime-stepping system matrix in eq. (3.22), causing the matrix to change every stepand destroying the main benefit of direct methods. Requiring a factorization atevery time-step is too slow to be of practical use. For example, in one sampleproblem (the gradient array example discussed later in this paper with an additionalsource added) with two long grounded wire sources on a mesh of approximately114000 cells with 80 time-steps and 6 unique step sizes, non-dispersive modelling,using 6 factorizations, took an average of 285 s on a high end laptop. Stretchedexponential modelling took 2415 s factoring at each time-step.Since it is now not possible to achieve the high reuse of matrix factors achieved48for non-dispersive simulations, I consider iterative Krylov subspace methods forsolving the time-stepping system of equations. These methods do not manipulatethe entries of the system matrix directly. They create sequences of vectors thathopefully converge to the solution of the system, requiring matrix-vector productsat each iteration but treating the matrix as a black box.The matrix in eq. (3.22) is symmetric positive-definite and so I turn to thepreconditioned conjugate gradient method, the most effective iterative method forsuch systems [81]. Convergence of the conjugate gradient algorithm depends onthe distribution of eigenvalues of the system matrix [81]. The algorithm will con-verge quickly if the eigenvalues are tightly clustered away from zero. Unfor-tunately the stretched exponential time-stepping matrix is highly ill-conditioned,making it difficult to use the preconditioned conjugate gradient algorithm withstandard preconditioners to solve the system. The ill-conditioning is a consequenceof the non-trivial the null-space of the curl operator and the fact that as a discretiza-tion of a differential operator the largest eigenvalue remains relatively constant,while the smallest eigenvalue shrinks with refinement of the mesh. Adding themass term eliminates the zero eigenvalues associated with the null space of the curlbut its values are generally not large enough to make the system well-conditioned.Thankfully, the explicit time dependence creates only a small perturbation in thesystem matrix away from the non-dispersive case. Furthermore, as the change inthe mass term will be relatively uniform across the mesh, it seems likely that thechange in the mass term at each time-step will cause a shift in the eigenvaluesrather than increasing or decreasing their spread.Motivated by the facts of the previous paragraph, I have developed a precon-ditioner for the stretched exponential time-stepping system matrix based on occa-sional matrix factorizations. Let us denote the system matrix at time t0 as A0 andthe matrix at time tn = t0 + n∆t by An. I emphasize again that these matrices de-pend on both the step size and the current time. While the step size ∆t remainsconstant, the matrix A0 is a good approximation of the matrices An at subsequenttimes t0+n∆t. I can writeAn = A0+δAn. (3.24)Now, consider solving systems involving An with the preconditioned conjugate49gradient method. If δAn is not too large and A−10 can be cheaply applied then A0should be an excellent preconditioner for the conjugate gradient method appliedto systems involving An. If A0 is factored and the factors are stored then forwardand back substitution can be used to cheaply and accurately apply A−10 to arbitraryvectors, as is required when using A0 as a preconditioner.This motivates the use of a hybrid direct-iterative approach to updating thestretched exponential electric field. At the first time-step, and subsequently whenthe step-size changes, the system matrix is factored. Otherwise the system is solvedusing the block preconditioned conjugate gradient method (PCG) [69] with thefactored matrix as a preconditioner. I use the block variant rather than standardPCG to improve convergence in systems with multiple right hand sides. This givesthe forward modelling algorithm shown in algorithm 1.Algorithm 1 Electric field update algorithmInput:Number of times steps nList of step sizes ∆tInitial electric fieldDo:Factor initial matrix A0 = LLTSolve for e by forward and back substitutionfor j = 2→ n doConstruct matrix A jif ∆t(n) = ∆t(n−1) thenSolve A jX = B with block PCGelseFactor A jSolve using forward and back substitutionend ifend forThis approach yielded a significant decrease in forward modelling run-times. Iobserved a small degradation in accuracy, which was negligible except after a largenumber of time steps. I tested the algorithm by comparing run-times and computedelectric fields when factoring at each time-step and when using the hybrid direct-5010−210−1100t (s)10−710−510−310−1101relative error normFigure 3.3: Error in PCG computed example problem electric fields for looseand tight PCG relative residual stopping tolerances.iterative algorithm. I considered the fields computed by factoring at each step(edir) to be correct and computed the error in the fields computed with the hybridalgorithm (eiter) asδe =‖edir− eiter‖2‖edir‖2 . (3.25)The highly ill-conditioned nature of the system makes it difficult to compute theelectric field accurately, with very small PCG residuals not necessarily correspond-ing to accurate electric fields. These errors accumulate over the course of the time-stepping process since fields and residual currents computed at a given time stepare used in the right hand side of future update linear systems. The relative errornorm of the electric field for PCG residual stopping tolerances of 1× 10−6 and1×10−14 is shown in fig. 3.3. The fields were computed from the example modelmentioned above, with 80 times steps, covering a time range 0.0025-3.1 s, using sixunique step sizes. Hybrid forward modelling on that model using a PCG relativeresidual stopping tolerance of 1×10−14 took an average of 702 s, 2.5 times slowerthan non-dispersive modelling of the same example but only 29% of run-time whenfactoring at every step. The number of PCG iterations per time step varied betweentwo and four. The time increase relative to non-dispersive modelling is mainly dueto the PCG iterations requiring multiple applications of forward and back substitu-51tion during application of the preconditioner but also due to the cost of assemblingmass matrices at each iteration and computation of the residual current, both ofwhich are not required in non-dispersive modelling.In my algorithm, A0 is a simple Neumann series preconditioner. For any squarematrix A for which the sequence {Ak}∞k=0 converges to zero, the following identityholds(I−A)−1 =∞∑k=0Ak (3.26)and the series is called the Neumann series of A. Defining A=A0+δA, it followseasily from the Neumann series identity thatA−1 =∞∑k=0(−A−10 δA)kA−10 (3.27)as long as {(A0δA)k}∞k=0 converges to zero. The eigenvalues of A−10 are boundedso convergence will occur if δA is small in a suitable norm. Using a truncationor other approximation to this series as a preconditioner is called Neumann seriespreconditioning [81]. My algorithm uses zeroeth order Neumann preconditioning.If δA is too large then the Neumann series may converge slowly or not at all.In the sample problems we tested, zeroeth order Neumann series preconditioningworked extremely well when refactoring the system matrix each time the step sizewas changed. I believe that the method would continue to work for all time-scales,models and meshes of interest to applied geophysicists. However, in regimes withlarger perturbations, a higher order Neumann series preconditioner could easily beimplemented to accelerate convergence.3.2.3 BDF2As in the non-dispersive case of chapter 2, higher order time-stepping is desirable,both to improve the accuracy of sequential time-stepping and to allow the use ofthe parallel time-stepping algorithm of section 2.2. The second order time-steppingscheme for stretched exponential forward modelling is based on the second orderbackward derivative approximation. We follow the same sequence of steps as inthe backward Euler discretization, discretizing the time derivatives in eqs. (3.1a),52(3.1b) and (3.12) with BDF2 approximations. Assuming a constant step-size, thisgives the semi-discrete system32bn+1−2bn+ 12bn−1 =−∆t∇× en+1 (3.28a)∇×µ−1bn+1−σ0en+1 = jsn+1− rn+1 (3.28b)32rn+1−2rn+ 12rn−1+ tβ−1n+1 ∆tθ−1rn+1 =ασ0(32en+1−2en+ 12en−1).(3.28c)As with the backward Euler scheme, eliminating explicit dependence on b anddiscretizing in space leads to a system of linear algebraic equations for en+1 at eachtime step in terms of fields at earlier times and the source current(CT M fµ−1C+32∆t−1Meσ0(1−γα))en+1 =32∆t−1[sn+1− 43sn+13sn−1+Meσ0(1−γα)(43en− 13en−1)−M1−γ(43rn−13rn−1)],(3.29)where ∆t is now a constant step-size and γ2 = (1+(2/3)tβ−1n+1 ∆tθ−1)−1.3.2.4 ImplementationThese algorithms were implemented in the Julia programming language. My codewas designed to be compatible with jInv [80], a Julia software framework for partialdifferential equation (PDE) constrained parameter estimation problems, which wasbriefly described in section 2.1.1. In the context of this chapter, in which I am onlydiscussing forward modelling, this is a minor point but it will make it simple toincorporate this simulation code into an inversion software package that will allowone to invert for intrinsic IP parameters.I verified the correctness of the code using the method of manufactured solu-tions [82]. In this technique, one defines an electric field that meets the boundaryconditions of the problem, then substitute it into Maxwell’s equations and analyti-cally calculates the magnetic flux density, residual current and source current that53match the manufactured electric field such that Maxwell’s equations are satisfied.Then these exact formulae are used to compute discretized initial electric field,residual current and source currents. These values may be fed into the software tocompute the approximate electric field and residual currents at later times. The rateof convergence of solutions in space and time can then be assessed by varying themesh cell and time-step sizes.3.3 Synthetic modelling examples3.3.1 Grounded source exampleI will now illustrate the capabilities of stretched exponential modelling of coupledelectromagnetic induction and induced polarization (EMIP) effects using two syn-thetic examples where both EM induction and IP effects are significant—a gradientarray survey and a coincident loop airborne time-domain electromagnetic (ATEM)sounding. The gradient array survey earth model was a homogeneous halfspacecontaining two anomalous bodies, one chargeable and the other conductive. TheATEM sounding model consisted of a buried conductive block in a resistive half-space with chargeable overburden.First I consider the gradient array example. Gradient arrays allow for efficientcollection of easily interpretable direct current induced polarization (DCIP) dataover large areas but they have significant deficiencies when viewed from the DCIPperspective. The major deficiencies are a lack of depth resolution and significantEM coupling. The EMIP approach can significantly address these deficiencies.The survey layout and model for our example is shown in fig. 3.4. A 2 km longstraight grounded wire is laid out over the survey area and an array of electric dipolereceivers parallel to the transmitter dipole are laid out in a grid covering a 500 ×500 m rectangular area centred at the centre of the transmitter wire. The transmittercurrent waveform is modelled as a step-off. In the transmitter on-time, galvanic DCcurrents allow polarization to build up in the presence of chargeable bodies. Whenthe transmitter current is shut off, the relaxation of these polarizations, as wellas electromagnetic induction, will generate transient electric currents in the earth.It can be useful to think of the EM induction aspect of an EMIP gradient array54TxA BRx(a)−1000 −500 0 500 1000x(m)−500−2500250500y(m)(b)Figure 3.4: Gradient array example survey layout. a) Cartoon showing trans-mitter, receiver array bounding box, and sample cartoon receiver. b)Plan view layout of transmitter, receiver array and two block syntheticmodel.survey as being analogous to a time-domain CSAMT survey for long groundedwire sources but near field effects can be important so the analogy is not exact.Away from the transmitter electrodes, in our area of interest near the centre of thewire, the primary magnetic field will be approximately solenoidal in the transmitteron-time. Transmitter shutoff will induce currents to flow parallel to the directionof the wire path.The plan-view locations of the two anomalous blocks in the earth model areshown, along with the transmitter wire path, in fig. 3.4b. Both blocks are 150 ×150 × 75 m rectangular prisms buried 50 m below the earth’s surface. The lefthand block had the same conductivity as the background halfspace (σ∞ = 0.01),chargeability η = 0.25, IP time constant τ = 0.5 s and stretched exponential decayexponent β = 0.5. Recall that the stretched exponential parameter α is relatedto η by the relation 1−α = 1/(1−η) and that the parameter θ is related to τby θ = τ(1− η). Recall also that chargeability is equivalent in the Cole-Coleand stretched exponential models. The time constants and exponents play roughlyequivalent roles in both models but of course because the models are not fullyequivalent we cannot in general choose a set of stretched exponential parametersthat will lead to an electric field decay exactly equal to some Cole-Cole decay. Withthose reminders of the roles of the parameters I move to the right hand block whichwas one order of magnitude more conductive than the background (σ∞ = 0.1) andwas not chargeable. Air conductivity was set to 1×10−8 S/m.55Figure 3.5: Visualization of OcTree mesh discretization of the conductivitymodel.The model was discretized on an OcTree mesh of 113772 cells, yielding alinear system of 319653 unknowns to be solved to compute the electric field onmesh edges at each time-step. Cubic cells with side length 25 m were used in thecore region of the mesh, including the area of interest and the transmitter wire path,with cell size expanding away from the area of interest and transmitter. A verticalcross section of the mesh showing cell size and conductivity in the y = 0 planeis pictured in fig. 3.5. Starting from a DC initial electric field, the simulation ranto a final time of 3.11 s after current shutoff, using 80 time-steps with 6 uniquestep sizes ranging from 2.5×10−3 s to 0.1 s. Using Backward Euler time-stepping,solving linear systems using the iterative hybrid direct algorithm described above,this simulation had an average run-time of 620 s, including the time to solve theDC problem for the initial electric field.The inline component of the DC electric field is shown in fig. 3.6. The plotshows a clear anomaly from the conductive block but only a small distortion ofthe field contours due to the chargeable block. The maximum discrepancy in theelectric field between this model and an otherwise equal model with the chargeableblock removed is approximately 5%. The difficulty of distinguishing between theIP and non-dispersive induction responses in the off-time electric field data depends56−250 0 250x(m)−400−2000200400y(m)1.7e-052.5e-053.2e-054.0e-054.7e-055.5e-056.3e-057.0e-05ex(V/m)Figure 3.6: DC inline electric field at the earth’s surface for two block ex-ample model. The locations of the buried blocks are shown by blacksquares.on several factors including the depth to conductive targets, depth to chargeable tar-gets, and the time scale of IP responses. In scenarios with very good conductorsand small IP time constants, traditional EM decoupling procedures may performpoorly. The ability to fully model the coupled EMIP response in 3D should allowfor better isolation of the IP response. More generally, having such modelling ca-pability allows inductively generated electric fields to be used as signal in groundedsource IP, not just noise. Early-time transient electric fields provide additional in-formation on ground conductivity, which can then be used to improve recovery ofIP parameters [52].The distinct characters of the decays of the electric fields near the conductiveand chargeable blocks in this example model are illustrated in fig. 3.7. It is difficultto distinguish the IP response of the chargeable block at early times when inductiondominates. The stretched exponential and non-dispersive fields are nearly identi-cal before 1× 10−2 s. After that the IP response appears clearly in the surfaceelectric field near the chargeable block, as shown in figs. 3.7b and 3.7c. By con-trast, fig. 3.7d shows that the stretched exponential surface electric field over theconductive block is almost identical to the non-dispersive field, as expected.57−250 0 250x(m)−400−2000200400y(m)-7.5e-06-7.0e-06-6.5e-06-6.0e-06-5.5e-06-5.0e-06-4.5e-06-4.0e-06-3.5e-06-3.0e-06ex(V/m)(a)−250 0 250x(m)−400−2000200400y(m)-1.5e-07-1.2e-07-9.0e-08-6.0e-08-3.0e-080.0e+003.0e-086.0e-089.0e-08ex(V/m)(b)10−210−1100t(s)10−910−710−5ex(V/m)(c)10−210−1100t(s)10−910−710−5ex(V/m)(d)Figure 3.7: In-line component of surface electric field for gradient array ex-ample. a) Plan-view in-line field at 0.01 s after current shutoff. Left andright black squares show the boundaries of the chargeable and conduc-tive blocks, respectively. b) Plan view field at 0.14 s. c) Electric fielddecays above the centre of the chargeable block, with stretched expo-nential field in blue and corresponding non-dispersive field (η = 0) inred. Dashed lines indicate negative values. d) Electric fields above thecentre of the conductive block, showing that the stretched exponentialand non-dispersive responses are almost identical.3.3.2 Inductive source exampleThe following example illustrates the applicability of stretched exponential mod-elling to the simulation of airborne inductive source induced polarization. The phe-nomenon of negative transients in coincident and concentric loop EM systems hasbeen known to exist for a long time. Weidelt [99] showed that these sign reversalscould not occur with non-dispersive electrical conductivity and magnetic perme-ability, suggesting that IP effects could be responsible. There has been renewedinterest in the phenomenon in recent years, as improvements in data acquisition58−100 0 100x(m)−100−500Depth (m)σ∞=0.01, η=0.3, τ=1×10−3σ∞=0.005σ∞=0.1Figure 3.8: Section view of block in a half-space with chargeable overbur-den model. The overburden depth is 10 m and the block extends from40-90 m depth with a plan view cross-section of 100×100 m. The back-ground halfspace conductivity was 0.005 S/m and the conductivity ofthe block was 0.1 S/m. The overburden had σ∞ = 0.01 S/m η = 0.3,τ = 1×10−3 s and β = 0.8have led to these signals appearing with increasing frequency in field data—seee.g. [61]. I illustrate how negative transients can be modelled with stretched ex-ponential polarization with a simple model of a single sounding of a coincidenthorizontal loop ATEM transmitter with step off source over a resistive half-spacecontaining a conductive rectangular prism buried 40 m beneath the earth’s surfaceand 10 m of chargeable overburden. The earth model is illustrated and the physicalproperty values are given in fig. 3.8. The transmitter was a square loop with 10 mside length located 30 m above the surface, centred over the conductive block. Thischoice of IP parameters is representative of the type of fine-grained near surfacematerial thought to be most commonly responsible for IP signals in ATEM data[60]. In addition to smaller time-constants, these materials tend to have sharperdecays than typical ground IP targets, leading to β values closer to 1.The effect of the chargeable overburden on the ATEM response of our examplemodel is illustrated in fig. 3.9. It shows the total EMIP responses alongside thecorresponding responses computed when neglecting the chargeable overburden.We used the stretched exponential algorithm with backward Euler time stepping tocompute the transient electric field over the time interval [1× 10−5,0.002] s. Thevertical component of the magnetic flux density (dbz/dt) was then computed from5910−410−3t(s)10−2100102104−dbz/dt(nT/m)(a)10−410−3t(s)10−2100102104−dbz/dt(nT/m)(b)10−410−3t(s)10−2100102104−dbz/dt(nT/m)(c)10−410−3t(s)10−2100102104−dbz/dt(nT/m)(d)Figure 3.9: Effect of chargeable overburden on ATEM response from buriedconductor. All plots show the vertical component of db/dt. The back-ground halfspace response is shown in green in each plot. Dashed linesindicate negative values. a) Block with no overburden. The responsewith the conductive block at a depth of 40 m is shown in blue. Themagenta line shows the response obtained when the top of the blockis moved to a depth of 70 m. b) Stretched exponential response with10 m of overburden (σ∞ = 0.005, η = 0.3, τ = 1× 10−3, β = 0.8) isshown in black, along with the overburden free and background re-sponses. c) Stretched exponential response with the top of the blockat 70 m in black, along with the overburden free and background re-sponses. d) Stretched exponential response with chargeable overburdenand no conductive block in black, along with background EM response.60the electric field values. When the chargeable overburden is neglected (fig. 3.9a)the response from the block at 40 m depth is clearly distinguishable from the back-ground response. The signal from the block remains but is much smaller when itis lowered to 70 m depth. The EM induction signal from the shallower block isstrong enough and emerges early enough that both it and the IP response of theoverburden are distinguishable in the EMIP data—fig. 3.9b. The EMIP responsenever changes sign due to the cancellation of the EM induction and IP signals butthe decay is clearly indicative of an IP response. It is important to note here that amore subtle IP response could be difficult to distinguish from a resistive body with3D geometry. I believe this is an important argument in favour of the use of 3Dmethods when considering the presence of IP effects in ATEM data.When the top of the block is moved to a depth of 70 m (fig. 3.9c), the responseof the conductor is completely masked by the IP response. However, as fig. 3.9dshows, the IP response generated by the overburden is much larger and is mani-fested in the data earlier in time without the presence of the block. Modelling thefull physics of an ATEM system with IP effects included allows for quantitativesimulations of the interactions of non-dispersive EM and IP responses. The abilityto model these effects in 3D is very useful in its own right and the key step towarddeveloping a 3D inversion algorithm for conductivity and IP parameters.3.4 ConclusionI have developed a novel computationally efficient forward modelling algorithmfor 3D time-domain electromagnetic modelling in the presence of chargeable ma-terials, based on stretched exponential relaxation. The Cole-Cole model is still byfar the most common framework for modelling and classifying chargeable materi-als. Stretched exponential relaxation is not equivalent to Cole-Cole relaxation butit is capable of modelling a similar range of phenomena. It is a significant com-fort to have shown that the two models do produce equivalent results at late times,showing that they can definitely model the same phenomena if the appropriatetime interval is considered. This work also identified the Cole-Cole and stretchedexponential responses within a more general context. When electric current maybe represented as a convolution with the electric field, the two approaches simply61imply different convolution kernels.I applied stretched exponential modelling to two synthetic examples possess-ing the key characteristics of two applied problems where both EM induction andIP effects are important. In the case of grounded source IP surveys with significantEM coupling, I believe that the ability to quantitatively simulate the transient EMand IP behaviour of the surveys opens up exciting possibilities for significantlyincreasing the amount of information about subsurface geology that they provide.In the case of ATEM surveying, the importance of treating IP effects in the anal-ysis of field data is increasingly clear. I have developed an algorithm capable ofefficiently modelling these effects. This algorithm forms the foundation of thethree-dimensional (3D) inversion algorithm presented in the next chapter, whichis to my knowledge the first algorithm to directly invert transient geophysical EMdata for intrinsic IP parameters.62Chapter 4Stretched exponential inversionThis chapter discusses my preliminary work into the inversion of transient elec-tromagnetic data for stretched exponential parameters. I have developed routinesfor computing the sensitivity of electromagnetic data to these parameters and per-formed simple proof of concept synthetic inversions. These preliminary results areencouraging but much more work needs to be done in order to understand how tobest extract induced polarization (IP) related information from transient electro-magnetic (EM) data in the stretched exponential approach.This is a multi-parameter inverse problem with a great deal of non-uniqueness.How to best handle that non-uniqueness is an open question. The best workflowis probably problem dependent. For example, should one invert for non-dispersiveconductivity and the other stretched exponential parameters simultaneously, or is itnecessary to have a good conductivity model before attempting to recover charge-ability and the other parameters? Is it always possible to get a good conductivitymodel before at all considering IP effects?In some problems it is possible a priori to separate a transient EM decay intothree sections: an early time range where EM induction dominates, a late timerange where IP effects dominate and a middle section where the two effects occuron similar scales and interact. In other problems this may not be possible andsimultaneous inversion for conductivity and IP parameters might be the only wayforward. In this chapter I focus on the former, simpler case. I demonstrate thecapabilities of stretched exponential inversion on the gradient array model problem63TxA BRx(a)−1000 −500 0 500 1000x(m)−500−2500250500y(m)(b)Figure 4.1: Gradient array example survey layout. a) Cartoon showing trans-mitter, receiver array bounding box, and sample cartoon receiver. b)Plan view layout of transmitter, receiver array and two block syntheticmodel.encountered in chapter 3.4.1 The gradient array exampleRecall the synthetic gradient array survey from section 3.3.1. Two perpendicularlong grounded wire sources are positioned such that they both pass through thecentre point of an array of 25 m electric dipole receivers oriented parallel to thetransmitters that measure the direct current (DC) and transient voltages in the earthdue to the transmitters. Thus, for the transmitter oriented parallel to the x coor-dinate axis only data from the x-oriented receivers were used. Similarly, only yreceivers were used for the y-axis oriented transmitter. The survey layout for thex-oriented receiver is shown in fig. 4.1. The model consists of two 150 × 150× 75 m blocks buried 50 m below the surface of a homogeneous half-space withconductivity σ∞ = 0.01 S/m. The true conductivity model is shown in fig. 4.2. Thewest block (left hand side of fig. 4.2a) has σ∞ = 0.018 S/m, chargeability η = 0.25,time constant τ = 0.5 s and exponent β = 0.5. The east block is conductive but notchargeable. It has σ∞ = 0.1 S/m.As in the forward modelling example in chapter 3, both transmitters used step-off waveforms to excite the earth. The initial DC voltages and transient voltages at85 times ranging from 5×10−4-1.128 s were computed. Standard deviations of 5%of the value plus a constant floor were assigned to each datum. The true noise-free64− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)10− 22 × 10− 23 × 10− 24 × 10− 26 × 10− 210− 1(S/m)(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.2: Selected slices of true conductivity model. Left hand block hadσ∞= 0.018 S/m and right block σ∞= 0.1 S/m. Left block had η = 0.25,β = 0.5 and τ = 0.5 s. Right block was not chargeable. a) z = −85 mb) y = 75 m c) y =−75 m.synthetic data were used in the inversions. I also committed the further “inversecrime” of using the same mesh and time-stepping scheme in the inversions that wasused to generate the synthetic data. As this work represents just an initial proof ofconcept for stretched exponential inversion, commission of the inverse crime wasdeemed acceptable.4.2 Inversion methodologyRecall that the stretched exponential forward modelling algorithm computes theelectric field on the cell edges of a computational mesh at a set of discrete times.The synthetic voltages are computed by approximating the line integrals of theelectric fields over the receiver wire paths. In this example the receiver wire pathswere always straight line segments coinciding with mesh edges.I seek to recover the non dispersive conductivity σ∞, the chargeability η , thetime constant τ and the exponent β from the DC and voltage data. The parametervalues are discretized onto the computational mesh, approximated as constant ineach cell. I seek to recover values of the parameters in all subsurface mesh cells. Iuse the following four stage inversion procedure:1. invert early time data, which is assumed to be negligibly effected by IP, forσ∞652. Hold σ∞ fixed and invert for all three IP parameters simultaneously3. Hold σ∞ and η fixed, invert for IP parameters τ and β simultaneously.4. Refine IP parameter estimates by inverting for all of them (η , τ , β ) simulta-neouslyStages 2-4 are initialized using the recovered models from the previous steps. Notethat DC data are included in the stage 1 inversion despite the fact that they areaffected by IP. However, even for DC, IP is a second order effect and it was judgedthat for this initial inversion it was acceptable to invert the DC data alongside theearly off-time data for σ∞, ignoring chargeability.At each stage the inverse problem is posed as a regularized least-squares pa-rameter estimation problem, using Tikhonov regularization [93]. I seek the modelm that minimizes the following objective function subject to bound constraints onthe parameter valuesminimizemΦ(m) =12‖W(F (m)−d)‖22+ γR(m)subject to ml < m < mh,(4.1)where d is the observed data, W is a diagonal matrix holding the standard devi-ations of the data, F is the forward modelling operator, R is the regularizationoperator, and γ is the scalar regularization parameter that balances the weights ofthe misfit and regularization terms in the objective function. For single physicalproperty inversions the regularization operator is the discretization of‖Wi∇m‖22 (4.2)where Wi is a weighting operator that allows smoothness to be enforced morestrongly in some regions of the model than others. The discrete regularizationoperator isR(m) = mT GT WiGm (4.3)where G is a cell-centered discretization of the gradient operator. It approximatesthe gradient of a cell-centre discretized scalar function on cell faces. Wi is an66interface weighting matrix that disretizes Wi on cell faces. For multi-parameterinversions the model is simply a concatenation of the models for each parameter.The regularization operator isRmulti(m) =∑jmT PTj GT Wi( j)GP jm, (4.4)where P j is a matrix that extracts the model for the jth parameter from the fullmodel and Wi( j) is the interface weighting matrix for the jth parameter.The optimization problem eq. (4.1) was solved using the Gauss-Newton algo-rithm with backtracking Armijo line search. I solved for the search direction usingthe projected preconditioned conjugate gradient algorithm with loose stopping tol-erance.There are several ways in which this procedure could be improved but they arebeyond the scope of this thesis. For example, structural constraints such as cross-gradient regularization (see e.g. [33]) could be imposed to penalize differences inspatial structure between the different physical property models. Or in cases wherethere is reason to believe that the IP signal in a transient EM dataset comes froma geometrically simple subsurface body, the techniques of chapter 5 of this thesiscould be used.4.3 Inversion resultsFirst I performed the conductivity only inversion, inverting the DC data and tran-sient data from 8 times ranging from 0.002-0.0055 s. Those times correspond tothe 5th to 12th steps of the forward modelling time-stepping. The first four stepscould not be used because they contained artefacts from the approximation of thestep-off waveform. In this synthetic example it was possible of course to knowexactly when IP effects start to become distinguishable in the data. In this exam-ple I tried to choose a cutoff time just before IP effects could be clearly seen in thedata without comparing it to the synthetic data generated from the true conductivitymodel without IP effects.Previous work by my supervisor and I [10] has shown that the use of early timetransient data can greatly improve the depth resolution of gradient array surveys.67However, in this work the focus is on studying IP and the interplay of EM andIP effects. In order to allow simulations to run to late enough times to see theIP effects overwhelm the EM effects without using a number of time steps thatwould render simulation too computationally demanding, I chose not to modelearly enough times to allow the conductivity inversion to recover the depths ofthe blocks correctly without the use of interface weighting. No sensitivity depthweighting was employed but the interface weights in the first subsurface layer ofthe mesh were set to 100 (they were 1 everywhere else) to penalize structures inthe model being pushed to the surface. This worked very well. The recoveredconductivity model is shown in fig. 4.3. The horizontal locations, depths, and− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)10− 22 × 10− 23 × 10− 24 × 10− 26 × 10− 210− 1(S/m)(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.3: Selected slices of recovered conductivity model from DC andearly off-time inversion. Thick black lines show true block boundaries.True conductivities are 0.01 S/m for the background, 0.018 S/m for theleft block and 0.1 S/m for the right block. a) z = −85 m b) y = 75 m c)y =−75 m.relative conductivities of both blocks were recovered well. This inversion was runto completion, meaning it was able to fit the data, according to the χ2 goodnessof fit criterion, stopping the inversion when the root mean squared (RMS) misfitdropped below 1. The misfit convergence curve is shown in figure fig. 4.4Next the inversion for η , τ , and β was performed with σ∞ held fixed to themodel recovered from the previous inversion. This inversion was haulted after 5Gauss Newton iterations. At that point the inversion began to struggle to decreasethe data misfit and inspection of the models at each iteration showed that the inver-sion had significantly improved the chargeability model but made only very small680 2 4 6 8 10iteration012345root mean squared errorFigure 4.4: Misfit convergence curve for σ∞ inversion.changes to the time constant and exponent models. I have not performed a sen-sitivity analysis but I believe the data are much more sensitive to changes in ηthan to the other IP parameters, especially at early iterations when, indeed, the τand β sensitivities depend on the chargeability model. The recovered chargeabilitymodel from this inversion is shown in fig. 4.5. The RMS misfit was approximately− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)0.000.050.100.150.200.25η(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.5: Selected slices of recovered chargeability model from DC andmid to late off-time data inversion. Thick black lines show true blockboundaries. Left hand block has η = 0.25. Right hand block is notchargeable. a) z =−85 m b) y = 75 m c) y =−75 m.3.8 when the inversion was stopped. So far the inversion has captured the location69of the chargeable block well but underestimated its chargeability and depth extent.Next I held both σ∞ and η fixed and inverted for the time parameters τ and β .This inversion was able to fit the data, in the χ2 goodness of fit sense, however, indoing so it created τ and β models with a great deal of spurious structure. The fi-nal recovered models from this stage 3 inversion are shown in figs. 4.6 and 4.7 the− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)0.50.60.70.80.91.0τ(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.6: Selected slices of recovered time constant model from late off-time data inversion. Thick black lines show true block boundaries. Lefthand block has true τ = 0.5. Right hand block is not chargeable. a)z =−85 m b) y = 75 m c) y =−75 m.− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)4 × 10− 16 × 10− 1100(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.7: Selected slices of recovered exponent model from late off-timedata inversion. Thick black lines show true block boundaries. Left handblock has true β = 0.5. Right hand block is not chargeable. a) z =−85 m b) y = 75 m c) y =−75 m.70spurious structure is clearly seen. Instead of using this as the final model, the mod-els from an earlier iteration of this inversion can be used to start another inversionfor all three IP parameters η , τ , and β . The early iterates from the stage 3 τ andβ inversion significantly improve on the results from the stage 2 inversion whileremaining relatively smooth. Using these early stage 3 iterates and the recoveredη model from the stage 2 inversion as starting models in a final inversion for allthree IP paramaters worked well. The recovered models from this final inversionare shown in figs. 4.8 to 4.10− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)0.000.050.100.150.200.25η(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.8: Selected slices of recovered chargeability model from stage 4 midto late off-time data inversion. Thick black lines show true block bound-aries. Left hand block has true η = 0.25. Right hand block is not charge-able. a) z =−85 m b) y = 75 m c) y =−75 m.The recovery of η is improved. The true depth extent is still not captured butthe planview geometry and chargeability of the recovered block are improved. Thetime constant model is also improved. There are many artefacts outside the charge-able region but this is not a concern. The recovered τ and β models should only beinterpreted in regions with non-neglible chargeability. These artefacts could likelybe significantly reduced by penalizing spatial inconsistencies between the η , τ , andβ models by a method such as cross-gradient. The β model contains artefacts onthe boundary of the chargeable region but has good recovery in the centre of theblock. I do not know the cause of these artefacts.Images of the true and predicted data from the stage 4 inversion for the x-directed receivers due to the x-directed transmitter for times 0.02 s and 0.17 s after71− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)0.50.60.70.80.91.0τ(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.9: Selected slices of recovered time constant model from stage 4mid to late off-time data inversion. Thick black lines show true blockboundaries. Left hand block has true τ = 0.5. Right hand block is notchargeable. a) z =−85 m b) y = 75 m c) y =−75 m.− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000100200y (m)5 × 10− 16 × 10− 17 × 10− 18 × 10− 19 × 10− 1(a)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(b)− 300 − 200 − 100 0 100 200 300 400x (m )− 200− 1000z (m)(c)Figure 4.10: Selected slices of recovered exponent model from stage 4 mid tolate off-time data inversion. Thick black lines show true block bound-aries. Left hand block has true β = 0.5. Right hand block is not charge-able. a) z =−85 m b) y = 75 m c) y =−75 m.transmitter shutoff are shown in figs. 4.11 and 4.12, respectively. The fit is quitegood at early times, as can be seen in fig. 4.11. The fit is not as tight at late times,with the predicted data showing an anomaly that is geometrically equivalent to butstronger than that seen in the true data. The observed and predicted voltage decaysat the centre of the chargeable block are shown in fig. 4.13.72−250 0 250x(m)−400−2000200400y(m)7.0e-078.3e-079.5e-071.1e-061.2e-061.3e-061.5e-061.6e-061.7e-061.8e-06ex(V/m)(a) true data−250 0 250x(m)−400−2000200400y(m)7.0e-078.3e-079.5e-071.1e-061.2e-061.3e-061.5e-061.6e-061.7e-061.8e-06ex(V/m)(b) Predicted dataFigure 4.11: Voltage data for the x-directed receivers due to the x-directedtransmitter at 0.02 s after source shutoff.−250 0 250x(m)−400−2000200400y(m)-3.6e-07-3.0e-07-2.4e-07-1.8e-07-1.2e-07-5.5e-086.3e-096.7e-081.3e-071.9e-07ex(V/m)(a) true data−250 0 250x(m)−400−2000200400y(m)-3.6e-07-3.0e-07-2.4e-07-1.8e-07-1.2e-07-5.5e-086.3e-096.7e-081.3e-071.9e-07ex(V/m)(b) Predicted dataFigure 4.12: Voltage data for the x-directed receivers due to the x-directedtransmitter at 0.02 s after source shutoff.4.4 ConclusionsMore work needs to be done in order to improve the stretched exponential inversionworkflow but this chapter gives a simple proof of concept example showing thepotential of the technique. The ease with which chargeability information wasrecovered from the transient EM data was encouraging, especially as this was thefirst work to recover the intrinsic chargeability directly from transient EM data inthree-dimensional (3D). More work needs to be done to understand whether it ispossible to achieve better recovery of the time parameters τ and β .A potential approach to reducing the non-uniqueness of the problem is to useparametric methods. The next chapter develops a stochastic parametric inversionalgorithm. Although IP is not addressed, the work of the chapter lays the ground-work for applying such an approach to stretched exponential inversion in the future.7310 1 100t (s)10 710 6ex(V/m)Figure 4.13: true and predicted voltage data above the centre of the charge-able block from 0.02 s after source shutoff onward plotted on log scale.Blue line shows observed data and red line shows predicted data.Dashed lines indicate negative values.74Chapter 5Randomized backgroundstochastic inversion5.1 IntroductionIn many cases, the goal of inverse problems is to estimate the physical propertiesof relatively simple objects that can be described by small numbers of parameters.Unfortunately, to solve such problems it is often necessary to also estimate a po-tentially large number of parameters that consist of the background to this model.Such background can be viewed as nuisance parameters. The quality of these nui-sance parameter estimates may then adversely affect the recovered estimates of theparameters of interest.In order to discretize the simple object and estimate its properties, in this paper,we take the view of parametric level-set methods. These methods fall under thegeneral setting of inverse shape reconstruction methods, which seek to estimatethe position and shape of an object embedded in some background medium fromindirect measurements such as geophysical or medical imaging data. Level setmethods have been a popular approach to shape reconstruction inverse problems.See [19] for a survey of early work on the subject and, e.g., [50, 73] for examplesof more recent work. In inverse level set methods, the boundary of the objectof interest is represented by the zero level-set of some function. The level-setfunction may be represented and manipulated directly, in which case it must be75represented cell-wise over the computational domain [19]. This is a quite flexibleapproach but creates a high-dimensional parameter space to optimize over. Thiscan lead to the inverse problem becoming highly ill-posed, requiring methods suchas regularization of the level-set function to address the ill-posedness [28].Alternatively, the level set function can be represented parametrically, usingsome sort of basis function expansion or geometric shape to represent the object ofinterest. For example, Aghasi et al. [1] used radial basis functions to reconstructgeneral 2D shapes. In our approach, we represent the object of interest as a skewedGaussian ellipsoid. This is based on the method of McMillan et al. [64], who usedGaussian ellipsoids to represent high contrast subsurface bodies in electromagneticgeophysical inversions.In the case of geophysical imaging, by which we are primarily motivated, thebackground medium in which the object of interest is hosted is usually heteroge-neous, and for an accurate reconstruction of the object, good prior knowledge ofthe background is necessary. Typically, the object of interest exhibits a high con-trast with the background but its response cannot be completely decoupled fromthe background response. Therefore, incorrect background estimation will intro-duce error into the estimate of the target object. A small amount of work has beendone on quantification of uncertainties in shape reconstruction (see e.g. [21, 50])but to our knowledge no work has specifically examined or tried to mitigate theeffect of uncertainty in background estimation.There are important practical applications in which it is reasonable to assumethat one might have some reliable but imperfect background information based ondata other than that used in the shape reconstruction problem. An example fromgeophysical imaging is time-lapse inversion of hydraulic fractures in the earth’ssubsurface. The physical properties of the background rocks before fracturing maybe imaged by electromagnetic or seismic geophysical surveying. Then, after in-jection of conductive fluid into the earth during the hydraulic fracturing process itis possible to estimate the volume of rock effectively fractured by electromagneticimaging, modelling the fractured zone as a thin conductive rectangular prism em-bedded in the previously estimated heterogeneous background medium [48]. Inboth these examples there will be uncertainty in the estimate of the backgroundused in the shape reconstruction.76In shape reconstruction with uncertain background information, the parametersdescribing the background model are effectively nuisance parameters. They arenot of direct interest but the accuracy of their estimation affects the accuracy of theshape reconstruction. To our knowledge, there has been very little work done ingeneral on estimating large space nuisance parameters in inverse problems. Theone notable example of which we are aware, which is very different in characterfrom the present contribution, is [5]. They use variable projection to eliminatetuning parameters from regularized least squares inverse problems. Their techniqueonly applies when the number of nuisance parameters is very small, relative tothe number of parameters of interest. We are interested in the opposite regime,where the object of interest is described by a small number of parameters and thebackground is a cell by cell discretization of a continuous function onto some sortof computational mesh, making the number of nuisance parameters very large.In this paper we present a stochastic optimization based method for improvingthe quality of shape recovery in parametric level-set shape reconstruction problemsin the presence of uncertainty in the background parameters. In section 5.2 we willdescribe the method. In section 5.3 we present a simple linear inverse problem thatwe use as a first application. This example illustrates the properties of our algo-rithm and examine the relative strengths and weaknesses of the two optimizationtechniques we propose for its solution. Section 5.4 discusses an application to anon-linear three-dimensional inverse problem, before concluding and summarizingthe paper.5.2 MethodIn general, we are concerned with the inverse problem of estimating the location,shape, and material properties of a homogeneous object embedded in a hetero-geneous background medium in two-dimensional (2D) or three-dimensional (3D)space from a noisy set of indirect observations d ∈ Rd . We assume that the dataare sensitive to some physical property θ of the material being probed. The dataare related to θ through the forward problemd =F (θ)+ ε, (5.1)77where ε is additive noise and F is the forward modelling operator, which is typi-cally a functional of the solution of a partial differential equation (PDE).In level-set shape reconstruction methods, the boundary of the object of interestis represented by the zero level-set of a function τ(x). The overall physical propertymodel can then be writtenθ(x) = θ0(x)+H(τ(x))(θ1−θ0(x)) , (5.2)where H is the Heaviside step function θ1 is the value of θ inside the object ofinterest and θ0 is the heterogeneous background physical property model. In tra-ditional level-set methods, the level-set function is discretized onto a grid and theinverse problem attempts to reconstruct it directly.In parametric methods, the level-set function is represented with a low dimen-sional parametrization. For example, in the application problem in section 5.3 ofthis paper we model the object of interest as an ellipse, which can be described byits centre position, radii, and angle of orientation. The level set function for thismodel can be writtenτ(x) = 1− (x−x0)T M(x−x0), (5.3)where x0 is the the location of the centre of the ellipse and M is a rotation andscaling matrix that depends on the radii of the ellipse and its angle of rotationabout the horizontal axis of the mesh.The goal of the parametric inverse problem is to estimate the shape and locationparameters that describe the object of interest as well as the value of the diagnosticphysical property inside the object. It should be noted here that in practice, in bothtraditional and parametric level-set methods, the Heaviside function in eq. (5.2) isreplaced by a smooth approximation in order to facilitate optimization. This willbe discussed further in the application sections.We formulate the inverse problem as the least squares optimization problemminimizemΦ(m) =12‖W(F (m)−d)‖22subject to ml < m < mh.(5.4)78We seek the parametric model m that minimizes the discrepancy between mea-sured data and simulated data computed by solving the forward problem—subjectto bound constraints on the parameter values. W in eq. (5.4) is a diagonal ma-trix containing the inverse standard deviations of the data. We assume that ourparametrization of the level-set function is sufficiently simple that it will have aregularizing effect and that therefore explicit regularization will not be required.This is a well known and well studied approach—see e.g. [1, 58, 64, 83]. Theoptimization problem eq. (5.4) can be solved by standard techniques. We havehad success in solving such problems with the projected Gauss-Newton method,solving for the search direction using the preconditioned conjugate gradient (PCG)algorithm with large stopping tolerance.Of course, perfect knowledge of the background is unlikely. To account foruncertainty in knowledge of the background, we model it as a random field. Then,since the objective function in eq. (5.4) is a function of the background model, wetreat it as a random variable. This recasts the problem of estimating the model pa-rameters of interest m as a stochastic optimization problem. Rather than minimiz-ing the objective function for a fixed background, we seek the m that minimizes theexpected value of the objective function with respect to the stochastic backgroundmodel.Let the background physical property model be a random field θ0(x,ω), whereω is an element of the sample space Ω of a probability space (Ω,F,P). For fixed ω(i.e. a realization of the random field), θ0 is a deterministic function of space. Wecan write the expected value of our now random least squares objective function asa Lebesgue integral with respect to the probability measure PE [Φ(m)] =∫Φ(m)dP(θ0). (5.5)We assume that P is continuous and that we can efficiently compute realizationsof the background. Having defined the expectation we can define our stochasticshape-reconstruction problemminimizemE [Φ(m)]subject to ml < m < mh.(5.6)79Stochastic optimization has been widely used in operations research for per-forming optimization under uncertainty in the formulation of or input to optimiza-tion problems—see e.g. [27, 95]. To our knowledge, it has not been used for thesolution of inverse problems. A version of stochastic optimization has been usedto reduce the computational complexity of solving large scale inverse problems in-volving many PDEs—e.g. [42, 78]. We are not aware of its use in improving thequality of the solution for an inverse problem.The high-dimensional integral in the expectation eq. (5.5) cannot be evaluatedanalytically for most problems and so it must be approximated. We use the sam-ple average approximation (SAA) to estimate this expectation and form a tractableoptimization problem. The SAA approximates the expectation using a sample av-erage [86],E[Φ(m)]≈ En[Φ(m)] = 1nn∑i=1Φ(m,θ (i)0 ), (5.7)where Φ(m,θ (i)0 ) denotes the misfit function evaluated with respect to a particularrealization of the background model θ (i)0 . Applying this approximation to eq. (5.6)gives the optimization problemminimizemEn[Φ(m)] =1nN∑i=1Φ(m,θ (i)0 ). (5.8)Intuitively, the misfit function is now a sample average of the individual misfitscomputed with different realizations of the background. The variance in the es-timate of the parameters of interest will decrease as the number of backgroundsamples increases.Let m∗n be the minimizer of the SAA problem eq. (5.8) and m∗ the minimizerof the true problem. The theory of stochastic programming tells us [86] that m∗nconverges in probability to m∗ such that‖m∗n−m∗‖=O(1/√n). (5.9)Despite the slow theoretical convergence rate, we have found that in practice wecan improve results in shape reconstruction problems using a number of samplesthat is computationally tractable.80Given a set of background samples {θ (i)0 }, eq. (5.8) may be solved using meth-ods of deterministic optimization. We have used the Gauss-Newton (GN) and ac-celerated mini-batch stochastic gradient descent (SGD) algorithms. In the nextsection we will analyze the performance of these two optimization approaches ona simple linear inverse problem. We note that, computationally, each backgroundsample requires separate forward problem, gradient, and possibly Hessian vectorproduct evaluations at each iteration of the optimization procedure but that thecomputations for each background are trivially parallel. Looking ahead, we havefound that given sufficient computing resources, the GN approach provides an effi-cient and robust solution while SGD can still give good results while using limitedcomputational resources.5.2.1 Computing Background SamplesBefore turning to an analysis of the performance of the SAA approach on exam-ple problems, we consider the problem of choosing background model probabilitydistributions and generating sample background models. Ideally the distributionswould be determined empirically, based on real examples from past experiments.This is not typically possible in geophysical imaging. As an alternative, we modelthe background as a Gaussian random field (GRF) with mean and covariance func-tions that encode our prior knowledge of the background. The stochastic opti-mization methods used in this paper do not depend on the background distributionbeing Gaussian. We use GRFs in our examples because of their flexibility and thesimplicity of sampling from them.We use a truncated Karhunen-Loe`ve expansion (KLE) to approximately sam-ple from GRFs. The KLE provides a way to represent a random field via a spec-tral expansion of its covariance function [35]. Consider a Gaussian random fieldm0(x,ω) with mean m0 and covariance function κ(x,y), where x ∈ D ⊆ Rs is thespatial variable and ω is an element of the sample space in some probability space(Ω,F,P). The KLE of m0 ism0(x,ϑ) = m0+∞∑n=0ξn(ω)√λnϕn(x), (5.10)81where the λn and ϕn are the eigenvalues and eigenfunctions of the covariance func-tion and the ξn(ω) are normally distributed random variables. For covariance func-tion κ(x,y), the eigenvalues and eigenfunctions solve the eigenvalue problem∫Ωκ(x,y)ϕn(y)dy = λnϕn(x). (5.11)In practice, except for special covariance functions, this must be solved approxi-mately. We use the Ny¨strom method [14] to approximate the KLE eigenvalues andeigenfunctions. This method approximates the integral eq. (5.11) by a quadraturerule. This leads to an nλ × nλ matrix eigenvalue problem that gives the first nλeigenvalues and the values of the corresponding eigenfunctions at the quadraturepoints. The eigenfunctions can then be evaluated at points of interest by interpo-lation. See [14] for a detailed description. For smaller problems the quadraturepoints may be evenly distributed over the domain and the eigenfunctions evalu-ated for each cell in the computational domain. However, for larger 3D problems,distribution of quadrature points over the full domain becomes prohibitively expen-sive. We address this problem in the discussion of our non-linear 3D applicationproblem in section 5.4.5.2.2 ImplementationThe software required to perform the numerical experiments in the remainder ofthis paper was developed in the Julia programming language [15], using the jInvframework [80]. jInv is a set of open source Julia software packages designed tosimplify the process of building software for PDE constrained parameter estima-tion problems by providing modular building block routines commonly needed inthese problems, including parallelized Gauss-Newton optimization and tools for fi-nite volume PDE discretization. We used these building blocks to develop routinesto solve our specific applications problems. To perform the Karhunen-Loe`ve basisfunction computations we used the open source Julia package GaussianRandom-Fields [77], modified to our needs.82Transmitter holeReceiverholeFigure 5.1: Cartoon schematic of straight ray tomography experiment show-ing straight ray paths from the transmitters to the receivers.5.3 Application to a Simple Model ProblemWe use two-dimensional (2D) straight ray cross-hole seismic tomography as anexample problem that is computationally simple enough to permit extensive exper-imentation. Straight ray tomography is a linearization of seismic tomography inwhich seismic waves (pressure waves moving through the earth’s subsurface) areassumed to propagate in straight lines from the source of an excitation to a receiverpoint [18]. Transmitters are placed in a borehole in the earth and receivers in a sep-arate borehole. The experiment then measures the times required for the signal topropagate from each transmitter to each receiver. The speed of wave propagationdepends on the material properties of the propagation medium, allowing seismictomography to be used as an imaging method. See fig. 5.1 for an illustration of atypical survey setup.The straight-ray tomography forward problem ist = As. (5.12)The data t are the travel times between each source and receiver pair. s is theinverse velocity, or slowness of each mesh cell, and each entry of A, ai j, givesthe distance travelled by the straight ray through cell j along the path connectingtransmitter-receiver pair i.83We attempt to recover an elliptical anomaly embedded in a heterogeneousbackground from travel time data. Our synthetic experiment is performed on aunit domain in R2, with 50 transmitters spaced evenly along the left edge of thedomain and 50 receivers spaced evenly along the right edge.The background slowness model is described by a GRF with mean and covari-ance based on the geological scenario of a horizontally layered earth with seismicwave velocity increasing with depth—a scenario commonly encountered in seis-mic imaging applications [66]. The mean slowness model is composed of fivehomogeneous layers. It is shown in fig. 5.2a. The covariance function is0 20 40 60 80x (m)020406080y (m)234567slowness (s/m)(a) mean background model0 20 40 60 80x (m)020406080y (m)−1.0−0.50.00.51.0slowness (s/m)(b) zero mean perturbation0 20 40 60 80x (m)020406080y (m)234567slowness (s/m)(c) true background model0 20 40 60 80x (m)020406080y (m)2345678slowness (s/m)(d) total true modelFigure 5.2: True model used for straight ray tomography inversions. Themean model a) consists of homogeneous horizontal layers. The ran-dom perturbation shown in b) is added to the mean model to form c),used as the background in the true model shown in d).κ(x,y) = σ exp(−(x−y)T A(x−y)), (5.13)84withA =(l−1x 00 l−1y), (5.14)for characteristic length scales lx and ly. σ is the standard deviation, which isconstant, giving the overall scale of fluctuations from the mean background. Inthis example we use σ = 1, lx = 0.3 and ly = 0.1, encoding our assumption that theslowness should vary more smoothly in the horizontal direction than the vertical.The slowness model is completed by adding a homogeneous ellipse to the back-ground model. The boundary of the ellipse is given by the parametric level-setfunctionτ(x) = 1− (x−x0)T M(x−x0), (5.15)where x0 is the the location of the centre of the ellipse and M is a rotation and scal-ing matrix that depends on the radii of the ellipse and its angle of rotation about thehorizontal axis of the mesh. For more details on the level-set parametrization, see[64]. In order to use gradient based optimization techniques, the slowness modelmust be differentiable with respect ellipse parameters across the entire domain. Toachieve this we use a smooth approximation to the Heaviside function to describethe transition from the ellipse to the background medium. The total slowness any-where in the domain is given (for fixed background slowness model s0) bys(x,s0) = s0+12(1+ tanh(aτ(x)))(se− s0(x)), (5.16)where se is the slowness inside the ellipse and a is a scaling factor that controls thesharpness of the transition between the ellipse and the background.The model used to generate synthetic data for our inversion consisted of anellipse in the centre of the domain embedded in a randomly chosen slowness modelfrom the background GRF. We discretized the domain onto a 100× 100 regulargrid of square cells. A 176 term KLE was used to compute a sample from thezero mean GRF with covariance eq. (5.13). This was then added to the meanbackground to generate a sample from the background model GRF. The zero meanperturbation and resulting background sample are shown in figs. 5.2b and 5.2c,respectively. To this background model we added an ellipse with an anomalously85high slowness of 8 s/m. The ellipse was centred at the centre of the domain with anangle of 30◦ from vertical and radii of 0.1 m and 0.2 m. This total model is shown infig. 5.2d. In the rest of this section this will be referred to as the true model. Usinga random sample from the background GRF, rather than the mean background asthe true model in our inversion study implies that the sample backgrounds usedduring inversion are being drawn from a distribution with the correct covariancefunction but wrong mean. In practice both the mean and covariance could only beapproximately known. We have not attempted to assess the success of our methodwhen the background covariance function is only approximately known.Generally, the accuracy of the background model mean and covariance func-tions required for our SAA technique to be effective in practice will depend onspecific characteristics of the application problem, and in particular, on how wellseparated the material properties of the object of interest are from the background.When the vast majority of the signal in the data comes from the object of interest,then adequate shape reconstruction might be possible without particularly precisecharacterization of the background. On the opposite end of the spectrum, if thebackground generates signals similar to that of the object of interest, then our ap-proach will likely fail and shape reconstruction will only be possible with a veryaccurate description of the background. It is in the middle ground where we believeour method can be most useful.5.3.1 SAA inversionWe conducted several numerical experiments to study the behaviour of our method.All these experiments attempted to fit the same synthetic data, which was gener-ated by applying the straight ray tomography forward modelling operator to thetrue model and then adding artificial independent Gaussian noise to each mea-surement. The added noise on each datum had standard deviation equal to 1% ofits magnitude. We attempted to recover the orientation angle, centre position andslowness of the anomalous ellipse. For simplicity we did not attempt to recover theradii of the ellipse. They were fixed to their true values. The starting guesses andbounds on the active inversion parameters were the same in all inversions. Theyare listed in table 5.1.86parameter starting guess lower bound upper boundφ 1.5 0 90x 0.2 0.1 0.9y 0.4 0.1 0.9s 4 0.01 15Table 5.1: Starting guesses and parameter bounds for straight ray tomographyinversion parameters. φ is the ellipse orientation angle in , x and y arethe horizontal and vertical coordinates of the centre position, and s theslowness inside the ellipse.Before studying the effectiveness of the SAA on the straight ray example prob-lem, we wanted to characterize the variability of results in single background inver-sions caused by uncertainty in background estimation. Assuming our backgroundmodel random field accurately captures the uncertainty in our knowledge of thebackground, this can be done by conducting many inversions, each using a singlerandom background. We performed 100 such inversions. For each we solved thedeterministic optimization problem eq. (5.4) using the projected Gauss-Newton al-gorithm with backtracking Armijo line-search. Inversions could end after attainingthe desired data fit as determined by a χ2 statistic, when the magnitude of the gra-dient was reduced below a threshold value relative to its initial magnitude, or whenthe line search failed to find a step that reduced the objective function. Resultsfrom all these categories were included in assessing the variance in the recoveredparameters as a function of the background. The results are summarized in fig. 5.3,which shows histograms and kernel density estimates of the recovered values ofeach inversion parameter. There were a wide range of results. Excellent recoverywas achieved for a few sample backgrounds but the variance in results was largeand many inversions achieved very poor recovery of the anomalous ellipse. Ori-entation angle was particularly poorly resolved, with the angle being pushed to itslower bound being the mode result. This seems to show that the angle is poorly re-solved by this particular straight ray tomography survey, rather than to give insightinto the effect of uncertainty in the background model. Overall, these results showthat the level of uncertainty in our knowledge of the background slowness model istoo high for parametric inversion using a fixed background to be reliably effective8740 20 0 20 40 60 80 100 120angle (degrees)0.0000.0050.0100.0150.0200.0250.030Normalized frequencytrue value = 30 degrees(a) rotation angle0.0 0.2 0.4 0.6 0.8 1.0x (m)0.00.51.01.52.02.53.0Normalized frequencytrue value = 0.5 m(b) horizontal position0.0 0.2 0.4 0.6 0.8 1.0y (m)0.00.51.01.52.02.53.0Normalized frequencytrue value = 0.5 m(c) vertical position0 5 10 15slowness (s/m)0.0000.0250.0500.0750.1000.1250.1500.175Normalized frequencytrue value = 8.0 s/m(d) slownessFigure 5.3: Results of 100 independent single background inversions. Therecovered values for each inversion parameter are plotted as histogramswith accompanying kernel density estimates. The red vertical linesshow the true parameter values. Note that the estimated densities mayhave support outside the parameter bounds. They are included in orderto give general visual impressions of the spread of inversion results andthese effects are not of concern.88ns 1 10 100φ 25.2 64.0 88.8x 0.38 0.50 0.51y 0.51 0.53 0.54s 6.96 7.97 8.36(a) sample meanns 1 10 100φ 28.2 31.2 8.29x 0.17 0.051 0.0091y 0.16 0.040 0.0073s 3.03 1.11 0.286(b) sample standard deviationTable 5.2: Marginal statistics of SAA inversion results. The sample means ofthe recovered values of each parameter are shown for the 1, 10, and 100background inversion experiments in a), with the corresponding samplestandard deviations shown in b). Recall that the true parameter values areφ = 30◦, x = 0.5, y = 0.5, and s = 4. See table 5.1 for parameter boundsand initial guesses.for this particular problem.To show the effect of the SAA approach in improving ellipse recovery in thepresence of uncertainty in background estimation, we performed repeated multiplebackground SAA inversions and studied the variability of the inversion results. Foreach of two SAA sample sizes ns = 10,100 we performed 100 inversions, with nsnew independent background models sampled from the background GRF. Kerneldensity plots showing the range of recovered parameter values for the SAA andsingle background inversions are plotted together in fig. 5.4. It is clear from theseresults that the SAA approach did not help recovery of the orientation angle—inline with our expectations—but that it did a great deal to reduce the variance in therecovered values of the other parameters. The density plots also show that thereare small but distinct biases in the distributions of recovered values. This stemsfrom the fact that we don’t know the true background model distribution and arethus not drawing our sample backgrounds from the correct distribution. Thus as theSAA sample size is increased we see the sample means of the recovered parametersconverging to values slightly different than the true values.The sample means and standard deviations of the recovered parameter valuesfor each SAA sample size are shown in table 5.2. The values in the table showthe O(1/√ns) convergence in standard deviation expected from the SAA. Overallthese statistics show that a single SAA inversion with a large enough sample sizecan be expected to give much better parameter recovery than a single one back-8925 0 25 50 75 100 125angle (degrees)0.0000.0050.0100.0150.0200.025Normalized frequency1 background10 backgrounds SAA100 backgrounds SAAtrue parameter value(a) rotation angle0.0 0.2 0.4 0.6 0.8 1.0x (m)0.02.55.07.510.012.515.017.5Normalized frequency(b) horizontal position0.0 0.2 0.4 0.6 0.8 1.0y (m)0246810121416Normalized frequency(c) vertical position0 5 10 15slowness (s/m)0.00.10.20.30.40.5Normalized frequency(d) slownessFigure 5.4: Kernel density plots showing inversion results for single back-ground, 10 background SAA, and 100 background SAA inversions. 100background angle results are not shown because they took the upperbound value of 90◦ for all but two of these inversions rendering thecorresponding kernel density plot meaningless. The vertical axes of theplots were clipped in order to show the single and 10 background resultsat reasonable scales.90ground inversion.5.3.2 Optimization by stochastic gradient descentA major downside of the SAA method described in this paper, when using tra-ditional deterministic optimization methods, is that the computational cost growslinearly with the number of samples ns, while the parameter estimates convergeonly at a rate of O(1/√ns). Each term in the SAA objective function requiresseparate forward modelling, gradient and hessian computations. When large scaleparallel computing resources are available, a large number of terms can be handledin parallel, allowing the approach to remain efficient. But when such resourcesare not available, the computational cost per iteration becomes prohibitive as thenumber of samples increases.In order to gain the statistical benefit of using a large number of samples inthe SAA technique without incurring a prohibitively large computational cost, weturn to stochastic gradient methods. These algorithms converge much more slowlythan the Gauss-Newton method but the per-iteration cost is drastically reduced andthey can perform SAA inversions with large numbers of background samples whenonly modest computational resources are available since only a small number ofbackground models are accessed at each iteration.Stochastic gradient descent (SGD) and related algorithms such as Polyak’s mo-mentum method [72], have been used in stochastic programming for a long time[86]. Recently, these methods, along with other variants, have gained great pop-ularity in the machine learning community due to their ability to work with smallsubsets of large datasets and their surprising performance on difficult non-convexmachine learning optimization problems—see e.g. [17, 57]. We are interested inSGD and related methods because of their ability to use information from a largenumber of sample background models without needing to process many of thesemodels at each iteration.We use the ADAM accelerated stochastic gradient descent algorithm [55]. TheADAM method combines ideas from momentum methods and gradient descentwith adaptive step size control. The method is shown in algorithm 2. The modelis updated by an exponentially weighted average of past gradients, with the indi-91Algorithm 2 The ADAM method, algorithm presentation adapted from algorithm1 of [55]Choose hyperparameters α , β1, β2, ε .Choose initial model m0 (Abuse of notation, do not confuse with backgroundmodel)Initialize v0 = 0, w0 = 0, k = 0Draw a set Ξ= {ξi} of ns samples from the background model distributionwhile not converged doCompute a random permutation of the samples and divide into mini-batchesfor j in mini-batches dok← k+1g← Φ[m,Ξ( j)] (gradient of objective function w.r.t. samples in mini-batch j).vk← (1−β t1)−1(β1vk−1+(1−β1)g)wk← (1−β t2)−1(β2wk−1+(1−β2)(g◦g))mk←mk−1−αvk/(√wk + ε)end forend whilereturn mkvidual components scaled by an exponentially weighted average of the magnitudesof past gradients. After reviewing several accelerated gradient descent algorithms,we chose ADAM because of its proven performance on many difficult large-scaleproblems and the fact that in our experience it seems to perform well with lesshyperparameter tuning than other related methods. In order to achieve satisfactoryperformance it was still necessary to tune the step-size. After some experimen-tation, we settled on a value of α = 0.002, which we used in the experimentsdiscussed below. We set the other parameters to the default values recommendedin [55]. These are β1 = 0.9, β2 = 0.99, and ε = 1×10−8.At each iteration of the optimization, the approximate gradient of the SAA ob-jective function is computed using a subset of the available background samples,called a mini-batch. We studied the performance of the ADAM method in solvingthe SAA optimization problem eq. (5.8) for the straight ray tomography problemby studying the variability of the recovery of ellipsoid parameters obtained withindependent sets of background samples, as we did with GN. We studied the per-92ns 100 500φ 24.1 26.3x 0.49 0.50y 0.51 0.52s 9.82 9.77(a) sample meanns 100 500φ 1.10 0.515x 0.0068 0.0038y 0.0077 0.0038s 0.288 0.139(b) sample standard deviationTable 5.3: Marginal statistics of SAA inversion results using the ADAM op-timization algorithm. Recall that the true parameter values are φ = 30◦,x= 0.5, y= 0.5, and s= 4. See table 5.1 for parameter bounds and initialguesses.formance of ADAM using SAA sample sizes of ns = 100,500, with mini-batchsizes of 1 and 10. We performed 100 inversions in each of these configurations.The optimizations generally converged faster with the larger mini-batch size, asexpected. The parameter recovery was equivalent in both cases. The summarystatistics of the results (computed from the 1 mini-batch runs) are shown in ta-ble 5.3. Kernel density estimates for each parameter value, for both SAA samplesizes are plotted in fig. 5.5. The standard deviations of the ellipsoid position andslowness from the 100 background inversions was very close to what was observedin the 100 background GN experiments. Unexpectedly, the orientation angle wasmuch better resolved by the ADAM inversions. The biases of the x and y posi-tion variables were similar to the GN results while the ADAM inversions tended tooverestimate the slowness more than the GN inversions.We cannot definitively say why the biases in the ADAM results were differentfrom the GN biases. It may be because the ADAM algorithm is better able to escapefrom local optima and saddle surfaces, or that because of its much smaller stepsizes and the scaling of the problem, ADAM is able to better explore the interiorof the parameter space and find reasonable values for the position and slownessparameters without quickly pushing the angle to its upper bound value.5.4 Application to non-linear 3D problemWe illustrate the application of our method to a larger scale non-linear problemusing another example from geophysics, direct current resistivity (DCR) imaging.9325 0 25 50 75 100 125angle (degrees)0.00.10.20.30.4Normalized frequency100 backgrounds500 backgroundstrue value(a) rotation angle0.0 0.2 0.4 0.6 0.8 1.0x (m)01020304050607080Normalized frequency(b) horizontal position0.0 0.2 0.4 0.6 0.8 1.0y (m)010203040506070Normalized frequency(c) vertical position0 5 10 15slowness (s/m)0.00.20.40.60.81.01.21.41.6Normalized frequency(d) slownessFigure 5.5: Kernel density plots of inversion results for 100 background and500 background ADAM inversions with mini-batch size 1. Horizontalaxis limits are the same as in fig. 5.494In DCR imaging, a steady electric current is injected into the ground and resultingelectrical potential differences due to electric currents flowing in the ground aremeasured at a small number of locations on the earth’s surface. These voltages canbe used to estimate the electrical conductivity (inverse of resistivity) of the earth’ssubsurface. Mathematically, DCR surveys are governed by the elliptic boundaryvalue problem∇ · (θ∇φ) =−∇ · js ∈ (Ω⊂ R3) (5.17)∇φ · nˆ ∈ ∂Ω1 (5.18)φ = 0 ∈ ∂Ω2 (5.19)where φ is the electrical potential in the earth, θ the conductivity, here assumedto be a scalar function of space, and js, the electric current driving the experi-ment, which is divergence free except at a finite number of isolated points wherethe source of current is connected to the ground. ∂Ω1 denotes the portion of thedomain boundary that represents the surface of the earth and nˆ is the unit vectornormal to ∂Ω. ∂Ω2 is the rest of the domain boundary. The Dirichlet conditionsmodel the fact that the potential tends toward zero far from sources of current.In order to approximately meet this condition in a finite domain, we must makethe domain large enough that the non-surface boundaries are indeed far from anysources of current. We discretize eq. (5.19) on locally refined rectilinear meshescalled OcTrees (see e.g. [45]), which makes it simple to enlarge the domain with-out significantly increasing the computational cost of solving the boundary valueproblem.The potential difference measurements are linear functionals of the potentialφ , a measurement being simply the difference in the potential between two lo-cations on the earth’s surface. The goal of this example problem is to use thesemeasurements to characterize a plate-like highly conductive body embedded in aheterogeneous vertically layered background medium. We model the conductivebody as a rectangular prism. The prism is incorporated with the background con-ductivity model using the parametric level-set approach. Similarly to the previous95example the overall conductivity model isθ(x) = θ0(x)+ H˜(x)(θp−θ0(x)), (5.20)where θ0 is the heterogeneous background conductivity model, θp is the conduc-tivity of the prism, and H˜ is a smooth approximation to an indicator function thattakes value 1 inside the prism and 0 outside. The indicator function is a prod-uct of one-dimensional (1D) indicator functions in the three Cartesian coordinatedirectionsH˜(x) = σx(R(x−x0))σy(R(x−x0))σz(R(x−x0)), (5.21)where x0 is the position of the centre of the prism and R is a rotation matrix thattransforms x− x0 to a coordinate frame whose axes are the principal axes of theprism. The 1D indicator function for coordinate i isσi(r) =11+ exp(−(ri+hi)/ai) −11+ exp(−(ri−hi)/ai) (5.22)where ri is the i component of r, hi is the side length of the prism in the i directionand ai the scaling factor that controls the sharpness of the prism boundary transitionin the i direction.The background conductivity θ0 is modelled as a log-normal random field. Themean of the underlying GRF is constant in the y direction and piecewise constant inthe x and z directions. It is illustrated in fig. 5.6 by a slice in the x-z plane showingthe core region of the model. The covariance function of log(θ0) isκ(x,y) = σ(x)σ(y)exp(−(x−y)T A(x−y)), (5.23)whereA =l−1x 0 00 l−1y 00 0 l−1z , (5.24)with lx = ly = 2.5lz. The standard deviation σ(x) is piecewise constant, being 50%of the mean value in each constant region of the mean model.We take approximate samples of the log(θ0) GRF by using the Karhunen-96− 400 − 300 − 200 − 100 0 100 200 300 400x (m )− 300− 200− 1000z (m)10− 210− 1(S/m)Figure 5.6: Mean background conductivity in the x-z plane at y = 0 in thecore region of the mesh.Loe`ve expansion. In order to make the computation of the Karhunen-Loe`ve eigen-functions tractable and to make sampling efficient, we sample the background ona subset of the full computational domain. We use a large enough volume suchthat the lack of randomization outside it does not significantly effect the data. Wethen simply take the exponential of these GRF samples to get samples from thebackground random field θ0.The 3D conductivity model used to generate synthetic data for our experiment(which we’ll call the true model below) was taken to be a thin dipping rectangularprism embedded in a random sample from the background field. It is illustrated byorthogonal 2D slices in fig. 5.7.With our computational resources, using the Gauss-Newton SAA inversiontechnique with many sample backgrounds was computationally intractable for thisproblem. However, we found that the stochastic gradient technique introduced inthe previous section worked well. As in the last section we performed a set ofinversions using single random samples from the background distribution. Dueto the computational cost, we did not perform a series of SAA inversions. Weperformed a single SAA inversion using the ADAM stochastic optimization algo-rithm, with 400 background models and a mini-batch size of 4. This mini-batchsize was chosen simply because this was the number of model realizations that wecould process concurrently.In this case the anomalous parametric body was better separated from the back-ground than in the 2D example. Because of this, the single background DCR inver-97− 400 − 300 − 200 − 100 0 100 200 300 400y (m )− 300− 200− 1000z (m)10− 210− 1(S/m)(a) Slice through model at x =−25 m− 400 − 300 − 200 − 100 0 100 200 300 400x (m )− 300− 200− 1000z (m)(b) Slice through model at y = 0 m− 300 − 200 − 100 0 100 200 300x (m )− 300− 200− 1000100200300y (m)(c) Slice through model at z =−60 mFigure 5.7: True model slices98sions gave better results than in the 2D tomography case. However, there was stilla significant variability in the results. This is illustrated in fig. 5.8, which showsthe recovered parameter values for the single background and SAA inversions. Thesingle background inversions were able to somewhat reliably recover the positionof the anomalous block but showed much larger variability in the recovery of theother parameters. None of the inversions were able to recover the thickness (z ex-tent) well, suggesting that the data may not be very sensitive to this parameter. Thesingle background inversions were able to reduce misfit by tweaking the extents ofthe block and varying its orientation angle but they struggled to capture its highlyconductive nature. In most cases they increased the signal from the block only byincreasing its thickness, and not by increasing the conductivity. The SAA inversionstill pushed the thickness to its upper bound but, perhaps because it could not finddescent directions by fitting small signals from background features, it was bet-ter able to capture the essential features of the anomalous block, leading to goodrecovery of its position, orientation, and conductivity.5.5 ConclusionsIn parametric shape reconstructions problems, estimating the properties of thebackground medium can be a major stumbling block to accurate reconstruction.We have shown that the effect of uncertainty in the background properties can bemitigated by the use of stochastic optimization methods. This can be thought ofas a means to improve the shape recovery possible given a specified characteriza-tion of the background or as a justification to limit the effort one might expend onachieving a very accurate representation of the background.The efficiency of this method could likely be improved in future. In this workwe used ADAM, a first order stochastic optimization method, to make our approachtractable for large-scale problems. There may be an opportunity to use StochasticNewton methods to improve the efficiency of optimization.990.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #255075100125150175200hx (m)(a) x extent0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #255075100125150175200hy (m)(b) y extent0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #1020304050hz (m)(c) z extent0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #−100−75−50−250255075100xc(m)(d) x centre position0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #−100−75−50−250255075100yc(m)(e) y centre position0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #−160−140−120−100−80−60−40zc(m)(f) z centre position0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #010203040ϕy(m)(g) rotation angle0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0sample #246810σ(S/m)(h) conductivityFigure 5.8: Results from single background and SAA ADAM inversions ofDC data. Blue dots show the recovered values of each parameter fromthe 20 single background inversions and the stars show the parametervalues recovered by the SAA inversion. In each plot, the vertical axislimits are the parameter bounds, the dashed black line shows the initialguess and the red line the true value of the parameter.100Chapter 6ConclusionsThis thesis has focused on the development of computational methods relatedto forward modelling and inversion of time-domain geophysical electromagnetic(EM) data in the presence of induced polarization effects. The work led to new,more efficient numerical methods for time discretization of the quasi-static Maxwellequations, a new forward modelling algorithm for simulating coupled EM induc-tion and induced polarization (IP) effects, and a corresponding inversion algorithmfor recovering three-dimensional (3D) conductivity and IP parameter earth modelsfrom transient EM data. The final main chapter concerned stochastic parametricinversion. It did not address IP directly but contributed a novel inversion algorithmthat is well suited for application to IP inversion problems in the future.This work was motivated by interest in induced polarization in the geophysicalcommunity and by a lack of appropriate computational tools for the investigationof coupled EM induction and IP effects. The first attempts to model these coupledeffects in 3D directly in the time-domain have emerged only in the last five years,in the work of Marchant et al. [63] and Commer et al. [25]. This work sought toimprove on the computational efficiency of these methods and to develop a purelytime-domain EM and IP forward modelling algorithm scalable to large industrialproblems that is amenable to use in inversion for IP parameters.Efficient modelling of coupled EM induction and IP effects clearly requiresefficient modelling of electromagnetic induction. Therefore, before addressingIP, Chapter 2 of the thesis examined the discretization in time of the quasi-static101Maxwell equations. Time-discretization in geophysical electromagnetics has tra-ditionally used quite simple methods. As a parabolic problem, solutions tend todisplay simple smoothly decaying behaviour in time but short time-scale effectsmakes the problem very stiff. In order to maintain stability one is forced to use dis-cretization methods that possess the property of stiff decay. In addition to stability,these methods have the advantage that they can use large time-steps to recover thelate time behaviour of a solution without resolving the details of its early time be-haviour. This key property allowed the development of the parallel in time forwardmodelling algorithm discussed in chapter 2. This simple modification eliminatesthe bottleneck in speeding up forward modelling caused by the inherently sequen-tial nature of time-stepping. A key difference between the parallel time-steppingapproach and lower level optimizations such as parallelization of the solution of thelinear systems at each time-step is its scalability. The tests performed in chapter 2showed that the approach could scale to multiple nodes of a computer cluster andachieve significantly faster forward modelling times than traditional methods.Chapter 3 detailed the development of the coupled EM and IP forward mod-elling algorithm based on stretched exponential relaxation. The algorithm was im-plemented using finite volume spatial discretization on OcTree meshes. BackwardEuler and second order backward differentiation formula (BDF2) time-steppingschemes were examined. Using the stretched exponential relaxation to model theIP effect was a departure from the Cole-Cole model that has become the de-factostandard model of IP in the geophysical community. These models are both phe-nomenological. The analysis of chapter 3 showed that they are not fully equivalentbut can be used to model similar phenomena and are equivalent in some cases.The examples in chapter 3 showed that the stretched exponential model is capableof simulating electromagnetic fields qualitatively similar to those observed in thefield by practicing geophysicists. However, the ability of the stretched exponentialapproach to model real data remains untested. This is an important area for fu-ture work. The stretched exponential forward modelling algorithm and the paralleltime-stepping algorithm described in chapter 2 were published in [11].Ultimately, the usefulness of the stretched exponential approach should bejudged by its ability to generate useful earth models that fit experimental transientEM data through the process of geophysical inversion. Chapter 4 showed how102stretched exponential modelling could be used in an inversion algorithm to recoverintrinsic IP parameters from grounded source transient EM data. The work waspreliminary and tested only on a simple synthetic example. It leaves several openquestions on how to best extract information on the IP properties of the earth fromtransient EM data. Some of these questions can be addressed by further syntheticinversion studies, such as the question of if the non-uniqueness of the inversioncan be mitigated through regularization that couples the three IP parameters? Thecircumstances, if any, under which conductivity and chargeability can be recoveredsimultaneously, without a good starting conductivity model could also be investi-gated in synthetic studies? However, any such synthetic studies should be coupledwith applications to field data.A further avenue for mitigating the non-uniqueness in the IP inverse problemis to use parametric methods, which pose the inverse problem as that of recoveringa small number of parameters that characterize a homogeneous body embeddedin some background medium. This can be particularly useful in the common casewhere there is a priori reason to believe that the IP signal observed in a geophysicaldata-set is the result of a relatively homogeneous compact body. This scenario wasexplored in chapter 5 of the thesis, through the development of a stochastic para-metric level-set inversion algorithm. In this class of methods, the domain of inves-tigation into the homogeneous body of interest and a heterogeneous backgroundmedium that is held fixed in the inversion. The work in this thesis made two maincontributions to this problem. It showed that better recovery of the homogeneousobject of interest could be achieved by addressing the uncertainty in the model ofthe background medium. This was done by modelling the background as a ran-dom field and posing the inverse problem as a stochastic programming problem.The initial implementation of this approach was quite computationally intensive.It was shown that the algorithm could be adapted to make it scalable to large, in-dustrial scale problems by employing stochastic gradient descent in the solution ofthe stochastic programming problem. The work in chapter 5 has been submittedfor publication in the SIAM/ASA Journal on Uncertainty Quantification.The stochastic parametric level-set inversion approach could be used to helpthe stretched exponential inverse problem in multiple ways that could be exploredin future work. One potential use is to model the earth conductivity as part of the103background model when it is not possible to get a precise estimate of the conduc-tivity before inverting for IP parameters. This would hopefully reduce the qualityof the starting conductivity model needed to produce good models of the IP param-eters. One could also use the stochastic approach to model a possibly polarizablebackground medium that has not been well characterized, when a compact conduc-tive target is of interest.Overall the problem of inverting transient EM data for IP parameters is stillnew and under explored. This thesis has developed computational methods thatcan now be used to investigate the recovery of IP information from real data. IPeffects seem to be showing up more commonly in recent inductive source EMdata and whether the target of interest in a particular survey is polarizable or ifa conductive target is simply masked by IP signals from background media, it isimport to understand how the presence polarizable bodies is manifest in EM dataand to incorporate this understanding into the modelling and inversion tools thatare used to analyze the data. The study of IP in grounded source electrical data ismore established but there is potential for much more information to be extractedfrom such data through the use of modelling and inversion tools, such as thosedeveloped in this thesis, that fully account for the behaviour of these data in time.104Bibliography[1] A. Aghasi, M. Kilmer, and E. L. Miller. Parametric Level Set Methods forInverse Problems. SIAM Journal on Imaging Sciences, 4(2):618–650,2011. doi:10.1137/1008002081.[2] F. Alvarez, A. Alegria, and J. Colmenero. Relationship between thetime-domain Kohlrausch-Williams-Watts and frequency-domainHavriliak-Negami relaxation functions. Physical Review B, 44(14):7306–7312, 1991.[3] P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L’Excellent. A fullyasynchronous multifrontal solver using distributed dynamic scheduling.SIAM Journal of Matrix Analysis and Applications, 23(1):15–41, 2001.[4] S. Ansari and C. G. Farquharson. 3D finite-element forward modeling ofelectromagnetic data using vector and scalar potentials and unstructuredgrids. Geophysics, 79(4):E149–E165, 2014. ISSN 0016-8033.doi:doi:10.1190/geo2013-0172.1.[5] A. Y. Aravkin and T. Van Leeuwen. Estimating nuisance parameters ininverse problems. Inverse Problems, 28, 2012. ISSN 02665611.doi:10.1088/0266-5611/28/11/115016.[6] U. Ascher and L. Petzold. Computer Methods for Ordinary DifferentialEquations and Differential-Algebraic Equations. SIAM, Philadelphia, PA,1998.[7] U. M. Ascher. Numerical Methods for Evolutionary Differential Equations.SIAM, 2008. ISBN 9780898716528.[8] U. M. Ascher and C. Greif. A First Course in Numerical Methods. SIAM,Philadelphia, PA, 2011. ISBN 3717089528.doi:10.1137/1.9780898719987.105[9] P. Belliveau. Parallelizing the 2.5D Airborne Electromagnetic InversionProgram ArjunAir. M.sc. thesis, Memorial University of Newfoundland,2014.[10] P. Belliveau and E. Haber. Achieving depth resolution with gradient arraysurvey data through transient electromagnetic inversion. In SEG TechnicalProgram Expanded Abstracts 2016, number 1, pages 857–861, 2016.[11] P. Belliveau and E. Haber. Coupled simulation of electromagneticinduction and induced polarization effects using stretched exponentialrelaxation. Geophysics, 83(2):WB109–WB121, 2018. ISSN 0016-8033.doi:10.1190/geo2017-0494.1. URLhttps://library.seg.org/doi/10.1190/geo2017-0494.1.[12] K. M. Bennett, K. M. Schmainda, R. B. (Tong), D. B. Rowe, H. Lu, andJ. S. Hyde. Characterization of continuously distributed cortical waterdiffusion rates with a stretchedexponential model. Magnetic Resonance inMedicine, 50(4):727–734, 2003.[13] M. N. Berberan-Santos, E. N. Bodunov, and B. Valeur. Mathematicalfunctions for the analysis of luminescence decays with underlyingdistributions: 1. Kohlsrausch decay function (stretched exponential).Chemical Physics, 315(1-2):171–182, 2005. ISSN 03010104.doi:10.1016/j.chemphys.2005.05.026.[14] W. Betz, I. Papaioannou, and D. Straub. Numerical methods for thediscretization of random fields by means of the Karhunen-Loe`veexpansion. Computer Methods in Applied Mechanics and Engineering,271:109–129, 2014. ISSN 00457825. doi:10.1016/j.cma.2013.12.010.URL http://dx.doi.org/10.1016/j.cma.2013.12.010.[15] J. Bezanson, A. Edelman, S. Karpinski, and V. Shah. Julia: A FreshApproach to Numerical Computing. SIAM Review, 59(1):65–98, 2017.URL https://arxiv.org/abs/1411.1607.[16] R.-U. Bo¨rner, O. G. Ernst, and S. Gu¨ttel. Three-dimensional transientelectromagnetic modelling using Rational Krylov methods. GeophysicalJournal International, 202(3):2025–2043, 2015. ISSN 0956-540X.doi:10.1093/gji/ggv224. URLhttp://gji.oxfordjournals.org/lookup/doi/10.1093/gji/ggv224.[17] L. Bottou. Large-Scale Machine Learning with Stochastic GradientDescent. In Proceedings of COMPSTAT’2010. Physica-Verlag HD, 2010.106[18] N. D. Bregman, R. C. Bailey, and C. Chapman. Crosshole seismictomography. Geophyics, 54(2):200–215, 1989.[19] M. Burger and S. J. Osher. A survey on level set methods for inverseproblems and optimal design. European Journal of Applied Mathematics,16:263–301, 2005. ISSN 09567925. doi:10.1017/S0956792505006182.[20] H. Cai, X. Hu, B. Xiong, and M. S. Zhdanov. Finite-element time-domainmodeling of electromagnetic data in general dispersive medium usingadaptive Pade´ series. Computers & Geosciences, 109(January):194–205,2017. ISSN 00983004. doi:10.1016/j.cageo.2017.08.017. URLhttp://linkinghub.elsevier.com/retrieve/pii/S0098300417300079.[21] M. Cardiff and P. K. Kitanidis. Bayesian inversion for facies detection: Anextensible level set framework. Water Resources Research, 45(10):1–15,2009. ISSN 00431397. doi:10.1029/2008WR007675.[22] K. S. Cole and R. H. Cole. Dispersion and absorption in dielectrics I.Alternating current characteristics. The Journal of Chemical Physics, 9(4):341–351, 1941. ISSN 00219606. doi:10.1063/1.1750906.[23] M. Commer and G. Newman. A parallel finitedifference approach for 3Dtransient electromagnetic modeling with galvanic sources. Geophysics, 69(5):1192–1202, 2004. ISSN 0016-8033. doi:10.1190/1.1801936. URLhttp://library.seg.org/doi/abs/10.1190/1.1801936.[24] M. Commer, G. A. Newman, K. H. Williams, and S. S. Hubbard. 3Dinduced-polarization data inversion for complex resistivity. Geophyics, 76(3):F157–F171, 2011.[25] M. Commer, P. V. Petrov, and G. A. Newman. FDTD modelling of inducedpolarization phenomena in transient electromagnetics. GeophysicalJournal International, 209:387–405, 2017. doi:10.1093/gji/ggx023.[26] T. A. Davis, S. Rajamanickam, and W. Sid-lakhdar. A survey of directmethods for sparse linear systems. Technical report, Texas A&MUniversity, 2016.[27] U. Diwekar. Optimization Under Uncertainty. In Introduction to AppliedOptimization. Boston, MA, 2008.[28] K. V. D. Doel and U. Ascher. On level set regularization for highlyill-posed distributed parameter estimation problems. Journal of107Computational Physics, 216:707–723, 2006.doi:10.1016/j.jcp.2006.01.022.[29] V. Dolean, M. Gander, and L. Gerardo-Giorda. Optimized SchwarzMethods for Maxwell’s Equations. SIAM Journal on Scientific Computing,31(3):2193–2213, 2009. ISSN 1064-8275. doi:10.1137/080728536. URLhttp://epubs.siam.org/doi/abs/10.1137/080728536.[30] V. Druskin and L. Knizhnerman. Spectral approach to solvingthree-dimensional Maxwell’s diffusion equations in the time and frequencydomains. Radio Science, 29(4):937–953, 1994.[31] M. E. Everett. Near-Surface Applied Geophysics. Cambridge UniversityPress, Cambridge, 2013.[32] G. Fiandaca, E. Auken, A. V. Christiansen, and A. Gazoty.Time-domain-induced polarization: Full-decay forward modeling and 1Dlaterally constrained inversion of Cole-Cole parameters. Geophysics, 77(3):E213–E225, 2012. ISSN 00168033. doi:10.1190/geo2011-0217.1.[33] L. A. Gallardo and M. A. Meju. Joint two-dimensional DC resistivity andseismic travel time inversion with cross-gradients constraints. Journal ofGeophysical Research, 109(B3):B03311, 2004. ISSN 0148-0227.doi:10.1029/2003JB002716. URLhttp://doi.wiley.com/10.1029/2003JB002716.[34] M. J. Gander and S. Gu¨ttel. PARAEXP: A Parallel Integrator for LinearInitial-Value Problems. SIAM Journal on Scientific Computing, 35(2):C123–C142, 2013. ISSN 1064-8275. doi:10.1137/110856137. URLhttp://epubs.siam.org/doi/abs/10.1137/110856137{%}5Cnhttp://dx.doi.org/10.1137/110856137http://epubs.siam.org/doi/abs/10.1137/110856137{%}7B{%}25{%}7D5Cnhttp://dx.doi.org/10.1137/110856137.[35] R. G. Ghanem and P. D. Spanos. Stochastic Finite Elements: A SpectralApproach. Springer-Verlag, 1991. ISBN 9781461277958.[36] N. I. M. Gould, J. A. Scott, and Y. Hu. A numerical evaluation of sparsedirect solvers for the solution of large sparse symmetric linear systems ofequations. ACM Transactions on Mathematical Software, 33(2), 2007.doi:10.1145/1236463.1236465. URL http://doi.acm.org/10.1145.[37] E. Haber. Computational Methods in Geophysical Electromagnetics.SIAM, Philadelphia, PA, 2015. ISBN 9781611973792.108[38] E. Haber and S. Heldmann. An octree multigrid method for quasi-staticMaxwell’s equations with highly discontinuous coefficients. Journal ofComputational Physics, 223(2):783–796, 2007. ISSN 00219991.doi:10.1016/j.jcp.2006.10.012.[39] E. Haber and L. Ruthotto. A multiscale finite volume method forMaxwell’s equations at low frequencies. Geophysical JournalInternational, 199:1268–1277, 2014. ISSN 0956-540X.doi:10.1093/gji/ggu268.[40] E. Haber, U. M. Ascher, and D. W. Oldenburg. Inversion of 3Delectromagnetic data in frequency and time domain using an inexactallatonce approach. Geophysics, 69(5):1216–1228, 2004. ISSN0016-8033. doi:10.1190/1.1801938. URLhttp://library.seg.org/doi/abs/10.1190/1.1801938.[41] E. Haber, D. W. Oldenburg, and R. Shekhtman. Inversion of time domainthree-dimensional electromagnetic data. Geophysical JournalInternational, 171:550–564, 2007.[42] E. Haber, M. Chung, and F. Herrmann. An Effective Method for ParameterEstimation with PDE Constraints with Multiple Right-Hand Sides. SIAMJournal on Optimization, 22(3):739–757, 2012.[43] E. Hairer, S. N{\o}rset, and G. Wanner. Solving Ordinary DifferentialEquations I: Nonstiff Problems, volume 32. Springer, 1993. ISBN354078862X. doi:10.1016/0378-4754(87)90083-8. URL https://link.springer.com/content/pdf/10.1007{%}2F978-3-540-78862-1.pdf.[44] R. Hilfer. H-function representations for stretched exponential relaxationand non-Debye susceptibilities in glassy systems. Physical Review E -Statistical, Nonlinear, and Soft Matter Physics, 65(6):9–13, 2002. ISSN15393755. doi:10.1103/PhysRevE.65.061510.[45] L. Horesh and E. Haber. A Second Order Discretization of Maxwell’sEquations in the Quasi-Static Regime on OcTree Grids. SIAM Journal onScientific Computing, 33(5):2805–2822, 2011.[46] G. Horton and S. G. Vandewalle. A Space-Time Multigrid Method forParabolic Partial Differential Equations. SIAM Journal on ScientificComputing, 16(4):848–864, 1995. ISSN 1064-8275. doi:10.1137/0916050.URL http://epubs.siam.org/doi/abs/10.1137/0916050.109[47] G. Horton, S. G. Vandewalle, and P. Worley. An Algorithm with PolylogParallel Complexity for Solving Parabolic Partial Differential Equations.SIAM Journal on Scientific Computing, 16(3):531–541, 1994.[48] G. M. Hoversten, C. Schwarzbach, P. Belliveau, E. Haber, andR. Shekhtman. Borehole to Surface Electromagnetic Monitoring ofHydraulic Fractures. In 79th EAGE Conference and Exhibition, 2017.doi:10.3997/2214-4609.201700853.[49] J. M. Hyman and M. Shashkov. Mimetic Discretizations for Maxwell’sEquations. Journal of Computational Physics, 151:881–909, 1999.[50] M. A. Iglesias, Y. Lu, and A. M. Stuart. A Bayesian level set method forgeometric inverse problems. Interfaces and Free Boundaries, 18(2):181–217, 2016. ISSN 14639963. doi:10.4171/IFB/362.[51] S. G. Johnson. Wikimedia-YeeCube. URLhttps://commons.wikimedia.org/wiki/File:Yee-cube.svg.[52] S. Kang and D. Oldenburg. Revisiting the time domain inducedpolarization technique, from linearization to inversion. In AGU FallMeeting, 2015.[53] S. Kang and D. W. Oldenburg. On recovering distributed IP informationfrom inductive source time domain electromagnetic data. GeophysicalJournal International, 207:174–196, 2016. ISSN 0956-540X.doi:10.1093/gji/ggw256.[54] A. Kemna, A. Binley, A. Ramirez, and W. Daily. Complex resistivitytomography for environmental applications. Chemical EngineeringJournal, 77(1-2):11–18, 2000. ISSN 13858947.doi:10.1016/S1385-8947(99)00135-7.[55] D. P. Kingma and J. L. Ba. ADAM: A Method for Stochastic Optimization.In Proceedings of the 3rd International Conference on LearningRepresentations (ICLR 2015), 2015.[56] T. Kratzer and J. C. Macnae. Induced polarization in airborne EM.Geophysics, 77(5):E317–E327, 2012. ISSN 0016-8033.doi:10.1190/geo2011-0492.1.[57] Y. Lecun, Y. Bengio, and G. Hinton. Deep learning. Nature, 521(7553):436–444, 2015. ISSN 14764687. doi:10.1038/nature14539.110[58] Y. Li and D. W. Oldenburg. 3-D inversion of magnetic data. Geophysics,61:394–408, 1996.[59] J.-l. Lions, Y. Maday, and G. Turinici. A parareal in time discretization ofPDE’s. Comptes Rendus de l’Acade´mie des Sciences - Series I -Mathematics, 1(2):661–668, 2001.[60] J. Macnae. Quantifying Airborne Induced Polarization effects in helicoptertime domain electromagnetics. Journal of Applied Geophysics, 2015. ISSN09269851. doi:10.1016/j.jappgeo.2015.10.016.[61] J. Macnae. Fitting superparamagnetic and distributed Cole-Coleparameters to airborne electromagnetic data : A case history from Quebec.Geophyics, 81(6):B211–B220, 2016. doi:10.1190/geo2016-0119.1.[62] D. Marchant. Induced Polarization Effects in Inductive SourceElectromagnetic Data. Ph.d. thesis, University of British Columbia, 2015.[63] D. Marchant, E. Haber, and D. W. Oldenburg. Three-dimensional modelingof IP effects in time-domain electromagnetic data. Geophyics, 79(6), 2014.[64] M. S. Mcmillan, C. Schwarzbach, E. Haber, and D. W. Oldenburg. Threedimensional parametric hybrid inversion of airborne time-domainelectromagnetic data. Geophyics, 80(6):K25–K36, 2014.doi:10.1190/geo2015-0141.1.[65] P. Monk. Finite Element Methods for Maxwell’s Equations. OxfordUniversity Press, 2003.[66] O. Okko. Vertical increase in seismic velocity with depth in shallowcrystalline bedrock. Journal of Applied Geophysics, 32(4):335–345, 1994.ISSN 09269851. doi:10.1016/0926-9851(94)90032-9.[67] D. W. Oldenburg and Y. Li. 3-D inversion of induced polarization data.Geophysics, 59(9):1327–1341, 1994. ISSN 00168033.doi:10.1190/1.1444877.[68] D. W. Oldenburg, E. Haber, and R. Shekhtman. Three dimensionalinversion of multisource time domain electromagnetic data. Geophysics, 78(1):E47–E57, 2013. ISSN 0016-8033. doi:10.1190/geo2012-0131.1. URLhttp://apps.webofknowledge.com/full{ }record.do?product=WOS{&}search{ }mode=GeneralSearch{&}qid=4{&}SID=3D9g64O1pAi9CCJ9N1a{&}page=1{&}doc=3.111[69] D. O’Leary. The block conjugate gradient algorithm and related methods,1980. URLhttp://www.sciencedirect.com/science/article/pii/0024379580902475.[70] G. J. Palacky. Resistivity Characteristics of Geologic Targets. InElectromagnetic Methods in Applied Geophysics, Volume 1, Theory, pages53–130. Society of Exploration Geophysicists, 1987.[71] W. Pelton, S. Ward, P. Hallof, W. Sill, and P. Nelson. MineralDiscrimination and Removal of Inductive Coupling with MultifrequencyIP. Geophyics, 43(3):588–609, 1978.[72] B. T. Polyak and A. B. Juditsky. Acceleration of Stochastic Approximationby Averaging. SIAM Journal on Control and Optimization, 30(4):838–855,1992. ISSN 0363-0129. doi:10.1137/0330046. URLhttp://epubs.siam.org/doi/10.1137/0330046.[73] K. Prieto and O. Dorn. Sparsity and level set regularization for diffuseoptical tomography using a transport model in 2D. Inverse Problems, 33,2017.[74] I. T. Rekanos and T. G. Papadopoulos. FDTD Modeling of WavePropagation in Cole-Cole Media With Multiple Relaxation Times.Antennas and Wireless Propagation Letters, IEEE, 9:67–69, 2010. ISSN1536-1225. doi:10.1109/lawp.2010.2043410.[75] Z. Ren, T. Kalscheuer, S. Greenhalgh, and H. Maurer. Afinite-element-based domain-decomposition approach for plane wave 3Delectromagnetic modeling. Geophysics, 79(6):E255–E268, 2014. ISSN0016-8033. doi:10.1190/geo2013-0376.1. URLhttp://library.seg.org/doi/10.1190/geo2013-0376.1.[76] A. Revil and N. Florsch. Determination of permeability from spectralinduced polarization in granular media. Geophysical Journal International,181:1480–1498, 2010. doi:10.1111/j.1365-246X.2010.04573.x.[77] P. Robbe. GaussianRandomFields, 2017.[78] F. Roosta-khorasani, K. V. A. N. D. E. N. Doel, and U. R. I. Ascher.Stochastic Algorithms for Inverse Problems Involving PDEs and ManyMeasurements. SIAM Journal on Scientific Computing, 36(5):S3–S22,2014.112[79] P. S. Routh and D. W. Oldenburg. Electromagnetic coupling infrequency-domain induced polarization data: a method for removal.Geophysical Journal International, 145:59–76, 2001. ISSN 0956-540X.doi:10.1046/j.1365-246x.2001.00384.x.[80] L. Ruthotto, E. Treister, and E. Haber. jInvA Flexible Julia Package forPDE Parameter Estimation. SIAM Journal on Scientific Computing, 39(5):702–722, 2017.[81] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia,PA, 2 edition, 2003. ISBN 9780898715347.[82] K. Salari and P. Knupp. Code Verification by the Method of ManufacturedSolutions. Technical report, 2000.[83] F. Santosa. A Level-Set Approach for Inverse Problems InvolvingObstacles. ESAIM: Control, Optimisation and Calculus of Variations, 1(January):17–33, 1996. doi:10.1051/cocv:1996101.[84] O. Schenk and K. Gartner. Solving unsymmetric sparse systems of linearequations with PARDISO. Journal of Future Generation ComputerSystems, 20(3):475–487, 2004.[85] H. O. Seigel. Mathematical Formulation and Type Curves for InducedPolarization. Geophyics, XXIV(3):547–565, 1959.[86] A. Shapiro, D. Dentcheva, and A. Ruszczyski. Lectures On StochasticProgramming: Modeling and Theory. SIAM, Philadelphia, PA, 2009.ISBN 089871687X. doi:http://dx.doi.org/10.1137/1.9780898718751.[87] R. S. Smith. Induced polarization effects in airborne electromagnetic data:estimating chargeability from shape reversals. In SEG Technical ProgramExpanded Abstracts 2016, pages 2211–2217, 2016.[88] K. Sørensen and E. Auken. SkyTEM a new high-resolution helicoptertransient electromagnetic system. Exploration Geophysics, 35:191–199,2004.[89] D. E. Stewart. Dynamics with Inequalities: Impacts and Hard Constraints.SIAM, Philadelphia, PA, 2011. ISBN 978-1-61197-070-8.[90] D. M. Sullivan. Z-Transform Theory and the FDTD Method. IEEETransactions on Antennas and Propagation, 44(1):28–34, 1996.113[91] Y. Takayama and W. Klaus. Reinterpretation of the auxiliary differentialequation method for FDTD. IEEE Microwave and Wireless ComponentsLetters, 12(3):102–104, 2002. ISSN 15311309. doi:10.1109/7260.989865.[92] W. M. Telford, L. P. Geldart, and R. E. Sheriff. Applied Geophysics.Cambridge University Press, second edition, 1990.[93] A. Tikhonov and Arsenin. Solutions ofIll-posed Problems. Winston &Sons, Washington, D.C., 1977.[94] G. VanVoorhis, P. Nelson, and T. Drake. Complex resistivity spectra ofporphyry copper mineralization. Geophyics, 38:49–60, 1973.[95] S. W. Wallace and W. T. Ziemba, editors. Applications of StochasticProgramming. SIAM, Philadelphia, PA, 2005.[96] S. H. Ward. Resistivity and Induced Polarization Methods. In Geotechnicaland Environmental Geophysics: Volume 1: Review and Tutorial, pages147–190. Society of Exploration Geophysicists, 1990.doi:10.1190/1.9781560802785.ch6.[97] S. H. Ward and G. W. Hohmann. Electromagnetic theory for geophysicalapplications. In M. N. Nabighian, editor, Electromagnetic Methods inApplied Geophysics, Volume 1, Theory, pages 131–311. Society ofExploration Geophysicists, 1987.[98] W. Weedon and C. Rappaport. A general method for FDTDmod- elingofwave propagation in arbitrary frequency-dispersive media. IEEETransactions on Antennas and Propagation1, 45(3):401–410, 997.[99] P. Weidelt. Response characteristics of coincident loop transientelectromagnetic systems. Geophysics, 47(9):1325–1330, 1982. ISSN1070485X. doi:10.1190/1.1441393. URLhttp://geophysics.geoscienceworld.org/content/47/9/1325.abstract.[100] G. A. Wilson, A. P. Raiche, and F. Sugeng. 2.5D inversion of airborneelectromagnetic data. Exploration Geophysics, 37:363–371, 2006.[101] M. G. Wismer and R. Ludwig. An Explicit Numerical Time DomainFormulation to Simulate Pulsed Pressure Waves in Viscous FluidsExhibiting Arbitrary Frequency Power Law Attenuation. IEEETransactions on Ultrasonics, Ferroelectrics, and Frequency Control, 42(6):1040–1049, 1995. ISSN 08853010. doi:10.1109/58.476548.114[102] Z. Xu and M. S. Zhdanov. Three-Dimensional Cole-Cole Model Inversionof Induced Polarization Data Based on Regularized Conjugate GradientMethod. IEEE Geoscience and Remote Sensing Letters, 12(6):1180–1184,2015. ISSN 1545598X. doi:10.1109/LGRS.2014.2387197.[103] K. Yee. Numerical solution of initial boundary value problems involvingMaxwell’s equations in isotropic media. IEEE Transactions on Antennasand Propagation, 14(3):302–307, 1966.[104] M. Zaslavsky and V. Druskin. Solution of time-convolutionary Maxwell’sequations using parameter-dependent Krylov subspace reduction. Journalof Computational Physics, 229(12):4831–4839, 2010. ISSN 0021-9991.doi:10.1016/j.jcp.2010.03.019. URLhttp://dx.doi.org/10.1016/j.jcp.2010.03.019.[105] M. S. Zhdanov. Generalized effective-medium theory of inducedpolarization. Geophysics, 73(5):F197–F211, 2008. ISSN 00168033.doi:10.1190/1.2973462.[106] K. Zonge, J. Wynn, and S. Urquart. Resistivity, Induced Polarization, andComplex Resistivity. In Investigations in Geophysics No 13, page 265.2005.[107] K. L. Zonge and J. C. Wynn. Recent advances and applications in complexresistivity measurements. Geophysics, 40(5):851–864, 1975.115
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Forward modelling and inversion of time-domain electromagnetic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Forward modelling and inversion of time-domain electromagnetic geophysical surveys in the presence of… Belliveau, Patrick 2019
pdf
Page Metadata
Item Metadata
Title | Forward modelling and inversion of time-domain electromagnetic geophysical surveys in the presence of chargeable materials |
Creator |
Belliveau, Patrick |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | The work described in this thesis examines transient geophysical electromagnetic forward modelling and inversion in the presence of induced polarization (IP) effects. The thesis introduces a new method of modelling IP using stretched exponential relaxation. A three-dimensional (3D) forward modelling algorithm taking full account of the coupling of IP and electromagnetic induction is developed. The stretched exponential modelling algorithm has been implemented using efficient numerical methods that allow it to tackle large-scale problems and are amenable to use in inversion. In particular, a parallel time-stepping technique has been de- veloped that allows transient electric fields at multiple time steps to be computed simultaneously. The behavior of the stretched exponential model is demonstrated by applying it to synthetic numerical examples that simulate a grounded source IP survey with significant electromagnetic induction effects and a concentric-loop airborne electromagnetic sounding over a polarizable body. An inversion algorithm using the stretched exponential model was developed that is able to recover the 3D structure of physical properties of the earth related to IP from transient geophysical electromagnetic data. The method is tested on a simple synthetic example problem. The thesis finishes with the development of a novel stochastic parametric level-set inversion algorithm, which could be useful in applying stretched exponential inversion to real world problems in the future. The algorithm addresses some of the shortcomings of the simple inversion approach used for stretched exponential inversion earlier in the thesis. The stochastic parametric inversion algorithm is used to solve shape reconstruction inverse problems in which the object of interest is embedded in a heterogeneous background medium that is known only approximately. Shape reconstruction is posed as a stochastic programming problem, in which the background medium is treated as a random field with a known probability distribution. It is demonstrated that by using accelerated stochastic gradient descent the method can be applied to large- scale problems. The capabilities of the method are demonstrated on a simple 2D model problem and in a more demanding application to a 3D inverse conductivity problem in geophysical imaging. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-04-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0377797 |
URI | http://hdl.handle.net/2429/69457 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2019-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2019_may_belliveau_patrick.pdf [ 2.16MB ]
- Metadata
- JSON: 24-1.0377797.json
- JSON-LD: 24-1.0377797-ld.json
- RDF/XML (Pretty): 24-1.0377797-rdf.xml
- RDF/JSON: 24-1.0377797-rdf.json
- Turtle: 24-1.0377797-turtle.txt
- N-Triples: 24-1.0377797-rdf-ntriples.txt
- Original Record: 24-1.0377797-source.json
- Full Text
- 24-1.0377797-fulltext.txt
- Citation
- 24-1.0377797.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0377797/manifest