Induced Polarization Effects inInductive Source ElectromagneticDatabyDavid MarchantB.Sc., The University of Western Ontario, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Geophysics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2015c© David Marchant 2015AbstractInduced polarization (IP) surveys are commonly conducted to map the dis-tribution of electrical chargeability, a diagnostic physical property in mineralexploration and in many environmental problems. Although these surveyshave been successful in the past, the galvanic sources required make theirapplication labour intensive and prevents them from being applied in somegeological settings. The ability to detect chargeability with a geophysicaltechnique that employs inductive sources, eliminating the need to injectcurrent directly into the ground, would provide a valuable tool to appliedgeophysicists.In this work, two aspects of inductive source induced polarization areexamined. First, a new methodology for processing inductive source fre-quency domain EM data to identify IP effects is presented. The methodmakes use of the asymptotic behaviour of the secondary magnetic fields atlow frequency. A new quantity, referred to as the ISIP datum, is defined sothat it equals zero at low frequencies for any frequency-independent (non-chargeable) conductivity distribution. Thus, any non-zero response in theISIP data indicates the presence of chargeable material. Numerical simu-lations demonstrate that the method can be applied even in complicatedgeological situations. A 3D inversion algorithm is developed to recover thechargeability from the ISIP data and the inversion is demonstrated on syn-thetic examples.Understanding the impacts of IP effects on time-domain electromagneticdata requires the ability to simulate common survey techniques while takingchargeability into account. Most existing techniques preform this modellingin the frequency domain prior to transforming their results to the time do-main. Application of those techniques on three dimensional problems caniiAbstractbe computationally limiting. In the second part of this thesis, three newtechniques for forward modelling the time-domain electromagnetic responseof chargeable materials directly in the time domain are developed. Thefirst considers the convolution in the time domain directly, while the othersuse auxiliary differential equations to analytically transform the governingequations into the time domain. The resulting methods are verified bycomparing their results with analytic solutions. The potential applicationof the method was demonstrated by modelling the occurrence of negativetransients in airborne time-domain electromagnetic data.iiiPrefaceThis thesis contains the results of original research preformed while studyingat the Geophysical Inversion Facility at the University of British Columbia.The research has resulted in two publications and four expanded conferenceproceedings. An additional publication containing the remaining materialis in preparation.The idea of looking at scaled combinations of fields came originally fromconversations with Dr. Eldad Haber. This idea forms the basis of the ISIPdatum which is presented in Chapter 2. I carried out the subsequent deriva-tions, numerical tests and manuscript preparation under the supervisionof Dr. Doug Oldenburg. An early version of this work was published inMarchant et al. (2013a), and various parts of were adapted from two con-ference proceedings (Marchant et al., 2012b,c).Dr. Oldenburg encouraged me to look into the inversion of the ISIPdatum. The development of this routine is presented in Chapter 3. I carriedout the derivation of the sensitivities (with some advice from Dr. Haber),wrote all of the code, and produced the synthetic examples. The bulk of thetext in this chapter is adapted from Marchant et al. (2013a).I carried out all of the all of the work involved in the derivation of theconvolution method presented in Chapter 4. I performed all of the program-ming (with some technical advice from Dr. Haber), numerical testing, andmanuscript preparation. This work is currently being prepared for publica-tion.The auxiliary differential equation method presented in Chapter 5 wasan extension of unpublished work started previously by Dr. Laurens Beranand Dr. Oldenburg. I performed all of the programming and testing of thismethod, along with the preparation of the manuscript. An early versionivPrefaceof this work appeared in Marchant et al. (2014) and a pair of conferenceproceedings (Marchant et al., 2012a, 2013b).Dr. Haber provided the suggestion of using a rational function approx-imation to represent the conductivity spectra. This resulted in the Pade´approximation technique present in Chapter 6. I performed the subsequentderivation, implementation and testing of the idea. An early version of thiswork also appeared in Marchant et al. (2014).Appendix A contains material adapted from the course note (MATH522 and EOSC 555B) prepared by Dr. Haber. He also provided the idea ofusing cylindrical meshes for 1D problems which resulted in the material inAppendix B. The reduced model for grounded source experiments (AppendixC) is adapted from a section of Marchant et al. (2014).vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research motivation . . . . . . . . . . . . . . . . . . . . . . . 11.2 Electrical surveying techniques . . . . . . . . . . . . . . . . . 31.3 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . 41.4 Dispersive electrical conductivity . . . . . . . . . . . . . . . . 51.4.1 Chargeability . . . . . . . . . . . . . . . . . . . . . . . 51.4.2 Practical importance of chargeable materials . . . . . 71.4.3 The Cole-Cole model . . . . . . . . . . . . . . . . . . . 81.4.4 Cole-Cole response of earth materials . . . . . . . . . 101.5 The state of the art of chargeability detection . . . . . . . . . 121.5.1 The induced polarization method . . . . . . . . . . . . 12viTable of Contents1.5.2 The magnetic induced polarization method . . . . . . 131.5.3 Difficulties with IP & MIP . . . . . . . . . . . . . . . 141.5.4 IP effects on inductive source electromagnetic surveys 151.6 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . 18I Inductive Source Induced Polarization . . . . . . . . . . . 202 Inductive Source Induced Polarization . . . . . . . . . . . . . 212.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.2 A mathematical model of low-frequency electromagnetics . . 232.3 Inductive source IP . . . . . . . . . . . . . . . . . . . . . . . . 252.4 Synthetic examples . . . . . . . . . . . . . . . . . . . . . . . . 262.4.1 Example #1: Blocks in a uniform half-space . . . . . 262.4.2 Example #2: Blocks in a three-dimensional background 312.5 Sources of error . . . . . . . . . . . . . . . . . . . . . . . . . . 342.5.1 Noisy magnetic fields and error propagation . . . . . . 342.5.2 Low frequency assumption . . . . . . . . . . . . . . . . 362.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393 Inversion of Inductive Source Induced Polarization Data . 423.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.2 Linearization of the ISIP sensitivities . . . . . . . . . . . . . . 433.3 Discretization of the ISIP data . . . . . . . . . . . . . . . . . 443.4 Solving the linear inverse problem . . . . . . . . . . . . . . . 493.4.1 Data misfit . . . . . . . . . . . . . . . . . . . . . . . . 493.4.2 Regularization . . . . . . . . . . . . . . . . . . . . . . 503.4.3 Solving the optimization problem . . . . . . . . . . . . 503.4.4 ISIP specifics . . . . . . . . . . . . . . . . . . . . . . . 523.5 Synthetic Examples . . . . . . . . . . . . . . . . . . . . . . . 533.5.1 Example #1: True sensitivities . . . . . . . . . . . . . 563.5.2 Example #2: Sensitivity approximated with half-space 583.5.3 Example #3: Sensitivity approximated with recov-ered conductivities . . . . . . . . . . . . . . . . . . . . 59viiTable of Contents3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61II Three Dimensional Modelling of IP Effects in TimeDomain Electromagnetic Data . . . . . . . . . . . . . . . . . . 624 Convolutionary Forward Modelling of Time Domain Elec-tromagnetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.2 Convolutionary Ohm’s law . . . . . . . . . . . . . . . . . . . . 644.3 The impulse response of the Cole-Cole model . . . . . . . . . 654.4 Discretization of Ohm’s law in time . . . . . . . . . . . . . . . 704.5 Discretization of Maxwell’s equations . . . . . . . . . . . . . . 754.6 Implementation and validation . . . . . . . . . . . . . . . . . 774.6.1 Comparison to analytic solutions . . . . . . . . . . . . 784.6.2 Three dimensional example . . . . . . . . . . . . . . . 824.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875 Modelling Debye Dispersions Using an Auxiliary Differ-ential Equation Approach . . . . . . . . . . . . . . . . . . . . . 885.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.2 Auxiliary differential equation approach . . . . . . . . . . . . 895.3 Discretization of Maxwell’s equations . . . . . . . . . . . . . . 925.4 Implementation and validation . . . . . . . . . . . . . . . . . 935.4.1 Comparison to 1D analytic solutions . . . . . . . . . . 935.4.2 Comparison to 3D convolution result . . . . . . . . . . 965.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 976 Modelling Cole-Cole Dispersions Using Pade´ Approxima-tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006.2 Rational function approximation of Ohm’s law . . . . . . . . 1016.3 Pade´ approximation of the Cole-Cole model . . . . . . . . . . 1066.3.1 Calculation of Pade´ coefficients . . . . . . . . . . . . . 106viiiTable of Contents6.3.2 Approximation order and central frequency . . . . . . 1086.3.3 The TEM Response of a Pade´ Model . . . . . . . . . . 1126.4 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 1136.5 Implementation and Validation . . . . . . . . . . . . . . . . . 1176.5.1 Comparison to analytic solutions . . . . . . . . . . . . 1176.5.2 Comparison to 3D convolution result . . . . . . . . . . 1216.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1247 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1257.1 Inductive source induced polarization . . . . . . . . . . . . . . 1257.1.1 Future research opportunities . . . . . . . . . . . . . . 1277.2 Three dimensional modelling of IP effects in time domainelectromagnetic data . . . . . . . . . . . . . . . . . . . . . . . 1277.2.1 Future research opportunities . . . . . . . . . . . . . . 129Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132A Finite Volume Modelling of Electromagnetic Methods . . . 145A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145A.2 Discretization in space . . . . . . . . . . . . . . . . . . . . . . 146A.2.1 Fields and physical properties . . . . . . . . . . . . . . 146A.2.2 The curl operator . . . . . . . . . . . . . . . . . . . . . 149A.2.3 Discretization of inner products . . . . . . . . . . . . . 149A.3 Maxwell’s equations in the frequency domain . . . . . . . . . 156A.3.1 Governing equations . . . . . . . . . . . . . . . . . . . 156A.3.2 Weak formulation . . . . . . . . . . . . . . . . . . . . 157A.3.3 The discrete system . . . . . . . . . . . . . . . . . . . 158A.4 Maxwell’s equations in the time Domain . . . . . . . . . . . . 159A.4.1 Governing equations . . . . . . . . . . . . . . . . . . . 159A.4.2 Discretization in time . . . . . . . . . . . . . . . . . . 160A.4.3 Discretization in space . . . . . . . . . . . . . . . . . . 161B Cylindrical Meshes for 1D Problems . . . . . . . . . . . . . . 163B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163ixTable of ContentsB.2 Discretization of fields and physical properties . . . . . . . . . 164B.3 The curl operator . . . . . . . . . . . . . . . . . . . . . . . . . 165B.4 Inner products . . . . . . . . . . . . . . . . . . . . . . . . . . 166B.4.1 Face variables . . . . . . . . . . . . . . . . . . . . . . . 166B.4.2 Edge variables . . . . . . . . . . . . . . . . . . . . . . 168C A Reduced Model for Grounded Source Experiments . . . 170C.1 A reduced model for low frequencies . . . . . . . . . . . . . . 170C.2 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 172C.3 Synthetic example: Effects of EM-coupling on induced polar-ization data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173xList of Tables1.1 Range of measured Cole-Cole parameters from various charge-able deposits as published in Pelton et al. (1978). . . . . . . . 104.1 Execution times for the c = 0.25 1D half-space example usingthe convolution algorithm. . . . . . . . . . . . . . . . . . . . . 824.2 Execution times for the c=0.5 3D example. . . . . . . . . . . 875.1 Execution times for first 1D half-space example (Figure 5.1a)using the ADE algorithm. . . . . . . . . . . . . . . . . . . . . 965.2 Execution times for 3D example using the ADE algorithm. . 986.1 Execution times for the c = 0.25 1D half-space example usingthe Pade´ algorithm. . . . . . . . . . . . . . . . . . . . . . . . 1216.2 Execution times for the c = 0.5 3D example using the Pade´algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122xiList of Figures1.1 Common IP mechanisms . . . . . . . . . . . . . . . . . . . . . 61.2 Influence of η on the Cole-Cole conductivity model . . . . . . 111.3 Influence of τ on the Cole-Cole conductivity model . . . . . . 111.4 Influence of c on the Cole-Cole conductivity model . . . . . . 111.5 Signal waveforms that could be expected in IP surveys. . . . 131.6 Negative transients observed in a VTEM survey . . . . . . . . 172.1 Model and survey geometry for large loop ISIP examples . . . 272.2 Magnetic fields and ISIP data from example #1 using a largeloop survey layout . . . . . . . . . . . . . . . . . . . . . . . . 282.3 Model and survey geometry for Slingram ISIP examples . . . 292.4 Magnetic fields and ISIP data from example #1 using a Slin-gram survey layout. . . . . . . . . . . . . . . . . . . . . . . . 302.5 Conductivity model for example #2. . . . . . . . . . . . . . . 312.6 Magnetic fields and ISIP data from example #2 using a largeloop survey layout. . . . . . . . . . . . . . . . . . . . . . . . . 322.7 Magnetic fields and ISIP data from example #2 using a Slin-gram survey layout. . . . . . . . . . . . . . . . . . . . . . . . 332.8 ISIP data calculated from magnetic fields that have been con-taminated with random noise . . . . . . . . . . . . . . . . . . 362.9 The linear approximation of =(Hz) compared to the true an-alytic response . . . . . . . . . . . . . . . . . . . . . . . . . . 382.10 ISIP data calculated with frequencies with increasing magni-tudes for the large loop example . . . . . . . . . . . . . . . . 40xiiList of Figures2.11 ISIP data calculated with frequencies with increasing magni-tudes for the slingram example . . . . . . . . . . . . . . . . . 413.1 Example of an orthogonal tensor mesh . . . . . . . . . . . . . 453.2 ISIP data compared to ISIP approximations . . . . . . . . . . 483.3 The layout of the survey used in the inversion examples . . . 543.4 True σ0 model used in synthetic inversion examples. . . . . . 553.5 True <[δσ] model used in synthetic inversion examples. . . . 553.6 Synthetic ISIP data used for inversions . . . . . . . . . . . . . 573.7 Recovered <(δσ) model for example #1 . . . . . . . . . . . . 583.8 Recovered <(δσ) model for example #2 . . . . . . . . . . . . 593.9 True and recovered conductivity models . . . . . . . . . . . . 603.10 Recovered <(δσ) model for example #3 . . . . . . . . . . . . 614.1 The impulse response of the Debye model . . . . . . . . . . . 664.2 The impulse response of the Cole-Cole model . . . . . . . . . 674.3 Early time approximation of the impulse response of the Cole-Cole model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.4 True and transformed impulse responses for a Cole-Cole model 694.5 Impulse response of the Cole-Cole model for different valuesof c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.6 Cole-Cole impulse responses compared to the low frequencyapproximation . . . . . . . . . . . . . . . . . . . . . . . . . . 714.7 Cartoon depiction of the terms involved in the evaluation ofEquation 4.18. . . . . . . . . . . . . . . . . . . . . . . . . . . 724.8 Comparison of the magnetic fields simulated using the con-volution algorithm and the transformed analytic expression . 814.9 The layout and core region discretization of the 3D example . 834.10 Vertical component of the calculated b-field data from the 3Dexample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.11 Time decay of the magnetic field for the 3D example . . . . . 854.12 The y-component of current densities for the 3D example . . 86xiiiList of Figures5.1 Magnetic fields simulated using the ADE algorithm and thetransformed analytic expression . . . . . . . . . . . . . . . . . 955.2 Calculated magnetic field data for the 3D example using theADE method . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.3 Decay of the magnetic field calculated using the ADE method 986.1 Pade´ approximations to s0.5 of varying orders . . . . . . . . . 1096.2 True and approximate sc curves for different values of c . . . 1116.3 Cole-Cole dispersions compared to Pade´ approximations fordifferent values of c . . . . . . . . . . . . . . . . . . . . . . . . 1116.4 Cole-Cole dispersions compared to Pade´ approximations fordifferent values of ω0 . . . . . . . . . . . . . . . . . . . . . . . 1126.5 Response of a vertical magnetic dipole at the surface of achargeable half-space with c = 0.75 calculated using the Pade´method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1146.6 Response of a vertical magnetic dipole at the surface of achargeable half-space with c = 0.25 calculated using the Pade´method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1156.7 Comparison of the magnetic fields simulated using the Pade´algorithm and the analytic expression . . . . . . . . . . . . . 1206.8 Calculated magnetic field data for the 3D example using thePade´ method . . . . . . . . . . . . . . . . . . . . . . . . . . . 1226.9 Decay of the magnetic field calculated using the Pade´ method 123A.1 An example of a non-uniform orthogonal mesh. . . . . . . . . 146A.2 Discretization of a cell . . . . . . . . . . . . . . . . . . . . . . 147A.3 Indexing convention of the mesh . . . . . . . . . . . . . . . . 148A.4 Quantities required to calculate the discrete approximationof the curl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150A.5 Sparsity structure of the discrete curl operator . . . . . . . . 151A.6 Field components involved in discretizing an inner productover a single cell . . . . . . . . . . . . . . . . . . . . . . . . . 152A.7 Sparsity patterns of the averaging matrices . . . . . . . . . . 154xivList of FiguresB.1 An example of a cylindrical mesh. . . . . . . . . . . . . . . . 163B.2 Cylindrical cell edges and face locations . . . . . . . . . . . . 164B.3 Elements involved in the approximation of the curl on a cylin-drical mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165B.4 Location of face elements for cylindrical meshes . . . . . . . . 166B.5 Location of edge elements for cylindrical meshes . . . . . . . 168C.1 Layout of synthetic gradient array experiments . . . . . . . . 174C.2 The x-component of the electric fields . . . . . . . . . . . . . 176C.3 Time evolution of the x-component of the electric fields . . . 177xvList of Algorithms4.1 Convolution based forward modelling . . . . . . . . . . . . . . 795.1 ADE based forward modelling . . . . . . . . . . . . . . . . . . 946.1 Pade´ forward modelling . . . . . . . . . . . . . . . . . . . . . 118xviAcknowledgementsThis thesis would never have been possible without the guidance, assistanceand support of a large group individuals.First and foremost, I would like to express my deepest gratitude to mysupervisor Prof. Doug Oldenburg for all of the guidance he has provided.Thank you for opening my eyes to the world of applied geophysics andinversion. Working with you has provided me with an understanding ofgeophysics that extends far beyond the scope of this thesis and for which Iam incredibly grateful.Next, I am heavily indebted to my committee member Prof. EldadHaber. Without the foundation in computational electromagnetics you gaveme, as well as the encouragement and guidance you provided, much of thework in this document would never have been completed.The funding of this work was generously provided from a variety ofsources. Thanks to University of British Columbia, the Department of Earthand Ocean Sciences, the Society of Exploration Geophysics and the NaturalSciences and Engineering Research Council of Canada for their support.Many thanks to all of my friends and colleagues at the Geophysical Inver-sion Facility and Computational Geosciences that I have shared time withover the years. Thanks to Rob Eso, Laurens Beran, Elliot Holtham, GregNash, Dikun Yang, Justin Granek, Gudni Rosenkjar, Livia Mahler, SarahDevriese, Kris Davis, Daniel Bild-Enkin, Mike McMillan, Mike Mitchell,Rowan Cockett, Seogi Kang, and Lindsey Heagy (to name a few) for mak-ing my time here so enjoyable. A very special thanks to Roman Shekhtmanfor all of your assistance with all things code and computer related.Finally, I would like to express my love and gratitude to Jenn Fohringfor all of the support and encouragement she provided throughout the highsand lows of this process.xviiTo my parents.xviiiChapter 1Introduction1.1 Research motivationMost of the world’s remaining undiscovered natural resources are buriedbeneath the surface of the Earth, and are thus hidden from traditional sur-face exploration methods. Locating these resources without the use of costlydrilling programs requires the application of geophysical imaging techniques.The data obtained from the application of these techniques can be inter-preted to gain a better understanding of the distribution of the physicalproperties of materials in the ground. For example, the ability to detectthe chargeability of a material, a physical property characterized by a fre-quency dependent electrical conductivity, has found numerous applicationsin applied geophysics and is of particular importance when exploring fordisseminated metallic mineralization.Existing techniques for detecting chargeable material rely on a groundedgalvanic source to generate an induced polarization response. While thesetechniques have been applied successfully for years, it would be beneficialto use inductive transmitters in place of the galvanic sources. The use ofinductive sources removes the necessity to inject current directly into theground which is limited in environments with extremely resistive overbur-dens. An inductive source IP method could also potentially lead to a sig-nificant decrease in the cost and effort involved in conducting IP surveys.The difficulty associated with the application of inductive source electromag-netic techniques arises when attempting to identify the induced polarizationeffects in the data.Electromagnetic experiments can be carried out either in the time or thefrequency domain. In the frequency domain, the magnitude and phase of11.1. Research motivationelectric or magnetic fields are recorded as a function of transmitter frequency.There are some situations where induced polarization effects are easilyidentifiable in time domain experiments. Notably, the presence of a signreversal in coincident loop time domain data can be attributed to chargeablematerial. For other configurations of sources and receivers there may beno obvious manifestation of induced polarization effects. A procedure foridentifying the presence of induced polarization in frequency domain dataremains unknown.Simulating the response of a geophysical experiment is an invaluable toolwhen trying to understand and interpret observations. Simulating the elec-tromagnetic response of chargeable materials in the frequency domain is, inprinciple, straight forward. Modelling the response of a frequency dependentconductivity in the time domain is computationally challenging. Existingtechniques rely on frequency to time domain transforms. These techniqueswork well when working on small one or two dimensional problems, butbecome computationally limiting when extended to 3D.In this thesis I consider two open questions in the area of inductive sourceinduced polarization.1. Can the effects of chargeability be identified in frequency domain in-ductive source electromagnetic data?2. Can we efficiently simulate the electromagnetic response of chargeablematerials in the time domain?Part 1 of this thesis addresses the first question. In it, I propose anacquisition and processing technique that allows the identification of inducedpolarization effects in frequency domain electromagnetic data.In Part 2 I consider the second question. Here, a suite of new forwardmodelling techniques is developed to simulate the time domain response ofthree dimensional distributions of chargeable material.The remainder of this introductory chapter is devoted to a review ofelectrical geophysical techniques, the causes of chargeability, and the stateof the art in chargeability detection. I conclude this chapter with an outlineof the thesis.21.2. Electrical surveying techniques1.2 Electrical surveying techniquesElectrical methods are commonly applied to detect the distribution andstructure of the electrical conductivity of the subsurface. All electrical meth-ods measure the electromagnetic response of the Earth to external stimula-tion, but the way in which this is achieved varies greatly.The Earth can be illuminated using either galvanic or inductive trans-mitters. Galvanic transmitters operate by injecting electrical current di-rectly into the ground through pairs of grounded electrodes. The currentcan be injected at a constant rate to carry out direct current resistivity (DCresistivity) (Kunetz, 1966) or magnetometric resistivity (MMR) (Edwardset al., 1978) experiments, or varied as a function of time as is the case incontrolled source audio-frequency magnetotellurics (CSAMT) (Zonge andHughes, 1991).Alternatively, the Earth can be illuminated inductively using a timevarying magnetic field. One can make use of naturally occurring fields,as is done in the magnetotelluric (MT) (Simpson and Bahr, 2005) or Z-axis tipper electromagnetic (ZTEM) (Lo and Zang, 2008) techniques, or themagnetic field can be created by passing a time varying current though a wireloop. This final class of experiments are commonly referred to as controlledsource electromagnetics (CSEM) (Nabighian, 1988). CSEM experiments areconducted with either the transmitter loop placed on the ground, or towedfrom a moving platform such as a helicopter.Similarly, data can be recorded in a variety of ways. For instance, thedifference in the electrical potential between two points can be measuredby recording the voltage across a pair of grounded electrodes, or the timederivative of the magnetic flux density can be sampled by observing thevoltage induced in wire loop. The magnitude or the components of themagnetic field itself could also be measured using a magnetometer.Regardless of the transmitter or type of observation, it is possible toinfer details of the distribution of conductive and resistive materials in thesubsurface from the resulting information. The understanding of the con-ductivities that are provided by these methods can often be linked to other31.3. Maxwell’s equationsproperties of interest, such as the location and extent of various rock units,the depth to the water table, or the location of a hydrocarbon reservoir.1.3 Maxwell’s equationsAll electromagnetic processes are described by Maxwell’s equations. Takingan eiωt convention for the Fourier transform, the frequency domain form ofMaxwell’s equations are~∇× ~E + iω ~B = 0 (1.1a)~∇× ~H + iω ~E = ~J + ~Js (1.1b)∇ · ~E = ρ (1.1c)∇ · ~B = 0 (1.1d)where ~E is the electric field, is the dielectric permittivity, ~H is the magneticfield, ~J is the current density, ~Js is the current density of the source, and ~Bis the magnetic flux density.The fields and fluxes are related through the constitutive relationships~B = µ ~H (1.2a)~J = σ ~E (1.2b)where µ is the magnetic permeability and σ is the electrical conductivity.The electrical resistivity ρ is often considered instead of the conductivity.The two are related by ρ = 1/σ.Most geophysical electromagnetic surveys operate at relatively low fre-quency (< 105Hz) and usually involve materials that are relatively conduc-tive (> 10−5S/m). As the permittivity of most earth materials is quite small( 10−12F/m) it is common to neglect the displacement currents present inAmpe´re’s law when analyzing the results of electromagnetic surveys (Wardand Hohmann, 1988). To this end, it is assumed that iω ~E ≈ 0 allowing usto work with the simplified magnetoquasistatic form of Maxwell’s equations41.4. Dispersive electrical conductivity~∇× ~E + iω ~B = 0 (1.3a)~∇× ~H = ~J + ~Js (1.3b)∇ · ~B = 0 (1.3c)Transforming these equations to the time domain results in the followingsystem of equations~∇× ~e+ ∂~b∂t= 0 (1.4a)~∇× ~h = ~j +~js (1.4b)∇ · ~b = 0 (1.4c)1.4 Dispersive electrical conductivityThe physical properties appearing in Equations 1.2a and 1.2b are not nec-essarily constant with respect to frequency. Instead, in some circumstances,the physical properties depend on the frequency at which it is being mea-sured.1.4.1 ChargeabilityWhen the frequency dependence, or dispersion, of the electrical conductiv-ity of a material is significant, the material is considered to be chargeable.Chargeable materials act like capacitors and become electrically polarizedwhen exposed to an electric field.In general, the conductivity of a chargeable material will increase withincreasing frequency. Causality conditions require that any change in theamplitude of a materials properties be accompanied by a corresponding non-zero phase in accordance to the Kramers-Kronig relations (Fuller and Ward,1970). Thus, a frequency dependent conductivity will have both real andimaginary parts.51.4. Dispersive electrical conductivityThe causes of chargeability are complex and vary depending on the typeof material, but all share some common characteristics (Everett, 2013). Asan electric current flows through the chargeable material, electrochemicaleffects cause charges to accumulate non-uniformly within the media, mak-ing it become polarized. The polarization then impedes the flow of currentthrough the material. When the imposed current is switched off, the accu-mulated charges relax back to their original state. This results in a net flowof current in the opposite direction of the original current flow.There are three simplified models that are commonly discussed in theliterature: polarization by constriction of pores, membrane polarization, andelectrode polarization (Scho¨n, 2011; Snyder et al., 1977).(a) (b) (c)Figure 1.1: Three common mechanisms of IP (adapted from Scho¨n (2011));(a) Polarization by constriction of pores (b) Membrane Polarization and (c)Electrode PolarizationConceptually, the simplest of these mechanisms is constrictivity polar-ization (Bo¨rner, 2006; Titov et al., 2004). The bulk of current flow throughporous media is electrolytic, that is, the movement of negative and positiveions in groundwater. A constriction in the pore space of the backgroundmedia can sometimes result in the variable mobility of anions and cations.This results in the build up of charge at the boundaries of these constric-tions (Figure 1.1a). This mechanism is the primary source of a frequency61.4. Dispersive electrical conductivitydependent conductivity in media free of clay or metallic material.Membrane polarization (Marshall and Madden, 1959; Vinegar and Wax-man, 1984) occurs when a charge builds up in non-metallic mineralization(such as clay) partially blocking a pore (Figure 1.1b). In the equilibriumstate, the negatively charged clay causes positively charged cations to accu-mulate in the pore surrounding the clay particle. When an external currentis applied this concentration of cations selectively impedes the ability ofother ions to move through the pore space. The result is a net buildup ofcharge at either end of the pore.Finally, electrode polarization (Marshall and Madden, 1959; Wong, 1979)occurs at the boundary of an electrolytic solution (an ionic conductor, suchas groundwater) and a conductive material (Figure 1.1c). This mechanismis also commonly called over-voltage. In the presence of an external current,electrolysis occurs at the metal electrolyte boundary, with electrons beinggained by the conductor at one boundary and lost at the other. Currentflows faster within the conductive material than in the electrolyte, resultingin accumulations of ions in the electrolyte at the boundary of the conductor.1.4.2 Practical importance of chargeable materialsAs a result of the over-voltage effect, the presence of polarizable materialsin the ground is an excellent proxy for the distribution of metallic minerals.Thus, surveys sensitive to the presence of chargeable materials are routinelyused in mineral exploration. These surveys are of particular value when ex-ploring for disseminated mineralization where chargeability is commonly themost diagnostic physical property. Numerous examples of this applicationcan be found in and Pelton et al. (1978) and Fink et al. (1990).While mineral exploration was the original focus, recent work has ex-panded to include a wide variety of other applications including the obser-vation of IP effects to hydrocarbon exploration (Burtman et al., 2014; Con-nell and Key, 2013; Davydycheva et al., 2006; Heenan et al., 2013; Veekenet al., 2012, 2009a), environmental and groundwater studies (Ho¨rdt et al.,2007; Kozhevnikov, 2012; Revil et al., 2012; Slater and Glaser, 2003), and71.4. Dispersive electrical conductivitynumerous other environmental and engineering applications (Beran and Old-enburg, 2008; Kemna et al., 2012, 2004; Okay et al., 2013; Weller et al., 2000;Yuval and Oldenburg, 1996).1.4.3 The Cole-Cole modelThere are a number of parametric models that are used to describe thenature of the frequency dependence of a material’s electrical conductivity.An overview and comparison of the different models is provided by Dias(2000). In the geophysics community, the most commonly used model is theCole-Cole model. The model was first proposed to fit observations of thedispersion spectra of dielectric permittivity (Cole and Cole, 1941)(ω) = +0 − ∞1 + (iωτ)c(1.5)It was later observed that the same model could also be used to fit observeddispersion in electrical resistivities (Pelton et al., 1978). Substituting ρ for in 1.5 and definingη =ρ0 − ρ∞ρ0(1.6)gives an expression for a frequency dependent resistivityρ(ω) = ρ0(1− η(1− 11 + (iωτ)c))(1.7)This expression was referred to by Pelton et al. (1978) as the Cole-Colemodel.A recent paper (Tarasov and Titov, 2013) suggested that a substitutionfor σ for would have been more appropriate. However, this approach doesnot impact the model’s ability to fit experimental data. Due to the extensivevolume of literature existing using Pelton’s form of the Cole-Cole equationthe work in this thesis will focus on materials whose frequency dependenceis described by Equation 1.7.In Equation 1.7, ρ0 is the electrical resistivity measured at zero fre-81.4. Dispersive electrical conductivityquency. The parameter η is commonly referred to as the chargeability, τ isa characteristic time constant, and c is the frequency dependence. In thehigh frequency limit Equation 1.7 reduces toρ∞ = (1− η)ρ0 (1.8)Taking the reciprocal of Equation 1.7 results in the expression for theCole-Cole conductivity model in terms of σ0σ(ω) = σ0(1 + η((iωτ)c1 + (1− η)(iωτ)c))(1.9)Making use of the high frequency limit of this expressionσ0 = (1− η)σ∞ (1.10)results in a Cole-Cole conductivity model in terms of σ∞σ(ω) = σ∞(1− η1 + (1− η)(iωτ)c)(1.11)The parameter η is the fractional decrease in the conductivity betweenthe high and low frequency asymptotes given byη = 1− σ0σ∞(1.12)Typically, η can take values between 0 and 1 with larger numbers corre-sponding to a larger difference between σ0 and σ∞ (Figure 1.2).The time constant τ influences the frequency at which the frequencydependence of the real part of the conductivity is maximized (Figure 1.3).This also corresponds to the maximum magnitude of the imaginary part ofthe conductivity.The frequency dependence, c, ranges from 0 to 1 and determines howrapidly the transition of the real part of the resistivity takes place (Figure1.4). A low value of c results in a very broad dispersion curve, whereas ahigh c value transitions rapidly from the low frequency to high frequency91.4. Dispersive electrical conductivityasymptotes. In the special case when c is equal to 1, the Cole-Cole modelreduces to the Debye model (Debye and Huckel, 1923). In this case, thereis a very rapid transition from the low to the high frequency asymptote. Ithas been observed that most chargeable materials of economic interest tothe mining community have a frequency dependence between 0.1 and 0.5(Pelton, 1977; Wong, 1979).1.4.4 Cole-Cole response of earth materialsTo gain an understanding of the Cole-Cole response of earth materials Pel-ton et al. (1978) performed in-situ measurements of the conductivity spectraof common occurrences of chargeable mineralization. The goal of the studywas to determine the range of values for the Cole-Cole parameters typical ineach deposit type with the ultimate goal of providing a means to discrimi-nate between economic and non-economic IP responses. The range of valuesdetermined in this study is summarized in table 1.1. The study concludedthat the properties of the dispersion related more to the concentration andthe grain size of the sulphides present, and fairly little to do with the com-position of the sulphides. The exception to this was the comparison betweenmassive sulphides and graphite, where the range in observed parameters wasconcluded to be such that rudimentary classification may be possible.ρ0 (Ωm) η τ (s) cPorphyry (dry) 1× 101 - 1× 103 0.1 - 0.5 1× 10−5 - 1× 100 0.1 - 0.6Porphyry (wet) 1× 100 - 1× 104 0.1 - 0.8 1× 10−2 - 7× 104 0.1 - 0.5Magnetite 1× 101 - 1× 103 0.1 - 1.0 8× 10−4 - 3× 100 0.1 - 0.6Pyrrhotite 1× 100 - 1× 103 0.3 - 0.8 3× 100 - 1× 105 0.1 - 0.5Massive Sulfide 1× 10−2 - 1× 103 0.6 - 0.95 8× 10−4 - 2× 100 0.1 - 0.4Graphite 1× 10−2 - 1× 103 0.75 - 0.98 3× 101 - 8× 103 0.1 - 0.5Table 1.1: Range of measured Cole-Cole parameters from various chargeabledeposits as published in Pelton et al. (1978).101.4. Dispersive electrical conductivity10-210-1100101102103ω (rad/s)0.0250.0500.0750.100ℜ(σ) (S/m)σ∞=0.1S/mτ=1.0sc=0.5η=0.25 η=0.50η=0.75(a)10-210-1100101102103ω (rad/s)0.000.010.020.03ℑ(σ) (S/m)σ∞=0.1S/mτ=1.0sc=0.5η=0.25 η=0.50η=0.75(b)Figure 1.2: Influence of η on the (a) real and (b) imaginary part of theCole-Cole conductivity model.10-210-1100101102103ω (rad/s)0.050.060.070.080.090.10ℜ(σ) (S/m)σ∞=0.1S/mη=0.5c=0.5τ=0.1 τ=1.0 τ=10.0(a)10-210-1100101102103ω (rad/s)0.0000.0020.0040.0060.0080.0100.0120.0140.0160.018ℑ(σ) (S/m)σ∞=0.1S/mη=0.5c=0.5τ=0.1 τ=1.0 τ=10.0(b)Figure 1.3: Influence of τ on the (a) real and (b) imaginary part of theCole-Cole conductivity model.10-210-1100101102103ω (rad/s)0.050.060.070.080.090.10ℜ(σ) (S/m)σ∞=0.1S/mη=0.5τ=1.0sc=0.25 c=0.50c=0.75c=1.0(a)10-210-1100101102103ω (rad/s)0.0000.0050.0100.0150.0200.025ℑ(σ) (S/m)σ∞=0.1S/mη=0.5τ=1.0sc=0.25 c=0.50c=0.75c=1.0(b)Figure 1.4: Influence of c on the (a) real and (b) imaginary part of theCole-Cole conductivity model.111.5. The state of the art of chargeability detection1.5 The state of the art of chargeability detectionTraditionally, there are two geophysical techniques that are applied to mapthe distribution of chargeable material: the induced polarization (IP) tech-nique (Bleil, 1953; Seigel, 1959) and the magnetic induced polarization(MIP) technique (Seigel, 1974).1.5.1 The induced polarization methodThe induced polarization method is an extension of the DC resistivity methodof mapping the electrical conductivity of the ground. In the DC experiment,a steady state current is injected into the ground through a pair of trans-mitter electrodes. The resulting differences in the electrical potential arethen measured across pairs of receivers. These data are easily inverted torecover the distribution of electrical conductivity in the ground (Ellis andOldenburg, 1994; Li and Oldenburg, 1994).The induced polarization technique makes use of the same equipmentas the DC resistivity experiment. Rather than measuring the steady statepotential differences, voltages are measured after the transmitter currenthas been switched off. The presence of a slowly decaying potential (in theabsence of significant electromagnetic induction) is indicative of the presenceof chargeable materials (Figure 1.5a).Alternatively, measurements can be taken in the frequency domain. Whenoperating at low enough frequencies (where negligible electromagnetic induc-tion occurs), IP effects are recognized as either variations in the magnitudeof the received voltage with respect to frequency, or as a phase shift betweenthe transmitted and received signal (Figure 1.5b).Both time domain and frequency domain IP measurements are com-monly inverted in two or three dimensions to recover models of the distri-bution of chargeable material (Li and Oldenburg, 2000; Oldenburg and Li,1994).A common extension of the frequency domain technique is the spectralinduced polarization technique (also referred to as the complex resistivitymethod) (Kemna et al., 2012). This technique attempts to detect the nature121.5. The state of the art of chargeability detectiontVTransmittedReceived(a) Time domain experimenttVTransmittedReceived(b) Frequency domain experimentFigure 1.5: Signal waveforms that could be expected in (a) time or (b)frequency domain induced polarization surveys.of the frequency dependence over a wide band of frequencies. Measuring theproperties of the dispersion curve could allow for discrimination betweendifferent chargeable materials (Pelton et al., 1978).An excellent summary of the historical development of the IP methodwas given by Seigel et al. (2007). The method has a proven track recordin mineral exploration and is widely considered the geophysical method ofchoice when looking for porphyry copper deposits (Fink et al., 1990).1.5.2 The magnetic induced polarization methodThe magnetic induced polarization (MIP) method was proposed by Seigel(1974) to make use of observations of magnetic fields rather than electricalpotentials. Just as the IP technique is closely related to the DC resistivitytechnique, the MIP technique is related to the magnetometric resistivity(MMR) technique.As in the DC resistivity experiment, a steady current is injected into the131.5. The state of the art of chargeability detectionground across a pair of transmitter electrodes. The magnetic fields producedfrom the currents flowing through the ground are measured and interpretedto recover the distribution of current flow and thus the distribution of con-ductivities (Chen et al., 2002).For MIP, the magnetic fields are observed after the current is interrupted.When electromagnetic coupling is neglected, these fields are directly at-tributed to the current flow in the ground as the accumulated charges relaxback to their equilibrium state.As magnetic fields are measured instead of electric fields, the time con-suming requirement of placing receiver electrodes is eliminated. The MIPtechnique also provides improved performance when operating where highlyconductive overburden exists. A 3D inversion technique for MIP data wasdeveloped by Chen and Oldenburg (2003).1.5.3 Difficulties with IP & MIPDespite the induced polarization technique’s widespread and successful ap-plication in the mining industry, it has its shortcomings. The IP techniquerequires the time consuming task of planting numerous electrodes into theground, often making the survey prohibitively expensive and preventing itsapplication to large scale reconnaissance. All of the electrodes must makesufficient electrical contact with the ground, and it must be possible to pushenough current into the ground through the transmit electrodes to excite anIP response. It is not always possible to meet these requirements in somegeologic environments, such as those with very resistive overburden. Someof these problems were solved with the development of the MIP technique.However it is still necessary to transmit current into the ground, and nu-merous transmitter pairs are required to gain a 3D understanding of thechargeability distribution.Another major problem with both the IP and the MIP experiments is thepresence of EM coupling (Dey and Morrison, 1973). Both techniques assumethat it is possible to work at late enough time (or low enough frequency)that the time derivative of the magnetic fields associated with the switching141.5. The state of the art of chargeability detectionoff of the transmit current can be neglected. While this is a valid assumptionin more resistive environments, it can render the IP data practically uselessin more conductive environments. A number of different methods have beenproposed to ’decouple’ the IP and EM portions of the response. Theseapproaches either try to design IP surveys in such a way as to minimize thepresence of EM effects (Coggon, 1973; Schmutz et al., 2014; White et al.,2003) or try to separate and the EM effects from the observed data (Coggon,1984; Dey and Morrison, 1973; Hallof, 1974; Hohmann, 1973; Routh andOldenburg, 2001). While these techniques can improve the resulting datain some situations, they are all approximations and will fail under certaincircumstances.1.5.4 IP effects on inductive source electromagnetic surveysInductive source electromagnetic techniques are also extensively used to mapthe distribution of conductivities (Nabighian, 1988). EM methods have anumber of advantages over the DC resistivity and MMR techniques. Thetransmitter operates by driving a current through a loop of wire, so it is notnecessary to force a current into the ground through electrodes. Observa-tions are typically made of the magnetic fields, or of the time derivative ofthe magnetic fields, which also requires no contact with the ground. Theseobservations are readily taken from moving platforms allowing for the rapidacquisition of large amounts of data.Electromagnetic surveys are also affected by chargeable materials. Un-fortunately, for most survey geometries, the effects are often impossible torecognize in the data. In these situations, the resulting IP effects are essen-tially noise, which hinder the interpretation of the data.For the particular case of coincident loop time-domain EM data, neg-ative transients - soundings with a reversal in sign of the received field -are diagnostic of chargeable material. While early papers speculated thatthis effect could be caused by particular conductivity distributions or mag-netic effects, Weidelt (1982) showed that for a coincident loop system witha step-off primary field, the measured secondary field must be non-negative151.5. The state of the art of chargeability detectionover non-polarizable ground regardless of the subsurface distribution of con-ductivity. Observed negative transients in coincident loop EM systems cantherefore only be attributed to polarization effects. In many cases, thisproperty can be extended to center loop systems, including many airborneplatforms Smith and Klein (1996). Negative transients are commonly ob-served in airborne TEM systems, such as AeroTEM or VTEM. An example,taken from a VTEM survey performed over the Mt. Milligan deposit inBritish Columbia, is shown in Figure 1.6.Despite the fact that negative transients can be related to the presenceof chargeable material, relatively little has been done to try and interpretthem directly. Smith et al. (2008) used the presence of negative transientsto map the presence of tailings around a mine site. Beran and Oldenburg(2008) applied an over-determined inversion routine to fit IP affected TEMdata to a simple layered model exhibiting Cole-Cole dispersion properties.Kozhevnikov and Antonov (2010) applied a similar methodology, solvingthe over-determined problem to recover Cole-Cole parameters of either auniform halfspace, or a uniform two-layer model.Compared to the large amount of literature dealing with negative tran-sients, very little has been published on effects of dispersive conductivitieson inductive source frequency domain electromagnetic data. One notableexception is the work presented in Hohmann et al. (1970, 1971). The fre-quency domain response of a two layer, chargeable earth was simulated for arange of model parameters at frequencies between 10 Hz and 1500 Hz. Theresults showed that the presence of a layer with a frequency dependent con-ductivity did have an effect on the data, but that the effects were very small.They went on to show that the effects could be explained by a different non-chargeable layered model. Field data from three test sites with known IPanomalies were also considered. The field tests only considered the ampli-tude response and all were carried out in relatively conductive environments.All of the observations could be explained by reasonable earth models thatcontained no IP effects and thus it was concluded that frequency domaininductive source experiments were not suitable for detecting the presence ofchargeable material under the conditions that were tested.161.5. The state of the art of chargeability detection0 200 400 600 800 1000 1200 1400 1600Easting (m)−10−50 10−510−410−310−2nTs(a)10-410-310-2Time (s)10-610-510-410-310-2nTs(b)Figure 1.6: Example of a negative transient observed in a VTEM survey atthe Mt. Milligan deposit in British Columbia (a) Line of data above thedeposit. (b) Single sounding at 850m along the line. Dashed lines indicatea negative value.171.6. Outline of the thesis1.6 Outline of the thesisThe content of this thesis is divided into two parts. Part 1 (Chapters 2and 3) considers the detection and interpretation of IP effects on inductivesource frequency domain electromagnetic data. In Chapter 2 I propose anew data processing technique for processing frequency domain EM data.This technique identifies the presence of IP effects in observations of themagnetic fields arising from an inductive source. The method makes use ofthe asymptotic behaviour of the secondary magnetic fields at low frequency.A new quantity, referred to as the ISIP datum, is defined so that it equalszero at low frequencies for any frequency-independent (non-chargeable) con-ductivity distribution.Chapter 3 considers the interpretation of this new type of data. Themodelling of the ISIP data is linearized, and the resulting linearization isemployed to develop a three dimensional inversion routine. This routine canbe used to recover a chargeability distribution model from observations ofISIP data.In Part 2 I address the problem of simulating the response of chargeablematerials directly in the time domain. Three different methods are presentedin Chapters 4, 5 and 6.Directly transforming Ohm’s law with a frequency dependent conductiv-ity to the time domain results in a convolution. Modelling this convolutiondirectly for three dimensional distributions of chargeable material is pre-sented in Chapter 4. This approach is shown to produce accurate resultsand is demonstrated with synthetic examples.While the convolution based approach produces accurate results, theevaluation of the convolution can require significant computational resources.In Chapter 5 I develop an additional modelling algorithm for Debye dis-persions that makes use of an auxiliary differential equation to avoid theevaluation of the convolution. This approach significantly reduces the com-putational expense.Chapter 6 uses an approximation to the frequency dependant conductiv-ity to derive an auxiliary differential equation similar to that used in Chapter181.6. Outline of the thesis5. While this approach does result in a slightly less accurate simulation thantreating the full convolution, it significantly decreases the required compu-tational resources while allowing any Cole-Cole conductivity to be modelled.Finally, Chapter 7 contains a discussion of the thesis contributions anddirections for future research.This material is accompanied by three appendices. The first two of theseprovide an overview of the finite volume modelling techniques employedthroughout this thesis. The final appendix presents a reduced model of theforward modelling routine developed in Chapter 5 that can be used to modelgrounded source experiments in the absence significant EM coupling.19Part IInductive Source InducedPolarization20Chapter 2Inductive Source InducedPolarization2.1 IntroductionCompared to the large body of literature dealing with negative transientsand induced polarization effects in time domain electromagnetic data, verylittle has been published on the effects of dispersive conductivities on induc-tive source frequency domain data.One notable exception is the work presented in Hohmann et al. (1970)and Hohmann et al. (1971), where the synthetic frequency domain responseof a chargeable two-layer model was computed for a variety of model pa-rameters at frequencies ranging from 10 Hz to 1500 Hz. The results showedthat the presence of a layer with a frequency dependent conductivity hasan effect on the data, but that those effects were very small. They went onto show that these effects could be explained by a different non-chargeablelayered model. Field data from two test sites with known IP anomalies werealso examined. The field tests only considered the amplitude response andall were carried out in relatively conductive environments. All observationscould be explained by reasonable earth models that contained no IP effects.It was concluded that frequency domain inductive source experiments werenot suitable for detecting the presence of chargeable material under the con-ditions of the tests. Furthermore, additional tests would be needed to assessthe potential of the method in high-resistivity environments and considermore than just amplitude data.In another study, Gasperikova (1999) and Gasperikova and Morrison212.1. Introduction(2001) considered mapping induced polarization using natural field electro-magnetic data. They proposed that any change in the response of a bodyat frequencies low enough to be in the DC limit must be indicative of afrequency dependent conductivity. This technique was shown to work onsimple synthetic examples, and then applied to field data sets collected overknown IP anomalies. However, it was observed that their technique couldnot work in a layered media as no DC limit would be reached and the pres-ence of a frequency-dependent conductivity could not be separated from thefrequency dependent EM response. Luo et al. (2003) examined the proposalof natural source induced polarization and concluded that, as the structureof the earth is generally layered, this non-layered earth requirement wouldnever be satisfied in reality. They also went on to note that, at the frequen-cies specified by Gasperikova and Morrison (2001), the IP effects would betoo small to extract for the observed fields.In this chapter, a new data collection and processing methodology isproposed which detects the presence of chargeable materials using induc-tive sources and observations of magnetic fields in the frequency domain.The asymptotic behaviour of the imaginary part of the magnetic fields atlow frequency is exploited to provide a direct indication of the presence ofchargeability. As long as the technique is performed at low enough frequencyto satisfy a low frequency assumption, the observed data will indicate thepresence of chargeable material regardless of the distribution of conductivi-ties.Section 2.2 examines the mathematics of electromagnetics at low fre-quency, and leads to the definition of the inductive source induced polariza-tion data in section 2.3. The idea is demonstrated on a series of syntheticexamples in section 2.4. Finally, sources of error and the feasibility of theidea are discussed in section 2.5.222.2. A mathematical model of low-frequency electromagnetics2.2 A mathematical model of low-frequencyelectromagneticsWith an e−iωt time dependence, Maxwell’s equations in the frequency do-main are~∇× ~E + iω ~B = 0 (2.1a)~∇× 1µ0~B − σ ~E = ~J0 (2.1b)where J0 is defined to be the current density flowing in a transmitter loop.The magnetic flux density ~B can be decomposed into primary ( ~B0) andsecondary ( ~Bs) parts such that~B = ~Bs + ~B0 (2.2)~B0 is defined such that it is equal to the magnetic flux density generated bythe transmitter loop operating at zero frequency~∇× 1µ0~B0 = ~J0 (2.3)Writing Equation 2.1 in terms of Bs and B0~∇× ~E + iω ~Bs = −iω ~B0 (2.4a)~∇× 1µ0~Bs − σ ~E = 0 (2.4b)and eliminating ~E gives~∇× ρ~∇× 1µ0~Bs + iω ~Bs = −iω ~B0 (2.5)where ρ is the electrical resistivity (ρ = 1σ ).Assuming a non-dispersive (frequency independent or non-chargeable)resistivity distribution and differentiating equation 2.5 with respect to the232.2. A mathematical model of low-frequency electromagneticsangular frequency, ω, yields~∇× ρ~∇× 1µ0∂ ~Bs∂ω+ iω∂ ~Bs∂ω= −i(~B0 + ~Bs)(2.6)Equation 2.6 is a partial differential equation for ∂~Bs∂ω . It has the same formas Equation 2.5 but has a different right hand side.As ~Bs is zero when ω = 0 by definition, at low enough frequency thesecondary magnetic field can be approximated by the first order Taylor series~Bs(w) ≈ ∂~Bs∂ω∣∣∣∣ω=0ω (2.7)Taking the limit of equation 2.6 as ω → 0 leaves~∇× ρ~∇× 1µ0∂ ~Bs∂ω= −i ~B0 (2.8)where ρ is assumed to be non-dispersive so that the operator ~∇× ρ~∇× ispurely real. Since ~B0 is also a real quantity by definition, the right handside of Equation 2.8 is therefore purely imaginary. Thus, for small ω thereal part of the first derivative of the magnetic field is equal to zero<[∂ ~Bs∂ω∣∣∣∣ω=0]≈ 0 (2.9)and the first order Taylor series defined in 2.7 is purely imaginary=[~Bs(ω)]≈ ∂~Bs∂ω∣∣∣∣ω=0ω (2.10)Therefore, in an environment with frequency independent conductivities theimaginary part of the magnetic fields will vary linearly with frequency.242.3. Inductive source IP2.3 Inductive source IPSo far it has been shown that the behaviour of the imaginary part of themagnetic field is easily predicted at low frequencies when no chargeablematerial is present. This property will now be used to detect chargeablematerial.Consider the magnetic response of a non-chargeable earth to forcingfrom an inductive source operating at two closely spaced, low frequencies,ω1 and ω2. If the frequencies are sufficiently low then the skin-depth willbe very large compared to the geometric decay of the source fields, and themeasured response will be sensitive to the same volume of earth, regardlessof frequency. From equation 2.10 it can be said that=[~Bs(ω)]ω≈ ∂~Bs∂ω∣∣∣∣ω=0(2.11)and that this quantity is constant with respect to frequency. The differencein this quantity measured at ω1 and ω2 would approximately equal zero.=[~Bs(ω1)]ω1−=[~Bs(ω2)]ω2≈ 0 (2.12)This observation motivates the introduction of a new quantity, dISIPwhich is defined as the inductive source induced polarization (ISIP) datadISIP = =[ ~Bs(ω2)]− ω2ω1=[ ~Bs(ω1)] (2.13)The quantity dISIP is obtained by taking a scaled linear combinationof two recorded magnetic fields. Thus, the ISIP data are secondary signalsthat are directly related to the chargeability of the earth. For any real,non-dispersive resistivity distribution the ISIP data, dISIP , should approxi-mately equal zero, while non-zero values indicate the presence of chargeablematerial.252.4. Synthetic examples2.4 Synthetic examplesThe sensitivity of the ISIP data to chargeable material will now be demon-strated using synthetic examples. The response from two different mod-els is tested. The first model consists of a changeable block buried in auniform half-space. The second will place the chargeable block in a non-uniform background made up of both resistive and conductive units. TheISIP data are simulated for a pair of electromagnetic surveys; a single largeloop transmitter with a roving receiver, and a Slingram survey with a con-stant transmitter-receiver offset.2.4.1 Example #1: Blocks in a uniform half-spaceThe first model features two blocks, buried in a uniform resistive back-ground. Block #1 is non-chargeable with a conductivity of 1 S/m. Block#2 has a zero-frequency conductivity (σ0) of 1 S/m, and is chargeable withCole-Cole parameters η = 0.1, τ = 0.1 and c = 0.5. Both blocks are 100m ona side, with the chargeable block centred at the point (250m, 500m) and thenon-chargeable block centred at the point (500m, 250m). The tops of theblocks are buried 125m below the surface of the half-space. The backgroundhalf-space has a conductivity of 10−3S/m.Large Loop EMIn this example, a large loop, electromagnetic survey is simulated. Thetransmitter is offset from the two conductive blocks, 200m on a side, andcentred at the point (200m, 200m). The transmitter location, and the loca-tion of the two blocks is shown in Figure 2.1. The three components of themagnetic field were simulated over the entire area.The magnetic fields are simulated at 1 and 2 Hz. This choice of frequen-cies and Cole-Cole parameters result in a change between the two frequenciesof 1.1× 10−2S/m in the real part and 7.5× 10−4S/m in the imaginary partof the resistivity of the chargeable block.The resulting magnetic fields and the calculated ISIP data are shown in262.4. Synthetic examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)TransmitterBlk #1Blk #2Figure 2.1: Plan view of the geometry of the two block model and the loca-tion of the transmitter wire for large loop examples. The tops of the blocksare 125m below the surface, and they extend to a depth of 225m. Block #1is conductive (1 S/m) but it is not chargeable. Block #2 is conductive andchargeable, with Cole-Cole parameters σ0 = 1S/m, η = 0.1, τ = 0.1 andc = 0.5. The dark black line shows the layout of the transmitter wire.272.4. Synthetic examplesFigure 2.2. The magnetic fields at the two frequencies are similar in formbut have different amplitudes. It is difficult to determine the presence ofthe blocks from the appearance of the magnetic fields, let alone determinewhether either of them are chargeable. Equation 2.13 is used to calculatethe ISIP data from these magnetic fields. The ISIP data are shown in thethird column in Figure 2.2. The existence and approximate location of thechargeable block is easily determined from looking at these plots.0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.9e-03-1.4e-03-9.6e-04-4.8e-040.0e+004.8e-049.6e-041.4e-031.9e-03(a) =( ~Bsx(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.6e-03-2.7e-03-1.8e-03-9.0e-040.0e+009.0e-041.8e-032.7e-033.6e-034.5e-03(b) =( ~Bsx(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-4.5e-06-3.6e-06-2.7e-06-1.8e-06-9.0e-070.0e+009.0e-071.8e-062.7e-063.6e-06(c) dISIPx0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.9e-03-1.4e-03-9.6e-04-4.8e-040.0e+004.8e-049.6e-041.4e-031.9e-03(d) =( ~Bsy(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.6e-03-2.7e-03-1.8e-03-9.0e-040.0e+009.0e-041.8e-032.7e-033.6e-034.5e-03(e) =( ~Bsy(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.1e-05-8.8e-06-7.0e-06-5.3e-06-3.5e-06-1.8e-060.0e+001.7e-06(f) dISIPy0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT4.0e-041.1e-031.8e-032.5e-033.2e-033.9e-034.6e-035.3e-036.0e-03(g) =( ~Bsz(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT8.0e-042.2e-033.6e-035.0e-036.4e-037.8e-039.2e-031.1e-021.2e-02(h) =( ~Bsz(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.2e-05-9.6e-06-7.2e-06-4.8e-06-2.4e-060.0e+002.4e-064.8e-067.2e-069.6e-06(i) dISIPzFigure 2.2: The x, y, and z components of the imaginary part of the magneticfields (nT) simulated using the large loop survey at 1 Hz and 2 Hz. Thecalculated ISIP data are shown in the third column. The centre of thechargeable block is located at (250m, 500m) and the centre of the conductiveblock is located at (500m, 250m).282.4. Synthetic examplesSlingram EMIn the next example, the response of a Slingram system was simulated.The transmitter was modelled to be a horizontal loop with a radius of 1m.The transmitter was moved throughout the area, and the fields from eachtransmitter location were recorded at a single receiver offset by 5m in the xdirection (figure 2.3).0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)Blk #1Blk #2Figure 2.3: Pan view of the geometry of the two block model and the locationof the transmitter-receiver pairs. Transmitters are shown with red dots,receivers with black dots.The simulated fields and the calculated ISIP data are shown in figure2.4. With this survey geometry it is easy to recognize the locations of thetwo blocks by looking at the fields. However, the ISIP data identifies onlythe location of the chargeable block.292.4. Synthetic examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.9e-07-1.7e-07-1.5e-07-1.3e-07-1.1e-07-8.7e-08-6.6e-08-4.5e-08-2.4e-08(a) =( ~Bsx(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.8e-07-3.4e-07-3.0e-07-2.6e-07-2.2e-07-1.7e-07-1.3e-07-9.0e-08-4.8e-08(b) =( ~Bsx(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.7e-09-3.0e-09-2.2e-09-1.5e-09-7.5e-100.0e+007.5e-101.5e-092.3e-093.0e-09(c) dISIPx0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-9.0e-08-7.2e-08-5.4e-08-3.6e-08-1.8e-080.0e+001.8e-083.6e-085.4e-087.2e-08(d) =( ~Bsy(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.7e-07-1.3e-07-8.4e-08-4.2e-080.0e+004.2e-088.4e-081.3e-071.7e-07(e) =( ~Bsy(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.0e-09-2.2e-09-1.5e-09-7.5e-100.0e+007.5e-101.5e-092.3e-093.0e-09(f) dISIPy0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT1.2e-061.3e-061.3e-061.4e-061.5e-061.5e-061.6e-061.6e-061.7e-061.8e-06(g) =( ~Bsz(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT2.4e-062.6e-062.7e-062.8e-062.9e-063.0e-063.2e-063.3e-063.4e-063.5e-06(h) =( ~Bsz(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT0.0e+002.4e-094.8e-097.2e-099.6e-091.2e-081.4e-081.7e-081.9e-082.2e-08(i) dISIPzFigure 2.4: The x, y, and z components of the imaginary part of the mag-netic fields (nT) simulated using the Slingram survey at 1 Hz and 2 Hz.The calculated ISIP data is shown in the third column. The centre of thechargeable block is located at (250m, 500m) and the centre of the conductiveblock is located at (500m, 250m).302.4. Synthetic examples2.4.2 Example #2: Blocks in a three-dimensionalbackgroundIn the second set of synthetic examples, the same two blocks were placed in acomplicated 3-D background and buried beneath a very resistive overburden.The conductivity of the background units ranged between 10−1S/m and10−4S/m. The overburden had a conductivity of 10−4S/m. Cross-sectionsthrough this conductivity model are shown in Figure 2.5.Large Loop EMThe resulting magnetic fields, and the calculated ISIP data are shown inFigure 2.6. The anomalous ISIP response from the chargeable block is clearlyevident, and is very similar to that of the previous example despite thefact that the chargeable body was buried in a host with a very differentconductivity structure.0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)σ010010−110−210−310−4(a) z = 150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)σ010010−110−210−310−4(b) y = 500mFigure 2.5: The real part of the complex three dimensional resistivity model.The background is not chargeable and the resistivity of the units varies from10−1S/m and 10−4S/m. The overburden is 50m thick and has a resistivity of10−4S/m. The two blocks have the same properties as they did in example#1. (a) Depth slice 150m below the surface. (b) Slice through the model aty = 500m.312.4. Synthetic examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-2.8e-03-2.1e-03-1.4e-03-7.0e-040.0e+007.0e-041.4e-032.1e-03(a) =( ~Bsx(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-5.6e-03-4.2e-03-2.8e-03-1.4e-030.0e+001.4e-032.8e-034.2e-03(b) =( ~Bsx(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-4.5e-06-3.6e-06-2.7e-06-1.8e-06-9.0e-070.0e+009.0e-071.8e-062.7e-063.6e-06(c) dISIPx0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.0e-03-2.3e-03-1.5e-03-7.5e-040.0e+007.5e-041.5e-032.2e-033.0e-033.8e-03(d) =( ~Bsy(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-6.0e-03-4.5e-03-3.0e-03-1.5e-030.0e+001.5e-033.0e-034.5e-036.0e-037.5e-03(e) =( ~Bsy(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.1e-05-8.8e-06-7.0e-06-5.3e-06-3.5e-06-1.8e-060.0e+001.7e-06(f) dISIPy0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT6.0e-041.4e-032.1e-032.9e-033.6e-034.4e-035.1e-035.9e-036.6e-037.4e-03(g) =( ~Bsz(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT1.3e-032.7e-034.3e-035.8e-037.3e-038.7e-031.0e-021.2e-021.3e-021.5e-02(h) =( ~Bsz(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.2e-05-9.6e-06-7.2e-06-4.8e-06-2.4e-060.0e+002.4e-064.8e-067.2e-06(i) dISIPzFigure 2.6: The x, y, and z components of the imaginary part of the mag-netic fields (nT) simulated using the large loop survey at 1 Hz and 2 Hz.The calculated ISIP data is shown in the third column. The centre of thechargeable block is located at (250m, 500m) and the centre of the conductiveblock is located at (500m, 250m).Slingram EMThe magnetic fields, and the resulting ISIP data are shown in Figure 2.7. Inthis example, the complicated background conductivity adds considerablestructure to the observed magnetic fields. The calculated ISIP data stillclearly shows the location of the chargeable block.322.4. Synthetic examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.2e-07-9.6e-08-7.2e-08-4.8e-08-2.4e-080.0e+002.4e-084.8e-087.2e-08(a) =( ~Bsx(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-2.4e-07-1.9e-07-1.4e-07-9.6e-08-4.8e-080.0e+004.8e-089.6e-081.4e-07(b) =( ~Bsx(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.7e-09-3.0e-09-2.2e-09-1.5e-09-7.5e-100.0e+007.5e-101.5e-092.3e-093.0e-09(c) dISIPx0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-7.2e-08-5.4e-08-3.6e-08-1.8e-080.0e+001.8e-083.6e-085.4e-087.2e-08(d) =( ~Bsy(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.4e-07-1.1e-07-7.2e-08-3.6e-080.0e+003.6e-087.2e-081.1e-071.4e-07(e) =( ~Bsy(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.0e-09-2.2e-09-1.5e-09-7.5e-100.0e+007.5e-101.5e-092.3e-093.0e-09(f) dISIPy0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT3.0e-074.1e-075.1e-076.2e-077.2e-078.3e-079.3e-071.0e-061.1e-06(g) =( ~Bsz(ω1))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT6.0e-078.1e-071.0e-061.2e-061.4e-061.7e-061.9e-062.1e-062.3e-06(h) =( ~Bsz(ω2))0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT0.0e+002.4e-094.8e-097.2e-099.6e-091.2e-081.4e-081.7e-081.9e-082.2e-08(i) dISIPzFigure 2.7: The x, y, and z components of the imaginary part of the mag-netic fields (nT) simulated using the Slingram survey at 1 Hz and 2 Hz.The calculated ISIP data is shown in the third column. The centre of thechargeable block is located at (250m, 500m) and the centre of the conductiveblock is located at (500m, 250m).332.5. Sources of error2.5 Sources of errorThe results in the last section show that ISIP data can be a direct indicatorof the presence of chargeable material. The important question remains as towhether or not the ISIP signal can actually be measured in field situations.There are two factors that must be considered. The first is the intrinsicnoise arising from the instrument and the measurement procedure. In orderto be useful, the ISIP signal must be larger than this noise. The secondfactor is the noise inherent in the ISIP derivation when the low frequencyassumption is violated.2.5.1 Noisy magnetic fields and error propagationSince an ISIP datum is a linear combination of two measured data, itsaccuracy is determined by the accuracy of each measurement used in itscalculation. Assuming that the observation of =(~Bs(ω1))is contaminatedwith noise with standard deviation ω1 , =(~Bs(ω2))is contaminated withnoise with standard deviation ω2 , and that the errors in the two observationsare uncorrelated, then the standard deviation in the resulting ISIP data isgiven byISIP =√2ω2 +(ω2ω1ω1)2(2.14)For the case where ω1 = ω2 = ω this simplifies toISIP =√1 +(ω2ω1)2ω (2.15)Thus, the standard deviation of the uncertainty in the ISIP data willalways be larger than the standard deviations in the measurements of themagnetic fields used to calculate the ISIP response. It is important to notethat while the magnitude of the standard deviations will only slightly in-crease for closely spaced frequencies, the magnitude of the ISIP data willbe significantly less than the magnitude of the original magnetic fields as342.5. Sources of errora result of equation 2.13. This results in an increase in the relative uncer-tainty of the ISIP data compared to the relative uncertainty in the originalmagnetic fields.The magnitude of the ISIP signal is controlled by many factors. Prin-cipally, it depends upon the size and geometry of the target, the targetslocation relative to the transmitter, the geometry of the transmitter, themagnitude of the current, the nature of the frequency dependence of the ma-terial, and the choice of frequencies. The relationship between most of theseparameters and the ISIP data is complicated. The exception is transmittercurrent; the ISIP data depend linearly on the magnitude of the current.Using modern SQUID magnetometers, resolutions of 20 fT/√Hz havebeen achieved at 1 Hz (Kawai et al., 1999). Thus, averaging the recordedsignal for 100 s with this type of sensor would result in an effective sensornoise of 2 fT. If the ISIP data were calculated from observations of themagnetic field at 1 and 2 Hz, then ISIP = 4.47 fT.To illustrate the fact that ISIP data can be large enough to be recordedwith this type of instrumentation, the large loop survey from example #1will be considered again. For a 1A transmitter current, the highest mag-nitude ISIP signal in Figure 2.2 occurs in the z-component and is approx-imately 12 fT. If the fields used to compute the ISIP data were contami-nated with random Gaussian noise with a standard deviation of 2 fT thanthe standard deviation of the uncertainty in the resulting ISIP datum wouldbe 4.47 fT. Figure 2.8 presents the noisy ISIP data that would be measuredwhen the transmitter current is 1 and 10 A. EM transmitters capable ofproducing currents up to (and above) 10A are commonly available on themarket today. With a 10A transmitter, and resolving power of 2fT, the ISIPresponse is easily visible in the data (Figure 2.8b).This analysis of noise is fairly simplistic and considers only the presencesensor noise in the observed data. The time varying geomagnetic fields ofthe Earth (Chwala et al., 2013) would introduce an additional source ofnoise potentially more significant than sensor noise. Accurate extractionof the imaginary component of the observed signal would also require anextremely precise understanding of the system geometry and transmitted352.5. Sources of error0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-1.4e+01-1.1e+01-7.2e+00-3.6e+000.0e+003.6e+007.2e+001.1e+011.4e+01(a) 1A0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-1.2e+02-9.6e+01-7.2e+01-4.8e+01-2.4e+010.0e+002.4e+014.8e+017.2e+019.6e+01(b) 10AFigure 2.8: The vertical component of the ISIP data (fT) calculated frommagnetic fields that had been contaminated with random Gaussian noisewith a standard deviation of 2 fT and a variable transmitter current. (a) 1Amp transmitter current (b) 10 Amp transmitter currentwaveform (Ley-Cooper et al., 2006). That being said, this example doesnot represent the best case scenario. Different transmitter and receiver ge-ometries, a shallower target, or a different Cole-Cole model could achieve ahigher magnitude response.2.5.2 Low frequency assumptionThe second area in which errors can contaminate the analysis occurs if thelow frequency assumptions are violated. Low frequencies are required in theISIP derivation for two reasons. First, while deriving the expressions for theISIP data, it was assumed that observations of the magnetic fields are madeat sufficiently low frequencies that the higher order terms of the expansion ofequation 2.10 can be dropped. The imaginary portion of the recorded mag-netic field then depends linearly upon frequency. The value of the inductionnumber is the determining factor in the validity of this assumption. As longas the induction number is much less than unity, the ISIP data will be equalto zero unless chargeable material is present. However, as the value of theinduction number increases, additional sources of signal become apparentin the data. Insight can be gained into the validity of this assumption byconsidering a simple analytic example.362.5. Sources of errorThe vertical magnetic fields measured by a horizontal co-planar loop atthe surface of a conductive half-space (Ward and Hohmann, 1988) is givenbyHz =m2pik2r5[9− (9 + 9ikr − 4k2r2 − ik3r3) e−ikr] (2.16)where m is the dipole moment of the transmitter loop, k is the wave number(k = (−iµσω)1/2), and r is the distance between the transmitter and thereceiver. By defining the quantity β asβ = ikr = −(1 + i)α (2.17)where α is the induction numberα =rδ= r(σµω2)1/2(2.18)Equation 2.16 can be rewritten asHz =m2piβ2r3[(9 + 9β + 4β2 + β3)e−β − 9](2.19)The behaviour of the function at low wavenumber can be approximatedby expanding e−β about β = 0 and substituting the result into equation2.19. This givesHz ≈ − m4pir3(1 +β24− 4β315+β48− 4β5105+ · · ·)(2.20)If a non-dispersive conductivity is assumed, this expression is easily sep-arated into real and imaginary parts with<[Hz] ≈ − m4pir3(1− 8α315− α42− 16α5105+ · · ·)(2.21a)=[Hz] ≈ − m4pir3(α22+8α315− 16α5105+ · · ·)(2.21b)For low induction numbers (α << 1 or r << δ) the higher order terms drop372.5. Sources of errorleaving=[Hz] ≈− mα28pir3≈− mσµω16pir(2.22)As expected, this expression is dominated by a term that varies as ω forlow induction numbers.This assumption will not hold for values of α that are not much less than1, thus, the ISIP definition will break down at high frequency, in conductiveenvironments, and for large transmitter receiver offsets. An example of thelinear approximation (equation 2.22) compared to the true response (equa-tion 2.16) is shown in figure 2.9 for a Slingram survey with a transmitter-receiver offset of 5m, above a 1S/m half-space.101102103104Frequency (Hz)10-810-710-610-510-410-3ℑ(Hz) (A/m)Analytic ResponseLinear Approximation10-210-1100Induction NumberFigure 2.9: The linear approximation of =(Hz) compared to the true analyticresponse 5m from a vertical magnetic dipole on the surface of a 1S/m half-space. The linear approximation matches the true response while α << 1.The validity of the linearity assumptions can also be examined usingthe synthetic examples shown in section 2.4. The model and survey ge-ometries from example #2 are used again, but the data will be computedusing frequencies from 5 − 55 Hz. The resultant z-component of the ISIP382.6. Conclusionsdata are shown in figures 2.10 and 2.11. ISIP data require measurementsof the magnetic field at two frequencies. The lower of the two frequenciesused are shown in the figure labels. The higher frequency is 5% larger.The chargeable block has been modified to have the ColeCole parametersη = 0.1, τ = 0.1 and c = 0.25. This was done so that the value of δρ remainsrelatively constant for each pair of frequencies.At low frequencies, the ISIP data show a clear response centred overthe chargeable block, with zero response elsewhere in the area. As thefrequencies increase, non-zero signals begin to appear in the ISIP data awayfrom the chargeable block. At these frequencies, starting at around 25 Hz,the magnetic fields no longer vary linearly with frequency. This result agreeswith the theoretical example shown in Figure 2.9, were deviation from thelinear approximation begins to be noticeable in the vicinity of this frequency.2.6 ConclusionsIn this work, a new methodology to identify the presence of chargeable mate-rial using frequency domain electromagnetics was proposed. The techniqueexploits the simple asymptotic behaviour of the fields at low frequencies and,if observed at low enough frequency, these data are identically zero if theconductivity is purely real and frequency independent. Thus any non-zerovalue of this datum is a direct indicator of chargeable material. Numeri-cal simulations demonstrate that this is true even in a complex geologicalenvironment.This technique provides a new methodology for detecting and mappingthe presence of chargeable material without needing to inject current intothe ground or place electrodes to measure potentials. By avoiding these re-quirements this technique may prove to be a useful tool in geological settingswhere traditional IP surveys are difficult to conduct.392.6. Conclusions0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-2.1e-06-1.0e-060.0e+001.1e-062.1e-063.1e-064.2e-065.2e-066.3e-06(a) 5Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT3.1e-053.4e-053.7e-054.0e-054.2e-054.5e-054.8e-055.1e-055.4e-055.6e-05(b) 15Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT4.6e-055.0e-055.5e-055.9e-056.3e-056.7e-057.1e-057.6e-058.0e-05(c) 25Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT1.3e-041.4e-041.4e-041.4e-041.4e-041.5e-041.5e-041.5e-041.6e-041.6e-04(d) 35Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT2.3e-042.3e-042.4e-042.4e-042.5e-042.5e-042.6e-042.6e-042.7e-04(e) 45Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT3.2e-043.3e-043.4e-043.5e-043.6e-043.8e-043.9e-044.0e-044.1e-044.2e-04(f) 55HzFigure 2.10: Vertical component of ISIP data for the grounded loop survey.The ISIP data is calculated for pairs of frequencies with increasing magni-tudes. The lower frequency is shown in the figure. The second frequency is5% higher.402.6. Conclusions0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT0.0e+001.1e-092.1e-093.2e-094.2e-095.2e-096.3e-097.4e-098.4e-09(a) 5Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.0e-090.0e+003.0e-096.0e-099.0e-091.2e-081.5e-081.8e-082.1e-082.4e-08(b) 15Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.4e-08-9.6e-09-4.8e-090.0e+004.8e-099.6e-091.4e-081.9e-082.4e-08(c) 25Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-3.8e-08-3.4e-08-2.9e-08-2.4e-08-1.9e-08-1.4e-08-9.6e-09-4.8e-090.0e+004.8e-09(d) 35Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-6.7e-08-6.3e-08-5.9e-08-5.6e-08-5.2e-08-4.9e-08-4.5e-08-4.1e-08-3.8e-08-3.4e-08(e) 45Hz0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-1.2e-07-1.1e-07-1.0e-07-9.3e-08-8.6e-08-7.8e-08-7.1e-08-6.3e-08-5.6e-08-4.8e-08(f) 55HzFigure 2.11: Vertical component of ISIP data for the Slingram survey. TheISIP data is calculated for pairs of frequencies with increasing magnitudes.The lower frequency is shown in the figure. The second frequency is 5%higher.41Chapter 3Inversion of Inductive SourceInduced Polarization Data3.1 IntroductionInverse modelling is a common and valuable tool for interpreting geophysicaldata. Solving the inverse problem involves finding a physical property modelthat can predict the observed data, subject to minimizing a user definedobjective function. The inversion process converts information containedin the observations from data space to a one, two, or three dimensionalmodel space. For many types of geophysical surveys, it is very difficult todraw physical interpretations directly from the data. Therefore, invertingthe data can considerably improve data interpretation and expose subtlefeatures that might otherwise have been missed.In this chapter an inversion scheme to recover the three dimensionalchargeability distribution from ISIP data is developed. The inversion method-ology employed is closely based on existing inversion schemes for invert-ing conventional induced polarization data developed by Oldenburg and Li(1994) and Li and Oldenburg (2000). Section 3.2 will show that ISIP dataare approximately linear with respect to the change in conductivity over thefrequencies being used. Section 3.3 presents the discretization of this result.Section 3.4 reviews the theory of linear inverse problems, and finally, Section3.5 presents a series of synthetic examples of the inversion of ISIP data.423.2. Linearization of the ISIP sensitivities3.2 Linearization of the ISIP sensitivitiesInverse modelling first requires an understanding of the connection betweenchanges in the Earth’s chargeability and changes in the ISIP data. Thepresence of chargeability will result in small perturbations in the conductiv-ity with respect to frequency at which they are measured. Let σ1 and σ2be the complex conductivities observed at frequencies ω1 and ω2. The twofrequencies are closely spaced such that σ2 can be written in terms of σ1plus a small perturbation,σ2 = σ1 + δσ (3.1)Now, consider a magnetic field ~Bs (ω, σ) at frequency ω existing in aconductivity distribution σ. An approximation for the magnetic field atω = ω2 and σ = σ2 can be found by expanding the field about σ = σ1,discarding higher order terms, and evaluating the expansion at ω2, to obtain~Bs (ω2, σ2) ≈ ~Bs (ω2, σ1) + ∂~Bs∂σ∣∣∣∣(ω2,σ1)δσ (3.2)The first term on the right hand side of this expression is the secondarymagnetic field that would be observed in the presence of the complex, fre-quency independent conductivity distribution σ1. As long as the frequenciesused are sufficiently low, then, from the results in Chapter 2=[~Bs (ω2, σ1)]≈ ω2ω1=[~Bs (ω1, σ1)](3.3)By considering only the imaginary part of Equation 3.2=[~Bs (ω2, σ2)]≈ =[~Bs (ω2, σ1)]+ =[∂ ~Bs∂σ∣∣∣∣(ω2,σ1)δσ](3.4a)≈ ω2ω1=[~Bs (ω1, σ1)]+ =[∂ ~Bs∂σ∣∣∣∣(ω2,σ1)δσ](3.4b)Finally, rearranging to make use of the definition of the ISIP datum433.3. Discretization of the ISIP dataleavesdISIP ≈ =[∂ ~Bs∂σ∣∣∣∣(ω2,σ1)δσ](3.5)Expanding this expression in terms of the real and imaginary parts givesdISIP ≈ <[∂ ~Bs∂σ∣∣∣∣(ω2,σ1)]= [δσ] + =[∂ ~Bs∂σ∣∣∣∣(ω2,σ1)]< [δσ] (3.6)This expression approximates the ISIP data using the sensitivities of themagnetic fields to conductivity perturbations and changes in conductivitywith frequency. Using this expression, a discrete linear inverse problem canbe formulated to recover the distribution of chargeable material.3.3 Discretization of the ISIP dataIn order to invert dISIP for the distribution of chargeable material, theproblem must first be discretized. The volume of interest is discretized ontoan orthogonal tensor mesh (Figure 3.1). The physical properties are placedat cell centres, electric fields are located on cell edges, and magnetic fieldsare located on cell faces.The discrete forms of Faraday’s law and Ampere’s law (Equation 2.4)are given byCE + iωBs = −iωB0 (3.7a)C>Mfµ−1Bs −MeσE = 0 (3.7b)where C is the discrete curl operator, which maps from cell edges to cellfaces. Mfµ−1 and Meσ are mass matrices containing the inverse magneticpermeability (averaged to cell faces) and the conductivity (averaged to celledges), respectively, and E and B are the discrete electric and magneticfields. Further details on the formation of discrete differential operators andthe mass matrices can be found in Appendix A.443.3. Discretization of the ISIP dataFigure 3.1: Example orthogonal tensor mesh of the type used in this chapter.The linearized forward model in Equation 3.6 requires the derivative of~Bs with respect to conductivity, evaluated at ω2 and σ1. Differentiating 3.7with respect to σ givesC∂E∂σ+ iω∂Bs∂σ= 0 (3.8a)C>Mfµ−1∂Bs∂σ−Meσ∂E∂σ− diag(E)Ace>diag(v) = 0 (3.8b)where Ace is an averaging operator, acting from cell edges to cell centres.Eliminating ∂E∂σ gives(Mfµ−1CMeσ−1C>Mfµ−1 + iωMfµ−1) ∂Bs∂σ= Mfµ−1CMeσ−1diag(E)AceTdiag(v)(3.9)This system can be solved for ∂B∂σ at ω2 and σ1 given E(ω2, σ1).Unfortunately, the system has a non-trivial null space (the gradient ofany scalar function) when ω is small and therefore must be stabilized. The453.3. Discretization of the ISIP datacondition∇ · ~Bs = 0 (3.10)is imposed to stabilize the system by subtracting Mfµ−1D>Meσ−1DMfµ−1∂B∂σfrom the left side of Equation 3.9, obtaining(Mfµ−1CMeσ−1C>Mfµ−1 −Mfµ−1D>Meσ−1DMfµ−1 + iωMfµ−1)∂Bs∂σ(3.11)= Mfµ−1CMeσ−1diag(E)Ace>diag(v)This guaranties a solution of the system even at very low frequency (Haberand Ascher, 2001).The sensitivity matrix J is defined byJ = Q∂Bs∂σ∣∣∣∣(ω2,σ1)(3.12)where Q is a projection matrix that projects the magnetic fields from cellfaces to receiver locations. The (i, j) element of J is a measure of how theith observation of ~Bs will be affected by a small perturbation in conductivityin the jth cell.The discrete form of the forward problem given by Equation 3.6 is nowdISIP = <[J]=[δσ] + =[J]<[δσ] (3.13)Note that dISIP is dependent on both <[δσ] and =[δσ]. While it would bepossible to use Equation 3.13 to invert for these two quantities, some simpleassumptions are made to linearize the problem.Beginning by examining the right hand side of Equation 3.9Mfµ−1CMeσ1−1diag(E)Ace>diag(v) (3.14)the first assumption assumes a low frequency such that the total magneticfield is dominated by the contribution of the primary field from the trans-mitter. Faraday’s law then implies that the resulting electric fields are dom-463.3. Discretization of the ISIP datainated by their imaginary part, or E ≈ =(E).The second assumption assumes relatively low chargeability values. Ifthis is the case, the elements of Meσ1 are dominated by its real part, and areapproximately equal to Meσ0 . Thus, the right hand side of Equation 3.9 isprimarily imaginary.Examining the matrix on the left hand side of Equation 3.9Mfµ−1CMeσ1−1C>Mfµ−1 −Mfµ−1D>Meσ1−1DMfµ−1 + iωMfµ−1 (3.15)the discretization of the differential terms are of order h−2µ−2σ−1 while thelast term is of order ωµ−1, where h is the length of one edge of the smallestcells in the mesh. If ω h−2σ−1µ−1 then the term involving ω can beneglected, and thus for low frequencies this matrix is primarily real.To this end, it can be concluded thatJ ≈ =[J] (3.16)orJ ≈ QA−1G (3.17)whereA = Mfµ−1CMeσ0−1C>Mfµ−1 −Mfµ−1D>Meσ0−1DMfµ−1 (3.18)andG = Mfµ−1CMeσ0−1diag(=(E))Ace>diag(v) (3.19)Since <[J] is very small, the term <[J]=[δσ] is dropped from Equation3.13 leaving a linear relationship between the ISIP data and the change inthe real conductivities.dISIP ≈ J<[δσ] (3.20)473.3. Discretization of the ISIP data0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-5.4e-18-4.5e-18-3.6e-18-2.7e-18-1.8e-18-9.0e-190.0e+009.0e-191.8e-182.7e-18(a) =[ ~Bs(ω2)]− ω2ω1=[ ~Bs(ω1)]0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-5.4e-18-4.5e-18-3.6e-18-2.7e-18-1.8e-18-9.0e-190.0e+009.0e-191.8e-182.7e-18(b) =[δσ]<[J] + <[δσ]=[J]0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)nT-5.4e-18-4.5e-18-3.6e-18-2.7e-18-1.8e-18-9.0e-190.0e+009.0e-191.8e-182.7e-18(c) <[δσ]=[J]Figure 3.2: The true ISIP data (a) compared to the approximations madein equation 3.13 (b) and equation 3.20 (c). The data is calculated for themodel and the survey layout shown in Figure 2.1.483.4. Solving the linear inverse problemThe true ISIP data, as well as the approximations from Equations 3.13and 3.20 are pictured in Figure 3.2 for the example previously presented inFigure 2.2. Both of the approximations provide an excellent result.3.4 Solving the linear inverse problemWhen solving an inverse problem, the goal is to recover a distribution ofphysical properties in the earth (the model, m) that can explain a set ofgeophysical observations (the data, dobs) to a desired level of accuracy. Todo this, the earth is discretized into cells (as in Figure 3.1), where each cellis assumed to have a constant physical property. Many cells are required toadequately represent the complicated geometry of the earth, and in threedimensional problems there are far more cells than there are data. This typeof problem, known as under-determined, has an infinite number of uniquemodels which reproduce the observations to the desired level of accuracy.This section will provide a brief overview of the topic of linear inverseproblems. Oldenburg and Li (2005) provide an excellent survey of the sub-ject in addition to what is included here.3.4.1 Data misfitThe data misfit between a given model and observed data is measured usingthe sum of the squares,φd =N∑i=1(F [m]i − dobsii)2= ||Wd(F [m]− dobs)||22 (3.21)In this expression, F [m] is the data predicted by the model m, i is theuncertainty in the ith datum, and Wd is a diagonal weighting matrix, wherethe diagonal elements contain the inverse of the data uncertainties. This def-inition of the misfit is convenient as it is equal to the number of observationswhen the observed and predicted data differ by .493.4. Solving the linear inverse problem3.4.2 RegularizationThe non-uniqueness of the problem requires a regularization such that aunique result can be obtained. This is achieved through the design of amodel objective function whose value is minimized for models exhibitingdesirable characteristics.In this work, the regularization is designed to favour smoothly varyingmodels that are close to a given reference model. This is a common choice inmany geophysical inverse problems. The model objective function is definedbyφm =αs||Ws (m−mref ) ||22 + αx||Wx (m−mref ) ||22+ αy||Wy (m−mref ) ||22 + αz||Wz (m−mref ) ||22=||Wm(m−mref )||22 (3.22)where Ws is a diagonal matrix, and Wx, Wy and Wz are discrete ap-proximations of the first derivative operator in the x, y and z directionsrespectively. The α’s are weighting parameters that balance the relativeimportance of producing small or smooth models.3.4.3 Solving the optimization problemThe inverse problem of estimating m is accomplished by solving the followingconstrained optimization problemminφ = φd(m) + βφm(m) (3.23)s.t. mlower ≤m ≤mupperwhere mlower and mupper are upper and lower bound constraints, set foreach cell in the model.The regularization or trade-off parameter β defines the relative impor-tance of minimizing the data misfit or minimizing the model objective func-tion. The regularization parameter is chosen using a cooling schedule whereβ is initially set at a very large value so that βφm dominates the objective503.4. Solving the linear inverse problemfunction. When the m that minimizes Equation 3.23 is identified, β is de-creased by a constant factor and the optimization problem is solved again.This continues until the desired level of data misfit is achieved.For each particular value of β, the constrained optimization problem issolved using a projected Gauss-Newton algorithm (Kelley, 1999). The modelis updated iteratively until a stopping criteria is met.The search direction, δm, is generated by solvingR˜δm = −g (3.24)where g is the gradient of the objective function with respect to m, and isgiven byg = J>W>d Wd(F [m]− dobs) + βW>mWm(m−mref ) (3.25)The matrix J is the sensitivity, ∂F [m]∂m , and R˜ is the reduced Hessian definedbyR˜ =H˜ij if i ∈ A(m) or j ∈ A(m)δij otherwise (3.26)H˜ = J>W>d WdJ + βW>mWm (3.27)The matrix A(m) is the set of active constraints, that is, the cells whosevalues are between, but not equal, to the upper and lower bounds.Upon computing δm, the updated model is given bymn+1 = Proj(mn + γδm) (3.28)where γ (0 < γ ≤ 1) is chosen by a simple backtracking line search suchthat mn+1 reduces the objective function. Proj is a projection operator thatprojects the updated model back to the bounded model space and thereforereassigns parameters to active and inactive sets.513.4. Solving the linear inverse problem3.4.4 ISIP specificsThe linear approximation for the ISIP data (Equation 3.20) is used to formu-late an inverse problem which is then solved to recover an estimate of <[δσ].This recovered quantity is not the chargeability in the sense of the Cole-Colemodel, however, it is a direct proxy for the distribution of chargeable ma-terials which exhibit a frequency dependence at the operating frequencies.The presence of a zero <[δσ] does not necessarily imply that the material isnot chargeable. It does however indicate that the dispersion of the materialis such that there is no measurable <[δσ] between the frequencies used.Prior to inverting the ISIP data, an approximation for the matrix =[J] isrequired. This matrix depends on the background conductivity distributionσ0, as well as the electric fields that would exist in the presence of this con-ductivity distribution. This requirement is very similar to the conventionalinduced polarization case, where an estimate of σ0 is required in order toinvert for chargeability. While σ0 is unknown, it is commonly possible to ob-tain a reasonable estimate. The estimated conductivity could be obtainedby performing a 3D inversion of data at one of the two frequencies, or itcould be generated using some other source of information. The importanceof obtaining a good estimate of the background conductivity distribution tothe success of the inversion will be considered in Section 3.5.Forward modelling the predicted data involves multiplying the sensitivitymatrix (Equation 3.17) by the model vector. This matrix is large (numberof observations times number of cells) and dense. For problems of evenmoderate size this matrix cannot be stored in memory. Rather than formthe matrix completely, the matrix vector product is carried out in threesteps. First, the model vector, m, is multiplied by the matrix G (Equation3.19),y = Gm (3.29)523.5. Synthetic ExamplesNext, the systemAx = y (3.30)is solved for x. Finally, the result is projected to the data locationsdISIP = Qx (3.31)The system 3.30 is solved using the MUMPs direct solver (Amestoy et al.,2001). The sensitivity matrix does not depend on m, and thus does notchange throughout the inversion. This means that the matrix A only needsto be factored once, and the factors can be stored. Once factored, A−1 iseasily and efficiently multiplied onto a vector with a forward and a backwardsubstitution.The projected Gauss-Newton equation (Equation 3.24) is solved itera-tively using the conjugate gradient algorithm. This requires only the abilityto multiply J and J> onto a vector. Multiplying J onto a vector is a forwardmodelling, and is accomplished in the same way as is described above. Sincethe matrix A is symmetric,J> = G>A−1Q> (3.32)The product of J> and a vector is again performed in three steps; first multi-plying by Q>, solving Ax = y, and finally multiplying by finally multiplyingby G>.3.5 Synthetic ExamplesThe previously described inversion routine will now be demonstrated on aseries of synthetic examples. All of examples include the same model andsurvey geometry.The survey geometry (Figure 3.3) consists of nine transmitters, laid outin a three by three grid. Each of the transmitters are 200m on a side. Thetransmitters frequencies are 1Hz and 2Hz, with a peak to peak transmitter533.5. Synthetic Examplescurrent of 30A. All three components of the magnetic field are recorded at169 receivers laid out in a 13 by 13 grid.0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)Figure 3.3: The layout of the survey used in the synthetic inversion exam-ples. The locations of the nine transmitters are shown with dashed lines,and the locations of the receivers shown with black dots. The extent of thechargeable block is shown in light grey, and the extent of the conductiveblock in dark grey.The conductivity model used in the examples consists of a pair of blockslocated in a variable background. Both of the blocks have a zero-frequencyconductivity of 1S/m, and the background units have conductivities rang-ing from 10−4S/m to 10−1S/m. Both of the blocks are 125m by 125m by100m and are buried 125m below the surface. Slices from the zero-frequencyconductivity model are shown in Figure 3.4.The block in the North-West corner of the model is chargeable withCole-Cole parameters of η = 0.1, τ = 0.1s, and c = 0.5. These parametersresult in a <[δσ] of 0.011S/m between 2Hz and 1Hz. Slices of the true <[δσ]model are shown in figure 3.5.543.5. Synthetic Examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)ℜ(σ1)10010−110−210−310−4(a) z = -150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)ℜ(σ1)10010−110−210−310−4(b) y = 500mFigure 3.4: True σ0 model used in synthetic inversion examples.0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)ℜ(δσ)0.000.01(a) z = -150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)ℜ(δσ)0.000.01(b) y = 500mFigure 3.5: True <[δσ] model used in synthetic inversion examples.553.5. Synthetic ExamplesThe z-component of the simulated ISIP data resulting from this surveygeometry and model is pictured in Figure 3.6. With this layout there are4563 unique data points. Prior to inversion, the data were contaminatedwith random Gaussian noise with a standard deviation of 45fT. This isapproximately the noise level one would expect from a receiver with a 20fTresolution.This data set was inverted three times, using a three different conduc-tivity models for the sensitivity calculations. The first example inverted thedata using the real part of the true conductivity model. The second andthird examples tested the impact of the conductivity model on the inversionresult. The second example used a simple half-space conductivity model,and the final inversion used the conductivity model recovered by invertingthe magnetic fields at one frequency while assuming a non-chargeable model.In each of the inversions a zero <[δσ] model was used for both the startingand reference models. The uncertainty of the data was set at 5% of themagnitude of the data values plus a floor of 45fT.3.5.1 Example #1: True sensitivitiesIn the first example, the real part of the true background conductivity modelwas used for sensitivity calculations. This represents an ideal, but not nec-essarily realistic case.A plan view and cross-section of the true and recovered <[δσ] modelsare shown in Figure 3.7. The depth and horizontal extents of the chargeablematerial are well located. As usual, with such inversions, there is someextension of the chargeability away from the boundaries of the true blockand the amplitude of the recovered chargeability is lower than the truevalue. However, this inversion has been successful in locating the chargeablematerial.563.5. Synthetic Examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-3.0e+020.0e+003.0e+026.0e+029.0e+021.2e+031.5e+031.8e+032.1e+032.4e+03(a)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-4.2e+020.0e+004.2e+028.4e+021.3e+031.7e+032.1e+032.5e+032.9e+03(b)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-4.5e+02-3.8e+02-3.0e+02-2.2e+02-1.5e+02-7.5e+010.0e+007.5e+011.5e+022.2e+02(c)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-3.6e+02-2.4e+02-1.2e+020.0e+001.2e+022.4e+023.6e+024.8e+026.0e+027.2e+02(d)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-3.0e+02-1.5e+020.0e+001.5e+023.0e+024.5e+026.0e+027.5e+029.0e+02(e)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-3.6e+02-3.0e+02-2.4e+02-1.8e+02-1.2e+02-6.0e+010.0e+006.0e+011.2e+02(f)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-3.9e+02-3.4e+02-2.8e+02-2.2e+02-1.7e+02-1.1e+02-5.6e+010.0e+005.6e+01(g)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-3.8e+02-3.4e+02-2.9e+02-2.4e+02-1.9e+02-1.4e+02-9.6e+01-4.8e+010.0e+004.8e+01(h)0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)fT-2.9e+02-2.5e+02-2.2e+02-1.8e+02-1.4e+02-1.1e+02-7.2e+01-3.6e+010.0e+003.6e+01(i)Figure 3.6: The vertical component of the synthetic ISIP data used forinversions. The dashed black squares indicate the location of the transmitterwire.573.5. Synthetic Examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)ℜ(δσ)0.0000.0010.0020.0030.0040.0050.0060.007(a) z = -150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)ℜ(δσ)0.0000.0010.0020.0020.0030.0040.0050.0060.0060.007(b) y = 500mFigure 3.7: Recovered <[δσ] model for example #1. The real part of the trueconductivity model was used to generate sensitivities for the <[δσ] inversion.3.5.2 Example #2: Sensitivity approximated withhalf-spaceIn reality, the true distribution of the conductivity in the earth is not knownand thus an estimate of the conductivity structure is required. Depending onthe method used to generate the estimate, the resulting conductivity modelmay not be a good representation of the true conductivities in the area ofinterest. In order to test the importance of the background conductivitymodel, the next example will use a 5× 10−3S/m half-space to generate thesensitivities.The resulting chargeability model obtained using sensitivities from a5×10−3S/m half-space is shown in figure 3.8. A chargeable body was clearlyrecovered, but the resolution is substantially reduced compared to that infigure 3.7. The body has moved toward the surface and is also spread out inthe horizontal direction. Nevertheless, the result provides useful informationas its maximum value coincides horizontally with the centre of the trueprism. This is a positive result and shows that knowledge of the background583.5. Synthetic Examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)ℜ(δσ)2.0e-044.0e-046.0e-048.0e-041.0e-031.2e-031.4e-031.6e-031.8e-03(a) z = -150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)ℜ(δσ)0.0e+005.0e-041.0e-031.5e-032.0e-03(b) y = 500mFigure 3.8: Recovered <[δσ] model for example #2. A 5 × 10−3S/m half-space was used to generate sensitivities for the <[δσ] inversion.is important but not critical to getting some valuable information from theISIP data. Moreover, this is a fairly extreme example. The true conductivityvaries between 1 and 10−1 S/m and it has been replaced by a uniform earthof 5× 10−3S/m.3.5.3 Example #3: Sensitivity approximated withrecovered conductivitiesIn the final example, noisy magnetic fields at 1Hz were first inverted usingthe FEM inversion code EH3D (Haber et al., 2004) to recover a real 3Dconductivity model. Plan view and cross-sections of the true and recoveredconductivity models are shown in Figure 3.9. As only a single frequency wasinverted, the resulting resistivity model differs substantially from the truemodel. The major features of the blocks are recognizable in the results, butthey are highly distorted.The inversion of the ISIP data using the recovered conductivities areshown in 3.10. The chargeable block is now well recovered and is nearly the593.5. Synthetic Examples0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)ℜ(σ1)10010−110−210−310−4(a) True conductivity at z = -150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)ℜ(σ1)10010−110−210−310−4(b) True conductivity at y = 500m0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)σ110010−110−210−310−4(c) Recovered conductivity at z = -150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)σ110010−110−210−310−4(d) Recovered conductivity at y = 500mFigure 3.9: True and recovered conductivity models. The magnetic fields at1Hz were inverted to generate the recovered model.603.6. Conclusionsame as that obtained when using the true conductivity model to calculatethe sensitivities.0 100 200 300 400 500 600 700Easting (m)0100200300400500600700Northing (m)ℜ(δσ)0.0000.0020.0030.0050.0060.0070.0090.011(a) z = -150m0 100 200 300 400 500 600 700Easting (m)−500−400−300−200−1000Depth (m)ℜ(δσ)0.0000.0020.0030.0050.0060.0070.009(b) y = 500mFigure 3.10: Recovered <[δσ] model for example #3. The recovered con-ductivity model (Figure 3.9) was used to generate sensitivities for the <[δσ]inversion.3.6 ConclusionIn this chapter an inversion scheme was developed for the ISIP data de-fined in Chapter 2. The scheme involved the development of a linearizedapproximation relating ISIP data to chargeability. The resulting sensitivitymatrix, which relates ISIP data to chargeability, was evaluated using a real,frequency independent conductivity model.Through synthetic inversion examples, I show that important informa-tion about the existence and location of the chargeable structure can beobtained. This is true even with a fairly poor knowledge about the back-ground conductivity is available.61Part IIThree DimensionalModelling of IP Effects inTime DomainElectromagnetic Data62Chapter 4Convolutionary ForwardModelling of Time DomainElectromagnetics4.1 IntroductionThe economic importance of chargeable materials has generated a large bodyof work focusing on the forward modelling of their electromagnetic response.The majority of the early work focused on modelling the response of geomet-rically simple polarizable bodies. This is accomplished by either transform-ing results from the frequency-domain to the time-domain (Bhattacharyya,1964; Flis et al., 1989; Hohmann and Newman, 1990; Lee, 1975, 1981; Leeand Thomas, 1992; Lewis and Lee, 1984; Morrison et al., 1969; Rathor, 1978;Wait and Debroux, 1984), or by treating the time domain convolution di-rectly (Smith et al., 1988). Zaslavsky and Druskin (2010) and Zaslavskyet al. (2011) developed a modelling technique based on the rational Krylovsubspace projection approach that can model the response of three dimen-sional distribution of chargeable material. However, their technique stillrequired a frequency-to-time domain transformation.Techniques making use of a frequency-to-time domain transformationscan easily incorporate dispersive conductivities, but are only efficient whenmodelling sources containing a small number of frequencies. Accuratelymodelling sources containing a broad band of frequency content, such as thesquare waves commonly used in time domain electromagnetic or inducedpolarization experiments, requires the solution of Maxwell’s equations at a634.2. Convolutionary Ohm’s lawlarge number of frequencies. For example, Newman et al. (1986) and Fliset al. (1989) reported requiring between 20 and 50 frequencies to accuratelymodel a step-off response. If iterative methods are employed, the response ateach frequency can be calculated in parallel. This makes frequency domainmethods efficient at solving problems with a limited number of sources,however, they quickly become computationally limiting as the number ofsources grows (da Silva et al., 2012; Streich, 2009).In this chapter, a technique for modelling the electromagnetic responseof a three dimensional distribution of chargeable material is developed. Thetechnique models the convolutionary nature of Ohm’s law directly, avoidingthe need for frequency to time domain transformations. The technique istested by comparing results to analytic half-space responses and is thendemonstrated for a three-dimensional example.4.2 Convolutionary Ohm’s lawMultiplication in the frequency domain results in a convolution when trans-formed to the time domain. For dispersive conductivities, Ohm’s law in thetime domain is therefore given by the following convolution,~j(t) =∞∫−∞σ(t− τ)~e(τ)dτ (4.1)where σ(t) is the impulse response of the conductivity. Due to the causalnature of the system, σ(t) must have the formσ(t) =0 : t < 0σ∞ : t = 0−σˆ(t) : t > 0(4.2)which can also be written asσ(t) = σ∞δ(t)− σˆ(t)u(t) (4.3)644.3. The impulse response of the Cole-Cole modelwhere u(t) is the Heaviside step function. Furthermore, assuming that theelectric field is zero for t < 0 results in the convolutionary form of Ohm’slaw for chargeable media~j(t) = σ∞~e(t)−t∫0σˆ(t− τ)~e(τ)dτ (4.4)To evaluate Ohm’s law for chargeable materials, an understanding of theimpulse response σ(t) is essential.4.3 The impulse response of the Cole-Cole modelA common description of the IP effect impulse response is given by thewell-known Cole-Cole relation (Pelton et al., 1978)σ(ω) = σ∞ − σ∞η1 + (1− η)(iωτ)c (4.5)The evaluation of Equation 4.4 requires the conductivity impulse response,σ(t). Analytic expressions for σ(t) only exist for c = 1 and c = 0.5 (Smithet al., 1988).For the Debye model (c = 1), the frequency dependent conductivity isgiven byσ(ω) = σ∞ − σ∞η1 + (1− η)(iωτ) (4.6)or by rearranging,σ(ω) = σ∞ − σ∞η(1− η)τ11(1−η)τ + iω(4.7)Taking the inverse Laplace transform and making use of the identity (Abramowitzand Stegun, 1964, equation 29.3.8),L−1[1a+ s]= e−atu(t) (4.8)654.3. The impulse response of the Cole-Cole model10-310-210-1100Time (s)10-1010-910-810-710-610-510-410-310-210-1100101ˆσ(t) (S/m)η=0.2η=0.4η=0.6η=0.8(a) Variable η, τ = 0.1s10-310-210-1100101102Time (s)10-1010-910-810-710-610-510-410-310-210-1100ˆσ(t) (S/m)τ=101τ=100τ=10−1τ=10−2(b) Variable τ , η = 0.1Figure 4.1: The impulse response of the Debye model (Equation 4.9) fordifferent values of η (a) and τ (b). In both figures σ∞ = 0.1S/mthe Debye model in the time domain is thenσ(t) = σ∞δ(t)− σ∞ητ(1− η)e− tτ(1−η)u(t) (4.9)Figure 4.1 shows the impulse response of the Debye model for variousvalues of η and τ . The early time limit of this function is well defined, withσˆ(0) =σ∞ητ(1− η) (4.10)A similar process can be used to obtain the impulse response of the Cole-Cole model with c=0.5. By making use of the identity (Abramowitz andStegun, 1964, equation 29.3.37)L−1[1a+√s]=(1√pit− aea2t erfc(a√t))u(t) (4.11)the impulse response of the Cole-Cole model for c = 0.5 is thenσ(t) = σ∞δ(t)− σ∞ηb(1√pit− beb2t erfc(b√t))u(t) (4.12)where b = 1(1−η)τ1/2 . Figure 4.2 shows the impulse response of the Cole-Cole model for c = 0.5, with various values of η and τ . In this case, the664.3. The impulse response of the Cole-Cole model10-410-310-210-1100101102103104105Time (s)10-1010-910-810-710-610-510-410-310-210-1100101102ˆσ(t) (S/m)η=0.2η=0.4η=0.6η=0.8(a) Variable η, τ = 0.1s10-410-310-210-1100101102103104105Time (s)10-1010-910-810-710-610-510-410-310-210-1100101ˆσ(t) (S/m)τ=101τ=100τ=10−1τ=10−2(b) Variable τ , η = 0.1Figure 4.2: The impulse response of the Cole-Cole model with c = 0.5(Equation 4.12) for different values of η (a) and τ (b). In both figuresσ∞ = 0.1S/mconductivity impulse response goes to ∞ as t goes to 0. However, at earlytimes (when t τ) the conductivity is approximated byσˆ(t) ≈ σ∞ηb(1√pit− b)(4.13)A comparison of this approximation to the true expression is pictured infigure 4.3.In order to consider Cole-Cole models with other values of c, the impulseresponse must be approximated numerically. Here, this approximation isachieved by transforming the Cole-Cole model to the time domain by ap-plying digital filters (Anderson, 1983; Guptasarma, 1982). An example ofthe results of this transformation is shown in Figure 4.4. The transform pro-vides an excellent approximation of the true impulse response, with relativeerrors below 1%.The resulting impulse responses for different values of c are shown inFigure 4.5. As with the c = 0.5 case, the impulse responses grow to ∞ as tgoes to 0. It is observed that the impulse response at early times (t τ)goes asσˆ(t) ≈ mtc−1 + d (4.14)674.3. The impulse response of the Cole-Cole model10-1010-810-610-410-2100102104Time (s)10-710-510-310-1101103105ˆσ(t) (S/m)True ˆσ(t)Approx. ˆσ(t)τFigure 4.3: The true impulse response of the Cole-Cole model with c =0.5 (Equation 4.12, black line) compared to the early time approximation(Equation 4.13, red line). The two results agree when t τwhere m and d are constants. For the special cases where c = 1 and c = 0.5the values of m and d can be determined analytically. For the Debye model(c = 1),m1 = 0 (4.15a)d1 =σ∞ητ(1− η) (4.15b)When c = 0.5, these values arem0.5 =σ∞η√piτ(1− η) (4.16a)d0.5 =σ∞ητ(1− η)2 (4.16b)For other values of c, the values of m and d can be determined by fittingthe early time behaviour of σˆ(t). The true behaviour of σˆ(t), as well as theearly time approximation from Equation 4.14 for different values of c areshown in Figure 4.6.684.3. The impulse response of the Cole-Cole model10-610-510-410-310-210-1100101102Time (s)10-510-410-310-210-1100101102103ˆσ(t) (S/m)True ˆσ(t)Transformed ˆσ(t)(a)10-610-510-410-310-210-1100101102Time (s)10-810-710-610-510-410-310-2Relative Error (S/m)(b)Figure 4.4: (a) The true (black line, Equation 4.12) and transformed (redline) impulse responses for a Cole-Cole model with c = 0.5. (b) The relativeerror of the transformed impulse shown above. The error is less than 1%over this time interval.694.4. Discretization of Ohm’s law in time10-510-310-1101103105Time (s)10-1010-810-610-410-2100102ˆσ(t) (S/m)c=0.20c=0.40c=0.60c=0.80Figure 4.5: The impulse response of the Cole-Cole model for different valuesof c. Each response has σ∞ = 0.1S/m, η = 0.1, and τ = 1s.In the next section these time-domain impulse responses will be usedto numerically evaluate convolutionary Ohm’s law in the form of Equation4.4.4.4 Discretization of Ohm’s law in timeIn this section I derive an approach to numerically approximate the convo-lutionary form of Ohm’s law for materials exhibiting a Cole-Cole frequencydependence. The goal is to approximate Equation 4.4 in terms of electricfields given at a set of discrete times. In this way, a version of Ohm’s lawcan then be included in an explicit time-stepping electromagnetic forwardmodelling routine.The time axis will be discretized intoN+1 discrete points [t(0), t(1), ..., t(N)]withδt(n) = t(n) − t(n−1) (4.17)Splitting the integral Equation 4.4 into time steps, gives704.4. Discretization of Ohm’s law in time10-1010-810-610-410-2100102104Time (s)10-810-710-610-510-410-310-210-1100101ˆσ(t) (S/m)True ˆσ(t) Approx. ˆσ(t)(a) c = 0.810-1010-810-610-410-2100102104Time (s)10-710-610-510-410-310-210-1100101102103ˆσ(t) (S/m)True ˆσ(t) Approx. ˆσ(t)(b) c = 0.610-1010-810-610-410-2100102104Time (s)10-710-610-510-410-310-210-1100101102103104105ˆσ(t) (S/m)True ˆσ(t) Approx. ˆσ(t)(c) c = 0.410-1010-810-610-410-2100102104Time (s)10-710-610-510-410-310-210-1100101102103104105106107ˆσ(t) (S/m)True ˆσ(t) Approx. ˆσ(t)(d) c = 0.2Figure 4.6: Transformed Cole-Cole impulse responses compared to the lowfrequency approximation from Equation 4.14. The coefficients in these ap-proximations are calculated by fitting the early time Cole-Cole response.The value of τ is denoted by the dotted grey line.714.4. Discretization of Ohm’s law in time~j(t(n)) = σ∞~e(t(n))−t(1)∫0σˆ(t(n) − τ)~e(τ)dτ−t(2)∫t(1)σˆ(t(n) − τ)~e(τ)dτ− · · ·−t(n)∫t(n−1)σˆ(t(n) − τ)~e(τ)dτ (4.18)τt(n−1)t(n−2)t(n−3)t(n−4)t(n−5)t(n−6)t(n−7)t(n−8) t(n)σˆ(τ − t(n))~e(τ)Figure 4.7: Cartoon depiction of the terms involved in the evaluation ofEquation 4.18.The first n − 1 terms in this expression are easily approximated usingthe trapezoid rule,t(k)∫t(k−1)σˆ(t(n) − τ)~e(τ)dτ ≈ δt(k)2(σˆ(t(n) − t(k−1))~e(k−1) + σˆ(t(n) − t(k))~e(k))(4.19)The final term requires additional consideration as σˆ(0) is undefinedwhen c 6= 1.724.4. Discretization of Ohm’s law in timeAssuming that the electric field is linear over the time interval δt,~e(t) ≈ t(n) − tδt(n)~e(n−1) +t− t(n−1)δt(n)~e(n) (4.20)allows the final integral in Equation 4.18 to be written ast(n)∫t(n−1)σˆ(t(n) − τ)~e(τ)dτ=~e(n−1)δt(n)t(n)∫t(n−1)σˆ(t(n) − τ)(t(n) − τ)dτ + ~e(n)δt(n)t(n)∫t(n−1)σˆ(t(n) − τ)(τ − t(n−1))dτ= 1δt(n)δt(n)∫0τ σˆ(τ)dτ~e(n−1) + δt(n)∫0σˆ(τ)dτ − 1δt(n)δt(n)∫0τ σˆ(τ)dτ~e(n)(4.21)As long as δt τ , the Cole-Cole impulse response can be approximatedby Equation 4.14 on the interval from 0 to δt. If this is the case, thenδt∫0σˆ(τ)dτ =mδtcc+ dδt (4.22a)δt∫0τ σˆ(τ)dτ =mδtc+1c+ 1+dδt22(4.22b)andt(n)∫t(n−1)σˆ(t(n) − τ)~e(τ)dτ ≈(mδtcc+ 1+dδt2)~e(n−1) +(mδtcc(c+ 1)+dδt2)~e(n)(4.23a)≈κ(δt)~e(n−1) + γ(δt)~e(n) (4.23b)734.4. Discretization of Ohm’s law in timewhere the coefficients κ(δt) and γ(δt) are defined byκ(δt) =mδtcc+ 1+dδt2(4.24a)γ(δt) =mδtcc(c+ 1)+dδt2(4.24b)Combining Equations 4.18, 4.19 and 4.23b results in the semi-discretizedequation~j(n) = σ∞~e(n) − δt(1)2(σˆ(t(n))~e(0) + σˆ(t(n) − t(1))~e(1))− δt(2)2(σˆ(t(n) − t(1))~e(1) + σˆ(t(n) − t(2))~e(2))− · · ·− δt(n−1)2(σˆ(t(n) − t(n−2))~e(n−2) + σˆ(t(n) − t(n−1))~e(n−1))− κ(δt(n))~e(n−1) − γ(δt(n))~e(n) (4.25)or~j(n) =(σ∞ − γ(δt(n)))~e(n) −~j(n−1)p (4.26)where~j(n−1)p =n−1∑k=1δt(k)2[σˆ(t(n) − t(k−1))~e(k−1)) + σˆ(t(n) − t(k))~e(k)]+ κ(δt(n))~e(n−1)(4.27)By using Equations 4.26 and 4.27 to approximate Ohm’s law for a Cole-Cole media in combination with Maxwell’s equations, electromagnetic ex-periments in the presence of chargeable materials can be simulated. To doso, Maxwell’s equations must be discretized in time and space.744.5. Discretization of Maxwell’s equations4.5 Discretization of Maxwell’s equationsMaxwell’s equations with a convolutionary Ohm’s law are~∇× ~e+ ∂~b∂t= 0 (4.28a)~∇× 1µ~b−~j = ~js (4.28b)~j = σ∞ −t∫0σˆ(t− τ)~e(τ)dτ (4.28c)Applying a backward Euler discretization in time,∂f∂t(t(n+1))≈ f(n+1) − f (n)δt(n+1)(4.29)and incorporating the results from Section 4.4 gives the following semi-discretized set of differential equations~∇× ~e(n+1) + ~e(n+1) − ~e(n)δt(n+1)= 0 (4.30a)~∇× 1µ~b(n+1) −~j(n+1) = ~j(n+1)s (4.30b)~j(n+1) =(σ∞ − γ(δt(n+1)))~e(n+1) −~j(n)p (4.30c)The earth is now discretized unto an orthogonal tensor mesh where phys-ical properties are located at cell centres, electric fields and current densitieslocated on cell edges, and magnetic fields are located on cell faces.Let b, e and j be grid functions that are the staggered discretization of~b,~eand~j. Using a standard staggered discretization of the differential operators,the following system of discretized differential equations is obtainedCe(n+1) +b(n+1) − b(n)δt(n+1)= 0 (4.31a)CMfµ−1b(n+1) −Mej(n+1) = Mej(n+1)s (4.31b)Mej(n+1) = MeAe(n+1) − j(n)p (4.31c)754.5. Discretization of Maxwell’s equationsIn these equations, C and C> are the discrete curl operators which mapfrom edges to faces or faces to edges, respectively. The vector jp is thediscrete representation of ~jp, defined byj(n)p =n∑k=1δt(k)2[Meσˆ(n,k−1)e(k−1)) + Meσˆ(n,k)e(k)]+ Meκe(n) (4.32)The matrices Me, Mfµ−1 , MeA, Meσˆ(n,k), and Meκ are mass matrices, givenbyMe = diag(Ace>v)(4.33a)Mfµ−1 = diag(Acf> (v ◦ µ−1)) (4.33b)MeA = diag(Ace> (v ◦ (σ∞ − γ)))(4.33c)Meσˆ(n,k) = diag(Ace>(v ◦ σˆ(t(n) − t(k))))(4.33d)Meκ = diag(Ace> (v ◦ κ))(4.33e)In these definitions, Ace and Acf are averaging operators, which map fromcell edges to cell centres and cell faces to cell centres, respectively. Thevector v contains the cell volumes, and ◦ denotes the Hadamard product.Details of the formation of the discrete differential operators and the massmatrices are given in Appendices A and B.In order to solve the system of equations 4.31, we first solve Equation4.31c for the current density, j(n+1)j(n+1) = Me−1MeAe(n+1) −Me−1j(n)p (4.34)Substituting this expression into Equation 4.31b, and solving 4.31b for e(n+1)givese(n+1) = MeA−1CMfµ−1b(n+1) + MeA−1j(n)p −MeA−1Mej(n+1)s (4.35)Eliminating e(n+1) in Equation 4.31a and rearranging, results in a linear764.6. Implementation and validationsystem relating the magnetic fields at the time t(n+1) to the magnetic fieldsat time t(n), physical properties, and current densities(Mfµ−1CMeA−1CMfµ−1 +1δt(n+1)Mfµ−1)b(n+1) =1δt(n+1)Mfµ−1b(n)−Mfµ−1CMeA−1j(n)p + Mfµ−1CMeA−1Mej(n+1)s(4.36)Prior to substitution, Equation 4.31a was multiplied by Mfµ−1 to ensure thatthe forward modelling matrix (the matrix in front of b(n+1)) is symmetric.Given the fields at previous time steps, Equation 4.36 is solved for the mag-netic fields at t(n+1) and Equations 4.34 and 4.35 are then used to calculatethe electric fields and current densities at time t(n+1).4.6 Implementation and validationThe forward modelling of the magnetic and electric fields was implementedin the Python programing language. The resulting algorithm is shown onpage 79. The linear system in Equation 4.36 is factored (Algorithm 4.1,line 4) and solved (Algorithm 4.1, line 11) using the direct method Mul-tifrontal Massively Parallel Solver (MUMPS) developed by the CERFACSgroup (Amestoy et al., 2001).The application of a direct solver has distinct strengths and weaknesseswhen compared to iterative methods. While the forward modelling matrixis sparse and requires relatively little memory, factoring this matrix comeswith significant memory requirements. However, once the matrix has beenfactored, the linear system can be solved quickly. The strength of directmethods comes when we consider problems with multiple transmitters anda uniform or semi-uniform discretization in time. The forward modelling ma-trix does not change as long as the time step, δt, remains constant, thereforeonly a single factorization is necessary as it can be reused for all transmittersand over multiple time steps. The use of direct methods for time domainelectromagnetic problems and its comparison to iterative methods is pre-774.6. Implementation and validationsented in detail in Oldenburg et al. (2013).The performance of this modelling algorithm can be assessed using twodifferent metrics; the amount of time it takes to perform the simulation, andthe maximum memory required.The execution time of Algorithm 4.1 is dominated by three operations.The first is the factorization of the forward modelling matrix (Line 4). Fora given physical property model, this matrix changes only when the size ofthe time step (δt) changes. If multiple steps are taken with the same valueof δt, the factorization of this matrix can be stored and reused.The second time consuming operation is the solution of the factoredlinear system (Line 11). Since simulating multiple transmitters results inmultiple right-hand-sides, this operation is carried out for the number oftime steps multiplied by the number of transmitters times. Finally, thecalculation of jp (Lines 6 to 9) is expensive for problems with fine timediscretization. For a problem with Nt time steps, Line 8 is executed Nt ·(Nt − 1) /2 times.There are two significant uses of memory in Algorithm 4.1. The firstcomes from storing the dense factors of the forward modelling matrices, andthe second comes form the storage of the electric fields at all previous timesteps. Storing the electric fields requires the storage of a dense array thathas dimensions of the number of cell edges by the number of time steps bythe number of transmitters.4.6.1 Comparison to analytic solutionsIt is possible to test the accuracy of the simulated time domain responseusing an analytic frequency domain expression for the magnetic fields. Thevertical component of the magnetic field measured at the surface of a uniformhalf-space, a distance r away from a vertical magnetic dipole (Ward andHohmann, 1988), is given byBz = − mµ02pik2r5[9− (9 + 9ikr − 4k2r2 − ik3r3) e−ikr] (4.37)784.6. Implementation and validationAlgorithm 4.1 Convolution based forward modellingRequire: Initial fields b0, e0 and j0Source current density jsProjection matrices QL and QTList of time steps, δtTotal number of time steps, Nt1: for n = 0 to Nt do2: if δt(n) 6= δt(n−1) or n = 0 then3: Form the forward modelling matrix AA←Mfµ−1CMeA−1CMfµ−1 +1δt(n)Mfµ−14: Compute factors of A5: end if6: Initialize jpjp ←Meκe(n−1)7: for k = 1 to n− 1 do8: Calculate the contribution of the kth interval to jpjp ← jp + δt(k)2(Meσˆ(n,k−1)e(k−1)) + Meσˆ(n,k)e(k))9: end for10: Calculate the RHSrhs← 1δt(n+1)Mfµ−1b(n) + Mfµ−1CMeA−1(Mej(n+1)s − jp)11: Solve Ab(n+1) = rhs12: Update ee(n+1) ←MeA−1CMfµ−1b(n+1) + MeA−1(jp −Mej(n+1)s)13: Project fields to receiver locationsdn+1 = QLb(n+1)14: end for15: Project results in timed = QTd794.6. Implementation and validationIn this expression, k = (−iσµ0ω)1/2 is the wavenumber and m is the dipolemoment of the transmitter.The time domain response of a frequency dependent σ is obtained bytransforming the output of this expression to the time domain. The frequency-domain to time-domain transform is accomplished through the use of digitalfilters (Anderson, 1983; Guptasarma, 1982).To test the implementation of the convolution-based forward modelling,the step-off responses of a vertical magnetic dipole source located at thesurface was simulated. The simulations were performed for various uniformhalf-space models, and the results were compared to the results obtainedfrom transforming the output of Equation 4.37.For the tests, the earth was discretized onto a cylindrical mesh. The coreregion or the cylindrical mesh was comprised of twenty five 2.5 m cells in theradial direction, and twelve, 2.5 m cells in the vertical direction. The coreregion was padded with an additional thirty five cells in the up, down andradial directions. These padding cells expanded with a fixed ratio of 1 to 3.The time axis was discretized into 400 segments, with 100 steps taken witheach of four values of δt = (10−5s, 5× 10−5s, 2.5× 10−4s and 1.25× 10−3s.)The receiver was placed 50 m from the transmitter which had a a dipolemoment of 1Am2. Both the transmitter and the receiver were located onthe surface of a uniform chargeable half-space. The half-space had Cole-Coleparameters of σ∞ = 10−2S/m, η = 0.75 and τ = 1s. Simulations were runfor four values of c = (1.0, 0.75, 0.5 and 0.25).The resulting magnetic fields, along with the analytic results are shownin Figure 4.8. An excellent agreement with the transformed analytic resultis obtained for all four of the tested models.The four half-space examples all had very similar execution times andtotal memory requirements. For these small 1D examples, memory require-ments were minimal. Storage of the electric fields at all times required15.2MB, and the factors of the forward modelling matrix used 7.5MB. Theexecution times for factorization, jp calculation, and solve for the exam-ple with c = 0.25 are shown on Table 4.1. The small size of the problemmakes the factorization and solve times almost negligible. Due to the fine804.6. Implementation and validation10-410-310-210-1Time (s)10-2110-2010-1910-1810-1710-1610-1510-1410-13bz (T)Analytic CalculatedNon-Chg. (σ∞) Non-Chg. (σ0)(a) c = 1.0010-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-1410-13bz (T)Analytic CalculatedNon-Chg. (σ∞) Non-Chg. (σ0)(b) c = 0.7510-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)Analytic CalculatedNon-Chg. (σ∞) Non-Chg. (σ0)(c) c = 0.5010-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)Analytic CalculatedNon-Chg. (σ∞) Non-Chg. (σ0)(d) c = 0.25Figure 4.8: A comparison of the magnetic fields simulated using the convo-lution algorithm and the transformed analytic expression. Excellent agree-ment between the pair of decays is obtained for all of the examples.814.6. Implementation and validationdiscretization in time, the jp calculation time is significant.Operation # of Calls Min. Time Max. Time Avg. Time Total TimeFactorization(Line 4)4 0.0552s 0.0624s 0.0574s 0.2296sCalculate jp(Lines 6 to 9)400 0.063s 4.338s 2.032s 812.863sSolve(Line 11)400 0.0019s 0.0034s 0.0020s 0.7854sTotal run time 815.564sTable 4.1: Execution times for the c = 0.25 1D half-space example using theconvolution algorithm.4.6.2 Three dimensional exampleThe convolution method will now be demonstrated on a synthetic threedimensional example. The physical property model for this test consistsof a chargeable block placed in a uniform, non-chargeable background. Theblock has dimensions 100m×100m×80m with the top of the block being 40mbelow the surface. The block has Cole-Cole parameters of σ∞ = 10−1S/m,η = 0.3 and τ = 10−1s. Simulations were performed for two different valuesof c, a Debye model with c = 1, and a Cole-Cole model with c = 0.5. Thebackground has a conductivity of 10−3S/m.Vertical point dipole transmitters with unit dipole moments are placedin a single line passing directly over the centre of the block at 30m intervals,30m above the surface of the model. Receivers are modelled to be co-locatedwith the transmitters, measuring the vertical component of the magneticfield.The conductivity model is discretized onto a regular tensor mesh withcore cells of 20m on a side. The core region consists of 31 × 11 × 20 cells.15 padding cells, growing by 30% with each cell, are then added in eachdirection. This results in a 61 × 41 × 50 cell mesh (125050 cells in total)modelling a 9.3km× 8.9km× 9.1km volume. The core region discretization,transmitter locations, and extent of the chargeable block are shown in Figure4.9. Time is discretized into 160 steps, using four values of δt (1 × 10−5,824.6. Implementation and validation−400 −200 0 200 400Easting (m)−200−1000100200Northing (m)Figure 4.9: Plan view of the layout and core region discretization of thethree-dimensional example. The extent of the chargeable block is shownin grey. The location of the transmitter-receiver pairs are shown by blackcircles.5 × 10−5, 2.5 × 10−4 and 1.25 × 10−3s) and 40 steps being taken for eachvalue.The resulting b-field data is shown in Figures 4.10 and 4.11. Figure 4.10shows the b-field response along the line as a bi-log plot. In this figure, eachline represents the magnitude of the magnetic field at a single time channel.The plotted lines represent 15 different times logarithmically distributedbetween 10−4s and 10−2s. On a bi-log plot, both positive and negativevalues are plotted logarithmically, with the transition zone (in this casebetween ±10−9nT) plotted linearly. Figure 4.11 shows the decay of themagnetic fields recorded directly over the centre of the chargeable block forboth values of c, as well as the decay of the non-chargeable σ∞ model. They-component of the computed current densities at a few times is shown inFigure 4.12.Details of the execution time of the c = 0.5 example are shown in Table4.2. In this 3D example, by including multiple transmitters, both the factor-ization and the solve times have become considerably more significant. Thefactorization and solve steps make up 4.4% and 17.6% of the total runtime,834.6. Implementation and validation−300 −200 −100 0 100 200 300Easting (m)−10−8−10−90 10−910−810−710−610−5bz (nT)(a) Debye dispersion (c = 1)−300 −200 −100 0 100 200 300Easting (m)−10−8−10−90 10−910−810−710−610−5bz (nT)(b) Warberg dispersion (c = 0.5)Figure 4.10: Bi-log plot of the vertical component of the calculated b-fielddata from the model depicted in Figure 4.9 for a Debye dispersion (a) anda Warberg dispersion (b). Values between ±10−9nT are plotted on a linearscale.844.6. Implementation and validation10-410-2Time (s)10-1110-1010-910-810-710-610-5bz (nT)Debye (c=1) Warberg (c=0.5)Non-Chg. (σ∞)Figure 4.11: Time decay of the magnetic field for Debye (Black), Warberg(Red), and non-chargeable (Gray) 3D models. Plotted decays come from thetransmitter receiver pair located directly over the centre of the chargeableblock. Negative values are shown as dashed lines.854.6. Implementation and validation−600 −400 −200 0 200 400 600Easting (m)−700−600−500−400−300−200−1000Depth (m)jy−10−11−10−13−10−15 0 10−15 10−13 10−11(a) t = 1.0× 10−4s−600 −400 −200 0 200 400 600Easting (m)−700−600−500−400−300−200−1000Depth (m)jy−10−11−10−13−10−15 0 10−15 10−13 10−11(b) t = 2.0× 10−4s−600 −400 −200 0 200 400 600Easting (m)−700−600−500−400−300−200−1000Depth (m)jy−10−11−10−13−10−15 0 10−15 10−13 10−11(c) t = 1.0× 10−3sFigure 4.12: The y-component of current densities produced from the trans-mitter located at x=0m and z=30m and the c=0.5 model. (a) Initially, astrong positive response is observed in the block and the background. (b)At later times, the direction of current flow has reversed in the block. Eventhough the direction of current flow in the block has reversed and is decay-ing back to zero, the decay of the remaining positive background responsestill dominates. (c) At late times the response at the receiver is negative,resulting from the decay of the polarization currents in the block.864.7. Conclusionsrespectively. The calculation of jp again makes up the majority of the totalrun time, contributing 75.4% of the total simulation time.Operation # of Calls Min. Time Max. Time Avg. Time Total TimeFactorization(Line 4)4 55.19s 56.03s 55.44s 221.77sCalculate jp(Lines 6 to 9)160 0.19s 46.37s 23.54s 3766.84sSolve(Line 11)160 5.42s 5.58s 5.49s 878.43sTotal run time 4997.18sTable 4.2: Execution times for the c=0.5 3D example.Memory requirements of this example were significant. Storage of theelectric fields at all time steps (a 390504 (# of edges) by 160 (# of time steps)by 21 (# of transmitters) dense array) required 9.77GB, with an additional5.42GB used to store the factors of the forward modelling matrix.Runtime and memory requirements for the c = 1.0 example were verysimilar to those of the c = 0.5 example presented here.4.7 ConclusionsIn this chapter I have developed a new technique for modelling time domainelectromagnetic experiments in the presence of chargeable materials. Thistechnique works directly in the time domain by numerically approximatingthe convolution that appears when transforming Ohm’s law for frequencydependent conductivities back to the time domain.The technique was shown to accurately model the half-space responsesof chargeable models by comparing the simulated results to analytic expres-sions. The technique was then demonstrated on a 3D example involvingmaterial for a pair of values of c.Despite the accuracy of the methods results, its usefulness may be limitedin some situations by the extensive memory requirements and large run timesthat occur for problems involving large number of cells or time steps. Thenext two chapters will develop methods to approximate jp with significantlysmaller run times and memory footprints.87Chapter 5Modelling Debye DispersionsUsing an AuxiliaryDifferential EquationApproach5.1 IntroductionIn the previous chapter a methodology was developed to simulate the elec-tromagnetic response of a three dimensional distribution of chargeable ma-terial. This was accomplished by discretizing the convolutionary form ofOhm’s law and incorporating it into an implicit time stepping routine. Themethod was shown to produce accurate results when compared to analyticsolutions, however, the numerical evaluation of the convolution integral re-quired significant computational resources.For materials exhibiting Debye dispersions (c = 1) Ohm’s law can betransformed to the time domain while avoiding the computationally expen-sive convolution. This property has been applied in various engineeringapplications (Teixeira, 2008). Recent work has focused on the developmentof a method for modelling time-domain wave propagation in media featur-ing dispersive electrical permittivities using an explicit finite difference timedomain (FDTD) scheme. However, with the resulting explicit techniquesdeveloped, stability requirements limit the size of the time step makingthem prohibitively computationally expensive for application to geophysicalproblems.885.2. Auxiliary differential equation approachIn this chapter I develop a similar forward modelling methodology tosimulate the response of Debye materials using an implicit time steppingmethodology better suited for geophysical problems. The technique suffersno significant loss in accuracy when compared to the convolution approachdeveloped previously, while significantly decreasing the required computa-tional resources.5.2 Auxiliary differential equation approachAssuming materials that only exhibit Debye dispersion (Cole-Cole modelwith c = 1) Ohm’s law becomes~J = σ∞(1− η1 + (1− η)(iωτ))~E (5.1)or, by rearranging this expression~J + τ(1− η)iω ~J = σ∞(1− η) ~E + σ∞τ(1− η)iω ~E (5.2)Making use of the Fourier transform pairF(∂f(t)∂t)= iωF (ω) (5.3)equation 5.2 is easily transformed to the time domain,~j + τ(1− η)∂~j∂t= σ∞(1− η)~e+ σ∞τ(1− η)∂~e∂t(5.4)Discretizing this differential equation in time using a backward Eulermethod~j(n+1) + τ(1− η)~j(n+1) −~j(n)δt(n+1)= σ∞(1− η)~e(n+1) + σ∞τ(1− η)~e(n+1) − ~e(n)δt(n+1)(5.5)895.2. Auxiliary differential equation approachand rearranging gives,~j(n+1) = σ∞(1− η)(δt+ τ)δt+ τ(1− η) ~e(n+1) − σ∞τ (1− η)δt+ τ(1− η)~e(n) +τ (1− η)δt+ τ(1− η)~j(n)(5.6)and finally,~j(n+1) = σ∞~e(n+1) − δtσ∞ηδt+ τ(1− η)~e(n+1) − σ∞τ (1− η)δt+ τ(1− η)~e(n) +τ (1− η)δt+ τ(1− η)~j(n)(5.7)As with the convolution approach from Chapter 4, this equation can bewritten in the form~j(n) =(σ∞ − γ(δt(n)))~e(n) −~j(n−1)p (5.8)withγ(δt(n)) =δt(n)σ∞ηδt(n) + τ(1− η) (5.9)and~j(n−1)p =σ∞τ (1− η)δt+ τ(1− η)~e(n−1) − τ (1− η)δt+ τ(1− η)~j(n−1) (5.10)Equations 4.26 and 5.7 accomplish the same thing. They provide a wayto calculate the current densities at a given time from the past history ofthe system. The convolution approach (Equation 4.26) requires knowledgeof the electric fields at the present time as well as all previous time steps.The ADE approach (Equation 5.7) requires the electric field at the presenttime step, as well as the electric field and current densities at one precedingtime step.To make equation 5.7 look more like a convolution, it can be applied toitself recursively to eliminate ~j(n). This results in an expression that, like905.2. Auxiliary differential equation approachEquation 4.26, includes a sum over all previous electric fields~j(n) =(σ∞ − δtσ∞ηδt+ τ(1− η))~e(n) − δtσ∞ηδt+ τ(1− η)(n−1∑k=1(1− δtδt+ τ(1− η))n−k~e(k))(5.11)Although Equations 5.7 and 5.11 are equivalent, Equation 5.7 is far easierto implement. Its evaluation requires the storage and summation of onlytwo fields.The convolution result (Equation 4.26) for a Debye medium while as-suming a constant δt becomes~j(n) =(σ∞ − δtσ∞η2τ(1− η))~e(n) − δtσ∞ητ(1− η)n−1∑k=1ζn−k~e(k) (5.12)whereζ = e− δt(1−η)τ (5.13)and δt is assumed to be small (δt τ). In this way, ζ can be approximatedby its first order Taylor series. By replacing ζ with its first order Taylorseries expansion about δt = 0 gives~j(n) =(σ∞ − δtσ∞η2τ(1− η))~e(n) − δtσ∞ητ(1− η)n−1∑k=1(1− δtτ(1− η))n−k~e(k)(5.14)Equations 5.11 and 5.14 are very similar in form, but not identical. Thisdifference is to be expected as the derivation of each expression involvedthe application of a different method to discretize in time. The derivationof Equation 5.11 used the backward Euler method (order O(δt)) whereasthe derivation of Equation 5.14 used the trapezoid rule (order O(δt2)). Asexpected, the difference between Equation 5.11 and 5.14 can be easily shownto be O(δt).915.3. Discretization of Maxwell’s equations5.3 Discretization of Maxwell’s equationsThe electric and magnetic fields, the current densities and the physical prop-erties will be discretized in the same way as was described in Section 4.5.This results in the discrete system of differential equationsCe(n+1) +b(n+1) − b(n)δt(n+1)= 0 (5.15a)CMfµ−1b(n+1) −Mej(n+1) = Mej(n+1)s (5.15b)Mej(n+1) = MeAe(n+1) − j(n)p (5.15c)In this case, the discrete form of ~jp isj(n)p = MeEe(n) −MeJ j(n) (5.16)The mass matrices appearing in 5.15 and 5.16 are defined to beMe = diag(Ace>v)(5.17a)Mfµ−1 = diag(Acf> (v ◦ µ−1)) (5.17b)MeA = diag(Ace> (v ◦ (σ∞ − γ)))(5.17c)MeE = diag(Acf>(v ◦ σ∞τ (1− η)δt+ τ(1− η)))(5.17d)MeJ = diag(Acf>(v ◦ τ (1− η)δt+ τ(1− η)))(5.17e)Eliminating e(n+1) and j(n+1) from Equations 5.15b using Equations5.15c and 5.15a results in the linear system for b(n+1) given by,(Mfµ−1CMeA−1CMfµ−1 +1δt(n+1)Mfµ−1)b(n+1) =1δt(n+1)Mfµ−1b(n)−Mfµ−1CMeA−1j(n)p + Mfµ−1CMeA−1Mej(n+1)s(5.18)925.4. Implementation and validationUpdates to e and j are calculated usinge(n+1) = MeA−1CMfµ−1b(n+1) + MeA−1j(n)p −MeA−1Mej(n+1)s (5.19a)j(n+1) = Me−1MeAe(n+1) −Me−1j(n)p (5.19b)The expressions in 5.18 and 5.19 are the same as those that appeared inChapter 4, but with different definitions of MeA and jp.5.4 Implementation and validationAs with the convolution based method presented in Chapter 4, the algorithmwas implemented in the Python programming language, making use of theMUMPs direct solver to solve the linear system 5.18.The application of this methodology is summarized in Algorithm 5.1.The matrix factorization (Line 4) and system solve (Line 8) operations re-main significant with regard to execution time, but the calculation of jp(Line 6) is now much less expensive.The performance of this algorithm will again be tested using two differ-ent approaches; the simulated responses of a uniform chargeable half-spacecompared the analytic responses of the equivalent model, and the responseof the three-dimensional model shown in Chapter 4 compared to the convo-lution based results.5.4.1 Comparison to 1D analytic solutionsFor the test the earth was discretized onto a cylindrical mesh with a coreregion comprised of 25, 5 m cells in the radial direction, and 10, 5 m cellsin the vertical direction. The core region was padded with an additional 35cells in the up, down and radial directions. The padding cells were expandedwith a fixed ratio of 1.3. The time axis was discretized into 400 segments,with 100 steps taken for each four values of δt = 10−5s, 5×10−5s, 2.5×10−4sand 1.25×10−3s. Both the spatial and temporal discretizations are identicalto those employed in the half-space examples presented in Section 4.6.1.935.4. Implementation and validationAlgorithm 5.1 ADE based forward modellingRequire: Initial fields b0, e0 and j0Source current density jsProjection matrices QL and QTList of time steps, δtTotal number of time steps, Nt1: for n = 0 to Nt do2: if δt(n) 6= δt(n−1) or n = 0 then3: Form the forward modelling matrix AA←Mfµ−1CMeA−1CMfµ−1 +1δt(n)Mfµ−14: Compute factors of A5: end if6: Calculate jpjp ←MeEe(n−1) −MeJ j(n−1)7: Calculate the RHSrhs← 1δt(n+1)Mfµ−1b(n) + Mfµ−1CMeA−1(Mej(n+1)s − jp)8: Solve Ab(n+1) = rhs9: Update ee(n+1) ←MeA−1CMfµ−1b(n+1) + MeA−1(jp −Mej(n+1)s)10: Project fields to receiver locationsdn+1 = QLb(n+1)11: end for945.4. Implementation and validationThe simulated receiver was placed 50 m from a vertical magnetic dipoletransmitter with a dipole moment of 1Am2. Both the transmitter and thereceiver were placed on the surface of a uniform chargeable half-space. Thetest was run for four different sets of Debye parameters, with σ∞ being either10−2S/m or 1S/m and τ being either 10−2s or 1s. For all tests η was set at0.5.The resulting magnetic fields, along with the analytic results are shownin Figure 5.1. Excellent agreement with the transformed analytic result isobtained for all four of the tested models.10-410-310-210-1Time (s)10-2110-2010-1910-1810-1710-1610-1510-1410-13bz (T)Analytic CalculatedNon-Chg. (σ0) Non-Chg. (σ∞)(a) σ∞ = 10−2S/m, η = 0.5, τ = 1s10-410-310-210-1Time (s)10-2210-2110-2010-1910-1810-1710-1610-1510-1410-13bz (T)Analytic CalculatedNon-Chg. (σ0) Non-Chg. (σ∞)(b) σ∞ = 10−2S/m, η = 0.5, τ = 10−2s10-410-310-210-1Time (s)10-1910-1810-1710-1610-1510-1410-1310-12bz (T)Analytic CalculatedNon-Chg. (σ0) Non-Chg. (σ∞)(c) σ∞ = 1S/m, η = 0.5, τ = 1s10-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-1410-1310-12bz (T)Analytic CalculatedNon-Chg. (σ0) Non-Chg. (σ∞)(d) σ∞ = 1S/m, η = 0.5, τ = 10−2sFigure 5.1: A comparison of the magnetic fields simulated using the ADEalgorithm and the transformed analytic expression. All of the models areuniform Debye (c = 1) half-spaces, with Cole-Cole parameters as shownin the labels. The analytic responses of non-chargeable half-spaces withconductivities of σ0 and σ∞ are also shown for comparison.955.4. Implementation and validationOperation # of Calls Min. Time Max. Time Avg. Time Total TimeFactorization(Line 4)4 0.0553s 0.0632s 0.0576s 0.2305sCalculate jp(Lines 6)400 0.0001s 0.0027s 0.0001s 0.0553sSolve(Line 8)400 0.0018s 0.0031s 0.0020s 0.7887sTotal run time 2.0537sTable 5.1: Execution times for first 1D half-space example (Figure 5.1a)using the ADE algorithm.Details of the time taken to perform the significant steps of Algorithm 5.1for the simulation of the response shown in Figure 5.1a are given on Table5.1. As the discretization was the same as those used in the convolutionexample, the factor and solve times remained relatively unchanged fromthose given on Table 4.1.The time required to calculate jp decreased from 812.863s to 0.0553s.This improvement results in a significantly improved total run-time.5.4.2 Comparison to 3D convolution resultTo test the three-dimensional ADE result the calculated response of a 3Dmodel containing a chargeable block is compared to the computed responseobtained using the convolution approach developed in Chapter 4.The model consists of a chargeable block (100m× 100m× 80m) located40m below the surface of a uniform 10−3S/m half-space. The block hasCole-Cole parameters of σ∞ = 10−1S/m, η = 0.3 and τ = 10−1s.The discretizations in space and time are identical to the examples inSection 4.6.2. The spatial mesh has 61×41×50 cells, including a 31×11×20core region of uniform cells that are 20m on a side. Time is discretized into160 steps, using 4 values of δt = 1×10−5, 5×10−5, 2.5×10−4 and 1.25×10−3sand 40 steps are taken for each δt.The modelled bz response using the ADE and convolution algorithmsare shown in Figures 5.2 (the response along the line) and 5.3 (the verticalcomponent of the b-field measured over the centre of the block). A negativeresponse as clearly visible over the chargeable block, and the results of the965.5. Conclusiontwo approaches are in excellent agreement with each other.−300 −200 −100 0 100 200 300Easting (m)−10−8−10−90 10−910−810−710−610−5bz (nT)ADE Conv.Figure 5.2: Bi-log plot of the vertical component of the calculated b-fielddata from the model depicted in Figure 4.9 computed using the ADE (black)and convolution (red) methods. Values between ±10−9nT are plotted lin-early.Details of the execution time of the example run with the ADE algo-rithm are shown in Table 5.2. As the discretization was identical, both thefactorization and solve times are very similar to those seen with the convo-lution algorithm in Section 4.6.2. However, as observed in the 1D examples,the time required to calculate jp has decreased significantly. This results ina speedup of the entire simulation of approximately 75%.The factorization memory requirements remained constant at 5.42GB,while past field storage needed for the calculation of jp decreased to just131Mb. The total memory footprint of the simulation decreased to 5.55GB.This is a 63% reduction in total memory over the simulation run with theconvolution algorithm.5.5 ConclusionThe auxiliary differential equation technique presented in this chapter cansignificantly decrease the computational difficulty when solving the electro-975.5. Conclusion10-410-310-2Time (s)10-1110-1010-910-810-710-610-5bz (nT)ADE Conv.Non-Chg. (σ∞)Figure 5.3: The decay of the vertical component of the magnetic fieldcalculated using the ADE (black) and convolution (red) methods for thetransmitter-receiver pair located at (0,0). Negative values are denoted by adashed line. The computed response of the non-chargeable model (σ∞) isshown in grey for reference.Operation # of Calls Min. Time Max. Time Avg. Time Total TimeFactorization(Line 4)4 54.54s 56.32s 55.37s 221.48sCalculate jp(Line 6)160 0.0004s 0.0034s 0.0005s 0.0824sSolve(Line 8)160 5.45s 5.57s 5.48s 878.86sTotal run time 1174.01sTable 5.2: Execution times for 3D example using the ADE algorithm.985.5. Conclusionmagnetic forward problem when all of the chargeable material being mod-elled is limited to Debye dispersion. The examples presented in this chapterhave shown that this technique is capable of producing results as accurateas those with the convolution algorithm presented in Chapter 4 in a fractionof the time and with a significantly lower memory requirement.The limitation to Cole-Cole models with c = 1 is problematic. It iswidely accepted that the majority of chargeable material of economic interestexhibits dispersion values of c that are much less than 1 (Pelton et al., 1978;Wong, 1979). This shortcoming will be addressed in the next chapter, wherethe ADE algorithm will be adapted to handle a wider variety of dispersionmodels.99Chapter 6Modelling Cole-ColeDispersions Using Pade´Approximations6.1 IntroductionThe auxiliary differential equation methodology developed in Chapter 5showed significant performance improvements over the convolution-basedalgorithm developed in Chapter 4, but it lacked the flexibility to deal witha wide variety of dispersion models.Similar limitations have been addressed by the FDTD community whensimulating the response of materials with various permittivities. For exam-ple, Rekanos and Papadopoulos (2010) and Rekanos (2012) used a rationalseries approximation of frequency dependent physical properties to obtaina higher order auxiliary differential equation for modelling. Furthermore,they went on to develop an explicit time stepping algorithm for simulatingthe response of materials with frequency dependent permittivities.As mentioned previously, explicit techniques are not well suited to geo-physical electromagnetic problems. In this chapter, I will develop an im-plicit forward modelling routine that makes use of Pade´ approximations toapproximate the frequency dependence of the chargeable materials. Thetechnique will be demonstrated on 1D and 3D geophysical examples, andthe results will be compared to those generated with the methods presentedin the previous chapters.1006.2. Rational function approximation of Ohm’s law6.2 Rational function approximation of Ohm’slawFor materials that exhibit a Cole-Cole like dispersion, Ohm’s law in thefrequency domain is given by~J = σ∞(1− η1 + (1− η)(iωτ)c)~E (6.1)which can be rearranged as~J + τ c(1− η)(iω)c ~J = σ∞(1− η) ~E + σ∞τ c(1− η)(iω)c ~E (6.2)For c < 1, the inverse Fourier transformation of this expression to the timedomain will result in fractional derivatives. The fractional derivative of orderc for 0 < c < 1 is defined by∂c∂tcf(t) =1Γ(1− c)∂∂tt∫0f(x)(t− x)cdx (6.3)While it is possible to numerically approximate this expression, the pres-ence of the integration of f over all previous times would lead us to expectlittle difference in performance over the convolution algorithm presented inChapter 4.Rather than treat the fractional derivatives directly, the frequency de-pendent portion of equation 6.2 is approximated by a rational function. Thisis similar to the approach taken by Weedon and Rappaport (1997) whenmodelling the effects of dispersive media on wave propagation in biologicalmaterials.By introducing a characteristic frequency ω0 (Rekanos and Papadopou-los, 2010), Equation 6.2 becomes~J + (τω0)c(1− η)(iωω0)c~J = σ∞(1− η) ~E + σ∞(τω0)c(1− η)(iωω0)c~E(6.4)1016.2. Rational function approximation of Ohm’s lawwhere the frequency dependent term(iωω0)cis now approximated by therational function RN,M(iωω0)(iωω0)c≈ RN,M(iωω0)=N∑n=0Pn(iωω0)n1 +M∑m=1Qm(iωω0)m (6.5)Ohm’s law now becomes(1 +M∑m=1Qm(iωω0)m)~J + (τω0)c(1− η)N∑n=0Pn(iωω0)n~J =σ∞(1− η)(1 +M∑m=1Qm(iωω0)m)~E + σ∞(τω0)c(1− η)N∑n=0Pn(iωω0)n~E(6.6)By definingPk = 0, if k > M (6.7a)Qk = 0, if k > N (6.7b)with Q0 = 1 and K = max (M,N), Equation 6.6 is rewritten asA0 ~J +K∑k=1Akωk0(iω)k ~J = σ∞(1− η)B0 ~E + σ∞(1− η)K∑k=1Bkωk0(iω)k ~E (6.8)whereAk =Qk + (τω0)c(1− η)Pk (6.9a)Bk =Qk + (τω0)cPk (6.9b)Since the powers of iω are now exclusively integer order, this expressioncan easily be transformed to the time-domain yielding the k-order, ordinary1026.2. Rational function approximation of Ohm’s lawdifferential equationA0~j +K∑k=1Akωk0∂k~j∂tk= σ∞(1− η)B0~e+ σ∞(1− η)K∑k=1Bkωk0∂k~e∂tk(6.10)In order to simplify Equation 6.10, and reduce it to a system of firstorder differential equations, the dummy variables ~fi and ~gi are defined suchthat~fi =∂i~j∂ti=∂ ~fi−1∂t(6.11a)~gi =∂i~e∂ti=∂~gi−1∂t(6.11b)Equation 6.10 can now be rewritten as the first order systemA0~j +K∑k=1Akωk0~fk = σ∞(1− η)B0~e+ σ∞(1− η)K∑k=1Bkωk0~gk (6.12)~f1 =∂~j∂t~f2 =∂ ~f1∂t...~fK =∂ ~fK−1∂t~g1 =∂~e∂t~g2 =∂~g1∂t...~gK =∂~gK−1∂tDiscretizing in time using the backward Euler method, yields the semi-1036.2. Rational function approximation of Ohm’s lawdiscretized system of equationsA0~j(n+1) +K∑k=1Akωk0~f(n+1)k = σ∞(1− η)B0~e(n+1) + σ∞(1− η)K∑k=1Bkωk0~g(n+1)k(6.13)~f(n+1)1 =~j(n+1) −~j(n)δt~f(n+1)2 =~f(n+1)1 − ~f (n)1δt...~f(n+1)K =~f(n+1)K−1 − ~f (n)K−1δt~g(n+1)1 =~e(n+1) − ~e(n)δt~g(n+1)2 =~g(n+1)1 − ~g(n)1δt...~g(n+1)K =~g(n+1)K−1 − ~g(n)K−1δtApplying the definitions of ~fk and ~gk recursively, ~f(n+1)k and ~g(n+1)k areeliminated from the left hand sides of the differential equations in 6.13 suchthat ~f(n+1)k and ~g(n+1)k can be written as~f(n+1)k =~j(n+1)δtk−k−1∑i=0~f(n)iδtk−i(6.14a)~g(n+1)k =~e(n+1)δtk−k−1∑i=0~g(n)iδtk−i(6.14b)where ~f(n)0 =~j(n) and ~g(n)0 = ~e(n). Substituting this observation into thefirst equation of 6.13A0~j(n+1)+K∑k=1Akωk0(~j(n+1)δtk−k−1∑i=0~f(n)iδtk−i)= σ∞(1− η)B0~e(n+1) + σ∞(1− η)K∑k=1Bkωk0(~e(n+1)δtk−k−1∑i=0~g(n)iδtk−i)(6.15)1046.2. Rational function approximation of Ohm’s lawand collecting like indices of ~j, ~e, ~f and ~g gives the simplified relation(K∑k=0Akωk0δtk)~j(n+1) −K−1∑k=0(K∑i=k+1Aiωi0δti−k)~f(n)k= σ∞(1− η)(K∑k=0Bkωk0δtk)~e(n+1) − σ∞(1− η)K−1∑k=0(K∑i=k+1Biωi0δti−k)~g(n)k(6.16)Dividing this expression by the coefficient leading ~j(n+1) results in the fa-miliar form~j(n+1) = (σ∞ − γ)~e(n+1) −~j(n)p (6.17)whereγ = σ∞ηK∑i=0Qi(ω0δt)K−iK∑i=0Ai(ω0δt)K−i(6.18)and~j(n)p =K−1∑k=0κ[k]g ~g(n)k −K−1∑k=0κ[k]f~f(n)k (6.19)with the coefficients κg and κf defined byκ[k]g =δtkK∑i=k+1Bi(ω0δt)K−iK∑i=0Ai(ω0δt)K−i(6.20a)κ[k]f =δtkK∑i=k+1Ai(ω0δt)K−iK∑i=0Ai(ω0δt)K−i(6.20b)1056.3. Pade´ approximation of the Cole-Cole modelIn order to simulate electromagnetic fields given this definition of currentflow in a chargeable media that is approximated by a rational function, thecoefficients of the rational function are required. In the following section amethod for calculating these coefficients is presented.6.3 Pade´ approximation of the Cole-Cole modelThe method described in Section 6.2 requires a rational function approxi-mation for(iωω0)cof the formRN,M =N∑n=0Pn(iωω0)n1 +M∑m=1Qm(iωω0)m (6.21)One way of calculating RN,M is to use a Pade´ approximation. A Pade´approximation is a rational function whose Taylor series expansion agreeswith the Taylor series of function it is approximating to the highest degreepossible (Baker and Gammel, 1961; Baker and Graves-Morris, 1996). Tothis end, the Taylor series is used to compute the necessary rational functioncoefficients, or Pade´ coefficients.6.3.1 Calculation of Pade´ coefficientsThe method for calculating the Pade´ coefficients Pn, Qm follows that de-scribed in Press et al. (2007) and employed in Rekanos and Papadopoulos(2010).Let s = iωω0 , then since c < 1, sc is not differentiable at s = 0 and thus itsMaclaurin series does not exist. Instead, the Taylor series expansion abouts = 1 is used. The K order Taylor series expansion of sc about the points = 1 is given byTK =K∑k=0αk (s− 1)k (6.22)1066.3. Pade´ approximation of the Cole-Cole modelwhereαk =(ck)(6.23)To simplify the calculation of the coefficients Pn and Qm we will first com-pute P˜n and Q˜m such thatRN,M =N∑n=0P˜n (s− 1)n1 +M∑m=1Q˜m (s− 1)m(6.24)and then use these values to calculate Pn and Qm.To calculate P˜n and Q˜m, RN,M (s) is equated to the Taylor series ap-proximation of order N + M . Multiplying both sides of the equations bythe value in the denominator of Equation 6.21 gives(N+M∑k=0αk (s− 1)k)(1 +M∑m=1Q˜m(s− 1)m)=N∑n=0P˜n(s− 1)n (6.25)Expanding this expression and equating terms with like powers of (s−1)up to (s− 1)N+M results in a linear system1 0 · · · 0 0 · · · 00 1 · · · 0 −α0 · · · 0....... . ........ . ....0 0 · · · 1 −αN−1 · · · −α00 0 · · · 0 −αN · · · −α1...............0 0 · · · 0 −αN+M−1 · · · −αNP˜0...P˜NQ˜1...Q˜M=α0α1...αN+M (6.26)This system is easily solved for P˜n and Q˜m. Finally, expanding the summa-tions in Equation 6.24 and collecting like powers of s allows Pn and Qm to1076.3. Pade´ approximation of the Cole-Cole modelbe calculated usingPn =N∑k=n(−1)k−n(kn)P˜kM∑k=0(−1)kQ˜k(6.27a)Qm =M∑k=m(−1)k−m( km)P˜kM∑k=0(−1)kQ˜k(6.27b)This section has outlined the method used to calculate the values ofthe coefficients that appear in the rational function RN,M . The followingsections present the implementation of this method when approximating aCole-Cole conductivity spectra.6.3.2 Approximation order and central frequencyThree parameters must be provided when using this algorithm: the orderof the polynomial appearing in the numerator (N) and the denominator(M) of the Pade´ approximation, and the central frequency (ω0) used in theevaluation of Equation 6.17. The values of these parameters impact theperformance (the run-time and the memory requirements) and the accuracyof the modelling algorithm in a variety of ways.Both the runtime and the memory requirements increase as the order ofthe approximation increases. This is clear from the definition of the quantity~j(n)p (Equation 6.19). Evaluating this expression requires a sum from 0 toK− 1 (where K = max(M,N)) which will take longer to evaluate for largervalues of K, and also requires the storage of 2K field derivative vectors(~g(n)k and~f(n)k ) from the previous time step. It is important to note thatthe runtime and the memory requirements vary only with changes in K,and not with N or M individually. Thus, there is no benefit to using anapproximation with differing orders in the numerator and denominator froma purely performance perspective.The ability of the Pade´ approximation to represent sc improves as N1086.3. Pade´ approximation of the Cole-Cole model10-310-210-1100101102103s10-210-1100101102s0.5R2,2(s) R2,1(s) R1,2(s)(a) Approximations with K = 210-310-210-1100101102103s10-1210-1010-810-610-410-2100102Relative ErrorR2,2(s) R2,1(s) R1,2(s)(b) Relative error with K = 210-310-210-1100101102103s10-210-1100101102s0.5R3,3(s) R3,2(s) R2,3(s)(c) Approximations with K = 310-310-210-1100101102103s10-1210-1010-810-610-410-2100102Relative ErrorR3,3(s) R3,2(s) R2,3(s)(d) Relative error with K = 310-310-210-1100101102103s10-210-1100101102s0.5R4,4(s) R4,3(s) R3,4(s)(e) Approximations with K = 410-310-210-1100101102103s10-1210-1010-810-610-410-2100102Relative ErrorR4,4(s) R4,3(s) R3,4(s)(f) Relative error with K = 2Figure 6.1: Pade´ approximations to s0.5 of varying orders from K = 2 toK = 4 (left column) and the relative errors of the approximation (rightcolumn).1096.3. Pade´ approximation of the Cole-Cole modeland M increase. Pade´ approximations of s0.5 for different values of N andM are shown in Figure 6.1. Plots in the left column of Figure 6.1 show thetrue value of sc for values of s from 10−3 to 103 (black lines) along withthe resulting Pade´ approximations using equal numerator and denominatororder (green lines), increased numerator order (red lines) and increased de-nominator order (blue lines). Plots in the right column show the relativeerror of the Pade´ approximations. It is clear that as K increases, the widthof the range of values of s for which the Pade´ approximation provides agood approximation for sc increases. Also, the Pade´ approximation calcu-lated with N = M provides a better approximation over the entire range ofs than either of the cases with the same value of K but with N 6= M .The value of the frequency dependence, c, also plays an important rolein the success of the Pade´ based modelling algorithm. While the value of chas little impact on the range of frequencies over which the Pade´ approx-imation provides an adequate fit to sc (Figure 6.2), the Cole-Cole modelbecomes harder to represent for smaller values of c. Figure 6.3 shows theCole-Cole model for two different values of c along with the conductivityspectrum generated when the frequency dependent part of the Cole-Colemodel,(iωω0)c, is replaced by a Pade´ approximation,σ(ω) = σ∞1− η1 + (1− η)(τω0)cRN,M(iωω0) (6.28)All of the examples were generated with σ∞ = 1 S/m, τ = 1 s, η = 0.5,and ω0 = 1 rad/s. The top pair of figures show results for c = 0.9 and thebottom two for c = 0.5. It is clear from these figures that approximationsof equal order provide a better approximation of a Cole-Cole model with alarger value of c.Finally, the choice of ω0 impacts the range over which the approximationgenerates a good representation of the conductivity spectra. Figure 6.4shows the Cole-Cole conductivity spectra for σ∞ = 1 S/m, τ = 1s, η = 0.5,and c = 0.5 along with the Pade´ approximation of the spectra calculatedwith N = M = 4 and three different values for ω0. The frequency where1106.3. Pade´ approximation of the Cole-Cole model10-310-210-1100101102103s10-210-1100101102s0.25s0.5s0.75(a)10-310-210-1100101102103s10-1410-1210-1010-810-610-410-2100102Relative Errors0.25s0.50s0.75(b)Figure 6.2: (a) True and approximate sc curves for different values of c. Truecurves are shown as solid lines and the N = M = 4 Pade´ approximationsshown as dashed lines of the same colour. (b) Relative error of the Pade´approximations. The value of c has little impact on the frequency rangewhere the approximation is valid.10-310-210-1100101102103Freq. (hz)0.50.60.70.80.91.0ℜ(σ)True N=2 N=4 N=6(a) <(σ) with c = 0.910-310-210-1100101102103Freq. (hz)0.000.050.100.150.200.25ℑ(σ)True N=2 N=4 N=6(b) =(σ) with c = 0.910-310-210-1100101102103Freq. (hz)0.50.60.70.80.91.0ℜ(σ)True N=2 N=4 N=6(c) <(σ) with c = 0.510-310-210-1100101102103Freq. (hz)0.000.020.040.060.080.100.12ℑ(σ)True N=2 N=4 N=6(d) =(σ) with c = 0.5Figure 6.3: Real (a & c) and imaginary (b & d) Cole-Cole dispersions com-pared to Pade´ approximations of the conductivity spectra. Approximationsof the same order do better for the model with the lower value of c. Theapproximations were all calculated with N = M and ωo = 1.1116.3. Pade´ approximation of the Cole-Cole modelthe minimum error occurs is approximately ω0.10-210-1100101102103ω (Rad/s)0.50.60.70.80.91.0ℜ(σ)True ω0=1 ω0=10 ω0=100(a)10-210-1100101102103ω (Rad/s)0.000.020.040.060.080.100.12ℑ(σ)True ω0=1 ω0=10 ω0=100(b)10-210-1100101102103ω (Rad/s)10-610-510-410-310-210-1100Relative Errorω0=1 ω0=10 ω0=100(c)10-210-1100101102103ω (Rad/s)10-510-410-310-210-1100Relative Errorω0=1 ω0=10 ω0=100(d)Figure 6.4: Real (a) and imaginary (b) parts of a Cole-Cole model withσ∞ = 1 S/m, τ = 1 s, η = 0.5, and c = 0.5 compared with the N = M = 4Pade´ approximation with different values of ω0. The relative error of eachof the approximations is shown in (c) and (d).6.3.3 The TEM Response of a Pade´ ModelTo test the accuracy of the time domain response of a Pade´ model throughthe use of a synthetic model requires two approximations; first, a frequencydomain analytic expression for the vertical component of the magnetic fieldmeasured at the surface of a uniform half-space (Equation 4.37), and sec-ond, the approximated frequency dependent conductivities (Equation 6.28).A step-off time domain response can then be calculated form these approx-imations through the use of digital filters (Anderson, 1983; Guptasarma,1982).1126.4. DiscretizationResults of this test are shown in Figures 6.5 and 6.6. The calculatedresponse of a uniform half-space exhibiting a Cole-Cole dispersion is com-pared to the response of Pade´ approximations of different orders. All ofthe examples were run with ω0 = 100Hz. All approximations adequatelyreproduced the early time results, but the lower order approximations didnot provide accurate results at late times. For the larger value of c (Figure6.5) a third or forth order approximation reproduces the response over thistime range. For the smaller value of c (Figure 6.6) more terms are requiredto achieve an acceptable result.In this section I have explained the method used to calculate the coef-ficients of a rational function approximation of frequency dependence, andhave demonstrated its ability to successfully reproduce the time domain re-sponse of a time domain survey. I will now combine the results of Section6.2 with Maxwell’s equations in the time-domain to derive a discrete systemof equations in which the results of this section can be applied.6.4 DiscretizationFollowing the same discretization methodology as was applied in Sections4.5 and 5.3 results in the discrete form of Maxwell’s equationsCe(n+1) +b(n+1) − b(n)δt(n+1)= 0 (6.29a)CMfµ−1b(n+1) −Mej(n+1) = Mej(n+1)s (6.29b)Mej(n+1) = MeAe(n+1) − j(n+1)p (6.29c)with the discrete definition of ~j(n+1)p given byj(n)p =K−1∑k=0(Meκg(k)g(n)k −Meκf (k)f(n)k)(6.30)1136.4. Discretization10-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(a) N=M=110-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(b) N=M=210-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(c) N=M=310-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(d) N=M=4Figure 6.5: Step off response 100m from a vertical magnetic dipole sourceplaced at the surface of a uniform chargeable half-space with Cole-Coleparameters of σ∞ = 10−2 S/m, τ = 0.1s, η = 0.5, and c = 0.75. Theanalytic response of the true model is shown in black and the responseof the Pade´ approximation is shown in red. The Pade´ coefficients wherecalculated using a Taylor series approximation about 100 Hz. Dashed linesindicate a negative response.1146.4. Discretization10-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(a) N=M=110-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(b) N=M=210-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(c) N=M=310-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(d) N=M=410-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)True PadeNon-Chg. (σ∞) Non-Chg. (σ0)(e) N=M=5Figure 6.6: Step off response 100m from a vertical magnetic dipole sourceplaced at the surface of a uniform chargeable half-space with Cole-Coleparameters of σ∞ = 10−2 S/m, τ = 0.1s, η = 0.5, and c = 0.25. Theanalytic response of the true model is shown in black and the response of aPade´ approximation is shown in red. The Pade´ coefficients where calculatedusing a Taylor series approximation about 100 Hz. Dashed lines indicate anegative response.1156.4. Discretizationand where the mass matrices in Equations 6.29 and 6.30 areMe = diag(Ace>v)(6.31a)Mfµ−1 = diag(Acf> (v ◦ µ−1)) (6.31b)MeA = diag(Ace> (v ◦ (σ∞ − γ)))(6.31c)Meκg(k) = diag(Acf>(v ◦ κ[k]g))(6.31d)Meκf (k) = diag(Acf>(v ◦ κ[k]f))(6.31e)with the definitions of γ, κ[k]g and κ[k]f as defined in Equations 6.18 and 6.20.Eliminating e(n+1) and j(n+1) from 6.29 results in the linear system forb(n+1) given by(Mfµ−1CMeA−1CMfµ−1 +1δt(n+1)Mfµ−1)b(n+1) =1δt(n+1)Mfµ−1b(n)−Mfµ−1CMeA−1j(n+1)p + Mfµ−1CMeA−1Mej(n+1)s(6.32)Updates to e and j are calculated usinge(n+1) = MeA−1CMfµ−1b(n+1) + MeA−1j(n+1)p −MeA−1Mej(n+1)s (6.33a)j(n+1) = Me−1MeAe(n+1) −Me−1j(n+1)p (6.33b)The expressions in 6.32 and 6.33 are the same as those that appeared pre-viously, again with different definitions of MeA and jp.1166.5. Implementation and ValidationFinally, g and f are calculated and stored for the next time stepg(n+1)0 = e(n+1) (6.34a)g(n+1)k =e(n+1)δtk−k−1∑i=0g(n)iδtk−i(6.34b)f(n+1)0 = j(n+1) (6.34c)f(n+1)k =j(n+1)δtk−k−1∑i=0f(n)iδtk−i(6.34d)Given the derivation of the discrete form of Maxwell’s equations whichmodel Cole-Cole conductivities approximated by rational functions, the ap-plication of this methodology is presented using 1D and 3D examples in thefollowing section.6.5 Implementation and ValidationAs with the previously presented methods, the Pade´ method was imple-mented in the Python programming language making use of the MUMPsdirect solver to solve the linear system in Equation 6.32.A summary of the application of the Pade´ methodology is outlined inAlgorithm 6.1. Once again, the performance of the Pade´ method is expectedto be worse than the Debye ADE approach due to the loop in the calculationof jp. However, the loop is over a much smaller range of time steps than theconvolution algorithm, so there is still a significant improvement in run-timescompared with the ADE method.6.5.1 Comparison to analytic solutionsThe Pade´ approximation approach is tested using the same two approachesemployed in Section 5.4. First, the response of a uniform chargeable half-space is simulated on a cylindrical mesh and compared to analytic solutions.Second, the three dimensional example first presented in Section 4.6.2 is1176.5. Implementation and ValidationAlgorithm 6.1 Pade´ forward modellingRequire: Initial fields b0, e0 and j0Source current density jsProjection matrices QL and QTList of time steps, # of steps δt, Nt1: Initialize g and fg0 = e0f0 = j02: for k = 1 to K − 1 do3: gk = 0fk = 04: end for5: for n = 0 to Nt do6: if δt(n) 6= δt(n−1) or n = 0 then7: Form the forward modelling matrix AA←Mfµ−1CMeA−1CMfµ−1 +1δt(n)Mfµ−18: Compute factors of A9: end if10: Initialize jpjp ← 011: for k = 0 to K − 1 do12: Calculate the contribution of the kth interval to jpjp ← jp + Meκg(k)g(n)k −Meκf (k)f(n)k13: end for14: Calculate the RHSrhs← 1δt(n+1)Mfµ−1b(n) + Mfµ−1CMeA−1(Mej(n+1)s − jp)15: Solve Ab(n+1) = rhs16: Update fieldse(n+1) ←MeA−1CMfµ−1b(n+1) + MeA−1(jp −Mej(n+1)s)j(n+1) = Me−1MeAe(n+1) −Me−1j(n+1)pg(n+1)0 = en+1f(n+1)0 = jn+117: for k = 1 to K − 1 do18: g(n+1)k =en+1δtk−k−1∑i=0g(n)iδtk−if(n+1)k =jn+1δtk−k−1∑i=0f(n)iδtk−i19: end for20: end for1186.5. Implementation and Validationrepeated and the results will be compared to those obtained using the con-volution approach.To test the Pade´ approach, the earth was discretized onto the samecylindrical mesh used in Section 4.6.1. This mesh had a core region of 25,5 m cells in the radial direction, and 10, 5 m cells in the vertical direction.The core region was padded with an additional 35 cells in the up, down andradial directions with the padding cells expanding with a constant expansionratio of 1.3. The time axis was discretized into 400 segments, with 100 stepstaken with each of four values of δt = 10−5s, 5 × 10−5s, 2.5 × 10−4s and1.25× 10−3s.For the simulations, the receiver was located 50 m from a vertical mag-netic dipole transmitter with a dipole moment of 1 Am2. Both the trans-mitter and the receiver are located on the surface of a uniform chargeablehalf-space. The test was run for three different sets of Cole-Cole parameters,with σ∞ = 10−2 S/m, τ = 1s, η = 0.75. The simulation was run for threedifferent values of c = (0.75, 0.5, 0.25). All three examples used N = M = 5Pade´ approximations with ω0 = 250 rad/s. These parameters were chosenbecause they provided good results when compared to the transformed ana-lytic solution. The same model, discretization and survey layout were usedto generate Figures 4.8b, 4.8c and 4.8d.The simulated bz response, along with the corresponding analytic solu-tion, is shown in Figure 6.7. For the largest value of c (Figure 6.7a) thesimulated response agrees closely with the analytic solution. As the valueof c decreases, the difference between the analytic solution and the simu-lated response grows, particularly at early times. That being said, even forc = 0.25, the simulated response agrees with the analytic solution over themajority of the time range.The slight decrease in accuracy of the simulated response comes witha significant decrease in the run-time and memory requirements. Whereasthe convolution algorithm required the storage of the electric fields at allprevious time steps (a # of time steps by # of edges dense matrix) tofacilitate the calculation of jp, the Pade´ algorithm requires just the storage ofg and f (a 2K by # of edges dense matrix). For this simple 1D example, total1196.5. Implementation and Validation10-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)Analytic CalculatedNon-Chg. (σ∞) Non-Chg. (σ0)(a) c = 0.7510-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)Analytic CalculatedNon-Chg. (σ∞) Non-Chg. (σ0)(b) c = 0.5010-410-310-210-1Time (s)10-2010-1910-1810-1710-1610-1510-14bz (T)Analytic CalculatedNon-Chg. (σ∞) Non-Chg. (σ0)(c) c = 0.25Figure 6.7: A comparison of the magnetic fields simulated using the Pade´algorithm and the transformed analytic expression.1206.5. Implementation and Validationfield storage decreases to 380kB, down from 15.2MB with the convolutionalgorithm. Run-time is also significantly improved (Table 6.1), with thetotal time decreased by over two orders of magnitude.Operation # of Calls Min. Time Max. Time Avg. Time Total TimeFactorization(Line 8)4 0.0553s 0.0583s 0.0564s 0.2256sCalculate jp(Lines 10 to 13)400 0.0007s 0.0422s 0.0012s 0.4640sSolve(Line 15)400 0.0018s 0.0033s 0.0019s 0.7759Total run time 2.6739sTable 6.1: Execution times for the c = 0.25 1D half-space example using thePade´ algorithm.6.5.2 Comparison to 3D convolution resultThe calculated response of a three-dimensional model containing a charge-able block with c = 0.5 will now be compared to the computed responseobtained using the convolution approach developed in Chapter 4.The conductivity model is discretized onto a regular tensor mesh withcore cells 20m on a side. The core region consists of 31× 11× 20 cells, with15 padding cells, growing by 30% with each cell, are added in each direction.This results in a 61 × 41 × 50 cell mesh (125050 cells in total) modelling a9.3km× 8.9km× 9.1km volume. The core region discretization, transmitterlocations, and extent of the chargeable block are shown in Figure 4.9. Thetime axis is discretized into 160 steps, using 4 values of δt = (1 × 10−5,5× 10−5, 2.5× 10−4 and 1.25× 10−3)s with 40 steps for each.The bz responses computed using the two approaches are shown in Fig-ures 6.8 (the response along the line) and 6.9 (the decay of the magneticfield measured over the centre of the block). A negative response is clearlyvisible over the chargeable block, and the results of the two approaches arein excellent agreement with each other.Details of the execution time of the c = 0.5 Pade´ example are shown inTable 6.2. As the discretization was identical, both the factorization andsolve times are very similar to those seen with the convolution algorithm in1216.5. Implementation and Validation−300 −200 −100 0 100 200 300Easting (m)−10−8−10−90 10−910−810−710−610−5bz (nT)ADE Conv.Figure 6.8: Bi-log plot of the vertical component of the calculated b-fielddata from the model depicted in Figure 4.9 computed using the Pade´ (black)and convolution (red) methods. Values between ±10−9nT are plotted lin-early.Section 4.6.2. However, as was seen in the 1D examples, the time requiredto calculate jp has decreased significantly. The entire simulation took ap-proximately 25% of the time required by the convolution algorithm.Total memory requirements also decreased significantly. The factoriza-tion storage requirements remained constant at 5.42GB, while past fieldstorage decreased to 656Mb. The total memory footprint of the simula-tion decreased 6.08GB. This is a 60% reduction in total memory over theconvolution algorithm.Operation # of Calls Min. Time Max. Time Avg. Time Total TimeFactorization(Line 8)4 55.21s 55.89s 55.64s 222.56sCalculate jp(Lines 10 to 13)160 0.0025s 0.0576s 0.0038s 0.6086sSolve(Line 15)160 5.48s 5.52s 5.50s 880.65sTotal run time 1196.96sTable 6.2: Execution times for the c = 0.5 3D example using the Pade´algorithm.1226.5. Implementation and Validation10-410-310-2Time (s)10-1110-1010-910-810-710-610-5bz (nT)ADE Conv.Non-Chg. (σ∞)Figure 6.9: The decay of the vertical component of the magnetic fieldcalculated using the Pade´ (black) and convolution (red) methods for thetransmitter-receiver pair located at (0, 0). Negative values are denoted bya dashed line. The computed response of the non-chargeable model (σ∞) isshown in grey for reference.1236.6. Conclusion6.6 ConclusionThe Pade´ modelling technique developed in this chapter has considerablymore flexibility than the ADE method developed in Chapter 5, while en-joying a similar boost in performance over the convolution based algorithmfrom chapter 4. This improvement does come with a decrease in accuracy(especially for very low values of c and simulations involving many decadesof time). However this slight decrease does not limit the methods potentialusefulness as a tool when considering either very large problems where theconvolution algorithms application will be limited by memory requirementsor small 1D problems where obtaining a solution rapidly is important.124Chapter 7ConclusionsThis thesis began by identifying two questions, the answers to which couldprovide a better understanding of induced polarization effects on inductivesource electromagnetic data. Those questions were:1. Can the effects of chargeability be recognized in frequency domaininductive source electromagnetic data?2. Can we efficiently simulate the electromagnetic response of chargeablematerials in the time domain?The results of my work on these issues have been detailed in the precedingchapters of this thesis. In this chapter I will revisit the conclusions of mywork and provide directions for future research on the subject.7.1 Inductive source induced polarizationIn Chapter 2, a new methodology to identify the presence of chargeablematerial using frequency domain electromagnetics was proposed. Using thesimple asymptotic behaviour of the fields at low frequencies a new datumwas introduced, called the inductive source induced polarization, or ISIP,datum. At low enough frequencies these data are identically zero in thepresence of an entirely non-chargeable conductivity distribution. Thus anynon-zero value in the ISIP datum is a direct indicator of the presence ofchargeable material.Numerical simulations were presented to demonstrate the potential ap-plication of the method to exploration projects. Two different survey geome-tries were considered; a large loop survey and fixed offset Slingram experi-1257.1. Inductive source induced polarizationment. The resulting synthetic ISIP data were able to identify the locationof the chargeable target in both examples.The synthetic tests revealed that for these examples, the magnitude ofthe ISIP data are small and very accurate measurements of the magneticfields would be necessary in order for it to be of use. However, a simpleanalysis of the propagation of error when calculating the data suggests thatthe required level of accuracy should be achievable using modern equipment.The derivation of the ISIP data required that the frequencies employedin their calculation be ‘low enough’. The definition of ‘low enough’ wasexamined using a special case for which an analytic expression exists. Itwas possible to show that value of the induction number, which dependson the product of the frequency and the conductivity of the surroundingmaterial, is the quantity of interest. As long as the induction number ismuch less than unity, the assumptions used when deriving the ISIP datashould hold true. This conclusion was tested on the synthetic examples bysimulating the ISIP data for frequencies of increasing value.Chapter 3 focused on the interpretation of the ISIP data. An inversionscheme was developed to produce three dimensional models of the distribu-tion of chargeable material from ISIP data. This process began by linearizingthe ISIP equations in terms of changes in the real and imaginary parts ofthe conductivity with frequency. Through this process, it was possible toshow that to first order, the ISIP data depend only on the change in thereal part of the conductivity between the chosen frequencies. Thus, the ISIPdata are not sensitive to the overall nature of the conductivity frequency de-pendence. Zeros in the data tell us that there is no frequency dependenceover the frequency range being considered.The resulting sensitivity matrix requires an estimate of the backgroundconductivity distribution in order to invert for a chargeability distribution.This is similar to commonly employed electrical induced polarization inver-sion schemes. Synthetic examples examined the importance of the accuracyof the conductivity model. These results showed that information about theexistence and approximate location of chargeable bodies could be obtainedwith a smooth estimate of the true conductivity or even a representative half1267.2. Three dimensional modelling of IP effects in time domainelectromagnetic dataspace. This lack of sensitivity to the true conductivity is important becauseit is unlikely that this distribution is known in true field situations.Although it does seem as though it would to be difficult to collect, theISIP data is the first proposed methodology for the extraction of induced po-larization information from frequency domain inductive source CSEM data.Theory suggests that it should work in moderately resistive environmentswith currently available equipment. By avoiding the requirements of tra-ditional induced polarization surveying techniques, the ISIP method mayprove to be a valuable tool in geological settings where traditional IP sur-veys are difficult to carry out.7.1.1 Future research opportunitiesAt the time of writing, the ISIP method is still just an idea, demonstratedonly on synthetic cases. Turning this idea into a useful geophysical tool willrequire addition work.The feasibility analysis presented in this work provides a first pass proofof concept of the ISIP method. A more in-depth study will need to be carriedout before moving forward. This study will require input from instrumenta-tion and data acquisition experts to ensure that important possible sourcesof error in the observations have been considered.The only true proof that the ISIP data can be employed as a usefulgeophysical technique would be the successful application of the technique.This test would require a known chargeable target in a suitably resistive en-vironment. Due to the high degree of accuracy required in the observations,a successful trial would likely require cutting edge equipment operated bydata acquisition experts.7.2 Three dimensional modelling of IP effects intime domain electromagnetic dataThe second part of this thesis focused on the forward modelling of timedomain electromagnetic data influenced by the presence of chargeable ma-1277.2. Three dimensional modelling of IP effects in time domainelectromagnetic dataterial. The goal of this work was to develop techniques for modelling theeffects of a frequency dependent conductivity model directly in the time do-main. The resulting methods eliminate the need of Fourier transforms andonly require the solution of real, symmetric, positive definite systems ratherthan complex, non-Hermitian systems that are present in frequency-domainmodelling.In Chapter 4, Ohm’s law, with a frequency dependent conductivity, wastransformed to the time domain. This results in the convolution of theimpulse response of the conductivity and the electric fields at all previoustimes. The evaluation of the convolution was numerically approximated,and this approximation was used within a finite volume discretization inspace and a backward Euler discretization in time. The methodology wasverified by comparing its results to the analytic solution for the responsefor a chargeable half-space for various Cole-Cole dispersions. Finally, themethodology was used on a three dimensional example of an airborne timedomain TEM survey.The analytic comparisons showed that the method produced accurateresults, however, its usefulness may be limited in some situations by theextensive memory requirements and large run times for problems involvinga large number of cells or time steps. Run times may be improved througha better implementation, but the large memory footprint is impossible toavoid. These potential limitations motivated the auxiliary differential equa-tion approach presented in the Chapters 5 and 6.When conductivity spectra are limited to the Debye model (Cole-Colemodel with a frequency dependence of c = 1) it is possible to manipulateOhm’s law so that it can be transformed to the time domain analytically.This results in a first order auxiliary ordinary differential equation that canbe considered along with the remainder of Maxwell’s equations to simulatethe response of these materials. This approach was considered in Chapter5. The resulting algorithm requires the electric fields and current densitiespresent at only one previous time step to be stored resulting in a significantreduction in the both the run time and the memory footprint. This algo-rithm was compared to the analytic response of a chargeable half space and1287.2. Three dimensional modelling of IP effects in time domainelectromagnetic dataalso compared to the response of the three dimensional model calculatedusing the convolution algorithm. These tests showed that the method wasable to produce results in a fraction of the time taken using the previousapproach with no significant loss of accuracy.Unfortunately, the restriction of the method to Debye chargeabilitieslimits the application of this approach in real world situations. Existingliterature (Pelton et al., 1978; Wong, 1979) widely agrees that the majorityof chargeable material of economic interest exhibits dispersions that arebetter represented by Cole-Cole models with values of c that are much lessthan 1.Approximating the frequency dependence of a Cole-Cole model with arational function again also allows Ohm’s law to be analytically transformedto the time domain. The result is an auxiliary differential equation withorder greater than one. This approach was taken in Chapter 6. The coeffi-cients in the rational function approximations were calculated with a Pade´approximation.The resulting algorithm requires the electric fields and current densitydistributions for a previous number of time steps equal to the order of therational function employed. The Pade´ modelling technique developed in thisChapter has considerably more flexibility than the ADE method developedin Chapter 5, while enjoying a similar boost in performance over the convo-lution based algorithm from chapter 4. This improvement does come witha decrease in accuracy (especially for very low values of c and simulationsinvolving many decades of time). However, this slight decrease does notlimit the potential usefulness of this method as a tool when considering ei-ther very large problems where the convolution algorithms application willbe limited by memory requirements, or small 1D problems where obtaininga solution rapidly is important.7.2.1 Future research opportunitiesThe forward modelling methodologies included in this thesis have focusedon the development of modelling tools for synthetic examples. There is a1297.2. Three dimensional modelling of IP effects in time domainelectromagnetic datawide variety of the potential applications for these methods. The modellingtools can be used to provide a better understanding of observations of IPeffects in TEM data. For example, a suite of simulations for varying phys-ical property distributions could provide an understanding of the extent orphysical property range that gives rise to observations of negative transientsin central loop survey data. These methods can also be used to provide abetter understanding of the manifestation of IP effects in non-centre loopTEM survey designs. The simulated results could also be used to test newideas on how to recognize or interpret IP effects.The application of this work is not limited to inductive source electro-magnetics. It is also a valuable tool to simulate the electromagnetic couplingeffects in traditional time domain induced polarization and magnetic inducedpolarization surveys. A simple example of this application was presented inMarchant et al. (2014) and has been included in Appendix C of this the-sis. Full waveform simulations of the response of these surveys will openthe door for better survey design, waveform selection and choice of integra-tion windows in order to maximize the value of IP data for an explorationproject.The material in this thesis restricted itself to conductivity spectra de-scribed by the Cole-Cole model. This choice was made as a result of thepredominance of the Cole-Cole model in geophysical literature. Other mod-els have been suggested that are parameterized in a way that better de-scribes the physical characteristics of the rock (for example, the GEMTIPmodel proposed in Zhdanov (2008)). Modelling dispersions not described bya Cole-Cole model may be desirable in some situations. Nevertheless, themethodologies and ideas presented here can likely form a good foundationfor how such alterations can be made computationally tractable.Finally, a forward modelling routine opens the door for inverse modelling.Once a sound understanding of the sensitivity of TEM responses to featuresof the dispersion curve is obtained, the modelling strategies developed in thiswork could be used to formulate an inversion routine to attempt to recoverinformation. Such work would look to recover the distribution of chargeablematerial, or try to estimate properties of the frequency dependence with the1307.2. Three dimensional modelling of IP effects in time domainelectromagnetic datagoal of mineral discrimination. Preliminary work on this subject was startedin Marchant et al. (2013b) and Kang et al. (2014).131BibliographyAbramowitz, M. and Stegun, I. A. (1964). Handbook of mathematical func-tions with formulas, graphs, and mathematical tables. Dover Publications.Amestoy, P., Duff, I., L’Excellent, J., and Koster, J. (2001). A FullyAsynchronous Multifrontal Solver Using Distributed Dynamic Schedul-ing. SIAM Journal on Matrix Analysis and Applications, 23(1):15–41.Anderson, W. L. (1983). Fourier cosine and sine transforms using lagged con-volutions in double-precision (Subprograms DLAGFO/DLAGF1). Tech-nical report, U.S. Geological Survey.Ascher, U. M. and Petzold, L. R. (1998). Computer Methods for OrdinaryDifferential Equations and Differential-Algebraic Equations. Society forIndustrial and Applied Mathematics, Philadelphia.Baker, G. A. and Gammel, J. (1961). The Pade´ approximant. Journal ofMathematical Analysis and Applications, 2(1):21–30.Baker, G. A. and Graves-Morris, P. (1996). Pade´ Approximants. CambridgeUniversity Press, Cambridge, 2nd edition.Beran, L. and Oldenburg, D. W. (2008). Estimation of Cole-Cole parame-ters from time-domain electromagnetic data. In SEG Technical ProgramExpanded Abstracts, pages 1–5.Bhattacharyya, B. (1964). Electromagnetic fields of a small loop antennaon the surface of a polarizable medium. Geophysics, 29(5):814–831.Bleil, D. F. (1953). Induced Polarization: A Method of GeophysicalProspecting. Geophysics, 18(3):636–661.132BibliographyBo¨rner, F. (2006). Complex conductivity measurements. In GroundwaterGeophysics, pages 119–153. Springer-Verlag, Berlin.Burtman, V., Fu, H., and Zhdanov, M. S. (2014). Experimental Study ofInduced Polarization Effect in Unconventional Reservoir Rocks. Geoma-terials, 04(04):117–128.Chen, J., Haber, E., and Oldenburg, D. W. (2002). Three-dimensional nu-merical modelling and inversion of magnetometric resistivity data. Geo-physical Journal International, 149(3):679–697.Chen, J. and Oldenburg, D. W. (2003). 3-D inversion of magnetic inducedpolarization data. In ASEG Extended Abstracts, number February, pages1–12.Chwala, A., Kingman, J., Stolz, R., Schmelz, M., Zakosarenko, V.,Linzen, S., Bauer, F., Starkloff, M., Meyer, M., and Meyer, H. (2013).Noise characterization of highly sensitive SQUID magnetometer systemsin unshielded environments. Superconductor Science and Technology,26(3):035017.Coggon, J. (1984). New three-point formulas for inductive coupling removialin induced polarization. Geophysics, 49(3):307–309.Coggon, J. H. (1973). A Comparison Of IP Electrode Arrays. Geophysics,38(4):737–761.Cole, K. and Cole, R. (1941). Dispersion and absorption in dielectrics I.Alternating current characteristics. The Journal of Chemical Physics,9:341–351.Connell, D. and Key, K. (2013). A numerical comparison of time andfrequency-domain marine electromagnetic methods for hydrocarbon ex-ploration in shallow water. Geophysical Prospecting, 61(1):187–199.da Silva, N. V., Morgan, J. V., MacGregor, L., and Warner, M. (2012). A fi-nite element multifrontal method for 3D CSEM modeling in the frequencydomain. Geophysics, 77(2):E101–E115.133BibliographyDavydycheva, S., Rykhlinski, N., and Legeydo, P. J. (2006). Electrical-prospecting method for hydrocarbon search using the induced-polarization effect. Geophysics, 71(4):G179–G189.Debye, P. and Huckel, E. (1923). The theory of electrolytes. I. Lowering offreezing point and related phenomena. Physikalische Zeitschrift, 24:185–206.Dey, A. and Morrison, H. F. (1973). Electromagnetic Coupling in Fre-quency and Time-Domain Induced-Polarization Surveys Over a Multilay-ered Earth. Geophysics, 38(2):380–405.Dias, C. A. (2000). Developments in a model to describe lowfrequencyelectrical polarization of rocks. Geophysics, 65(2):437–451.Druskin, V. and Knizhnerman, L. (1994). Spectral approach to solvingthree-dimensional Maxwell’s diffusion equations in the time and frequencydomains. Radio Science, 29(4):937–953.Edwards, R. N., Lee, H., and Nabighian, M. N. (1978). On the theory ofmagnetometric resistivity (MMR) methods. Geophysics, 43(6):1176–1203.Ellis, R. G. and Oldenburg, D. W. (1994). The pole-pole 3-D DC-resistivityinverse problem: a conjugategradient approach. Geophysical Journal In-ternational, 119:187–194.Everett, M. E. (2013). Near-Surface Applied Geophysics. Cambridge Uni-versity Press, Cambridge.Fink, J. B., McAlister, E. O., Sternberg, B. K., Wieduwilt, W. G., andWard, S. H., editors (1990). Induced Polarization: applications and casehistories, volume 4 of Investigations in Geophysics. Society of ExplorationGeophysicists, Tulsa, Oklahoma.Flis, M., Newman, G. A., and Hohmann, G. W. (1989). Induced-polarization effects in time-domain electromagnetic measurements. Geo-physics, 54(4):514–523.134BibliographyFullagar, P., Zhou, B., and Bourne, B. (2000). EM-coupling removal fromtime-domain IP data. Exploration Geophysics, 31:134–139.Fuller, B. D. and Ward, S. H. (1970). Linear System Description of the Elec-trical Parameters of Rocks. IEEE Transactions on Geoscience Electronics,GE-8(1):7–18.Gasperikova, E. (1999). Mapping of Induced Polarization Using NaturalFields. PhD thesis, University of California Berkeley.Gasperikova, E. and Morrison, H. F. (2001). Mapping of induced polariza-tion using natural fields. Geophysics, 66(1):137–147.Guptasarma, D. (1982). Computation of the timedomain response of apolarizable ground. Geophysics, 47(11):1574–1576.Haber, E. and Ascher, U. M. (2001). Fast Finite Volume Simulation of3D Electromagnetic Problems with Highly Discontinuous Coeddicients.SIAM Journal on Scientific Computing, 22(6):1943–1961.Haber, E., Ascher, U. M., Aruliah, D., and Oldenburg, D. W. (2000). FastSimulation of 3D Electromagnetic Problems Using Potentials. Journal ofComputational Physics, 163(1):150–171.Haber, E., Ascher, U. M., and Oldenburg, D. W. (2004). Inversion of 3Delectromagnetic data in frequency and time domain using an inexact al-latonce approach. Geophysics, 69(5):1216–1228.Haber, E. and Heldmann, S. (2007). An octree multigrid method for quasi-static Maxwells equations with highly discontinuous coefficients. Journalof Computational Physics, 223(2):783–796.Hallof, P. (1974). The IP phase measurement and inductive coupling. Geo-physics, 39(October):650–665.Heenan, J., Porter, A., Ntarlagiannis, D., Young, L. Y., Werkema, D. D.,and Slater, L. D. (2013). Sensitivity of the spectral induced polariza-135Bibliographytion method to microbial enhanced oil recovery (MEOR) processes. Geo-physics, 78(5):E261–E269.Hohmann, G. W. (1973). Electromagnetic Coupling Between GroundedWires At The Surface Of A Twolayer Earth. Geophysics, 38(5):854–863.Hohmann, G. W., Kintzinger, P. R., Van Voorhis, G. D., and Ward, S. H.(1970). Evaluation of the measurement of induced electrical polarizationwith an inductive system. Geophysics, 35(5):901–915.Hohmann, G. W., Kintzinger, P. R., Van Voorhis, G. D., and Ward, S. H.(1971). On: ”Evaluation of the Measurement of Induced Electrical Polar-ization with an Inductive System”. Geophysics, 36(2):427–429.Hohmann, G. W. and Newman, G. A. (1990). Transient electromagneticresponses of surficial, polarizable patches. Geophysics, 55(8):1098–1100.Ho¨rdt, A., Blaschek, R., Kemna, A., and Zisser, N. (2007). Hydraulic con-ductivity estimation from induced polarisation data at the field scale theKrauthausen case history. Journal of Applied Geophysics, 62(1):33–46.Hyman, J. and Shashkov, M. (1997). Natural discretizations for the di-vergence, gradient and curl on logically rectangular grids. Computers &Mathematics with Applications, 33(4):81–104.Hyman, J. M. and Shashkov, M. (1999). Mimetic Discretizations forMaxwell’s Equations. Journal of Computational Physics, 151(2):881–909.Kang, S., Oldenburg, D. W., Yang, D., and Marchant, D. (2014). On recov-ering induced polarization information from airborne time domain EMdata. In SEG Technical Program Expanded Abstracts 2014, number 4,pages 1785–1789. Society of Exploration Geophysicists.Kawai, J., Uehara, G., and Kohrin, T. (1999). Three axis SQUID magne-tometer for low-frequency geophysical applications. IEEE Transactionson Magnetics, 35(5):3974–3976.136BibliographyKelley, C. T. (1999). Iterative Methods for Optimization, volume 18. SIAMFrontiers in Applied Mathematics.Kemna, A., Binley, A., Cassiani, G., Niederleithinger, E., Revil, A., Slater,L. D., Williams, K. H., Orozco, A. F., Haegel, F. H., Ho¨rdt, A., Kr-uschwitz, S., Leroux, V., Titov, K., and Zimmermann, E. (2012). Anoverview of the spectral induced polarization method for near-surface ap-plications. Near Surface Geophysics, 10:453–468.Kemna, A., Binley, A., and Slater, L. D. (2004). Crosshole IP imaging forengineering and environmental applications. Geophysics, 69(1):97–107.Kozhevnikov, N. (2012). Fast-decaying inductive IP in frozen ground. Rus-sian Geology and Geophysics, 53(4):406–415.Kozhevnikov, N. and Antonov, E. (2010). Inversion of IP-affected TEM re-sponses of a two-layer earth. Russian Geology and Geophysics, 51(6):708–718.Kunetz, G. (1966). Principles of direct current resistivity prospecting. Ge-bruder Borntraeger, Berlin.Lee, T. J. (1975). Sign Reversals in the Transient Method of ElectricalProspecting (OneLoop Version). Geophysical Prospecting, 23(4):653–662.Lee, T. J. (1981). Transient electromagnetic response of a polarizableground. Geophysics, 46(7):1037–1041.Lee, T. J. and Thomas, L. (1992). The Transient Electromagnetic Responseof a Polarizable Sphere in a Conducting Half Space. Geophysical prospect-ing, 40:541–563.Lewis, R. J. G. and Lee, T. J. (1984). The Detection of Induced Polar-ization with a Transient Electromagnetic System. IEEE Transactions onGeoscience and Remote Sensing, GE-22(1):69–80.Ley-Cooper, Y., Macnae, J., Robb, T., and Vrbancich, J. (2006). Identi-fication of calibration errors in helicopter electromagnetic (HEM) data137Bibliographythrough transform to the altitude-corrected phase-amplitude domain.Geophysics, 71(2):G27.Li, Y. and Oldenburg, D. W. (1994). Inversion of 3-D DC resistivity data us-ing an approximate inverse mapping. Geophysical Journal International,116(3):527–537.Li, Y. and Oldenburg, D. W. (2000). 3-D inversion of induced polarizationdata. Geophysics, 65(6):1931–1945.Lo, B. and Zang, M. (2008). Numerical modeling of ZTEM (airborne AF-MAG) responses to guide exploration strategies. In SEG Technical Pro-gram Expanded Abstracts 2008, volume 27, pages 1098–1102. Society ofExploration Geophysicists.Luo, Y., Zhang, S., and Xiong, B. (2003). Feasibility of Natural SourceInduced Polarization. Chinese Journal of Geophysics, 46(1):169–178.Madden, T. and Mackie, R. (1989). Three-dimensional magnetotelluric mod-elling and inversion. Proceedings of the IEEE, 77(2).Marchant, D., Haber, E., Beran, L., and Oldenburg, D. W. (2012a). 3Dmodeling of IP effects on electromagnetic data in the time domain. InSEG Technical Program Expanded Abstracts, pages 1–5.Marchant, D., Haber, E., and Oldenburg, D. W. (2012b). Inductive sourceinduced polarization. In SEG Technical Program Expanded Abstracts,pages 1–5.Marchant, D., Haber, E., and Oldenburg, D. W. (2012c). Inductive SourceInduced Polarization. In ASEG Extended Abstracts, volume 2012, pages1–4.Marchant, D., Haber, E., and Oldenburg, D. W. (2013a). Inductive sourceinduced polarization. Geophysical Journal International, 192(2):602–612.138BibliographyMarchant, D., Haber, E., and Oldenburg, D. W. (2013b). Recovery of 3DIP distribution from airborne time-domain EM. Proceedings of the 23rdInternational Geophysical Conference and Exhibition, pages 1–4.Marchant, D., Haber, E., and Oldenburg, D. W. (2014). Three-dimensionalmodeling of IP effects in time-domain electromagnetic data. Geophysics,79(6):E303–E314.Marshall, D. J. and Madden, T. R. (1959). Induced Polarization, A Studyof its Causes. Geophysics, 24(4):790–816.Morrison, H. F., Phillips, R., and O’Brien, D. (1969). Quantitative in-terpretation of transient electromagnetic fields over a layered half space.Geophysical prospecting, 14:82–101.Nabighian, M. (1988). Electromagnetic methods in applied geophysics, vol-ume 2. Society of Exploration Geophysicists.Newman, G. A. and Alumbaugh, D. L. (1995). Frequency-domain modellingof airborne electromagnetic responses using staggered finite differences.Geophysical Prospecting, 43(April):1021–1042.Newman, G. A., Hohmann, G. W., and Anderson, W. L. (1986). Transientelectromagnetic response of a three-dimensional body in a layered earth.Geophysics, 51(8):1608–1627.Okay, G., Cosenza, P., Ghorbani, A., Camerlynck, C., Cabrera, J., Florsch,N., and Revil, A. (2013). Localization and characterization of cracks inclay-rocks using frequency and time-domain induced polarization. Geo-physical Prospecting, 61(1):134–152.Oldenburg, D. W., Haber, E., and Shekhtman, R. (2013). Three dimensionalinversion of multisource time domain electromagnetic data. Geophysics,78(1):E47–E57.Oldenburg, D. W. and Li, Y. (1994). Inversion of Induced Polarization Data.Geophysics, 59(9):1327–1341.139BibliographyOldenburg, D. W. and Li, Y. (2005). Inversion for Applied Geophysics: ATutorial. In Butler, K., editor, Near-Surface Geophysics, chapter 5, pages89–150. Society of Exploration Geophysicists.Pelton, W. (1977). Interpretation of Induced Polarization and ResistivityData. PhD thesis, University of Utah.Pelton, W., Ward, S. H., Hallof, P., Sill, W., and Nelson, P. (1978). MineralDiscrimination and Removal of Inductive Coupling with MultifrequencyIP. Geophysics, 43(3):588–609.Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.(2007). Numerical Recipes: The art of scientific computing. CambridgeUniversity Press, Cambridge, 3rd edition.Rathor, B. (1978). Transient electromagnetic field of a polarizable half-spacedue to various current pulses. Geophysical Prospecting, 26:337–351.Rekanos, I. (2012). FDTD schemes for wave propagation in Davidson-Coledispersive media using auxiliary differential equations. IEEE Transactionson Antennas and Propogation, 60(3):1467–1478.Rekanos, I. and Papadopoulos, T. (2010). An auxiliary differential equationmethod for FDTD modeling of wave propagation in Cole-Cole dispersivemedia. IEEE Transactions on Antennas and Propogation, 58(11):3666–3674.Revil, A., Karaoulis, M., Johnson, T., and Kemna, A. (2012). Review:Some low-frequency electrical methods for subsurface characterization andmonitoring in hydrogeology. Hydrogeology Journal, 20(4):617–658.Routh, P. S. and Oldenburg, D. W. (2001). Electromagnetic coupling infrequency-domain induced polarization data: a method for removal. Geo-physical Journal International, 145(1):59–76.Schmutz, M., Ghorbani, A., Vaudelet, P., and Blondel, A. (2014). Cable ar-rangement to reduce electromagnetic coupling effects in spectral-inducedpolarization studies. Geophysics, 79(2):A1–A5.140BibliographyScho¨n, J. (2011). Physical Properties of Rocks: A Workbook.Schwarzbach, C. and Haber, E. (2013). Finite element based inversion fortime-harmonic electromagnetic problems. Geophysical Journal Interna-tional, 193(2):615–634.Seigel, H. O. (1959). Mathematical Formulation and Type Curves for In-duced Polarization. Geophysics, 24(3).Seigel, H. O. (1974). The magnetic induced polarizaton (MIP) method.Geophysics, 39(3).Seigel, H. O., Nabighian, M., Parasnis, D. S., and Vozoff, K. (2007). Theearly history of the induced polarization method. The Leading Edge,26(3):312–321.Simpson, F. and Bahr, K. (2005). Practical magnetotellurics. CambridgeUniversity Press.Slater, L. D. and Glaser, D. R. (2003). Controls on induced polarization insandy unconsolidated sediments and application to aquifer characteriza-tion. Geophysics, 68(5):1547–1558.Smith, J. T. (1996). Conservative modeling of 3-D electromagnetic fields,Part I: Properties and error analysis. Geophysics, 61(5):1308–1318.Smith, R. S., Cheng, L., and Chouteau, M. (2008). Using reversed polarityairborne transient electromagnetic data to map tailings around mine sites.The Leading Edge, (November):1470–1478.Smith, R. S. and Klein, J. (1996). A special circumstance of airborneinduced-polarization measurements. Geophysics, 61(1):66–73.Smith, R. S., Walker, P., Polzer, B., and West, G. F. (1988). The time-domain electromagnetic response of polarizable bodies: an approximateconvolution algorithm. Geophysical Prospecting, 36(April):772–785.141BibliographySnyder, D. D., Merkel, R. H., and Williams, J. T. (1977). Complex formationresistivity - the forgotten half of the resistivity log. In SPWLA 18thAnnual Logging Symposium Transactions, pages 1–39, Houston, Texas.Strang, G. and Fix, G. (2008). An Analysis of the Finite Element Method.Wellesley-Cambridge, Wellesley, MA, 2nd editio edition.Streich, R. (2009). 3D finite-difference frequency-domain modeling ofcontrolled-source electromagnetic data: Direct solution and optimizationfor high accuracy. Geophysics, 74(5):F95–F105.Tarasov, a. and Titov, K. (2013). On the use of the Cole-Cole equations inspectral induced polarization. Geophysical Journal International, pages352–356.Tavlove, A. (1995). Computational Electrodynamics: The Finite-DifferenceTime-Domain Method. Artech House, Norwood, MA.Teixeira, F. (2008). Time-domain finite-difference and finite-element meth-ods for Maxwell equations in complex media. IEEE Transactions on An-tennas and Propogation, 56(8):2150–2166.Titov, K., Kemna, A., Tarasov, A., and Vereecken, H. (2004). Induced Po-larization of Unsaturated Sands Determined through Time Domain Mea-surements. Vadose Zone Journal, 3:1160–1168.Veeken, P. C. H., Kudryavceva, E. O., Putikov, O. F., Legeydo, P. Y.,and Ivanov, S. A. (2012). Modelling induced polarization effects due topyrite in geochemical alteration zones above hydrocarbon accumulations.Petroleum Geoscience, 18(1):59–72.Veeken, P. C. H., Legeydo, P. J., Davidenko, Y. A., Kudryavceva, E. O.,Ivanov, S. A., and Chuvaev, A. (2009a). Benefits of the induced po-larization geoelectric method to hydrocarbon exploration. Geophysics,74(2):B47–B59.142BibliographyVeeken, P. C. H., Legeydo, P. J., and Pesterev, I. (2009b). Geoelectric mod-elling with separation between electromagnetic and induced polarizationfield components. First Break, 27(December):25–36.Vinegar, H. J. and Waxman, M. H. (1984). Induced polarization of shalysands. Geophysics, 49(8):1267–1287.Wait, J. and Debroux, P. (1984). Induced Polarization in ElectromagneticInductive Schemes. Geophysical prospecting, 32:1147–1154.Wang, T. and Hohmann, G. W. (1993). A finite-difference, time-domainsolution for three-dimensional electromagnetic modeling. Geophysics,58(6):797–809.Ward, S. H. and Hohmann, G. W. (1988). Electromagnetic Theory for Geo-physical Applications. In Nabighian, M., editor, Electromagnetic Methodsin Applied Geophysics, volume 1, chapter 4, pages 130–331. Society ofExploration Geophysicists, Tulsa, Oklahoma.Weedon, W. H. and Rappaport, C. M. (1997). A General Method for FDTDModeling of Wave Propagation in Arbitrary Frequency-Dispersive Media.IEEE Transactions on Antennas and Propogation, 45(3):401–410.Weidelt, P. (1982). Response characteristics of coincident loop transientelectromagnetic systems. Geophysics, 47(9):1325–1330.Weller, A., Frangos, W., and Seichter, M. (2000). Three-dimensional in-version of induced polarization data from simulated waste. Journal ofApplied Geophysics, 44:67–83.White, R., Collins, S., and Loke, M. (2003). Resistivity and IP arrays, opti-mised for data collection and inversion. Exploration Geophysics, 34(4):229.Wong, J. (1979). An electrochemical model of the induced-polarization phe-nomenon in disseminated sulfide ores. Geophysics, 44(7):1245–1265.143Yee, K. (1966). Numerical solution of initial boundary value problems in-volving maxwell’s equations in isotropic media. IEEE Transactions onAntennas and Propogation, 14(3):302–307.Yuval and Oldenburg, D. W. (1996). DC resistivity and IP methods inacid mine drainage problems: results from the Copper Cliff mine tailingsimpoundments. Journal of Applied Geophysics, 34(3):187–198.Zaslavsky, M. and Druskin, V. (2010). Solution of time-convolutionaryMaxwells equations using parameter-dependent Krylov subspace reduc-tion. Journal of Computational Physics, 229(12):4831–4839.Zaslavsky, M., Druskin, V., and Knizhnerman, L. (2011). Solution of 3Dtime-domain electromagnetic problems using optimal subspace projection.Geophysics, 76(6):F339–F351.Zhdanov, M. (2008). Generalized effective-medium theory of induced polar-ization. GEOPHYSICS, 73(5):F197–F211.Zonge, K. and Hughes, L. (1991). Controlled source audio-frequency mag-netotellurics. In Nabighian, M. N., editor, Electromagnetic Methods inApplied Geophysics, pages 713–809. Society of Exploration Geophysicists.144Appendix AFinite Volume Modelling ofElectromagnetic MethodsA.1 IntroductionA large portion of this thesis involves the simulation of the response ofthe subsurface to electromagnetic survey techniques. The new simulationtechniques developed in the thesis were built using many aspects of exist-ing numerical methodologies. This appendix provides the reader with anoverview of the numerical techniques that were utilized in the thesis.The overview will limit itself to orthogonal, tensor meshes of the typeused to run all of the three dimensional examples throughout the thesis.Cylindrical meshes, used for the one dimensional examples, will be discussedseparately in Appendix B.I first provide an overview of the discrete approximations used whenmodelling, followed by the derivation of the discrete forms of Maxwell’sequations in both the frequency domain and time domain. The majorityof the content presented is adapted from course notes and lecture materialprovided by Professor Eldad Haber, however several additional valuable ref-erences include Yee (1966), Madden and Mackie (1989), Wang and Hohmann(1993), Druskin and Knizhnerman (1994), Newman and Alumbaugh (1995),Tavlove (1995), Smith (1996), Hyman and Shashkov (1997), Hyman andShashkov (1999), Haber et al. (2000), Haber and Ascher (2001), Haber et al.(2004), Strang and Fix (2008) and Schwarzbach and Haber (2013).145A.2. Discretization in spaceFigure A.1: An example of a non-uniform orthogonal mesh.A.2 Discretization in spaceThis section will derive the discrete spacial approximations required for thesimulation of electromagnetic geophysical experiments in both the time andfrequency domain.A.2.1 Fields and physical propertiesAs mentioned previously, this appendix will be limited to considering non-uniform orthogonal meshes. An example of such a mesh is shown on FigureA.1. While other meshes are potentially better suited to large geophysicalproblems, as they are able to resolve fine detail with a smaller number ofcells, (for example, ocTree meshes used in Haber and Heldmann (2007)),orthogonal meshes are by comparison very simple to work with.Vector fields and physical properties are discretized on the mesh by way146A.2. Discretization in spaceyˆxˆzˆEzExEyByBxBzFigure A.2: Discretization of a cell. The electric field, ~E, is placed on celledges and the magnetic flux density, ~B, is placed on cell faces. Physicalproperties (such as σ and µ) are placed at cell centres.of Yee’s method (Yee, 1966). Physical properties are defined at cell centres,and are assumed to be constant throughout each individual cell. Fields arediscretized on a staggered grid (Yee, 1966) with ~B located on cell boundariesand ~E located at cell edges. This results in different components of the fieldsbeing located at different locations on the cells. The x component of a fielddiscretized onto cell edges is located at the centre of edges parallel to thex axis. The x component of a field discretized onto cell faces is located atthe centre of faces whose normal is parallel to the x axis. The locationsof the different components of the magnetic and electric fields on a cell arepictured in Figure A.2.Indexing cells, edges and facesWhen storing either a discrete cell-centred physical property, or a discreteedge or face based vector field, it is convenient to collapse the three dimen-sional information into a one dimensional vector. The ordering scheme usedis shown in Figure A.3.The index of cells (Figure A.3a) increases first in the x, then the y, and147A.2. Discretization in spacefinally the z direction. Cell number one is located at the bottom, south-westcorner of the mesh. Cell edges (Figure A.3b) and faces (Figure A.3c) arecounted similarly, with the x directed edges and faces counted first, followedby y directed and the z directed edges and faces.1 2 336910 11 1212151819 20 2121242719 20 2122 23 2425 26 27(a) Cell numbering1 2 313 14 1525 26 2737 38 3940 41 4243 44 4546 47 4852566064687276808485 86 87 8889 90 91 9293 94 95 9697 98 99 100104108112113 114 115 116120124128129 130 131 132136140144(b) Edge numbering481216202428323637 38 3949 50 5161 62 63100 101 102103 104 105106 107 108(c) Face numberingFigure A.3: Indexing convention of the mesh (a) cells, (b) edges and (c)faces. In (b) and (c) x directed edges and faces are numbered in red, ydirected edges and faces numbered in blue, and z directed edges and facesnumbered in grey.148A.2. Discretization in spaceA.2.2 The curl operatorThe curl of a vector field ~F is defined by(~∇× ~F)· nˆ ≡ limA→0(1|A|∮C~F · d~r)(A.1)where nˆ is a unit vector orthogonal to the surface A, C is the boundary ofA, and |A| is the area of A. The definition of the curl is approximated by(~∇× ~F ) · nˆ ≈ 1|A|∮C~F · d~r (A.2)Since the equations considered require only the curl of a field defined on celledges, a natural choice for the surface A in equation A.2 is the cell faces.Applying the midpoint rule to approximate the integral, the x componentof the curl of the vector field ~F (depicted in Figure A.4) is given by(~∇× ~F )x ≈(F+z − F−z )∆z − (F+y − F−y )∆y∆y∆z(A.3)Repeating this discretization for each cell face in the mesh produces a sparsematrix C with dimensions consisting of the number of faces by the number ofcells. The matrix C can then be used to compute the discrete approximationof the curl of any edge based vector field. The sparsity structure of C forthe simple three by three mesh used in Figure A.3 is shown in Figure A.5.A.2.3 Discretization of inner productsThe inner product of two vector fields is given by the following integral(~u,~v) =∫Ω~u · ~v dV (A.4)The subsequent discretization of Maxwell’s equations will require the evalu-ation of the inner product for both face and edge based vector fields, as wellas the case where one of the vector fields is multiplied by a scalar physicalproperty distribution.149A.2. Discretization in space∆y∆zF−yF+yF−zF+z(~∇× ~F )xFigure A.4: The quantities required to calculate the discrete approximationof the curl of the edge based vector field ~F at one cell face.Face variablesConsider the inner product of two face variables A and B over a single cellvolume. Splitting the integral into three components(~A, ~B)=∫ΩAxBx dV +∫ΩAyBy dV +∫ΩAzBz dV (A.5)and applying the trapezoid rule to approximate the volume integrals, gives(~A, ~B)≈ 12v(A−xB−x +A+xB+x ) +12v(A−y B−y +A+y B+y )+12v(A−z B−z +A+z B+z)(A.6)150A.2. Discretization in space20 40 60 80 100 120 140Edge Index20406080100Face IndexFigure A.5: Sparsity structure of the discrete curl operator for the 3× 3× 3mesh shown in Figure A.3.151A.2. Discretization in spacewhere v is the volume of the cell, and the values of A and B correspond tothose shown in Figure A.6a.Multiplying ~A by a physical property γ, which is constant over the cell,results in(γ ~A, ~B)≈ 12γv(A−xB−x +A+xB+x ) +12γv(A−y B−y +A+y B+y )+12γv(A−z B−z +A+z B+z ) (A.7)A+y , B+yA−y , B−yA+x , B+xA−x , B−xA+z , B+zA−z , B−z(a) Face variablesC+−x , D+−xC−−x , D−−xC++x , D++xC−+x , D−+xC−−y , D−−y C+−y , D+−yC++y , D++yC−+y , D−+yC+−z , D+−zC−−z , D−−zC−+z , D−+z C++z , D++z(b) Edge variablesFigure A.6: Location of field components involved in discretizing an innerproduct over a single cell.The volume of integration is then expanded to the entire mesh by sum-ming the approximation of the integral over each cell, resulting in the fol-lowing matrix-vector expressions(~A, ~B)≈ v>Acf (A ◦B) (A.8a)(γ ~A, ~B)≈ (v ◦ γ)>Acf (A ◦B) (A.8b)The vectors A and B contain the discrete form of the vector fields ~A and152A.2. Discretization in space~B, v is a vector containing the volume of each cell in the mesh, and γ isa vector containing the value of γ in each cell. The symbol ◦ denotes theHadamard product, or element-wise multiplication. The averaging matrixAcf is defined asAcf =[Acfx Acfy Acfz](A.9)where Acfx, Acfy and Acfz average from cell faces to cell centres in x,y, and zrespectively. Acf is a sparse matrix with dimensions consisting of the numberof cells by number of faces, where every nonzero element in Acf is equal to12 . The sparsity pattern of Acf is shown in Figure A.7a.Equations A.8a and A.8b are then rearranged to the following convenientform, (~A, ~B)≈ B>MfA (A.10a)(γ ~A, ~B)≈ B>MfγA (A.10b)where, Mf and Mfγ are mass matrices defined asMf = diag(Acf>v)(A.11a)Mfγ = diag(Acf>(v ◦ γ))(A.11b)Edge variablesThe discrete inner product of two edge variables is obtained in a similarfashion to that of the face variables. Beginning with a single cell and apply-ing the trapezoid rule to approximate the volume integral for the cell shown153A.2. Discretization in space20 40 60 80 100Face Index1020Cell Index(a)20 40 60 80 100 120 140Edge Index1020Cell Index(b)Figure A.7: Sparsity patterns of the (a) Acf and (b) Ace averaging matrices.The value of every non-zero in Acf is equal to 1/2, while the value of everynon-zero in Ace is 1/4154A.2. Discretization in spacein Figure A.6b) gives(~C, ~D)≈ 14v(C−−x D−−x + C−+x D−+x + C+−x D+−x + C++x D++x)+14v(C−−y D−−y + C−+y D−+y + C+−y D+−y + C++y D++y)+14v(C−−z D−−z + C−+z D−+z + C+−z D+−z + C++z D++z)(A.12)When a cell centered uniform physical property is assigned to the cell, theinner product is given by(γ ~C, ~D)≈ 14γv(C−−x D−−x + C−+x D−+x + C+−x D+−x + C++x D++x)+14γv(C−−y D−−y + C−+y D−+y + C+−y D+−y + C++y D++y)+14γv(C−−z D−−z + C−+z D−+z + C+−z D+−z + C++z D++z)(A.13)Once again we sum over all cells in the mesh to approximate the innerproduct over the entire domain(~C, ~D)≈ v>Ace (C ◦D) (A.14a)(γ ~C, ~D)≈ (v ◦ γ)>Ace (C ◦D) (A.14b)Here, the averaging matrix Ace is given byAce =[Acex Acey Acez](A.15)where Acex, Acey and Acez average one component from the edge location tothe cell centre. Every nonzero element in Ace is equal to14 , and its sparsitypattern is shown in Figure A.7b.Finally, Equations A.14a and A.14b are written in the convenient form(~C, ~B)≈ B>MeC (A.16a)(γ ~C, ~B)≈ B>MeγC (A.16b)155A.3. Maxwell’s equations in the frequency domainwhere Me and Meγ are mass matrices defined byMe = diag(Ace>v)(A.17a)Meγ = diag(Ace>(v ◦ γ))(A.17b)Given the derivations of the curl operator, and edge and face inner prod-ucts, the discretization of Maxwell’s equations are presented in the followingsections.A.3 Maxwell’s equations in the frequency domainModelling of electromagnetic surveys can be carried out either in the fre-quency domain or the time domain. This section will cover the frequencydomain formulation and its discretization. The method described in thissection was used to carry out the modelling presented in Chapters 3 and 2of the thesis.A.3.1 Governing equationsI will take the convention of the Fourier transform to be eiωt. Assuming thisconvention, the inverse transform is given byf(t) =12pi∞∫−∞F (ω)eiωtdω (A.18)The magnetoquasistatic form of Maxwell’s equations in the frequencydomain are~∇× ~E + iω ~B = 0 (A.19a)~∇× ~H = ~J + ~Js (A.19b)∇ · ~B = 0 (A.19c)where ~E is the electric field, ~H is the magnetic field, ~J is the current density,~Js is the current density of the source, and ~B is the magnetic flux density.156A.3. Maxwell’s equations in the frequency domainFor linear, non-dispersive media, the fields and fluxes are related throughthe following constitutive equations~B = µ ~H (A.20a)~J = σ ~E (A.20b)where µ is the magnetic permeability and σ is the electrical conductivity.Computing a particular solution to a system of differential equationsrequires a set of boundary conditions. When considering only the controlledsource experiments of the type discussed in this thesis, one simple choiceof boundary conditions is to have fields go to zero on the boundary of thedomain. This is a reasonable assumption as long as the boundary of thedomain is defined to be far away form the source. There are many ways toexpress this condition, but one that will prove useful later on is1µ~B × nˆ = 0∣∣∣∣δΩ(A.21)where nˆ is normal to the boundary of the domain Ω.A.3.2 Weak formulationPrior to discretization, we will rewrite the system of equations in so-calledweak form (Strang and Fix, 2008). This will provide some advantages overdiscretizing A.19 directly, in particular, it will remove the need to evaluatethe curl in the functional space of ~H.Let ~F and ~W be arbitrary, smooth test functions where ~F is defined inthe same space as ~B, and ~W is defined in the same space as ~E. MultiplyingA.19a by ~F and A.19b (written in terms of ~B and ~E) by ~W and integratingover the domain gives (~∇× ~E, ~F)+ iω(~B, ~F)= 0 (A.22a)(~∇× 1µ~B, ~W)−(σ ~E, ~W)=(~Js, ~W)(A.22b)157A.3. Maxwell’s equations in the frequency domainIntegrating the first term of A.22b by parts and applying Green’s The-orem yields∫Ω~∇× 1µ~B · ~W dV =∫Ω1µ~B · ~∇× ~W dV −∫∂Ω(1µ~B × nˆ)· ~W dS(A.23)Making use of the boundary condition defined in Equation A.21, thesystem simplifyies to (~∇× ~E, ~F)+ iω(~B, ~F)= 0 (A.24a)(1µ~B, ~∇× ~W)−(σ ~E, ~W)=(~Js, ~W)(A.24b)Note that the equations in A.22 require the curl of ~E (a field defined oncell edges) and ~B (a field defined on cell faces), whereas the equations inA.24 require the curl of ~E and ~W , which are both discretized on cell edges.Thus the need to evaluate the curl in the functional space of ~B has beeneliminated.A.3.3 The discrete systemApplying the discrete approximation of the curl operator, and the discreteevaluation of the inner product (Equations A.10 and A.16) Equations A.24aand A.24b are then written in their discrete form as follows,F>MfCE + iωF>MfB =0 (A.25a)W>C>Mfµ−1B−W>MeσE =W>MeJs (A.25b)As these expressions are true for any choice of ~F and ~W , F>Mf canbe dropped from Equation A.25a and W> dropped from Equation A.25b,leavingCE + iωB =0 (A.26a)C>Mfµ−1B−MeσE =MeJs (A.26b)158A.4. Maxwell’s equations in the time DomainThroughout this thesis, I have chosen to work with a B-formulation of thisdiscrete system. That is, Equation A.26b is used to eliminate E from Equa-tion A.26a, leaving a linear system that can be solved to obtain a solutionfor B. Solving Equation A.26b for E givesE = Meσ−1C>Mfµ−1B−Meσ−1MeJs (A.27)and substituting this result into Equation A.26a leaves(Mfµ−1CMeσ−1C>Mfµ−1 + iωMfµ−1)B = Mfµ−1CMeσ−1MeJs (A.28)Prior to the substitution, Equation A.26a was multiplied by Mfµ−1 to makethe matrix multiplying B symmetric. There are approaches other than theB-formulation that can be taken to solve A.26. Rather than eliminatingE, B could be eliminated instead (the approach taken in Schwarzbach andHaber (2013)), or the system could be redefined in terms of potentials (asin Haber et al. (2000)).A.4 Maxwell’s equations in the time DomainThis section presents the discretization of Maxwell’s equations in the timedomain. The methodology forms the basis of the work presented in Chapters4, 5 and 6 of the thesis.A.4.1 Governing equationsIn the time domain, the magnetoquasistatic form of Maxwell’s equations aregiven by~∇× ~e+ ∂~b∂t= 0 (A.29a)~∇× ~h = ~j (A.29b)∇ · ~b = 0 (A.29c)159A.4. Maxwell’s equations in the time DomainFor linear, non-dispersive physical properties, these are accompanied by theconstitutive relationships~b = µ~h (A.30a)~j = σ~e (A.30b)with boundary and initial conditions1µ~b× nˆ = 0∣∣∣∣δΩ(A.31a)~b(0, ~x) = ~b0 (A.31b)~e(0, ~x) = ~e0 (A.31c)When considering inductive sources starting in a steady state, ~e0 = 0.A.4.2 Discretization in timeIn the presence of a conductive media, Maxwell’s equations are stiff. Thatis, numerical methods for solving the differential equations become unstableunless very small time steps are taken. In order to ensure the stability of thesolution, the backward Euler method is applied to discretize in time (Ascherand Petzold, 1998). With this approach, the derivative of the function fevaluated at time t(n+1) is approximated as∂f∂t∣∣∣∣t=t(n+1)≈ f(n+1) − f (n)δt(n+1)(A.32)where f (n+1) denotes the function f evaluated at time t(n+1) andδt(n+1) = t(n+1) − t(n) (A.33)Applying this approach to Equations A.29 (in terms of ~b and ~e) and A.31160A.4. Maxwell’s equations in the time Domainresults in the semi-discretized system of partial differential equations~∇× ~e(n+1) +~b(n+1) −~b(n)δt(n+t)= 0 (A.34a)~∇× 1µ~b(n+1) − σ~e(n+1) = ~j(n+1)s (A.34b)1µ~b(n) × nˆ = 0∣∣∣∣δΩ(A.34c)~b(0) = ~b0 (A.34d)~e(0) = ~e0 (A.34e)We can now follow the approach taken in Section A.3 to discretize the systemA.34 in space.A.4.3 Discretization in spaceThe discretization of A.34 in space follows an identical procedure to thatapplied to the frequency domain case. Multiplying A.34a by the test function~f and A.34b by the test function ~w and integrating over all space gives(~∇× ~e(n+1), ~f)+1δt(n+1)(~b(n+1), ~f)=1δt(n+1)(~b(n), ~f)(A.35a)(~∇× 1µ~b(n+1), ~w)−(σ~e(n+1), ~w)=(~j(n+1)s , ~w)(A.35b)or, by making use of the boundary condition in Equation A.34c,(~∇× ~e(n+1), ~f)+1δt(n+1)(~b(n+1), ~f)=1δt(n+1)(~b(n), ~f)(A.36a)(1µ~b(n+1), ~∇× ~w)−(σ~e(n+1), ~w)=(~j(n+1)s , ~w)(A.36b)Given the discrete curl operator, and the previously derived approximationof the inner products, this system can be written in discrete form asf>MfCe(n+1) +1δt(n+1)f>Mfb(n+1) =1δt(n+1)f>Mfb(n) (A.37a)w>C>Mfµ−1b(n+1) −w>Meσe(n+1) = w>Mej(n+1)s (A.37b)161A.4. Maxwell’s equations in the time DomainSince this formulation is true for any choice of ~f and ~w, choose ~f = 1 = ~wto obtain,Ce(n+1) +1δt(n+1)b(n+1) =1δt(n+1)b(n) (A.38a)C>Mfµ−1b(n+1) −Meσe(n+1) = Mej(n+1)s (A.38b)As with the frequency domain case, a B-formulation of the system ischosen. Solving Equation A.38b for e(n+1)e(n+1) = Meσ−1C>Mfµ−1b(n+1) −Meσ−1Mej(n+1)s (A.39)and using this result to eliminate e(n+1) from Equation A.38a results in thefollowing linear system(Mfµ−1CMeσ−1C>Mfµ−1 +1δt(n+1)Mfµ−1)b(n+1)=1δt(n+1)Mfµ−1b(n) + Mfµ−1CMeσ−1Mej(n+1)s (A.40)Given the source current distribution and the magnetic field at the previoustime step, the linear system can then be solved to recover the magnetic fieldat the next time step.162Appendix BCylindrical Meshes for 1DProblemsB.1 IntroductionFor problems exhibiting symmetry, different meshing schemes can be em-ployed to decrease the computational difficulty. Specifically, if we limit our-selves to problems where the conductivity varies only as a function of depthand whose source term, ~js, is cylindrically symmetric and is exclusively inthe tangential direction, a cylindrical mesh significantly decreases the sizeof the problem. For this class of problems, the electric fields and currentdensities point exclusively in the tangential direction, and the tangentialcomponent of the magnetic fields are equal to zero (Ward and Hohmann,1988).Figure B.1: An example of a cylindrical mesh.163B.2. Discretization of fields and physical propertiesB.2 Discretization of fields and physicalpropertiesThe spatial domain is discretized into a series of rings surrounding a stackof cylinders (Figure B.1). Each ring (Figure B.3a) or cylinder (Figure B.3b)is an individual cell. Cells, edges and faces are counted first in the radialdirection, then in the vertical direction beginning with the bottom inner-most ring.Physical property values are located at cell centres and are assumed tobe constant throughout the cell. The electric fields and current densitiesare located on the edges, all of which point in the tangential direction (Fig-ure B.2a). Magnetic fields are located at the centre of cell faces (FigureB.2b), pointing either in the radial or the vertical direction. In this way,the dimension of the forward problem is reduced from three (x, y and z)components to either two components if solving for magnetic fields (r, z)or one component (θ) if solving for electric fields. This allows for the samevolume of earth to be modelled with a significant reduction in the size ofthe linear system.(a) (b)Figure B.2: (a) Cell edge locations. All cell edges point in the tangentialdirection. (b) Radial (green) and vertical (red) cell face locations.The final discrete expressions (Equation A.28 for the frequency domainor Equation A.40 for the time domain) derived in Appendix A are indepen-dent of the spacial discretization employed. All that is needed to apply theresults of Appendix A to a cylindrical mesh is a new definition of the discretecurl operator and the averaging matrices involved in the mass matrices.164B.3. The curl operatorB.3 The curl operatorTo discretize the curl operator two cases need to be considered; the ringcells that occur away from the mesh’s centre, and the cylindrical cells in themiddle. Using the discrete approximation of the curl (Equation A.2), andthe definitions shown in Figure B.3a, the radial and vertical components ofthe curl of an edge based vector field ~f on a ring-like cell are given by(~∇× ~f)· rˆ ≈ f+− − f++2pir+∆z(B.1a)(~∇× ~f)· zˆ ≈ 2 (r+f++ − r−f−+)r2+ − r2−(B.1b)Applying the definitions shown in Figure B.3b, the components of the curlfor the central cylindrical cells are given by(~∇× ~f)· rˆ ≈ f− − f+2pir∆z(B.2a)(~∇× ~f)· zˆ ≈ 2f+r(B.2b)f+−f++f−+r−r+∆z(~∇× f) · zˆ(~∇× f) · rˆ(a)rf+f−(~∇× f) · zˆ(~∇× f) · rˆ∆z(b)Figure B.3: Elements involved in the approximation of the curl of an edgevariable evaluated on the faces of a (a) ring cell or (b) cylinder cell.165B.4. Inner productsB.4 Inner productsInner products need to be considered for both face and edge variables, andeach of the two types of cells need to be considered separately.B.4.1 Face variablesFace variables have non-zero components in the radial and the vertical di-rections so the inner product of vector fields ~A and ~B are thus given by(~A, ~B)=∫ΩArBrdV +∫ΩAzBzdV (B.3)rA+z , B+zA−z , B−zA+r , B+r∆z(a)r−r+∆zA+z , B+zA−z , B−zA−r , B−r A+r , B+r(b) Location of face elements for ring cellsFigure B.4: (a) Location of face elements for cylindrical cells (b) Locationof face elements for ring cellsCylindrical cellsI will start by approximating the integral in Equation B.3 over a singlecylindrical cell (Figure B.4a). With this volume of integration, EquationA.6a becomes(~A, ~B)=∆z∫02pi∫0r∫0ArBrrdrdθdz +∆z∫02pi∫0r∫0AzBzrdrdθdz (B.4)166B.4. Inner productsAs only one radially directed face surrounds the cylindrical cell, the quan-tity ArBr is taken to be constant throughout the cell. The quantity AzBzis assumed to vary linearly in the z-direction. Using these assumptions,Equation B.4 is approximated by(~A, ~B)≈ vA+r B+r +12v(A−z B−z +A+z B+z)(B.5)where v is the volume of the cell, v = pir2∆z.Ring cellsTaking the volume of integration to be a single ring like cell (Figure B.4b)Equation A.6a can be written as(~A, ~B)=∆z∫02pi∫0r+∫r−ArBrrdrdθdz +∆z∫02pi∫0r+∫r−AzBzrdrdθdz (B.6)Assuming ArBr and AzBz vary linearly in the radial and vertical directionsrespectively, Equation B.6 is then approximated by(~A, ~B)≈13v(2r− + r+r− + r+A−r B−r +r− + 2r+r− + r+A+r B+r)+12v(A−z B−z +A+z B+z)(B.7)where v is the volume of the cell, v = pi(r2+ − r2−)∆z.Extending to the entire meshEquation B.3 is now extended to the entire mesh by summing the approxi-mated integral over each cell, resulting in the following relation(~A, ~B)≈ v>Acf (A ◦B) (B.8)where v is a vector containing the volume of each cell, A and B are vectorscontaining the discrete values of ~A and ~B, and Acf is a sparse averaging167B.4. Inner productsmatrix that averages from cell faces to cell centers. The values of the non-zero entries of Acf are obtained from Equations B.5 and B.7. Equation B.8has the same form as that obtained when using a 3D tensor mesh, with theonly difference being the definition of the averaging matrix, Acf .B.4.2 Edge variablesEdge variables have non-zero components only in the tangential directionsso the inner product of vector fields ~C and ~D can be written as(~C, ~D)=∫ΩCθDθdV (B.9)rC+θ , D+θC−θ , D−θ∆z(a)C+−θ , D+−θC++θ , D++θC−+θ , D−+θC−−θ , D−−θr−r+∆z(b) Location of edge elements for ringcellsFigure B.5: (a) Location of face elements for cylindrical cells (b) Locationof edge elements for ring cellsCylindrical cellsIntegrating over the single cylindrical cell pictured in Figure B.5a, the inte-gral in Equation B.9 becomes(~C, ~D)=∆z∫02pi∫0r∫0CθDθrdrdθdz (B.10)168B.4. Inner productsOnly two edges border cylindrical cells. Assuming Cθ and Dθ vary linearlyin the z direction, and because cylindrical cells have only two borderingedges, the integral is approximated by(~C, ~D)≈ 12v(C−θ D−θ + C+θ D+θ)(B.11)where v = pir2∆z is the volume of the cell.Ring cellsFor a single ring cell (Figure B.5b), the integral in Equation B.9 is given by(~C, ~D)≈∆z∫02pi∫0r+∫r−CθDθrdrdθdz (B.12)Assuming Cθ and Dθ vary linearly throughout the ring, and becauserings cells have four bordering edges, the integral is approximated by(~C, ~D)≈ 16v(2r− + r+r− + r+(C−−θ D−−θ + C−+θ D−+θ)+r− + 2r+r− + r+(C+−θ D+−θ + C++θ D++θ))(B.13)where v = pi(r2+ − r2−)∆z is the volume of the cell. .Extending to the entire meshEquation B.3 is now extended to the entire mesh by summing the approxi-mated integral over each cell. This results in(~C, ~D)≈ v>Ace(C ◦D) (B.14)The averaging matrix Ace averages from cell edges to cell centres, and thevalues of its non-zero entries are obtained from Equations B.11 and B.13.This expression is again the same as that obtained for the tensor mesh, witha different definition of Ace.169Appendix CA Reduced Model forGrounded SourceExperimentsIP phenomena arising from traditional galvanic survey techniques are rarelymodelled by the full Maxwell’s equations.Modelling is generally performed in the frequency domain, using thestatic Maxwell’s equation while assuming a frequency dependent conductiv-ity. The common frequency dependent model is~J = σ(ω)~∇φ (C.1a)∇ · ~J = ~Js (C.1b)A similar simplified model can be used in time in order to generate IP decaycurves.C.1 A reduced model for low frequenciesIn this appendix I will use the auxiliary differential equation approach forDebye media discussed in Chapter 5 to represent the chargeable response.However, the convolution or Pade´ methods from Chapters 4 and 6 couldalso be used.The magnetoquasistatic form of Maxwell’s equations in the time domain170C.1. A reduced model for low frequenciesare given by~∇× ~e+ ∂~b∂t= 0 (C.2a)~∇× ~h = ~j (C.2b)~j + τ(1− η)∂~j∂t= σ∞(1− η)~e+ σ∞τ(1− η)∂~e∂t(C.2c)Note that this system has two time dependent equations; first, Faraday’slaw (C.2a) contains the derivative ∂~b∂t , and second, Ohm’s law (C.2c) containstime derivatives of both the electric fields and the current densities.Assuming a finite time source and, that after some period of time thetime derivative of the magnetic fields approximately equals zero, while thederivatives of the electric fields and the current densities do not, then~∇× ~e ≈ 0 (C.3)Thus, the electric field can be written as the gradient of a scalar potential,φ~e = ~∇φ (C.4)By taking the divergence of Ampere’s law to eliminate the ~∇ × ~h, asimplified time dependent system is obtained that can be solved for theelectric potential φ and the flux density ~j:∇ · ~j = ∇ · ~js (C.5a)~j + τ(1− η)∂~j∂t= σ∞(1− η)~e+ σ∞τ(1− η)∂~e∂t(C.5b)This system is the time-dependent equivalent of the usual simplifiedfrequency modelling of IP phenomena when inductive electromagnetic effectsare neglected. It assumes that the IP effects are still present at times when∂~b∂t is small enough to be neglected.171C.2. DiscretizationClearly, this assumption breaks down in the presence of large conductors.The contamination of induced polarization surveys by significant inductionis often referred to as electromagnetic coupling. Therefore, using the reducedsystem as a modelling and inversion tool should be done carefully.C.2 DiscretizationUsing the methodology discussed in Appendix A we can discretize the re-duced model. Setting e = Gφ where G is a discretization of the nodalgradient (see Haber and Ascher (2001) for details of such discretization)and following the same procedure used previously results in the discretesystemG>j(n+1) = G>j(n+1)s (C.6a)Mej(n+1) = MeAGφ(n+1) − j(n)p (C.6b)In these equations,j(n)p = MeEGφ(n) −MeJ j(n) (C.7)and the mass matrices MeA MeE and MeJ are identical to those defined inSection 5.3.Eliminating j(n+1) from these equations results in a linear systemG>Me−1MeAGφ(n+1) = G>Me−1j(n)p + G>j(n+1)s (C.8)that can be solved for φ(n+1). Once this has been solved, the current den-sities at the next time step are calculated usingj(n+1) = Me−1MeAGφ(n+1) −Me−1j(n)p (C.9)172C.3. Synthetic example: Effects of EM-coupling on induced polarizationdataC.3 Synthetic example: Effects of EM-couplingon induced polarization dataNeglecting the effects of induction is often reasonable in induced polarizationsurveys, but this assumption can break down in more conductive environ-ments (Dey and Morrison, 1973; Hallof, 1974). Attempts have been made tomitigate these issues by estimating and removing the inductive effects fromcollected data (Fullagar et al., 2000; Routh and Oldenburg, 2001; Veekenet al., 2009b), or designing arrays to minimize their effects (White et al.,2003). Nevertheless the effects of EM coupling still pose a problem for manysurveys.In this appendix, I have developed a procedure for calculating the IPresponses when the the effects of induction are neglected. Comparing theseresults to the modelled response obtained while including induction allowsfor the evaluation of the effects of EM coupling on IP data.As an example, I simulate the responses of a typical gradient array in-duced polarization survey. A single transmitter, oriented in the East-Westdirection with an A-B offset of 1225m is simulated. The transmitter wireis run just south of the area of interest. This survey geometry is shown inFigure C.1.In the first example, a single chargeable block (dark gray square in FigureC.1) is buried in a uniform background. The block exhibits Debye dispersionwith σ∞ = 0.1 S/m, τ = 1s, and η = 0.1 and the background has a conduc-tivity of 0.005 S/m. The top of the block is located 100m below the surfaceand extends to a depth of 200m. Three simulations were run; the electro-magnetic response of the model with no chargeability present, the reducedmodel of the IP response in the absence of electromagnetic induction, andthe total response including both IP effects and electromagnetic induction.Figure C.2a shows the x-component of the electric fields from the three sim-ulations observed at the point (0, 0). Due to the very resistive backgroundmodel, induction effects decay rapidly, and the reduced model accuratelypredicts the response from approximately 2× 10−2 seconds onwards.In a second example, 50m of conductive overburden is placed at the173C.3. Synthetic example: Effects of EM-coupling on induced polarizationdata−500 0 500Easting (m)−5000500Northing (m)Figure C.1: Layout of synthetic gradient array experiments. Thick lineshows the path of the transmitter wire, and dots show receiver locations.Dark gray box shows the extents of the chargeable target and the light greyboxes show the location of additional conductive blocks present in the secondexample.174C.3. Synthetic example: Effects of EM-coupling on induced polarizationdatasurface of the model shown in Figure C.1. Three additional conductiveblocks are also placed directly under the overburden, extending to a depthof 150m below the surface. The overburden and conductive blocks have aconductivity of 1 S/m. The chargeable block and the background have thesame properties as in the first example.In each of the examples, the conductivity was discretized onto a recti-linear mesh with 55× 55× 54 cells in the x, y and z directions, respectively.This includes a uniform core region of 25× 25× 24 cells that are 25m on aside as well as 10 padding cells in each direction. The padding cells expandaway from the core region with an expansion rate of 1.3. This results ina core region that is 625m × 625m × 600m with approximately 5.4km ofpadding in every direction. Time was discretized into 4 intervals of constantδt, (1×10−4s, 5×10−4s, 2.5×10−3s, and 6.25×10−3s), with 65 steps takenfor each δt.The x-component of the electric field at (0,0) is shown in Figure C.2b.This time, the inductive effects take much longer to decay due to the presenceof the conductive units. The reduced model still accurately predicts the latetime response of the system, but in this case, inductive effects dominate untilapproximately 2×10−1s in this location. Figure C.3 shows the time evolutionof the x-component of the electric fields during the simulation. Immediatelyfollowing the switch off of the transmitter current (t = 1.0× 10−3s) electricfields are induced directly below the path of the transmitter wire. Thesefields completely mask the presence of the chargeable target and diffusethrough the area, interacting with the other conductive units. It is not untillate times (t = 2.5 × 10−1s) that the response of the chargeable block iseasily identified. At the end of the simulation (t = 0.1s) the electric fieldsresemble those of a dipole, centred at the location of the chargeable block,and they agree with the response of the reduced model.175C.3. Synthetic example: Effects of EM-coupling on induced polarizationdata10-310-210-1100Time (s)10-810-710-610-510-4exTotal responseEM responseIP response(a)10-310-210-1100Time (s)10-1610-1210-810-4exTotal responseEM responseIP response(b)Figure C.2: The x-component of the electric fields recorded at (0, 0) for (a)a chargeable block in a resistive half-space and (b) a chargeable block in aresistive half-space with a conductive overburden.176C.3. Synthetic example: Effects of EM-coupling on induced polarizationdata−400 0 400Easting (m)−4000400Northing (m)ex-1.4e-04-1.1e-04-7.2e-05-3.6e-050.0e+003.6e-057.2e-051.1e-041.4e-04(a) t = 1.0× 10−3s−400 0 400Easting (m)−4000400Northing (m)ex-2.4e-05-1.8e-05-1.2e-05-6.0e-060.0e+006.0e-061.2e-051.8e-052.4e-05(b) t = 5.0× 10−3s−400 0 400Easting (m)−4000400Northing (m)ex-4.2e-06-3.1e-06-2.1e-06-1.0e-060.0e+001.1e-062.1e-063.2e-064.2e-06(c) t = 1.7× 10−2s−400 0 400Easting (m)−4000400Northing (m)ex-1.9e-07-1.4e-07-9.6e-08-4.8e-080.0e+004.8e-089.6e-081.4e-071.9e-07(d) t = 6.4× 10−2s−400 0 400Easting (m)−4000400Northing (m)ex-4.8e-10-3.6e-10-2.4e-10-1.2e-100.0e+001.2e-102.4e-103.6e-104.8e-10(e) t = 2.5× 10−1s−400 0 400Easting (m)−4000400Northing (m)ex-2.2e-10-1.7e-10-1.1e-10-5.6e-110.0e+005.6e-111.1e-101.7e-102.2e-10(f) t = 1× 100sFigure C.3: Time evolution of the x-component of the electric fields recordedabove the model shown in Figure C.1.177
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Induced polarization effects in inductive source electromagnetic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Induced polarization effects in inductive source electromagnetic data Marchant, David 2015
pdf
Page Metadata
Item Metadata
Title | Induced polarization effects in inductive source electromagnetic data |
Creator |
Marchant, David |
Publisher | University of British Columbia |
Date Issued | 2015 |
Description | Induced polarization (IP) surveys are commonly conducted to map the distribution of electrical chargeability, a diagnostic physical property in mineral exploration and in many environmental problems. Although these surveys have been successful in the past, the galvanic sources required make their application labour intensive and prevents them from being applied in some settings. The ability to detect chargeability with a geophysical technique that employs inductive sources, eliminating the need to inject current into the ground, would provide a valuable tool to applied geophysicists. In this work, two aspects of inductive source induced polarization are examined. First, a new methodology for processing inductive source frequency domain EM data to identify IP effects is presented. The method makes use of the asymptotic behaviour of the secondary magnetic fields at low frequency. A new quantity, referred to as the ISIP datum, is defined so that it equals zero at low frequencies for any frequency-independent (non-chargeable) conductivity distribution. Thus, any non-zero response in the ISIP data indicates the presence of chargeable material. Numerical simulations demonstrate that the method can be applied even in complicated geological situations. A 3D inversion algorithm is developed to recover the chargeability from the ISIP data and the inversion is demonstrated on synthetic examples. Understanding the impacts of IP effects on time-domain electromagnetic data requires the ability to simulate common survey techniques while taking chargeability into account. Most existing techniques preform this modelling in the frequency domain prior to transforming their results to the time domain. Application of those techniques on three dimensional problems can be computationally limiting. In the second part of this thesis, three new techniques for forward modelling the time-domain electromagnetic response of chargeable materials directly in the time domain are developed. The first considers the convolution in the time domain directly, while the others use auxiliary differential equations to analytically transform the governing equations into the time domain. The resulting methods are verified by comparing their results with analytic solutions. The potential application of the method was demonstrated by modelling the occurrence of negative transients in airborne time-domain electromagnetic data. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-10-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
IsShownAt | 10.14288/1.0135704 |
URI | http://hdl.handle.net/2429/52685 |
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 | 2015-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_may_marchant_david.pdf [ 15.41MB ]
- Metadata
- JSON: 24-1.0135704.json
- JSON-LD: 24-1.0135704-ld.json
- RDF/XML (Pretty): 24-1.0135704-rdf.xml
- RDF/JSON: 24-1.0135704-rdf.json
- Turtle: 24-1.0135704-turtle.txt
- N-Triples: 24-1.0135704-rdf-ntriples.txt
- Original Record: 24-1.0135704-source.json
- Full Text
- 24-1.0135704-fulltext.txt
- Citation
- 24-1.0135704.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:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0135704/manifest