Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

On recovering distributed induced polarization information from time-domain electromagnetic data Kang, Seogi 2018

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

Item Metadata

Download

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

Full Text

On recovering distributed induced polarizationinformation from time-domain electromagnetic databySeogi KangGeosystem Engineering, Hanyang University, 2010Geophysics, Hanyang University, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POST DOCTORALSTUDIES(Geophysics)The University of British Columbia(Vancouver)March 2018c Seogi Kang, 2018AbstractThe electrical conductivity of earth materials is frequency-dependent. This is dueto a phenomenon known as induced polarization (IP), wherein electrical chargesbuild-up under the application of an electric field. Macroscopically, rocks may beconsidered chargeable, as they act like electric capacitors. The goal of this the-sis is to show how IP data can be extracted from time-domain electromagnetic(TEM) data, then inverted to recover information about chargeable targets. Al-though both frequency and time-domain electromagnetic (EM) surveys measureIP signals, this dissertation will focus solely on TEM. To recover chargeability in-formation, the following TEM-IP inversion workflow is developed. (1) Extract abackground conductivity model that is assumed free of IP signals. (2) Decouplethe TEM and IP signals by subtracting the fundamental responses estimated usingthe background conductivity. (3) Invert the resultant IP data to recover pseudo-chargeabilities at multiple times for a set of 3D volumes. This is used to inferthe location and dimensions of chargeable targets. (4) Carry out further analysesof pseudo-chargeabilities at multiple times to estimate intrinsic parameters such asCole-Cole chargeability and its associated time constant. For grounded sources, theworkflow is implemented for a synthetic gradient array example. Results show thatthe early time signals, which are often discarded, can be used to estimate the back-ground conductivity. Applying the workflow to inductive sources such as airborneEM (AEM) is more challenging, as steady-state electric fields are not produced.This was overcome by developing an IP function which (1) accurately character-izes how electric fields from inductive sources behave in the earth and (2) allowsthe recovery of a 3D chargeability by solving a linear inverse problem. The efficacyof the aforementioned approach is validated using field AEM surveys over the Mt.iiMilligan porphyry deposit in British Columbia and Tli Kwi Cho kimberlite depositin the Northwestern Territories. For the kimberlite deposit, the recovered charge-ability information is able to distinguish two distinct kimberlite units. To validatethe approach, a 3D rock model for Tli Kiw Cho is constructed using the recov-ered chargeability and background conductivity. This model is compared againstgeological models obtained through drilling and shows good agreement.iiiLay SummaryWhen materials are excited by electromagnetic (EM) sources, they may exhibitsignificant capapcitive properties due to the buildup of electrical charges; a phe-nomenon known as induced polarization (IP). Rocks exhibiting IP are said to bechargeable. Signals measured during EM geophysical surveys can include IP ef-fects from chargeable rocks. This has been observed in mining, oil and gas, andgroundwater applications. The strength of the IP effects, and their impact on field-collected data, depends upon the concentration and distribution of chageable mate-rials below the surface. In this thesis, I develop a workflow capable of recoveringa 3D chargeability distribution from EM geophysical data. The workflow is vali-dated using both synthetic and field examples. In particular, I focus on the Tli KwiCho kimberlite complex in NWT, Canada.ivPrefaceThe following document presents original research I completed at the Geophysi-cal Inversion Facility (GIF) in the Department of Earth, Ocean and AtmosphericSciences at the University of British Columbia (UBC), Vancouver, Canada. Sig-nificant portions of this research was used to produce three peer-reviewed journalarticles and five expanded conference proceedings. Work as part of this dissertationhas also been presented at seven conferences.Chapters 2 and 3 contain the text and figures from the published paper by Kangand Oldenburg [2016] in Geophysical Journal International. The development ofa TEM-IP inversion workflow which linearizes inductive source IP (ISIP) and ac-counts for inductive responses was originally proposed by Dr. Douglas Olden-burg. Under the supervision of Dr. Douglas Oldenburg, I produced all subsequentderivations, numerical modeling codes, numerical tests and manuscript prepara-tions. Numerical modeling codes were completed under the supervision of Dr.Elded Haber. These codes were written within the framework of an open-sourcegeophysical simulation and inversion package known as SIMPEG [Cockett et al.,2015]. Successful integration of original coding into this framework was done withthe assistance of Mr. Rowan Cockett and Ms. Lindsey Heagy. Original numericalmodeling codes developed for this dissertation were validated against the EMTDIPcode developed by Dr. David Marchant [Marchant et al., 2014].The workflow I developed from Dr. Douglas Oldenburg’s original concept isapplied to airborne IP data in Chapter 4. This chapter is divided into three sections.Section 4.1 summarizes content presented in the SEG expanded abstract [Kang andOldenburg, 2015]. Section 4.2summarizes content presented in the SEG expandedabstract [Kang et al., 2014]. Here, the conductivity inversion is based upon PhDvwork by Dr. Dikun Yang [Yang et al., 2014]. Additional numerical experimentswere done using pre-existing SimPEG modeling packages. Sections 4.1 and 4.2were edited with the help of Dr. Douglas Oldenburg. Original work summarized inSections 4.1 and 4.2 are currently in preparation to be submitted as a journal arti-cle. Section 4.3 contains text and figures from the published paper by Kang et al.[2017] in Interpretation. This work was done in collaboration with several mem-bers of the UBC Geophysical Inversion Facility (UBC-GIF) group, most notably,Mr. Dominique Fournier. Mr. Dominique Fournier made significant contributionsto the conductivity inversion code which was used; see Fournier et al. [2017] inInterpretation. All other numerical experiments are original work. Manuscriptwith editing was done with the help of Mr. Dominique Fournier and Dr. DouglasOldenburg.Chapter 5 is a revised version of the journal paper by Kang and Oldenburg[2017], which was published in Geophysical Prospecting. The idea was initiatedfrom Dr. Douglas Oldenburg’s suggestion to apply the workflow to DC-IP data. Iperformed numerical experiments with SIMPEG EM and DC-IP packages, and Iwrote the journal paper with editing help from Dr. Douglas Oldenburg.Appendices A, B, C contain materials relevant to Chapters 2 and 3.This thesis manuscript is entirely original work. The manuscript was com-pleted with editing help from Dr. Douglas Oldenburg, Dr. Eldad Haber and Dr.Randolph Enkin.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Polarization mechanisms . . . . . . . . . . . . . . . . . . . . . . 81.3 Complex Resistivity and Conductivity . . . . . . . . . . . . . . . 151.4 Lab-scale IP Measurements . . . . . . . . . . . . . . . . . . . . . 161.5 Modelling Maxwell’s Equations . . . . . . . . . . . . . . . . . . 201.6 TEM-IP Inversion Workflow . . . . . . . . . . . . . . . . . . . . 321.7 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37vii2.1 Pseudo-chargeability . . . . . . . . . . . . . . . . . . . . . . . . 382.2 Linear IP Function . . . . . . . . . . . . . . . . . . . . . . . . . 422.3 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . 482.3.1 IP Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 482.3.2 Polarization Currents . . . . . . . . . . . . . . . . . . . . 532.3.3 IP Currents . . . . . . . . . . . . . . . . . . . . . . . . . 572.3.4 Validations of Linearization . . . . . . . . . . . . . . . . 632.4 Effective Pseudo-chargeability . . . . . . . . . . . . . . . . . . . 652.5 Positivity on Time History of Electric Field . . . . . . . . . . . . 672.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 773 3D IP Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 793.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 793.2 3D IP Inversion with a Linear Equation . . . . . . . . . . . . . . 793.3 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . 813.3.1 Incorrect Conductivity . . . . . . . . . . . . . . . . . . . 823.3.2 Extracting Intrinsic IP Parameters . . . . . . . . . . . . . 883.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 894 Airborne Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 924.1 Synthetic Example . . . . . . . . . . . . . . . . . . . . . . . . . 994.1.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 994.1.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1014.1.3 Conductivity Inversion . . . . . . . . . . . . . . . . . . . 1054.1.4 EM-decoupling . . . . . . . . . . . . . . . . . . . . . . . 1074.1.5 IP Inversion . . . . . . . . . . . . . . . . . . . . . . . . . 1124.1.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 1164.2 Field Example: Mt Milligan . . . . . . . . . . . . . . . . . . . . 1174.2.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1174.2.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1214.2.3 Conductivity Inversion . . . . . . . . . . . . . . . . . . . 1244.2.4 EM-decoupling . . . . . . . . . . . . . . . . . . . . . . . 1254.2.5 IP Inversion . . . . . . . . . . . . . . . . . . . . . . . . . 128viii4.2.6 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . 1314.2.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 1354.3 Field Example: DO-27/18 kimberlites . . . . . . . . . . . . . . . 1364.3.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1364.3.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1394.3.3 Conductivity Inversion . . . . . . . . . . . . . . . . . . . 1434.3.4 EM-decoupling . . . . . . . . . . . . . . . . . . . . . . . 1464.3.5 IP Inversion . . . . . . . . . . . . . . . . . . . . . . . . . 1494.3.6 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . 1554.3.7 Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . 1584.3.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 1605 Grounded Source Example . . . . . . . . . . . . . . . . . . . . . . . 1625.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1625.2 TEM-IP Inversion Workflow for Grounded Source . . . . . . . . . 1665.3 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1685.4 Conductivity Inversion . . . . . . . . . . . . . . . . . . . . . . . 1725.5 EM-decoupling . . . . . . . . . . . . . . . . . . . . . . . . . . . 1745.6 IP Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1765.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1785.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1806 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1816.1 Development of the TEM-IP Inversion Workflow . . . . . . . . . 1826.2 Application to Airborne IP Data . . . . . . . . . . . . . . . . . . 1836.3 Application to DC-IP Data . . . . . . . . . . . . . . . . . . . . . 1856.4 Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . 1866.5 Concluding comments . . . . . . . . . . . . . . . . . . . . . . . 188Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191A Handling Multiple Sources for AEM Surveys . . . . . . . . . . . . . 199B Extracting Intrinsic IP Parameters . . . . . . . . . . . . . . . . . . . 205ixC Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207C.1 Steady-state Maxwell’s Equations . . . . . . . . . . . . . . . . . 207C.2 Linearized Kernel for IP Responses . . . . . . . . . . . . . . . . 208xList of TablesTable 1.1 Cole-Cole parameters from different rocks (obtained from Pel-ton et al. [1978]). . . . . . . . . . . . . . . . . . . . . . . . . 7Table 2.1 Amplitudes of decomposed IP currents at two marked points(white stars) shown in Fig. 2.10(b). Units in A/m2 . . . . . . . 58Table 4.1 Cole-Cole parameters of four anomalous bodies (A1-A4). . . . 100Table 4.2 Parameters of conductivity inversion. . . . . . . . . . . . . . . 105Table 4.3 Parameters of IP inversion. . . . . . . . . . . . . . . . . . . . 113Table 4.4 Recovered Cole-Cole parameters of four chargeable blocks (A1-A4). Values in parenthesis indicates true values. . . . . . . . . 114Table 4.5 Expected physical property contrast for kimberlite deposits inthe Lac de Gras region [Power and Hildes, 2007]. . . . . . . . 138Table 4.6 Mean and standard deviation of estimated t (µs) and h for A1-A4. Used values of c for fitting are given. . . . . . . . . . . . 153Table 4.7 Petrophysical domains built from inversions of airborne EMdata sets. Here s , h˜E , and h˜L correspondingly conductivity,and pseudo-chargeability at 130 and 410µs. . . . . . . . . . . 156Table 5.1 Conductivity (s•) and resistivity (r•) at infinite frequency, andCole-Cole chargeability (h) values for five units: A1-A4 andhalf-space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169Table 5.2 Parameters for the TEM and DC inversions. See Section 3.2for explanation of parameters. shal f indicates conductivity ofthe homogeneous half-space, which has a value of 0.01 S/m. . . 172xiList of FiguresFigure 1.1 Cole-Cole resistivity as a function of frequency [Pelton et al.,1978]. Black and red lines indicate real and imaginary partof the Cole-Cole resistivity. Used Cole-Cole parameters are:r0=1.25 Wm, h=0.2, t=0.1 s, and c=0.5. . . . . . . . . . . . . 6Figure 1.2 Conceptual diagram of (a) half-cycle input current waveformand (b) measured IP responses at surface electrodes. . . . . . 6Figure 1.3 Polarization phenomenon at a pore throat: (a) Initial state, (b)While electric field is applied, and (c) Right after removal ofthe electric field. (d) Input current. (e) Overvoltage at on-time.(f) Overvoltage at off-time. . . . . . . . . . . . . . . . . . . . 9Figure 1.4 Microscope images of rocks. (a) Illite (a clay mineral) withsurface area-to-volume ratio of 100 m2/g(1000 times greaterthan sandstone). (b) Quartz overgrowths in sandstone with sur-face area-to-volume ratio of 0.1 m2/g. Modified from [GPG,2018]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10Figure 1.5 Electrode polarization: (a) Initial state, (b) While electric fieldis applied, and (c) Right after removal of the electric field.Modified from Revil et al. [2015]. . . . . . . . . . . . . . . . 12xiiFigure 1.6 Electrical double layer at the surface of insulating minerals.X indicates active sites at mineral surface, which attracts metal-lic cations M+. M+ is attached to the active sites, and boundedby the Outer Helmholtz Plane (OHP). The OHP separates thestern layer from the diffuse layer. A- indicates anion. Obtainedfrom Revil [2013] . . . . . . . . . . . . . . . . . . . . . . . . 13Figure 1.7 Polarization mechanisms related to Electrical Double Layer(EDL): (a) Initial state, (b) While electric field is applied, and(c) Right after removal of the electric field. Modified from Re-vil [2013]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13Figure 1.8 Maxwell-Wagner polarization: (a) Initial state. For a rock con-taining three different phases: solid, water, and air that havedifferent conductivity and dielectric permittivity. (b) Polariza-tion charge can be built when an electric field is applied. (c)Charges diffuse away when the electric field is removed, gen-erating opposite direction of the current flow to the appliedelectric field. Modified from Chen and Or [2006]. . . . . . . . 14Figure 1.9 Lab scale IP measurement system in Geological Survey of Canada(GSC). The photo is taken from Cowan [2015]. . . . . . . . . 18Figure 1.10 Complex impedance measurements using the Geological Sur-vey of Canada system [Enkin et al.]. . . . . . . . . . . . . . . 19Figure 1.11 Cole-Cole response in frequency domain (a) and time (b) do-main. The Cole-Cole parameters are s• = 102 S/m, h = 0.5,t = 0.1 s, and c=1. The arrow indicates a Dirac-delta function(s•d (t)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 1.12 3D chargeability model. (a) Plan and (b) Section views. Whitesolid line on (a) shows the wire path, and the arrow indicatesthe direction of current in the wire when current is turned on.THe half-space conductivity is 0.05 S/m. . . . . . . . . . . . 25xiiiFigure 1.13 Time curves of the normalized potential difference at a receiverlocation (marked as black solid line with dots in Fig. 1.12(a)).Black, blue and red lines correspondingly indicate observed,fundamental and IP data. Solid and dashed lines distinguishpositive and negative values. Black vertical dashed lines atthree time channels indicate times at which the electric fielddistributions in Fig. 1.14 are provided. The halfspace conduc-tivity is 0.05 S/m. . . . . . . . . . . . . . . . . . . . . . . . 26Figure 1.14 3D electric field distributions at three different times in the on-time. Top, middle and bottom panel correspondingly show theelectric fields at 6 ms, 77 ms, and 1717 ms. (a) Plan view at thesurface and (b) Section view at Northing 0 m. The halfspaceconductivity is 0.05 S/m. . . . . . . . . . . . . . . . . . . . . 27Figure 1.15 Time decaying curves of the normalized potential differencewhen the current is turned off. Black, blue, and red lines re-spectively indicate observed, fundamental and IP data. Dottedand solid lines distinguish positive and negative values. Halfs-pace conductivity is 0.05 S/m. . . . . . . . . . . . . . . . . . 28Figure 1.16 3D electric field distribution at three different times when thecurrent is turned off. Top, middle and bottom panel corre-spondingly show the electric fields at 6 ms, 715 ms, and 2116ms. (a) Plan view at the surface and (b) Section view at Nor-thing 0 m. The halfspace conductivity is 0.05 S/m. . . . . . . 29Figure 1.17 Observed potential difference at three different off times. (a) 6ms, (b) 77 ms, and (c) 1717 ms. Halfspace conductivity is 0.05S/m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 1.18 Time curves of the normalized potential difference at a receiverlocation (marked as black solid line with dots in Fig. 1.12(a)).Halfspace conductivity is 0.5 S/m. . . . . . . . . . . . . . . . 30Figure 1.19 Time decaying curves of the normalized potential differencewhen current is turned off. Halfspace conductivity is 0.5 S/m. 31xivFigure 1.20 Observed potential difference at three different off times. (a) 6ms, (b) 77 ms, and (c) 1717 ms. Halfspace conductivity is 0.5S/m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Figure 2.1 Cole-Cole response in frequency domain (a) and time (b) do-main. The Cole-Cole parameters are s• = 102 S/m, h = 0.5,t = 0.01 s, and c=1. The arrow shown in Fig. 2.1b indicates adelta function (s•d (t)). . . . . . . . . . . . . . . . . . . . . 42Figure 2.2 Convolution of 4s and ~e resulting in polarization current,~jpol . (a) DC source. (b) Inductive source. 4s and h˜ I arerespectively defined in eqs. ( 2.8) and ( 2.7). . . . . . . . . . . 43Figure 2.3 Conceptual diagram for the amplitude of the fundamental elec-tric fields. (a) DC source. (b) Inductive source. . . . . . . . . 43Figure 2.4 Plan (a) and section (b) views of the IP model. The solid linein (a) delineates the boundary of the IP body. Solid circles in(a) denote the sounding locations. In (b) the conductivity s2 isvariable so that canonical, conductive and resistive blocks canbe examined . . . . . . . . . . . . . . . . . . . . . . . . . . 50Figure 2.5 Time decaying curves of the observations (d; black line), fun-damental (dF ; blue line) and IP (dIP; red line) responses. Allthree cases: (a) canonical, (b) conductive and (c) resistive arepresented. Right and left panels show bz and  ∂bz∂ t . The verti-cal black dotted line indicates the time at which the polariza-tion field reaches its maximum value. The flight height of thecollocated transmitting and receiving loop is 30 m above thesurface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 2.6 Interpolated maps of observed (left panel), fundamental (mid-dle panel) and IP (right panel) responses. Two time channelsat (a) 0.86 ms and (b) 6.7 ms are presented. The white linecontours a zero-crossing in the observed response. . . . . . . 52xvFigure 2.7 Maps of reference currents: (a) canonical and (b) conductivemodels. Left and right panels show plan and section views at -125 m-depth and 0 m-easting, respectively. A source is locatedat (-200 m, 0 m, 30 m). Black arrows and colored backgroundrespectively indicate the direction and amplitude of the current.Black solid lines outline the boundary of the chargeable body. 54Figure 2.8 Maps of polarization currents: (a) canonical and (b) conductivemodels at 0.86 ms. Left and right panels show plan and sectionviews at -125 m-depth and 0 m-easting, respectively. A sourceis located at (-200 m, 0 m, 30 m). Black arrows and shadedvalues respectively indicate the direction and amplitude of thecurrent. Black solid lines outlin boundary of the surface or thechargeable body. . . . . . . . . . . . . . . . . . . . . . . . . 55Figure 2.9 Maps of polarization currents: (a) canonical and (b) conduc-tive models at 6.7 ms. Left and right panels show plan andsection views at -125 m-depth and 0 m-easting, respectively.A source is located at (-200 m, 0 m, 30 m). Black arrows andshaded values indicate the direction and amplitude of the cur-rent, respectively. Black solid outlines boundary of the surfaceor the chargeable body. . . . . . . . . . . . . . . . . . . . . . 56Figure 2.10 Decomposition of the IP currents as~jpol (left panel),s•~—f IP(middle panel), ands•~aIP (right panel) at 0.86 ms. Plan viewmaps of the currents at -125 m-depth are shown: (a) canonical,(b) conductive, and (c) resistive cases. . . . . . . . . . . . . . 59Figure 2.11 Comparisons of contributions of ~jpol , s•~—f IP, and s•~aIPto the observed IP magnetic field. Solid line indicates true bIPzresponses. Stars, rectangles, and circles correspondingly indi-cate each IP response generated by applying the Biot-Savartlaw to ~jpol , s•~—f IP, and s•~aIP. Empty and solid markersrepresent positive and negative values, respectively. . . . . . . 60Figure 2.12 Interpolated maps of (a) true and (b) approximate IP currentsat 0.86 ms. Left and right columns respectively show plan andsection view maps at -125 m-depth and 0 m-easting. . . . . . 61xviFigure 2.13 Interpolated maps of (a) true and (b) approximate IP currentsat 6.7 ms. Left and right columns respectively show plan andsection view maps at -125 m-depth and 0 m-easting. . . . . . 62Figure 2.14 Comparison of true and approximate IP responses (bIPz ). Black,blue, and red colors respectively indicate canonical, conduc-tive, and resistive cases. Solid lines indicate true bIPz , com-puted by subtracting the fundamental response from the obser-vation. The stars are the application of Biot-Savart to true IPcurrent and generate bIPz BS. Empty circles show our approxi-mate bIPz approx response. . . . . . . . . . . . . . . . . . . . . 64Figure 2.15 Conceptual diagram of EM induction process from a loop trans-mitter when a vertical conductor is embedded in a homoge-neous half-space. After the transmitter current is turned-offfrom the loop, currents are induced in the subsurface and dif-fuse away like a smoke ring. This smoke ring currents will bechanneled into the vertical conductor (dashed lines), and therewill be vortex currents rotating (red dashed circle) in the con-ductor due to the time-varying magnetic field. In the conduc-tor, at the top voxel (yellow) the half-space channeled currentand the vortex current are in the same direction (in to the page),whereas at the bottom voxel they are in the opposite direction. 70Figure 2.16 Plan and section views of a 3D conductivity model is shown inright and left panels, respectively. A vertical conductor (0.005S/m) is embedded in a homogeneous half-space (0.001 S/m).Green solid circle indicates the horizontal location of the trans-mitter loop having 13 m-radius, and located 30 m above thesurface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71xviiFigure 2.17 The current densities at three different times due to a 13 m-radius loop located 30 m above the surface (marked as greensolid circle) with a conductivity model having a vertical con-ductor shown in Fig. 2.16. Plan and section views of thecurrent density at three different times after the current is off:(a) 0.01 ms, (b) 0.07 ms, and (c) 1.17 ms. The time history ofthe electric field at the top and bottom voxels in the conductor(yellow and green crosses, respectively) are presented in Fig.2.19. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 2.18 The peak current (or reference current) selected from the timehistory of the current density shown in Fig. 2.17. Plan andsections views are shown in right and left panels, respectively.The time history of the electric field at the top and bottom vox-els in the conductor (yellow and green crosses, respectively)are presented in Fig. 2.19. . . . . . . . . . . . . . . . . . . . 73Figure 2.19 The time history of the electric field, wre f , at the two voxelsshown in Figs. 2.17 and 2.18. Yellow and green distinguishwre f for the top and bottom voxels, respectively. (a) wre f at theboth voxels in semi-log scale (log(t)wre f ). (b) wre f at the topvoxel only showing later times after 0.1 ms with linear scale intime. The negative values of wre f exist at the very late-time forthe bottom voxel (green) as shown in (b). . . . . . . . . . . . 74Figure 2.20 (a) Intrinsic pseudo-chargeability, h˜ I as a function of time.Used Cole-Cole parameters are denoted in (a). (b) Pseudo-chargeability, h˜ , calculated by convolving h˜ I and wre f at thebottom voxel shown in Fig. 2.19. Solid line and cross marksdistinguish h˜ computed with and without positive projectionon wre f , P0[wre f ]. . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 2.21 Computed IP data, dIP, using linear IP function. IP data arecomputed with and without positive projection on wre f anddenoted as dIP and dIP0 , respectively; red dashed line and redcrosses are distinguishing dIP and dIP0 , respectively. . . . . . . 76xviiiFigure 3.1 Effect of depth weighting in 3D IP inversion. (a) True pseudo-chargeability model on vertical section at 0 m-northing. Re-covered pseudo-chargeability models (b) without depth weight-ing and (c) with depth weighting. . . . . . . . . . . . . . . . 83Figure 3.2 IP responses on a profile line at 0 m-northing. IP responses arecomputed from perturbed s• models. Halfspace conductivity(s1) is perturbed: 2⇥s1 and 0.5⇥s1, which respectively re-sulting in overestimated (dotted line) and underestimated (dashedline) IP responses. The solid line shows the true IP response. 85Figure 3.3 Recovered pseudo-chargeability sections from 3D IP inversionsat 0 m-northing. (a) dIP with true s1. (b) dIP with 2⇥s1. (c)dIP with 0.5⇥s1. (d) dIP with 0.5⇥s1 and the positivity con-straint on the pseudo-chargeability. White dashed lines contourzero-crossing lines. . . . . . . . . . . . . . . . . . . . . . . . 86Figure 3.4 Recovered pseudo-chargeability sections from the 3D IP in-versions at 0 m-northing. (a) True and (b) incorrect s• is usedto compute sensitivity function. For the incorrect sensitivity Iused a halfspace conductivity s1. . . . . . . . . . . . . . . . . 87Figure 3.5 Section views of recovered: (a) time constant and (b) charge-ability. Any region where the pseudo-chargeability shown inFig. 3.4a is smaller that 0.001 is ignored in this analysis, andblanked. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89Figure 3.6 Comparisons of the observed pseudo-chargeability at a singlepixel in a chargeable body from the 14 inversions at 14 timechannels, and the predicted pseudo-chargeability. The emptycircles and solid line respectively indicate predicted and ob-served pseudo-chargeability. The estimated time constant andchargeability are respectively expressed as test and hest . Thetrue values for t and h are respectively 0.005 s and 0.2. . . . . 90Figure 4.1 Geometry of the VTEM system (Geotech). . . . . . . . . . . 95xixFigure 4.2 Observed VTEM data over the Mt Milligan porphyry deposit,BC, Canada. (a) Profile line at Northing 6109000 m, whichcrosses Harris fault and Rainbow faults. (b) Time decay atNorthing 6109000 m and Easting 433800 m. Here black andred lines distinguish positive and negative values. . . . . . . . 96Figure 4.3 Section of a chargeable cylinder embedded in a halfspace. Con-ductivity of the halfspace is s1=0.001 S/m. The conductivechargeable cylinder is parameterized with a Cole-Cole modeland has higher conductivity than that of halfspace: s•=0.01S/m. Other parameters are given in the diagram. . . . . . . . . 97Figure 4.4 Propagation of electric fields in time: 0.1 ms (left), 0.5 ms(middle), and 7.9 ms (right). ey is plotted in the x-z plane, andfundamental electric field is rotating in counter-clock-wise-direction. In the late-time, the reversed direction of ey is ob-served due to strong polarization effects. Boundaries of theair-earth interface and chargeable cylinder are shown as blacklines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Figure 4.5 Time decays of observed (d), fundamental (dF ), and IP (dIP)voltages over a chargeable cylinder. Solid and dashed linesdistinguish positive and negative values. (a) Only observedvoltage. (b) All three voltages with distinction of early, inter-mediate, and late times. Vertical grey lines indicate three times(0.1, 0.5, and 7.9 ms) when electric fields are plotted in Fig. 4.4. 98Figure 4.6 Sectional views of s•(x,y,z)model. Solid lines delineate bound-aries of four IP blocks. . . . . . . . . . . . . . . . . . . . . . 99Figure 4.7 2D map of the simulated data including both EM and IP effects(d). A1-A4 indicate corresponding anomalies due to four IPbodies. Black dotted lines indicate boundaries of the four IPbodies. Negative responses are shaded as white regions. Blackdots are locations of Fig. 4.9. Black dots are locations of Fig.4.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102xxFigure 4.8 Profile line of the observed data, d at Northing (a) 630 m and(b) -630 m. Black and red lines distinguish positive and nega-tive data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103Figure 4.9 Time decays of the observed data, d at four sounding locationsclose to A1-A4. Here solid and dashed lines distinguish posi-tive and negative data. Error bars are shown from 0.1 to 1 ms,which indicate the uncertainty used in conductivity inversion.Three grey vertical liens indicate time channels shown in Fig.4.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104Figure 4.10 Sectional views of the recovered conductivity model (sest).Solid lines delineate boundaries of the four IP blocks. . . . . . 106Figure 4.11 Comparison of true and estimated fundamental responses at0.2, 1.8 and 5.6 ms. (a) true fundamental response: F [s•] and(b) estimated fundamental response: F [sest ]. . . . . . . . . . 108Figure 4.12 Time decays of the d, F [s•], and F [sest ] at four sounding loca-tions close to A1-A4. Here solid and dashed lines distinguishpositive and negative datum. Error bars are shown from 0.1 to1 ms, which indicate uncertainty used in conductivity inversion. 109Figure 4.13 Comparison of true and estimated IP responses at 0.2, 1.8 and5.6 ms. (a) true IP response: dIP[s•] and (b) estimated IPresponse: dIP[sest ]. . . . . . . . . . . . . . . . . . . . . . . . 110Figure 4.14 Time decays of the d, dIP[s•], and dIP[sest ] at four soundinglocations close to A1-A4. Here solid and dashed lines or solidand empty circles distinguish positive and negative datum. . . 111Figure 4.15 Regional removal procedure. Left, middle and right panelscorrespondingly indicate before removal, estimated regionalfields, and after removal of dIP[sest ] at 1.8 ms. Black dots andempty white circles indicate all stations and selected stationsto estimate regional field. . . . . . . . . . . . . . . . . . . . . 111Figure 4.16 Sectional views of the recovered dh˜/dt models at 1.8 ms.(a) Without regional removal. (b) With regional removal. . . . 114xxiFigure 4.17 Comparison of observed and predicted pseudo-chargeabilityfor four chargeable anomalies: A1-A4. Solid lines and emptycircles distinguish observed and predicted pseudo-chargeability.115Figure 4.18 Location of Mt Milligan and QUEST survey area. . . . . . . . 118Figure 4.19 Geology and VTEM survey at Mt Milligan porphyry depositin British Columbia, Canada [Yang and Oldenburg, 2012]. . . 119Figure 4.20 Normalized current waveform of the VTEM survey (2007) atMt Milligan deposit. . . . . . . . . . . . . . . . . . . . . . . 119Figure 4.21 Topography of the VTEM survey of the Mt Milligan region.Geology indicated with white lines. . . . . . . . . . . . . . . 120Figure 4.22 2D map of the observed VTEM data at eight different timesranging from 4.6-7.1 ms. Solid white lines show boundaries ofdifferent gelogical units and the white areas indicate where theresponse has gone negative. . . . . . . . . . . . . . . . . . . 122Figure 4.23 2D map of the VTEM data at 2.7 ms (left panel) and time de-cays at three stations close the negative anomalies at A1-A3(right panel). Solid white lines on the right panel show bound-aries of different geological units, and white dashed lines indi-cate the approximate locations of where there are. . . . . . . . 123Figure 4.24 3D Conductivity model of Mt Milligan porphyry deposit ob-tained by Yang et al. [2014]. . . . . . . . . . . . . . . . . . . 124Figure 4.25 EM-decoupling of the VTEM data over Mt Milligan region.2Dmaps of d, F [sest ], dIP[sest ] are shown at (a) 2.7, (b) 0.7 and(c) 0.2 ms. F [sest ] is subtracted from d to obtain dIP[sest ]. A4-A5 are IP anomalies recognized by EM-decoupling, whereasA1-A3 were shown as negatives in d. A1-A3 and A4-A5 aremarked as grey and black dashed circles, respectively. . . . . 126Figure 4.26 Time decays of d, F [sest ], dIP[sest ] at the five IP anomalies:A1-A5. Top left panel showed 2D map of dIP[sest ] at 2.7 ms. 127xxiiFigure 4.27 Plan and section views of the recovered pseudo-chargeability:(a) 2.7, (b) 0.7, and (c) 0.2 ms. Five anomalies recognizedfrom dIP[sest ] are marked as dashed circles on plan view (leftpanels). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129Figure 4.28 Comparison of observed and predicted IP data for the IP inver-sion at three time channels: 2.7, 0.7 and 0.2 ms. . . . . . . . . 130Figure 4.29 Plan and section views of 3D conductivity and pseudo-chargeabilityover Mt Milligan area. . . . . . . . . . . . . . . . . . . . . . 132Figure 4.30 The 3D rock model obtained from both conductivity and pseudo-chargeability models. Red and black dots indicates fault struc-tures, and three mineralized zones, respectively. R3 shows therock units obtained from the airborne IP data. . . . . . . . . . 133Figure 4.31 The 3D rock model obtained from both conductivity and pseudo-chargeability models. Red and black dots indicates fault struc-tures, and three mineralized zones, respectively. R2 shows therock units obtained from the airborne IP data. . . . . . . . . . 134Figure 4.32 Location map for the Tli Kwi Cho (TKC) kimberlites, NWT.DO-18 and DO-27 are two main kimberlite pipes at TKC re-gion. PK, VK, and HK correspondingly indicate pyroclas-tic kimberlite (PK), volcaniclastic kimberlite (VK), and hy-pabyssal kimberlite (HK). . . . . . . . . . . . . . . . . . . . 137Figure 4.33 Schematic diagram of a kimberlite pipe in the Lac de Grasregion (Modified from Devriese et al. [2017]). A lake may bepresent after glaciation and is often used as a first indicator ofa possible kimberlite. Transverse lines are from the DIGHEMsurvey (1992). . . . . . . . . . . . . . . . . . . . . . . . . . 138xxiiiFigure 4.34 Plan maps of four EM data sets at TKC: (a) DIGHEM (56000Hz), (b) NanoTEM (77 µs), (c) AeroTEM II (26 µs), (d) VTEM(90 µs). For TEM data sets, a smaller region (red box) closeto DO-18 and -27 is presented. The black line is a contour lineof the negative anomaly (-8 nT/s) from AeroTEM data at 26µs. The white line shows the boundary of the lakes. Nega-tive anomalies: A1-A4 are correspondingly marked as purple,yellow, red, and green solid circles; A1-A3 showed strong neg-atives for all TEM data at an early-time as shown in Fig. 4.35. 140Figure 4.35 Plan maps of four EM data sets at TKC: (a) DIGHEM (7200Hz), (b) NanoTEM (603 µs), (c) AeroTEM II (534 µs), (d)VTEM (680 µs). For TEM data sets, a smaller region (redbox) close to DO-18 and -27 is presented. The black line is acontour line of the negative anomaly (-8 nT/s) from AeroTEMdata at 26 µs. The white line shows the boundary of the lakes.Negative anomalies: A1-A4 are correspondingly marked aspurple, yellow, red, and green solid circles; A1-A3 showedstrong negatives for all TEM data at an early-time as shown inFig. 4.35. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141Figure 4.36 Transient curves of NanoTEM and VTEM data. (a) NanoTEMsoundings away from the pipe and representative of background(blue lines) and over the DO-18 pipe (black lines). They aremarked as blue and black solid circles in Fig. 4.34(b). (b)VTEM soundings at A1-A4 ( correspondingly purple, yellow,red, and green lines). They are marked as purple, yellow, red,and green solid circles in Fig. 4.34 (d) and Fig. 4.35(d). Solidand dashed lines distinguish positive and negative observations. 142Figure 4.37 Plan and section views of the recovered conductivity modelfrom the cooperative inversion of the VTEM and DIGHEMdata sets: (a) VTEM and (b) DIGHEM. The white outlinesdelineate boundaries of the lake. The green outlines show theextent of DO-27 and DO-18 at the surface, based on drilling. . 144xxivFigure 4.38 Observed and estimated fundamental responses at TKC. Foursounding locations at A1-A4 are marked as solid circles. Re-gions having negative values are shown as white. . . . . . . . 145Figure 4.39 Time decaying curves of the d (black line), F [sest ] (blue line),and dIP[sest ] (red line) data at (a) A1, (b) A2, (c) A3, and (d)A4. Solid and dashed lines distinguish positive and negativevalues. Vertical black dashed line indicate 130 and 410 µs,respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 147Figure 4.40 Plan maps of d, F [sest ] and dIP[sest ] at (a) 130 and (b) 410µs. Left, middle, and right panels correspondingly show the d,F [sest ] and dIP[sest ]. . . . . . . . . . . . . . . . . . . . . . . 148Figure 4.41 Plan and section views of the recovered pseudo-chargeabilitymodel from the 3D IP inversion of the raw IP at (a) 130 µsand (b) 410 µs, respectively. Left panel shows plan map at99 m below surface. Top and bottom of right panels show A-A’ and B-B’ sections, respectively. Solid and dashed red linesdelineate contours of the recovered pseudo-chargeability at 50s1 (130 µs) and 10 s1 (410 µs). The black outlines delineateboundaries of the lake. The green outlines show the extent ofDO-27 and DO-18 at the surface, based on drilling. . . . . . . 150Figure 4.42 Cole-Cole conductivity as a function of frequency [Pelton et al.,1978]. Black and red lines indicate real and imaginary part ofthe Cole-Cole conductivity. Used Cole-Cole parameters are:s•=1 S/m, h=0.2, and c=0.5; two t values are considered:103 s and 104 s. Solid and dashed line indicates complexconductivity generated with t = 103 s and t = 104 s, re-spectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152Figure 4.43 Comparison of the observed (lines) and predicted (solid cir-cles) pseudo-chargeability at cells close to A1-A4 (correspond-ingly purple, yellow, red, and green colors). (a) All time decay-ing curves of the observed and predicted pseudo-chargeability.(b) Median values of the observed and predicted pseudo-chargeabilityat each time. . . . . . . . . . . . . . . . . . . . . . . . . . . 153xxvFigure 4.44 A cross plot between the estimated time constant (t) and charge-ability (h) at cells close to A1-A4. Solid circles shaded aspurple, yellow, red, and green colors correspondingly indicateA1-A4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154Figure 4.45 Comparative sections through the TKC kimberlites. (a) Over-laid anomalies of three physical property models: conductiv-ity (dark grey), pseudo-chargeability at 130 (red), and 410 µs(green). (b) Generated rock model from three anomalies. Thewhite outlines delineate boundaries of the lake. The greenoutlines show the extent of DO-27 and DO-18 at the surface,based on drilling. . . . . . . . . . . . . . . . . . . . . . . . . 157Figure 4.46 Comparison of 3D geology and rock model. (a) 3D geologicmodel obtained from the known geology and drilling resultsat the TKC area (b) 3D petrophysical model obtained fromairborne EM. The black outlines delineate boundaries of thelake. The green outlines show the extent of DO-27 and DO-18at the surface, based on drilling. . . . . . . . . . . . . . . . . 159Figure 5.1 Conceptual diagram of a ground-based grounded source withhalf-duty cycle current waveform. . . . . . . . . . . . . . . . 164Figure 5.2 A typical example of overvoltage effects in electric IP data. . . 164Figure 5.3 Observed voltage with EM induction effects. EM effects dom-inate the early on and off-time data. Inset figure shows the en-larged off-time voltage and it demonstrates early EM inductioneffects (negative) and late-time IP effects (positive). Dashedand solid lines distinguish positive and negative voltages. . . . 165Figure 5.4 Plan and section views of the 3D mesh. Black solid lines showthe boundaries of four blocks (A1-A4). Only A2 and A3 (redlabels) are chargeable. Arrows indicate a wire path for thegrounded source. Black dots indicate potential electrodes. . . 169Figure 5.5 Plan maps of the observed DC data: (a) Voltage and (b) Ap-parent resistivity. . . . . . . . . . . . . . . . . . . . . . . . . 170xxviFigure 5.6 Plan maps of the observed TEM data at (a) 5 ms, (b) 80 ms (c)130, and (d) 500 ms. Dashed and solid contours differentiatenegative and positive data. . . . . . . . . . . . . . . . . . . . 170Figure 5.7 Time decaying curves of the observed (d) and fundamental(dF ) data. Four sounding locations close to (a) A1, (b) A2, (c)A3, and (d) A4 are presented. Blue and black color indicatesd and dF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171Figure 5.8 Recovered conductivity models from: (a) TEM and (b) DCinversions by inverting off-time (EM) and on-time data (DC),respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . 173Figure 5.9 Plan maps for 80 ms. (a) observed data, (b) estimated funda-mental using sEMest , and (c) IP. . . . . . . . . . . . . . . . . . 175Figure 5.10 True and estimated fundamental responses at 80 ms. (a) Truefundamental response; (b), (c), and (d) correspondingly indi-cate estimated fundamental responses using sEMest , sDCest , andshal fest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175Figure 5.11 True and estimated IP responses at 80 ms. (a) True dIP; (b),(c), and (d) correspondingly indicate estimated IP responsesusing sEMest , sDCest , and shal fest . . . . . . . . . . . . . . . . . . . 175Figure 5.12 Recovered chargeability models from IP inversions. IP datacomputed using an estimated conductivity model from (a) TEMand (b) DC inversions, respectively. . . . . . . . . . . . . . . 177Figure 5.13 True and estimated IP responses at 130 ms. (a) True dIP; (b),(c), and (d) correspondingly indicate estimated IP responseswith use of EM, DC, and half-space conductivity. . . . . . . . 179Figure 5.14 Decay curves of the observed, fundamental and IP data at A3.(a) True and estimated fundamental responses with the ob-served data. (b) True and IP responses with the observed data. 179Figure A.1 Normalized weights for the conductive case for all source lo-cations. A single pixel located at (0 m, 0 m, -75 m) is used. . 202xxviiFigure A.2 Time decays of we(t) and wˆ(t) for the conductive case. Asingle pixel located at (0 m, 0 m, -75 m) is used. Solid lineand dashed lines correspond to we(t) and wˆk(t) for all sources(k= 1, . . . , nTx); wˆk at the center source located at (0 m, 0 m,30 m) is marked as solid circles. A number of we(t) curves areoverlaid due to the symmetric position of source locations tothe conductive block. . . . . . . . . . . . . . . . . . . . . . 203Figure A.3 Comparison of true and approximate bIPz responses at 0.86 mson a plan view map. . . . . . . . . . . . . . . . . . . . . . . 204xxviiiAcknowledgmentsThis thesis work has built upon a number of people’s previous efforts, and wouldhave not been completed without the guidance and support from a large group ofindividualstod. So, I would to take this as an opportunity to express my gratitudeto them.First of all, I would like to express my deepest gratitude to my supervisor Prof.Doug Oldenburg for his academical, emotional, and financial supports for the pastfive and half years. When I first came to Vancouver, Canada from South Korea, Iwas eager to see and experience how electromagnetic geophysics has been used inthe bigger world, and learn the state-of-art technology. Doug not only taught meabout academic knowledge about the applied geophysics, but also opened up thedoors that I can be connected to people from various disciplines which allowed meto experience how I can help solve real world problems. Further, his moral supportwas a fuel for me whenever I was down.Next, I would like to thank my committee members, Prof. Eldad Haber andDr. Randy Enkin. The foundation in computational electromagnetics that I havelearned from Eldad was the basis for most of computational aspects in the thesis.Randy’s guidance on understanding mechanisms of induced polarization effectsmade significant impact for me to interpret the field examples shown in the thesis.In addition, I would like thank my external examiner Prof. Richard Smith for hiscritical comments suggestions on my thesis.The funding of this work was generously provided from a variety of sources.Thanks to University of British Columbia, the Department of Earth and OceanSciences, the Mathematics of Information Technology and Complex Systems, andthe Canadian Exploration Geophysical Society for their support.xxixMany thanks to all of my friends and colleagues at Geophysical Inversion Fa-cility, Simulation and Parameter Estimation in Geophysics community, and Uni-versity of British Columbia that I have shared time with over the years. Thanks toDavid Marchant, Dikun Yang, Rowan Cockett, Lindsey Heagy, Mike Mitchell, Do-minique Fournier, Sarah Deveriese, Kris Davis, Mike McMillan, Gudni Rosenkjar,Luz Angelica, Thibaut Astic, Devin Cowan, Patrick Belliveau, Daniel Bild-Enkin,Roman Shekhtman, Gesa Myer, Mike Wathen, Julie Nutini, and Jilmarie Stephensfor making my Canadian life delightful.Finally, I would like to show my love and gratitude to Angela, for all of thesupport and love that she provided rain or shine and being my lifetime companion.xxxDedicationTo my parents and in the loving memory of my mother.xxxiChapter 1Introduction1.1 Research MotivationThe electrical conductivity of earth materials is usually frequency dependent withthe effective conductivity decreasing with decreasing frequency. The change inconductivity is due to the buildup of electric charges that occur under the applica-tion of an electric field. The rock is electrically polarized and thus it can be saidthat the rock is chargeable. This charge build-up phenomenon is called induced po-larization (IP). At a microscopic level the causes of the IP depends upon the veryfine scale structure of the rocks.Rocks can have different polarization mechanisms and this will result in dif-ferent IP characteristics as a function of frequency. There are many ways to pa-rameterize the complex behavior but a commonly used one is the Cole-Cole model[Pelton et al., 1978]:r(w) = r0h1h⇣1 11+(ıwt)c⌘i, (1.1)where r0 is resistivity at zero frequency, h is chargeability, t is a time constant,c is a frequency exponent. Fig. 1.1 shows real and imaginary part of complexresistivity as a function of frequency. Using eq. (1.1), Pelton et al. [1978] fittedmeasured complex impedances of different rock samples and provided a range ofCole-Cole parameters. These are summarized in Table 1.1. IP surveys have been1successfully conducted in various geoscience applications. For mining applica-tions, IP surveys are recognized as a principal geophysical technique for findingdisseminated sulfide or porphyry deposits. This is mainly because of the metal-lic polarization effects arising from metallic minerals such as pyrite. On the otherhand, non-metallic polarization effects (e.g. membrane polarization) can also gen-erate IP signals in various geological settings. IP surveys have been used in a broadspectrum of environmental problems including hydraulic conductivity estimation[Ho¨rdt et al., 2007] and hydrocarbon contaminant mapping [Kemna et al., 2012]. Inaddition, Veeken et al. [2009a,b] has used IP surveys to find offshore hydrocarbons.A common type of IP survey is an electrical IP (EIP) survey (or often calledDC-IP survey) where current electrodes, connected to a generator, are used to in-ject currents into the earth (grounded or galvanic source) and voltages as a functionof time are measured at potential electrodes. For time-domain measurements, theinput current waveform has an on- and off-time as shown in Fig. 1.2. Typically,the on-time is of the order of 0.1-10 seconds. At the potential electrodes on the sur-face, an electric potential difference is measured while input currents are switchedon and off. Fig. 1.2 shows measured voltage, and here Vs indicates the secondaryvoltage (off-time),V0 andV• are potential differences at zero and infinite frequency,respectively. Assuming there are no electromagnetic (EM) induction effects, anysignal in the off-time is the result of IP phenomena. Hence, any datum that cap-tures some aspect of this secondary voltage is considered to be an IP datum. Acommonly used form of the time-domain IP data, dIP can be written asdIP =1V0Z t2t1Vs(t) (1.2)Here t1 and t2 are arbitrary times included in the off-time with t1 < t2. This measurecorrespond to to the area of the shaded region in Fig. 1.2, and when t is measuredin milliseconds (ms), then its unit is ms. However, each time channel of secondarypotential can be considered as an IP datum, then the units of dIP can be mV/V. Inthe frequency domain, since measured voltage is complex-valued, both amplitudeand phase are considered as IP data. Assuming that the voltage is measured attwo different frequencies, percentage frequency effect (PFE), and phase difference2(mrad), are also commonly used as IP data. When only one frequency is used,phase is considered as IP data, as when there is no chargeability, the phase is zero.Interpretation of the IP data is most commonly carried out by taking one ofthe many definitions of IP data, assume that the ground chargeability is small,linearize the problem [Seigel, 1959], and invert to find 2D or 3D volumes of achargeability [Oldenburg and Li, 1994]. This IP inversion is commonly carriedout using a two-stage inversion procedure: (a) first recover a resistivity distributionat zero frequency from DC data (in practical terms, the latest on-time data or thelowest frequency), then (b) recover a chargeability distribution by inverting dIPdata using a sensitivity function from the previous DC inversion.Not all chargeable volumes are the economic resources sought and it is desir-able to distinguish between the different chargeable units. Determining different IPcharacteristics by interpreting multiple time and frequency channels has thus beena long term goal. The generic term for this study is spectral IP (SIP). For instance,one may want to differentiate between a porphyry and a magnetite embedded inthe earth. Both of them have significant chargeability, but different time constantsas shown in Table 1.1. By inverting SIP data, and recovering both chargeabilityand time constants, these two different rocks can possibly be distinguished.For extracting spectral IP information from IP data, I consider two main ap-proaches. The first approach is simple and most often used. I invert each timeor frequency IP data separately, and then interpret the recovered models to distin-guish different IP sources using a complex resistivity model [Yuval and Oldenburg,1997, Kemna et al., 2004, Ho¨rdt et al., 2006]. The second is more sophisticated,and this inverts multiple time or frequency IP data simultaneously by parameteriz-ing the complex resistivity model with a few parameters, e.g. Cole-Cole [Kemnaet al., 2004, Fiandaca et al., 2012, 2013]. Either approach can provide distributedIP information that is insightful. A principal challenge for inverting SIP data in3D space is the increased dimensionality due to frequency- or time-dependent re-sistivity. For instance, in the frequency domain, the output of the inversion can ber(x,y,z;w). This is four-dimensional function. Even if I chose a parameterizationsuch as Cole-Cole [Pelton et al., 1978], which reduces the dimensionality fromthe number of frequencies to four representative values, I still have four times thenumber of parameters in the model that must be solved. With the first approach3where each time or frequency is inverted separately, this dimensionality issue maynot be a critical problem. However, for the second approach where all parametersare found at once, one needs an effective strategy to handle the increased dimen-sionality to obtain a successful 3D inversion. Considering this generic challengefor inverting SIP data, I chose the first approach to recover 3D distribution of IPparameters.A significant assumption behind the interpretation of DC-IP data is that EMinduction effects can be ignored. EM induction effects arise from the time varyingmagnetic fields generated by the transmitter. These general electric fields are super-posed on our desired IP signal. Depending on both the resistivity of the earth andmeasured times or frequencies, there can be situations where this assumption is notreasonable. For several decades, this issue has been treated using EM-decouplingtechniques; that is, use a method to remove EM induction effects from the obser-vations and thereby extract IP signals. There have been various EM-decouplingtechniques suggested based on approximate representations of EM induction ef-fects using a halfspace or layered-earth resistivity model [Wynn and Zonge, 1975,Routh and Oldenburg, 2001]. Although they have been effective for some casehistories, I believe that the EM-decoupling issue has to be revisited using the latest3D EM forward modelling and inversion capabilities.On the other hand, whenever an electric field is applied to the earth, IP sig-nals can be generated from a chargeable body [Macnae, 1988]. Thus, not onlygrounded sources, but also inductive sources such as airborne loop sources canexcite chargeable bodies. However, the electric field generated by an inductivesource is solely based on EM induction, and hence it is fundamentally differentfrom the grounded source excitation. This difference needs to be considered inorder to proceed to the interpretations of inductive source IP (ISIP) data. The air-borne IP problem in time-domain is important in geophysics research. Negativetransients measured in various time-domain airborne EM (AEM) systems are com-monly observed [Smith and Klein, 1996, Jansen and Witherly, 2004, Kang et al.,2017, Kaminski and Viezzoli, 2017]. In theory, a chargeable earth is required toexplain such negatives [Weidelt, 1982] particularly from a coincident loop system.I call either AEM survey or data including IP effects as an airborne IP. Invertingthese airborne IP data, and recovering IP parameters, may yield information about4near surface structure that is relevant to mineral exploration and/or environmentalproblems such as permafrost [Smith and Klein, 1996, Kang et al., 2017]. Theseairborne IP data have typically been inverted by invoking a halfspace or layered-earth assumption [Kratzer and Macnae, 2012, Viezzoli et al., 2017]. However, IPsignals are small and 3D analysis is likely necessary for many problems. Carryingout this analysis for airborne IP data is an important topic.In this thesis I consider three open questions in the area of inductive source(particularly airborne IP) and grounded source (DC-IP).1. How can we understand the fundamentals of inductive source IP, and can thetime-domain airborne IP data be linearized similar to DC-IP data?2. Can distributed IP information be recovered from the airborne IP data?3. Can EM effects in DC-IP data be effectively removed to generate high qual-ity IP data?I develop an effective workflow that can separate EM and IP portions in the obser-vation, and recover distributed IP information in 3D space. To apply the workflowto airborne IP data, the first question needs to be answered, and this is treated inChapter 2. Based upon this linearization of IP data, 3D IP inversion code is de-veloped in Chapter 3; this provides a complete workflow that can be applied toboth airborne IP and DC-IP data. Chapter 4 considers the second question. Here,the workflow is applied to both synthetic and field examples of airborne IP data.In Chapter 5, the third and last question is treated by applying the workflow to asynthetic DC-IP example.In this introductory Chapter 1, I will present mechanisms of IP phenomenon, asimple mathematical model for IP, how the IP effect can be simulated in TEM datawith Maxwell’s equations, and the workflow to restore polarization informationfrom EM data including IP effects.5Figure 1.1: Cole-Cole resistivity as a function of frequency [Pelton et al.,1978]. Black and red lines indicate real and imaginary part of the Cole-Cole resistivity. Used Cole-Cole parameters are: r0=1.25 Wm, h=0.2,t=0.1 s, and c=0.5.Figure 1.2: Conceptual diagram of (a) half-cycle input current waveform and(b) measured IP responses at surface electrodes.6Table 1.1: Cole-Cole parameters from different rocks (obtained from Peltonet al. [1978]).r0(Wm) h t (s) c (0< c<1)Porphyry (dry) 1 ⇥ 101 - 1 ⇥ 103 0.1 - 0.5 1⇥ 105 - 1⇥ 100 0.1 - 0.6Porphyry (wet) 1 ⇥ 100 - 1 ⇥ 104 0.1 - 0.8 1⇥ 102 - 7⇥ 104 0.1 - 0.5Magnetite 1 ⇥ 101 - 1 ⇥ 103 0.1 - 1.0 8⇥ 104 - 3⇥ 100 0.1 - 0.6Pyrrhotite 1 ⇥ 100 - 1 ⇥ 103 0.3 - 0.8 3⇥ 100 - 1⇥ 105 0.1 - 0.5Massive sulfide 1 ⇥ 102 - 1 ⇥ 103 0.6 - 0.95 8⇥ 104 - 2⇥ 100 0.1 - 0.4Graphite 1 ⇥ 102 - 1 ⇥ 103 0.75 - 0.98 3⇥ 101 - 8⇥ 103 0.1 - 0.571.2 Polarization mechanismsIn this section, I want to provide some insight about polarization mechanismssince this thesis is devoted to recovering distributed polarization parameters of theground from various TEM data. Polarization mechanisms are complex and I willnot investigate them in detail nor do the specific details of mechanisms really mat-ter to my thesis.Induced polarization (IP) is caused by the transport and accumulation of chargecarriers (ions and electrons) in earth materials (e.g. solid rock or unconsolidatedsoils) due to an applied electrical field. From an abstract point of view, the accumu-lation of charge carriers will oppose the current flow in the material. This createsan electric dipole that results in an overvoltage in potential measured at the surfaceand is also associated with a reversed current (in the opposite direction to the ap-plied electric field). Fig. 1.3(e) illustrates this induced polarization phenomenon.Consider a pore throat with an ionic fluid. At the initial state, the system (rocks andfluids) is electrically neutral. When an electrical field is applied (in the case by abattery or a generator), ions will move through the rock, and accumulate at the porethroat, these charges accumulate and oppose the current flow in the rock. If voltageis measured as a function of time, it will increase in time, and reach to the steadystate as shown in Fig. 1.3 (d). After the applied electric field is removed, polar-ized ions will revert back to the initial state, and this current flow is in the reverseddirection of the applied electric field (see how the positive charges move). Thiseffectively generates a decaying voltage as shown in Fig. 1.3. At the microscopiclevel, intuition can come from thinking about the surface area-to-volume ratio ofpore material. As the surface area-to-volume ratio increases, so does the numberof microscopic pores that can accumulate charges. This results in stronger inducedpolarization effects. Fig. 1.4 (a) and (b) show microscope images of clay mineraland sandstone. The surface area-to-volume ratio of the clay minerals is much big-ger than that of the sandstone; this will results in much stronger polarization effectsfor the clay minerals when an electric field is applied.Having a simple intuition is good; however, recognizing various polarizationmechanisms are still important because this will allow us to identify potential ap-plications of IP. Thus, I provide a brief overview of some important IP mechanisms8Figure 1.3: Polarization phenomenon at a pore throat: (a) Initial state, (b)While electric field is applied, and (c) Right after removal of the electricfield. (d) Input current. (e) Overvoltage at on-time. (f) Overvoltage atoff-time.below 1 MHz. These can be classified as five different categories:1. Electrode polarization2. Stern layer polarization3. Diffuse layer polarization4. Membrane polarization9Figure 1.4: Microscope images of rocks. (a) Illite (a clay mineral) with sur-face area-to-volume ratio of 100 m2/g(1000 times greater than sand-stone). (b) Quartz overgrowths in sandstone with surface area-to-volume ratio of 0.1 m2/g. Modified from [GPG, 2018].5. Maxwell-Wagner polarizationElectrode polarization can be observed in the presence of disseminated conduc-tive minerals such as pyrite, chalcopyrite, pyrrhotite, and magnetite [Pelton et al.,1978, Wong, 1979, Revil et al., 2015]. This has been very useful model for min-ing applications such as finding porphyry deposits. Fig. 1.5 illustrates before and10after the electric field is applied to a metallic particle in a solution. The electrodepolarization will mainly depends upon type, volume and geometry (e.g. diameterand shape) of the metallic particle, and it can explain both low and high frequencypolarization effects [Pelton et al., 1978, Wong, 1979, Revil et al., 2015]. In mostmining applications, low frequency polarization effects are observed in DC-IP data(< 10 Hz).Stern (2)[Reference], Diffuse layer (3)[Reference], and Membrane (4)[Refer-ence] polarization are related to the electrical double layer (EDL) structure; theyexplain polarization mechanisms of porous media, but in particular for clays. Fig.1.6 shows the EDL composed of Stern layer (fixed) and Diffuse layer. An insulat-ing grain is in the the electrolyte containing a dissolved salt and hence both cationsand anions exist in the solution. Cations are bounded in the Stern layer, whereasboth cations and anions can move in the Diffuse layer. Fig. 1.7 illustrates thethree polarization mechanisms:(2)-(4). Both the Stern and the Diffuse layer canbe polarized with the applied electric field, and after removal of the electric field,polarized charges diffuse [Marshall and Madden, 1959, Revil, 2013]. These threemechanisms can contribute to low frequency polarization effect (< 10 Hz; Revil[2013]). In Fig. 1.7(b) due to greater mobility of cations than anions, they can beaccumulated on the right hand side as a consequence of the applied electric field.After removal of the electric field, accumulated charges diffuse backward. Thesetwo phenomenon acting together is sometimes called membrane polarization, andthis can contribute to both low and high frequency polarization effects [Bu¨cker andHo¨rdt, 2013]. All three mechanisms will depends upon structure of the pore, andhow cations can be exchanged on the surface of the grain. Hence, both cation ex-change capacity (CEC) and the surface area-to-volume ratio will be key parametersfor these mechanisms.The Maxwell-Wagner polarization [Chen and Or, 2006] is an interfacial polar-ization due to the discontinuity of displacement currents in a multiphase systemwith discontinuities of the dielectric permittivity and/or electrical conductivity atthe interface between the different phases [Kemna et al., 2012]. Fig. 1.8 illustratesthe Maxwell-Wagner polarization. When an electric field is applied to a mediumhaving three different phases: solid (rock), water, and air. This mechanism is con-trolled by the tortuosity of the different phases, their volume fractions, and the11conductivity and permittivity of the different phases. The Maxwell-Wagner polar-ization is mainly responsible for polarization phenomena at frequencies above 10Hz [Chen and Or, 2006].To summarize, from a conceptual view point, the polarization effect is ex-plained as buildup of polarization charge which disturbs flow of the electrical cur-rent. When the applied electric field is removed, the accumulated polarizationcharges diffuse backward, which generates a current that has opposite direction tothe electric field. Five different mechanisms of IP are introduced here: (1) Elec-trode polarization, (2) Stern layer polarization, (3) Diffuse layer polarization, (4)Membrane polarization, and (5) Maxwell-Wagner polarization. (1)-(4) were re-lated to electromigration of ions or electrons, whereas (5) was mainly due to in-terfacial discontinuities in both conductivity and dielectric permittivity. (1)-(4) cancontribute to low frequency polarization effects (< 10 Hz), whereas (5) mostlycontribute to high frequency polarization effect (> 10 Hz). (1) and (4) can alsocause high frequency polarization effects.Figure 1.5: Electrode polarization: (a) Initial state, (b) While electric field isapplied, and (c) Right after removal of the electric field. Modified fromRevil et al. [2015].12Figure 1.6: Electrical double layer at the surface of insulating minerals. Xindicates active sites at mineral surface, which attracts metallic cationsM+. M+ is attached to the active sites, and bounded by the OuterHelmholtz Plane (OHP). The OHP separates the stern layer from thediffuse layer. A- indicates anion. Obtained from Revil [2013]Figure 1.7: Polarization mechanisms related to Electrical Double Layer(EDL): (a) Initial state, (b) While electric field is applied, and (c) Rightafter removal of the electric field. Modified from Revil [2013].13Figure 1.8: Maxwell-Wagner polarization: (a) Initial state. For a rock con-taining three different phases: solid, water, and air that have differentconductivity and dielectric permittivity. (b) Polarization charge can bebuilt when an electric field is applied. (c) Charges diffuse away whenthe electric field is removed, generating opposite direction of the currentflow to the applied electric field. Modified from Chen and Or [2006].141.3 Complex Resistivity and ConductivityMathematical representations of IP effects are often defined in the frequency do-main with the resistivity represented as a complex function. An extremely simple,but a most used model, is the Cole-Cole model. One representation due to Peltonet al. [1978] can be expressed asZ(w) = R0h1h⇣1 11+(ıwt)c⌘i, (1.3)where R0 is the resistance at zero frequency, h is the chargeability, t is the timeconstant, c is the frequency exponent. As frequency goes to infinity, the above com-plex impedance value (Z(w)) asymptotes to R• = R01h . This yields an expressionof the chargeability ash = R0R•R0. (1.4)For a sample rock with area A (m2) and length l (m), the complex resistivity, r(w)can be simply obtained byr(w) = Z(w)Al. (1.5)Pelton’s Cole-Cole model can also be expressed in resistivity form:r(w) = r0h1h⇣1 11+(ıwt)c⌘i, (1.6)and the same equation can be written in conductivity form:s(w) = r(w)1 = s•s• h1+(1h)(ıwt)c , (1.7)where subscripts 0 and • respectively stands for a value at zero and inifinte fre-quency. Here, the chargeability can be defined as either using either resistivity andconductivity:h = r0r•r0=s•s0s•. (1.8)Note that eqs. (1.6) and (1.7) is essentially the same as the equation for re-sistivity when the dielectric permittivity is complex [Cole and Cole, 1941]. For15conductivity, the original Cole-Cole model can be written as:scc(w) = s•+s0s•1+(ıwtcc)ccc= s•s• h1+(ıwtcc)ccc . (1.9)I shall refer to eq. (1.7) a Pelton’s model (PM) and eq. (1.9) a Cole-Cole model(CCM) following Tarasov and Titov [2013]. Here, tcc and ccc differ from t and c ofPM. Noticeable difference lies in the denominator: PM includes (1h), but CCMdo not. There are some preferences among different researchers for CCM and PM,and Tarasov and Titov [2013] throughly investigated pros and cons of each model.However, it is not the focus of my PhD research, and I choose Pelton’s model,which extensively used in mining applications, throughout my thesis. Hence whenI am referring to the ‘Cole-Cole model’, that indicates PM.1.4 Lab-scale IP MeasurementsTo understand how IP effects can be measured, I consider a lab-scale IP measure-ment system used in Geological Survey of Canada [Enkin et al.] as shown inFig. 1.9. A core sample cut in a standardized size and shape is located betweentwo cylinder-shaped copper electrodes. To prevent possible charge build-up on thecontacts of electrodes and rock, two paper membranes, soaked with copper sul-phate, are attached between the rock sample and copper electrodes. A source and avolt-meter are attached to the copper electrodes. The source generates a sinusoidalcurrent ranging from 0.025Hz -1MHz, and the volt-meter measures the compleximpedance at each transmitting frequency.Fig. 1.10(a) shows the measured complex impedance from a core sampleobtained from a Highland Valley region at British Columbia, Canada. This is aporphyry-copper deposit. The black solid lines show the observed amplitude andphase of the measured impedance. The observations are assumed to be the sum ofa complex impedance associated with the rock as well as other instrument effectssuch as capacitive effects associated with the copper electrodes and electrode ef-fects associated with the contacts between the copper sulfate paper and rock. Theblack dots are the predicted data associated with all of these effects. The contribu-16tion to the complex impedance that arises from the rock is marked as stars. Thisis the response of a Cole-Cole model, with parameters shown on the plot. For fur-ther details about extracting the impedance of the rock from the observations, thereader is referred to Enkin et al.. The Cole-Cole impedance reproduces the lowerand mid-frequency parts of the observations. The chargeability of the porphyry isestimated to be h = 0.16 and the time constant t = 2.6⇥103s.As a second example, I present the complex impedance of a core sample ob-tained from the DO-27/18 kimberlites. It is shown in Fig. 1.10b. This sampleshows significantly different Cole-Cole parameters especially for the time con-stant (t). For the kimberlite t ' 4⇥ 105s which is very small compared to thatfrom the porphyry deposit. This difference may be due to different polarizationmechanisms. The first sample from a porphyry deposit contains significant amountof pyrites which are possibly more related to metallic polarization. On the otherhand, the sample from the kimberlite deposit (crater facies) contains a considerableamount of clay minerals with high porosity, and might be filled with water; morerelevant to the membrane polarization and Maxwell-Wagner polarization. Nev-ertheless, the exact understanding for the cause of the IP response is not crucialfor my research. An important point is that at the hand-specimen scale, these tworocks are associated with different IP parameters, and a geophysical experimentmeasured distinctly different Cole-Cole parameters may thus allow us to distin-guish between them.con17Figure 1.9: Lab scale IP measurement system in Geological Survey ofCanada (GSC). The photo is taken from Cowan [2015].18Figure 1.10: Complex impedance measurements using the Geological Sur-vey of Canada system [Enkin et al.].191.5 Modelling Maxwell’s EquationsIP responses in various types of time-domain electromagnetic (TEM) surveys aregoverned by Maxwell’s equations. When the electrical conductivity of the earth isfrequency-dependent, then current density, ~J(w), in the frequency domain is~J = s(w)~E, (1.10)where ~E is the electric field (V/m) and s(w) is the frequency-dependent conduc-tivity (S/m), which is complex-valued. The complex conductivity is plotted in Fig.1.11(a) as a function of frequency. In the time-domain, eq. (1.10) can be written as~j = s(t)⌦~e, (1.11)where⌦ indicates a causal convolution, and~e is the electric field. To compute TEMsignals including IP effects, this convolution between s(t) and~e needs to be takenaccount when solving Maxwell’s equations. I first illustrate complex conductivity.For problems in the time-domain, it is convenient to work in the Laplace domain.The complex conductivity model in the Laplace domain can be expressed ass(s) = s•+4s(s), (1.12)where s= ıw is Laplace transform parameter and s• is the conductivity at infinitefrequency. The inverse Laplace transform of this yields:s(t) =L 1[s(s)] = s•d (t)+4s(t), (1.13)where L [·] is a Laplace transform, d (t) is a Dirac-delta function, and 4s(t) =L 1[4s(s)], which look like Fig. 1.11(b). Note that s(t) includes a delta func-tion at t=0, and a decaying function after t=0 with negative amplitude. Its convo-lution with ~e is well-defined. Consider a special case for a Cole-Cole model (eq.1.9) with c=1,4s(s) can be expressed as4s(w) =s• h1+(1h)(st) , (1.14)20In the time-domain,4s(t) can be expressed as4s(t) =s• h(1h)t e t(1h)t . (1.15)Now I turn our attention to Maxwell’s equations. Here I first define EM fieldswithout and with IP effects, then introduce how IP data can be defined. The electro-magnetic (EM) fields without IP effects are referred to as fundamental fields [Smithand West, 1989], and the corresponding Maxwell’s equations can be written as~—⇥~eF =∂~bF∂ t, (1.16)~—⇥ 1µ~bF s•~eF = ~js. (1.17)Here the superscript F on EM fields in the above equations refers to fundamentalfields,~b is the magnetic flux density (Wb/m2), ~js is the current source (A/m2) andµ is the magnetic permeability (H/m). Observed data will include both EM and IPeffects, and corresponding Maxwells equations can be written as~—⇥~e=∂~b∂ t, (1.18)~—⇥ 1µ~b~j = ~js, (1.19)Again, the total current ~j is a convolution of s(t) and~e(t):~j(t) =L 1[s(s)~E(s)] = s(t)⌦~e(t) =Z t0s(tu)~e(u)du, (1.20)where ~E(s) is an electric field in the Laplace domain. Observed fields includingboth EM and IP effects can be defined as ~e = ~eF +~eIP, ~b =~bF +~bIP and ~j =~jF +~jIP, where the superscript IP is induced polarization.The IP response, dIP, can be computed by subtracting the fundamental responsedF from the observation d:dIP = F [s(s)]F [s•] = ddF (1.21)21where F [·] is Maxwell’s operator. For example, F [s•] outputs EM fields usinga conductivity model s•. Observed and fundamental responses correspond tod = F [s(s)] and dF = F [s•], respectively. By solving both Maxwells equations, Ican evaluate both the observed and the fundamental responses, which allows me tocompute dIP. The subtraction in eq. (1.21) correspond to the EM-decoupling [Pel-ton et al., 1978, Routh and Oldenburg, 2001]. Note that when chargeability is zero,s(s) = s• and dIP = 0. For these simulations I use EMTDIP code developed byMarchant et al. [2014]. This code is capable of solving the full Maxwell equationsand hence of simulating both inductive current flow and galvanic current flow.To understand how both fundamental and IP fields behave in a TEM survey,a simple grounded DC-IP survey, using a single source, is considered. The cur-rent waveform includes 2100 ms on- and off-time, in the transmitter wire, currentflows in counter-clockwise direction. Fig. 1.12 shows a source wire path andthe chargeability distribution in the core region of the 3D mesh (that is, paddingcells are omitted). Chargeability of the chargeable block is 0.1. The conductivitymodel at the infinite frequency is a half-space conductivity, shal f = 0.05 S/m. Twoforward simulations are performed to compute observed (d = F [s(s)]) and funda-mental (dF = F [s•]) responses. By subtracting dF from d, the IP response (dIP)is obtained. Computed potential difference at a receiver location (where the twoelectrodes are marked as black dots in Fig. 1.12(a) is shown in Fig. 1.13. Immedi-ately after the current is switched on, significant EM induction arises due to Fara-day’s law, then as time proceeds the EM induction effect decays and the observedresponse converges to the steady-state. Similarly after the current switch-off, sig-nificant EM induction rises, and this decays as time passes. Here, the normalizedpotential difference is provided, which means V (t)/V0, where V (t) is the observedpotential difference and V0 is the potential at the DC limit.Although only decay curves at a single receiver locations are provided in Fig.1.13, electric field exists everywhere in the domain and seeing that in higher di-mension can promote understanding of both EM and IP effects in this time-domainDC-IP survey. Figs. 1.14(a) and (b) respectively show plan and section views ofthe electric fields at 6 ms (the top panel), 715 ms (the middle panel), and 2116 ms(the bottom panel). These three times are marked as the black dashed lines in Fig.1.13. At 6 ms, I can clearly recognize the induced EM signals from the current22wire path; those then disappear as at later when the voltage asymptotes to the DClimit. In the on-time, the IP response (red line) is much smaller than the fundamen-tal response (blue line), whereas in the off-time the IP responses are large enoughto be identified at late times. Fig. 1.15 show only off-time signals in the form oflog-scaled normalized potential. Note that “off-time” indicates a time after currentis turned off, and this time is set to zero (0 ms) when presenting off-time data.Dashed and solid lines differentiate the sign of responses, negative and positiverespectively. At early times (⇠ 10 ms), observed and fundamental responses arealmost coincident because EM induction is dominant at this time. At later time, therelative strength of the IP increases, and is dominant at 500 ms. Fig 1.16 showsthe electric fields at three off-times (6, 77, and 1717 ms). Significant EM inductioneffects are observed at 6 ms near the current wires. At 77 ms, this EM inductioneffect smoothly diffuses away except for the chargeable region, where polarizationcharges have been been built. At 1717 ms EM induction is decayed significantlyand hence IP is dominant. Here the chargeable block acts as a current dipole, whichhas opposite direction to the DC currents shown in Fig. 1.14, and this is consis-tent with the conceptual model described in Fig. 1.3(c). Grid of electrodes aredeployed on the surface, and the potential differences of the shortest dipole in x-direction are shown as a map in Fig. 1.17. Three time channels at (a) 6 ms, (b)715 ms, and (c) 2116 ms are shown here, and the plotted potential differences canbe considered as electric fields in x-direction (~e=—f , where f is the potential).Conversely, most of time-domain DC-IP surveys assume that IP-dominant timesare measured. However, when earth medium is highly conductive this assumptioncould be invalid. For instance, consider a simple case where the halfspace conduc-tivity shal f is increased to 0.5 S/m (previous one was 0.05 S/m). Fig 1.18 showsthe normalized potential at both on- and off-time channels. Compared to the pre-vious example, EM induction decays much slowly resulting in slower convergesto the steady-state in on-time. In off-time, as shown in Fig. 1.19, d and dF arealmost coincident indicating EM is always dominant at all times. Fig 1.20 showsthe measured potential difference at the off times.To summarize, in observation d, which includes both EM and IP effects, naturalseparation of EM and IP responses in time exists as shown in Fig. 1.15:23• Early-time: EM-dominant• Intermediate time: Both EM and IP are considerable• Late-time: IP-dominantAt the EM-dominant time, the observation includes minor IP effect hence this canbe used to recover conductivity distribution. Conversely at the IP-dominant time,EM induction effects in the observation are minor, thus these data can be used torecover a chargeability distribution. At intermediate time, both EM and IP effectsare significant and here removal of either the EM or IP response (often called EM-or IP-decoupling) is required; this depends upon which information (conductivityor chargeability) is desired. Although a simple time-domain DC-IP example isshown here, a physical concept that an observed TEM signal shows natural separa-tion of EM and IP effects in time is general enough to be extended to various TEMsurveys (e.g. airborne EM).Figure 1.11: Cole-Cole response in frequency domain (a) and time (b) do-main. The Cole-Cole parameters are s• = 102 S/m, h = 0.5, t = 0.1s, and c=1. The arrow indicates a Dirac-delta function (s•d (t)).24Figure 1.12: 3D chargeability model. (a) Plan and (b) Section views. Whitesolid line on (a) shows the wire path, and the arrow indicates the direc-tion of current in the wire when current is turned on. THe half-spaceconductivity is 0.05 S/m.25Figure 1.13: Time curves of the normalized potential difference at a receiverlocation (marked as black solid line with dots in Fig. 1.12(a)). Black,blue and red lines correspondingly indicate observed, fundamental andIP data. Solid and dashed lines distinguish positive and negative val-ues. Black vertical dashed lines at three time channels indicate timesat which the electric field distributions in Fig. 1.14 are provided. Thehalfspace conductivity is 0.05 S/m.26Figure 1.14: 3D electric field distributions at three different times in the on-time. Top, middle and bottom panel correspondingly show the electricfields at 6 ms, 77 ms, and 1717 ms. (a) Plan view at the surface and (b)Section view at Northing 0 m. The halfspace conductivity is 0.05 S/m.27Figure 1.15: Time decaying curves of the normalized potential differencewhen the current is turned off. Black, blue, and red lines respectivelyindicate observed, fundamental and IP data. Dotted and solid lines dis-tinguish positive and negative values. Halfspace conductivity is 0.05S/m.28Figure 1.16: 3D electric field distribution at three different times when thecurrent is turned off. Top, middle and bottom panel correspondinglyshow the electric fields at 6 ms, 715 ms, and 2116 ms. (a) Plan viewat the surface and (b) Section view at Northing 0 m. The halfspaceconductivity is 0.05 S/m.29Figure 1.17: Observed potential difference at three different off times. (a) 6ms, (b) 77 ms, and (c) 1717 ms. Halfspace conductivity is 0.05 S/m.Figure 1.18: Time curves of the normalized potential difference at a receiverlocation (marked as black solid line with dots in Fig. 1.12(a)). Halfs-pace conductivity is 0.5 S/m.30Figure 1.19: Time decaying curves of the normalized potential differencewhen current is turned off. Halfspace conductivity is 0.5 S/m.Figure 1.20: Observed potential difference at three different off times. (a) 6ms, (b) 77 ms, and (c) 1717 ms. Halfspace conductivity is 0.5 S/m.311.6 TEM-IP Inversion WorkflowThe previous examples show how IP and EM induction phenomenon affect EMdata. A primary goal of this thesis is to recover distributed IP information in 3Dfrom observed TEM data, which include both EM and IP signals. An effectivestrategy is required to achieve this goal, and three main options can be considered:• Formulate an inverse problem to find a complex conductivity: s(x,y,z;w).For instance, in frequency domain, each frequency can be inverted separatelyor simultaneously, and recover a complex conductivity distribution [Kemnaet al., 2004, Commer et al., 2011].• Parameterize complex conductivity as a few parameters (e.g. Cole-Cole: s•,h , t , and c), then measured data can jointly be inverted for those multipleparameters [Fiandaca et al., 2012, Marchant et al., 2013, Xu and Zhdanov,2015].• Realize that the signal has an early (EM) part, a late (IP) part and an interme-diate zone of mixed signals as shown in Fig. 1.15. Isolate these by applyingEM-decoupling or IP-decoupling, and work with each independently. This iswhat has been done with most routine applications Oldenburg and Li [1994].For instance, a two-stage approach for DC-IP data inverts late on-time data(DC) to recover conductivity and uses that conductivity to invert early off-time data (IP) to recover a chargeability.The first and second options are advantageous given that they are taking ac-count both EM and IP effects simultaneously. However, in particular for a 3Dinversion, those choices will increase both non-uniqueness and computational costof the inversion dramatically, since model space is enlarged at least four times (e.g.Cole-Cole) and the computational cost of solving Maxwell’s equations is increaseddue to time-dependent conductivity [Marchant et al., 2014]. The third option isrobust and cost-effective, but this will have issues when earth medium includesconductive structures resulting in significant EM-coupling even at the latest timechannel, or there are very early-time IP effects that make conductivity estimationdifficult.32Each option has pros and cons, but basically depends on how observed data,d are represented. Here d can be any types of EM fields or fluxes that can bemeasured. Two choices ared = F [s(x,y,z;w)], (1.22)d = F [s•]+dIP. (1.23)eq. (1.22) indicates that the observation can be simulated as a non-linear function,which takes a complex conductivity s(x,y,z;w); here both EM and IP effects arecoupled into the observation. In contrast as described in eq. 1.23, the observationcan also be separated as EM and IP, and treated separately. For this case, theIP datum, dIP can be linearized as described in Seigel [1974], Oldenburg and Li[1994]: dIP = Jh˜ , which yieldsd = F [s•]+ Jh˜ , (1.24)where J is a sensitivity function and h˜ is a pseudo-chargeability 1; its definitionwill be described in Chapter 2. I prefer the second view point (eq. 1.24), whichdecomposes two different physical effects, and provides a way to linearize dIP.From a learning point of view, separating EM and IP parts in the observation,investigating how each physical effect behaves, and even linearizing the IP part forvarious TEM surveys, promotes physical understanding of each phenomenon.Therefore, I choose the third option and modify the two-stage approach so thatit could handle EM-coupling in the observation. This new workflow is named theTEM-IP inversion workflow, and its procedures are as follows:• The first step is to estimate a background conductivity. To achieve this, thelate on-time or the early off-time data can be inverted by respectively usingthe DC or TEM inversion algorithms to recover conductivity, sest . Eitherway, a selection of data for inversion assumes that any IP effects in the dataare insignificant.• Using the estimated conductivity sest , a TEM response can be computed;1I note that pseudo-chargeability, h˜ is different from Cole-Cole chargeability, h = s•s0s•33this is an approximate fundamental response: F [sest ]. By subtracting thisfrom the observation, the IP data can be defined:dIP[sest ] = dF [sest ],and this corresponds to EM-decoupling.• Then develop a linearized function of dIP:dIP = Jh˜and this linear equation is used to recover a pseudo-chargeability by invertingthe obtained IP data, dIP[sest ].341.7 Thesis OutlineThe main goals of this thesis are:• Develop a methodology by which IP information about the 3D earth’s sub-surface can be extracted from inductive and grounded sources.• Contend with important aspects of EM-coupling.• Apply the technique to field data sets from TEM surveys.As introduced in Section 1.6, the TEM-IP inversion workflow has three mainsteps: (a) Invert TEM data to recover conductivity distribution. (b) Subtract thepredicted EM response (obtained by the recovered conductivity) from the obser-vation, to obtain IP data. (c) Using a linearized IP function, invert the IP data torecover a distributed chargeability information. The application of the workflow isfocused on two types of TEM surveys: (a) AEM and (b) DC-IP surveys. Chapters2 and 3 presents how the IP signal from an inductive source is linearized as a func-tion of a pseudo-chargeability, and how a 3D IP inversion is developed to recover3D distribution of the pseudo-chargeability. The TEM-IP inversion workflow isapplied to both the AEM and DC-IP surveys. In Chapter 4, I concentrate upon theinversion of synthetic and field AEM data sets. Here, one synthetic example (Sec-tion 4.1) and two field examples: the Mt. Milligan porphyry deposit (Section 4.2)and Tli Kwi Cho (TKC) kimberlites (Section 4.3) are presented. The workflow isgeneral and hence it can also be applied to grounded sources. Chapter 5 presentsthe application of the workflow to the gradient array DC-IP survey, and emphasizesthe use of early-time TEM data to recover a better conductivity model, and its usefor the EM-decoupling for the late off-time IP data.Collectively, my specific contribution to science includes:1. Developing the TEM-IP inversion workflow to recover both conductivity andchargeability information in 3D from various TEM surveys.2. Deriving and verifying a linear IP functional for inductive sources and inparticular for AEM surveys.353. Recovering both conductivity and chargeability information from AEM datathat include negative transients. In particular, my work on a kimberlite de-posit (TKC) shows how 3D conductivity and chargeability information canbe an important component of generating a geophysical (and alternativelygeologic) rock model from AEM data.4. Removing EM-coupling in the DC-IP data to obtain better 3D chargeabilityby using early-time TEM signals, which are often discarded.36Chapter 2LinearizationEarth materials are chargeable, hence polarization charges and currents can be gen-erated whenever an electric field is applied. Hence any EM data set can include IPeffects, and effectively interpreting these IP data has been an important issue. Fora DC source (grounded source without EM-coupling), a linearized form of IP datais well-known [Seigel, 1959], and successfully used as a forward kernel for IP in-versions [Oldenburg and Li, 1994, Li and Oldenburg, 2000]. Seigel’s basic idea isthat polarization effects can be considered as a small perturbation in conductivity,which can be written ass = s•(1 h˜). (2.1)Here s includes polarization effects, but s• does not; recall, s(s) = s• (withs = ıw) when h = 0. Here w is the angular frequency. Pseudo-chargeability h˜ isoften assumed to be small (h˜ ⌧ 1). Note that the pseudo-chargeability is not thesame as the chargeability (h = s•s0s• ), but it is similar in the sense that both ofthem can be thought as a small fraction of conductivity change. Now based uponeq. ( 2.1) and following Seigel [1959], IP data can be linearized using Taylor’sexpansiondIP = Fdc[s•(1 h˜)]Fdc[s•]⇡ d Fdc[s•]d log(s•)h˜ = Jdc[s•]h˜ , (2.2)37where Fdc[·] is a DC operator, which takes conductivity and computes a DC re-sponse, and Jdc[s•] = d Fdc[s•]d log(s•) is a sensitivity function.Extension of this linearization to an inductive source is not available, and thislimits the capability to apply the TEM-IP inversion workflow to inductive sourceIP (ISIP) data to restore polarizable information. The challenge posed by the useof inductive sources is that steady-state electric fields are not established inside theearth as they are for the DC source. At any location in the earth the electric fieldwill increase to a maximum value and then decrease as the electromagnetic (EM)wave diffuses through. The EM fields at any position and time depending upon theconvolution of the electric field with the time-dependent conductivity of the rock.Unraveling these complexities, and linearizing ISIP data are issues that I address inthis Chapter. Systematic numerical verifications will support validity of the deriva-tion, and these will be focused upon the airborne EM (AEM) geometry. Althoughthe focus is linearizing the ISIP, the developed IP function will be general enoughto handle IP data from various TEM surveys (e.g. AEM and DC-IP surveys).2.1 Pseudo-chargeabilityThe goal of this chapter is to linearize ISIP data as a function of pseudo-chargeability.Here I define what the pseudo-chargeability is, and build a foundation for the fol-lowing sections.Polarization phenomenon is effectively described by a complex conductivitymodel, and its general form (with s= ıw) is a good starting point:s(s) = s•+4s(s). (2.3)Here w is the angular frequency. Depending upon the choice of a complex con-ductivity model,4s(s) can be defined differently; for instance4s(s) of Pelton’sCole-Cole model is:4s(s) = s•h1+(1h)(st)c . (2.4)Fig. 2.1(a) shows an example Cole-Cole model as a function of frequency. Thetime-domain form of complex conductivity can be obtained using inverse Laplace38transform:s(t) =L 1[s(s)] = s•d (t)+4s(t), (2.5)where d (t) is Dirac delta function, andL 1[·] is inverse Laplace transform opera-tor. Note that only a causal function is considered here, which is only defined whent  0.The chargeability h ish = 1s•limt!•L1[4s(s)s] =s•s0s•. (2.6)Here 4s(s) is divided by s to obtain a step-function in time. It is convenient todefine an time-dependent impulse reponse, h˜ I(t), ash˜ I(t) =4s(t)s•. (2.7)which transforms to4s(t) =s•h˜ I(t), (2.8)Note that the intrinsic chargeability, h , is not time-dependent but the impulsepseudo-chargeability, h˜ I(t) is time-dependent. The Cole-Cole response in time-domain is shown in Fig. 2.1(b). The arrow at t=0 s indicates s•d (t), which is adelta function, and after t=0 s, s(t) =4s(t). When c = 1, the inverse Laplacetransform can be evaluated and the corresponding h˜ I(t) is:h˜ I(t) = h(1h)t e t(1h)t . (2.9)Ohm’s law describes how currents are generated when an electric field is ap-plied to a conductive medium, and in Laplace domain it can be written as~J(s) = s(s)~E(s) (2.10)The inverse Laplace transform of this yields~j(t) =L 1[~J(s)] = s(t)⌦~e(t), (2.11)39where ⌦ is a causal convolution.I define the IP currents, ~jIP, as the total current (~j) minus the fundamentalcurrent (~jF ):~jIP = ~j~jF . (2.12)Using eqs. (2.5) and (2.11) I obtain~jIP = s•~e IP+~jpol, (2.13)where the polarization current (~jpol) is~jpol(t) =4s(t)⌦~e(t). (2.14)If the electric field has different characteristics for inductive and groundedsources this will generate different features in the polarization current. Considertwo cases: (1) a DC source with grounded electrodes and (2) an inductive source.The first case corresponds to the usual approache for interpreting DC-IP data [Seigel,1959, Oldenburg and Li, 1994], and the second is associated with ISIP. The polar-ization current is a convolution between 4s and ~e (eq. 2.14), and Fig. 2.2(a)illustrates this process for case (1) with known~e. Here the electric field is instanta-neously turned on then off, which emulates an electric field from case (1) withoutEM induction effects. This electric field is convolved with 4s , and outputs the~jpol , which has the opposite sign of the electric field (due to negative 4s ). Asshown in Fig. 2.2(b), a similar process is applied for case (2), and here ~e hasa different time behavior because case (2) excites the earth inductively. It startsfrom zero, increases to a peak then decays. The resultant ~jpol shows a similartime behavior to the applied electric field. Comparison of ~jpol due to cases (1) and(2) illustrates different polarization buildups when the earth’s subsurface is exciteddifferently. The main difference in between cases (1) and (2) is that the absenceof steady-state electric field for case (2) resulting in a more dynamic polarizationprocess for case (2).To capture this difference in a linearized IP function, pseudo-chargeability,40h˜(t) is defined ash˜(t) =~jpol(t)~j re f, (2.15)where the reference current, ~jre f is defined as~j re f = s•~e re f . (2.16)Here ~e re f is the reference electric field, and a choice of ~e re f is described below.The pseudo-chargeability defined in eq. (2.15) is the ratio of the polarization cur-rent to the reference current. The pseudo-chargeability is a small quantity and thisplays an essential role in our linearization. To evaluate the pseudo-chargeability,a reference current or reference electric field ~e re f is required to be identified; inmy definition (eq. 2.15) this is independent of time. For DC-IP, the value of theelectric field achieved when there is no IP is chosen, that is the value shown in Fig.2.3(a). For the inductive source the peak electric field is selected as shown in Fig.2.3(b).Each pixel in the earth has its own reference electric field and time that thepeak value occurs thus both ~e re f and t re f have a 3D distribution. For both DC-IP and ISIP cases, the choice of the reference electric field can mathematically bepresented as~e re f =~eF(t)⌦d (t t re f ). (2.17)The reference time for the DC-IP case can be any time in the on-time. By rearrang-ing eq. (2.15), I obtain~jpol =~j re f h˜(t). (2.18)This states that the polarization current has an opposite direction to the referencecurrent, and is proportional to the pseudo-chargeability. This reversed direction ofthe current in a chargeable medium results from the negative values of the time-dependent conductivity when t >0 as shown in Fig. 2.1(b). This conceptual modelabout the polarization current shown in eq. (2.18) is consistent with Seigel [1959]’sresult. Note that for any pixel, even if~e re f attains the same value for an ISIP surveyas for an DC-IP survey, the pseudo-chargeability resulting from an ISIP surveywill be less than that from a DC-IP survey. To illustrate, let us assume e0 and4s41shown in Fig. 2.2(a) and (b) are the same, but the area for the DC source willalways be greater than the inductive source resulting in greater ~jpol; that is greaterpseudo-chargeability. Considering the linearizations for DC-IP problems workedwell when h˜ is small (< 1), I can infer from this that linearization techniquesshould be successful in ISIP problems, which have even smaller h˜ compared toDC-IP problems.Figure 2.1: Cole-Cole response in frequency domain (a) and time (b) domain.The Cole-Cole parameters are s• = 102 S/m, h = 0.5, t = 0.01 s, andc=1. The arrow shown in Fig. 2.1b indicates a delta function (s•d (t)).2.2 Linear IP FunctionFollowing from the methodologies in DC-IP, our goal is to express the IP response,dIP as a linear function of the pseudo-chargeability, h˜(t). That is I wish to writedIP(t) = Jh˜(t), where J is a yet to be determined sensitivity function which isindependent of time. In doing this first consider a general EM system which isapplicable to grounded or inductive sources. For any volume pixel in the earth theamplitude and direction of the electric field can vary dramatically in time and thusthe IP charging process can be complicated. However, if substantial polarizationcurrents are developed, an assumption can be made that there was a sufficientlylarge electric field in a predominant direction to generate them. Although the di-rection of the electric field is constant the amplitude varies with time.Let~e(t) be approximated as~e(t)⇡~e re f wˆ(t), (2.19)42Figure 2.2: Convolution of 4s and ~e resulting in polarization current, ~jpol .(a) DC source. (b) Inductive source. 4s and h˜ I are respectively definedin eqs. ( 2.8) and ( 2.7).Figure 2.3: Conceptual diagram for the amplitude of the fundamental electricfields. (a) DC source. (b) Inductive source.where wˆ(t) is defined as:wˆ(t) = P0[wre f (t)]. (2.20)Here a projection P0[·] of an arbitrary time function, f (t) isP0[ f (t)] =(f (t) f (t) 00 if f (t)< 0,(2.21)43andwre f (t) =~eF(t) ·~e re f~e re f ·~e re f . (2.22)wre f (t) is a dimensionless function that prescribes the time history of the electricfield at each location along the direction of the chosen reference electric field,~e re f .According to this definition, negative values of wre f (t) are set to zero in accordancewith our conceptual model that polarization currents have an opposite direction tothe reference current (eq. 2.18). The pseudo-chargeability is redefined ash˜(t) = h˜ I(t)⌦ wˆ(t). (2.23)The polarization current, ~jpol shown in eq. (2.14), can be approximated with eq.(2.7) as~jpol(t)⇡h˜ I(t)⌦ wˆ(t)~j re f . (2.24)Substituting into eq. (2.13) yields~jIP(t)⇡ s•~e IP(t) h˜ I(t)⌦ wˆ(t)~j re f (2.25)and this yields~jIP(t)⇡ s•~e IP(t)~j re f h˜(t). (2.26)The second term,~j re f h˜(t) corresponds to polarization currents. In particularfor ISIP, the first term, s•~e IP(t) is usually omitted [Smith et al., 1988]. This wasbecause Smith et al. [1988] were mostly interested in chargeable targets that weresignificantly conductive compared to the background. However, if the conductivityof the chargeable target is similar to that of the background, the first term couldbe important. That first term is included here and the conditions in which it isimportant will be explored.Because the reference current is static, any time-dependence in the polarizationcurrents is encapsulated in the pseudo-chargeability. The buildup and decrease ofpolarization currents is a slow process and I assume therefore that this process doesnot produce induction effects ( ∂~bIP∂ t ⇡ 0) and hence the IP electric field is assumedthat it can be written as~e IP ⇡~e IPapprox =~—f IP, (2.27)44where f IP is the electrical potential for IP. By taking the divergence of eq. (2.26),I obtain— ·~jIP(t)⇡ — ·s•~e IP(t)— ·~j re f h˜(t). (2.28)Here — ·~jIP = 0 because — ·~j = — · (~jF +~jIP) = — ·~js, where — · jF = — ·~js.Substituting~e IP with eq. (2.27) in eq. 2.28 and carrying out some linear algebra,I obtainf IP(t)⇡[— ·s•~—]1— ·~j re f h˜(t). (2.29)By applying the gradient I obtain~e IPapprox = ~—⇣[— ·s•~—]1— ·~j re f h˜(t)⌘. (2.30)Thus, the electric field due to the IP effect can be expressed as a function of h˜(t)in time. This form is also applicable to the DC-IP case, and the same form used forthe conventional case (eq. 2.2) can be obtained from eq. ( 2.29):f IP ⇡ d fs•d log(s•)h˜ , (2.31)with d fs•d log(s•) = [— ·s•~—]1⇣— ·s•~—fs•⌘and fs• is a DC potential with s•.Here the reference electric field ~e re f is set to ~—fs• . For the DC case, the timehistory of the electric field is basically the same as the input current waveform.Hence, when a step-off function: 1 u(t) is used for the current waveform, thepseudo-chargeability ish˜(t) = h˜ I⌦ (1u(t)). (2.32)For a Cole-Cole model with c=1, this can be explicitly written assamh˜(t) = het(1h)t , (2.33)which is a similar result to Ho¨rdt et al. [2006] and Yuval and Oldenburg [1997].For inductive sources, often data are either~b or its time derivative and hence~bIPor its time derivative needs to be computed. For this, I first compute ~jIP then usethe Biot-Savart law. By substituting eq. ( 2.30) into eq. ( 2.26), the approximated45IP current density, ~jIPapprox can be expressed as~jIP(t)⇡ ~jIPapprox = S¯~j re f h˜(t), (2.34)where S¯ is a matrix operator comprised of to terms:S¯= s•~—[— ·s•~—]1— ·I¯ (2.35)and I¯ is an identity tensor. Applying the Biot-Savart law I have:~bIPapprox(~r; t) =µ04pZWS¯~j re f (~rs)⇥ rˆ|~r~rs|2 h˜(t)d~rs, (2.36)where~rs indicates a vector for a source location, and rˆ = ~r~rs|~r~rs| . If s•~eIP is omit-ted in ~jIP then the tensor, S¯, becomes I¯. In this situation, the IP current is thesame as the polarization current, and it always has an opposite direction to the ref-erence current. This reversed current, along with the Biot-Savart law, provides aphysical understanding about the negative transients in AEM data when the earthis chargeable.Observed data are often the time derivative of~b, hence by taking time derivativeof eq. (2.36), I obtain ∂~bIPapprox∂ t(~r; t) =µ04pZWS¯~j re f (~rs)⇥ rˆ|~r~rs|2⇣ ∂ h˜(t)∂ t⌘d~rs. (2.37)Here I have chosen to keep the minus signs in eq. (2.37) so that  ∂ h˜(t)∂ t is positivewhen h˜(t) is decaying in time. Accordingly, the IP datum is given by  ∂~bIP∂ t . Notethat for inductive sources the time history of the electric field wˆ(t) is not the sameas the input current waveform. Thus, h˜ defined in eq. ( 2.23) is dependent upon thetime history of electric field diffused in the earth medium for ISIP, whereas it onlydepends polarization parameters and the current waveform for DC-IP as shown ineq. ( 2.32).The IP fields shown in eqs (2.30), (2.36) and (2.37) are linear functionals of h˜46and the equations for a single time channel can be discretized in space asdIP = Jh˜ , (2.38)where J is the corresponding sensitivity matrix. In particular when the observeddatum is the time derivative of~b, the linear relationship can be written asdIP = J(∂ h˜∂ t⌘. (2.39)A detailed description for the discretization of the linearized kernel is shown inAppendices C.1 and C.2. The representation in eq. (2.38) is valid for bothgrounded and inductive sources but the two assumptions: (a) ~e ⇡ ~e re f wˆ(t) and(b)~e IP ⇡~—f IP need to be tested numerically for the case of inductive sources.472.3 Numerical ExperimentsFor numerical experiments I concentrate upon the AEM survey. This choice ismade because of the observed negative transients that are direct indicators of po-larization effects [Weidelt, 1982], and the extensive use of this survey by industry[Smith and Klein, 1996, Kratzer and Macnae, 2012, Kaminski and Viezzoli, 2017,Kang et al., 2017]. I begin with a simple IP model composed of a chargeable blockin a halfspace as shown in Figure 2.4. Cole-Cole parameters of the block are h= 0.2, t = 0.005 s and c = 1. The conductivity of the halfspace, s1, is 103 S/m,whereas the conductivity of the chargeable body, s2, will be assigned different val-ues; s• is thus a 3D distribution. I consider three cases: (a) canonical (s2 = s1),(b) conductive (s2 = 102⇥s1) and (c) resistive models (s2 = 102⇥s1). The 3Dspace is discretized with 50⇥ 50⇥ 50 m core cells and the number of cells in thedomain is 41⇥41⇥40. The size of the chargeable body is 250⇥250⇥200 m andthe top boundary is located 50 m below the surface. The EMTDIP code [Marchantet al., 2014] is used to simulate EM responses that include IP effects. The survey,consisting of 11 soundings along each of 11 lines, is shown in Fig. 2.4(a). Dataare from a coincident-loop system and the flight height is 30 m above the surface;the radius of the loop is 10 m. A step-off source waveform is used and the rangeof the observed time channels is 0.01-60 ms. The observed responses can be thevertical component of~b or ∂~b∂ t . In these numerical experiments, first the observedresponses and the total currents are decomposed into fundamental and IP portionsto aid in the basic understanding of IP effects in AEM data. Then, the linearized IPfunction is systematically validated by computing the approximate IP currents andIP responses, and comparing these with the true values.2.3.1 IP DataUsing the EMTDIP code and carrying out two simulations: d = F [s(s)] and dF =F [s•], dIP is obtained by subtracting dF from d. Fig. 2.5 shows the d, dF , and dIPat a sounding location above the center of the chargeable body for (a) canonical,(b) conductive and (c) resistive models. Both bz and  ∂bz∂ t data are shown. TheIP effects are most noticeable for the conductive body and I turn attention to thisexample first. The IP response starts to significantly affect the observations near 0.648ms and the observed responses show a sign reversal near 1 ms. Beyond that timethe signal is dominated by the IP. The dashed line in Fig. 2.5(b) shows that afterturning off the source current, the IP current increases (as inferred by the magnitudeof the bz field) until about 1 ms and then decreases. I interpret this in terms ofcharging and discharging phases for the chargeable body and a vertical dashed linein the figure defines the two phases. In the charging phase at early times the EMeffects dominate and IP signals are not expected to be observed. In the dischargingphase, which occurs at later time, the IP effects may eventually dominate the EMeffects. The maximum of the bIPz corresponds to the zero crossing for  ∂bIPz∂ t butthe times at which the IP signal becomes dominant are delayed compared to bIPz .By comparing the observations with the fundamental fields I see that the IP signalcould be recognized in the bz data near 0.7 ms and near 2.0 ms in the  ∂bz∂ t data.The plots for the canonical and resistive bodies show that the time that separatescharging and discharging occurs earlier than that for the conductive body. This isa reflection that the fundamental currents reside for a longer time in a conductor.For the canonical body, a significant difference between the measured responsesand the fundamental fields occur about 0.9 ms for bz and about 2 ms for  ∂bz∂ t .The amplitudes of the IP responses are significantly smaller than those for theconductor. Lastly, there is little IP signal for the resistive body; the IP signal ismuch smaller than the fundamental response throughout the given time range. Thisis a consequence of the small fundamental currents in the resistor.The decay curves from a sounding location provide insight about the IP re-sponse but more is gleaned by looking at data from all sounding locations in theAEM survey. I focus on bIPz for the conductive block at selected time channels.Fig. 2.6 shows interpolated maps of the d, dF and dIP at (a) 0.86 ms and (b) 6.7ms which are respectively included in the charging and discharging times. For theconductive block, 0.86 ms is close to the peak time when the transition from charg-ing to discharging occurs, but it is still included in the charging time. At this time,the d are dominated by the dF and no negative values, which are the signature ofthe IP effect, are observed. Subtracting the dF however, yields a residual dIP datamap that has a strong negative. This example shows that our EM-decoupling pro-cedure can work satisfactorily for this simple model. At 6.7 ms, obtaining goodIP data is easier because the d already show negative values. There is still a weak49Figure 2.4: Plan (a) and section (b) views of the IP model. The solid line in(a) delineates the boundary of the IP body. Solid circles in (a) denotethe sounding locations. In (b) the conductivity s2 is variable so thatcanonical, conductive and resistive blocks can be examinedfundamental field and the subtraction process improves the dIP response. The dIPdata at 0.86 ms and 6.7 ms shown in Fig. 2.6 are of sufficient quality to be inverted.50Figure 2.5: Time decaying curves of the observations (d; black line), funda-mental (dF ; blue line) and IP (dIP; red line) responses. All three cases:(a) canonical, (b) conductive and (c) resistive are presented. Right andleft panels show bz and  ∂bz∂ t . The vertical black dotted line indicatesthe time at which the polarization field reaches its maximum value. Theflight height of the collocated transmitting and receiving loop is 30 mabove the surface.51Figure 2.6: Interpolated maps of observed (left panel), fundamental (middlepanel) and IP (right panel) responses. Two time channels at (a) 0.86 msand (b) 6.7 ms are presented. The white line contours a zero-crossing inthe observed response.522.3.2 Polarization CurrentsTo evaluate the polarization current shown in eq. (2.14) for the linear functional,I assumed ~e(t) ⇡ ~e re f we(t) and defined our reference current as ~jre f = s•~e re f .That yielded the polarization current to be ~jpol(t)⇡~j re f h˜(t). This requires thatthe polarization current has a direction antiparallel to the reference current, and thedirection is the same for all times. With this approximation the time dependencefor the polarization currents only occurs through the scalar h˜(t). The validity ofthe approximation is investigated by evaluating both reference and polarizationcurrents numerically. From eq. (2.17), a reference current can be considered asthe maximum fundamental current that occurred throughout the time history. Toevaluate polarization currents I rearrange eq. (2.13) as ~jpol = ~jIPs•~e IP.Here canonical and conductive blocks are used to compute ~jpol . Figs 2.7(a)and (b) show reference currents for the canonical and conductive blocks, respec-tively. A source is located at (-200 m, 0 m, 30 m) and marked as a white solid circlein the figure, where (·, ·, ·) refers to a point at (easting, northing, depth). Referencecurrents, ~j re f for the canonical block are circular, centered on the source location,and decay with distance. For the conductive block, additional vortex currents areinduced. These reference currents are compared with the polarization currents.Fig. 2.8 shows the plan and section view maps of the polarization currents at 0.86ms. Comparisons of Figs 2.7 and 2.8 clearly show that polarization currents forboth canonical and conductive blocks are oppositely aligned with respect to theirreference current. This was the hypothesized outcome. Fig. 2.9 shows that thedirection of ~jpol at 6.7 ms is similar to those at 0.86 ms. Thus both for the canonicaland conductive blocks, the direction of ~jpol after 0.86 ms is constant in time.Of particular interest is the difference in character of the ~jpol for the canoni-cal and conductive bodies. For the canonical body the ~jpol look like anomalousgalvanic currents that would be expected from an DC-IP survey. The resultantmagnetic fields will be similar to the magnetic fields obtained from an electricdipole. For the conductive case however, the currents are circular and they reflectthe vortex nature of the induced currents. The resultant magnetic fields are thoseassociated with a magnetic dipole. The ~jpol inside a body are therefore compli-cated by the fact that they are a mixture of galvanic (charge buildup) and inductive53Figure 2.7: Maps of reference currents: (a) canonical and (b) conductivemodels. Left and right panels show plan and section views at -125 m-depth and 0 m-easting, respectively. A source is located at (-200 m, 0m, 30 m). Black arrows and colored background respectively indicatethe direction and amplitude of the current. Black solid lines outline theboundary of the chargeable body.processes. Our choice of ~j re f effectively incorporates this complexity.54Figure 2.8: Maps of polarization currents: (a) canonical and (b) conductivemodels at 0.86 ms. Left and right panels show plan and section views at-125 m-depth and 0 m-easting, respectively. A source is located at (-200m, 0 m, 30 m). Black arrows and shaded values respectively indicate thedirection and amplitude of the current. Black solid lines outlin boundaryof the surface or the chargeable body.55Figure 2.9: Maps of polarization currents: (a) canonical and (b) conductivemodels at 6.7 ms. Left and right panels show plan and section viewsat -125 m-depth and 0 m-easting, respectively. A source is located at(-200 m, 0 m, 30 m). Black arrows and shaded values indicate the di-rection and amplitude of the current, respectively. Black solid outlinesboundary of the surface or the chargeable body.562.3.3 IP CurrentsThe IP currents, as provided in eq. (2.13), are given as~jIP = s•~e IP+~jpol. (2.40)In most analyses, (e.g. Smith et al. [1988]), the term s•~e IP is neglected. I haveincluded this term but with an approximation that~e IP ⇡—f IP (eq. 2.27). Here Iinvestigate these approximations, and under what circumstances they hold.The electric field,~eIP, can be evaluated with forward modelling. It can be bro-ken into galvanic galvanic and inductive parts using the Helmholtz decomposition.So ~eIP = ~—f IP ~aIP. Taking the divergence and making use of the fact that— ·~aIP = 0 (vector potential is divergence free). I obtain— ·~eIP =— ·~—f IP. (2.41)Then ~aIP = ~eIP+~—f IP. The effects from the scalar potential are included, butin my approximate formulation (justified below) I have neglected any contributionfrom the vector potential. I look at the contributions of each of these terms forthe three cases of canonical, conductive and resistive bodies. Fig. 2.10 respec-tively shows plan view maps of ~jpol , s•~—f IP, and s•~aIPfor (a) canonical, (b)conductive, and (c) resistive models at 0.86 ms.Inside the body, the ~jpol has the greatest strength and the strength of thesecurrents is largest in the conductive body and smallest in the resistive body. Inall cases, ~jpol is the largest contribution to ~jIP. The second column in Fig. 2.10 isrelated to the scalar potential for the electric field, or effectively to the galvanic cur-rents, s•~—f IP. These exist both inside and outside the chargeable body. Again,these are largest for the conductive body. Note that inside the body, s•~—f IP hasa direction that is opposite to the ~jpol . The third column is associated with the vec-tor potential and is associated with vortex currents: s•~aIP. The effects of thesecurrents have not been included in our linearized approximations. These currentsare quite small for the canonical and resistive models but their amplitude starts tobe comparable to the galvanic portion for the conductive model.I evaluate~jIP and its components at two locations in the body for the conductive57Table 2.1: Amplitudes of decomposed IP currents at two marked points(white stars) shown in Fig. 2.10(b). Units in A/m2Division |~jIP| |~jpol| |s•~—f IP| |s•~aIP|Left 1.5⇥1010 2.5⇥1010 7.6⇥1011 1.9⇥1012Right 5.4⇥1011 1.2⇥1010 3.5⇥1011 3.3⇥1011model. These are denoted by white stars in the figures. For both locations, ~jpol hasthe greatest strength ands•~aIP is smaller thans•~—f IP. The ~jIP is smaller than~jpol mostly because thes•~—f IP is in the opposite direction compared to the ~jpol .The results are tabulated in Table 2.1.The above figures provide insight about the three contributions to ~jIP but ofultimate interest is the effect of these currents on the measured data. Therefore theBiot-Savart law is applied to each current. It suffices to work with the conductivecase. Fig. 2.11 shows dIP computed from the polarization current (stars), galvanic(rectangles) and inductive portions (circles) of the IP current. Here solid and emptymarkers show negative and positive signs, respectively. The polarization currenthas the major contribution to the IP response although it is larger than the truevalue. This overshoot is primarily negated by the galvanic portion of IP responsesand further reduced because of the vortex currents. I notice that the contribution ofthe galvanic currents is generally larger than those due to the vortex currents exceptnear 0.4 ms. At 6.7 ms, the amplitude of the IP response due to the polarizationcurrent is about 130 percent of the true polarization current, while the galvanicportion is 30 percent. These results show that the assumption by Smith et al. [1988]is reasonable, but incorporation of the galvanic portion to the IP datum is significantat later times. The inductive portion of the IP responses is small compared tothe galvanic portion except for the time before 0.2 ms, and hence ignoring this isjustified for the three cases we examined.58Figure 2.10: Decomposition of the IP currents as ~jpol (left panel), s•~—f IP(middle panel), and s•~aIP (right panel) at 0.86 ms. Plan view mapsof the currents at -125 m-depth are shown: (a) canonical, (b) conduc-tive, and (c) resistive cases.59Figure 2.11: Comparisons of contributions of ~jpol , s•~—f IP, and s•~aIPto the observed IP magnetic field. Solid line indicates true bIPz re-sponses. Stars, rectangles, and circles correspondingly indicate each IPresponse generated by applying the Biot-Savart law to ~jpol ,s•~—f IP,and s•~aIP. Empty and solid markers represent positive and negativevalues, respectively.60Figure 2.12: Interpolated maps of (a) true and (b) approximate IP currents at0.86 ms. Left and right columns respectively show plan and sectionview maps at -125 m-depth and 0 m-easting.61Figure 2.13: Interpolated maps of (a) true and (b) approximate IP currentsat 6.7 ms. Left and right columns respectively show plan and sectionview maps at -125 m-depth and 0 m-easting.622.3.4 Validations of LinearizationRecall that the linearization of the IP function shown in eqs. ( 2.38) and ( 2.39)basically requires obtaining approximate IP currents and then applying the Biot-Savart law. I validate this linearization procedure by numerically comparing boththe approximate currents and dIP with true ones. First, for the comparison of IPcurrents, only the conductive model is considered since this is the most challengingcase for linearization. Fig. 2.12 compares the true and approximate IP currentsat 0.86 ms. The approximate IP currents match well, both in direction and ampli-tude, with the true IP currents both inside and outside the body. As shown in Fig.2.13 the agreement improves as time increases (see the directions of the true andapproximate IP currents at (0,0,-350) on the right panels of Figs 2.12 and 2.13).Second, the validity of computing linear dIP is tested. True dIP is obtained by sub-tracting the dF from d. Approximate dIP is obtained by applying the Biot-Savartlaw to approximate IP currents shown in Fig. 2.13(b). In Fig. 2.14, comparisonof true (blue solid line) and approximate (blue empty circles) dIP are presented.Approximate dIP show lower amplitude after 0.2 ms and differ by 33 percent at theextreme. The difference decreases with increasing time. Overall the two curvesare in reasonable agreement, thus validating our linearized forward modeling (eq.2.38). Also, the discrete Biot-Savart operator is tested by applying it to true IPcurrents and evaluating dIP. As shown in Fig. 2.14 this dIP (blue stars) is almostcoincident with the true one except before 0.01 ms.In addition, the same analysis of comparing true and approximate dIP data wascarried out for the canonical and resistive models, and presented in Fig. 2.14. Thetrue and approximate dIP for both cases show good agreements. Note however, thatdespite the fact that the linear IP function reasonably explains dIP for the resistivecase, the amplitude of dIP is very small compared to EM signal, and this IP signallikely cannot be identified in practice.63Figure 2.14: Comparison of true and approximate IP responses (bIPz ). Black,blue, and red colors respectively indicate canonical, conductive, andresistive cases. Solid lines indicate true bIPz , computed by subtractingthe fundamental response from the observation. The stars are the ap-plication of Biot-Savart to true IP current and generate bIPz BS. Emptycircles show our approximate bIPz approx response.642.4 Effective Pseudo-chargeabilityThe IP datum is linearized as a function of the pseudo-chargeability, eq. ( 2.38),and this was developed for a single source. The pseudo-chargeability was definedas a convolution between the impulse pseudo-chargeability, h˜ I and the time historyof electric field, wˆ, as shown in eq. ( 2.23). Considering a TEM survey includ-ing multiple sources, wˆ is different for each source resulting in different pseudo-chargeabilities. This is problematic especially when IP inversion is considered,since only a single pseudo-chargeability model is desired to be restored from thedIP data. This important issue should be handled, and it requires combining thepseudo-chargeabilites that arise from individual sources into a source-independenteffective pseudo-chargeability. Hence, a desired result is that a dIP datum for anysource takes the form: 266664dIP1 (t)dIP2 (t)...dIPnTx(t)377775=266664J1J2...JnTx377775hh˜(t)i, (2.42)where dIPk (t) and Jk indicates the IP datum and sensitivity matrix at k-th source.Here h˜(t) stands for an effective pseudo-chargeability, which represents pseudo-chargeability from all sources. Hence, for a given effective pseudo-chargeability Ican compute IP responses at all sources.How to obtain this effective pseudo-chargeability will differ for different typesof TEM surveys. I consider three types of TEM surveys:• Large ground loop survey: a single loop source with multiple receivers• Airborne survey: a number of a source and receiver loop pairs• Grounded source survey: multiple current (source) and potential (receiver)electrodes connected by current wiresThe first case poses no problem, since the survey includes a single loop sourceresulting in only a single pseudo-chargeability. The second case is the most prob-lematic since it includes many sources. This issue is addressed in Appendix A.65This requires computing and combining the individual time histories of the electricfields due to each source into an effective time history. The survey geometry ofthird case is the same as a DC-IP survey, but the effects of EM induction inducedfrom current wire paths is not ignored, but considered. Similar to the DC-IP sur-vey, however, steady-state can be achieved for polarization charge buildup at thelate on-time, then the polarization simply decays in the off-time (see red dIP curvein Fig. 1.13). That is, for the grounded sources, a steady-state electric field canbe established, and hence assuming that wˆ is equivalent to the input current wave-form is reasonable. This makes pseudo-chargeabilities of all sources in the surveyequivalent. For the above three surveys types, IP data which are linearly related toeffective pseudo-chargeability (eq. 2.42) has been obtained. Now for a given h˜ ,dIP can simply be computed as a matrix-vector product when the sensitivity matrix,J is formed.662.5 Positivity on Time History of Electric FieldThe time history of the electric field, wre f , was defined in eq. ( 2.22), and there Iignored negative values of wre f by using a positive projection, P0[·], as shown in eq.( 2.20). Here, I discuss why negative portions of wre f can be ignored and discusspotential errors on the assumption. The first motivation for using the positive pro-jection was to enforce the sign of the pseudo-chargeability, h˜ , to be positive. Intrin-sic chargeability is positive and, as in the case of grounded IP, pseudo-chargeabilityis also a positive quantity. For the inverse problem, forcing positivity on h˜ is animportant constraint. I show this in Chapter 3.Because of the complexities of earth structure and the inductive source field, Irecognize that there can be situations where wre f is negative. For instance, considera vertical conductor offset from a transmitter loop as shown in Fig. 2.15. After thetransmitter current is turned off, the currents induced due to the half-space (smokering) result in a current that is into the page at all depths, but the vortex currentinduced in the conductor is positive (into the page) at the top and negative (out ofthe page) at the bottom as shown in Fig. 2.15. So, at the top (yellow voxel) thehalf-space channeled current and the vortex current are in the same direction, butat the bottom (green voxel) the two currents are in opposite directions. Hence, thebottom voxel has a potential to have negative wre f values based upon eq. ( 2.22).In order to address this potential problem, I first carry out a TEM simulationusing SIMPEG-EM package [Heagy et al., 2017]. A thin vertical conductor (0.005S/m) is embedded in a homogeneous half-space earth (0.001 S/m) as shown inFig. 2.16. A transmitter loop with radius 13 m is located 30 m above the surface,and 150 m away from the conductor (marked as a green solid dot in Fig. 2.16);step-off waveform is used. A receiver which measures the vertical magnetic fieldis co-located with the transmitter. Fig. 2.17 shows the current density at threedifferent times after the transmitter current is turned off. At 0.01 ms, the half-spacechanneled currents and the vortex currents are both evident. At the top part ofthe conductor (see yellow cross mark), the two currents are in the same direction,which increases the net current. However, at the bottom of the conductor (see yel-low cross mark), the two currents are in opposition so the net current is decreased.As time passes, the vortex currents decay away (Fig. 2.17b) and eventually only67the channeled currents (Fig. 2.17c) remain. The direction of the currents in theconductor changes with time, and hence wre f can have negative value at some vox-els.Now, suppose the vertical conductor is chargeable and its Cole-Cole parame-ters are h=0.5, t=104 s, and c=1. The question arises regarding how much impactthose negative wre f have on the IP data. To answer the question I carry out follow-ing procedures:• Step 1: Compute the reference electric field, e re f :~ere f =~e⌦~t re f .Here t re f is the time when the electic field,~e, has a peak amplitude in time.• Step 2: Evaluate wre f using eq. ( 2.22) and P0[wre f ].wre f (t) =~e(t) ·~e re f~e re f ·~e re f .• Step 3: Compute h˜ (eq. 2.23) with wre f and P0[wre f ]:h˜ = h˜ I⌦wre fandh˜0 = h˜ I⌦P0[wre f ],respectively.• Step 4: Calculate IP data, dIP, with h˜ by using the linear IP function:dIP = Jh˜ ,and dIP0 with h˜0:dIP0 = Jh˜0.As long as dIP and dIP0 are close enough, my assumption is valid in spite ofnegative wre f values.68Fig. 2.18 shows the peak current computed from the time history of the cur-rents shown in Fig. 2.17, and it effectively captures both the channeled currentsand the vortex currents. The time history of the electric field (wre f ) for two voxels ispresented in Fig. 2.19(a); the time is in log scale. Here I also provided P0[wre f ] forcomparison (cross marks). At the top voxel (yellow line), wre f is always positive,so there is no problem. However, at the bottom voxel (green line), wre f changesfrom positive to negative around 0.15 ms. These negative values are clearly shownin Fig. 2.20(b) where they are plotted with a linear scale in time. The resultantpseudo-chargeability (eq. 2.23) shows negative values after 0.4 ms.IP data (dIP) at the receiver location (evaluated with dIP = Jh˜) is presented inFig. 2.21(a). The same forward modelling is carried out using h˜0 and the resultantdIP0 is provided (red cross marks). Note that dIP and dIP0 are almost coincident.Thus even though some pixels will be modelled incorrectly, the fact that an IP da-tum is the cumulative response due to all of the pixels, most of which are correctlymodelled, means that our assumption is likely quite robust. An additional reasonto give confidence in my assumption is that an AEM survey uses multiple trans-mitters. Each transmitter will excite the earth differently. The pixel at the bottomof the plate, while having a negative wre f when the transmitter was displaced fromthe plate, may have a positive wre f when the transmitter moves location. The ef-fective pseudo-chargeability derived from all transmitters used in an AEM survey,basically averages h˜s from different transmitters and this may ameliorate negativeeffects caused by specific coupling geometries (See Appendix 2.4).Lastly, I comment on the effects of neglecting negative wre f in the inverse prob-lem. Our basic equation is dIP = Jh˜ . Negative values of h˜ will reduce the value ofdIP compared to that obtained with our assumption (although the plot show barelyno change). After the inverse problem is solved the recovered h˜ will therefore havevalues lower than they should; that is the recovered pseudo-chargeability will beunderestimated.69Figure 2.15: Conceptual diagram of EM induction process from a loop trans-mitter when a vertical conductor is embedded in a homogeneous half-space. After the transmitter current is turned-off from the loop, cur-rents are induced in the subsurface and diffuse away like a smoke ring.This smoke ring currents will be channeled into the vertical conductor(dashed lines), and there will be vortex currents rotating (red dashedcircle) in the conductor due to the time-varying magnetic field. Inthe conductor, at the top voxel (yellow) the half-space channeled cur-rent and the vortex current are in the same direction (in to the page),whereas at the bottom voxel they are in the opposite direction.70Figure 2.16: Plan and section views of a 3D conductivity model is shown inright and left panels, respectively. A vertical conductor (0.005 S/m)is embedded in a homogeneous half-space (0.001 S/m). Green solidcircle indicates the horizontal location of the transmitter loop having13 m-radius, and located 30 m above the surface.71Figure 2.17: The current densities at three different times due to a 13 m-radius loop located 30 m above the surface (marked as green solidcircle) with a conductivity model having a vertical conductor shown inFig. 2.16. Plan and section views of the current density at three differ-ent times after the current is off: (a) 0.01 ms, (b) 0.07 ms, and (c) 1.17ms. The time history of the electric field at the top and bottom voxelsin the conductor (yellow and green crosses, respectively) are presentedin Fig. 2.19.72Figure 2.18: The peak current (or reference current) selected from the timehistory of the current density shown in Fig. 2.17. Plan and sectionsviews are shown in right and left panels, respectively. The time his-tory of the electric field at the top and bottom voxels in the conductor(yellow and green crosses, respectively) are presented in Fig. 2.19.73Figure 2.19: The time history of the electric field, wre f , at the two voxelsshown in Figs. 2.17 and 2.18. Yellow and green distinguish wre f forthe top and bottom voxels, respectively. (a) wre f at the both voxels insemi-log scale (log(t)wre f ). (b) wre f at the top voxel only showinglater times after 0.1 ms with linear scale in time. The negative values ofwre f exist at the very late-time for the bottom voxel (green) as shownin (b).74Figure 2.20: (a) Intrinsic pseudo-chargeability, h˜ I as a function of time. UsedCole-Cole parameters are denoted in (a). (b) Pseudo-chargeability, h˜ ,calculated by convolving h˜ I and wre f at the bottom voxel shown inFig. 2.19. Solid line and cross marks distinguish h˜ computed withand without positive projection on wre f , P0[wre f ].75Figure 2.21: Computed IP data, dIP, using linear IP function. IP data arecomputed with and without positive projection on wre f and denoted asdIP and dIP0 , respectively; red dashed line and red crosses are distin-guishing dIP and dIP0 , respectively.762.6 ConclusionsThe linearization of dIP with respect to a pseudo-chargeability is derived for theISIP data. Pseudo-chargeability is defined as the ratio of the polarization currentto a reference current. Unlike the DC-IP case, the electric fields for an inductivesource do not achieve steady-state and hence neither do the polarization currents.To address this important difference I evaluate the fundamental electric field at eachlocation in the earth and generate a reference electric field that has the direction andmagnitude of the field at the time when the fundamental field reaches its maximumvalue. The pseudo-chargeability at a point in the earth thus depends upon the po-larization parameters (e.g. Cole-Cole: h , t , and c), the reference electric field,and the time history of the electric field. The situation becomes more complicatedwhen data from many sources such as an AEM survey are to be inverted simultane-ously because the time history of the electric field at a point in the earth is differentfor each source. This is handled by defining an effective pseudo-chargeability andan associated reference electric field that accommodates, in a least squares fashion,the effects of all sources acting on a single cell.To have confidence in when, and under what circumstances, the approxima-tions are sufficiently valid, a number of rigorous tests are carried out as a syntheticexample of an AEM survey. First, three test models consisting of a chargeableblock in a non-chargeable halfspace are introduced. The block can be conductive,canonical, or resistive with respect to the halfspace conductivity. The evaluationsshow that: (a) the choice of reference electric field and its time history producesa good estimate of the ~jpol; (b) the ~jIP is dominated by the ~jpol , which is an as-sumption that has been made. However, the galvanic and vortex currents arisingfrom the scalar and vector potentials in the Helmholtz decomposition of ~e IP canbe significant in some circumstances. The galvanic currents are the second mostimportant contribution to the IP currents and, in the body, they have a direction thatopposes the direction of the polarization currents. For the linear IP function, I haveincluded the galvanic currents and neglected the vortex currents which are almostalways smaller than the galvanic currents; (c) the dIP data can be evaluated usingthe Biot-Savart law which provides quite accurate results; (d) with the approximate~jIP, the predicted responses are in reasonably good agreement with true values al-77though they are underestimated at early times for the highly conductive example.These results lead us to infer that our linearized formulation dIP(t) = Jh˜(t) is a vi-able representation for the forward modelling at late times when the IP effects aresubstantial compared to the EM effects. (e) For the multi-source case I derivedan effective pseudo-chargeability which is a linear combination of the pseudo-chargeability of each source. These were forward modelled with the linearizedformulation and compared to the true responses. The values were underestimatedby 33 percent at the extreme for the conductive model but were almost identical forthe canonical and resistive models.Although the numerical experiments shown here focuses on the ISIP method,the linearization itself is general enough to handle dIP from various TEM surveys.This is a valuable contribution in two main perspectives:• The linearized formulation will be an effective forward simulator to computeIP effects in various TEM surveys.• The linearized formulation will be used as a forward engine to invert dIP torecover the distribution of h˜ .78Chapter 33D IP Inversion3.1 IntroductionFrom the previous chapter, dIP is linearized with respect to the pseudo-chargeability:dIP(t) = Jh˜(t) for various TEM surveys. Using this linear equation, I develop a 3DIP inversion algorithm, which separately inverts each time channel of dIP. Solvingthis linear inverse problem is common in applied geophysics so only an essentialsummary is provided in Section 3.2. After h˜ at multiple time channels are ob-tained, a small inverse problem is solved to extract intrinsic polarization parame-ters for each cell. The 3D IP inversion algorithm is tested with synthetic examplesusing an AEM survey geometry with a simple conductive and chargeable blockembedded in a halfspace (Section 3.3).3.2 3D IP Inversion with a Linear EquationThe linear inverse problem to recover chargeability is straightforward and is de-scribed in Oldenburg and Li [1994]. I rewrite eq. (2.38) asdpred = Jm, (3.1)where J is the sensitivity matrix of the linear problem, which corresponds to the Jshown in eq. (2.38). Here, dpred represents IP responses at a single time channel,mdenotes model parameters, which can be either h˜ or ∂ h˜∂ t . The important positivity79constraint results because the intrinsic chargeability h is restricted to the range[0,1).The solution to the inverse problem is the modelm that solves the optimizationproblemminimize f = fd(m)+bfm(m) (3.2)s.t.m 0,where fd is a measure of data misfit, fm is a model objective function and theb is a regularization or trade-off parameter.I use the sum of the squares to measure data misfitfd = kWd(Amdobs|)k22 =NÂj=1(dpredj dobsje j)2, (3.3)where N is the number of the observed data and Wd is a diagonal data weightingmatrix which contains the reciprocal of the estimated uncertainty of each datum(e j) on the main diagonal, dobs is a vector containing the observed data, dpred isa vector containing calculated data from a linear equation given in eq. (3.1). Themodel objective function, fm, is a measure of the amount of structure in the modeland upon minimization this will generate a smooth model which is close to a ref-erence model, mre f . I define fm asfm = Âi=s,x,y,zaikWiW(mmre f )k22, (3.4)where W is a model weighting matrix, which will be defined below, Ws is a diag-onal matrix containing volumetric information of prisms, andWx,Wy andWz arediscrete approximations of the first derivative operator in x, y and z directions, re-spectively. The a’s are weighting parameters that balance the relative importanceof producing small or smooth models [Tikhonov and Arsenin, 1977].Depending on the geometry of TEM surveys, additional weighting could benecessary. For instance in an AEM survey, because only a single datum for eachsource is available, intrinsic depth resolution is lacking. A similar lack of depthresolution occurs when a gradient array (a single grounded source and multiple re-80ceivers) is used for a DC-IP survey. These are also the same circumstances encoun-tered when inverting magnetic data [Li and Oldenburg, 1996]. Correspondingly Iapply a depth weighting through the model weighting matrix (W):W= diag(z z0)1.5, (3.5)where z and z0 are discretized depth locations and reference depth in the 3D do-main. Here the expnent 1.5 arises from the geometric decay of magnetic field⇠ 1r3 .Although I use the linear form of dIP data (eq. 2.38), the inverse problem isnonlinear because of imposed positivity on m. An initial model, m0 is requiredto start the inversion. This constrained optimization problem is solved by usinga projected Gauss-Newton (GN) method [Kelley, 1999]. The trade-off parameter,b , is determined using a cooling technique where b is progressively reduced fromsome high value. The inversion is stopped when the tolerance is reached; the mostobvious stopping criteria is a target misfit, fd and it is set to the number of data (cf.Oldenburg and Li [2005]). For the implementation of the IP inversion algorithm,I use an open source python package for simulation and gradient-based parameterestimation in geophysics called SIMPEG [Cockett et al., 2015].3.3 Numerical ExperimentsUsing the 3D IP inversion algorithm, which uses the linear IP equation as a forwardengine (eq. 2.38), I now invert dIP, and recover pseudo-chargeability. As a test, thesame model and the AEM survey setup used in the previous section are used exceptthat only the conductive and chargeable block model is used here. Conductivity ofthe block is s2=0.1 S/m, and the halfspace conductivity is s1=0.001 S/m. Thesurvey, consisting of 11 soundings along each of 11 lines, is shown in Fig. 2.4(a).Data are from a coincident-loop system and the flight height is 30 m above thesurface; the radius of the loop is 10 m. A step-off source waveform is used and therange of the observed time channels is 0.01-60 ms. The observed responses are thevertical component of~b.When computing the sensitivity matrix, a 3D distribution of s•(x,y,z) is re-quired, and an estimate of this distribution can be achieved in the first step of theTEM-IP workflow: inverting TEM data to recover conductivity. This estimated81conductivity, sest will always include some errors, and hence this will have someimpact on the sensitivity computation. To start, I use the true s• to calculate thesensitivity matrix, but the issue of incorrect conductivity is addressed in Section3.3.1. Using the true s•, the sensitivity is computed, dIP data at successive timechannels are inverted, and 3D pseudo-changeability at multiple times are recov-ered. The 3D IP inversion is based upon Oldenburg and Li [1994] and Li andOldenburg [2000], and it requires some choices for inversion parameters.For data uncertainties, I use one percent of the maximum amplitude of theobserved data (0.01max(|dobs|)). Coefficients for smallness and smoothness areset to as = 105 and ax = ay = az = 1, respectively. The reference model is zero,which means the pseudo-chargeability of every cell is zero, and I applied a depthweighting. The need for a depth weighting arises because the sensitivity functionJ is primarily controlled by a 1/r3 decay associated with the Biot-Savart kernelsand IP data map at a single time is inverted to recover a 3D distribution of thepseudo-chargeability. Thus an AEM data set is like a magnetic data set where it iswell established that a depth weighting is required to image objects at depth. Thefollowing example illustrates this.I first generate IP responses at a single time using the linear function and spec-ifying that the pseudo-chargeability is unity inside the body and zero outside, asshown in Fig. 3.1(a). Fig. 3.1(b) shows the recovered pseudo-chargeability withoutdepth weighting. The recovered anomalous pseudo-chargeability is concentratednear the surface and the magnitude of the pseudo-chargeability is underestimated;it is ⇠0.2 rather than unity. By using the depth weighting shown in eq. (3.5), theIP body is imaged closer to its true depth (Fig. 3.1b). Also, the magnitude of therecovered pseudo-chargeability (⇠0.6) is closer to the true value than the resultwithout depth weighting. Based on this analysis, I use the same depth weightingfor following examples.3.3.1 Incorrect ConductivityThe 3D distribution of s•(x,y,z) plays a central role in the TEM-IP inversion work-flow. It is used in the EM-decoupling process to obtain the IP data: dIP = dF [s•]and it is also needed to compute the linearized sensitivity, J[s•], for the IP inver-82Figure 3.1: Effect of depth weighting in 3D IP inversion. (a) True pseudo-chargeability model on vertical section at 0 m-northing. Recoveredpseudo-chargeability models (b) without depth weighting and (c) withdepth weighting.sion. Therefore, estimating the 3D distribution of s• is an essential step, andinverting early TEM signals having minor IP effects can be an effective option. Inthis section I do not focus on estimating s•. I do however appreciated that this willnever be known exactly. This issue is addressed in more detail with both syntheticand field AEM examples in Chapter 4, but here some consequences of having anincorrect sest are addressed, in particular with regard to the sensitivity function. Ireturn to the conductive block in a halfspace and evaluate the dIP data when thehalfspace conductivity is the true value (s1 = 103 S/m) as well as a factor of twotoo large (2⇥103 S/m) and a factor of two too small (5⇥104 S/m). Here theconductivity of the chargeable block is fixed to s2 = 0.1 S/m. The data along asurvey line are plotted in Fig. 3.2.I invert these three IP responses, and provide sections of the recovered pseudo-chargeability at 0 m-northing. Fig. 3.3(a), (b) and (c) correspondingly show therecovered pseudo-chargeability when the conductivity is: the true value, too high,or too low. With the correct conductivity the geometry of the IP body is reasonablyrecovered. When the halfspace conductivity is too high, the dIP have a negativebias that results in larger pseudo-chargeabilities and positive-valued artifacts nearthe IP body (Fig. 3.3b). When the halfspace conductivity is too small, the IPdata have a positive bias and this produces negative-valued pseudo-chargeabilityartifacts either side of the IP body (Fig. 3.3c). White dotted contours shown in Fig.3.3(c) show zero-crossing lines, which delineate those negative-valued artifacts.83However, based on the definition of the pseudo-chargeability shown in eq. (2.23),the sign of the pseudo-chargeability should be positive. By incorporating positivityas a constraint in the inversion, and re-inverting the IP data that have a positive bias,I obtain the result in Fig. 3.3(d). Considering there are no negative values in therecovered pseudo-chargeability, this is a better result than Fig. 3.3(c). This resultsshows that the positive constraints prevent fitting positive residual fields. I shalluse this positivity constraint for following 3D IP inversion examples.A distribution of s•(x,y,z) is required when computing the sensitivity func-tion. An incorrect conductivity will affect the sensitivity function. In order totest this, the sensitivity matrix is computed using a halfspace conductivity model(s1 = s2 = 0.001 S/m). Fig. 3.4 compares the recovered pseudo-chargeability fromthe 3D IP inversion of the IP datum at 0.86 ms with the true and incorrect sensi-tivity function using the halfspace conductivity. There is not a large differencebetween the two inversions and this suggests that an approximate conductivity forthe chargeable and conductive body may still provide sensitivities that are adequatefor inversion. This parallels results from inverting DC-IP data where even an ap-proximate conductivity can still yield good results when inverting the data. Thusthere is robustness in our sensitivity function with respect to an incorrect conduc-tivity. These results suggest that even if one cannot generate a good estimate ofs•, a halfspace conductivity might produce an adequate sensitivity function, andhence an inversion can provide some indication of a chargeable body.84Figure 3.2: IP responses on a profile line at 0 m-northing. IP responses arecomputed from perturbed s• models. Halfspace conductivity (s1) isperturbed: 2⇥s1 and 0.5⇥s1, which respectively resulting in overesti-mated (dotted line) and underestimated (dashed line) IP responses. Thesolid line shows the true IP response.85Figure 3.3: Recovered pseudo-chargeability sections from 3D IP inversionsat 0 m-northing. (a) dIP with true s1. (b) dIP with 2⇥s1. (c) dIPwith 0.5⇥s1. (d) dIP with 0.5⇥s1 and the positivity constraint on thepseudo-chargeability. White dashed lines contour zero-crossing lines.86Figure 3.4: Recovered pseudo-chargeability sections from the 3D IP inver-sions at 0 m-northing. (a) True and (b) incorrect s• is used to computesensitivity function. For the incorrect sensitivity I used a halfspace con-ductivity s1.873.3.2 Extracting Intrinsic IP ParametersBy applying the 3D inversion to each time channel of dIP data separately, 3D dis-tributions of pseudo-chargeability at multiple times can be recovered. The pseudo-chargeability at each time carries different information about the state of polariza-tion and these can be used to recover information about intrinsic IP parameters suchas Cole-Cole parameters. Diverse time-dependent conductivity models such as theCole-Cole model and the Stretched-exponential can be used for this interpretation[Kohlrausch, 1854, Cole and Cole, 1941, Pelton et al., 1978, Tarasov and Titov,2013].The Cole-Cole model with c= 1 is used here. Pseudo-chargeability at a singlepixel is parameterized in terms of chargeability and time constant, and recoveredby solving a small inverse problem (see Appendix B for further details). In pre-vious works about this task for the DC-IP problem [Yuval and Oldenburg, 1997,Ho¨rdt et al., 2006], the convolution shown in eq. (B.1) was not explicitly men-tioned because wˆ(t) is a step-off or step-on function and it does not change fordifferent cells and transmitters. This allowed an explicit equation for a step-off orstep-on response of the pseudo-chargeability to be derived. However, in our work,convolution plays a fundamental role and needs to be explicitly addressed whenextracting intrinsic IP parameters. Also, the details regarding how I defined theeffective pseudo-chargeability (eq. A.8) need to be included. Except for this addi-tional complexity related to the convolution, our approach parallels that of Yuvaland Oldenburg [1997] and Ho¨rdt et al. [2006].As an example, I use the conductive and chargeable block presented in the pre-vious section and invert 14 time channels of data ranging from 1-10 ms. The EMdata are forward modelled using EMTDIP code and the true s• model is used toevaluate the IP datum and compute the sensitivity function. The recovered pseudo-chargeability from one of the 14 inversions is shown in Fig. 3.4(a). In that pseudo-chargeability model, I select cells that have a pseudo-chargeability value greaterthan 0.001, and then use the data from all 14 inversions to carry out the nonlin-ear inversion to estimate the time constant, t , and chargeability, h , for each cellseparately. The forward modelling for this inversion is shown in eq. (A.8), whichrequires we(t) (eq. A.9). The we(t) for a pixel in the block is shown in Fig. A.2.88Figure 3.5: Section views of recovered: (a) time constant and (b) chargeabil-ity. Any region where the pseudo-chargeability shown in Fig. 3.4a issmaller that 0.001 is ignored in this analysis, and blanked.Fig. 3.5(a) and (b) correspondingly show the estimated time constants andchargeability as section maps. The estimated time constants show good agreementwith the true value t= 0.005 s. There is less agreement about chargeability forwhich the true value is h = 0.2. Recovered values range from about 0.04-0.2 somost values are significantly underestimated. In Fig. 3.6, I also provide time de-cays of the observed and predicted pseudo-chargeabilities at a single pixel markedas a black rectangle in Fig. 3.5. The estimated time constant, test , and chargeabil-ity, hest , for this pixel are 0.0046 and 0.09, respectively. These results imply thereis greater stability in recovering the time constant than in recovering chargeabilitywith our approach. Again, similar experiments were carried out for the canoni-cal and resistive bodies and the conclusions were also that the time constant wasadequately recovered with better fidelity than was the chargeability.3.4 ConclusionsThe 3D inversion algorithm to invert time-domain IP data using the linearized for-mulation is developed and tested for the synthetic AEM survey example. IP datafor the AEM survey have only one receiver for each transmitter and a data mapat a single time channel is essentially a potential field. The data do not have in-trinsic depth resolving power and hence, as in magnetics or gravity inversions, I89Figure 3.6: Comparisons of the observed pseudo-chargeability at a singlepixel in a chargeable body from the 14 inversions at 14 time channels,and the predicted pseudo-chargeability. The empty circles and solid linerespectively indicate predicted and observed pseudo-chargeability. Theestimated time constant and chargeability are respectively expressed astest and hest . The true values for t and h are respectively 0.005 s and0.2.attempt to counteract this by introducing a depth weighting. When this is done, the3D IP inversion recovers a reasonable geometric shape and location of the charge-able body but the amplitude is underestimated. For the inversion it is assumedthat a good estimate of s• is available. An incorrect s• has two effects in theinversion. Firstly it can generates errors in the dIP because the estimated funda-mental field, which is subtracted from the observations, is incorrect. To obtaininsight I looked at the effects when sest was too low or two high. This respectivelyyielded positive or negative residual fields in the IP response. A positivity con-straint on the pseudo-chargeability (similar to that used in DC-IP surveys) greatly90ameliorated the effects of the positive residuals. The other avenue by which anincorrect s• can affect the inversion is through the sensitivity matrix J. I showedthat, even with an approximate conductivity, I recovered important informationabout the chargeable body such as geometric shape and location. An inversionof the data at a particular time channel provides information about the effectivepseudo-chargeability for each pixel. Inversions carried out at multiple time chan-nels therefore generates a pseudo-chargeability as a function of time for each pixel.The pseudo-chargeability for pixels that had significant chargeability were subse-quently fit to a Cole-Cole model to estimate t and h by assuming c = 1. Theestimated t was close to the true value whereas h was underestimated and less ro-bust. This suggests that there is a possibility to extract intrinsic IP parameters fromthe recovered pseudo-chargeability from AEM surveys if the c values are knownor assummed91Chapter 4Airborne ExamplesMost airborne EM (AEM) surveys use wire loops for the source and receiver. Thispair of source and receiver loops is attached to either a helicopter or a plane, andcan cover a large area quickly. Fig. 4.1 illustrates the geometry of an exampleAEM system. In particular for time-domain systems, current in the source loopis abruptly turned off after the on-time. This changing current generates a time-varying magnetic field, which induces electrical fields in the earth. These electricfields in the earth diffuse away. Recalling that there will always be polarizationeffects whenever an electric field is applied to a chargeable medium, induced po-larization (IP) effects can occur due to a changing current induced by an airbornesystem, and IP responses could be measured in the air as well. These IP effects areoften observed to have negative transients for late-time channels. Weidelt [1982]showed that for a coincident-loop system, negative transients can be obtained onlyif the earth medium is chargeable. Historically, negatives have been measuredfrom a number of field AEM surveys [Smith and Klein, 1996, Kratzer and Macnae,2012, Kaminski and Viezzoli, 2017], and today it is common to measure negativesespecially because of the recent developments in instrumentation; specifically lowtransmitter and receiver heights, large inducing currents and low noise levels. Forinstance, Fig. 4.2(a) shows a profile line of data obtained from a VTEM survey(2007) over the Mt Milligan porphyry deposit. Around the easting of 433800 m,negatives are observed (red lines). The time decay at a sounding around that loca-tion are provided in Fig. 4.2(b), and this shows negative transients at time channels92(> 1 ms).To illustrate the physics behind negative transients, a simple forward simulationis carried out using EMTDIP code. Fig. 4.3 shows a cylindrical chargeable bodyembedded in a halfspace earth. The conductivity of the halfspace is 0.001 S/m,and the conductivity of the chargeable medium is 0.01 S/m. TEM simulation isconducted with a loop source with a step-off current waveform; the radius of theloop is 13 m, and it is located 30 m above the surface. Voltage data is measuredat the receiver loop when the current is turned off. I first present the electricalfield in the subsurface, then show the measured voltages as a function of time.Fig. 4.4 shows the x-z section of the electrical field in the y-direction at threedifferent times. At 0.1 ms, large electric fields are induced in a chargeable medium.Note that the electric field is rotating in the counter-clock-wise direction. As timeproceeds, these induced electric fields in the conductor diffuse away (0.5 ms), andat later time (7.9 ms), the electric field is reversed (clock-wise direction); this isdue to the polarization current. Fig. 4.5(a) shows the observation at the receiverloop, and negatives are apparent after 0.9 ms. Clearly, these negatives are dueto polarization currents flowing in the reverse direction. Another simulation isperformed without IP effects (dF = F [s•]), and dIP is obtained by subtracting dFfrom d. Time decays of all three data: d, dF and dIP are plotted in Fig. 4.5(b). Inthe early-time (< 0.03 ms) EM is dominant (d ' dF ), whereas in the late-time (>3 ms) IP is dominant (d ' dIP). Intermediate time exists when both EM and IP aresignificant. Thus, AEM data can include IP effects and an effective methodologyto handle these effects is necessary.There can be two perspectives on the IP effects in AEM data: (a) They arenoise, which disturbs conductivity signals and (b) IP effects are signals, and theyinclude chargeability information of the subsurface. Both perspectives are valid,and depend upon the goal of an AEM survey and the sources of IP. For miningapplications such as seeking kimberlites, IP effects due to permafrost can be con-sidered as geologic noise, which contaminates the conductivity signals from thekimberlites. On the other hand, if kimberlites are chargeable, then IP effects in thedata may be considered as signal. The main interest of this Chapter is in the latter,and hence I use the TEM-IP inversion workflow to recover distributed polarizationinformation. This workflow is based upon natural separation of EM and IP effects93in time (Fig. 4.5b), and includes three main steps:• Step 1: Conductivity inversion. Here I need to invert data that are not con-taminated with IP effects to recover the conductivity distribution, sest(x,y,z).Negatives are obvious, and those can be removed. However positive datacould include significant IP effects (see 0.1 ms in Fig. 4.5b), and hencedetermining the latest time at which EM is dominant is a challenge. Forthe work here I ignored all negative data, and also the first few positive databefore the sign reversal.• Step 2: EM-decoupling. In this step, EM effects are considered as the noise.The fundamental response is estimated by evaluating F [sest ], and subtractingit from the observation to yielddIP[sest ] = dF [sest ] (4.1)The obtained IP data, dIP[sest ] will always have some errors because sestis not the same as the true s•(x,y,z). EM-decoupling will be effective inthe intermediate time when both EM and IP effects are significant. For theEM-dominant time it will not be effective, since the IP is too small; EM-decoupling is not required for the late-time.• Step 3: IP inversion. Each time channel of the obtained dIP is separately in-verted, and pseudo-chargeability: h˜(x,y,z; t) is recovered. Pseudo-chargeabilityitself can provide some qualitative information about the distributed charge-ability. Further, intrinsic polarization information such as Cole-Cole param-eters, can be estimated by solving a small inverse problem using pseudo-chargeability of a single cell at multiple times.I now apply this TEM-IP inversion workflow to three examples. The first is asynthetic example (Section 4.1), and here I test the capability of the workflow tohandle airborne IP data. The second example shows an application of the workflowto field AEM data over the Mt Milligan porphyry deposit (Section 4.2), whichshowed clear negative transients. A 3D distribution of chargeability informationis recovered, and this reveals highly altered zones surrounding a resistive stock.94Figure 4.1: Geometry of the VTEM system (Geotech).The third example presents much stronger negatives obtained over a kimberlitecomplex at Tli Kwi Cho (TKC) region, and here I show how both conductivityand chargeability information obtained from AEM data can help distinguishingdifferent kimberlites in the ground.95Figure 4.2: Observed VTEM data over the Mt Milligan porphyry deposit,BC, Canada. (a) Profile line at Northing 6109000 m, which crossesHarris fault and Rainbow faults. (b) Time decay at Northing 6109000m and Easting 433800 m. Here black and red lines distinguish positiveand negative values.96Figure 4.3: Section of a chargeable cylinder embedded in a halfspace. Con-ductivity of the halfspace is s1=0.001 S/m. The conductive chargeablecylinder is parameterized with a Cole-Cole model and has higher con-ductivity than that of halfspace: s•=0.01 S/m. Other parameters aregiven in the diagram.Figure 4.4: Propagation of electric fields in time: 0.1 ms (left), 0.5 ms (mid-dle), and 7.9 ms (right). ey is plotted in the x-z plane, and fundamen-tal electric field is rotating in counter-clock-wise-direction. In the late-time, the reversed direction of ey is observed due to strong polarizationeffects. Boundaries of the air-earth interface and chargeable cylinderare shown as black lines.97Figure 4.5: Time decays of observed (d), fundamental (dF ), and IP (dIP) volt-ages over a chargeable cylinder. Solid and dashed lines distinguish pos-itive and negative values. (a) Only observed voltage. (b) All three volt-ages with distinction of early, intermediate, and late times. Vertical greylines indicate three times (0.1, 0.5, and 7.9 ms) when electric fields areplotted in Fig. 4.4.984.1 Synthetic Example4.1.1 SetupTo interrogate how the workflow could be applied to airborne IP data, a syntheticCole-Cole model is composed assuming c=1. The s•(x,y,z) model is shown inFig. 4.6. Conductivities of the resistive host rock and conductive overburden onthe west side are 103 and 102 S/m, respectively. Four IP blocks, which namedA1-A4 are embedded in the subsurface. IP blocks have different s• and h values asshown in Table 4.1. A3 and A4 are conductors, but A4 is buried at a greater depth.A1 has the same conductivity value as the overburden, and A2 is a resistor. Exceptfor the four IP blocks, h values are set to zero, and t is set to a constant value, 0.005s. To compute synthetic AEM data from this Cole-Cole model, the EMTDIP codeis used [Marchant et al., 2014]. The survey geometry includes 21 lines as shownin the top-left panel of Fig. 4.6 as black dots. We use a coincident-loop systemand the loop is located 30 m above the surface; the radius of the loop is 13 m. Astep-off current waveform is used and the range of measured time channels is 0.01-10 ms. The measured responses is voltage, which is proportional to the verticalcomponent of  ∂~b∂ t . To test the capability of the TEM-IP inversion workflow forairborne IP data, I generate synthetic airborne IP data using the above setup, andapply the workflow to them. The following section presents the synthetic airborneIP data.Figure 4.6: Sectional views of s•(x,y,z)model. Solid lines delineate bound-aries of four IP blocks.99Table 4.1: Cole-Cole parameters of four anomalous bodies (A1-A4).Division A1 A2 A3 A4s• (S/m) 0.01 0.001 0.1 0.1h 0.5 0.5 0.4 0.8t (s) 0.005 0.005 0.005 0.0051004.1.2 DataA forward simulation is performed using the setup described in Section 4.1.1, andhere I visualize this synthetic data in three different forms:• 2D map (at a single time, but all stations)• Profile line (at multiple times)• Time decay (at a station)Fig. 4.7 shows the 2D map of the observed voltage at 0.2 ms, 1.8 ms and 5.6ms. Here, negative values are presented as white on the color map. At 0.2 ms,EM induction is dominant. The boundary of the conductive overburden and resis-tive earth clearly show up around easting 250 m. Both conductive and resistiveanomalies appears near A3 and A2, respectively. The deep conductor, A4 does notshow a high anomaly at this time. At later time, the effect of polarization increases,and the first negatives are shown at 1.8 ms near A3; at 5.6 ms, stations around allfour blocks show negatives, hence IP is dominant. For visualizing airborne IP data,often a profile line of data at multiple times are plotted in double-log plot for bet-ter representation of negative IP signatures. Fig. 4.8 shows the observed data attwo profile lines, which cross the four IP blocks. Near A1 and A3, negatives areclearly shown at late-time channels (marked as red lines in Fig. 4.8 a). Early-timeresponses of A1 and A3 are different due to their different conductivity structures;A3 shows a high anomaly, whereas A1 is flat. This trend changes with time: bothA1 and A3 show low or negative anomaly data as time passes. Profile line dataaround A2 and A4 are shown in Fig. 4.8(b). Compared to A1 and A3, muchsmaller negatives are shown near A2 and A4, because A2 and A4 are correspond-ingly a resistor and a deep conductor; they generate only small polarization effects.Resistors do not have enough induction, and the deep conductor has weaker po-larization effects by geometric decay (1/r3). Effects of conductive overburden aremuch clearer on this line (⇠ 400 m easting). Fig. 4.9 shows four time decays overthe four A1-A4 blocks. All four decays show negative transients at later times, butthe sign reversals occur at different times due to different conductivity and charge-ability structures. In the following sections, I apply the workflow to these syntheticairborne IP data.101Figure 4.7: 2D map of the simulated data including both EM and IP effects(d). A1-A4 indicate corresponding anomalies due to four IP bodies.Black dotted lines indicate boundaries of the four IP bodies. Negativeresponses are shaded as white regions. Black dots are locations of Fig.4.9. Black dots are locations of Fig. 4.9.102Figure 4.8: Profile line of the observed data, d at Northing (a) 630 m and (b)-630 m. Black and red lines distinguish positive and negative data.103Figure 4.9: Time decays of the observed data, d at four sounding locationsclose to A1-A4. Here solid and dashed lines distinguish positive andnegative data. Error bars are shown from 0.1 to 1 ms, which indicatethe uncertainty used in conductivity inversion. Three grey vertical liensindicate time channels shown in Fig. 4.7.1044.1.3 Conductivity InversionConductivity inversion is the first step of the TEM-IP inversion workflow. Themain question in this step is: how do I ignore any IP effects when inverting theobserved data to obtain a conductivity model. Here, I simply choose time chan-nels from 0.1 to 1 ms (20 channels) for all sounding locations, which do not showany negatives. Using UBC TDoctree code [Haber and Schwarzbach, 2014], theseearly 20 channels of the observations are inverted, and a 3D conductivity model isrecovered. Parameters of the conductivity inversion are summarized in Table 4.2.The estimated conductivity model, sest(x,y,z), is shown in Fig. 4.10. The conduc-tive overburden, resistor (A1), and conductor (A2) are all recovered, but not a deepconductor (A4). This poor conductivity recovery of the deep conductor shows ageneral challenge of handling late-time conductivity signals when the observationis significantly contaminated by IP. Overall, the estimated conductivity is reason-able enough, and this will be used in the following steps: EM-decoupling and IPinversion.Table 4.2: Parameters of conductivity inversion.Parameters Valuesax = ay = az 1as 105f ⇤d 441 ⇥ 20 = 8820Uncertainty(e j) 101 |dobsj |m0 log(103)mre f log(103)Model weighting N/ABounds constraint N/A105Figure 4.10: Sectional views of the recovered conductivity model (sest).Solid lines delineate boundaries of the four IP blocks.1064.1.4 EM-decouplingThe conductivity model (sest) in Fig. 4.10 is used to estimate the fundamental re-sponses, F [sest ]. 2D maps of F [s•] and F [sest ] at three different times: 0.2 ms, 1.8ms, and 5.6 ms are shown in Fig. 4.11. At 0.2 ms, F [sest ] effectively estimated alinear edge of from the conductive overburden, and both anomalies at A2 (resistor)and A3 (conductor). Even at later times: 1.8 and 5.6 ms, F [s•] and F [sest ] showreasonable agreements. Note that the range of time used for the TEM inversion is0.1-1 ms, that is 1.8 and 5.6 ms are not included in the conductivity inversion. Fig.4.12 presents time decays of F [s•] and F [sest ] at four sounding locations close toA1-A4. The overall fit is good not only for times used in the inversion, but also atlater times: 1-10 ms. All four stations close to A1-A4 show a good match betweenF [s•] and F [sest ], and this will ensure effective EM-decoupling.The sought data dIP[sest ] can be obtained can be obtained bydIP[sest ] = dF [sest ] (4.2)The obtained IP data however, will include some errors since F [sest ] is not sameas F [s•]. To examine these more closely we writedIP[sest ] = dF [sest ] = dIP[s•]+4d+n, (4.3)where n is additive noise, and 4d (= F [s•]F [sest ]) is the error caused becauseof poor estimate of s•. Undoubtedly there are situations where the errors on theright hand side can become larger than dIP. This will always occur at early-timechannels where the IP response is small and EM induction response is large. Thusthere is an earliest time channel that can be used for this analysis. The secondissue concerns 4d. This is more difficult to quantify and needs to be treated on acase-by-case basis. I postpone this discussion until later in this section.I now proceed with EM-decoupling. As this is a synthetic example, both F [sest ]and F [s•] will be used, and the obtained dIP[sest ] and dIP[s•] are shown in Fig.4.13(a) and (b), respectively. A comparison shows that the relative strength of4dget smaller as time increases (e.g. 1.8 and 5.6 ms), whereas this error is large atearlier time (0.2 ms). Although there are significant errors, the IP anomaly at A3107is recognized even at 0.2 ms, and this was not shown as a negative anomaly ind (see Fig. 4.7a). At 1.8 ms d did not show negatives except for A3, but withthe EM-decoupling all four IP anomalies are recognized. At 5.6 ms dIP[s•] anddIP[s•] are almost identical. Fig. 4.14 compares time decays of dIP[s•] anddIP[sest ]; they have converged by 2 ms. To summarize, the EM-decoupling showsgood performance at the intermediate and late-times (see Fig. 4.5) when the IPsignal is significant, whereas it shows poor performance at early times when EMinduction is dominant.(a) True(b) EstimatedFigure 4.11: Comparison of true and estimated fundamental responses at 0.2,1.8 and 5.6 ms. (a) true fundamental response: F [s•] and (b) estimatedfundamental response: F [sest ].108Figure 4.12: Time decays of the d, F [s•], and F [sest ] at four sounding loca-tions close to A1-A4. Here solid and dashed lines distinguish positiveand negative datum. Error bars are shown from 0.1 to 1 ms, whichindicate uncertainty used in conductivity inversion.Ideas for Contending with4dAlthough dIP[sest ] = F [s•]F [sest ] data are obtained, still they include errors:4d. Here I comment about 4d, and suggest an idea that can reduce its impactof 4d. It is least important when dealing with resistive bodies and hosts, andmost problematic as the bodies and hosts become more conductive. There are afew items of note. Firstly, if sest is incorrect by a scale factor then this shifts thedIP data. Away from chargeable bodies the dIP response should be zero. Assum-ing these locations can be recognized, then the regional shift can be estimated andapplied to dIP[sest ]. The same idea is applicable to long wavelength spatial compo-nents of regional fields (4d). Any corrective procedure, which is akin to removalof regional fields in potential fields processing, relies on identification of areas inthe model believed to be free of IP responses.For example, I choose dIP[sest ] at 1.8 ms shown in Fig. 4.15a, and applya regional removal. Based on identified IP anomalies at A1-A4, some stations109(a) True(b) EstimatedFigure 4.13: Comparison of true and estimated IP responses at 0.2, 1.8 and5.6 ms. (a) true IP response: dIP[s•] and (b) estimated IP response:dIP[sest ].away from those area are selected to estimate regional fields. Selected locationsare marked as empty white circles in Fig. 4.15(a). A polynomial function wasused to fit dIP[sest ] at selected locations, and Fig. 4.15(b) shows the estimatedregional field. This is subtracted from dIP[sest ] to produce dIPRM[sest ] as shown inFig. 4.15(c). The dIPRM[sest ] shows a better match with dIP[s•]. Hence, a properregional removal could improve the quality of dIP[sest ]. Importantly, however extracare should be put on this step in practice since one could remove IP signals indIP[sest ] data. These are preliminary ideas at this stage, but the example showspromise and this is worthwhile pursuing.110Figure 4.14: Time decays of the d, dIP[s•], and dIP[sest ] at four soundinglocations close to A1-A4. Here solid and dashed lines or solid andempty circles distinguish positive and negative datum.Figure 4.15: Regional removal procedure. Left, middle and right panels cor-respondingly indicate before removal, estimated regional fields, andafter removal of dIP[sest ] at 1.8 ms. Black dots and empty white cir-cles indicate all stations and selected stations to estimate regional field.1114.1.5 IP InversionThe EM-decoupling method in the previous subsection has been used to generateIP responses, dIP. The next step inverts these data to recover a 3D distributionof pseudo-chargeability. 3D IP inversions applied to the dIP[sest ] (before regionalremoval) and the dIPRM[sest ] (after regional removal) at 1.8 ms are shown in Fig.4.15(a) and (c). Parameters of the IP inversion used are shown in Table 4.3.The recovered pseudo-chargeability models are shown in Fig. 4.16(a) and (b).Without the regional removal, the four chargeable blocks are reasonably imaged,but there is a large-scale rectangular artefact contaminating the image. This is dueto errors in the EM-decoupling caused by an inaccurate sest . These artifacts areeffectively removed when the dIP[sest ] are processed to remove a regional field.This result indicates that the recovered pseudo-chargeability can provide usefuldistributed IP information of the subsurface. The inversion results in Fig. 4.16provide insight about the existence and the geometry of the chargeable rocks. How-ever, there are two important aspects that should be remembered when interpretingthe resultant images:• Relative magnitude of h˜ . The recovered h˜ is different from h , and higherh˜ does not necessarily indicate higher chargeability.• Depth resolution of the IP inversion. A depth weighting is used to com-pensate for lack of intrinsic depth resolution. As shown in Fig. 4.16, relativedepth information can be extracted from the IP inversion. This is the sameissue as in magnetic inversion, and numerous successes of magnetic inver-sion supports that the inversion with depth weighting still provides relativedepth information.Extracting intrinsic IP parametersIn the above section I have concentrated upon IP inversion at a single time channel.The same IP inversion is now applied to each time channel of dIP[sest ], separately,to recover pseudo-chargeability at multiple times. Here pseudo-chargeability is ∂ h˜∂ t because the data are voltages, which is similar to  ∂~b∂ t and hence its dimen-sion is 1/s. From this pseudo-chargeability:  ∂ h˜∂ t (x,y,z; t), I extract Cole-Cole pa-112rameters by solving a small inverse problem as described in Appendix B. Pseudo-chargeability at a single cell is parameterized in terms of h and t by assuming cis known and is set to c=1. A representative cell, which has a maximum pseudo-chargeability value, for each of the four anomalies is selected, and correspondingdecays are shown in Fig. 4.17. Those four time decays of the pseudo-chargeabilityare data for a parametric inversion to recover h and t for each decay. Comparisonof observed and predicted pseudo-chargeability at the four cells are shown in Fig.4.17. Recovered values and selected locations are summarized in Table 4.4. Therecovered t’s for cells corresponding to A1-A3 are in fair agreement with true val-ues but A4 differs by a factor of four. Poor recovery of t for the A4 may be causedby poor estimation of conductivity for this deep conductor. The h values howeverare all underestimated (see Table 4.1). This result is consistent with the exampleshown in Section 3.3.2. I conclude that there is the potential to extract intrinsicCole-Cole parameters from AEM data but the time constant seems to be the mostrobust parameter to be extracted.Table 4.3: Parameters of IP inversion.Parameters Valuesax = ay = az 1as 105Uncertainty(e j) 102 max(dobs)m0 0mre f 0Model weighting 1/z3Bounds constraint m> 0113(a)(b)Figure 4.16: Sectional views of the recovered dh˜/dt models at 1.8 ms. (a)Without regional removal. (b) With regional removal.Table 4.4: Recovered Cole-Cole parameters of four chargeable blocks (A1-A4). Values in parenthesis indicates true values.Division A1 A2 A3 A4Cell location (-575, 625, -25) (-575, -525, -25) (625, 575, -75) (-625, -575, -25)h 0.39 (0.5) 0.033 (0.5) 0.006 (0.4) 0.01 (0.8)t 0.006 (0.005) 0.0035 (0.005) 0.0036 (0.005) 0.0019 (0.005)114Figure 4.17: Comparison of observed and predicted pseudo-chargeability forfour chargeable anomalies: A1-A4. Solid lines and empty circles dis-tinguish observed and predicted pseudo-chargeability.1154.1.6 ConclusionsThe TEM-IP workflow is applied to a synthetic AEM example, and each of threesteps is carefully analyzed and tested. The first step is conductivity inversion, andearly-time channels which do not have negative values are used to recover 3Dconductivity, sest . Overall the conductivity structures are imaged well, except forthe deep conductor, A4. Using the obtained conductivity, the second step is EM-decoupling. By subtracting F [sest ] from d, I obtain dIP[sest ]. The EM-decouplingshows great performance in the intermediate and late-times, but as expected it givespoor results in the early-time. Four IP anomalies are obtained, even at 1.8 ms whenthe observational data are positive for some blocks. To illustrate an approach toreduce the error of EM-decoupling characterized at early times, a regional removalis applied by assuming the error is a long-wavelength spatial component of a re-gional field. This regional removal procedure helped and it may be worthy of fol-low up research. The third step is IP inversion. Each time channel of the obtaineddIP[sest ] is inverted to generate 3D pseudo-chargeability models at multiple times;they provide useful information about the distributed polarization of the subsur-face. Furthermore, by interpreting the recovered pseudo-chargeability, Cole-Coleparameters: h and t are extracted for the four anomalies: A1-A4. Recovery of twas reasonable, but h is mostly underestimated. To conclude, this synthetic ex-ample demonstrates the capability of the TEM-IP inversion workflow to handleairborne IP data, and recover distributed polarization information in 3D.1164.2 Field Example: Mt Milligan4.2.1 SetupMt Milligan is located in BC, Canada (Fig. 4.18). It is a copper and gold por-phyry deposit within the Early Mesozoic Quesnel Terrane, a Late Triassic to EarlyJurassic magmatic arc complex that lies along the western North American con-tinental margin. Previous geologic and geophysical work has illuminated that themineralization is associated with monzonite stocks that intruded into basaltic vol-caniclastic rocks Oldenburg et al. [1997], Yang and Oldenburg [2012][Reference].A VTEM survey (2007) was flown over the Mt Milligan porphyry deposit, andit uses a helicopter-borne system with co-located source loop and receiver loop,which is similar to the coincident-loop system. This VTEM survey is a part of theQUEST project, which covers a much broader region; roughly 514 km in lengthby 124 km in width as shown in Fig. 4.18. The Mt Milligan VTEM survey has13 lines, each 2.7 km in length resulting in a total of 37.5 line-km. Voltages aremeasured at 27 time channels ranging from 0.099 to 9.3 ms; the current waveformis shown in Fig. 4.20. Survey lines of the VTEM survey are shown in Fig. 4.19as grey lines. That figure also presents the geology of the region and structuralfeatures. The topography of the area is shown in Fig. 4.21 and it ranges from 1044to 1419 m. The flight height of the EM loops was about 30 m.An effective 3D TEM inversion was developed by Yang et al. [2014], and theyinverted the VTEM data, and recovered a 3D conductivity model. The most promi-nent intrusive monzonite stock, MBX, was successfully recovered as a resistivetarget. A conductive alteration zone surrounding this resistive stock, often calleda halo structure, was recovered. This was an important example for the industrybecause it showed the necessity of 3D AEM inversion. 1D inversion could notrecover this 3D conductivity structure, and it generated a misleading conductivefeature at the location of the MBX stock [Yang and Oldenburg, 2012]. Negativetransients were observed at some stations (see Fig. 4.2), and these stations werenot included in the 3D AEM inversion. These negative transients suggest the exis-tence of chargeable materials in the Mt Milligan area, and this was also supportedby DC-IP inversion conducted by Oldenburg et al. [1997]. Based upon the TEM117inversion done by Yang et al. [2014], I apply the TEM-IP inversion workflow to theVTEM data over the Mt Milligan area. This is the first application of the TEM-IPinversion workflow to a field example, and an attempt to recover 3D chargeabilityinformation. The focus of this field AEM example will be EM-decoupling and theIP inversion since the 3D conductivity model has already been obtained by Yanget al. [2014].Figure 4.18: Location of Mt Milligan and QUEST survey area.118Figure 4.19: Geology and VTEM survey at Mt Milligan porphyry deposit inBritish Columbia, Canada [Yang and Oldenburg, 2012].Figure 4.20: Normalized current waveform of the VTEM survey (2007) atMt Milligan deposit.119Figure 4.21: Topography of the VTEM survey of the Mt Milligan region. Ge-ology indicated with white lines.1204.2.2 DataFrom the VTEM survey 14,760 stations of the TEM data are collected and eachstation has 27 time channels. Fig. 4.22 shows a 2D map of the observed data ateight time channels ranging from 0.2 to 2.7 ms. At 0.2 ms, high anomalies areshown in the central western side and the north east side of the map. High voltagesat the east side of the region are due to highly conductive Tertiary sediments. Thecentral anomaly may be due to a conductive alteration zone around the resistiveMBX stock. At later times, this high positive anomaly changes to negative valuesat 1 ms. Other negative anomalies are shown at later times. In the left panel ofFig. 4.23, three negative anomalies are marked as white dashed circles and namedA1-A3. A2 is the main negative anomaly, which is the positive high anomaly inthe central west discussed above. Corresponding time decays are shown on theright panel of Fig. 4.23. Decays at A2 show the greatest negative values indicat-ing strong IP effects. Both A1 and A3 show negatives, whereas their amplitudesare quite small: ⇠ 104 pV/A-m2, which is about the noise level for this survey;and hence it may be questionable to consider them as negative anomalies due tochargeable materials. However, the A2 anomaly is definitely strong enough to beconsidered as an IP signal. We will test this by applying the TEM-IP workflow tothis VTEM data set, and recover 3D chargeability information over the MtMilligandeposit.121Figure 4.22: 2D map of the observed VTEM data at eight different timesranging from 4.6-7.1 ms. Solid white lines show boundaries of dif-ferent gelogical units and the white areas indicate where the responsehas gone negative.122Figure 4.23: 2Dmap of the VTEM data at 2.7 ms (left panel) and time decaysat three stations close the negative anomalies at A1-A3 (right panel).Solid white lines on the right panel show boundaries of different geo-logical units, and white dashed lines indicate the approximate locationsof where there are.1234.2.3 Conductivity InversionBecause detailed information about the 3D conductivity inversion of the VTEMis described in Yang et al. [2014], only a brief summary of the inversion is pro-vided here. In the total 14,362 soundings (except for some soundings showingnegatives (around A1)), and 8 time channels ranging from 0.2 to 2.7 ms are usedin the 3D TEM inversion. Five earlier time channels are ignored due to lack ofknowledge about the exact current waveform; after 2.7 ms the signals are noisy.For data uncertainty, 10 percent of the data and a floor of 1013 (V/A) floor areused. Fig. 4.24 shows plan and sections views of the estimated conductivity, sest .The resistive MBX stock and the surrounding conductive halo structure are suc-cessfully recovered. Note that the region where the conductive halo structure isrecovered is where negatives are observed: A1-A3 anomalies. Using the obtainedsest , EM-decoupling will be applied to obtain IP data.Figure 4.24: 3D Conductivity model of Mt Milligan porphyry deposit ob-tained by Yang et al. [2014].1244.2.4 EM-decouplingAlthough negatives are measured at later time channels (¿ 1 ms), the observationsstill include EM effects, and these need to be removed to recognize IP signals inthe data. For this, the EM portion is estimated by forward modelling TEM datausing sest . The dIP[sest ] = dF [sest ].Fig. 4.25 illustrates the EM-decoupling applied to the observed data, d atthree different time channels: 2.7 ms, 0.7 ms, and 0.2 ms. In the figure, left,middle and right panels respectively show d, F [sest ], and dIP[sest ]. At 2.7 ms,F [sest ] reasonably estimates d except for the region where negatives are shown.The dIP[sest ] plotted in the right panel of Fig. 4.25(a) show five anomalies namedA1-A5. Note that A4-A5 did not show any negatives in the observations, whereasA1-A3 did. The EM-decoupling shows good performance even at earlier times(0.7 ms and 0.2 ms), as shown in Fig. 4.25(b) and (c). Fig. 4.26 shows timedecays at A1-A5. Overall, d and F [sest ] show reasonable agreements at earlytimes, but the discrepancy between them increases with time (see A2 in Fig. 4.26).The EM-decoupling converts this to an IP signal. However, the obtained dIP[sest ]can include errors originating from discrepancies between s•(x,y,z) and sest , andhence some caution is necessary; I postpone this issue to the following section.To summarize, the EM-decoupling process is effectively applied to the Mt Mil-ligan VTEM data, and dIP[sest ] are obtained, which illuminated five IP anomalies:A1-A5. In the following section, the obtained dIP[sest ] will be inverted to recovera 3D pseudo-chargeability in the Mt Milligan area.125(a) 2.7 ms(b) 0.7 ms(c) 0.2 msFigure 4.25: EM-decoupling of the VTEM data over Mt Milligan region. 2Dmaps of d, F [sest ], dIP[sest ] are shown at (a) 2.7, (b) 0.7 and (c) 0.2ms. F [sest ] is subtracted from d to obtain dIP[sest ]. A4-A5 are IPanomalies recognized by EM-decoupling, whereas A1-A3 were shownas negatives in d. A1-A3 and A4-A5 are marked as grey and blackdashed circles, respectively.126Figure 4.26: Time decays of d, F [sest ], dIP[sest ] at the five IP anomalies: A1-A5. Top left panel showed 2D map of dIP[sest ] at 2.7 ms.1274.2.5 IP InversionFrom the EM-decoupling, eight time channels of the dIP[sest ] (from 0.2 ms-2.7ms) are obtained. Each time channel of the dIP[sest ] is inverted separately, and a3D distribution of pseudo-chargeability is recovered at each of eight channels. Pa-rameters of the inversion used here are the same as those used in Section 4.1.5 asdescribed in Table 4.3. Fig. 4.27(a)-(c) shows plan and section views of the recov-ered pseudo-chargeability at three different time channels: 2.7 ms, 0.7 ms, and 0.2ms, respectively. At 2.7 ms anomalous pseudo-chargeability volumes are recov-ered around A1-A5, and in particular the pseudo-chargeability around A2 and A3show strong amplitude. Most anomalies are imaged 100 m below surface. Regionalremoval can be applied (as described in Section 4.1.4), and it could enhance thequality of the obtained dIP[sest ]. In contrast, it could remove the IP signals, whichindicates that extra care is necessary for the regional removal. Hence, I omittedthis step, and to prevent fitting those errors in the dIP[sest ], I imposed the positivityconstraint on the IP inversion. This positivity constraint effectively ignores mostof positive bias in the dIP[sest ] as shown in Fig. 4.28. In this figure, the observedand predicted data of the IP inversion are compared at three time channels. Earliertime shows greater misfits, which is a consequence of the greater the errors in theEM-decoupling process.128(a) 2.7 ms(b) 0.7 ms(c) 0.2 msFigure 4.27: Plan and section views of the recovered pseudo-chargeability:(a) 2.7, (b) 0.7, and (c) 0.2 ms. Five anomalies recognized fromdIP[sest ] are marked as dashed circles on plan view (left panels).129Figure 4.28: Comparison of observed and predicted IP data for the IP inver-sion at three time channels: 2.7, 0.7 and 0.2 ms.1304.2.6 InterpretationUsing the obtained conductivity and chargeability, in this section, 3D rock modelswill be constructed. A pseudo-chargeability at 2.7 ms is used for this interpretation;it is the latest time that has the greatest quality of dIP[sest ]. Fig. 4.29 showsplan and section views of the 3D conductivity and pseudo-chargeability with somegeological references (e.g., MBX stock and Rainbow fault). The MBX stock isimaged as an isolated resistor, and conductive halo structures are recovered aroundtheMBX stock . Tertiary sediments on the eastside of the inversion are well imagedas a highly conductive medium. High pseudo-chargeability anomalies, A1-A5, areimaged on the edges of three mineralized zones: MBX, DWBX, and SoutherneStar. In particular around the Rainbow fault, high pseudo-chargeability anomaliesare imaged and they show good alignment with the Rainbow fault. Note that ahigher degree of alteration is expected around faults surrounding the stocks, andthey could potentially have strong polarization effects.With the conductivity and pseudo-chargeability models, a 3D rock model hav-ing three units: R1-R3 is constructed and presented in Fig. 4.30. Each of rockunits can be interpreted as:• R1: Highly conductive medium, which correspond to both alteration zonesand Tertiary sediments• R2: Isolated resistive volumes at depth, which indicates resistive monzonitestocks (e.g MBX stock)• R3: Highly chargeable medium, which may corresponds to highly alteredzones especially near the surface.The above rock model is reasonable, but there can be multiple ways to constructa rock model. For instance, if the emphasis is on the resistive stock and alterationzone, one could obtain a different rock model that looks like Fig. 4.30. Here R1corresponds to the resistive monzonite stock, and R2 is the alteration zone comingfrom the obtained chargeability information.131(a) Conductivity(b) Pseudo-chargeabilityFigure 4.29: Plan and section views of 3D conductivity and pseudo-chargeability over Mt Milligan area.132Figure 4.30: The 3D rock model obtained from both conductivity andpseudo-chargeability models. Red and black dots indicates fault struc-tures, and three mineralized zones, respectively. R3 shows the rockunits obtained from the airborne IP data.133Figure 4.31: The 3D rock model obtained from both conductivity andpseudo-chargeability models. Red and black dots indicates fault struc-tures, and three mineralized zones, respectively. R2 shows the rockunits obtained from the airborne IP data.1344.2.7 ConclusionsThe VTEM data over the Mt Milligan porphyry deposit have negative anomalies atlate-times, which are due to chargeable materials. Three negative anomalies: A1-A3 are recognized in the observations. By applying the TEM-IP inversion work-flow, 3D distributions of both conductivity and pseudo-chargeability are obtained.Negatives in the observations are ignored when inverting observations to recover3D conductivity distribution. The obtained 3D conductivity, sest is used to carryout EM decoupling to yield dIP[sest ] at eight time channels. Two more IP anoma-lies: A4-A5 around the Rainbow fault are recognized in the obtained dIP[sest ];these were not identified in the observations. Using the IP inversion, each timechannel of the dIP[sest ] is separately inverted, and pseudo-chargeability is recov-ered at eight time channels. Highly chargeable volumes are successfully imagedon the edges of three mineralized zones, and in particular three of anomlies showedgood alignment with the Rainbow fault. These chargeable volumes are interpretedas highly altered zones. Using the conductivity and pseudo-chargeability model at2.7 ms, 3D rock models are constructed. Importantly, the pseudo-chargeability pro-vided valuable information to characterize highly altered zones of the monzonitestock.To conclude, the TEM-IP inversion workflow is successfully applied to theVTEM data over Mt Milligan resulting in both conductivity and chargeability mod-els in 3D. This result demonstrates the capability of the workflow to handle fielddata. Moreover, I believe this was the first time that a 3D chargeability distributionwas recovered from field AEM data using a 3D voxel inversion, and the resultantpolarization information is promising to be used for characterizing a porphyry de-posit. The final step of the workflow, extracting intrinsic IP information, has notbeen applied here yet but will be applied in the near future.1354.3 Field Example: DO-27/18 kimberlites4.3.1 SetupThe Tli Kwi Cho (TKC) kimberlites were identified from a DIGHEM survey in1992 (Fig. 4.32). The kimberlites are located 28 km southeast of the Diavik minein the Lac de Gras region, Northwest Territories, Canada. The initial discoverytargeted two anomalies, named DO-18 and DO-27. Following the initial discovery,several generations of EM systems have been deployed over the TKC area in aneffort to characterize the kimberlites. In 1999, the first TEM survey was carriedout using the AeroTEM I system [Boyko et al., 2001]. Negative transients weremeasured, in particular at DO-18, although it was not clear whether these weretrue signals from the earth or instrumental noise. Surveys with new generations ofequipment, AeroTEM II (2003) and VTEM (2004), reaffirmed the negatives. Inaddition, a ground loop NanoTEM (1993) survey showed negatives at the DO-18pipe [Jansen and Doyle, 2000]. Airborne TEM systems and NanoTEM have sim-ilar geometry and can be considered to be coincident loop systems and hence thenegatives are indicative of chargeable material [Weidelt, 1982]. From the perspec-tive of kimberlite exploration however, the existence of an IP signal is not neces-sarily significant. Ice and near surface clays are known to be chargeable. Theirpresence distorts EM signals and impede interpretation [Smith and Klein, 1996,Kozhevnikov and Antonov, 2012]. As such, the existence of negative transientsis usually considered to be “noise” and it is commonly referred to as IP contam-ination. Recent studies however, have suggested that negative transients could beattributed to more interesting geological features and thus the negative transientsare “signal” [ElKaliouby and Eldiwany, 2004, Flores and Peralta-Ortega, 2009,Kang et al., 2014]. It is this potential that I wish to pursue in this example.Kimberlite pipes in the Lac de Gras region are generally excellent geophysicaltargets as they exhibit high physical property contrasts with the granitic host rocks:higher magnetic susceptibility (with high remanence), lower density and higherconductivity [Power and Hildes, 2007]. The standard geological model adoptedhere for kimberlites consists of three different kimberlitic rocks: hypabyssal kim-berlite (HK), volcaniclastic kimberlite (VK), and pyroclastic kimberlite (PK) as136shown in Fig. 4.33 and summarized in Table 4.5. Since each kimberlite unitshows different physical property character, recovering 3D distributions of density,susceptibility, and conductivity will help characterize different kimberlites. In ad-dition, adding chargeability could help this characterization in the TKC region, andthis is the motivation to recover 3D chargeability information from the AEM data.For instance, PK and VK are very similar in terms of density, susceptibility, andconductivity. Hence, if chargeability could help distinguish PK and VK, that willbe valuable information for diamond exploration.In this example, I concentrate upon extracting both conductivity and charge-ability information from the VTEM data set. The TEM-IP inversion workflow hasthe three steps: (1) Conductivity inversion, (2) EM-decoupling, (3) 3D IP inver-sion. In a final stage I compare our petrophysical interpretation, based solely onAEM data, to the extensive drilling data available over the deposit.DO−27DO−187 13 30 007 13 40 007 13 50 00N or t hi n g ( m)556000 557000 558000 559000Easting (m)TKCNTNUBC AB SK MBFigure 4.32: Location map for the Tli Kwi Cho (TKC) kimberlites, NWT.DO-18 and DO-27 are two main kimberlite pipes at TKC region. PK,VK, and HK correspondingly indicate pyroclastic kimberlite (PK),volcaniclastic kimberlite (VK), and hypabyssal kimberlite (HK).137Figure 4.33: Schematic diagram of a kimberlite pipe in the Lac de Gras re-gion (Modified from Devriese et al. [2017]). A lake may be presentafter glaciation and is often used as a first indicator of a possible kim-berlite. Transverse lines are from the DIGHEM survey (1992).Table 4.5: Expected physical property contrast for kimberlite deposits in theLac de Gras region [Power and Hildes, 2007].Rock type Density Susceptibility ConductivityGlacial till moderate none moderate-highHost rock moderate none lowHypabyssal kimberlite (HK) low-moderate high low-moderateVolcaniclastic kimberlite (VK) low low-moderate moderate-highPyroclastic kimberlite (PK) low low-moderate moderate-high1384.3.2 DataNotable features about the TKC kimberlites can be identified, by simple visualinspection of the EM data. From the DIGHEM data at 7200 Hz, shown in Fig.4.35(a), positive anomaly highs are observed over the location of DO-18 and DO-27; this indicates that both pipes are more conductive than the host granitic rocks.Fig. 4.34 (b), (c), and (d) correspondingly show NanoTEM (77 µs), AeroTEMII (26 µs), and VTEM (90 µs) data. Here, all three TEM surveys show a positiveanomaly near DO-27. I also identify three important anomalies that were not cap-tured by the DIGHEM system and where the data are negative: A1 near DO-18,A2 between DO-18 and DO-27, and A3 near DO-27. Fig. 4.35 (b)-(d) shows allthree TEM data sets at later times. Both NanoTEM (603 µs) and AeroTEM II data(534 µs) are significantly noisy. In the VTEM data (680 µs) however the area thatwas previously positive at early times within DO-27 has switched to negative. Thisis referred as the A4 anomaly and it suggests that DO-27 has chargeable material.Data quality, and the time range for which data are sampled, vary across EMsystems, hence the EM data sets should show some differences. Fig. 4.36 (a)and (b) respectively provide transients of NanoTEM and VTEM data at severalsounding locations at DO-18. Soundings taken away from the pipe are referredto as ‘background’. In NanoTEM data, all transients show negative values, butthe negative transients from the background decay faster than those over the DO-18 pipe. The IP signal in the background soundings is likely due to surface glacialsediments (including ice and clays). These background negative transients were notidentifiable in the VTEM data, likely because the VTEM system does not extendas early in time and also the survey equipment is higher off of the ground.Since the VTEM data set includes most of the important IP features observedat TKC, while showing less noise at later time channels than other TEM data, Ifocus our analysis on the VTEM data. From those data, four negative anomalies ofinterest (A1-A4) are identified, and they appear to have different decaying rates asshown in Fig. 4.36 (b).139Figure 4.34: Plan maps of four EM data sets at TKC: (a) DIGHEM (56000Hz), (b) NanoTEM (77 µs), (c) AeroTEM II (26 µs), (d) VTEM (90µs). For TEM data sets, a smaller region (red box) close to DO-18and -27 is presented. The black line is a contour line of the nega-tive anomaly (-8 nT/s) from AeroTEM data at 26 µs. The white lineshows the boundary of the lakes. Negative anomalies: A1-A4 are cor-respondingly marked as purple, yellow, red, and green solid circles;A1-A3 showed strong negatives for all TEM data at an early-time asshown in Fig. 4.35.140Figure 4.35: Plan maps of four EM data sets at TKC: (a) DIGHEM (7200Hz), (b) NanoTEM (603 µs), (c) AeroTEM II (534 µs), (d) VTEM(680 µs). For TEM data sets, a smaller region (red box) close to DO-18 and -27 is presented. The black line is a contour line of the nega-tive anomaly (-8 nT/s) from AeroTEM data at 26 µs. The white lineshows the boundary of the lakes. Negative anomalies: A1-A4 are cor-respondingly marked as purple, yellow, red, and green solid circles;A1-A3 showed strong negatives for all TEM data at an early-time asshown in Fig. 4.35.141Figure 4.36: Transient curves of NanoTEM and VTEM data. (a) NanoTEMsoundings away from the pipe and representative of background (bluelines) and over the DO-18 pipe (black lines). They are marked as blueand black solid circles in Fig. 4.34(b). (b) VTEM soundings at A1-A4 ( correspondingly purple, yellow, red, and green lines). They aremarked as purple, yellow, red, and green solid circles in Fig. 4.34(d) and Fig. 4.35(d). Solid and dashed lines distinguish positive andnegative observations.1424.3.3 Conductivity InversionThe first step in the workflow is to estimate the background conductivity, s• fromthe TEM data. Detailed information about this step is described in Fournier et al.[2017] and hence only brief summary is provided here. All of the DIGHEM data,and VTEM data that were positive, were cooperatively inverted. In some areas nearDO-18, even the earliest time channel was negative so only DIGHEM data couldbe used there. The neglect of potential IP contamination in the DIGHEM data,and the likelihood of IP contamination in VTEM time channels that were positive,probably contributed to the difficulty we had in obtaining a single conductivitydistribution that fit both data sets. Nevertheless, the two models found through thecooperative inversion were very similar. For my purposes here, where I concentrateon extracting information from the VTEM data, I use the inversion model obtainedfrom the last step of the cooperative inversion where the starting and referencemodel was from the DIGHEM data but the data to be fit were the VTEM data.Fig. 4.37(a) and (b) shows the estimated conductivity model from the VTEMand the DIGHEM data, respectively. The two conductive pipes are imaged at depth.The conductive pipe for DO-27 extends deeper than the pipe for DO-18. Fig. 4.38shows plan maps of the observed (d) and estimated fundamental responses (F [sest ])at 130 µs. Regions in white correspond to negative data and these soundings werenot used in the cooperative inversion. Therefore, the conductivity structure nearA1-A3 shown in Fig. 4.37 is mostly coming from the DIGHEM data. The majorregion where both the VTEM and DIGHEM data contributed to the final conduc-tivity is near A4; d and F [sest ] in Fig. 4.38 show a good match there.143Figure 4.37: Plan and section views of the recovered conductivity modelfrom the cooperative inversion of the VTEM and DIGHEM data sets:(a) VTEM and (b) DIGHEM. The white outlines delineate boundariesof the lake. The green outlines show the extent of DO-27 and DO-18at the surface, based on drilling.144Figure 4.38: Observed and estimated fundamental responses at TKC. Foursounding locations at A1-A4 are marked as solid circles. Regions hav-ing negative values are shown as white.1454.3.4 EM-decouplingBased upon the estimated conductivity, sest from the previous conductivity inver-sion (Fig. 4.37a), EM-decoupling is applied: dIP[sest ] = dF [sest ]. I expectthat the EM-decoupling step will only be effective at intermediate times when bothEM and IP effects are considerable (see Fig. 4.5). Note that at early times (EM-dominant; positive data), dIP is too small to be obtained, and at late-times (IP-dominant; negative datum) EM-decoupling will have a minor impact. Fig. 4.39shows time decaying curves of d (black), F [sest ] (blue), and dIP[sest ] (red) at A1-A4. For A1-A3, as shown in Fig. 4.39 (a)-(c), even the earliest observation has anegative sign (IP-dominant time). The magnitude of F [sest ] is much smaller thandIP[sest ] except for the few earliest times, hence the EM-decoupling is not nec-essary. However, at A4 as shown in Fig. 4.39 (d), a full suite of EM-dominant,intermediate, and IP-dominant times are captured. The EM-decoupling is mosteffective in the intermediate time (200-1000 µs).To illustrate performance of the EM decoupling I focus on two time channels:130 and 410 µs. Plan view maps of d, F [sest ], and dIP[sest ] for the times areshown in Fig. 4.40 (a) and (b). At 130 µs, near A4 the positive high anomaly(from the conductive DO-27 pipe) is effectively removed resulting in some smalllow amplitude IP features. Near A1-A3 the EM-decoupling results in strongernegatives. At 410 µs, around A4, the EM-decoupling makes a greater impact,and it converts a small area with positive observations within negatives to a largernegative amplitude area. Having separated the EM and IP signals in the VTEMdata, the obtained dIP[sest ] at each time channel can now be inverted to recover a3D pseudo-chargeability. The inversion will be carried out for all time channels.146Figure 4.39: Time decaying curves of the d (black line), F [sest ] (blue line),and dIP[sest ] (red line) data at (a) A1, (b) A2, (c) A3, and (d) A4.Solid and dashed lines distinguish positive and negative values. Verti-cal black dashed line indicate 130 and 410 µs, respectively.147Figure 4.40: Plan maps of d, F [sest ] and dIP[sest ] at (a) 130 and (b) 410 µs.Left, middle, and right panels correspondingly show the d, F [sest ] anddIP[sest ].1484.3.5 IP InversionEach time channel of dIP[sest ] is separately inverted to recover a pseudo-chargeabilitymodel. The recovered pseudo-chargeabilities at two time channels: 130 µs and 410µs are shown in Figs 4.41(a) and (b), respectively. At 130 µs three chargeable bod-ies are recovered: A1-A3; A1 and A3 are imaged at DO-18 and the north-easternpart of DO-27, respectively. Dashed contours (red) shown in Fig. 4.41(a) delineatethree chargeable bodies. No chargeable volume is recovered around the southernpart of the DO-27 (see green lines at DO-27 shown in Fig. 4.41a). At 410 µshowever, a chargeable body is imaged at this southern part of the DO-27 (A4);solid contour (red) indicates this chargeable volume. Other chargeable structuresare imaged at 410 µs around A1 and A2, but their amplitudes are much smallerthan A4. This reflects the different time decaying features of the IP signals: A1-A3decay faster than A4 (see Fig. 4.36).The recovered pseudo-chargeability at 130 µs and 410 µs provides some im-portant information about the different kimberlites in the region. However, thatIP information is still qualitative and hence it motivates the following quantitativeanalysis.149Figure 4.41: Plan and section views of the recovered pseudo-chargeabilitymodel from the 3D IP inversion of the raw IP at (a) 130 µs and (b)410 µs, respectively. Left panel shows plan map at 99 m below sur-face. Top and bottom of right panels show A-A’ and B-B’ sections,respectively. Solid and dashed red lines delineate contours of the re-covered pseudo-chargeability at 50 s1 (130 µs) and 10 s1 (410 µs).The black outlines delineate boundaries of the lake. The green outlinesshow the extent of DO-27 and DO-18 at the surface, based on drilling.150Extracting intrinsic IP parametersA distribution of pseudo-chargeability values at multiple times has been recoveredand I now wish to use those results to extract intrinsic information about the po-larization parameters of the kimberlites. Pseudo-chargeabilities at representativecells close to A1-A4 anomalies are chosen, and plotted as a function of time inFig. 4.43 (a). Using the procedure described in Appendix B, I fit  ∂ h˜∂ t (t) usingCole-Cole model, which can be expressed ass(w) = s• s•h1+(1h)(ıwt)c , (4.4)where w is angular frequency (rad/s), s• is the conductivity at infinite frequency(S/m), h is the chargeability, t is the time constant (s), c is the frequency exponent.Fig. 4.42 shows an example Cole-Cole model as a function of frequency. EachCole-Cole parameter determine different characters of s(w):• Magnitude of the s(w) is dominated by s• at high frequencies.• h controls the difference of s(w) between low and high frequency asymp-tote.• t mainly controls the peak frequency (see imaginary part in Fig. 4.42), andthis effectively controls the characteristic rate of decay in time. For instance,as t decreases the peak frequency increases and the polarization responsedecays faster in time.• c controls the spread in time decays.In these analyses, two Cole-Cole parameters: s• and c are assumed to beknown, hence I only estimated h and t in this inversion. For s•, the sest ob-tained from the conductivity inversion is used. I empirically used c=1 for cellsclose to A1-A3 and c=0.5 for those close to A4. Justfication of this was trialand errors. I could not fit A4 using c=1, and vice versa. Solid circles in Fig.4.43(a), show the predicted pseudo-chargeability; they match well with the ob-served pseudo-chargeability (lines). Median values of the pseudo-chargeability ateach time channel are shown in Fig. 4.43(b), and they demonstrate the different151rates of decay between A1-A3 and A4. A summary of the estimated t and h ateach cell is presented in Table 4.6. Fig. 4.44 (a) shows a cross-plot of the esti-mated t and h for A1-A4. Clustering of A1-A4 indicates distinctions between thedifferent kimberlite units based upon the estimated t and h :• A4 can easily be distinguished from others on the basis of t• A1 and A3 can be differentiated by h and perhaps by t• The distinction between A1 and A2 is subtle, but it may be possible basedupon t valuesFigure 4.42: Cole-Cole conductivity as a function of frequency [Pelton et al.,1978]. Black and red lines indicate real and imaginary part of the Cole-Cole conductivity. Used Cole-Cole parameters are: s•=1 S/m, h=0.2,and c=0.5; two t values are considered: 103 s and 104 s. Solid anddashed line indicates complex conductivity generated with t = 103 sand t = 104 s, respectively.152Figure 4.43: Comparison of the observed (lines) and predicted (solid circles)pseudo-chargeability at cells close to A1-A4 (correspondingly purple,yellow, red, and green colors). (a) All time decaying curves of theobserved and predicted pseudo-chargeability. (b) Median values of theobserved and predicted pseudo-chargeability at each time.Table 4.6: Mean and standard deviation of estimated t (µs) and h for A1-A4.Used values of c for fitting are given.Division Mean t Mean h Std. t Std. h cA1 110 0.4 24 0.13 1A2 50 0.4 9 0.11 1A3 40 0.12 9 0.07 1A4 1600 0.15 380 0.1 0.5153Figure 4.44: A cross plot between the estimated time constant (t) and charge-ability (h) at cells close to A1-A4. Solid circles shaded as purple,yellow, red, and green colors correspondingly indicate A1-A4.1544.3.6 InterpretationFrom the application of the TEM-IP inversion workflow to the VTEM data overthe TKC region, both conductivity and chargeability information are recovered in3D. Using these, along with basic geological information given for kimberlites (e.g.Table 4.6), I can build a 3D rock model for the TKC kimberlites. For this, I will useDIGHEM conductivity model shown in Fig. 4.37(b), and the pseudo-chargeabilitymodels at 130 and 410 µs (Fig. 4.40). The pseudo-chargeability model for 130µs and 410 µs are correspondingly referred as the early pseudo-chargeability, h˜E ,and the late pseudo-chargeability, h˜L.Compared to the host rock, kimberlites are supposed to have higher conduc-tivity (Table 4.5), and hence conductivity can provide the volumes of kimberliticmaterials; chargeability information about kimberlites is not well-known. Fromthe conductivity model shown in Fig. 4.37(b), I extracted moderately conduc-tive volumes (< 1250 W m), and overlaid them in to Fig. 4.45(a) as grey regions;they indicate overall kimberlites at DO-18 and DO-27. From the early pseudo-chargeability (h˜E >55 s1), I obtained two main chargeable volumes close to A1and A3 (red regions); they have small time constant (t ' 75 µs). The late pseudo-chargeability (h˜L >8 s1) provides another chargeable volume at the southern partof the DO-27 (green region); it has a large time constant (t ' 1600 µs). Note thatall chargeable anomalies are within the conductive volume (grey region).By using the above anomalous volumes overlaid in the Fig. 4.45(a), a 3D rockmodel having four different rocks types (R0-R3) is constructed and presented inFig. 4.45(b). Each rock unit has different character with regard to conductivity andchargeability.• R0: Low conductivity and non-chargeable• R1: Moderate-High conductivity and non-chargeable• R2: Moderate-High conductivity and chargeable (small t)• R3: Moderate-High conductivity and chargeable (large t)This 3D rock model provides some important information about the TKC kimber-lites. First, the DO-27 pipe is embedded at depth, but the DO-18 pipe is exposed155to the surface. Second, the DO-18 and DO-27 kimberlite pipes are different basedupon distinction between R2 and R3; this supports that the two pipes are distinctevents. Third, at the DO-27 pipe there are at least two different kimberlites (see R2and R3 shown in the right-bottom panel of the Fig. 4.45b).One remaining question here is how do we interpret rock units obtained geo-physically as lithological units. From Table 4.5, there are three different kimberliteunits: PK, VK, and HK. HK does have good susceptibility contrast, but not conduc-tivity contrast from the host rock, and hence I ignore HK unit in this interpretation.PK and VK have moderate-high conductivity, and hence all rock units except R0can correspond to them. R2 and R3 can be different kimberlite units, and they havedifferent polarization character: R3 has greater t than R2. However, both R2 andR3 can be either PK or VK and hence distinction of PK and VK is not clear. PK is asubclass of VK and is deposited after an explosive event. Both units can be highlyweathered resulting in significant clay content which is conductive and chargeable.Because of its explosive origin, PK likely has greater pore size than VK; the poremay be filled with water. The pore size of the rock is strongly correlated with thetime constant: greater pore size result in larger time constant [Pelton et al., 1978,Revil et al., 2014]. Based upon that, I interpret R2 (small t) and R3 (large t) asVK and PK, respectively. Here, polarization mechanisms of PK and VK can berelated to both membrane polarization and Maxwell-Wagner polarization. Finalinterpretation is summarized in Table 4.7.Table 4.7: Petrophysical domains built from inversions of airborne EM datasets. Here s , h˜E , and h˜L correspondingly conductivity, and pseudo-chargeability at 130 and 410µs.Rock Unit s h˜E h˜L t InterpretationR0 Low Low Low N/A Host RockR1 Mod.-High Low Low N/A KimberliteR2 Mod.-High High Low Small VKR3 Mod.-High Low High Large PK156Figure 4.45: Comparative sections through the TKC kimberlites. (a) Over-laid anomalies of three physical property models: conductivity (darkgrey), pseudo-chargeability at 130 (red), and 410 µs (green). (b) Gen-erated rock model from three anomalies. The white outlines delineateboundaries of the lake. The green outlines show the extent of DO-27and DO-18 at the surface, based on drilling.1574.3.7 SynthesisExtensive drilling programs have been carried out at TKC, from which four differ-ent kimberlite units have been identified [Doyle et al., 1999, Eggleston and Brise-bois, 2008]:• Dipping sheets of hypabyssal kimberlite (HK) facies limited to the north-eastpart of DO-27.• Shale-rich volcaniclastic kimberlite (VK) crater facies in the northern part ofDO-27.• Resedimented pyroclastic kimberlite (PK) crater facies infilling the core ofDO-27; Diamondiferous unit.• Xenocryst-rich volcaniclastic kimberlite (XVK) crater facies mostly restrictedto the core of DO-18.Fig. 4.46(a) and (b) compares the final petrophysical model obtained frominterpretation of the AEM data sets to the geology based upon drilling results froman extensive drilling program. The agreement is quite good, particularly regardingthe geometric confinement of the pipes. For the DO-27 pipe, interpretation of R2and R3 as respectively VK, and PK, agrees with the ground truth. From the com-parative sections at DO-27 (B-B’), the upper part of the PK and VK units are wellimaged with R2 and R3, respectively. The deeper part of the PK unit is mappedby R1 (obtained from s ), but the bottom boundary is not well distinguished. Thismay due to lack of resolution for the AEM survey for the bottom boundary of aconductor since the EM field significantly decays in the conductor. In addition, ourinterpretation that the DO-18 pipe is VK, is reasonable. XVK is a subunit of VKand was identified through drilling. The HK unit is not distinguished since it doesnot show any anomalous conductivity or chargeability; however, this can be distin-guished by susceptibility which is often available through the magnetic survey thatare usually flown with an AEM survey. Overall, our analysis has clearly demon-strated impact of the obtained IP information from the AEM surveys to characterizevarious kimberlite units at the TKC region.158Figure 4.46: Comparison of 3D geology and rock model. (a) 3D geologicmodel obtained from the known geology and drilling results at theTKC area (b) 3D petrophysical model obtained from airborne EM.The black outlines delineate boundaries of the lake. The green out-lines show the extent of DO-27 and DO-18 at the surface, based ondrilling.1594.3.8 ConclusionsConsistent negative transients have been observed at TKC with various TEM sur-veys: NanoTEM, AeroTEM II, and VTEM. These are due to chargeable materials.Focusing on the most recent VTEM data, four distinct IP anomalies showing dif-ferent decaying features all delineated. The TEM-IP inversion workflow is appliedto the VTEM data. Even the earliest time channel of the VTEM data showed nega-tives and this presented a significant challenge to undertake conductivity inversion.This was overcome by cooperatively inverting the VTEM and DIGHEM data, andgenerated a 3D distribution of sest . dIP[sest ] are obtained by implementing theEM-decoupling procedure. The dIP[sest ] were inverted to recover the 3D pseudo-chargeability at multiple time channels. Then by using representative cells fromthe four chargeable bodies, pseudo-chargeabilities were fitted using a Cole-Colemodel to recover h and t . Four chargeable bodies are imaged at depth, but indifferent time channels: chargeable bodies near A1-A3 are seen at 130 µs; charge-ability near A4 is seen at 410 µs. The estimated t for cells close to A4 is 1160 µsand this is much greater than t (⇠70 µs) for A1-A3. The DO-18 pipe (A1) canthus be differentiated from the southern part of DO-27 (A4), which supports thegeologic hypothesis that the two pipes are distinct intrusive events.Finally, using both conductivity and chargeability information, a 3D rock modelis constructed, which includes four different rocks: R0-R3. Three of them were re-lated to the kimberlites: R1-R3. The moderate conductivity at the upper part of thetwo pipes suggested that this was a highly weathered kimberlitic rock that mightbe PK or VK, but we were not able to differentiate between those units. The ad-dition of IP information however, enabled us to make this distinction because thetwo rocks have significantly different t values. Our final rock model was comparedto the 3D geological model built up from an extensive drilling program; the twomodels are quite similar. No explicit information regarding the known geologyhas been used to constrain both conductivity and IP inversions. Only the generalkimberlite model, the component rock types and their relative physical propertycontrasts were used. Despite that, a rock model is obtained whose major featureswere representative of the geologic model from drilling. The AEM data used inthe analysis were obtained from airborne which are far easier, and less costly, to160collect than ground data.161Chapter 5Grounded Source Example5.1 IntroductionControlled-source electromagnetic (EM) methods excite the earth using either agrounded source (a generator attached to two grounded electrodes) or an inductivesource (currents flowing in a wire loop). For both cases, an electric field will begenerated in the earth, and polarization charges will be accumulated if chargeablematerial exists. Effectively, the rock is electrically polarizable, and this is under-stood by allowing the electrical conductivity to be frequency- or time-dependent[Pelton et al., 1978, Tarasov and Titov, 2013, Revil et al., 2015]. A typical elec-trical induced polarization (EIP; often called DC-IP) survey [Seigel, 1974] usinga grounded source is shown in Fig. 5.1. It consists of grounded electrodes car-rying a current waveform (like the square wave with half-duty cycle shown) andelectrodes to measure voltage differences. When the ground is chargeable the re-ceived voltage looks like that in Fig. 5.2. The decay in the current off-time isthe IP effect. To interpret the observed IP data, a two-stage inversion is usuallydeployed [Oldenburg and Li, 1994]. The first step is to invert late on-time data (V0)using a DC inversion to obtain the background conductivity. The second step is touse the obtained conductivity to generate a sensitivity function, and then invert lateoff-time data (Vs); this is often called DC-IP inversion.Although application of this method has been successful, especially for miningapplications [Fink et al., 1990, Oldenburg et al., 1997, Wait and Gruszka, 1986], a162main concern is the second step. The time-decaying fields are assumed to be purelythe result of IP phenomena and any EM induction effects in the data are ignored.This assumption can be violated when the earth has a significant conductivity andEM effects remain stable and flat even in the late off-time [Veeken et al., 2009a,b]Removing EM induction responses from the measured data is referred to as EM-decoupling and it has been a focus of attention for many years. Most analyseshave used simple earth structures, a half-space or a layered earth, to ameliorate itseffects [Routh and Oldenburg, 2001, Wynn and Zonge, 1975, Kratzer and Macnae,2012]. However, with recent capability to handle 3D TEM forward modellingand inversion [Commer et al., 2011, Haber, 2014, Oldenburg et al., 2013], I canuse more complicated conductivity models. In addition, full 3D TEM forwardmodelling using a complex conductivity model is currently available. Results canbe obtained directly in the time-domain [Marchant et al., 2014] or by transformingfrequency domain responses to the time-domain [Hohmann and Newman, 1990,Flis et al., 1989]. With these advances I can revisit the challenge of decoupling EMeffects and extracting better information about conductivity and chargeability.A TEM-IP workflow illustrated in Section 1.6 has successfully been appliedto airborne IP data using inductive sources (Chapter 4). The workflow includesthree main steps: 1) inverting early-time TEM data (those not affected by IP) torecover a 3D conductivity model, 2) EM-decoupling (forward modelling the EMresponse and then subtracting it from the observations), and 3) IP inversion torecover pseudo-chargeability distribution at each time channel. The problem ofinverting IP data using grounded sources follows the same workflow, but someaspects are greatly simplified because DC-IP measures data when electric fields,and charge accumulations, have reached a steady state. This provides another dataset from which information about the electrical conductivity can be extracted. Amajor difference between a conventional DC-IP inversion and our approach is theuse made of early-time channels in the DC-IP data. In conventional work thesehave been considered as “noise” being corrupted by EM coupling effects and hencethrown away. However, I consider these as “signal” to recover conductivity. In thisstudy, the TEM-IP inversion workflow is applied to a synthetic example using agradient array with a grounded source. I illustrate the steps in the workflow listedabove but the first step is altered so that one can optionally invert the DC data163Figure 5.1: Conceptual diagram of a ground-based grounded source withhalf-duty cycle current waveform.Figure 5.2: A typical example of overvoltage effects in electric IP data.and use that for decoupling. This allows me to compare the effects of differentestimates of the conductivity model on the final chargeability models.164Figure 5.3: Observed voltage with EM induction effects. EM effects domi-nate the early on and off-time data. Inset figure shows the enlarged off-time voltage and it demonstrates early EM induction effects (negative)and late-time IP effects (positive). Dashed and solid lines distinguishpositive and negative voltages.1655.2 TEM-IP Inversion Workflow for Grounded SourceThe TEM-IP inversion workflow has applied to airborne EM examples in Chapter4, and it is applicable here for the grounded sources, but there are some differ-ences to be taken care of. The workflow has three steps. The first step is to invertthe TEM data to recover the 3D conductivity model. As in the inductive sourcework, I use only early-time data, which are not IP-contaminated. Note that theseearly-time data have previously been considered as “noise in conventional DC-IPanalyses and hence have been thrown away. However, here I consider these as “sig-nal” and use them to recover a better conductivity model. Another possibility forobtaining a background conductivity is to use the steady-state fields (stable plateaulevel) just prior to switching the current off. These are the potentials that are tra-ditionally used in DC-IP inversion. Inversion of these data yields a conductivitythat is s0 = s•(1h) (when on-time pulse is long enough to be fully charged up)but if h is small then this will be a reasonable approximation to s•. The inversionof DC data is analogous to inverting only one frequency in a frequency-domaindata set. Hence it might be expected that inverting data at multi-times (equivalentto multi-frequencies) would produce a better result. Our experience verifies this.Nevertheless, the DC fields are valuable and I wish to use them. The options areto invert the DC and TEM data together, or treat them as two separate data sets.For the present I have chosen the latter since I then do not have to contend with theissue that the DC fields are really s0, but TEM data at early off-time are close tos•. To investigate how much information can be extracted from each data set, hereI separately invert DC and TEM data and recover both DC and EM conductivities.The second step of the workflow is EM-decoupling. The estimated conduc-tivity model, sest , from step 1, is used to generate estimated IP data, dIP[sest ],according todIP[sest ] = dF [sest ] (5.1)where F [·]Maxwell’s operator taking 3D conductivity as an argument, and F [sest ]is the estimated fundamental data. Note that F [sest ] might be different from F [s•]because sest is not identical to s•. The differences will be most pronounced atearly times when induction responses are large. Thus potential errors in dIP[sest ]will be largest at early times. EM-decoupling will be effective in the intermediate166times when both EM and IP signals are considerable; at late times (IP-dominant),EM-decoupling may not be required.The final step in the process is to carry out the IP inversion. I adopt the conven-tional IP inversion approach [Oldenburg and Li, 1994], which uses a linear form ofIP responses written asdIPi = Jh˜ i (5.2)where J is the sensitivity matrix, and h˜ i is a column vector for the pseudo-chargeabilityat the i-th time channel. The conductivity model sest is required to generate thesensitivity matrix. From eq. ( 2.23), pseudo-chargeability can be defined ash˜(t) =Z t•4s(tu)s•we(u)du (5.3)where4s(t) =L 1[4s(s)], andL 1[·] indicates the inverse Laplace transform;we(t) is the time history of the electric field. In particular, for the EIP case with thestep-off current, I let we(t)⇡ 1uon(t), where uon(t) is a Heaviside step function.For instance, with the Cole-Cole model shown in eq. ( 1.7), integration in eq. ( 5.3)can simply be evaluated when c=1,h˜(t) = he1(1h)t t (5.4)Here (1h) in the power of the denominator arises from the conductivity formu-lation of Pelton’s Cole-Cole model [Pelton et al., 1978], whereas it is absent in theresistivity formulation. I invert each time channel of IP data separately, and re-cover pseudo-chargeability at multiple times. The recovered pseudo-chargeabilitywill be a 4D property: h˜(x,y,z; t). Interpreting this recovered pseudo-chargeabilityto extract intrinsic IP information such as h , t , and c is possible [Yuval and Old-enburg, 1997, Ho¨rdt et al., 2006, Fiandaca et al., 2012], but I do not attempt that inthis grounded source example.1675.3 DataAs an example, I use a grounded source and multiple receivers which measurevoltages as shown in Figure 5.4; this is often called a gradient array. Four blocks(A1-A4) presented in Figure 5.4 have different s• and h values (see Table 5.1); allblocks have t=0.5 s and c=1 (Debye model). Only blocks A2 and A3 are charge-able. The length of the source wire is 6 km and the potential differences betweentwo electrodes along easting lines are measured at 625 locations (200-m lengthfor the potential electrodes). Simulations were performed using the EMTDIP code[Marchant et al., 2014] with a step-off waveform. The initial condition of thissimulation corresponds to solving a DC problem for on-time data using a DC con-ductivity, s0, then time-stepping proceeds to compute TEM responses at off-times.Two separate simulations are performed calculating d = F [s(s)] and dF = F [s•].For the latter simulation h was set to zero everywhere. Note that s(s) = s• whenh = 0.Fig. 5.5 shows the observed DC voltage during the on-time and its conver-sion to apparent resistivity. Lower apparent resistivity anomalies are reflective ofthe A1 and A3 bodies. After current switch-off, voltages are measured at sixtylogarithmic-based time channels ranging from 1-600 ms. Computed responses at5, 80, 130, and 500 ms are shown in Figure 5.6. To investigate EM effects atdifferent times, I present the time decaying curves of d and dF at A1-A4 in Figure5.7. The four vertical black lines correspond to the times for the data maps shownin Figure 5.6.At 5 ms EM induction effects are dominant and all data are negative. At 80 msboth EM and IP effects are considerable but still, all data are negative. At 130 ms asign reversal has occurred at A2 and A3 bodies resulting in positive data. A1 alsohas smaller positive values. Note that A2 and A3 are chargeable but A1, which isconductive, is not. Therefore it is difficult to differentiate chargeability and con-ductivity anomalies just by looking at observed data at 80 and 130 ms. At 500 mshowever, EM induction effects here significantly decayed and IP signals are domi-nant; only A2 and A3 show positive anomalies. Thus, depending on the measuredtime window, background conductivity, and IP parameters of chargeable bodies,data could be dominated by IP effects or not. Hence, whenever our measured time168window is not late enough to be considered as IP-dominant time, EM-decouplingis a crucial step. Note that the A1 anomaly at 80 and 130 ms could have beenmisinterpreted as a chargeable response.Table 5.1: Conductivity (s•) and resistivity (r•) at infinite frequency, andCole-Cole chargeability (h) values for five units: A1-A4 and half-space.Division A1 A2 A3 A4 Half-spaces• (S/m) 1 0.01 0.1 0.001 0.01r• (Wm) 1 100 10 1000 100h 0 0.1 0.1 0 0Figure 5.4: Plan and section views of the 3D mesh. Black solid lines showthe boundaries of four blocks (A1-A4). Only A2 and A3 (red labels) arechargeable. Arrows indicate a wire path for the grounded source. Blackdots indicate potential electrodes.169Figure 5.5: Plan maps of the observed DC data: (a) Voltage and (b) Apparentresistivity.Figure 5.6: Plan maps of the observed TEM data at (a) 5 ms, (b) 80 ms (c)130, and (d) 500 ms. Dashed and solid contours differentiate negativeand positive data.170Figure 5.7: Time decaying curves of the observed (d) and fundamental (dF )data. Four sounding locations close to (a) A1, (b) A2, (c) A3, and (d)A4 are presented. Blue and black color indicates d and dF .1715.4 Conductivity InversionTo recover s•(x,y,z), the first six channels of the TEM data (1-6 ms) are used;these have only minor contamination from IP. In addition, DC data which containIP effects, but have minor EM induction effects can be used. I first invert the DCdata and recover the DC conductivity, sDCest . Note that this theoretically representss0 rather than s•. For DC inversion, the DC-IP package of SIMPEG is used[Cockett et al., 2015]. The gradient array data used here have only a single sourceand hence they can be represented as a potential field and have no inherent depthresolution. Without a depth weighting the resultant conductivity from the inversioncan be a thin layer at the surface. To overcome this, I incorporate a depth weightingsimilar to that used in magnetic inversion [Li and Oldenburg, 1996]. The TEM dataare inverted with our H3DTD code [Oldenburg et al., 2013]. Parameters for bothinversions are summarized in Table 5.2 and further details are given in Section3.2. The recovered conductivity models from the 3D TEM and DC inversions arerespectively shown in Figure 5.8(a) and (b). The conductive blocks A1 and A3are better imaged with the TEM inversion. Note the TEM data have inherent depthresolution and hence no depth weighting was required.Table 5.2: Parameters for the TEM and DC inversions. See Section 3.2 forexplanation of parameters. shal f indicates conductivity of the homoge-neous half-space, which has a value of 0.01 S/m.Parameters TEM inversion DC inversionax=ay=ax 1 1as 104 104f ⇤d 625 ⇥ 6 = 3750Uncertainty(e j) 0.02|dobsj | 0.05|dobsj |+0.01m0 log(shal f ) log(shal f )mre f log(shal f ) log(shal f )Model weighting N/A (z z0)3Bounds constraint m> 0 N/A172Figure 5.8: Recovered conductivity models from: (a) TEM and (b) DC inver-sions by inverting off-time (EM) and on-time data (DC), respectively.1735.5 EM-decouplingThe next step is EM-decoupling. Eq. ( 5.1) is implemented by using sest from theTEM inversion (Fig. 5.8). In Fig. 5.9, I present d, F [sEMest ] and dIP[sEMest ] data at80 ms. At this time, both EM and IP effects are considerable. Our EM-decouplingprocedure removes EM effects due to conductivity especially for regions close toA1 (not chargeable) and A3 (chargeable). Removing the conductive anomaly at A1is especially important, because this could have been misinterpreted as chargeableanomaly.A crucial aspect of our EM-decoupling procedure is the effect of the back-ground conductivity. To show this I consider three other candidates, namely trues•, sDCest , and half-space conductivity, shal fest . I compare the performance of theEM-decoupling procedure for all four conductivity models. Fig. 5.10 shows thetrue and estimated fundamental responses. sEMest does the best job of predicting thefundamental response followed by sDCest . The halfspace response is very smoothand is a poor approximation to the fundamental response.The discrepancies between the true and estimated fundamental responses aredirectly observed when the true and estimated dIP data are compared as in Fig.5.11. dIP[sEMest ] are the closest representation of dIP, followed by dIP[sDCest ]. Notethat the latter still has some positive response over A1. As expected the dIP[shal fest ]results in a poor representation of the IP response. It indicates that A1 is charge-able, and that the data at A3 are higher than those at A2. If these data are input toa 3D IP inversion, they produce strong artifacts from which incorrect conclusionscan be drawn.174Figure 5.9: Plan maps for 80 ms. (a) observed data, (b) estimated fundamen-tal using sEMest , and (c) IP.Figure 5.10: True and estimated fundamental responses at 80 ms. (a) Truefundamental response; (b), (c), and (d) correspondingly indicate esti-mated fundamental responses using sEMest , sDCest , and shal fest .Figure 5.11: True and estimated IP responses at 80 ms. (a) True dIP; (b), (c),and (d) correspondingly indicate estimated IP responses using sEMest ,sDCest , and shal fest .1755.6 IP InversionTo recover 3D pseudo-chargeability, I invert the IP data sets at 80 ms obtained us-ing the two estimated conductivity models: sEMest and sDCest . These conductivities areused to generate the linearized sensitivities as outlined in Section 2.2. The unit ofthe pseudo-chargeability is basically dimensionless, but I use mV/V. The linear sys-tem is inverted with the added constraint of positivity on the pseudo-chargeability[Oldenburg and Li, 1994]. Inversions are done using a single time channel and asingle source so there is no depth resolution for the IP inversion. Thus, similar tothe previous DC inversion, a depth weighting needs to be incorporated. For the IPinversion the reference model, mre f is set to zero, and the uncertainty is set to 5%of the maximum amplitude of IP data (0.05 max(dIP)) at 80 ms. The recovered3D pseudo-chargeability models from dIP[sEMest ] and dIP[sDCest ] are shown in Fig-ure 5.12. The two true chargeable bodies, A2 and A3, are well imaged althoughthe pseudo-chargeability model from dIP[sEMest ] is of better quality than that fromdIP[sDCest ]: a) it has fewer artifacts near A1 as expected from dIP[sDCest ] and b) therecovered chargeable blocks are more compact. In addition, considering the truepseudo-chargeability value of both A2 and A3 is 83 mV (obtained from the eq.5.4), the recovered h˜ values of A2 and A3 blocks (⇠30 mV/V) are in a reasonablerange.176Figure 5.12: Recovered chargeability models from IP inversions. IP datacomputed using an estimated conductivity model from (a) TEM and(b) DC inversions, respectively.1775.7 DiscussionConductivity inversion is a crucial part of the TEM-IP inversion workflow, and forthe time-domain DC-IP case there are two types of data: a) DC (late on-time) andb) TEM (early off-time). The DC data can be inverted directly to give an estimate ofs0 = s•(1h). Often h is small (⌧1) hence s0 ⇡ s•. In the example presentedhere using a gradient array however, it was required to incorporate a depth weight-ing. With that, the recovered conductivity was sufficiently good to allow a goodEM-decoupling and subsequent inversion of the IP data. The recovered charge-ability nicely delineated the two chargeable bodies but a weak chargeable artefactwas observed at the location of conductor. In the synthetic example the early-timeTEM data produced a better conductivity model, without the need for any depthweighting. This illustrated the resolving power that is inherent in working withmultiple time channels. The pseudo-chargeability recovered from the inversionwas improved compared to that using the DC data. Notwithstanding that result,the challenge with inverting the early-time TEM data is the need to work with datathat do not exhibit significant IP effects. IP contamination in the data will bias re-covered conductivities towards increased resistivity and thereby generate artefacts.Deciding how to recognize when data are IP-coupled is an area for future research.So too is the need to invert both the DC and early-time TEM data simultaneously.Both issues need to be tested for a more realistic model (e.g. complicated back-ground model) and different arrays of DC-IP surveys (e.g. pole-dipole) to be usedin practice.The accuracy to which s• needs to be estimated depends upon the informationthat is ultimately desired. If the goal is for detection of a body then a less re-solved conductivity may suffice. As time increases, the dIP data obtained throughEM-decoupling is less adversely affected by inaccuracies in sest . For example,as shown in Fig. 5.13, performance of the EM-decoupling with the sDCest at 130ms was good. At 80 ms (Fig. 5.11) the performance is still generally good butartefacts are beginning to be more prominent (for example around A1). For thepurpose of detection, both of these results may suffice. But, although the detectionof a chargeable body can be accomplished with a moderately good approximationto the true conductivity, the need to obtain a high quality conductivity becomes178more important when spectral information is sought. This may require qualityIP data at early times. Fig. 5.14 shows time decaying curves of the observed,estimated fundamental and IP data at A3 (conductive and chargeable). dIP[sDCest ]converges to dIP[s•] after 130 ms. However, dIP[sEMest ] converges to dIP[s•] after20 ms resulting in a broader time band of the IP data compared to dIP[sDCest ]. Thiswill be a crucial factor for recovering intrinsic IP parameters because capturingIP information in a greater time band will reduce non-uniqueness in the inverseproblem for multiple parameters such as h(x,y,z), t(x,y,z), and c(x,y,z) [Lajaunieet al., 2016].Figure 5.13: True and estimated IP responses at 130 ms. (a) True dIP; (b),(c), and (d) correspondingly indicate estimated IP responses with useof EM, DC, and half-space conductivity.Figure 5.14: Decay curves of the observed, fundamental and IP data at A3.(a) True and estimated fundamental responses with the observed data.(b) True and IP responses with the observed data.1795.8 ConclusionsIn this chapter, the TEM-IP inversion workflow has been applied to a syntheticgrounded source TEM example. First, I inverted on-time DC data and recovereda 3D conductivity. Then, I inverted six of the earliest time channels of TEM data,which have minor IP-contamination, and recovered an another estiamte of the 3Dconductivity. These early TEM data have often been thrown away because theyare considered as “noise”. However, by considering them as “signal” and invertingthem, I recovered a better conductivity model in the sense of depth resolution. Sec-ond, the recovered conductivity, sEMest was used in our EM-decoupling procedureto generate IP data. The procedure was effective for removing EM induction inthe observations, especially for regions close to A1 and A3, which had significantconductivity responses. Third, I inverted the IP data set generated from the TEMconductivity model using conventional 3D IP inversion. The recovered pseudo-chargeability successfully imaged two true chargeable anomalies. This demon-strates that our TEM-IP inversion workflow can be effective for recovering a goodestimate of electrical conductivity, for removing EM-coupling from IP data, andfor obtaining a 3D distribution of pseudo-chargeability.180Chapter 6ConclusionsEarth materials are chargeable, and hence EM data can include IP effects. The lat-est developments of instrumentation and processing techniques provide high qual-ity EM data in which small IP signals can be observed. For instance, measuringnegative transients in AEM surveys are now commonly observed and regarded assignal whereas that was not the case a decade or more ago. Therefore, developingan effective methodology to handle EM data that includes both EM and IP signals,is timely. My goal is to develop physical understanding and computational pro-cedures to recover 3D conductivity and chargeability from EM data, and to applythe technique to field data sets. This requires separation of the EM and IP por-tions in the EM data, and thus contending with important aspects of EM-couplingis another essential task. I have designed this thesis to contribute to the scientificcommunity by developing a workflow, which extracts 3D conductivity and charge-ability from time-domain EM data. The procedure has been applied to syntheticand field examples. For the development of the workflow, and its application, Ifocus on answering three important and unsolved questions:1. How can we understand the fundamentals of inductive source IP, and can thetime-domain airborne IP data be linearized similar to DC-IP data?2. Can distributed IP information be recovered from the airborne data showingIP effects?3. Can EM effects in DC-IP data be effectively removed to generate high qual-181ity IP data?The following sections summarize my work on these questions. Of centralimportance is the development of the workflow to handle TEM-IP data.6.1 Development of the TEM-IP Inversion WorkflowThe TEM-IP inversion workflow includes three main steps: a) Conductivity inver-sion, b) EM-decoupling, and c) IP inversion. In the first step, TEM data at early-time channels (where IP effects are minor) are used to recover an estimated 3Dconductivity model, sest . In the second step, estimated fundamental data, F [sest ],are simulated and subtracted from the observations to generate IP data, dIP; thisis the EM-decoupling procedure. In the third step, the obtained IP data are lin-earized as a function of the pseudo-chargeability: dIP = Jh˜ . Then by using thatrelationship, the IP data can be inverted to generate a 3D distribution of pseudo-chargeability. This inversion can be done for multiple time channels. Lastly thepseudo-chargeability of a cell at different times can be used as data in a parametricmodel for estimating parameters to characterize the chargeability. Here I estimateCole-Cole parameters of chargeability and time constant. Linearizing IP data forinductive sources is more challenging than for grounded sources, because of thedynamic polarization buildup that arises from nature of the diffusive electric fields.This polarization character needs to be captured to develop a linear IP function.Chapter 2 is related to the first question stated at the beginning of this conclu-sion section. I present linearization of inductive source IP (ISIP) data. A core ideaof linearizing ISIP data is based upon the time behavior of the electric field from aninductive source. Consider a single pixel in the earth. When the current is turnedoff, the electric field at the pixel increases to a peak and then decays. I assumed thatthe major polarization effect is characterized by the electric field at the time whenit reaches the maximum. Based upon that hypothesis, the IP current is expressedas a function of the pseudo-chargeability. This IP current captures both vortex andgalvanic currents which contribute to the IP signals. With the application of theBiot-Savart law to the IP current, both the magnetic field and its time derivativeare obtained. This provides a desired linear form: dIP = Jh˜ for ISIP data. The as-sumptions made to linearize the IP function are thoroughly tested with numerical182experiments. Particular emphasis is placed upon the airborne IP case.Chapter 3 constructs a foundation to answer the second question. By using thelinearized IP function, I develop a 3D IP inversion algorithm, and test that with asynthetic airborne IP data. Depth weighting is introduced to compensate for a lackof intrinsic depth resolution in the airborne IP data. With this, the 3D IP inversionrecovers a reasonable geometric shape and location of the chargeable body. Forthe inversion I assumed that the true s•(x,y,z) is available. However, in practicesest obtained in the first step is always different from s•. This incorrect conductiv-ity has two effects on the IP inversion. First it generates errors in the dIP[sest ]. Apositivity constraint imposed on the pseudo-chargeability greatly ameliorates theseerrors particularly where errors are positive. The other avenue by which an incor-rect conductivity can affect the inversion is through the sensitivity matrix, which isideally J[s•]. However the impact of an incorrect conductivity in the sensitivity ismuch less than the errors obtained by an incorrect fundamental response F [sest ]. Asynthetic example is introduced and each time channel of the dIP is inverted, andpseudo-chargeability at multiple times are obtained. The pseudo-chargeabilities atmultiple times can be inverted to obtain Cole-Cole parameters h , t . I have foundthat the estimated t was close to the true value, whereas h was poorly estimated.6.2 Application to Airborne IP DataChapter 4 addresses the second question by applying the TEM-IP inversion work-flow to airborne EM (AEM) data. This chapter is composed of three sections,which present a synthetic AEM example, and two field airborne examples overmineral deposits: the Mt Milligan porphyry deposit and the Tli Kwi Cho kimber-lites.In Section 4.1, I generate the synthetic AEM data, and apply the workflow.Each of three steps in the workflow is carefully tested. For the conductivity inver-sion, early-time channels which do not have negative values are used to recovera 3D conductivity, sest . Overall, the conductivity structures are imaged well. Inthe second step, dIP[sest ] is obtained by subtracting F [sest ] from the observations.This EM-decoupling shows good performance in the intermediate and late times.However, the EM-decoupling shows poor performance in the early times, which183is expected. The third step is the IP inversion. By inverting each time channelof the obtained dIP[sest ], a pseudo-chargeability at multiple times is recovered.This demonstrates that an IP inversion can provide location and geometric shape ofchargeable bodies embedded in the subsurface as long as appropriate depth weight-ing is applied. By interpreting the recovered pseudo-chargeability, Cole-Cole pa-rameters, h and t , are extracted at four selected cells, that corresponded spatiallyto the four chargeable blocks. Recovery of t is reasonable, whereas h is mostlyunderestimated. From this result, I conclude that there is the potential to extractintrinsic Cole-Cole parameters from airborne IP data but the time constant seemsto be the most robust parameter to be extracted.Section 4.2 presents application of the workflow to field VTEM data over theMt Milligan porphyry deposit. Clear negative anomalies are shown in the VTEMdata. By applying the workflow, a 3D distribution pseudo-chargeability is recov-ered. Negative transients are ignored for the conductivity inversion. The obtained3D conductivity, sest , is used in an EM-decoupling procedure, and dIP[sest ] ateight time channels are obtained. Two IP anomalies, which are not visible as neg-ative voltages in the observation, are detected around the Rainbow fault. Usingthe IP inversion, each time channel of the dIP[sest ] is separately inverted, and a3D pseudo-chargeability at each of the eight time channels is recovered. Highlychargeable volumes are successfully imaged on the edges of three mineralizedzones; they are interpreted as highly altered zones. Using both conductivity andpseudo-chargeability models, 3D rock models are constructed. Importantly, highlyaltered parts of the monzonite stock is recognized from the pseudo-chargeability,and this may have potential to help characterize a porphyry deposit. I believe thatthis was the first time that a voxel inversion has been used to recover a 3D charge-ability from AEM data.In Section 4.3, the workflow is applied to the VTEM data over the Tli Kwi Chokimberlites where consistent negative transients have been observed with variousTEM surveys. A main question here is whether observed negatives are geologicalnoise or signal that can be useful for kimberlite exploration. I focus on the lat-est VTEM data, which shows four distinct IP anomalies, and apply the workflow.Even at the earliest time channel, some of the VTEM data are negatives, and thispresents challenges for the conductivity inversion. This is overcome by coopera-184tively inverting the VTEM and DIGHEM data, and a 3D distribution conductivityis recovered. Two conductive pipes are respectively imaged at DO-18 and DO-27pipes; the conductive pipe at DO-27 is embedded at depth, whereas that at DO-18extends to the surface. The estimated conductivity is used for EM-decoupling pro-cedure resulting in dIP[sest ] at multiple time channels. dIP[sest ] at 130 µs showsthree IP anomalies at DO-18 and in the north-eastern part of the DO-27. At 410µs, an IP anomaly at the southern part of the DO-27 is revealed. Hence four IPanomalies are identified in total. By inverting the obtained dIP[sest ], 3D pseudo-chargeabilities at multiple time channels are recovered. Four chargeable bodies areimaged at different time channels. By using the recovered pseudo-chargeabilityat multiple times, h and t values of representative cells for the four chargeablebodies are extracted. The recovered t from the chargeable body at southern part ofthe DO-27 is greater than that from the other three chargeable bodies showing thatthere is at least two distinct rock units. The pseudo-chargeability at 130 µs and 410µs along with the conductivity are used to generate a 3D rock model. Three differ-ent rock units related to kimberlites are identified. Particularly, two rock units areinterpreted as PK and VK units. Distinction between PK and VK is obtained fromthe different polarization character considering their pore size: PK is characterizedby large pore size and hence large t , and VK is characterized by small pore sizeand hence small t . The PK unit is the diamondiferous unit. In conclusion, the 3Drock model that is obtained from the interpretation of the AEM data shows majorfeatures that are a representation of the geologic model obtained from drillings. Ibelieve this demonstrates how the obtained polarizable information could help akimberlite exploration.6.3 Application to DC-IP DataChapter 5 deals with the third question. The workflow is applied to a syntheticgrounded source example. Emphasis is on EM-decoupling. This requires an es-timate of the background conductivity. This is first obtained by inverting the DCdata to generate sDCest . Then several early-time channels of TEM data are invertedwith a 3D TEM inversion is used to recover sEMest . These early TEM data haveoften been discarded because they are considered as “noise”. However, by consid-185ering them as “signal” and inverting them, a better conductivity model than thatfrom the DC data is recovered The sEMest is used in the EM-decoupling procedure.This was effective for removing EM induction in the observations, especially forthe conductive blocks embedded in the earth. The resultant IP data are invertedto recover a 3D pseudo-chargeability. Chargeable blocks are successfully imaged.These results demonstrate that the TEM-IP workflow can be effective for recov-ering a good estimate of electrical conductivity, for removing EM signals from IPdata, and for obtaining a 3D distribution of pseudo-chargeability from groundedsource TEM data.6.4 Future ResearchNotwithstanding the successful recovery of a 3D chargeability from synthetic andfield TEM data sets using the developed workflow, there are some potential weak-nesses in the workflow, and they require further investigations. Hence, I will elab-orate on them and also suggest potential approaches for improvement.A crucial challenge arises from an incorrect conductivity obtained in the firststep, and its impact on subsequent EM-decoupling. Let4d = F [s•]F [sest ]. Theobtained dIP[sest ] always includes residuals,4d as well as other additive noise.dIP[sest ] = dIPtrue+4d+noise. (6.1)I suggested applying a regional removal technique to reduce the effect of4d, andpresented a synthetic example. Basically I selected some regions of the data whichwere felt to have minimal IP effects. So dIP[sest ] '4d. I fitted those points witha polynomial and used the surface as an estimate of4d. This was subtracted fromdIP[sest ]. This procedure is similar to the removal of a regional field in gravity andmagnetics. It suffers form the same limitations but as in potential field inversions.However it may be an effective way to handle an unwanted background signal.An alternative and less subjective approach was to put a positivity constraint onh˜ when inverting dIP[sest ] to prevent positive bias, since negative h˜ usually gen-erate postive IP data for the coincident-loop airborne system. I used the positivityconstraint in the field examples to prevent fitting positive4d, so negative4d canstill make artefacts in a recovered pseudo-chargeability model. When that is the186case, developing an effective regional removal technique to obtain clean IP datamight be important, and it is worthy of further investigation.For the conductivity inversion, data which have minimal IP effects are required.I choose several early-time channels, which do not have negative values. Depend-ing upon how one selects time channels to obtain EM-dominant data, the quality ofan estimated conductivity can vary. A more robust approach needs to be developed.Three options are considered. The first option is using the on-time data [Smithand Klein, 1996] because the on-time data will be dominated by the primary fieldsand EM induction effects. This is the simplest and effective solution, but in realityobtaining the high quality on-time data for a time-domain AEM system is chal-lenging due to the inaccuarate measure of the primary field. The second option isto modify the TEM-IP inversion workflow. To improve the quality of the conduc-tivity inversion, an additional step is required for the original workflow to removeIP effects in observations when inverting for a conductivity. Accordingly, the orig-inal workflow can be altered as (a) Conductivity inversion, (b) EM-decoupling, (c)IP inversion, and (d) IP-decoupling. (a) and (b) can be same as the original work-flow, but (c) needs to be modified to estimate dIP for subsequent IP-decouplingprocedure. For this, rather than separately inverting for each time channel of dIPto recover a pseudo-chargeability at multiple time channels, I suggest parameter-izing the pseudo-chargeability with a representative function such as Cole-Coleor Stretched exponential [Kohlrausch, 1854, Pelton et al., 1978, Tarasov and Titov,2013], and inverting all time channels of dIP together to recover 3D distributions ofthe parameters; this can be considered as a spectral IP inversion similar to Fiandacaet al. [2013]. For instance, considering a Cole-Cole model for the parameterization,h(x,y,z), t(x,y,z), and c(x,y,z) can be recovered and then dIP can be estimated.The estimated dIP can be subtracted from observations to yield cleaner EM data;this is an IP-decoupling step. The four steps of the workflow can be iterated untilit reaches to a stopping criteria. This iterative workflow still require the approachto select EM-dominant data when inverting for conductivity, but the impact of thisselection is much less than the original workflow since the quality of EM data willbe enhanced in the iterative procedure by the IP-decoupling.The thrid option is inverting both EM and IP signals together. As mentionedpreviously there are various ways to parameterize a complex conductivity. The187TEM data with a Cole-Cole representation can be written asd[s•,h ,t,c]. (6.2)Four Cole-Cole parameters are defined in 3D space, and hence this increases thenumber of unknowns compared to the inversion used in the workflow; this in-creases non-uniqueness of the inversion. Moreover, the problem is now non-linearrather than linear. This makes the solution more complicated to achieve. Lastly,evaluation of the convolution between time-dependent conductivity and the elec-tric field is required for each transmitter to compute forward responses, and thisdramatically increases the computational cost of the inversion. However, the pos-sibility exists to make this option robust and computationally efficient, and it isworthwhile to pursue. For the non-uniqueness issue, rather than starting from ahomogeneous conductivity, one can invert early-time channels of TEM data to re-cover a 3D conductivity, then use the recovered conductivity as a starting modelfor the non-linear inversion.When linearizing IP response from an inductive source (Section 2.2), I as-sumed a step-off current waveform, and based upon that an electric field at max-imum was selected as the reference electric field to compute the sensitivity func-tion. However in practice, TEM surveys use various types of current waveforms.Although the linearization steps suggested in the thesis are general enough to han-dle this, more complicated waveforms still needs to be taken into account, andtested; this can be important especially when quantitative polarization informationis desired.Finally, the workflow has been applied to the synthetic grounded source exam-ple, but not to field examples. Therefore, the workflow must be applied to fieldexample to demonstrate its capability.6.5 Concluding commentsBefore I finish this thesis, I want to remind readers about why extracting polariza-tion information of the subsurface is important. Most of earth rocks are chargeable,and they have different polarization characteristics. This polarization informationwill be valuable information to characterize the subsurface in various applications:188mining, oil and gas, groundwater, and geotechnical problems. In addition, IP sig-nals, which have often been acquired from the ground using DC-IP surveys, cannow readily be obtained from the air using AEM surveys with the latest improve-ments in instrumentation, although AEM surveys has limited depth of investigationfor chargeable bodies compared to DC-IP surveys. In this perspective, the IP tech-nique is just now starting to be used, and I believe there will be more areas that IPcan make impact. I illustrate a few important areas below.First is airborne IP. This allows polarization information of a broad region tobe quickly mapped. It is still controversial about whether measured IP data are ge-ological noise or signals that can be useful for mineral explorations. Nonetheless,there are IP signals in AEM data, and hence some polarization information canbe extracted, and this may be useful for various applications. For instance in thethesis, I showed an example when airborne IP can be used to distinguish differentkimberlites in the ground. In addition, clays and permafrost are also chargeable,and hence they can be delineated from the airborne IP data. A second applica-tion is estimating hydraulic permeability of rocks using grounded IP surveys. Anumber of lab-scale and some field-scale research, which estimate permeability ofrocks from IP data, have been performed and showed promising results [Slater andLesmes, 2002, Revil and Florsch, 2010, Weller et al., 2010, Ho¨rdt et al., 2007].This can make an impact on groundwater studies, which uses hydraulic permeabil-ity as an input for their simulations. A third application is characterizing landfillor waste dumps. Some contaminated materials in the landfills can be chargeable,and hence grounded IP surveys can be used to characterize them, particularly forold landfill areas, which do not have sealing facility to prevent leaking leachate[Gazoty et al., 2012]. Further, depending upon the situation, contaminated groundneeds to be remediated, and for this bioremediation can be an option. If this re-mediation process removes chargeable contaminants, then the chargeability of theground may decrease, and hence these changes can be monitored by IP.To conclude, obtaining polarization information of the subsurface is an impor-tant task for geoscience applications, and further there will be more applicationsin which IP can play a crucial role. The workflow, and its application, presentedin this thesis provide a general framework that can extract polarization informationfrom EM data. The results obtained in this thesis show that important information189about chargeability can be obtained from EM data and this opens the scope forpotential new applications.190BibliographyW. Boyko, N. R. Paterson, and K. Kwan. AeroTEM: system characteristics andfield results. The Leading Edge, 20(10):1130–1138, 2001. ! pages 136M. Bu¨cker and A. Ho¨rdt. Long and short narrow pore models for membranepolarization. Geophysics, 78(6):E299–E314, 2013. ! pages 11Y. Chen and D. Or. Geometrical factors and interfacial processes affectingcomplex dielectric permittivity of partially saturated porous media. WaterResources Research, 42(6):1–9, 2006. ! pages xiii, 11, 12, 14R. Cockett, S. Kang, L. J. Heagy, A. Pidlisecky, and D. W. Oldenburg. SimPEG:An open source framework for simulation and gradient based parameterestimation in geophysical applications. Computers & Geosciences, 85, Part A:142–154, 2015. ! pages v, 81, 172, 208K. S. Cole and R. H. Cole. Dispersion and Absorption in Dielectrics I.Alternating Current Characteristics. The Journal of Chemical Physics, 9(4),1941. ! pages 15, 88M. Commer, G. A. Newman, K. H. Williams, and S. S. Hubbard. 3Dinduced-polarization data inversion for complex resistivity. GEOPHYSICS, 76(3):F157–F171, 2011. ! pages 32, 163D. Cowan. personal communication, 2015. ! pages xiii, 18S. G. R. Devriese, K. Davis, and D. W. Oldenburg. Inversion of airbornegeophysics over the DO-27/DO-18 kimberlites Part 1: Potential fields.Interpretation, 5(3):T299–T311, 2017. ! pages xxiii, 138B. J. Doyle, K. Kivi, and B. H. S. Smith. The Tli Kwi Cho (DO27 and DO18)Diamondiferous Kimberlite Complex, Northwest Territories, Canada. InProceedings of the VIIth International Kimberlite Conference, pages 194–204,1999. ! pages 158191T. Eggleston and K. Brisebois. DO-27 Diamond Project, Northwest Territories,Canada. Ni 43-101 report, AMEC, aug 2008. ! pages 158H. ElKaliouby and E. Eldiwany. Transient electromagnetic responses of 3Dpolarizable body. GEOPHYSICS, 69(2):426–430, 2004. ! pages 136R. J. Enkin, D. Cowan, J. Tigner, A. Severide, D. Gilmour, A. Tkachyk,M. Kilduff, and J. Baker. Physical property measurements at the GSCpaleomagnetism and petrophysics laboratory, including Electric ImpedanceSpectrum methodology and analysis. ! pages xiii, 16, 17, 19G. Fiandaca, E. Auken, A. V. Christiansen, and A. Gazoty. Time-domain-inducedpolarization: Full-decay forward modeling and 1D laterally constrainedinversion of Cole-Cole parameters. Geophysics, 77(3):E213, 2012. ! pages 3,32, 167G. Fiandaca, J. Ramm, A. Binley, A. V. Christiansen, and E. Auken. polarizationdata through 2-D inversion. pages 631–646, 2013. ! pages 3, 187J. Fink, E. McAlister, B. Sternberg, W. Wieduwilt, and S. Ward. InducedPolarization Applications and Case Histories. Society of ExplorationGeophysicists, 1990. ! pages 162M. Flis, G. Newman, and G. Hohmann. Induced olarization effects in timedomainelectromagnetic measurements. GEOPHYSICS, 54(4):514–523, 1989. !pages 163C. Flores and S. A. Peralta-Ortega. Induced polarization with in-loop transientelectromagnetic soundings: A case study of mineral discrimination at El Arcoporphyry copper, Mexico. Journal of Applied Geophysics, 68(3):423–436,2009. ! pages 136D. Fournier, S. Kang, M. S. McMillan, and D. W. Oldenburg. Inversion ofairborne geophysics over the DO-27/DO-18 kimberlites Part 2:Electromagnetics. Interpretation, 5(3):T313–T325, 2017. ! pages vi, 143A. Gazoty, G. Fiandaca, J. Pedersen, E. Auken, A. V. Christiansen, and J. K.Pedersen. Application of time domain induced polarization to the mapping oflithotypes in a landfill site. Hydrology and Earth System Sciences, 16(6):1793–1804, 2012. doi:10.5194/hess-16-1793-2012. ! pages 189GPG. Geophysics for practicing geoscientists, 2018. URL https://gpg.geosci.xyz/content/induced polarization/induced polarization physical properties.html. !pages xii, 10192E. Haber. Computational Methods in Geophysical Electromagnetics. Society forIndustrial and Applied Mathematics, Philadelphia, PA, 2014. ! pages 163,207, 208E. Haber and C. Schwarzbach. Parallel inversion of large-scale airbornetime-domain electromagnetic data with multiple OcTree meshes. InverseProblems, 30(5):055011, may 2014. ! pages 105L. J. Heagy, R. Cockett, S. Kang, G. K. Rosenkjaer, and D. W. Oldenburg. Aframework for simulation and inversion in electromagnetics. Computers &Geosciences, 107:1–19, 2017. ISSN 0098-3004.doi:https://doi.org/10.1016/j.cageo.2017.06.018. ! pages 67, 209G. W. Hohmann and G. A. Newman. Transient electromagnetic responses ofsurficial, polarizable patches. GEOPHYSICS, 55(8):1098–1100, 1990. !pages 163A. Ho¨rdt, T. Hanstein, M. Ho¨nig, and F. M. Neubauer. Efficient spectralIP-modelling in the time domain. Journal of Applied Geophysics, 59(2):152–161, jun 2006. ! pages 3, 45, 88, 167A. Ho¨rdt, R. Blaschek, A. Kemna, and N. Zisser. Hydraulic conductivityestimation from induced polarisation data at the field scale the Krauthausencase history. Journal of Applied Geophysics, 62(1):33–46, 2007. ! pages 2,189J. Jansen and K. Witherly. The Tli Kwi Cho kimberlite complex, NorthwestTerritories, Canada: A geophysical case study, chapter 289, pages 1147–1150.2004. ! pages 4J. C. Jansen and B. J. Doyle. The Tli Kwi Cho Kimberlite Complex, NorthwestTerritories, Canada: A Geophysical Post Mortum. Northwest MiningAssociation Practical Geophysics (III), 2000. ! pages 136V. Kaminski and A. Viezzoli. Modeling induced polarization effects in helicoptertime-domain electromagnetic data: Field case studies. GEOPHYSICS, 82(2):B49–B61, 2017. ! pages 4, 48, 92S. Kang and D. Oldenburg. Recovering IP information in airborne-time domainelectromagnetic data. ASEG Extended Abstracts, (1):1–4, jan 2015. ! pages vS. Kang and D. W. Oldenburg. TEM-IP: Extracting more IP information fromgalvanic source time domain EM data. Geophysical propecting, 2017. ! pagesvi193S. Kang, D. Oldenburg, D. Yang, and D. Marchant. On recovering inducedpolarization information from airborne time domain EM data. In SEGTechnical Program Expanded Abstracts 2014, pages 1785–1789. Society ofExploration Geophysicists, aug 2014. ! pages v, 136S. Kang, D. Fournier, and D. Oldenburg. Inversion of airborne geophysics overthe DO-27/DO-18 kimberlites Part 3: Induced polarization. Interpretation,pages T327–T340, mar 2017. ! pages vi, 4, 5, 48C. T. Kelley. Iterative Methods for Optimization. Society for Industrial andApplied Mathematics, 1999. ! pages 81A. Kemna, A. M. Binley, and L. Slater. Cross-borehole IP imaging forengineering and environmental applications. Geophysics, 69(1):97–107, 2004.! pages 3, 32A. Kemna, A. Binley, G. Cassiani, E. Niederleithinger, A. Revil, L. Slater, K. H.Williams, A. F. Orozco, F. H. Haegel, A. Ho¨rdt, S. Kruschwitz, V. Leroux,K. Titov, and E. Zimmermann. An overview of the spectral inducedpolarization method for near-surface applications. Near Surface Geophysics, 10(6):453–468, 2012. ! pages 2, 11R. Kohlrausch. Theorie des elektrischen Ru¨ckstandes in der Leidener Flasche.Annalen der Physik, 167(2):179–214, 1854. ! pages 88, 187, 206N. Kozhevnikov and E. Antonov. Fast-decaying inductively induced polarizationin frozen ground: A synthesis of results and models. Journal of AppliedGeophysics, 82:171–183, jul 2012. ! pages 136T. Kratzer and J. Macnae. Induced polarization in airborne EM. Geophysics, 77(5):E317—-E327, 2012. ! pages 5, 48, 92, 163M. Lajaunie, P. Maurya, and G. Fiandaca. Comparison of Cole-Cole and ConstantPhase Angle modeling in time-domain induced polarization. 4th IP workshop,B(4):1192–1202, 2016. ! pages 179Y. Li and D. W. Oldenburg. 3-D inversion of magnetic data. Geophysics, 61(2):394–408, 1996. ! pages 81, 172Y. Li and D. W. Oldenburg. 3-D inversion of induced polarization data.Geophysics, 65(6):1931–1945, nov 2000. ! pages 37, 82J. C. Macnae. Geophysical Prospecting with Electric Fields from an Inductive EMSource: Ph.D dissertation. PhD thesis, University of Toronto, 1988. ! pages 4194D. Marchant, E. Haber, and D. Oldenburg. Recovery of 3D IP distribution fromairborne time-domain EM. ASEG Extended Abstracts, 144:1–4, 2013. ! pages32D. Marchant, E. Haber, and D. W. Oldenburg. Three-dimensional modeling of IPeffects in time-domain electromagnetic data. Geophysics, 79(6):E303–E314,2014. ! pages v, 22, 32, 48, 99, 163, 168D. J. Marshall and T. R. Madden. INDUCED POLARIZATION, A STUDY OFITS CAUSES. GEOPHYSICS, 24(4):790–816, 1959. ! pages 11D. Oldenburg and Y. Li. Inversion of induced polarization data. Geophysics, 59(9):1327–1341, 1994. ! pages 3, 32, 33, 37, 40, 79, 82, 162, 167, 176D. W. Oldenburg and Y. Li. 5. Inversion for Applied Geophysics: A Tutorial. InNear-Surface Geophysics, chapter 5, pages 89–150. 2005. ! pages 81D. W. Oldenburg, Y. Li, and R. G. Ellis. Inversion of geophysical data over acopper gold porphyry deposit: A case history for Mt. Milligan. Geophysics, 62(5):1419–1431, 1997. ! pages 117, 162D. W. Oldenburg, E. Haber, and R. Shekhtman. Three dimensional inversion ofmultisource time domain electromagnetic data. GEOPHYSICS, 78(1):E47–E57, 2013. ! pages 163, 172W. Pelton, S. Ward, P. Hallof, W. Sill, and P. Nelson. MINERALDISCRIMINATION AND REMOVAL OF INDUCTIVE COUPLING WITHMULTIFREQUENCY IP. Geophysics, 43(3):588–609, 1978. ! pages xi, xii,xxv, 1, 3, 6, 7, 10, 11, 15, 22, 88, 152, 156, 162, 167, 187M. Power and D. Hildes. Geophysical strategies for kimberlite exploration innorthern Canada. In Proceedings of Exploration ’07: Fifth DecennialInternational Conference on Mineral Exploration, pages 1025–1031, 2007. !pages xi, 136, 138A. Revil. On charge accumulation in heterogeneous porous rocks under theinfluence of an external electric field. 78(4), 2013. ! pages xiii, 11, 13A. Revil and N. Florsch. Determination of permeability from spectral inducedpolarization in granular media. Geophysical Journal International, 181(3):1480–1498, 2010. ! pages 189A. Revil, N. Florsch, and C. Camerlynck. Spectral induced polarizationporosimetry. Geophysical Journal International, 198(2):1016–1033, 2014. !pages 156195A. Revil, G. Z. Abdel Aal, E. A. Atekwana, D. Mao, and N. Florsch. Inducedpolarization response of porous media with metallic particles Part 2:Comparison with a broad database of experimental data. Geophysics, 80(5):D539–D552, 2015. ! pages xii, 10, 11, 12, 162P. S. Routh and D. W. Oldenburg. Electromagnetic coupling in frequency-domaininduced polarization data: a method for removal. Geophysical JournalInternational, 145(1):59–76, 2001. ! pages 4, 22, 163H. Seigel. Mathematical formulation and type curves for induced polarization.Geophysics, 24(3):547–565, 1959. ! pages 3, 37, 40, 41H. Seigel. The magnetic induced polarization (MIP) method. Geophysics, 39(3):321–339, 1974. ! pages 33, 162L. Slater and D. P. Lesmes. Electrical-hydraulic relationships observed forunconsolidated sediments. Water Resources Research, 38(10):33–46, 2002. !pages 189R. S. Smith and J. Klein. A special circumstance of airborne inducedpolarizationmeasurements. Geophysics, 61(1):66–73, jan 1996. ! pages 4, 5, 48, 92, 136,187R. S. Smith and G. F. West. Field examples of negative coincidentloop transientelectromagnetic responses modeled with polarizable halfplanes.GEOPHYSICS, 54(11):1491–1498, 1989. ! pages 21R. S. Smith, P. W. Walker, B. D. Polzer, and G. F. West. The time-domainelectromagnetic response of polarizable bodies: an approximate convolutionalgorithm. Geophysical Prospecting, 36(April):772–785, 1988. ! pages 44,57, 58A. Tarasov and K. Titov. On the use of the Cole-Cole equations in spectralinduced: Polarization. i, 195(1):352–356, 2013. ! pages 16, 88, 162, 187A. N. Tikhonov and V. Y. Arsenin. Solutions of Ill-Posed Problems. W.H.Winston and Sons., 1977. ! pages 80P. Veeken, P. Legeydo, I. Pesterev, Y. Davidenko, E. Kudryavceva, and S. Ivanov.Geoelectric modelling with separation between electromagnetic and inducedpolarization field components. first break, 27(12):53–64, 2009a. ! pages 2,163196P. C. Veeken, P. J. Legeydo, Y. A. Davidenko, E. O. Kudryavceva, S. A. Ivanov,and A. Chuvaev. Benefits of the induced polarization geoelectric method tohydrocarbon exploration. Geophysics, 74:B47–B59, 2009b. ! pages 2, 163A. Viezzoli, V. Kaminski, and G. Fiandaca. Modeling induced polarization effectsin helicopter time domain electromagnetic data: Synthetic case studies.GEOPHYSICS, 82(2):E31–E50, 2017. ! pages 5J. R. Wait and T. P. Gruszka. On electromagnetic coupling ”removal” frominduced polarization surveys. Geoexploration, 24(1):21–27, 1986. ! pages162P. Weidelt. Response characteristics of coincident loop transient electromagneticsystems. GEOPHYSICS, 47(September):1325–1330, 1982. ! pages 4, 48, 92,136A. Weller, S. Nordsiek, and W. Debschutz. Estimating permeability of sandstonesamples by nuclear magnetic resonance and spectral-induced polarization.Geophysics, 75(6):E215, 2010. ! pages 189J. Wong. An electrochemical model of the induced-polarization phenomenon indisseminated sulfide ores. Geophysics, 44(7):1245, 1979. ! pages 10, 11J. C. Wynn and K. L. Zonge. EM coupling, its intrinsic value, its removal and thecultural coupling problem. Geophysics, 40(5):831–850, 1975. ! pages 4, 163Z. Xu and M. S. Zhdanov. Three-Dimensional Cole-Cole Model Inversion ofInduced Polarization Data Based on Regularized Conjugate Gradient Method.Geoscience and Remote Sensing Letters, IEEE, 12(6):1180–1184, jun 2015. !pages 32D. Yang and D. W. Oldenburg. Three-dimensional inversion of airbornetime-domain electromagnetic data with applications to a porphyry deposit.Geophysics, 77(2):B23–B34, 2012. ! pages xxii, 117, 119D. Yang, D. W. Oldenburg, and E. Haber. 3-D inversion of airborneelectromagnetic data parallelized and accelerated by local mesh and adaptivesoundings. Geophysical Journal International, 196(3):1492–1507, 2014. !pages vi, xxii, 117, 118, 124K. Yee. Numerical solution of initial boundary value problems involvingmaxwell’s equations in isotropic media. IEEE Transactions on Antennas andPropagation, 14(3):302–307, 1966. ! pages 207197Yuval and D. Oldenburg. Computation of ColeCole parameters from IP data.Geophysics, 62(2):436–448, 1997. ! pages 3, 45, 88, 167, 206198Appendix AHandling Multiple Sources forAEM SurveysThe linearization for inductive sources in Chapter 2 has been developed for a singlesource and 3D information about chargeability can be obtained if there are multiplereceivers. For ATEM data however, here is only a single receiver location for eachsource but There are multiple source locations. Our goal is to alter the problem towork with an effective pseudo-chargeability.In our linearized eq. (2.38), each source has its own sensitivity and pseudo-chargeabilty. For our airborne case the sensitivity for the k-th source is the k-throw of J and the pseudo-chargeability is h˜k. The corresponding IP datum isdIPk (t) =nCÂi=1Jk,ih˜ki (t), k = 1, . . . ,nTx, (A.1)where nTx is the number of sources, nC is the number of cells in the domain, andJk,i indicates an element of the Jacobian matrix for the k-th source and the i-th cell.I want to replace h˜ki with a single effective pseudo-chargeabilty h˜i and thereforewrite the IP datum asdIPk (t) =nCÂi=1Jk,ih˜i(t), k = 1, . . . ,nTx, (A.2)199The waveforms are different for each source and hence this representation cannotbe exact. To examine the implications of this it suffices to look at the contributionof any volumetric pixel. Each pixel contributes to all of the IP data but in differingamounts. The total contribution of the i-th pixel to the nTx data set at a single timeisqi =nTxÂk=1Jk,ih˜ki (t), i= 1, . . . ,nC. (A.3)Our goal is to find an effective chargeability that produces the same net effect onthe measured data. I search for a source-independent h˜i such thatqesti =nTxÂk=1Jk,ih˜i(t), i= 1, . . . ,nC. (A.4)Minimizing the least squares difference between eqs (A.3) and (A.4) yieldsh˜i(t) =SnTxk=1J2k,ih˜ki (t)SnTxk=1J2k,i=nTxÂk=1aki h˜ki (t), i= 1, . . . ,nC. (A.5)where the normalized weight (aki ) isaki =J2k,iSnTxk=1J2k,i, i= 1, . . . ,nC. (A.6)With the above understanding about how h˜i relates to the h˜ki from each sourceI can proceed as follows. Firstly, from eq. (2.23) I haveh˜ki (t) = h˜ I⌦ wˆki (t) (A.7)Substituting eqs (A.7) into (A.5) allows us to writeh˜i(t) = h˜ I(t)⌦wei (t), (A.8)where I define the effective time history of the electric field, wei (t) aswei (t) =nTxÂk=1aki wˆki (t), i= 1, . . . ,nC. (A.9)200The above equations shows that the pseudo-chargeabilty for any pixel recov-ered from the inversion is equal to the convolution of the impulse pseudo-chargeability,h˜ I(t), with an effective time history of the electric field we(t). Although it is some-what involved, the we(t) associated with each pixel can be evaluated by know-ing the electric fields associated with the fundamental EM problem. Ultimatelythis allows us to estimate the parameters associated with the impulse pseudo-chargeability in the same manner as outlined for the case with a single source.Our ability to evaluate the we(t) and test the validity of eq. (A.2) is treated below.For each pixel I have equation:h˜i(t) = h˜ Ii (t)⌦wei (t), (A.10)where h˜ Ii (t) is the impulse pseudo-chargeability associated with an individual pixel.The effective time history of the electric field, wei (t) is a linear combination of thefundamental electric fields due to the individual sources. I can calculate wei (t) andcarry out the convolution to evaluate the effective pseudo-chargeability. The IPdata can then be forward modelled using eq. (2.38). This allows us to validate eq.(A.2), which demonstrated the linear form of dIP data at all source locations, bycomparing results with the true IP data obtained via forward modelling. It is onlynecessary to apply this to the conductive model.The evaluation of the effective pseudo-chargeability is carried out on a cell bycell basis. For each cell I first evaluate we(t) (eq. A.9). This requires calculatingnormalized weights shown in eq. (A.6). Fig. A.1 shows these weights at a singlepixel located at (0 m,0 m,-75 m). These decay away from the center pixel becauseof the decay of the sensitivity functions. Because those are weights used to com-pute we(t), I could expect that the computed we(t) will be mostly affected by wˆkfrom a few stations close to the center. In Fig. A.2, I provide both wˆk (dashed lines)from all source locations and we(t) (eq. A.9; solid line). The we(t) is dominantlyaffected by the wˆ(t) at the center source location (solid circles)). Considering thatthe sources are 50 m apart, the decay of the sensitivity from center source locationto others is substantial (⇠ 1/r3). This results in the greatest normalized weight atthe center source location, and the observed result about we(t) is caused by this.we(t) is convolved with h˜ I(t) to compute the effective h˜(t) for that cell. When201Figure A.1: Normalized weights for the conductive case for all source loca-tions. A single pixel located at (0 m, 0 m, -75 m) is used.this is carried out for each cell then the approximate IP responses can be computedusing eq. (2.38). These can be compared with the true IP responses. Fig. A.3shows the comparisons at 0.86 ms. The images are nearly identical in shape butthe approximate IP responses are nearly a factor of two lower than the true values.This is not entirely unexpected. A similar effect was observed for IP responses fora single source shown in Fig. 2.14. At 0.86 ms, the approximate value was about70 percent of the true dIP. These results seem to be a worst case scenario. Thediscrepancy for a conductive body lessens as time increases and analyses for thecanonical and resistive bodies shows that the approximate and true IP data are invery good agreement (Section 2.3.4).202Figure A.2: Time decays of we(t) and wˆ(t) for the conductive case. A singlepixel located at (0 m, 0 m, -75 m) is used. Solid line and dashed linescorrespond to we(t) and wˆk(t) for all sources (k = 1, . . . , nTx); wˆk atthe center source located at (0 m, 0 m, 30 m) is marked as solid circles.A number of we(t) curves are overlaid due to the symmetric position ofsource locations to the conductive block.203Figure A.3: Comparison of true and approximate bIPz responses at 0.86 ms ona plan view map.204Appendix BExtracting Intrinsic IPParametersThe output of our IP inversion is a 3D distribution of the pseudo-chargeabilityat multiple time channels. As its name suggests, pseudo-chargeability is not anintrinsic IP parameter like chargeability, but it is a convolution of h˜ I(t) and wˆ(t):h˜(t) = h˜ I(t)⌦ wˆ(t), (B.1)using the definition of impulse pseudo-chargeability (eq. 2.7). We now use the h˜(t)as the data and recover intrinsic parameters such as h ,t,c in a Cole-Cole model.Assuming a Debye model (c=1), we obtainh˜ I(t) = h(1h)t e t(1h)t , (B.2)Since we have s• we can compute wˆ(t), which is the time history of the electricfield. Accordingly, we can unravel the recovered pseudo-chargeability to extractintrinsic IP parameters such as chargeability(h) and time constant (t). We usea gradient-based optimization and thus we need the sensitivity function for thepseudo-chargeability (eq. B.1) with respect to h and t . To simplify this procedure,we rewrite impulse pseudo-chargeability ash˜ I(t) = aebt , (B.3)205where a= h(1h)t and b=1(1h)t . Then we take the derivative of h˜(t) with regardto a and b:∂ h˜(t)∂a= ebt ⌦ wˆ(t), (B.4)∂ h˜(t)∂b=atebt ⌦ wˆ(t). (B.5)With these sensitivity functions, we can set up an inverse problem, and recover aand b. The chargeability and time constant can be obtained from a and b:h = ab, (B.6)t = 1(1a/b)b . (B.7)We apply this inversion separately to each cell in the recovered pseudo-chargeabilityin a manner similar to Yuval and Oldenburg [1997]. For a better alternative (rep-resentation) of time-dependent conductivity, a different parameterization such asstretched-exponential [Kohlrausch, 1854] or Cole-Cole model with variable c canbe implemented.206Appendix CDiscretizationC.1 Steady-state Maxwell’s EquationsAs shown in eq. (2.29), computation of our linearized kernel requires solvingsteady-state Maxwell’s equations. We discretize this system using a mimetic finitevolume (FV) method with weak formulation [Yee, 1966, Haber, 2014]. For thediscretization, we assume that the electric field~e is discretized by a grid function eon cell edges and the magnetic flux density~b is discretized by agrid function b oncell faces. The electrical potential f is discretized by a grid function f on the cellnodes. For a clear representation of the derivation, recall Maxwell’s equations insteady state are~j = s•~e=s•~—f , (C.1)— ·~j = — ·~js, (C.2)~j∂W · nˆ= 0, (C.3)where ∂W indicates a boundary surface of the system and nˆ is the normal vector ofthe boundary surface. The weak form of those equations can be written as(~j,~w)+(s•~—f ,~w) = 0, (C.4)(~j,~—y) = (~js,~—y). (C.5)207The inner products (~j,~w), (s•~—f ,~w), (~j,~—y) and (~js,~—y) are edge-based prod-ucts. Here we define the inner product as(~a,~b) =ZW~a ·~bdv, (C.6)where W is the volume of the system. By discretizing the ~— operator and the innerproduct in space, we obtainMej+Mes•Gf = 0, (C.7)GTMej=GTMejs, (C.8)whereMe performs volume averaging, andMes• is the mass matrix of conductivity(s•), which discretizes the edge based inner product. For further details on theformation of this matrix see Haber [2014].By substituting eq. (C.7) into (C.8), we haveAs•f = rhsDC, (C.9)where As• =GTMes•G and rhsDC =GTMejs. We use the SIMPEG’s tensor meshand solver classes to form and solve above the linear system [Cockett et al., 2015].C.2 Linearized Kernel for IP ResponsesTo obtain a linear form of eq. (2.38), we first discretize the Biot-Savart law shownin eqs (2.36) and (2.37). In our discretization ~jIP and h˜ are defined at the cellcenters, and those for each time channel are constant in a cell volume, whereas~e re f is defined on the cell edges. We define the number of cells and edges in 3Dspace as nC and nE, respectively. The discretized IP current density, jIPcc 2 R3nC1 , isdefined at the cell center. Since ~jIP has three components, we first discretize theintegration operator including cross product (Rv⇥rˆr2 dv) asGBiot =264eT 0 00 eT 00 0 eT375264 0 Sz SySz 0 SxSy Sx 0375 , (C.10)208whereSl = diag(v rl 1r2 ), l = x, y, zand the electric field at peak time, eFmax 2RnE1 is a column vector, diag(·) is the diag-onal matrix and is the Hadamard product. Computing the eFmax requries a forwardsimulation time-domain electromagnetic problem, and SIMPEG-EM package isused [Heagy et al., 2017]. Then we discretize ~jIP shown in eq. (2.28) asjIPcc (t) = Sdiag(eFmax)AecTdiag(v)diag(s•)h˜(t), (C.11)where Aec is a discrete averaging matrix from edge to cell center andS= AeccvMe1[Mes•GA1s•GT  I]diag(eFmax)AecTdiag(v)diag(s•). (C.12)HereAeccv is a discrete averaging matrix from edge to cell center with considerationof three component vector: 2R3nCnE . Thus, we can have a linear equation for a singletime channel asbIP =GBiotSh˜ ,Finally, by lettingJ=GBiotS, (C.13)we havebIP = Jh˜ , (C.14)where J is the Jacobian matrix of the linear equation, and since J is static, we alsoobtain ∂bIP∂ t= J(∂ h˜∂ t ). (C.15)209

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0364164/manifest

Comment

Related Items