Cooperative and parametric strategies for3D airborne electromagnetic inversionbyMichael S. G. McMillanB.Sc., The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Geophysics)The University of British Columbia(Vancouver)April 2017c©Michael S. G. McMillan, 2017AbstractAirborne electromagnetic (AEM) data is commonly collected for detecting buriednatural resources, and this technique is sensitive to subsurface electrical resistivitydistributions. The subsequent process of 3D AEM inversion constructs a physicalproperty model from this data in order to better understand the size and shape ofpotential hidden resources. This thesis is designed to develop practical strategiesto improve 3D AEM inversion accuracy in geologic settings where AEM data setsproduce inconsistent or unsatisfactory inversion results.In this research, two overarching problematic scenarios are examined. First,in regions where an AEM survey overlaps with other electromagnetic data sets, anovel cooperative approach is introduced. This method is first tested on syntheticdata where instead of producing an inversion model from each data set, the cooper-ative algorithm finds one consistent 3D resistivity model with improved resolution.The approach is then applied to field data over a high-sulfidation epithermal golddeposit where similar improvements are observed.The second scenario relates to improving 3D AEM inversions over thin conduc-tive anomalies, a common geophysical target for copper and gold deposits. A newparametric inversion is developed using skewed Gaussian ellipsoids to representtarget bodies. The approach is general but applied to frequency and time-domainAEM data with one or multiple anomalies. Combined with a voxel algorithm, theparametric inversion forms a hybrid approach capable of resolving thin conduc-tive targets with smooth surrounding features. This hybrid technique is tested onsynthetic data over conductive targets in a resistive background, and consistentlyproduces models that are easier to interpret compared to voxel inversions alone.Field examples from a volcanogenic massive sulfide and an orogenic gold depositiihighlight the practical nature of this method to image conductive mineralizationwith increased precision.The thesis concludes by analyzing a setting where multiple spatially overlap-ping AEM data sets exist over narrow conductive anomalies. Here, parametric,cooperative and voxel inversions are applied together to generate a consistent 3Dresistivity model with thin targets and smooth background features. This sectionincludes a discussion about potential pitfalls of such an approach when incompati-ble field measurements are encountered.iiiPrefaceThe 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. Por-tions of this research can be located in two peer-reviewed journal articles as wellas seven expanded conference proceedings, and I have presented this work at nineconferences worldwide.Chapter 2 contains text and figures from the published paper McMillan andOldenburg (2014) in Geophysics. I came up with the idea, carried out the numericalexperiments and wrote the manuscript. Dr. Douglas Oldenburg helped to edit themanuscript and Dr. Eldad Haber assisted with some of the numerical content. Thework first appeared in the SEG expanded abstract McMillan and Oldenburg (2012).Chapter 3 is a revised version of the published paper McMillan et al. (2015c)in Geophysics. The idea developed from discussions with Dr. Eldad Haber, Dr.Elliot Holtham and myself. I performed the numerical experiments with help fromDr. Christoph Schwarzbach, and I wrote the journal paper with editing help fromDr. Douglas Oldenburg. The work was first shown in the SEG expanded abstractMcMillan et al. (2014).Chapter 4 is currently in preparation for submission. The idea was my own, Idid the numerical experiments with help from Dr. Christoph Schwarzbach, and Iwrote the resulting chapter with edits from Dr. Douglas Oldenburg. The work firstappeared in the EAGE expanded abstract McMillan et al. (2015b).Chapter 5 is part of ongoing research that will lead to a future publication.I thought of the idea, and I completed the numerical trials with help from Dr.Christoph Schwarzbach and Dr. Eldad Haber. I wrote the chapter with edits fromivDr. Douglas Oldenburg, and the work was first shown in the EAGE expandedabstract McMillan et al. (2016).Appendix A and Appendix B contain material relevant to Chapter 4 and Chap-ter 5 and represents more of my original research.This thesis manuscript is entirely my own work, and I prepared it and wrote itwith editing help from Dr. Douglas Oldenburg and Dr. Eldad Haber.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Electromagnetic Survey Techniques . . . . . . . . . . . . . . . . 21.3 Maxwell’s Equations . . . . . . . . . . . . . . . . . . . . . . . . 41.4 Electromagnetic Modeling and Inversion . . . . . . . . . . . . . . 71.5 Airborne Electromagnetic Interpretation Techniques . . . . . . . . 91.5.1 Semi-quantitative methods . . . . . . . . . . . . . . . . . 91.5.2 Voxel inversion methods . . . . . . . . . . . . . . . . . . 91.5.3 Cooperative and joint inversion . . . . . . . . . . . . . . 101.5.4 Constrained inversion . . . . . . . . . . . . . . . . . . . 111.5.5 Parametric methods . . . . . . . . . . . . . . . . . . . . . 111.6 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12vi2 Cooperative Inversion of Multiple Electromagnetic Data Sets . . . . 142.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2 Antonio Geologic Background . . . . . . . . . . . . . . . . . . . 162.3 Antonio Geophysical Data Sets . . . . . . . . . . . . . . . . . . . 192.3.1 Time-domain AEM . . . . . . . . . . . . . . . . . . . . . 192.3.2 CSAMT . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3.3 DC resistivity . . . . . . . . . . . . . . . . . . . . . . . . 202.4 Inversion Methods . . . . . . . . . . . . . . . . . . . . . . . . . 202.4.1 Inversion preparation . . . . . . . . . . . . . . . . . . . . 202.4.2 Inversion methodology . . . . . . . . . . . . . . . . . . . 232.4.3 Individual Antonio inversions . . . . . . . . . . . . . . . 232.4.4 Cooperative inversion work flow . . . . . . . . . . . . . . 262.5 Synthetic Inversion Results . . . . . . . . . . . . . . . . . . . . . 282.5.1 Individual synthetic inversions . . . . . . . . . . . . . . . 282.5.2 Synthetic cooperative inversion . . . . . . . . . . . . . . 302.5.3 Synthetic constrained cooperative inversion . . . . . . . . 342.6 Antonio Inversion Results . . . . . . . . . . . . . . . . . . . . . . 362.6.1 Antonio cooperative inversion . . . . . . . . . . . . . . . 362.6.2 Antonio constrained cooperative inversion . . . . . . . . . 382.6.3 Geologic interpretation . . . . . . . . . . . . . . . . . . . 402.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423 Single Anomaly Parametric Inversion of AEM data . . . . . . . . . 433.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.2 Electromagnetic Background . . . . . . . . . . . . . . . . . . . . 463.3 Inversion Methodology . . . . . . . . . . . . . . . . . . . . . . . 463.4 Parametric Sensitivity and Initial Parameter Selection . . . . . . . 513.5 Synthetic Results . . . . . . . . . . . . . . . . . . . . . . . . . . 563.5.1 Thin plate responses . . . . . . . . . . . . . . . . . . . . 563.5.2 Synthetic dipping plate inversion results . . . . . . . . . . 573.6 Caber Case Study Results . . . . . . . . . . . . . . . . . . . . . . 633.6.1 Caber geology . . . . . . . . . . . . . . . . . . . . . . . 633.6.2 Caber time-domain AEM . . . . . . . . . . . . . . . . . . 64vii3.6.3 Caber parametric inversions . . . . . . . . . . . . . . . . 653.6.4 Caber hybrid parametric inversion . . . . . . . . . . . . . 663.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714 Multiple Anomaly Parametric Inversion of AEM data . . . . . . . . 734.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.2 Inversion methodology . . . . . . . . . . . . . . . . . . . . . . . 744.3 Synthetic Results . . . . . . . . . . . . . . . . . . . . . . . . . . 774.3.1 Synthetic model . . . . . . . . . . . . . . . . . . . . . . 784.3.2 Synthetic data . . . . . . . . . . . . . . . . . . . . . . . . 804.3.3 Synthetic voxel inversions . . . . . . . . . . . . . . . . . 814.3.4 Synthetic parametric inversions . . . . . . . . . . . . . . 814.4 West Plains Case Study Results . . . . . . . . . . . . . . . . . . . 874.4.1 West Plains geology . . . . . . . . . . . . . . . . . . . . 874.4.2 West Plains AEM data . . . . . . . . . . . . . . . . . . . 874.4.3 West Plains voxel inversions . . . . . . . . . . . . . . . . 894.4.4 West Plains parametric inversions . . . . . . . . . . . . . 894.4.5 West Plains hybrid parametric inversions . . . . . . . . . 934.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 995 Cooperative Multiple Anomaly Parametric Inversion of AEM Data 1005.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.2 Inversion Methodology . . . . . . . . . . . . . . . . . . . . . . . 1015.3 Synthetic Results . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.3.1 Cooperative parametric inversion . . . . . . . . . . . . . 1025.4 West Plains Case Study Results . . . . . . . . . . . . . . . . . . . 1055.4.1 Cooperative parametric inversions . . . . . . . . . . . . . 1055.4.2 Hybrid cooperative parametric inversion . . . . . . . . . . 1095.5 Cooperative Inversion Discussion . . . . . . . . . . . . . . . . . 1125.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1156 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166.1 Cooperative Inversion . . . . . . . . . . . . . . . . . . . . . . . . 1166.2 Parametric Inversion - Single Anomaly . . . . . . . . . . . . . . . 117viii6.3 Parametric Inversion - Multiple Anomalies . . . . . . . . . . . . . 1186.4 Cooperative Parametric Inversion . . . . . . . . . . . . . . . . . . 1196.5 Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . 120Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121A Multiple Sphere Parametric Inversion . . . . . . . . . . . . . . . . . 129A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129A.2 Sphere Parameterization . . . . . . . . . . . . . . . . . . . . . . 129A.3 Synthetic Example . . . . . . . . . . . . . . . . . . . . . . . . . 131A.4 Synthetic Result . . . . . . . . . . . . . . . . . . . . . . . . . . . 131B Airborne Electromagnetic Data Sampling . . . . . . . . . . . . . . . 133B.1 Total Horizontal Gradient Data Sampling . . . . . . . . . . . . . 133ixList of TablesTable 2.1 Error assignments and final data misfits for individual field andsynthetic inversions. . . . . . . . . . . . . . . . . . . . . . . . 25Table 2.2 Final data misfits for cooperative synthetic and field inversions. 33Table 2.3 Quantitative assessment of synthetic inversions using a residual(R) value. A lower residual refers to a more accurate recovery. 35Table 4.1 Synthetic data inversion summary. . . . . . . . . . . . . . . . 83Table 4.2 Field data inversion summary. . . . . . . . . . . . . . . . . . . 92Table 5.1 Synthetic data parametric inversion summary. . . . . . . . . . 105Table B.1 Synthetic inversions summary: down-sampling methods. . . . 136xList of FiguresFigure 1.1 Simplified helicopter airborne electromagnetic setup for a genericfrequency and time-domain system. ρ0 is the background re-sistivity and ρ1 represents the target resistivity. . . . . . . . . 5Figure 1.2 Generic frequency and time-domain waveforms used for anairborne electromagnetic survey. a) Frequency domain. b)Time domain. . . . . . . . . . . . . . . . . . . . . . . . . . . 6Figure 2.1 Map of Peru with location of Antonio deposit marked as ablack star. Modified from McMillan and Oldenburg (2014). . 16Figure 2.2 Yanacocha regional alteration map with location of major golddeposits. (modified from Hoschke (2011)). . . . . . . . . . . 17Figure 2.3 Antonio geology with lithologies, structures and near-surfacesilica alteration outlined. . . . . . . . . . . . . . . . . . . . . 18Figure 2.4 NEWTEM-I discretized waveform in red with measured timechannels marked as black crosses. . . . . . . . . . . . . . . . 19Figure 2.5 Antonio geophysics locations. Receiver sites for AEM (orangecircles), CSAMT (red stars), and DC Resistivity (green trian-gles) overlaid on topography with near-surface silica alterationoutlined in dashed red and geologic faults in thin blue. . . . . 21Figure 2.6 Observed and predicted field data from individual inversions.a) AEM plan view of ∂bz∂ t at 139 µs b) CSAMT plan view ofy-component electric field amplitudes at 16 Hz c) DC apparentresistivity pseudo-section from line 9230630. . . . . . . . . . 22xiFigure 2.7 Resistivity 3D inversions from field data at a 3870 m elevationslice with receiver locations shown. Near-surface silica alter-ation outline in dashed red and geologic faults in thin solidblue. a) AEM b) CSAMT c) DC Resistivity d) Cooperative. . 24Figure 2.8 Cooperative inversion work flow. σ (i)j = conductivity modelfrom the ith cooperative iteration for data set j, β = trade-offparameter, φd = data misfit, φ ∗d = target data misfit. . . . . . . 29Figure 2.9 Synthetic model resistivity at a 3870 m elevation slice. . . . . 31Figure 2.10 Resistivity 3D inversions from synthetic data at a 3870 m el-evation slice with receiver locations shown. True resistor out-line in black. a) AEM b) CSAMT c) DC Resistivity d) Coop-erative. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 2.11 Data misfit convergence curves with target misfit shown indashed black with final models circled in black. a) Syntheticdata without bounds b) Field data without bounds c) Syntheticdata with bounds d) Field data with bounds. . . . . . . . . . . 34Figure 2.12 Cooperative 3D inversions with bounds at a 3870 m eleva-tion slice a) Synthetic data with the true resistor outline in thinblack b) Field data with the outline of known near-surface sil-ica alteration in dashed red, interpreted silica alteration outlineat 3870 m elevation in solid black, geologic faults in thin blue,and conductive anomalies of interest numbered in yellow. . . . 36Figure 2.13 Resistivity 3D inversions from synthetic data at five elevationslices: 3895 m, 3870 m, 3845 m, 3820 m, 3795 m. a) Truemodel b) AEM c) CSAMT d) DC Resistivity e) Cooperative f)Cooperative with drilling bounds. . . . . . . . . . . . . . . . 37Figure 2.14 Linear regression curve, shown as a black line, for the resistiv-ity versus total sulfur relationship from 30 borehole samples. . 39Figure 2.15 Drilling information for model cells intersected by 78 bore-holes at Antonio. a) Alteration b) Total sulfur content c) Re-sistivity derived from a regression relationship between totalsulfur/alteration and resistivity. . . . . . . . . . . . . . . . . . 40xiiFigure 2.16 Resistivity 3D inversions from field data at five elevation slices:3895 m, 3870 m, 3845 m, 3820 m, 3795 m. a) Alteration fromboreholes b) AEM c) CSAMT d) DC Resistivity e) Coopera-tive. An example of a conductive anomaly potentially linkedto propylitic alteration is highlighted with a yellow star f) Co-operative with drilling bounds. . . . . . . . . . . . . . . . . . 41Figure 3.1 a) Analytic step-function with a = 10. b) Gaussian ellipsoidexample with a = 10, c = 1. c) Gaussian ellipsoid examplewith a = 10, c = 5. d) Analytic step-function with a = 0.5. e)Gaussian ellipsoid example with a = 0.5, c = 1. f) Gaussianellipsoid example with a = 0.5, c = 5. . . . . . . . . . . . . . 55Figure 3.2 Triangular discretized waveform in red with measured timechannels marked as black crosses. . . . . . . . . . . . . . . . 56Figure 3.3 Noise-free ∂bz∂ t responses from a coincident loop AEM systemflying 37.5 m over a thin 3 Ωm conductive plate in a 3000 Ωmbackground. a) Vertical plate buried 20 m below the surface.b) 45 degree dipping plate buried 20 m below the surface. c)Flat plate buried 20 m below the surface. . . . . . . . . . . . 57Figure 3.4 Plan view depth slices and cross-sections through the true model,recovered voxel-based inversion and recovered parametric mod-els. a) True model at z = -250 m. b) True model along y = 0 m.c) True model along x = -50 m. d) Voxel model at z = -250 m.e) Voxel model along y = 0 m. f) Voxel model along x = -50 m.g) Parametric model at z = -250 m for fixed ρ . h) Parametricmodel along y = 0 m for fixed ρ . i) Parametric model alongx = -50 m for fixed ρ . j) Parametric model at z = -250 m forvariable ρ . k) Parametric model along y = 0 m for variable ρ .l) Parametric model along x = -50 m for variable ρ . . . . . . . 59xiiiFigure 3.5 Synthetic dipping plate observed and predicted ∂bz∂ t data. a)Observed data at 1110 µs for fixed ρ , with a selected soundingmarked with a cross. b) Predicted data at 1110 µs for fixedρ , with a selected sounding marked with a cross. c) Fixed ρobserved and predicted data at selected sounding location. d)Fixed and variable resistivity data misfit progression. . . . . . 62Figure 3.6 a) Caber deposit location and geology, modified from Prikhodkoet al. (2010) and Adair (2011). b) Simplified deposit cross-section with drilling traces. . . . . . . . . . . . . . . . . . . . 64Figure 3.7 VTEM-35 (2012) discretized waveform in red with measuredtime channels marked as black crosses. . . . . . . . . . . . . 65Figure 3.8 Plan view depth slices and cross-sections through the initialguess and recovered Caber parametric and hybrid models. a)Initial guess at z = 142.5 m. b) Initial guess along y = 5513510m. c) Initial guess along x = 710146 m. d) Parametric modelat z = 42.5 m. e) Parametric model along y = 5513510 m. f)Parametric model along x = 710146 m. g) Hybrid model at z= 42.5 m. h) Hybrid model along y = 5513510 m. i) Hybridmodel along x = 710146 m. . . . . . . . . . . . . . . . . . . . 67Figure 3.9 Caber observed and predicted ∂bz∂ t data. a) Observed data at1010 µs with a selected sounding marked with a cross. b) Pre-dicted data at 1010 µs with a selected sounding marked witha cross. c) Observed and predicted data at selected soundinglocation. d) Data misfit progression. . . . . . . . . . . . . . . 68Figure 3.10 Caber observed and predicted data at 505 µs from the para-metric hybrid inversion. a) Observed data. b) Predicted data. . 70Figure 3.11 Caber deposit outline with inversion results. 0.4Ωm iso-surface(red) from parametric stage, and 0.4 Ωm iso-surface (hatchedgray) from hybrid inversion overlaid on Caber massive sulfidedeposit model from drilling (light gray outline) and Maxwellplate anomalies (multiple yellow sheets) from five central linesof AEM data (Prikhodko et al., 2012). a) Looking north-east.b) Looking north-west. . . . . . . . . . . . . . . . . . . . . . 71xivFigure 4.1 Cross section through a two-anomaly parametric model withfour unique resistivity values (ρ0 - ρ3). H1 and H2 representanomaly regions. . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 4.2 The West Plains area resides within the Archean CommitteeBay greenstone belt in Eastern Nunavut, Canada . . . . . . . 77Figure 4.3 Frequency and time-domain AEM synthetic voxel inversionsin plan view at z = 250 m (left panel) and cross-section acrosssolid white lines at y = 7333080 m and y = 7332780 m (rightpanel). a) True synthetic model. b) Frequency-domain voxelinversion. c) Time-domain voxel inversion. . . . . . . . . . . 79Figure 4.4 VTEM (2003) discretized waveform in red with measured timechannels marked as black crosses. . . . . . . . . . . . . . . . 81Figure 4.5 Synthetic parametric inversions in plan view at z = 250 m (leftpanel) and cross-section across solid white lines at y = 7333080m and y = 7332780 m (right panel). Dashed white lines repre-sent starting guess locations for parametric anomalies. a) Truesynthetic model. b) Frequency-domain parametric inversion.c) Time-domain parametric inversion. . . . . . . . . . . . . . 82Figure 4.6 Observed and predicted synthetic data from parametric inver-sions with locations shown as black dots. a) Observed data -real Hz at 450 Hz. b) Predicted parametric data - real Hz at 450Hz. c) Observed data - imaginary Hz at 450 Hz. d) Predictedparametric data - imaginary Hz at 450 Hz. e) Observed data -∂bz∂ t at 150 µs. f) Predicted parametric data -∂bz∂ t at 150 µs. . . 85Figure 4.7 Synthetic parametric data misfit progression. . . . . . . . . . 86Figure 4.8 Frequency and time-domain parametric inversion 5 Ωm iso-surfaces. a) True synthetic model. b) Starting guesses. c)Frequency-domain parametric inversion. d) Time-domain para-metric inversion. . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 4.9 Simplified West Plains geology map with major lithology unitsdefined and location of overlapping frequency and time-domainAEM surveys outlined in black. Conductive komatiite units ofinterest are numbered in yellow. . . . . . . . . . . . . . . . . 88xvFigure 4.10 Frequency and time-domain field data. a) Imaginary Hz at 385Hz. b) Imaginary Hz at 115000 Hz. c) ∂bz∂ t at 150 µs. d)∂bz∂ t at3180 µs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Figure 4.11 Frequency and time-domain West Plains voxel inversions inplan view at z = 190 m (left) and cross-section across solidwhite lines at y = 7333080 m and y = 7332780 m (right). a)Frequency domain. b) Time domain. . . . . . . . . . . . . . . 91Figure 4.12 West Plains data misfit progression. a) Voxel inversions. b)Parametric and hybrid parametric inversions. . . . . . . . . . 92Figure 4.13 West Plains parametric inversions in plan view at an elevationof z = 190 m (left) and cross-section across solid white lines aty = 7333080 m and y = 7332780 m (right). Dashed white linesrepresent starting guess locations for parametric anomalies. a)Frequency domain. b) Time domain. . . . . . . . . . . . . . . 93Figure 4.14 Observed and predicted field data from parametric and hybridinversions with locations shown as black dots. a) Observeddata - real Hz at 385 Hz. b) Predicted parametric data - real Hzat 385 Hz. c) Predicted hybrid parametric data - real Hz at 385Hz. d) Observed data - imaginary Hz at 385 Hz. e) Predictedparametric data - imaginary Hz at 385 Hz. f) Predicted hybridparametric data - imaginary Hz at 385 Hz. g) Observed data -∂bz∂ t at 150 µs. h) Predicted parametric data -∂bz∂ t at 150 µs. i)Predicted hybrid parametric data - ∂bz∂ t at 150 µs. . . . . . . . 94Figure 4.15 West Plains hybrid parametric inversions in plan view at an el-evation of z = 190 m (left) and cross-section across solid whitelines at y = 7333080 m and y = 7332780 m (right). a) Fre-quency domain. b) Time domain. . . . . . . . . . . . . . . . 95Figure 4.16 West Plains hybrid parametric inversions in plan view at anelevation of z = 310 m (surface). a) Frequency domain. b)Time domain. Note the colorbars are different to highlightnear-surface features. . . . . . . . . . . . . . . . . . . . . . . 96xviFigure 4.17 30 Ωm iso-surfaces from voxel and hybrid inversions fromWest Plains field data. 1 g/t gold shapes from drilling in theWest conductor are shown in yellow and correspond with theconductive komatiite unit. a) Frequency-domain voxel. b)Time-domain voxel. c) Frequency-domain hybrid parametric.d) Time-domain hybrid parametric. . . . . . . . . . . . . . . 98Figure 5.1 Synthetic cooperative parametric inversion data misfit progres-sion with the final model circled in black. The first three seg-ments are numbered for explanation purposes. . . . . . . . . . 103Figure 5.2 Synthetic parametric inversions in plan view at z = 250 m (leftpanel) and cross-section across solid white lines at y = 7333080m and y = 7332780 m (right panel). Dashed white lines rep-resent starting guess locations for parametric anomalies. a)True synthetic model. b) Frequency-domain parametric inver-sion. c) Time-domain parametric inversion. d) Cooperativeparametric inversion. . . . . . . . . . . . . . . . . . . . . . . 104Figure 5.3 Comparison of a practical waveform with a small decayingcurrent during early off-time ∂bz∂ t measurements compared toan ideal instantaneous shut-off waveform. a) Practical wave-form (red) and ideal waveform (black) in early off-times. b)∂bz∂ t responses at early times for the practical and ideal wave-forms modeled 30 m above a 1000 Ωm half-space. . . . . . . 107Figure 5.4 West Plains cooperative parametric inversion data misfit pro-gressions with varying numbers of Gauss-Newton iterationsper cooperative iteration. Final models circled in black. a)5 Gauss-Newton iterations. b) 10 Gauss-Newton iterations. c)25 Gauss-Newton iterations. d) Hybrid cooperative inversionwith 25 Gauss-Newton iterations, with each segment numberedfor explanation purposes. . . . . . . . . . . . . . . . . . . . . 108xviiFigure 5.5 Cooperative and hybrid cooperative parametric inversions atWest Plains. Dashed white lines represent starting guess loca-tions for parametric anomalies. Plan view slices at z = 190 m(left) and cross-section across solid white lines at y = 7333080m and y = 7332780 m (right). a) Cooperative parametric model.b) Hybrid cooperative parametric model. c) Hybrid coopera-tive parametric model with a plan view slice at surface. . . . . 111Figure 5.6 Final hybrid cooperative parametric inversion at West Plainswith 1 g/t gold shapes in yellow. a) 30 Ωm iso-surface. b)Closeup view of Western conductor with 30Ωm iso-surface. c)Closeup view of Western conductor with 280 m elevation slice.Orange circle highlights area of condensed mineralization thatcorresponds with a widening of the recovered conductor. d)Cross section at y = 7333230 m. . . . . . . . . . . . . . . . . 113Figure 5.7 Final hybrid cooperative parametric inversion at West Plains.Elevation slice at z = 190 m overlaid on simplified geology map. 114Figure A.1 Multiple sphere parametric inversion, frequency-domain AEMdata. Plan view at z = 250 m (left) and cross-section acrosssolid white lines at y = 7332880 m and y = 7332640 m (right).a) True model with cluster centers marked with black x’s. b)Initial guess. c) Parametric sphere inversion. . . . . . . . . . . 132Figure B.1 Time-domain AEM data sampling techniques. Data locationsshown as black dots. a) Full synthetic ∂bz∂ t data set at 150 µswith 1172 total source locations. b) Evenly down-sampledsynthetic ∂bz∂ t data set at 150 µs with 533 total source loca-tions. c) THG of ∂bz∂ t data at 150 µs with 522 selected sourcelocations based from the THG. . . . . . . . . . . . . . . . . . 134Figure B.2 a) Full sampling. b) Even sampling. c) THG sampling. . . . . 136xviiiAcknowledgmentsThis epic adventure known as a Ph.D. degree would not have been possible with-out the help of countless people along the way. My teachers and instructors haveprovided an amazing learning pathway for which I am forever grateful, from el-ementary school in Smithers, B.C. all the way to professors in graduate school atUBC. To Doug Oldenburg, who introduced me to applied geophysics back in 2007,and who continues to support me, I owe more than words can describe. To EldadHaber, who has given me the most unique type of nebulous guidance along withplenty of criticism throughout my Ph.D., I thank for making this process as enter-taining as it has been rewarding. Roman Shekhtman also needs special thanks forbeing the glue that holds the whole UBC-GIF group together.To my fellow graduate students at GIF who have gone through this adventurewith me, thank you for your cooperation, support, encouragement and friendship.We have shared laughter, pain, happiness, suffering, joy, despair, love, heartbreakand every other possible human emotion over these six years. I look forward tofollowing each of your successes in the future.Thank you to NSERC for sponsoring much of this research, and to NewmontMining and Computational Geosciences for allowing me to keep my fingers in theworld of applied geophysics while pursuing a doctorate degree.To my close friends, from Smithers, from UBC, from my travels near and far,I appreciate all the support that you guys have provided. In all my silly endeavors,I can always count on my huge network of amazing friends for help and advice.Finally to my parents and to my brother, thank you for always believing inmy abilities, it means the world to me and I could not have completed this degreewithout each of you.xixChapter 1Introduction1.1 Research MotivationThe subsurface of the earth, in particular the top kilometer of the crust, containsmost of the available natural resources that society will need over the course of thenext generation. Whether it be in the search for metal, water, oil or other com-modities of interest, accurately mapping these resources is vitally important forconservation, extraction and monitoring purposes. Imaging this upper crustal re-gion of the earth with geophysical techniques is an appealing cost effective solutioncompared to the likes of drilling. One such geophysical method, broadly known aselectromagnetics, is sensitive to electrical resistivity, which is a physical propertythat measures the degree to which a material opposes the flow of electric current.Variations in resistivity within the earth can indicate a change in rock type, alter-ation or in-situ fluid condition. These variations can in turn help locate and evencharacterize buried resources.The usefulness of electromagnetic and specifically airborne electromagnetic(AEM) data for mineral exploration and commodity detection has been known formany decades. Companies have steadily acquired data in various configurationsover the years, and as a result, an ongoing question is how to interpret the ensuingdata sets. Only within the past decade has the AEM problem become tractablefor 3D inversion, a process designed to reconstruct a feasible 3D earth model thatbest matches the collected electromagnetic observations. This is due to software1algorithms taking advantage of the increased computing power and novel massiveparallelization options not previously available. Consequently, decades of historicdata sets over a wide range of geologic targets are now being inverted with modern3D techniques.The process of inverting historic as well as modern data sets has opened upmany new questions with regards to best practices for 3D AEM inversion. Onescenario that commonly occurs due to the abundance of historic data sets, is whatto do when two surveys have been collected over the same area but at differenttimes. The question that arises is whether the information from each survey can becombined in a useful manner. The other common problem relates to the resolutionof 3D inversion models. Often exploration companies plan drill holes based on 3Dinversion models, and for certain targets such as narrow thin anomalies, conven-tional inversion results have trouble accurately representing this compact geometry.In this thesis, I introduce novel algorithms and develop practical strategies to ad-dress these questions with regards to improving the accuracy and capabilities of3D frequency and time-domain AEM inversion.This first introductory chapter will present a review on general electromagnet-ics and relevant equations, the airborne electromagnetic method, fundamentals ofvoxel, cooperative and parametric 3D inversion and will conclude with an outlineof the thesis.1.2 Electromagnetic Survey TechniquesElectromagnetic surveys detect contrasts in electrical resistivity, which can helpdistinguish rock types and alteration zones within the earth (Keller, 1988). But,these data can be collected in a multitude of ways, and selecting an appropriatesurvey design is critical to ensure that the target of interest can be imaged. Ageneral electromagnetic layout consists of a transmitter that generates a signal thatpenetrates the ground and a set of one or more receivers that measure the inducedresponse that emanates from the earth. When designing a survey, there are manychoices to consider in terms of the type of transmitter and receiver to use and how toplace these instruments during data acquisition. Some major distinctions betweensurvey configurations are listed below.21. Type of transmitter - inductive, galvanic or natural source.2. Type of receiver measurements - electric field (e), magnetic flux (b) or timerate of change of magnetic flux ( ∂b∂ t ).3. Spatial component of measurements - x, y, or z.4. Transmitter and receiver locations - airborne, ground or borehole.5. Domain of measurements - frequency or time domain.This list is not meant to be exhaustive, but it is clear that many choices are avail-able when designing an electromagnetic survey, and there are pros and cons to eachchoice depending on the goal of the particular project. In this thesis I will focusprimarily on AEM configurations (Macnae et al., 1991) both in the frequency andtime domain. AEM surveys have been collected since the late 1950’s (Palacky andWest, 1973), and deploy inductive sources, which excite the earth using a time-varying magnetic field. Figure 1.1 shows a cartoon of a generic frequency andtime-domain helicopter AEM configuration. A frequency-domain setup generallyconsists of multiple small transmitter loops that can be oriented in any direction.Figure 1.1 depicts a simplified version with one transmitter loop in a horizontalcoplanar orientation, also known as a vertical dipole. A time-domain system of-ten possesses a much larger transmitter loop, and due to the size of the loop itis usually oriented in a horizontal fashion as shown in Figure 1.1. The heaviertime-domain loop generally maintains a higher terrain clearance compared to thefrequency-domain counterpart, but the exact terrain clearance will vary depend-ing on the survey specifications and the amount of topography present in the area.Both frequency and time-domain AEM receivers are generally small loops locatedin the air that measure components of b or ∂b∂ t from secondary fields induced in theearth by the transmitted source. These secondary fields are dependent on the trueheterogeneous resistivity distribution of the earth, but for demonstration purposescan be thought of as due to a background resistivity ρ0 and an anomalous targetregion ρ1.In the frequency domain, the transmitted waveform oscillates continuously asdisplayed in Figure 1.2a, and the recorded data from the receiver is decomposed3into real (in-phase with the transmitter source) and imaginary (out-of-phase withthe transmitter source) components. Conversely in the time domain, the transmitterwaveform is usually turned on for a period of time in a pulse-like manner with alarge amount of current and then shut off as shown in Figure 1.2b. During this‘off-time’, the receiver records the secondary fields from the earth in the absenceof the transmitted primary signal. Note that many variations of waveforms areused in time-domain AEM systems, from square waves (Eaton et al., 2013) totriangular waves (Boyko et al., 2003) to half-sine waves (Mule` et al., 2012), withmany combinations in between. Some newer systems combine two transmitterpulses, one high and one low amplitude pulse, to suit both deep and near-surfacetargets (Chen et al., 2015).Other data sets mentioned throughout this work include controlled source au-dio frequency magnetotellurics (CSAMT) (Zonge and Hughes, 1991) and the directcurrent (DC) resistivity method (Kunetz, 1966). A CSAMT survey uses receiverson the ground to measure electric fields and magnetic fluxes in response to a gal-vanic transmitter many kilometers away that injects current directly into the earthat frequencies between 0.1 Hz and 10 kHz.In contrast, a DC configuration is often considered an electrical technique, butit can be thought of as a special case of the electromagnetic method that operates ata near-zero frequency. The DC method injects current with a galvanic transmitter,and often operates in the time domain using a half-duty cycle square wave elec-tric current with an eight second period. Frequency-domain DC surveys also exist,and they inject an oscillating electric current into the ground instead of a squarewave. In either the time or frequency domain, DC data consist of electric poten-tial differences measured between grounded receiver electrodes known as dipoles.CSAMT and DC data are sensitive to electrical conductivity, and as such they areboth helpful to reconstruct conductivity distributions through Maxwell’s equations.1.3 Maxwell’s EquationsElectromagnetic induction is governed by Maxwell’s equations, which are writtenin quasi-static form in the time domain as4TransmitterReceiverHorizontal coplanar conguration.ReceiverTime domainFrequency domain~ 35 m~ 45 m~ 30 m~ 30 mρ1ρ0TransmitterFigure 1.1: Simplified helicopter airborne electromagnetic setup for a genericfrequency and time-domain system. ρ0 is the background resistivity andρ1 represents the target resistivity.∇× e+µ ∂h∂ t= 0 (1.1)∇×h−σe = s (1.2)5TimeNormalized current (A)a) Timeb)Normalized current (A)Figure 1.2: Generic frequency and time-domain waveforms used for an air-borne electromagnetic survey. a) Frequency domain. b) Time domain.where e is the electric field, h is the magnetic field, µ is the magnetic permeabil-ity, σ is the electrical conductivity and s is a source vector. Electrical resistivityρ , the reciprocal of conductivity, is mentioned interchangeably with conductivitythroughout this dissertation.The same equations can be written in the frequency domain as∇×E+ iωµH = 0 (1.3)∇×H−σE = S. (1.4)whereω is the angular frequency and the capitalized notation represents the frequency-domain version of the fields and source term. The relationship between fields andfluxes is summarized by constitutive equationsB = µH (1.5)J = σE (1.6)where B is the magnetic flux density and J is the current density vector. In thiswork we presume that µ = µ0, the magnetic permeability of free space (4pi ×10−7 N/A2).From Maxwell’s equations, the skin depth (δ ) can be derived, which is an6informative metric that describes the depth at which point the amplitude of theelectromagnetic wave is reduced by a factor of 1 / e, and it is written asδ =√2µ0σω. (1.7)The equivalent in the time domain is the diffusion distance (d) that describes thedepth at which point the amplitude of the electromagnetic wave is at a maximumfor a fixed time t, and it is written asd =√2tµ0σ. (1.8)1.4 Electromagnetic Modeling and InversionThis work builds upon techniques for discretizing Maxwell’s equations for nu-merical modeling and inversion of electromagnetic data as shown in Haber et al.(2007b); Haber and Schwarzbach (2014); Oldenburg et al. (2013); Ruthotto et al.(2016). Namely that the forward problem is written asF(m) = d (1.9)where F is a forward operator of Maxwell’s equations, m is the relevant physicalproperty (i.e. electrical conductivity or resistivity) and d is the data. The choiceof discretization is a finite-volume approach, however finite-difference (Commerand Newman, 2004; Egbert and Kelbert, 2012; Weiss and Constable, 2006), finite-element (Key and Weiss, 2006; Li and Key, 2007; Schwarzbach and Haber, 2013)and integral equation (Cox et al., 2010; Zhdanov, 2010) techniques are also com-monly employed. The inverse problem is then written as the optimization ofargminmφ(m) = φd(m)+βφm(m) (1.10)7where φd is the data misfit, φm is the model regularization term, and β is a trade-offor regularization parameter. The data misfit is defined by a least-squares measureφd(m) = ‖Wd(dpred−dobs)‖22 (1.11)where Wd is a diagonal matrix containing the reciprocal of data error standarddeviations, dpred is the predicted data computed from the forward operator F(m),the observation data is given by dobs and ‖ ‖22 is the squared `2 norm. The defaultmodel regularization term is given byφm(m) = αs||Ws(m−m0)||22+αx||Wx(m)||22+αy||Wy(m)||22+αz||Wz(m)||22(1.12)where the α values are user-defined weights that control the proximity of the modelm to a reference model m0 and the overall smoothness of m. Diagonal weightingmatrices W act on particular model cells or cell boundaries and can be used toincorporate a priori information into the inversion process. In Equation 1.12, m0may also be included in the W(m) terms as W(m−m0) if desired. Although `2norms are common in the model regularization, `1 and approximations to `0 typenorms can be used as well to encourage compactness (Ekblom, 1973; Fournier,2015; Last and Kubik, 1983; Portniaguine and Zhdanov, 1999). The optimizationproblem is solved with an iterative Gauss-Newton approach, and the minimiza-tion of the objective function at the (i+ 1)th Gauss-Newton iteration requires thesolution of(JTWTd WdJ+βRm)δm =−JTWTd Wd(d(i)pred−dobs)−βRm(m(i)−m0) (1.13)where J is a Jacobian matrix of sensitivities, Rm is a regularization term com-posed of Ws, Wx, Wy and Wz matrices and δm is a model perturbation vector.The inversion model is updated through Equation 1.13 until a suitable step can nolonger be found or the target data misfit has been reached. In this thesis, the datamisfit will be normalized by the number of data points, meaning the target datamisfit will always be equal to unity.81.5 Airborne Electromagnetic Interpretation TechniquesAs AEM techniques have been around for decades, there have been a variety ofmethodologies developed to interpret and model AEM data sets, and these methodsare now discussed.1.5.1 Semi-quantitative methodsPrior to inversion, analysis of AEM data relied on semi-quantitative approaches,where approximations and simplifications to the structure of the earth are used.These techniques can produce valuable inferences and range from curve-fitting ofsimple geometric shapes such as vertical ribbons representing thin plates to de-termine the conductance of an anomaly (Palacky and West, 1973), to calculatinga best-fitting apparent resistivity value at each transmitter/receiver pair (Fraser,1978). Later on, groups started developing efficient conductivity-depth transformsfor imaging anomalies at depth (Eaton, 1998; Macnae et al., 1998; Macnae andLamontagne, 1987; Macnae et al., 1991) as well as analyzing the time-constantexponential rate of decay of the fields to gain knowledge of the conductivity andshape of the target (Macnae et al., 1998; Nabighian and Macnae, 1991). Whilesimple and fast, semi-quantitative methods have limitations due to their simplicityand various approximations, but they paved the way for electromagnetic inversion.1.5.2 Voxel inversion methodsAEM inversion began with 1D layered earth inversions (Constable et al., 1987;Farquharson and Oldenburg, 1993; Raiche et al., 1985), and moved to 2D (Wilsonet al., 2006), before 3D voxel AEM inversions became computationally feasible(Cox et al., 2010, 2012; Haber and Schwarzbach, 2014; Oldenburg et al., 2013;Ruthotto et al., 2016; Yang et al., 2014). The term 3D voxel inversion means thatthe earth is discretized into cells in a 3D mesh, where the conductivity value in eachcell is solved for during the inversion process. This technique, sometimes referredto as a pixel method, provides flexibility in the sense that any arbitrary anomalyshape or geologic structure that can be discretized on a mesh can be theoreticallyrecovered by the inversion.In a 3D AEM setting, the system of equations to solve can become excessively9large for a number of reasons. The vast area covered in an airborne survey meansthat a significant number of 3D mesh cells are required to discretize the expan-sive spatial domain. In addition, the number of source locations can easily exceed100,000, which makes for an expensive system to solve at many source locations.To tackle these numerical obstacles, Cox et al. (2012) introduces the concept ofa moving sensitivity footprint to reduce the spatial domain size of the forwardproblem. Alternatively, an adaptive mesh refinement can be used (Key and Weiss,2006), where the mesh is refined in areas that require a higher numerical accu-racy, while retaining a coarse solution in other areas. This approach is applied inHaber and Schwarzbach (2014), where the full spatial domain is modeled, but theproblem size is reduced by using decoupled forward and inverse adaptive octreemeshes, specifically for time-domain problems. Recently, Ruthotto et al. (2016)introduced the operators necessary to model the full spatial responses in both thefrequency and time domain through the Julia language framework (Bezanson et al.,2012). It is worth noting that a potential drawback of the voxel approach in largesurveys is that it can be difficult to accurately fit small-scale features, especiallythin anomalies with high resistivity contrasts compared to the background. Strate-gies to address this issue will be discussed in the parametric section.1.5.3 Cooperative and joint inversionThe inversion methodology becomes ambiguous when multiple data sets are presentover a common spatial area. If individual inversions produce inconsistent modelsthis can lead to the question of why the models are different, and which one shouldbe used for ensuing interpretation purposes. In a cooperative inversion, a commonphysical property model is sought by incorporating the inversion result from onedata set as a priori knowledge for another data set, either through a reference model,initial model or inversion weights (Commer and Newman, 2009; Lines et al., 1988;Oldenburg et al., 1997). A variation of the cooperative approach is the alternatingdirection method of multipliers (ADMM) (Wahlberg et al., 2012) where the objec-tive function is re-written as an optimization problem with a Lagrangian multiplierand a constraint that the resistivity models from each data set must be equal. Anapplication of ADMM to 3D hydrogeophysical inversion can be found in Steklova10and Haber (2016). These approaches differ from a joint approach, where multi-ple data sets are inverted simultaneously to produce one inversion model (Albouyet al., 2001; Gallardo and Meju, 2004; Oldenburg et al., 1997; Sosa et al., 2013;Vozoff and Jupp, 1975). This requires being able to forward model and computesensitivities for all data within a single code, and it also requires proper relativeassignments of uncertainties.1.5.4 Constrained inversionPhysical property values can be extracted from drilling or surface rock samples, andrepresents an efficient manner to incorporate geologic knowledge into the inversionas constraints. Constrained inversions have been studied at length and the readeris advised to peruse Lelie`vre and Oldenburg (2009); Lelie`vre et al. (2009); Li andOldenburg (2000); Williams et al. (2009) for a detailed account of various waysto add such information. In this thesis, geologic constraints will be added in theform of reference models and upper and lower physical property bounds when thisinformation is available.1.5.5 Parametric methodsParametric inversion methods can be advantageous when prior geometric or geo-logic knowledge suggests that the target anomaly has a particular shape or a sim-plistic geometry. Parametric inversions solve for a reduced set of parameters thatdescribe the physical property space instead of solving for the value in every meshcell. Early parametric methods incorporated a reduced-physics approach wherevarious assumptions and approximations to Maxwell’s equations were used. Areduced-physics approach was applied to the AEM problem by Keating and Cross-ley (1990) where a plate parameterization in free space is assumed, by Xiong andTripp (1995) with the software package Marco Air, where 3D prisms are placedin a layered-earth environment, and by Raiche (1998) with the package LeroiAir, where 3D plates are modeled in a layered earth setting. Currently, a com-monly used reduced-physics forward modeling and inversion parametric algorithmis Maxwell EMIT (2005), where the response from one or multiple plates is calcu-lated in free-air, in a half-space or in a layered-earth media using electrical ribbon11currents located on plate edges. Although these parametric programs are highlyuseful for modeling responses from simplified conductors, they contain significantnumerical approximations that can create issues when dealing with multiple tar-gets or complex geometries. One of the major contributions of this thesis is thatthe parametric algorithm that is introduced models the full 3D Maxwell’s equationsshown in Equations 1.1 and 1.2.1.6 Thesis OutlineThis thesis is composed of four main chapters, the first two have been publishedin peer-reviewed journals, the third is in preparation for submission and the fourthis part of ongoing research. The content is written under the common theme ofdeveloping practical cooperative and parametric strategies to improve the accuracyof 3D AEM inversion in a variety of geologic settings where conventional tech-niques face difficulties. Each chapter contains both synthetic and field examplesto demonstrate how the approaches are both theoretically and practically viable.These novel algorithms focus on two common scenarios, and the goals of this dis-sertation are listed below.1. To develop strategies to invert AEM data in 3D that spatially overlap withother electromagnetic surveys.2. To develop methods to invert AEM data in 3D to recover thin, high contrastanomalies.Chapter 2 focuses on the first topic where I propose a cooperative methodol-ogy for incorporating multiple electromagnetic data sets over a common area. InChapter 3 I look at the second goal where I introduce a new parametric method ofinverting AEM data for narrow, highly conductive targets. In Chapter 4 I expandthe method from Chapter 3 to allow for multiple anomalies, making the approachsuitable for a larger range of geologic scenarios. Finally in Chapter 5 I address thesituation where both topics apply, and here I modify the cooperative strategy fromChapter 2 to incorporate both parametric and voxel based inversion algorithms.Collectively, my specific contribution to science includes:121. Presenting a new cooperative algorithm for spatially overlapping electro-magnetic data sets that is also general for other types of geophysical data.2. Developing the first parametric algorithm for airborne electromagnetics thatmodels the full 3D Maxwell equations.3. Introducing the first parametric hybrid inversion for geophysical applica-tions where smooth features are incorporated along with compact parametricanomalies within one framework.13Chapter 2Cooperative Inversion of MultipleElectromagnetic Data SetsIn this chapter, I develop a cooperative methodology for inverting AEM data thatspatially overlaps with other electromagnetic surveys. The approach produces aconsistent 3D resistivity model with improved resolution compared to invertingeach data set independently, and adheres to geologic knowledge from boreholemeasurements through bound constraints. Synthetic and field data from the Anto-nio gold deposit in Peru are used to demonstrate the benefits of this technique.2.1 IntroductionConventionally, an individual AEM survey is inverted to create a single resistivityinversion model; however, when multiple spatially overlapping surveys, possiblyfrom different time periods, produce inconsistent inversion models, this can leadto difficulties in interpretation. These model discrepancies suggest that a joint orcooperative approach, where one consistent inversion model is sought, could bebeneficial.The advantage of a cooperative approach is that individual algorithms, tailoredto inverting a particular data type, can be used. This is beneficial because carry-ing out inversions of multiple data sets individually is generally much faster thaninverting them simultaneously. Cooperative inversion is also less sensitive to prob-14lematic or noisy measurements compared to joint inversion. Erroneous data withinthe entire suite of observations can cause a joint inversion and to a lesser degree acooperative approach to proceed very slowly or to stall, and may produce unwantedartifacts in the final model. By breaking up the inversion into small steps using onesingle data set at a time, as in the cooperative approach, the inversion has an easiertime finding a suitable model update step. In a joint inversion this computed stepis based on all data sets, each with a different noise signature, and can be moreproblematic to obtain.There are numerous ways to implement a cooperative inversion through the useof reference models, constraints, or weightings in the regularization term of the ob-jective function (Commer and Newman, 2009; Lines et al., 1988; Oldenburg et al.,1997). However, there is also the question of whether one such consistent modelexists due to the presence of induced-polarization effects or a frequency-dependentconductivity. Therefore, a cooperative inversion needs a set methodology or workflow to address specific issues in implementation. In this chapter, time-domainAEM data are cooperatively inverted with a joint CSAMT and DC resistivity dataset. Here, the DC data are modeled as a very low frequency EM survey in orderto model simultaneously with the CSAMT data in a generalized frequency-domaininversion code. This chapter focuses on three primary research objectives.• To develop a cooperative inversion method for AEM data that spatially over-lap with other geophysical EM data sets.• To prove, with a synthetic example, that a cooperative method increases theaccuracy of the resulting 3D resistivity model.• To apply a cooperative approach to field data to produce a consistent inver-sion model where additional geologic interpretations can be established.Three spatially overlapping data sets are first introduced over the Antonio high-sulfidation epithermal gold deposit: time-domain AEM, CSAMT, and DC resistiv-ity. As an initial step, each survey is inverted individually in 3D to estimate a re-sistivity structure sensitive to that particular survey. Similarities and differences inthe resulting models are then noted. The inversion results coupled with geologicalinsights regarding the deposit are subsequently used to construct a synthetic model15that emulates, as close as possible, the Antonio deposit. Simulated data sets, us-ing field measurement locations, are computed over the synthetic model. Syntheticdata are inverted individually to ascertain whether differences again occur betweenresulting models. At this point, a cooperative approach is applied, which aims toimprove the accuracy of the synthetic inversion model in a qualitative and quantita-tive sense. This cooperative work flow is then applied to the Antonio field data. Aconstrained cooperative inversion follows where bound constraints from boreholesmeasurements are used, which yields a consistent final model from which geologicinterpretations can be made.P E R UAntonioFigure 2.1: Map of Peru with location of Antonio deposit marked as a blackstar. Modified from McMillan and Oldenburg (2014).2.2 Antonio Geologic BackgroundThe Antonio gold deposit is located in the Andes mountains of Northern Peru,as shown in Figure 2.1, and resides within the larger Yanacocha high-sulfidationepithermal gold system. Newmont Mining Corporation owns the majority of thisactive mining and exploration project. The region experienced pervasive hydrother-mal alteration to form a zone of massive silicic alteration in the innermost zone,flanked by alunite, pyrophyllite, kaolinite and montmorillonite assemblages withan outermost halo of propylitic alteration (Teal and Benavides, 2010). This alter-ation zonation is characteristic of high-sulfidation deposits (Arribas, 1995). Fig-16ure 2.2 shows a regional view of the alteration zones at Yanacocha, and the corre-sponding relationship between high grade gold deposits in black and silica alter-ation in red.Figure 2.2: Yanacocha regional alteration map with location of major golddeposits. (modified from Hoschke (2011)).It is clear that the bulk of the gold mineralization resides within or near massivesilica, vuggy silica and granular silica units. The resistive nature of silica alterationcompared to the relatively conductive background makes it an applicable target forEM surveys (Goldie, 2000; Hoschke, 2011; Oldenburg et al., 2005). These quartz-rich areas of metasomatism are often found near faults where confined fluid flowoccurred (Teal and Benavides, 2010). Furthermore, the intersection of faults, of-ten conductive in nature, where hydrothermal breccias broke through overlayingvolcanic units is especially prospective for gold mineralization. These structuraltraps within favorable pyroclastic lithologies, such as ignimbrite beneath or prox-imal to flow dome complexes, are typical geologic hosts to gold deposits in theregion (Loayza and Reyes, 2010). Much of the surrounding clay alteration, which17is generally more conductive than silica alteration (Goldie, 2000), at Antonio hasalso been enriched with silica. For the sake of this study, any zone significantlyenriched with quartz is referred to as silica alteration. Figure 2.3 shows a geologicmap of the Antonio region with lithologies and structures marked. The primaryextent of near-surface (0 - 100 m depth) silica alteration is outlined in dashed red.X (m)Y (m)Figure 2.3: Antonio geology with lithologies, structures and near-surface sil-ica alteration outlined.182.3 Antonio Geophysical Data Sets2.3.1 Time-domain AEMIn 2003, a helicopter based time-domain AEM survey was collected using theNEWTEM I system (Eaton et al., 2013). Five East-West lines over the Antoniodeposit, extracted from the larger airborne survey, are analyzed in this chapter.The peak current of the transmitter was 245 A, with a transmitter loop area of 289m2. z-component responses ( ∂bz∂ t ), measured out to 6.35 ms after current shut-offwere recorded every 20 m with a line spacing of 200 m for a total of 268 trans-mitter positions. Time channels from 30 µs to 2000 µs are used for analysis, andthe discretized waveform is shown in Figure 2.4 with logarithmically spaced timechannels plotted for reference. The waveform closely approximates a square wave,which is advantageous for near-surface mapping of weakly conductive or resis-tive targets (Allard, 2007). Smaller discretized time steps are taken near the firstrecorded time channel for numerical accuracy purposes (Haber and Schwarzbach,2014). Due to mountainous terrain, the drape of the transmitter above the groundvaried from 32 m to 146 m, with a mean value of 62 m. Exact system locations areshown as orange circles in Figure 2.5.Time (s)Normalized current (A)Time (s)Figure 2.4: NEWTEM-I discretized waveform in red with measured timechannels marked as black crosses.192.3.2 CSAMTIn 2003, an asynchronous scalar CSAMT (Zonge and Hughes, 1991) survey wasacquired by Quantec Geoscience with a total of five East-West (EW) and eightNorth-South (NS) lines over the Antonio deposit. Two transmitters: an EW ori-ented transmitter 6.2 km to the south, and a NS oriented transmitter 5.9 km to theeast provided the source for EW and NS lines respectively. Having transmitterswith orthogonal orientations permitted the earth to be energized from two direc-tions. Line spacings varied between 150 m and 200 m, and stations were spread 50m apart. In-line electric field, and orthogonal magnetic field measurements from11 frequencies ranging between 2 Hz and 2048 Hz are used in this study. Receiverlocations are shown as red stars in Figure 2.5.2.3.3 DC resistivityIn 1998, a conventional in-line time-domain pole-dipole DC resistivity survey wasacquired with five EW lines spread 150 m to 200 m apart with 50 m spaced dipoles,and a transmitter/receiver separation (N spacing) of 1 to 6. In 2004, one additionalline of 150 m spaced dipoles situated 50 m south of the previous survey was col-lected. Receiver locations are shown as green triangles in Figure 2.5. Examples offield data from each survey can be seen in Figure 2.6.2.4 Inversion Methods2.4.1 Inversion preparationPrior to inverting field data, initial steps such as quality control of the data, prepar-ing suitable meshes, and assigning appropriate uncertainty values, need to be com-pleted, as well as removing field data below an estimated noise threshold. Eventhough inversions are carried out separately in a cooperative inversion, it is de-sirable to have a common mesh so that difficulties in transferring results fromone mesh to another are avoided. Here, because there are frequency and time-domain surveys, the design of the mesh is based on diffusion distance and skindepth (Nabighian and Macnae, 1991; Ward and Hohmann, 1988) shown in Equa-tions 1.7 and 1.8.20X (m)Y (m)Figure 2.5: Antonio geophysics locations. Receiver sites for AEM (orangecircles), CSAMT (red stars), and DC Resistivity (green triangles) over-laid on topography with near-surface silica alteration outlined in dashedred and geologic faults in thin blue.In order to minimize numerical accuracy issues, mesh creation in this chap-ter and throughout the thesis adheres to the guidelines below, which are based onexperience working with EM inversion codes. Additionally, when significant to-pography, such as at Antonio, is present, the topography discretization error mustbe taken into consideration when determining a finest mesh cell size.• The mesh contains padding cells extending to a minimum of twice the largestdiffusion distance, or twice the largest skin depth. For a 50 Ωm half-space21a) b) c)Predicted DataPredicted DataAmplitude EyPredicted DataObserved Data Observed DataObserved DataPredicted DataX (m)X (m)X (m)X (m)Y (m)Y (m)Y (m)Y (m)X (m)X (m)(ohm-m)(ohm-m)(V / Am)(V / Am)(V / Am2)(V / Am2)AEM CSAMTDCFigure 2.6: Observed and predicted field data from individual inversions. a)AEM plan view of ∂bz∂ t at 139 µs b) CSAMT plan view of y-componentelectric field amplitudes at 16 Hz c) DC apparent resistivity pseudo-section from line 9230630.this equates to roughly four kilometers of padding around the core area ofinterest for AEM, and five kilometers of padding around the CSAMT re-ceivers.• The smallest cell size in the mesh is at most half the minimum diffusiondistance or minimum skin depth. For a 50 Ωm half-space this produces aminimum cell size of 24 m and 40 m for AEM and CSAMT respectively.In an inversion, there is always a trade-off between accuracy and total size of themesh. Using the above guidelines, the final mesh is designed to be suitable for allthree surveys, and contains core cell sizes of 25 m × 50 m × 25 m in x,y,z. Toaccommodate CSAMT transmitters far away from the receiver area without greatlyincreasing the number of total cells, an adaptive ocTree mesh is used (Haber et al.,2007a) with a total number of inversion mesh cells of 86,942. For AEM inversions,forward and inverse meshes are decoupled for maximum efficiency as described inHaber and Schwarzbach (2014), and these forward meshes have roughly 10,000cells per mesh.The next task is to assign uncertainties to the data. The assigned error uncer-22tainty is made up of a percentage and a floor value. Percentages range between10% and 15%, depending on the noisiness of the data, and noise floors are selectedas shown in Table 2.2.2.4.2 Inversion methodologyAEM, CSAMT and DC 3D inversions follow algorithms outlined in (Haber et al.,2004; Haber and Schwarzbach, 2014; Oldenburg and Li, 1994). In this chapter,three Gauss-Newton iterations are computed for each trade-off value of β , knowncollectively as one β iteration, and the inversion terminates when the data misfitreaches its target level φ ∗d of unity.2.4.3 Individual Antonio inversions3D inversions for all individual field data sets are performed, and examples ofpredicted data results are shown in Figure 2.6. Resistivity plan maps at a constantelevation of 3870 m above sea-level through the resulting models are shown inFigure 2.7 with a consistent color scale. Due to rolling topography, a constantelevation slice of 3870 m corresponds to an average depth below surface of roughly75 m, although it ranges from 10 m to 150 m throughout the survey area. Additionalconstant elevations slices at 3895 m, 3845 m, 3820 m, and 3795 m are displayedin Figure 2.16. These inversions are unconstrained, and use a 50 Ωm half-spacereference model. Details about uncertainty assignments and final data misfits aresummarized in Table 2.1. Geologic faults and the outline of known near-surfacesilica alteration are plotted for reference. Figure 2.7a shows a slice through the3D AEM inversion model along with data locations, and the result recovers a largeuniform resistor in the center of the survey area. This anomaly agrees well with theknown resistive silica alteration outline and past studies (Oldenburg et al., 2005,2004), although the resistor extends beyond the known outline to the north-westfor approximately 200 m.23c)a) b)d)X (m) X (m) X (m) X (m) Y (m)Y (m)Y (m)Y (m)Figure 2.7: Resistivity 3D inversions from field data at a 3870 m elevationslice with receiver locations shown. Near-surface silica alteration out-line in dashed red and geologic faults in thin solid blue. a) AEM b)CSAMT c) DC Resistivity d) Cooperative.24Table 2.1: Error assignments and final data misfits for individual field andsynthetic inversions.Data Set Error (%) Error (floor) Final φdFieldAEM 15 4.2×10−10 V/Am2 0.98CSAMT 10 3.5×10−8 V/m (E) 0.962.7×10−8 A/m (H)DC 10 0.5 mV 0.91SyntheticAEM 10 4.2×10−10 V/Am2 1.00CSAMT 10 1.0×10−9 V/m (E) 0.521.0×10−9 A/m (H)DC 10 1.5 mV 0.7225Figure 2.7b portrays a 3870 m constant elevation slice through the 3D CSAMTinversion result along with receiver locations. Due to the asynchronous nature ofthe CSAMT survey, meaning the transmitter and receiver clocks are not synchro-nized, only the amplitude and not the phase of the electric and magnetic fields areinverted. The imaged resistor has some similarity to that from the AEM inver-sion, but it deviates more from the known silica alteration outline. Although therecovered resistivity magnitudes compare well with that from AEM, the center ofthe CSAMT anomaly is shifted to the west, and the resistor extends 100 m to 200m west of the mapped alteration zone. The CSAMT anomaly is broken up intomultiple pieces, unlike the cohesive AEM image, and some small conductive areasoccur within the resistive region. There is also a spurious anomaly in the extremenorth-west corner of Figure 2.7b, which is outside the data supported region andshould be ignored.The inversion of DC potential differences is shown as a 3870 m constant ele-vation slice in Figure 2.7c, along with receiver electrode locations. The recoveredmodel exhibits a similar curved north-northwest trending resistivity feature as theAEM survey, although the strongest resistor occurs outside the dashed outline tothe north-west. Within the known alteration zone, much of the area is modeled asconductive, although there is a thin resistive feature that extends down through themarked anomaly to the south-east portion of the survey. Resistivity magnitudes aregenerally weaker compared to AEM and CSAMT results, and the overall shape isnarrower and smaller compared to the other two models.2.4.4 Cooperative inversion work flowDiscrepancies between individual field inversion results suggest that joint and/orcooperative methods are warranted to produce one resistivity model that fits all thedata sets. Intrinsically this may be challenging. Each survey has its own noisesignature, and the geometric layout of transmitters and receivers may be sensitiveto different portions of the earth. In addition the process of homogenization, orrepresenting micro-scale features in the earth as a homogeneous physical propertyvalue within a mesh cell, may not yield equivalent results across all electromag-netic techniques (Caudillo-Mata et al., 2016). Hence, each survey will not neces-26sarily image the ground in the same manner. Furthermore, induced-polarizationeffects, anisotropy, data quality variations, acquisition location differences, mod-eling errors, and different source types (galvanic vs. inductive) can all potentiallycomplicate results. That being said, at Antonio, the frequency-dependent natureof conductivity within sulfides is thought to be minimal and there is no reason tobelieve that anisotropy is a major factor. Galvanic and inductive sources will im-age the ground in unique ways, but it is postulated that the combined informationfrom these data sets will be beneficial for inversion accuracy instead of detrimen-tal. Therefore, it is assumed that producing a consistent isotropic resistivity modelis possible, while the mesh generation guidelines above should minimize model-ing errors. First, the AEM data is cooperatively inverted with a joint CSAMT/DCdata set. In this joint CSAMT/DC data set, DC voltages are converted into elec-tric fields, and then treated as a 0.125 Hz frequency input in the CSAMT data.This frequency is sufficiently low that there are no inductive effects observed in thesimulated data.The cooperative work flow developed is proposed in Figure 2.8. Starting withan initial model σ0, the first task is to calculate a model update from data set d1using a starting beta, β (1)1 , where the subscript refers to data set (d1), and the super-script in parenthesis refers to the cooperative iteration number. The Gauss-Newtonsystem is solved multiple times for a given β , known as an inner iteration, to pro-duce a revised model σ (1)1 . This updated model becomes the initial and referencemodel for the first cooperative iteration for data set d2, and the output is σ(1)2 . Sub-sequently, the process repeats and σ (1)2 becomes the initial and reference model fora second cooperative iteration. The values of β for this next, and future, iterationsare reduced according to a schedule β (i+1)1 = γβ(i)1 and β(i+1)2 = γβ(i)2 , where γ ≤ 1.For the work here γ = 0.2, which is an aggressive β reduction scheme; however,experience with the code suggests that solving the Gauss-Newton system multipletimes per β makes this a feasible approach.This cooperative work flow continues up to a maximum number of iterationsimax, until the target misfit is reached for one or both data sets or until a suitablestep that reduces the data misfit can no longer be found. If a single model hits bothtarget misfits, φ ∗d1 for d1 and φ∗d2 for d2, each data set has been adequately fit andthe process stops.27If only one data set is fit, i.e. d1, then the emphasis shifts to continued Gauss-Newton iterations for d2. If the output model σ(i)2 is still compatible with d1 butstill not satisfactory for d2, then further iterations are performed with d2. Afteradditional work with d2, if the misfit for d1 is significantly increased, then thecooperative inversion cycle resumes with d1 starting with the last used β i1. At thispoint, if the misfit for d2 no longer decreases, the current model is accepted as thefinal result. The strategy is automated and requires only a few user parameters.2.5 Synthetic Inversion Results2.5.1 Individual synthetic inversionsIn synthetic modeling, a geophysical data set is computed over a pre-defined resis-tivity distribution, noise is then added to simulate errors encountered in field dataand this noisy synthetic data set is inverted in an attempt to recover the originalmodel. The synthetic model in this section is designed to encapsulate the primaryfeatures of the Antonio area. It includes a 225 m thick 1000 Ωm resistor placed ina uniform 50 Ωm background. Two conductive 10 Ωm blocks with dimensions of100 m × 150 m × 100 m are embedded. The northern conductive block is buried75 m below the surface while the southern block is exposed at the surface. To-pography for the Antonio area is used for the synthetic study. Figure 2.9 shows a3870 m constant elevation slice through this synthetic model. Additional elevationslices, are shown in Figure 2.13a.Forward modeling of AEM, CSAMT and DC surveys is carried out by keepingdata locations and other specifications equivalent to those of the field setups. 10%Gaussian noise is added to the measurements prior to inversion. Error assignmentsand final data misfit values are summarized in Table 2.1. Constant elevation slicesat 3870 m through the resulting inversions are shown in Figure 2.10, with addi-tional elevation slices presented in Figure 2.13. The images all display a resistivebody centered in the correct location, but the extent of the resistor, and the abilityto detect the conductive blocks varies between the three models. In Figure 2.10a,the AEM inversion is able to image the outline of the resistor, but is not able to de-tect either of the two conductive bodies. Figure 2.10b shows the CSAMT inversion28Figure 2.8: Cooperative inversion work flow. σ (i)j = conductivity model fromthe ith cooperative iteration for data set j, β = trade-off parameter, φd =data misfit, φ ∗d = target data misfit.29recovery. The result accurately images the overall resistor geometry, while clearlydetecting the southern conductor, and faintly identifying the northern anomaly. ForCSAMT data, real and imaginary electric and magnetic fields are inverted insteadof amplitude-only data as in the field example, and this produces twice the numberof data points. As field components (real and imaginary) represent a typical dataset collected by industry, these are used to showcase the maximum improvementthat can be obtained with the cooperative inversion. Experience has also shown thatusing field components can sometimes be better at recovering the true resistivitystructure compared to amplitude data only. Numerical computations are also easierwith field components as the equations are more linear compared to using ampli-tudes and phases, and this reduces the chance of being stuck in a local minimumduring optimization. Error floors for synthetic CSAMT data are also set consider-ably lower to prevent dramatic under fitting of the data. The DC inversion resultis displayed in Figure 2.10c, which depicts the overall geometry of the resistor andthe southern conductive block. Resistivity magnitudes are closer to those in thesynthetic model compared to the CSAMT and AEM results; however, the northernblock is not seen. The three inversions have each produced valuable information,but significant differences exist in the models. At this point, the cooperative workflow is tested on these data sets in the hopes of establishing a common resistivitymodel.2.5.2 Synthetic cooperative inversionThe cooperative work flow from Figure 2.8 is applied to the synthetic data sets.A common mesh is used for all inversions, and previous error assignments arekept from the individual synthetic inversions with the addition of an error floor of3.0×10−5 V/m for electric fields at 0.125 Hz converted from DC voltages. A con-vergence curve documenting the data misfit progression for each data set is plottedin Figure 2.11, and summarized in Table 2.2. For each cooperative iteration thereare a maximum of four data misfit evaluations, representing the initial data misfit,and the resulting misfit after each of the three inner iterations. It can be seen on Fig-ure 2.11a that the joint CSAMT/DC data, pictured as red x’s, hits the target misfit,shown by a black dashed line, on the third cooperative iteration. On the fourth co-3050 Ωm10 Ωm1000 ΩmX (m)Y (m)Figure 2.9: Synthetic model resistivity at a 3870 m elevation slice.operative iteration, the initial data misfit for the joint data set is still below the targetlevel, and thus no additional model update is performed. Eventually the process isstopped after iteration 7, because the AEM data misfit increases compared to thatin iteration 6, indicating the inversion is trending in the wrong direction. Therefore,the model after iteration 6 of the AEM inversion is chosen as the final result. At thatpoint, the CSAMT/DC and AEM final data misfits are very close to their desiredtarget misfits. Figure 2.10d shows a 3870 m constant elevation slice through thissequential synthetic inversion model. The resistor magnitude and shape is definedmore uniformly compared to the individual inversions, and the northern conductoris now better detected. The southern conductor is still clearly defined, althoughnot as strongly as in the individual CSAMT or DC inversion. Although the AEMindividual inversion does little to resolve the two conductive targets, it contributesto the cooperative inversion by helping to define the main resistive target. As AEMis an induction method, it should be sensitive to conductive targets; however, thenorthern conductor is buried, which masks the signal. Furthermore, the flight lines31b)c)a)d)X (m) X (m) X (m) X (m) Y (m)Y (m)Y (m)Y (m)Figure 2.10: Resistivity 3D inversions from synthetic data at a 3870 m eleva-tion slice with receiver locations shown. True resistor outline in black.a) AEM b) CSAMT c) DC Resistivity d) Cooperative.32Table 2.2: Final data misfits for cooperative synthetic and field inversions.Inversion Method Data Set Final φdSynthetic Cooperative AEM 1.01CSAMT/DC 1.02Field Cooperative AEM 1.62CSAMT/DC 1.15Synthetic Constrained Cooperative AEM 1.00CSAMT/DC 1.00Field Constrained Cooperative AEM 1.29CSAMT/DC 1.15do not directly pass over the southern conductor. In a resistive setting, the sig-nal level is also smaller, and the conductive responses may be washed out by theGaussian noise added, and by the relatively high uncertainty assignments.Visual interrogation of the models shows that the cooperative result is betterthan any of the individual inversions. In an attempt to quantify this numerically,a region of the model is extracted that consists of the resistor and its includedconductive blocks. The recovered models are numerically evaluated to determinehow close they are to the true model using a residual (R) as the metric of choice asshown in Equation 2.1.R =1N‖log10(m)− log10(mtrue)‖22 (2.1)where m and mtrue are the recovered and true resistivity model values respectively,and N is the number of cells in the volume of interest. A lower value of R refersto a smaller deviation from the true model, and hence a more accurate recovery.The residuals corresponding to the individual and cooperative inversions, and alsothe residual from the starting model, a 50 Ωm half-space, are shown in Table 2.3.The cooperative inversion performs the best, followed closely by the DC result,while the CSAMT and AEM models fare increasingly worse. As expected, allfour inversions recover a more accurate model compared to a uniform 50 Ωm half-space. This synthetic example demonstrates that a cooperative method improvesthe accuracy of the recovered anomalies. Because of its close association with the33a)c) d)b)Cooperative iterationData misfitData misfitCooperative iterationData misfitCooperative iterationData misfitCooperative iterationFigure 2.11: Data misfit convergence curves with target misfit shown indashed black with final models circled in black. a) Synthetic data with-out bounds b) Field data without bounds c) Synthetic data with boundsd) Field data with bounds.field data example, the same conclusions are anticipated there.2.5.3 Synthetic constrained cooperative inversionThus far, all inversions have been unconstrained, and all start with a 50 Ωm half-space initial and reference model. To test a constrained approach, the syntheticexample is revisited. Constraints are extracted from drilling information from 78synthetic boreholes which are replicas of the field boreholes over the Antonio de-posit. Each hole intersects multiple discretized cells in the synthetic 3D model, andthis information is used to produce upper and lower resistivity bounds. The upper34Table 2.3: Quantitative assessment of synthetic inversions using a residual(R) value. A lower residual refers to a more accurate recovery.Inversion Model Residual (R)50 Ωm half-space 1.62AEM 0.80CSAMT 0.66DC 0.55Cooperative 0.47Bounded Cooperative 0.26and lower bounds for each model cell intersected by a borehole are respectively setto 5% above and below the corresponding resistivity value in the synthetic model.The regularization of the inversion algorithm attempts to smooth this informationaway from constrained cells, and a global bound range of 1 Ωm to 5000 Ωm is ap-plied to all cells not intersected by boreholes. This global bound prevents extremeinversion artifacts from emerging, and allows the conductivity to vary freely withinthis range.The cooperative method from Figure 2.8 is implemented. The convergence plotfrom Figure 2.11c shows that after the third cooperative iteration, the CSAMT/DCdata set reaches its target value, and after the sixth AEM cooperative iteration,both data sets reach convergence. Therefore, the addition of bounds helps guidethe cooperative approach to a single solution that fits both data sets. A 3870 m con-stant elevation image of the constrained inversion result is shown in Figure 2.12awith additional elevation slices in Figure 2.13f. The outline of the resistor is muchimproved compared to previous results, and both the northern and southern con-ductors are clearly imaged. Even the small resistive zone at the southern tip of thesurvey is recovered, thanks primarily to a synthetic borehole. Quantitatively, theconstrained cooperative method produces a residual value of 0.26, which outper-forms all other synthetic inversions.35a)a) b) X (m)X (m)Y (m)Y (m)Figure 2.12: Cooperative 3D inversions with bounds at a 3870 m elevationslice a) Synthetic data with the true resistor outline in thin black b)Field data with the outline of known near-surface silica alteration indashed red, interpreted silica alteration outline at 3870 m elevation insolid black, geologic faults in thin blue, and conductive anomalies ofinterest numbered in yellow.2.6 Antonio Inversion Results2.6.1 Antonio cooperative inversionThe cooperative approach is applied to the field data, once again using the samecommon mesh and error assignments as before. A noise floor of 1.0× 10−5 V/mis placed on 0.125 Hz electric fields converted from DC data. The convergencecurves are plotted in Figure 2.11b and final data misfit values are described in Ta-ble 2.2. The joint CSAMT/DC data set reaches its target level, after cooperativeiteration 8, and the inversion terminates after iteration 9, when the AEM data mis-fit is greater than after iteration 8. The model after the eighth AEM iteration is36a) Resistivity: True model b) Resistivity: AEM c) Resistivity: CSAMTd) Resistivity: DC e) Resistivity: Cooperative f) Resistivity: Cooperative with bounds3895m3870m3845m3820m3795m775400X (m)923000092304009230800Y (m)Elevation (m)775800 7762003895m3870m3845m3820m3795m775400X (m)923000092304009230800Y (m)Elevation (m)775800 7762003895m3870m3845m3820m3795m775400X (m)923000092304009230800Y (m)Elevation (m)775800 7762003895m3870m3845m3820m3795m775400X (m)923000092304009230800Y (m)Elevation (m)775800 7762003895m3870m3845m3820m3795m775400X (m)923000092304009230800Y (m)Elevation (m)775800 7762003895m3870m3845m3820m3795m775400X (m)923000092304009230800Y (m)Elevation (m)775800 776200Figure 2.13: Resistivity 3D inversions from synthetic data at five elevationslices: 3895 m, 3870 m, 3845 m, 3820 m, 3795 m. a) True model b)AEM c) CSAMT d) DC Resistivity e) Cooperative f) Cooperative withdrilling bounds.chosen as the end result. The final AEM data misfit of 1.62 is above the targetvalue of 1.0, which demonstrates that the inversion has a more difficult time fittingAEM measurements to the assigned error levels, compared to CSAMT/DC data,where the final data misfit is only slightly higher than the target value. A constant3870 m elevation image of the result is displayed in Figure 2.7d. The recoveredresistive anomaly has a better agreement with the silica alteration outline comparedto the individual inversions and additional conductive features are clearly visiblewithin the deposit region. Comparable to the synthetic case, the AEM contributesmost to the cooperative result by mapping the extent of the large resistor, whilethe conductive features are largely due to the CSAMT/DC data. This may seem37counter intuitive at first, but the NEWTEM waveform is designed specifically tobe sensitive to near-surface resistors, and the CSAMT/DC surveys have more datapoints directly over the well-imaged southern conductor compared to AEM mea-surements. The magnitude of the cooperative resistor is slightly stronger than theindividual inversions, and the resistive anomaly extends past the mapped outline tothe north-west. Erroneous resistive anomalies in the extreme north-west and north-east corner are beyond the data supported region, and should be neglected. Aswould be expected, the cooperative model is similar to all three individual imagesin some aspects, but not identical to any, and is interpreted to be a more accuraterepresentation of the resistivity signature over the Antonio deposit.2.6.2 Antonio constrained cooperative inversionThe constrained cooperative method is now applied to field data, which requiresborehole resistivity information. However, boreholes at Antonio have alterationlogging and geochemical assay values, but only a limited number of resistivitymeasurements. Surprisingly, from the samples available there is a poor correlationbetween alteration or rock type with resistivity; thus, another proxy for resistivityinformation is needed. Nelson and Van Voorhis (1983) noted an inverse relation-ship between total weight percent sulfide and resistivity for in-situ rock measure-ments in a porphyry environment. Although Antonio is a high-sulfidation depositand not a porphyry, this concept is investigated by plotting in Figure 2.14 the totalweight percent sulfur against resistivity for 30 rock samples at Antonio, coloredby alteration type. The linear regression curve, represented by the black line inFigure 2.14, does not include propylitic samples, and has a resulting Pearson cor-relation coefficient of -0.73, indicating a statistical negative correlation. Nearly anidentical relationship is extracted when resistivity is compared to total weight per-cent sulfide, but total sulfur content is used because more of these measurementsare available.This regression relationship is applied to the alteration values from 78 fieldboreholes, whose locations are shown in Figure 2.15a, and total sulfur contentfrom 61 boreholes, displayed in Figure 2.15b. The product of the relationshipis a resistivity reference model, shown in Figure 2.15c. Propylitic samples are380 10 20 30101102103104 Total Sulfur (%) Resistivity (ohm−m) SilicaPropyliticRegression CurveFigure 2.14: Linear regression curve, shown as a black line, for the resistivityversus total sulfur relationship from 30 borehole samples.assigned a resistivity value of 23 Ωm in the reference model, equal to the mean ofthe three samples in Figure 2.15a. To create upper and lower bounds, the mean andstandard deviation of the resistivity values within each model cell are calculated,and the bounds are set respectively to one standard deviation above or below themean cell value. The incorporation of these bounds to a cooperative constrainedfield inversion produces the convergence curve shown in Figure 2.11d. The finaltarget AEM misfit of 1.29 is lower than the unconstrained cooperative inversion,but still above the target level, whereas the joint CSAMT/DC final misfit of 1.15is identical to the unconstrained case. A 3870 m constant elevation slice fromthis model is displayed in Figure 2.12b. This result has a strong agreement withthe known silica outline, and with previous studies (Oldenburg et al., 2005, 2004),although some additional features are present, and will be discussed further in thegeologic interpretation section.Caution must be placed when implementing constraints such as upper andlower resistivity bounds, because they could erroneously bias the inversion model.Field constraints are derived from a relationship between total sulfur content andresistivity based on 30 laboratory rock measurements. This represents a small sam-ple size from which to produce constraints for an entire model. But the validity ofthis relationship is corroborated by the similarity between the unconstrained and39a) c) b) 3895m3870m3845m3820m3795m3895m3870m3845m3820m3795m3895m3870m3845m3820m3795m775400 775800 776200X (m)92304009230800Y (m)9230000Elevation (m)775400 775800 776200X (m)92304009230800Y (m)9230000Elevation (m)775400 775800 776200X (m)92304009230800Y (m)9230000Elevation (m)Figure 2.15: Drilling information for model cells intersected by 78 boreholesat Antonio. a) Alteration b) Total sulfur content c) Resistivity derivedfrom a regression relationship between total sulfur/alteration and resis-tivity.constrained cooperative inversions. For additional confidence in the bounds, morepetro-physical measurements should be acquired, and this is recommended for fu-ture constrained cooperative studies.2.6.3 Geologic interpretationThe constrained cooperative inversion result shown in Figure 2.12b is consid-ered the final field inversion, and additional elevation slices are displayed in Fig-ure 2.16f. This final result compares well with the unconstrained cooperative modelin Figure 2.7d, which adds more confidence to the regression relationship thathelped calculate the upper and lower resistivity bounds. In all elevation slices,the final model maps the target silica alteration zone as a strong resistor. At anelevation of 3870 m, the resistive zone extends past the mapped alteration outlineto the north-west, which suggests that the silica alteration does as well. Conse-quently, an outline shown in solid black in Figure 2.12b represents the interpretedsilica alteration zone at 3870 m elevation based on the final field inversion.Within this resistive alteration region, small conductive areas are recovered,and two such anomalies are numbered 1 and 2 in Figure 2.12b. From borehole alter-ation logs, anomaly 1 represents a small area of propylitic alteration, and anomaly 2coincides with a large total sulfur anomaly within silica alteration, which suggests40a) Alteration b) Resistivity: AEM f) Resistivity: Cooperative with boundsd) Resistivity: DC3895m3870m3870m3845m3820m3795m775400 775800 776200X (m)92304009230800Y (m)Elevation (m)3895m3870m3845m3820m3795m3895m3870m3845m3820m3795m3895m3870m3845m3820m3795m3895m3870m3845m3820m3795m3895m3870m3845m3820m3795m775400923000092304009230800Elevation (m)775400 775800 776200923000092304009230800Elevation (m)775400 775800 776200X (m)923000092304009230800Y (m)Elevation (m)775400776200X (m)923000092304009230800Y (m)Elevation (m)775400 775800 776200X (m)923000092304009230800Y (m)Elevation (m)9230000775800775800 776200**X (m) X (m)Y (m)Y (m)e) Resistivity: Cooperativec) Resistivity: CSAMTFigure 2.16: Resistivity 3D inversions from field data at five elevation slices:3895 m, 3870 m, 3845 m, 3820 m, 3795 m. a) Alteration from bore-holes b) AEM c) CSAMT d) DC Resistivity e) Cooperative. An exam-ple of a conductive anomaly potentially linked to propylitic alterationis highlighted with a yellow star f) Cooperative with drilling bounds.extensive sulfide mineralization. This location of inferred sulfide mineralizationis potentially explained by its proximity to fault intersections directly to the northand west, where fluid flow would have been confined. Anomaly 2 also sits on achargeability high from induced-polarization data. Based on these characteristics,anomaly 2 is a prime target for gold mineralization. Borehole assays corroboratethis observation with anomalous gold values present within anomaly 2 but not inanomaly 1.At greater depths, many conductive anomalies within the confines of silicaalteration can be explained by the presence of propylitic alteration. This agreement41is noted even in the unconstrained cooperative inversion where borehole constraintshave not been applied. An excellent example is the large conductive zone at anelevation of 3795m as shown by yellow stars in Figures 2.16e and 2.16f. Whencompared to the region of propylitic alteration in Figure 2.16a there is a strongspatial correlation. Other smaller conductive anomalies within the silica-rich zonenot linked to propylitic alteration may indicate the presence of sulfides, and theseregions are considered prospective targets for gold mineralization.2.7 ConclusionsThe mineralization at Antonio is known to occur within areas of silica alteration,which is a resistive geophysical target. However, AEM data in addition to twoother geophysical surveys over the Antonio deposit image the resistivity signaturedifferently. The cooperative inversion method inputs all this information into onephysical property model, and images the resistive zone at Antonio with a high levelof agreement to the known silica alteration. This method is successful in ensuringa consistent inversion result, as demonstrated with a synthetic and field example.Through synthetic modeling, the cooperative method is shown to improve theaccuracy of 3D resistivity inversions. Moreover, adding geologic constraints to theinversion, in the form of upper and lower resistivity bounds for each 3D model cell,produces a result closer to the true synthetic model. Consequently, in a field setting,when reliable borehole physical property information exists, it should be incorpo-rated into the 3D inversion. The final constrained cooperative field model clearlyimages the large silica alteration zone, and its shape is in general correspondencewith drilling and previous studies of the area. An area of silica alteration from thefield inversion that extends beyond the previous known outline offers a new inter-pretation of the mapped alteration. Further analysis from this final model highlightspotential areas of sulfide and gold mineralization, within the silica alteration zone,in the form of small conductive anomalies.Collectively this study illustrates that a practical cooperative constrained inver-sion methodology is possible for AEM and other electromagnetic data sets, andthat exploration companies could benefit from such a technique for spatially over-lapping surveys with or without borehole constraints.42Chapter 3Single Anomaly ParametricInversion of AEM dataIn this chapter, I demonstrate a method to invert AEM data using a parametriclevel-set approach combined with a conventional voxel-based technique to forma parametric hybrid inversion. The approach is designed for situations where avoxel-based inversion alone may struggle, such as a conductive thin dipping platewith a large physical property contrast between itself and a highly resistive back-ground. Synthetic and field data from the Caber volcanogenic massive sulfide de-posit in Quebec, Canada, are shown to highlight the advantages of this technique.Although the approach is shown for time-domain data, the same functionality hasbeen developed for the frequency domain and will be showcased in later chapters.3.1 IntroductionIt has been established that AEM surveys are widely used as a mineral explorationtool for imaging subsurface electrical resistivity distributions to help identify rocktypes, mineralization, and alteration zones (Keller, 1988). In a conventional voxel-based electromagnetic inversion, the resistivity value in every active mesh cell issolved for, and in recent years, inversion algorithms have been modified to facil-itate thousands of source locations. This has resulted in numerous 3D inversioncodes tailored toward AEM data (Cox et al., 2012; Haber and Schwarzbach, 2014;43Oldenburg et al., 2013; Ruthotto et al., 2016).When an anomaly of interest has sharp boundaries and a large physical propertycontrast compared to the background, experience with voxel-based inversion codessuggests that accurately defining these abrupt resistivity contrasts can be challeng-ing. Moreover, the true resistivity of the target is often poorly estimated, especiallywith a low resistivity or high conductivity target. One option is to use minimum orcompact norm solutions for voxel-based problems to help achieve sharp anomalies(Farquharson and Oldenburg, 1998; Fournier, 2015; Last and Kubik, 1983; Port-niaguine and Zhdanov, 1999; Zhdanov and Cox, 2013); however, experience hasshown mixed results when applying these types of solutions to 3D AEM problems.Another approach that will be explored further in this chapter is the parametricmethod (Abubakar et al., 2006; Aghasi et al., 2011). Here, instead of solving forthe resistivity in every cell, only a few parameters are sought to describe the phys-ical property space, and the reduction in the number of variables can typically bemany orders of magnitude. Parametric inversions can also be coupled with suchmethods as level sets (Osher and Fedkiw, 2001; Osher and Sethian, 1988), to solvefor the resistivity and shape of a target of interest (Aghasi et al., 2011). Level-set methods have previously been applied to 2D electromagnetic tomography andinverse scattering (Dorn and Lesselier, 2006; Dorn et al., 2000), as well as othergeophysical inverse problems such as 2D travel time (Zheglova and Farquharson,2016), 2D DC resistivity (Ascher and Roosta-khorasani, 2016; van den Doel andAscher, 2006), 3D gravity (Isakov et al., 2011) and 3D gravity gradient data (Luand Qian, 2015).This research focuses on 3D AEM data, and specifically on parametric level-setinversion methods, which adds a level of stability to the inverse problem comparedto traditional level-sets by restricting the shape of the anomaly; thus, acting asa form of regularization (Aghasi et al., 2011). A parametric level-set approachenables the inversion to find optimally shaped anomalies with potentially sharpboundaries and a high resistivity contrast compared to the background with a mini-mal number of variables. This study adds to previous hybrid modeling approachesfor AEM data such as Leroi Air (Raiche, 1998), where thin sheet integral equa-tions are used with 3D plates in a layered earth, and Marco Air (Xiong and Tripp,1995), where volume integral equations are incorporated with 3D prisms in a lay-44ered earth.Finding an appropriate parameterization is critical, and this work chooses towork with a skewed Gaussian ellipsoid to represent the target anomaly, althoughother options are available such as radial basis functions or truncated Gaussian dis-tributions (Aghasi et al., 2011; Pidlisecky et al., 2011). The Gaussian ellipsoidshape is chosen for its flexibility, as it can easily discretize high-frequency featuresincluding plates and thin layers, as well as broader shapes such as large intrusions.It is recognized that the skewed Gaussian parameterization has limitations, espe-cially in the single-anomaly case. Therefore, to account for additional features inthe model space, a voxel-based stage is added to solve for smooth background fea-tures. The choice of parameterization is specific, but the method developed hereis general to any parameterization, and could be applied to other geophysical datasets such as ground EM, potential fields or induced-polarization. The data setspresented in this section are all in the time domain, but the method has also beendeveloped for frequency-domain observations.The three research goals in this chapter are listed below.• To develop a parametric inversion for AEM data with a skewed Gaussianellipsoid parameterization.• To show, with a synthetic and field example, that the novel parametric inver-sion can accurately recover a thin dipping plate; a pertinent target for mineralexploration.• To combine parametric and voxel-based methods to form a hybrid techniquethat can model a single target of interest with parametric inversion, whilefilling in remaining features with a voxel-based code.This chapter first discusses general electromagnetic theory before going intodetail regarding the parametric hybrid methodology. Simulated AEM data overa synthetic model composed of a thin dipping conductor in a resistive half-spaceis then introduced to provide a means to test the code. Following the syntheticexample, the hybrid inversion, which to our knowledge is the first algorithm tocombine smooth and parametric features for a geophysical application, is tested on45AEM field data over a volcanogenic massive sulfide (VMS) deposit. The results aresubsequently compared to previous work and to geologic knowledge from drilling.3.2 Electromagnetic BackgroundThis chapter focuses on AEM platforms where the receiver is contained in thecenter of the transmitter loop, known as coincident loop systems (Allard, 2007).The data collected are typically x,y,z component ∂b∂ t data, although components ofb alone can also be calculated.For conventional voxel-based time-domain AEM inversions, the approach dis-cussed in Haber and Schwarzbach (2014) is followed, where Gauss-Newton basedoptimization is used to solve quasi-static Maxwell’s equations in space and timesubject to boundary and initial conditionsn× (∇×h) = 0 (3.1)h(x,y,z, t = 0) = h0 (3.2)∇ ·µh0 = 0 (3.3)using a finite-volume discretization on ocTree meshes (Haber et al., 2007a). Heren is a normal vector, x,y,z are spatial observation coordinates and t is time. Forfurther information on voxel-based electromagnetic inversion theory, see Haberet al. (2004); Haber and Heldmann (2007); Haber et al. (2007b). Details regardingthe parametric hybrid inversion are now discussed.3.3 Inversion MethodologyConsider an inverse problem where the forward problem has the formF(m(x))+ ε = d (3.4)where F maps the function m(x), with position vector x, to the discrete data d, andε is the noise that is assumed to be Gaussian. After discretization of the model, the46discrete problemF(m)+ ε = d (3.5)is obtained where m is a discrete approximation to the function m(x). A maximumlikelihood approach would minimize the global misfit φ as shown in Chapter 2.minmφ(m) =12(F(m)−d)TWTd Wd(F(m)−d) (3.6)However, this problem is typically ill-posed, and there are many possible solutionsthat minimize the global misfit. To obtain a well-posed problem, two possibleroutes can be taken. First, if the model m does not have any particular form somesmoothness can be assumed and this results in a regularized least-squares approach.A second option is where some specific a priori information is available, and it isassumed that the model m, with n number of cells, can be expressed by a smallnumber of j parameters denoted as p. This can be expressed bym = f (p) (3.7)where f : R j→Rn is a known smooth function that is continuously differentiable.In some cases both assumptions about the model are valid. The model can bemade of a smooth background and an anomalous body that can be parametrized.That ism = ms+ f (p) (3.8)where ms is some smooth background and f (p) describes an anomalous conductiveor resistive body. This leads to the following regularized problem to be solved.minms,pφ(ms,p) =12(F(ms+ f (p))−d)>WTd Wd(F(ms+ f (p))−d)+βR(ms)(3.9)Here, R(·) is a regularization term that enforces smoothness on the backgroundmodel and β is a regularization parameter. Additional restrictions, such as boundson p or ms can also be invoked. Equation 3.9 is a discrete optimization problem47for the smooth background model ms and the parameters p. In general, it is non-convex and therefore care must be taken to obtain feasible solutions.A hybrid approach is proposed where a block coordinate descent (Gill et al.,1981) is used to fix ms and minimize over p in the first, or parametric stage, andthen fix p and minimize over ms in the second, or voxel-based stage. Anotheroption is to run the voxel-based stage first, prior to the parametric stage, but ex-perience with the code so far suggests that the former procedure produces betterresults. The proposed hybrid scheme runs through each stage once, but this processcould be iterated. The stopping criteria are when the inversion: a) has convergedto a target data misfit of unity b) has reached a maximum number of user-definediterations, or c) can not find a Gauss-Newton step that lowers the normalized datamisfit by more than 0.1%. The program may also be terminated if the inversionis deemed to be over-fitting the data by placing obvious resistivity artifacts neartransmitter or receiver locations.This hybrid approach has the advantage of scale separation, in that the para-metric inversion of f (p) typically affects data locally, whereas optimizing over msaffects data globally, and can fit large-scale and smooth features. In the first para-metric stage, the method searches for one anomaly of interest, either conductive orresistive, in a background ms. This background can either be a uniform half-spaceor a heterogeneous resistivity distribution from a priori information or previous in-version work. Once again, the time-dependent quasi-static Maxwell’s equations aresolved with initial and boundary conditions as shown in Equations 3.1 through 3.3.The parametric hybrid approach finds a best-fitting skewed Gaussian ellipsoid,by means of a finite-volume discretization on local and global ocTree meshes(Haber and Schwarzbach, 2014). The inversion requires an initial guess, whichis composed of the quantities, rx, ry, rz, φx, φy, φz, x0, y0, z0, ρ0, and ρ1. The valuesr and φ represent an estimate for the radius and rotation angle of the ellipsoid foreach Cartesian direction, while x0, y0 and z0 represent the center coordinates of theanomaly, and ρ0 and ρ1 are the background and anomalous resistivities. Resistiv-ity values can be fixed by the user, or alternatively, these resistivities can be set asactive parameters in the inversion. The initial guesses for radii and rotation anglesare multiplied together to give a transform matrix T as shown in Equations 3.10to 3.14. Equation 3.15 then forms a symmetric positive definite matrix M, com-48posed of stretching and skewing parameters m1 through m6, and more informationregarding rotation matrices can be found in Modersitzki (2003). All unscaled pa-rameters p˜ are scaled by element-wise division denoted by with the vector s, toimprove the conditioning of the system as seen in Equations 3.16 and 3.17. Thevector s is composed of an appropriate length scale L, and a characteristic resis-tivity ρˆ . In total, p contains 11 scaled parameters that are used in the parametricinversion.S =1rx0 00 1ry 00 0 1rz (3.10)Rx =1 0 00 cos(φx) −sin(φx)0 sin(φx) cos(φx) (3.11)Ry = cos(φy) 0 sin(φy)0 1 0−sin(φy) 0 cos(φy) (3.12)Rz =cos(φz) −sin(φz) 0sin(φz) cos(φz) 00 0 1 (3.13)T = SRxRyRz (3.14)M =m1 m4 m5m4 m2 m6m5 m6 m3= TTT (3.15)49p˜ =m1m2m3m4m5m6x0y0z0log(ρ0)log(ρ1)s =L−2L−2L−2L−2L−2L−2LLLlog(ρˆ)log(ρˆ)(3.16)p = p˜ s. (3.17)For any position x,y,z in the spatial domain Ω setx =xyz x0 =x0y0z0 (3.18)and the level-set function τ is introduced in each mesh cellτ = c− (x−x0)TM(x−x0) (3.19)where c represents a positive constant and will be discussed in more detail later inthis chapter. τ , ρ0, and ρ1 are used to generate the resistivity distribution throughan analytic step functionρ(τ,ρ0,ρ1) = ρ0+12(1+ tanh(aτ))(ρ1−ρ0) (3.20)50wherelimτ→−∞ρ(τ,ρ0,ρ1) = ρ0limτ→+∞ρ(τ,ρ0,ρ1) = ρ1ρ(τ = 0,ρ0,ρ1) =12(ρ0+ρ1).A hyperbolic tangent is chosen for the analytic step function, but other choices arepossible (Tai and Chan, 2004). The transition zone between ρ0 and ρ1 occurs whenτ = 0, also known as the zero level set (Osher and Sethian, 1988), and its width iscontrolled by the parameter a. For numerical stability, a minimum of two meshcells should occur within the transition zone to ensure a suitable Gauss-Newtonstep can be found. The optimization of the inversion follows a conventional Gauss-Newton procedure for time-domain AEM data (Haber and Schwarzbach, 2014),and a line search is used to determine an appropriate model update step within aminimum and maximum value. Other techniques such as the trust region approachintroduced by Sturler and Kilmer (2011) can be applied if an acceptable step lengthis not found with a conventional line search, although this issue has not been en-countered in this research. In the case where m = f (p) is sufficient to model theregion of interest, the parametric code can be used as a stand-alone algorithm wherems is fixed. Otherwise, the hybrid technique is achieved by setting the parametricinversion model as the initial and reference model when optimizing over ms.3.4 Parametric Sensitivity and Initial ParameterSelectionGauss-Newton optimization for minimizing Equation 3.9 requires a sensitivity ma-trix Ji, j, or ∂di∂ p j , where di is the ith data point and p j is the jth inversion parameter.For p1 through p9, the chain rule is used to calculate Ji, j as shown in Equation 3.21Ji, j = ∑α,β∂di∂ρα∂ρα∂ρβ∂ρβ∂τ∂τ∂ p j(3.21)51where ρα is the cell center resistivity for each local mesh cell, α is the local meshcell index, ρβ is the cell center resistivity for each global mesh cell, and β is theglobal mesh cell index. This domain separation into local and global meshes is usedfor maximum computational efficiency as described in Haber and Schwarzbach(2014). For p10 and p11, Ji, j is shown in Equation 3.22 where the derivative withrespect to τ is replaced with a derivative with respect to ρ0 or ρ1Ji, j = ∑α,β∂di∂ρα∂ρα∂ρβ∂ρβ∂ρ0,1∂ρ0,1∂ p j(3.22)The voxel-based framework calculates the derivatives ∂di∂ρα∂ρα∂ρβ, therefore, the para-metric algorithm requires the evaluation of ∂ρβ∂τ∂τ∂ p j or∂ρβ∂ρ0,1∂ρ0,1∂ p j . These derivativesare expressed below in Equation 3.23.52∂ρβ∂m1= −0.5a(ρ1−ρ0)(1− tanh2(aτ))(x− x0)2∂ρβ∂m2= −0.5a(ρ1−ρ0)(1− tanh2(aτ))(y− y0)2∂ρβ∂m3= −0.5a(ρ1−ρ0)(1− tanh2(aτ))(z− z0)2∂ρβ∂m4= −a(ρ1−ρ0)(1− tanh2(aτ))(x− x0)(y− y0)∂ρβ∂m5= −a(ρ1−ρ0)(1− tanh2(aτ))(x− x0)(z− z0)∂ρβ∂m6= −a(ρ1−ρ0)(1− tanh2(aτ))(y− y0)(z− z0)∂ρβ∂x0= a(ρ1−ρ0)(1− tanh2(aτ))(M(x− x0))∂ρβ∂y0= a(ρ1−ρ0)(1− tanh2(aτ))(M(y− y0))∂ρβ∂ z0= a(ρ1−ρ0)(1− tanh2(aτ))(M(z− z0))∂ρβ∂ρ0= 0.5(1− tanh(aτ))∂ρβ∂ρ1= 0.5(1+ tanh(aτ))(3.23)To calculate Ji, j, the user defined inversion parameters L, ρˆ,a,c,, from Equa-tions 3.16, 3.19, and 3.20 need to be chosen. Experience suggests to choose L suchthat it represents a typical length scale of the problem, meaning it should be on theorder of a few cell lengths, and in this thesis, L = 100. Select ρˆ such that it repre-sents a moderately low resistivity value in the inversion. With resistivities varyingby many orders of magnitude this can be difficult to select, but 10 Ωm is chosenthroughout this work. Based on initial tests, the exact choice of L or ρˆ does notsubstantially change the inversion results and should not be a critical choice.The parameters a and c collectively change the width of the transition zonebetween ρ1 and ρ0, as depicted in Figure 3.1 with initial parameters [rx,ry,rz] =53[25,50,25], [φx,φy,φz] = [0,0,0], [x0,y0,z0] = [0,0,0] and [ρ0,ρ1] = [0,1]. Select-ing the parameter a to be larger results in a smaller transition zone. For the ex-amples shown in this thesis, a = 10, such that a minimum of two mesh cells arepresent within the transition zone. It is suggested to choose a between 5 and 15for best results. The inversion is less sensitive to the parameter c, but it is recom-mended to be set to 1 to keep initial radius guesses rx, ry and rz consistent withthe mapped radius values of the ensuing ellipsoid after the transformations fromEquations 3.11 to 3.20. Figures 3.1a, 3.1b and 3.1c demonstrate a sharp step-offfunction followed by the resulting ellipsoids with a = 10 and the parameter c set to1 and 5 respectively. Figures 3.1d, 3.1e and 3.1f show a gradual step-off followedby the corresponding ellipsoids with a = 0.5 and the parameter c set to 1 and 5 re-spectively. Figure 3.1b displays the ellipsoid constructed with the suggested choiceof inversion parameters of a= 10 and c= 1.54a) b) c)d) e) f)y (m)x (m) x (m)y (m)x (m) x (m)y (m)y (m)Figure 3.1: a) Analytic step-function with a = 10. b) Gaussian ellipsoid ex-ample with a = 10, c = 1. c) Gaussian ellipsoid example with a = 10, c =5. d) Analytic step-function with a = 0.5. e) Gaussian ellipsoid examplewith a = 0.5, c = 1. f) Gaussian ellipsoid example with a = 0.5, c = 5.553.5 Synthetic Results3.5.1 Thin plate responsesPrior to inversion, it is worthwhile to look at ∂bz∂ t responses from a concentric looptime-domain AEM system over a thin, vertical or dipping conductive plate. Thediscretized triangular waveform used to calculate these curves is depicted in Fig-ure 3.2 with measured time channels marked in black for reference. This waveformhas a longer on-time pulse coupled with a slower shut-off decay compared to theNEWTEM-I waveform from Chapter 2, which is beneficial for detecting deep con-ductive targets (Allard, 2007).Time (s)Normalized current (A)Time (s)Figure 3.2: Triangular discretized waveform in red with measured time chan-nels marked as black crosses.Figure 3.3 displays three scenarios with a conductive 3Ωm thin plate of dimen-sions 20 x 300 x 300 m in x,y,z buried 20 m below a flat topographic surface in aresistive half-space of 3000 Ωm. The transmitter and receiver are located 37.5 mabove the earth, and the data consists of 19 time channels recorded between 10 µs- 7000 µs. Figure 3.3a demonstrates the classic symmetric double-peak responseover a vertical dipping plate. Figure 3.3b displays how the symmetry is brokenwhen the plate is tilted at a 45 degree angle, and Figure 3.3c shows a single broadhigh over a flat-lying thin conductive sheet.56y(m)x(m)a) b) c)x (m) x (m) x (m)δBz/dt (V / Am2 )δBz/dt (V / Am2 )δBz/dt (V / Am2 )Figure 3.3: Noise-free ∂bz∂ t responses from a coincident loop AEM systemflying 37.5 m over a thin 3 Ωm conductive plate in a 3000 Ωm back-ground. a) Vertical plate buried 20 m below the surface. b) 45 degreedipping plate buried 20 m below the surface. c) Flat plate buried 20 mbelow the surface.3.5.2 Synthetic dipping plate inversion resultsTo test the parametric inversion, a synthetic model of a buried plate-like dippingconductor is constructed in a resistive background. For simplicity, flat topographyis included in this example. Figures 3.4a, 3.4b, and 3.4c portray sections throughthe synthetic model in each Cartesian direction. A time-domain AEM survey issimulated at a height of 37.5 m above ground with 49 transmitter locations overthe dipping plate as shown by stars in Figure 3.4a. The data is composed of ∂bz∂ tdata, contaminated with 5% Gaussian noise, calculated at 19 time channels usingthe waveform and time gates shown in Figure 3.2. The model is discretized on anocTree mesh with core cells of 20 × 20 × 20 m for a total of 96,272 cells in theinversion mesh, and roughly 10,000 cells in the forward meshes. The resistivityof the plate is 3 Ωm in a uniform 3000 Ωm background. The dimensions of thedipping plate are 320× 60× 640 m in x,y,z respectively with a near vertical dip of80 degrees in the positive y direction (east). Figures 3.4d, 3.4e, and 3.4f show theresults from a voxel-based inversion alone (Haber and Schwarzbach, 2014) with acompact Ekblom regularization (Ekblom, 1973) on the smoothness component ofthe model norm, which demonstrates how a voxel algorithm has trouble imaging57this type of target. The resulting anomaly gets smeared, the resistivity magnitudeis much too high, and the inversion is unable to capture the true dip.For the parametric inversion, the initial guess consists of a 50 m radius spherelocated at [x,y,z] = [-120, 0, -150] in a uniform background. This is the same initialguess as the voxel-based inversion, although a similar voxel result ensues if a half-space is used as the starting guess. In the first example, the true anomalous andbackground resistivity is presumed known. In the second trial, incorrect values areassigned and the inversion calculates the resistivity parameters. Having a sphereas a starting guess provides minimal information with regards to the true size anddip of the anomaly, and it is a reasonable starting guess without having any a prioriknowledge. Note, in the case where a portion of the recovered ellipsoid residesabove ground, this region is set to a typical resistivity value of air.58-500 -400 -300 -200 -100 0 100 200 300 400 500-500-400-300-200-10001002003004005000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-500-400-300-200-100010020030040050033.053.13.153.23.253.33.353.43.453.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000--700-600-500-400-300-200-100033.053.13.153.23.253.33.353.4.453.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000--700-600-500-400-300-200-100033.053.13.153.23.253.33.353.4.453.5-500 -400 -300 -200 -100 0 100 200 300 400 500-500-400-300-200-10001002003004005000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-500-400-300-200-10001002003004005000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-1000-900-800-700-600-500-400-300-200-10000.511.522.533.5-500 -400 -300 -200 -100 0 100 200 300 400 500-500-400-300-200-10001002003004005000.511.522.533.5a) b) c)d) e) f)g) h) i)x (m) x (m) y (m)x (m) x (m) y (m)x (m) x (m) y (m)y (m)z (m)z (m)y (m)z (m)z (m)y (m)z (m)z (m)z = -250 m y = 0 m x = -50 mz = -250 m y = 0 m x = -50 mz = -250 m y = 0 m x = -50 mlog10(Ωm) log10(Ωm) log10(Ωm)log10(Ωm) log10(Ωm) log10(Ωm)log10(Ωm) log10(Ωm) log10(Ωm)j) k) l)x (m) x (m) y (m)y (m)z (m)z (m)z = -250 m y = 0 m x = -50 mlog10(Ωm) log10(Ωm) log10(Ωm)Figure 3.4: Plan view depth slices and cross-sections through the true model,recovered voxel-based inversion and recovered parametric models. a)True model at z = -250 m. b) True model along y = 0 m. c) True modelalong x = -50 m. d) Voxel model at z = -250 m. e) Voxel model alongy = 0 m. f) Voxel model along x = -50 m. g) Parametric model at z =-250 m for fixed ρ . h) Parametric model along y = 0 m for fixed ρ . i)Parametric model along x = -50 m for fixed ρ . j) Parametric model at z= -250 m for variable ρ . k) Parametric model along y = 0 m for variableρ . l) Parametric model along x = -50 m for variable ρ .59When the anomalous and background resistivities are fixed to the true values,the optimization over p finishes in 32 Gauss-Newton iterations. For the parametricstage, only a small number of poorly correlated parameters are sought, which initself acts as a form of regularization. As such, an extra regularization term is notrequired and β is set to zero. Based on testing, this approach has been successful,but different regularization schemes in parametric and level-set inversions can befound in Aghasi et al. (2011); Dorn and Lesselier (2006); van den Doel and Ascher(2006); Zheglova and Farquharson (2016). At this point, the voxel-based stagesolves for ms, but in this example there are no regional background features toresolve, and the sharp resistivity contrast between the homogeneous backgroundand the thin plate favors a parametric inversion. Consequently, the voxel-basedstage does not decrease the data misfit without adding spurious inversion artifacts.Although the parametric result does not reach the target data misfit, the voxel-basedstage is considered to over-fit the data and the model after the first parametric stageis chosen as the final answer.Although the concept of a target misfit is valuable for a stopping criteria, itis acknowledged that the Gaussian ellipsoid parameterization may not perfectlyrepresent the true model or anomaly of interest, and the inversion may not reachthe target misfit. As such, if the inversion terminates due to a stalled data misfitbut does not reach the target level of unity, this suggests that the inversion hasprogressed as far as possible based on the assigned uncertainty values. In thiscase, the resulting model should have a significant reduction in data misfit to beconsidered an acceptable final result.Sections through the recovered synthetic model are shown in Figures 3.4g, 3.4h,and 3.4i. These images show that the parametric algorithm is able to accurately re-cover the size, shape and dip of the anomaly compared to the true answer. The dipof the parametric model is roughly 74 degrees to the east, and this closely matchesthe true dip of 80 degrees to the east. One aspect not perfectly detected is the po-sition of the plate bottom, which sits around 600 m in the parametric recovery and800 m in the true model. This discrepancy can most likely be attributed to reducedsensitivity to model cells at depths greater than 600 m. Similar parametric resultsare produced when the initial guess is centered at z = -300 and z = -450, so therecovery is not highly sensitive to the exact depth of the initial guess.60The next example allows the resistivity to be variable, and incorrect anomalousand background values of 1 Ωm and 1000 Ωm are assigned. The choice of startingguess is open to the user, but naturally, the further the starting guesses are fromthe true resistivities, the more difficulties the program will have converging to thecorrect solution. With the chosen starting guesses, the parametric stage concludesafter 29 Gauss-Newton iterations. As before, the voxel-based stage adds little im-provement, and the model after the parametric stage is chosen as the final result.Cross-sections through the recovered variable resistivity inversion are displayed inFigures 3.4j, 3.4k, and 3.4l. The inversion defines the shape of the target withaccurate precision with an estimated dip of 73 degrees to the east. The depth extentof the anomaly is slightly closer to the true model compared to the case when thetrue resistivities are assigned. Recovered anomalous and background resistivitiesare 3.97 Ωm and 3173 Ωm respectively.A plan map of observed and predicted data at a mid-range time channel, 1110µs, is displayed in Figure 3.5a and 3.5b. The plot is for the fixed resistivity in-version where the true values are assigned, although predicted data are similar inthe variable resistivity case. An observed and predicted sounding from a selectedlocation, marked with a cross in Figure 3.5a, is plotted in Figure 3.5c. Collectively,these images demonstrate the high level of agreement between the observed andpredicted data. The initial and final data misfits are 70.4 and 3.8 for the fixed resis-tivity case and 2328.4 and 2.8 for the variable resistivity scenario. For both trials,the misfit at each Gauss-Newton iteration is summarized in Figure 3.5d. The misfitsummary shows that both parametric inversions make excellent strides in reducingthe data misfit but fall short of the target level.61-200 0 200-300-200-1000100200300observed data×10-120.511.522.533.5-200 0 200-300-200-1000100200300predicted data×10-120.511.522.533.50 5 10 15 20 25 30100101102103Fixed SigmaVariable SigmaTime (s)10-5 10-4 10-3Response10-1410-1210-10x,y,z = [ 100.00, 0.00, 37.50]observed datapredicted dataa) b)c)(V / Am2) (V / Am2)x xY (m)Y (m)X (m) X (m)observed ata predicted dataTime (s)ResponseGauss-Newton iteration(V / Am2)data misfitd)Figure 3.5: Synthetic dipping plate observed and predicted ∂bz∂ t data. a) Ob-served data at 1110 µs for fixed ρ , with a selected sounding markedwith a cross. b) Predicted data at 1110 µs for fixed ρ , with a selectedsounding marked with a cross. c) Fixed ρ observed and predicted data atselected sounding location. d) Fixed and variable resistivity data misfitprogression.62A possible explanation for the lower final misfit when using a variable resis-tivity instead of true values is that a skewed Gaussian ellipsoid can not perfectlyrecover the staircase nature of a discretized dipping plate. It is also worth remem-bering that within the transition zone between ρ0 and ρ1, the anomaly does notcontain the true resistivity value, but instead a weighted average of ρ0 and ρ1.Therefore, if the anomalous resistivity quantity is allowed to vary, it is possiblethat the inversion can find a shape and resistivity combination that fits the observeddata better than correct values of ρ0 and ρ1. As an additional check to the code, theinversion reaches the same solution and data misfit as shown in the variable resis-tivity example when the initial guess is composed of spheres with fixed resistivitiesthat match the previously achieved values of ρ0 = 3.97Ωm and ρ1 = 3173Ωm. Thisensures that the fixed and variable version of the code are acting consistently.Based on testing not presented, the parametric inversion achieves similar re-sults when dealing with smaller resistivity contrasts, one or two orders of mag-nitude difference, between the background and target resistivity. In general, theencouraging results from the synthetic dipping plate parametric inversion suggestthat similar success can be found with a more complicated field example.3.6 Caber Case Study Results3.6.1 Caber geologyThe parametric hybrid approach is now applied to data from the Caber volcanogenicmassive sulfide (VMS) deposit of western Quebec, Canada, as shown in Figure 3.6a.Within the Superior Province, the copper and zinc rich Caber deposit is part of theMatagami camp of the Abitibi greenstone belt (Carr et al., 2008). Geologically, theprominent McIvor fault separates Caber and accompanying gabbros, rhyolites andbasalts from a granodiorite unit to the north-east (Adair, 2011). The local geologyin plan view is depicted in Figure 3.6a and a simplified cross-section is displayedin Figure 3.6b (Prikhodko et al., 2010). The cross-section shows there is a thinconductive overburden layer above the ore-body that thickens to the north-east.The deposit itself is cigar-like in shape with a steep dip of 75 to 85 degrees to thesouth-west, and a strike direction of 125 degrees (Adair, 2011). The cross-section63also shows the near-vertical nature of a narrow shear zone proximal to the deposit,and nearby steeply dipping rhyolite and gabbro units.The thin, buried nature of the Caber deposit, coupled with its position belowconductive overburden, makes it a challenging target to detect with AEM tech-niques. Fortunately, the elevated conductivity of the deposit compared to surround-ing rock units produces an anomalous electromagnetic response that is measurablefrom the air (Prikhodko et al., 2010). A small shear zone next to the deposit, whichmay be more conductive relative to the background but more resistive than the de-posit, may contribute to the conductive response. For the purpose of this research,the response from the deposit and any contribution from the shear zone will beconsidered the target of interest.McIvor Fault0 1000metersQuebecCaber DepositDiabaseCaberRhyoliteGabbroGranodioriteSW NEShear zone100 m100 m200 m300 ma) b)3002001SulfidesFigure 3.6: a) Caber deposit location and geology, modified from Prikhodkoet al. (2010) and Adair (2011). b) Simplified deposit cross-section withdrilling traces.3.6.2 Caber time-domain AEMThe Caber AEM field data are inverted with the parametric hybrid approach. TheAEM data at Caber consist of eight lines of Versatile Time-Domain Electromag-netic (VTEM-35) data, collected in 2012 with a 35 m diameter transmitter loop64and a peak dipole moment of 1,300,000 Am2. See Allard (2007) for further infor-mation regarding the VTEM system. Due to relatively flat terrain, the drape of thetransmitter above the ground varied between 34 m to 46 m, with a mean value of 40m. The flight lines are spaced 50 m apart with a heading of 225 degrees. The dis-cretized modified-triangular VTEM waveform is displayed in Figure 3.7 with 19time channels ranging from 167 µs to 2021 µs marked with black x’s for reference.Once again, a slower transmitter current shut-off compared to NEWTEM-I is de-signed to achieve a deeper penetration of the signal, but it is modified from a puretriangular wave to achieve a slightly steeper turn-off slope to detect near-surfaceweak conductors (Allard, 2007).Time (s)Normalized current (A)Time (s)Figure 3.7: VTEM-35 (2012) discretized waveform in red with measuredtime channels marked as black crosses.3.6.3 Caber parametric inversionsFor the parametric inversion stage, a subsection of the total AEM survey over theCaber deposit is used, consisting of 102 transmitters of ∂bz∂ t data with 11 time chan-nels ranging from 505 µs to 2021 µs. The model is discretized on an ocTree meshwith core cells of 20 × 20 × 20 m for a total number of 63,876 cells in the in-version mesh. The initial guess is a 50 m radius sphere, buried 150 m below thesurface, with a resistivity of 0.2 Ωm positioned in the center of the observed ∂bz∂ tanomaly with a background of 1000 Ωm. In this field example, the true anomalous65and background quantities are unknown, and ρ0 and ρ1 are variables solved for inthe inversion. Data uncertainties of 5% plus a time-channel dependent noise floorare applied. The noise floor is set to one order of magnitude lower than a 1000Ωm half-space response. This varying noise floor is selected to weight each timechannel as equally as possible in the inversion. The optimization over p concludesafter 30 Gauss-Newton iterations and sections through the initial guess and result-ing model in each Cartesian direction are shown in Figures 3.8a through 3.8f. Therecovered model has a steep 80 degree dip to the south-west, which agrees with theknown dip of the deposit. The final resistivities of the anomaly and background are0.084 Ωm and 2118 Ωm respectively. As in the synthetic case, similar results areproduced if the initial guess is moved to various locations and depths near the cen-ter of the observed ∂bz∂ t anomaly. The achieved anomaly value of 0.084 Ωm is wellwithin the range of resistivities for massive sulfides (Palacky, 1988), and is similarto the value of 0.14Ωm over a 30 m thickness obtained by Maxwell plate modeling(EMIT, 2005) of the 2012 AEM data over Caber from previous work (Prikhodkoet al., 2012). Overall, the similarities between the two results are encouraging,and an exact match to a 30 m plate can not be expected with a Gaussian ellipsoidon a 20 m mesh discretization. Even with 15 m mesh cells, the results would notbe expected to be identical since the methods use different approximations andparameterizations.3.6.4 Caber hybrid parametric inversionA plan map of observed and predicted parametric data for a mid-range time chan-nel, 1010 µs, is displayed in Figure 3.9a and 3.9b. Overall there is a close agree-ment between the observed and predicted data, as they both exhibit a similar asym-metric double peaked response, indicative of a single dipping plate anomaly. Somediscrepancies between the data sets exist, such as the break in the observed datathat separates the southern lobe of the double peak that is not present in the pre-dicted data. The inability to resolve fine levels of detail, such as this break, is alimitation of the first parametric stage of the hybrid approach, but the primary pur-pose of finding a best-fitting skewed Gaussian ellipsoid that matches the overalltrend of the data is achieved. Observed and predicted data at a selected sounding66×1057.098 7.099 7.1 7.101 7.102 7.103 7.104 7.105 7.106×1065.51315.51325.51335.51345.51355.51365.51375.51385.5139-1-0.500.511.522.53×1065.5132 5.5134 5.5136 5.5138-10001002003002.42.62.833.23.43.6×1065.5132 5.5134 5.5136 5.5138-1000100200300-10123×1057.098 7.1 7.102 7.104 7.106-1000100200300-10123×1065.5132 5.5134 5.5136 5.5138-1000100200300-10123×1057.098 7.1 7.102 7.104 7.106-1000100200300-10123y (m)y = 5513510 mx = 710146 m×1057.098 7.099 7.1 7.101 7.102 7.103 7.104 7.105 7.106×1065.51315.51325.51335.51345.51355.51365.51375.51385.5139-1-0.500.511.522.53a) b) c)x (m)x (m)y (m)y (m)z (m)z (m)log10(Ωm) log10(Ωm)log10(Ωm)z = 142.5 my = 5513510 md) e) f)x (m)x (m)y (m)z (m)z (m)z = 42.5 mx = 710146 mlog10(Ωm) log10(Ωm)log10(Ωm)×1057.098 7.1 7.102 7.104 7.106-10001002003002.42.62.833.23.43.6×1057.0987.099 7.1 7.1017.1027.1037.1047.1057.106×1065.51315.51325.51335.51345.51355.51365.51375.51385.51392.42.62.833.23.43.6y (m)y = 5513510 mx (m)x (m)y (m)z (m)z (m)z = 42.5 mx = 710146 mlog10(Ωm) log10(Ωm)g) h) i)log10(Ωm)Figure 3.8: Plan view depth slices and cross-sections through the initial guessand recovered Caber parametric and hybrid models. a) Initial guess at z= 142.5 m. b) Initial guess along y = 5513510 m. c) Initial guess alongx = 710146 m. d) Parametric model at z = 42.5 m. e) Parametric modelalong y = 5513510 m. f) Parametric model along x = 710146 m. g)Hybrid model at z = 42.5 m. h) Hybrid model along y = 5513510 m. i)Hybrid model along x = 710146 m.in the center of the northern lobe, marked with a cross in Figure 3.9b are shownin Figure 3.9c, and a strong agreement is evident. The initial and final data misfitsare 58.6 and 5.7 respectively, meaning the parametric approach reduces the datamisfit by an order of magnitude, and the data misfit progression is summarized in67Figure 3.9d.0 5 10 15 20 25 300102030405060Time (s) ×10-30.5 1 1.5 2×10-12468101214 observed datapredicted data×1057.1 7.1005 7.101 7.1015 7.102 7.1025 7.103 7.1035×1065.51335.51335.51345.51345.51355.51355.51365.51365.51375.5137×10-122345678×1057.1 7.1005 7.101 7.1015 7.102 7.1025 7.103 7.1035×1065.51335.51335.51345.51345.51355.51355.51365.51365.51375.5137×10-122345678a) b)c) d)(V / Am2) (V / Am2)x xy (m)y (m)x (m) x (m)responseGauss-Newton iteration(V / Am2)observed data Predicted datadata misfitObserved datatime (s)Figure 3.9: Caber observed and predicted ∂bz∂ t data. a) Observed data at 1010µs with a selected sounding marked with a cross. b) Predicted data at1010 µs with a selected sounding marked with a cross. c) Observed andpredicted data at selected sounding location. d) Data misfit progression.For the second stage of the hybrid approach, the parametric result is placed asthe initial and reference model in order to optimize over ms, allowing the resistiv-ity in every mesh cell to vary. Data from 19 time channels ranging from 167 µs to2021 µs are included, and responses below a threshold of 1e-13 VAm2 are discardedbecause of potential noise concerns. In total, 727 transmitter locations from acrossthe entire survey area are inverted with the same error assignments as in the para-68metric stage. Data over the central anomaly are included to allow the resistivityand shape of the parametric result to adjust if needed.Over nine Gauss-Newton iterations, each composed of three inner iterationsfor a fixed β , ms is optimized and the inversion converges to a data misfit of 0.95.The regularization parameter β starts at a value of 1e-3 and decreases by a fac-tor of γ = 0.2 with each Gauss-Newton iteration. The starting β is much smallercompared to Chapter 2 because the parametric starting guess is a much more de-tailed initial model compared to a half-space, and the same cooling scheme asChapter 2 is kept. Figures 3.8g, 3.8h and 3.8i show sections through the final hy-brid parametric model. Once again the thin conductor is imaged, representing theCaber deposit. A different color bar is used for the hybrid result to show both theweakly conductive overburden and the strongly conductive target together. TheCaber anomaly steeply dips to the south-west with a dip of 80 degrees, and anoverburden unit over the conductor initially thickens to the north-east. The re-covered size, dip and presence of overburden corroborates geologic informationfrom Figure 3.6b, which adds confidence to this hybrid result. The shape and sharpboundaries of the dipping conductor are achieved through the parametric inversion,while the voxel-based inversion adds smooth features such as the overburden. Al-though the voxel-based stage can potentially alter the shape and resistivity of theCaber deposit, the anomaly only changes slightly from the parametric result.Observed and predicted ∂bz∂ t data from time channel 505 µs, are illustrated inFigure 3.10. The predicted data at this time channel closely resemble the observeddata, and clearly illustrate the response of the Caber deposit in addition to theconductive overburden in the north-east portion of the survey area. Holes in theobserved and predicted data represent areas of resistive terrain where ∂bz∂ t responsesdrop below 1e-13 VAm2. Data from much of the south-west survey area is not showndue to data below this noise threshold.To validate and ground truth the field inversion, Figure 3.11 displays a frontand side view of a 0.4 Ωm iso-surface from the parametric result in red, the hybridresult in hatched gray, the massive sulfide deposit outline from drilling (M. Allard,personal communication, 2014) in light gray, and the aforementioned plate modelsfor individual lines (Prikhodko et al., 2012) in yellow. The image demonstrateshow the depiction of the Caber anomaly is extremely similar after both the para-69×1057.1 7.105 7.11×1065.5135.51325.51345.51365.51385.5145.51425.51445.51465.5148-11.6-11.5-11.4-11.3-11.2-11.1-11-10.9-10.8-10.7Source : dataTD-pred Component : dBz Time : 0.000505×1057.1 7.105 7.11×1065.5135.51325.51345.51365.51385.5145.51425.51445.51465.5148-11.6-11.5-11.4-11.3-11.2-11.1-11-10.9-10.8-10.7Source : dataTD Component : dBz Time : 0.000505a) b)x (m) x (m)y (m)y (m)Observed data Predicted datalog10(V / Am2) log10(V / Am2)Figure 3.10: Caber observed and predicted data at 505 µs from the paramet-ric hybrid inversion. a) Observed data. b) Predicted data.metric and hybrid stages, however, the hybrid result is treated as the final model.The hybrid iso-surface accurately portrays the general size and dip of the targetanomaly, although it stretches further to the south-east compared to the depositmodel. Interestingly, previous AEM plate modeling also places conductive platesextending off the deposit. This suggests that additional conductive material to thesouth-east of the recorded sulfide zone is needed to explain the observed response.This anomalous material outside the deposit could be due to an unknown extensionof the economic mineralization or the mapped shear zone shown in Figure 3.6b, orsome combination of the two. As shown in the synthetic example, there is reducedinversion sensitivity at depth and caution must be taken when interpreting deepfeatures. However, the parametric hybrid inversion provides a new interpretationto the cause of the conductive response from the 2012 time-domain AEM data.The Gaussian ellipsoid parameterization has been shown to be successful forthin targets, but it also has the flexibility to map out larger anomalies such as in-trusions or kimberlites. McMillan et al. (2015a) provides an example not shownin this thesis, where the parametric inversion uses time-domain AEM data to accu-rately image the DO-27 diamondiferous kimberlite deposit at Tli-Kwi-Cho (TKC)in the North-West Territories, Canada. Lessons from this work will be discussed in70a) b)x (m)710100710150710200 710250710350 (m)710350710016710200710250710100710150x (m)-135165100500-50-135z (m)165100500-50-135z (m)551329655134505513500551355055136185513400y (m)y (m)551350055135505513296551340055134505513618710016710350Maxwell platesCaber deposit outlineParametric iso-surfaceHybrid iso-surfaceFigure 3.11: Caber deposit outline with inversion results. 0.4 Ωm iso-surface(red) from parametric stage, and 0.4 Ωm iso-surface (hatched gray)from hybrid inversion overlaid on Caber massive sulfide deposit modelfrom drilling (light gray outline) and Maxwell plate anomalies (mul-tiple yellow sheets) from five central lines of AEM data (Prikhodkoet al., 2012). a) Looking north-east. b) Looking north-west.further detail in Chapter 5.3.7 ConclusionsA parametric hybrid inversion has been developed for AEM data to recover a tar-get represented by a skewed Gaussian ellipsoid in a smooth background. Thisapproach is valid for frequency and time-domain surveys, but it is tested in thischapter on time-domain synthetic and field AEM data. Here the target is a thinconductive plate-like body with a large resistivity contrast between itself and theresistive background. Both examples recover models that agree well with eitherthe true synthetic answer or geologic information from past drilling respectively.The approach can be used as an alternative to a purely voxel-based inversion wheresharp boundaries and large resistivity contrasts may not be imaged accurately. Theparametric hybrid models are produced using a basic initial guess of a sphere with-out a priori information, which adds robustness to the algorithm. The results show71that the choice of parameterization is well suited for a thin, dipping plate target, butthe approach is also applicable for other exploration environments such as mappinglarge intrusions or kimberlites.For simplistic cases such as a synthetic dipping plate, the parametric stage doesan admirable job in recovering the narrow target accurately. In a more complexexample, such as the Caber case study, a hybrid approach is required, which is ableto recover both the target of interest while including smooth features such as theconductive overburden. This task is not feasible with the parametric stage alone,and highlights the benefits of a hybrid approach.For this single-anomaly version of the code, care must be taken in the paramet-ric stage because in a scenario where two nearby conductive or resistive anomaliesexist, the response may be modeled erroneously as a single anomaly. Thus, any apriori knowledge of the area coupled with a careful analysis of the observed andpredicted data is needed to see whether one or multiple bodies are warranted. Inthis example, a priori information from plate modeling and drilling suggested thatthe Caber deposit could be depicted appropriately as one anomaly. Unfortunately,no physical property data existed from these bore holes to include as constraintsinto the inversion.72Chapter 4Multiple Anomaly ParametricInversion of AEM dataIn this chapter, I demonstrate an extension of the single anomaly parametric in-version method to include multiple anomalies with distinct resistivity values forfrequency and time-domain AEM data sets. The approach is designed to gener-alize the parametric inversion technique, making it more flexible for a variety ofgeologic settings and targets. The algorithm is tested with synthetic and field datafrom the West Plains orogenic gold deposit in Nunavut, Canada.4.1 IntroductionNatural resource deposits come in all shapes and sizes. Some can be straightfor-ward in their geometric nature as shown in the previous chapter, but this is true foronly a handful of cases. A single body parametric code can be applied to simplisticexamples, but for multiple targets, additional tools are needed. A logical progres-sion, shown in this chapter, is to incorporate multiple parametric bodies into theinversion scheme, thus allowing the anomalies to interact and arrange themselvesin more complex patterns. There are three primary research goals of this chapterand they are listed below.• To extend the parametric inversion framework to include multiple anomaliesand unique resistivity values for frequency and time-domain AEM data using73a skewed Gaussian ellipsoid parameterization.• To show, with synthetic and field examples, that the multiple anomaly para-metric inversion can accurately recover several thin dipping conductors, astypically found in orogenic gold settings.• To apply the hybrid parametric approach to improve data fits and to fill inmissing features not captured with a parametric inversion alone.This chapter begins with a discussion regarding the extension of the parametricmethod from single to multiple anomalies. Synthetic inversions are then shownfrom frequency and time-domain AEM data over a model with three thin dippingconductors in a resistive half-space. Following the synthetic case, the hybrid in-version is tested on frequency and time-domain AEM field data from an orogenicgold deposit in the Canadian Arctic. The field results are subsequently evaluatedin areas where geologic and mineralization knowledge from drilling exists.4.2 Inversion methodologyThe multiple anomaly parametric inversion method solves for an integer number(n) targets in a smooth or heterogeneous background. Each anomaly is once againrepresented by a skewed Gaussian ellipsoid; however, further discussion regardingthe use of a multiple sphere parameterization can be found in Appendix A. Equa-tion 4.1 shows the symmetric matrix Mn composed of the stretching and skewingparameters for the nth ellipsoid.Mn =m1(n) m4(n) m5(n)m4(n) m2(n) m6(n)m5(n) m6(n) m3(n) (4.1)Observation locations (x) are defined as described in the single anomaly ver-sion, with the central position of the nth ellipsoid (x0(n)) written asx0(n) =x0(n)y0(n)z0(n) . (4.2)74The parametric function τn is now defined for the nth anomaly asτn = c− (x−x0(n))TMn(x−x0(n)) (4.3)where c represents a positive constant, which is set to unity for all inversions. Thezero level set is taken (Osher and Sethian, 1988) with a hyperbolic tangent analyticHeaviside function to split anomalous regions from the background. In each meshcell, the separation of the nth ellipsoid, denoted as Hn, from the background isdefined asHn = 1−0.5(1+ tanh(aτn)) (4.4)where a controls the steepness of the Heaviside function and is once again set to10 for all inversions. Hn = 1 refers to a background resistivity value (ρ0) and Hn= 0 equates to an anomalous resistivity value (ρa). To assign specific anomalousresistivity values, the approach from Tai and Chan (2004) is applied. For twoanomalies in a homogeneous background the resistivity in each mesh cell is set toρ =ρ0(H1)(H2)+ρ1(1−H1)(H2)+ρ2(H1)(1−H2)+ρ3(1−H1)(1−H2).(4.5)Figure 4.1 shows an example of a cross-section through a parametric model withtwo anomalies, H1 and H2, with distinct resistivities, ρ1 and ρ2, an overlappingregion with resistivity ρ3 and a homogeneous background set to ρ0.For more than two anomalies, the method can be easily modified. One ap-proach is to use an arbitrary number of parametric bodies, but to keep four uniqueresistivity values as in Equation 4.5. In this scenario, Heaviside functions fromeach anomaly Hn can be multiplied together (element-wise) to form two globalHeaviside functions, Ψ1 and Ψ2, which are then substituted in for H1 and H2 inEquation 4.5. Here, all anomalies within Ψ1 are assigned a value of ρ1 and allwithin Ψ2 are set to ρ2. The overlap region of the two global Heaviside functionsis set to ρ3. The second approach is to assign a unique resistivity to each para-metric body, however, with each additional anomaly, an overlap region must be75x (m)z (m)log10(Ωm)ρ0ρ1 ρ2ρ3H1 H2Figure 4.1: Cross section through a two-anomaly parametric model with fourunique resistivity values (ρ0 - ρ3). H1 and H2 represent anomaly re-gions.defined between itself and every other anomaly. This quickly increases the numberof parameters needed, but may be necessary for certain situations. All examplesshown in this thesis are suitable for the former approach, and this method is thefocus moving forward.The multiple anomaly parametric algorithm recovers a user-defined number ofanomalies, and requires a starting guess for the location and resistivity of eachtarget like the single anomaly version. Similar to Chapter 3, the inversion op-timization follows a Gauss-Newton approach with a line search to determine anappropriate model update step, and the inversion parameters are scaled in the samemanner. The sensitivity matrix Ji, j is also computed as in Chapter 3, except that∂ρβ∂H∂H∂τ replaces∂ρβ∂τ for the derivatives with respect to location and skewing param-eters as shown in Equation 4.6.Ji, j = ∑α,β∂di∂ρα∂ρα∂ρβ∂ρβ∂H∂H∂τ∂τ∂ p j(4.6)764.3 Synthetic ResultsTo test the multiple body parametric approach, a frequency and time-domain AEMsynthetic modeling study is performed, simulating an Arctic gold exploration set-ting. The motivation for the synthetic study comes from field data over the WestPlains area of the Committee Bay greenstone belt in Nunavut, Canada. The loca-tion of this auriferous region is shown in Figure 4.2. This greenstone belt containsgold mineralization within ultramafic units such as komatiites and banded iron for-mations, which exist as a series of mainly linear, thin, dipping conductors in aresistive background (Kerswill, 1996). Hence, the synthetic model will containmultiple thin conductors, representing these ultramafic targets within a resistivehost.Figure 4.2: The West Plains area resides within the Archean Committee Baygreenstone belt in Eastern Nunavut, Canada774.3.1 Synthetic modelThe true synthetic model is composed of three dipping conductors of 0.2 Ωm, 2.0Ωm and 2.0 Ωm respectively in a 1000 Ωm background with flat topography asshown in Figure 4.3a. The left panel displays a plan view at an elevation slice of250 m (60 m below surface), and the solid white lines represent two East-Westcross-sections shown in the right panel from top to bottom at y = 7333080 m andy = 7332780 m respectively. This same plan view and cross-section layout areused for all images in the synthetic study. Two of the conductors have a dip of65 degrees with respect to the horizontal, one to the East and one to the West,while the third conductor has a vertical dip. The white dashed lines in plan viewoutline the starting guess locations for the parametric inversions. The syntheticdata locations are taken from two overlapping AEM field data sets collected overthe West Plains area: a 2005 frequency-domain RESOLVE (Viezzoli et al., 2009)and a 2003 time-domain VTEM (Allard, 2007) survey. The RESOLVE survey has25 lines collected with a 310 degree heading with a variable line spacing of 60 mand 120 m. The drape of the transmitter varies between 22 m and 39 m with amean terrain clearance of 30 m. The VTEM survey has 18 lines, and is also flownin a 310 degree direction with a variable line spacing of 60 m and 120 m. Thetransmitter height fluctuates between 22 m and 36 m with a mean height of 26 m.This average instrument drape is extremely low for a time-domain system, and isonly possible due to essentially flat topography.78y(m)x(m)y(m)x(m)x(m)z(m)a)c)ΩmΩmz(m)x(m)x(m)z(m)z(m)x(m)True modelTime domainy(m)x(m)x(m)z(m)b)Ωmz(m)x(m) Frequency domainFigure 4.3: Frequency and time-domain AEM synthetic voxel inversions inplan view at z = 250 m (left panel) and cross-section across solid whitelines at y = 7333080 m and y = 7332780 m (right panel). a) True syn-thetic model. b) Frequency-domain voxel inversion. c) Time-domainvoxel inversion.794.3.2 Synthetic dataForward modeled data are contaminated with 3% Gaussian noise, and 3D inver-sions are completed with a conventional voxel-based inversion algorithm (Haberand Ascher, 2001; Oldenburg et al., 2013), and then with the multiple body para-metric approach. In the frequency domain, the data consists of 793 source locationswith 10 logarithmically spaced frequencies between 315 Hz and 6000 Hz of realand imaginary Hz data. In the time domain, 522 source locations are used with 20logarithmically spaced time channels between 150 µs and 3180 µs of ∂bz∂ t data. Thesynthetic data uses a transmitter waveform similar to that seen in Figure 3.7 fromChapter 3, but with a 7.5 ms on-time pulse as shown in Figure4.4 to match the 2003VTEM system. The inversion mesh consists of 289,979 cells with a finest dimen-sion size of 20 x 20 x 20 m in x,y,z respectively, and the forward meshes containroughly 10,000 cells per mesh. Inversion uncertainties are assigned as 3% of theobserved data plus a noise floor that is one order of magnitude below the responsefrom a 1000 Ωm half-space, to weight each time channel as equally as possible.In order to reduce the size of the system to speed up numerical computations, thedata are down-sampled using a total horizontal gradient (THG) method of sam-pling. This technique computes the horizontal gradient amplitude for a particularfrequency or time and uses this information to discard data in regions with slowlyvarying responses while preserving measurements where the response is rapidlychanging. More details regarding THG sampling are described in Appendix B. Tokeep the locations consistent between the field and synthetic data, the THG sam-pling method is based on the field data observations at 385 Hz and 150 µs, althoughAppendix B shows that synthetic results improve when THG sampling is appliedbased on the synthetic data.The time channels selected are from the 2003 VTEM system, and frequenciesare chosen to give a similar depth of penetration compared to these time channelsas per the skin depth and diffusion distance shown in Equations 1.7 and 1.8. Thismeans that both data sets should have similar sensitivities to the target conductors.The synthetic study incorporates only 10 frequencies compared to 20 time channelsto honor the reality that currently fewer frequencies are available on airborne sys-tems such as RESOLVE compared to the number of time channels on time-domain80Time (s)Normalized current (A)Time (s)Figure 4.4: VTEM (2003) discretized waveform in red with measured timechannels marked as black crosses.airborne platforms such as VTEM, NewTEM, AeroTEM, HeliTEM MegaTEM andSkyTEM (Allard, 2007; Boyko et al., 2003; Eaton et al., 2013; Mule` et al., 2012;Smith et al., 2003; Sørensen and Auken, 2004).4.3.3 Synthetic voxel inversionsConventional voxel inversion results for the frequency and time-domain syntheticstudy are shown in Figure 4.3b and 4.3c respectively. The true background ispresumed to be known for this study, and a 1000 Ωm half-space reference modelis used with an Ekblom norm to promote thin anomalies. In Figures 4.3b and 4.3c,voxel inversions give a relatively accurate spatial recovery of the conductors in planview but struggle to clearly identify the true dip and depth extent of the anomaliesin cross-section.4.3.4 Synthetic parametric inversionsFigure 4.5 shows the parametric inversion results using the same data and uncer-tainty values as the voxel trials, with the true model displayed in Figure 4.5a, thefrequency inversion in Figure 4.5b and the time inversion in Figure 4.5c. Onceagain the true background is incorporated and the initial guess consists of three81y(m)x(m)y(m)x(m)x(m)z(m)a)c)ΩmΩmz(m)x(m)x(m)z(m)z(m)x(m)True ModelTime domainy(m)x(m)x(m)z(m)b)Ωmz(m)x(m) Frequency domainFigure 4.5: Synthetic parametric inversions in plan view at z = 250 m (leftpanel) and cross-section across solid white lines at y = 7333080 mand y = 7332780 m (right panel). Dashed white lines represent start-ing guess locations for parametric anomalies. a) True synthetic model.b) Frequency-domain parametric inversion. c) Time-domain parametricinversion.200 m radius spheres of 1 Ωm centered over the anomalies of interest at a depthof z = 150 m. The starting guess locations are marked by dashed white lines inFigure 4.5a. Both the frequency and time-domain parametric results successfullyspatially reconstruct the anomalies in plan view and give a much more accurateapproximation of the respective dips in cross-section compared to the voxel in-versions. The recovered target resistivities are 0.17 Ωm and 2.62 Ωm for the fre-82Table 4.1: Synthetic data inversion summary.Inversion model ρ1 ρ2 Rc RAll # GN Final φdTrue model 0.20 Ωm 2.00 Ωm 0.0 0.0 - -Voxel (freq) Variable Variable 2.73 0.39 96 4.2Voxel (time) Variable Variable 2.58 0.75 45 17.2Parametric (freq) 0.17 Ωm 2.62 Ωm 2.18 0.26 38 19.6Parametric (time) 0.15 Ωm 1.99 Ωm 3.33 0.31 28 17.6quency domain and 0.15 Ωm and 1.99 Ωm for the time domain. These valuesare extremely close to the true resistivities of 0.20 Ωm and 2.00 Ωm respectively.The results for all synthetic inversions are summarized in Table 4.1. ExaminingFigure 4.3 and Figure 4.5 suggests that a parametric inversion provides increasedinterpretation accuracy over voxel counterparts, and differences between frequencyand time-domain inversions are minimal.In an attempt to objectively evaluate the accuracy of each synthetic inversion,the quantitative metric from 2.1 is employed, and these values are displayed inTable 4.1. When only the conductive cells are taken into account, the frequency-domain parametric model outperforms the voxel result with a lower residual valueof 2.18 compared to 2.73, but surprisingly the time-domain parametric model has ahigher residual value of 3.33 compared to the voxel inversion at 2.58. This raises aconcern with assessing synthetic inversions, because quantitatively the voxel modelis more accurate, but the human eye can interpret and glean more true informationin terms of the size and dip of the various conductors from the parametric model. Incontrast, when all of the cells in the model are taken into account, both parametricmodels show a better quantitative result. These metrics demonstrate a potential is-sue when interpreting purely quantitative assessments of synthetic models, becausea model may be deemed more accurate but could be less useful to the interpreter.This suggests that an improved model for human interpretation is not necessarilyalways quantitatively more accurate.83Observed and predicted Hz data at 450 Hz and ∂bz∂ t data at 150 µs from thesynthetic parametric inversions are displayed in Figure 4.6 along with survey loca-tions as black dots. All predicted data visually match closely with observed mea-surements, and data misfit curves are presented in Figure 4.7. Both data sets showa comparable decrease in misfit although the final data misfits of 19.6 and 17.6 arewell above the target misfit of unity. This is most likely due to the small 3% as-signed error percentage, and future models can be run with a higher assigned noisepercentage to achieve a smaller data misfit, however similar results are expected.Note that the time-domain voxel inversion terminated at a similar data misfit levelof 17.2, whereas the frequency-domain voxel inversion converged to a value of4.2. This small misfit level for the frequency voxel model is most likely due tothe placement of high-frequency near-surface conductive features, which may beartifacts.Figure 4.8 shows a 5 Ωm iso-surface view of the parametric synthetic inver-sions, which gives a way of looking at only the conductive cells of interest. Fig-ure 4.8a displays the true model with the starting guess model shown in Figure 4.8b.Figure 4.8c and Figure 4.8d display the frequency and time-domain parametric in-versions. Collectively Figure 4.8 visually demonstrates how a Gaussian ellipsoidparameterization geometrically matches a dipping plate target extremely well. Ofcourse this example is designed for this purpose, but the motivation stems from acase study that will now be discussed.84a)c)b)d)x (m) x (m)x (m) x (m)y (m)y (m)y (m)y (m)e)y (m)x (m) f )y (m)x (m)log10(A/m) log10(A/m)log10(A/m) log10(A/m)log10(V/Am2) log10(V/Am2)Observed data Predicted dataFigure 4.6: Observed and predicted synthetic data from parametric inversionswith locations shown as black dots. a) Observed data - real Hz at 450Hz. b) Predicted parametric data - real Hz at 450 Hz. c) Observed data -imaginary Hz at 450 Hz. d) Predicted parametric data - imaginary Hz at450 Hz. e) Observed data - ∂bz∂ t at 150 µs. f) Predicted parametric data- ∂bz∂ t at 150 µs.85Gauss-Newton iterationData mistFigure 4.7: Synthetic parametric data misfit progression.a)c) d)x(m)y(m)z(m)x(m)y(m)z(m)x(m)y(m)z(m)310-130310-130310-130b)x(m)y(m)z(m)310-130Figure 4.8: Frequency and time-domain parametric inversion 5 Ωm iso-surfaces. a) True synthetic model. b) Starting guesses. c) Frequency-domain parametric inversion. d) Time-domain parametric inversion.864.4 West Plains Case Study ResultsBased on the success of the synthetic study, the multiple body parametric approachis now applied to field AEM data from the West Plains project.4.4.1 West Plains geologyWest Plains is an active orogenic gold exploration project, located in the South-West portion of the Archean Committee Bay greenstone belt, which itself is part ofthe larger Rae geologic domain, in Eastern Nunavut, Canada (Skulski et al., 2003).The project is owned by Auryn Resources Inc., and the economic gold mineral-ization is hosted in conductive units within a resistive background, a scenario wellsuited for AEM data. A simplified geology map of West Plains is shown in Fig-ure 4.9 with the outline of overlapping frequency and time-domain AEM surveysin black. The geology consists of a series of ultramafic komatiite units shown indark purple, named the Prince Albert Group, that cut through a sediment pack-age shown in light purple (Carson et al., 2004). The komatiites and sediments arebounded by granodiorites to the North-West and tonalites to the South-East. Thekomatiites contain the majority of the gold mineralization, and garnering a betterunderstanding of their 3D geometry represents the focal point of this study.Geologic outcrops are scarce in the West Plains region, and the simplified ge-ology map is compiled mainly from regional magnetic data analysis. The geologicinterpretation from Figure 4.9 does not include information from AEM data, whichmeans there is plenty of room for improvement based off 3D AEM inversion re-sults. From the geology and drilling, it is believed that three conductive units exist,comprised of a Western conductor, numbered 1 on Figure 4.9, and two near-parallelor possibly intersecting conductors in the East, numbered 2 and 3 on Figure 4.9.4.4.2 West Plains AEM dataAEM field data locations and inversion meshes are kept the same as the syntheticcase and the assigned uncertainties are applied in the same manner. Hz data fromfive coplanar frequencies (385, 1500, 6200, 25,000 and 115,000 Hz) are used in thefrequency domain, and the same 20 time channels of ∂bz∂ t data from the syntheticsection are inverted in the time domain.87KomatiiteSedimentsGranodioriteTonaliteSurvey outlinex (m)y (m) 12 3Figure 4.9: Simplified West Plains geology map with major lithology unitsdefined and location of overlapping frequency and time-domain AEMsurveys outlined in black. Conductive komatiite units of interest arenumbered in yellow.Selected data maps from West Plains are shown in Figure 4.10. Imaginary com-ponent Hz data from the lowest frequency, 385 Hz, is portrayed in Figure 4.10a andfrom the highest frequency, 115,000 Hz, in Figure 4.10b. The low-frequency gridclearly highlights the conductors of interest and has a relatively quiet signal else-where in the survey area. Conversely, the high-frequency grid has a response fromthe conductors, but also detects near-surface features such as clay layers associatedwith lakes and rivers. Subsets of the data could be used in an attempt to isolate theconductors, however all the data are input to show the robustness of the parametricalgorithm.Time-domain ∂bz∂ t data from the earliest time, 150 µs, are shown in Figure 4.10cand from the latest time, 3180 µs, in Figure 4.10d. The earliest time channel does88an excellent job in mapping out the target conductors without much influence fromlakes. However, at later times the data become noisier and measurements below aperceived noise threshold of 5e-12 VAm2are removed for quality control purposes.This threshold is chosen by looking at the amplitude at which point the electromag-netic response no longer consistently decays monotonically and oscillates due tonoise. As can be seen in Figure 4.10d, only data points near the strong conductorshave amplitudes above this threshold at late times. This perceived noise level is anorder of magnitude higher compared to data over the Caber deposit described inChapter 3. The noisiness of VTEM data at West Plains compared to Caber, em-phasizes the improvements in noise reduction of time-domain AEM data between2003 and 2012.4.4.3 West Plains voxel inversionsConventional voxel inversion results are shown for frequency and time-domainAEM data at West Plains in Figure 4.11a and 4.11b with plan view slices at aconstant elevation of z = 190 in the left panel and cross-sections from y = 7333080m and y = 7332780 m in the right panel. Inversion parameters are selected to matchthose from the synthetic example. Data misfit progression curves from the voxelinversions are shown in Figure 4.12a.The voxel inversions recover a linear Western conductor and two Eastern con-ductors although the distinction between the two conductive units in the East isnot consistently clear. In cross-section, the two inversions bare few similarities,as the conductors in the frequency-domain are compact near-surface objects whilethe time-domain result produces anomalies that are larger and extend to depth.The presence of near-surface features in the frequency domain and not in the timedomain is not entirely surprising, since the skin depth from observed frequenciesranges from 47 m to 810 m, whereas the diffusion distance from observed timechannels varies from 487 m to 2247 m, using a 1000 Ωm background.4.4.4 West Plains parametric inversionsTo complement the voxel models, parametric inversions are computed to imagemore clearly the three conductive komatiite units at West Plains. Appropriately,89a)c)b)d)x (m) x (m)x (m) x (m)y (m)y (m)y (m)y (m)log10(V/Am2) log10(V/Am2)log10(A/m) log10(A/m)Figure 4.10: Frequency and time-domain field data. a) Imaginary Hz at 385Hz. b) Imaginary Hz at 115000 Hz. c) ∂bz∂ t at 150 µs. d)∂bz∂ t at 3180µs.the parametric algorithms are set to recover three conductors, with the Westernconductor chosen to have its own resistivity value and the Eastern units sharing abest-fitting resistivity. This is based on a priori knowledge that the two Easternconductors are similar in nature. Figure 4.13a shows the frequency-domain fieldparametric inversion, with a plan view slice through the model at z= 190 m in theleft panel. The solid white lines represent two East-West cross-sections displayedin the right panel at y = 7333080 m and y = 7332780 m. The dashed white linesindicate the spatial outlines of the 1 Ωm starting guesses centrally located at adepth of z = 150 m. The recovered target resistivities in the frequency-domainare 0.15 Ωm and 7.28 Ωm with a background of 12,609 Ωm. Comparatively, the90y(m)x(m)x(m)z(m)a)Ωmz(m)x(m) Frequency domainy(m)x(m)x(m)z(m)b)Ωmz(m)x(m) Time domainFigure 4.11: Frequency and time-domain West Plains voxel inversions inplan view at z = 190 m (left) and cross-section across solid white linesat y = 7333080 m and y = 7332780 m (right). a) Frequency domain. b)Time domain.time-domain model, which is shown in plan view and cross section in Figure 4.13b,recovers resistivities of 0.28Ωm and 0.11Ωm with a background of 495Ωm. WestPlains parametric inversion statistics, including recovered resistivities, number ofGauss-Newton iterations required and final data misfit values are summarized inTable 4.2. Data misfit curves are displayed in Figure 4.12b, where the initial partof the curve represents the parametric portion while the latter half is the hybridinversion discussed in the next section.In plan view, the frequency and time-domain models have many spatial sim-ilarities, although the absolute recovered resistivities are much lower in the time-domain for the Eastern conductors (ρ2) and the background (ρ0). This suggeststhat there is some disagreement between the two AEM data sets, and the discrep-ancy in ρ0 will be addressed in greater detail in Chapter 5. In cross-section, theWestern conductor has a similar shape and steep dip to the West in both inversions;91Gauss-Newton iterationData mista) Gauss-Newton iterationData mistb)Figure 4.12: West Plains data misfit progression. a) Voxel inversions. b)Parametric and hybrid parametric inversions.Table 4.2: Field data inversion summary.Inversion model ρ1 ρ2 ρ0 # GN Final φdVoxel (freq) variable variable variable 49 1.50Voxel (time) variable variable variable 60 2.92Parametric (freq) 0.15 Ωm 7.28 Ωm 12609 Ωm 96 6.27Parametric (time) 0.28 Ωm 0.11 Ωm 495 Ωm 175 8.49Hybrid (freq) variable variable variable 47 1.87Hybrid (time) variable variable variable 40 3.51however, the Eastern conductors appear different. In the frequency domain, theEastern conductors are rounder compared to the time-domain results, with a steepdip to the East. In contrast, the time-domain model has a thin central conductor thatsteeply dips to the West, and a small near-surface Eastern conductor. Rectangularblocky anomalies in the North-East corner of the parametric models are outside theobserved data area, within padding cells, and should be ignored.Observed and predicted field data are shown in Figure 4.14. Real-componentHz data at 385 Hz are shown in Figure 4.14a along with parametric predicted data inFigure 4.14b. Observed and parametric predicted imaginary Hz data are displayedin Figure 4.14d and Figure 4.14e respectively. Time-domain observed ∂bz∂ t data at150 µs are shown in Figure 4.14g with parametric predicted data in Figure 4.14h.92y(m)x(m)x(m)z(m)a)Ωmz(m)x(m) Frequency domainTime domainy(m)x(m)x(m)z(m)b)Ωmz(m)x(m)Figure 4.13: West Plains parametric inversions in plan view at an elevationof z = 190 m (left) and cross-section across solid white lines at y =7333080 m and y = 7332780 m (right). Dashed white lines representstarting guess locations for parametric anomalies. a) Frequency do-main. b) Time domain.The parametric predicted data from Figure 4.14 closely resemble the observeddata, but further improvements can be made such as fitting two amplitude highs inthe real-component image of Figure 4.14a that are missing in the predicted data.Additional modifications to the parametric models to fit detailed features is thebenefit of the hybrid parametric method.4.4.5 West Plains hybrid parametric inversionsAs shown in Chapter 3, a hybrid parametric inversion incorporates the parametricresult as an initial and reference model for a voxel inversion. In the voxel stage,all cells are kept active, meaning that the resistivity value can change in each meshcell. This process allows information from the parametric stage to be passed to thevoxel stage, which then inserts complexities into the model not possible with Gaus-sian ellipsoids only. The hybrid method is applied to parametric results from the93a)d)b) c)e) f )x (m) x (m) x (m)x (m) x (m) x (m)y (m)y (m)y (m)y (m)y (m)y (m)g)y (m)x (m) h)y (m)x (m) i)y (m)x (m)log10(A/m) log10(A/m) log10(A/m)log10(A/m) log10(A/m) log10(A/m)log10(V/Am2) log10(V/Am2) log10(V/Am2)Observed data Predicted data - parametric Predicted data - hybridFigure 4.14: Observed and predicted field data from parametric and hybridinversions with locations shown as black dots. a) Observed data - realHz at 385 Hz. b) Predicted parametric data - real Hz at 385 Hz. c)Predicted hybrid parametric data - real Hz at 385 Hz. d) Observed data- imaginary Hz at 385 Hz. e) Predicted parametric data - imaginary Hzat 385 Hz. f) Predicted hybrid parametric data - imaginary Hz at 385Hz. g) Observed data - ∂bz∂ t at 150 µs. h) Predicted parametric data -∂bz∂ t at 150 µs. i) Predicted hybrid parametric data -∂bz∂ t at 150 µs.previous section and hybrid frequency and time-domain AEM models from WestPlains are displayed in Figure 4.15. Figure 4.15a displays the frequency model, andin plan view it is evident that the voxel stage alters the resistivity within the threeconductors in an attempt to concentrate areas of high conductivity in certain places.In cross-section, this effect is even more apparent as the Eastern conductor now has94a concentrated near-surface conductivity high similar to the time-domain result.Figure 4.15b shows the hybrid time-domain inversion, which only has minor de-viations from the parametric result. The voxel inversion places diffuse conductivezones around the main anomalies, but does not change the resistivity distributionwithin the conductors to any large degree.y(m)x(m)x(m)z(m)a)Ωmz(m)x(m) Frequency domainTime domainy(m)x(m)x(m)z(m)b)Ωmz(m)x(m)Figure 4.15: West Plains hybrid parametric inversions in plan view at an ele-vation of z = 190 m (left) and cross-section across solid white lines aty = 7333080 m and y = 7332780 m (right). a) Frequency domain. b)Time domain.The hybrid inversion also adds near-surface conductive zones due to clay layersand lakes as seen in Figure 4.16 where plan view slices at the model surface aredisplayed from both data sets. The colorbars have been modified to best displaynear-surface features, which are more pronounced in the frequency-domain modelcompared to the time-domain, as would be expected due to a greater near-surfacesensitivity.Modified conductivity structures added from the voxel stage of the hybridmethod permit a closer data fit compared to parametric inversions alone as pre-sented in Figure 4.14. Here the hybrid predicted data closely matches the observed95y(m)x(m)a)Ωmb)Figure 4.16: West Plains hybrid parametric inversions in plan view at an ele-vation of z = 310 m (surface). a) Frequency domain. b) Time domain.Note the colorbars are different to highlight near-surface features.data. Specifically, the double amplitude high in the real-component Hz data is nowpresent in the hybrid predicted data in Figure 4.14c. Hybrid misfit curves are dis-played in the latter half of Figure 4.12b where the voxel stage lowers the achievedmisfit from 6.27 to 1.87 in the frequency domain, and from 8.49 to 3.51 in thetime domain. Note that final hybrid data misfit values are marginally higher thanvoxel inversion misfits, but the hybrid model is judged to give a clearer interpre-tation of the shape and dip of the conductive anomalies. Originally, the startingregularization parameter β for the voxel portion of the hybrid inversion was cho-sen to be 1e-3 with a cooling parameter of γ = 0.2 like in Chapter 3. However, thissetup had trouble reducing the data misfit and the model update simply put diffuseconductivity zones around the central conductors with little change to the actualconductor shape. As a result, the hybrid inversions in this section use β = 0 for allGauss-Newton iterations.To compare the inversion models to ground truth information, Figure 4.17 dis-plays conductive zones from the West Plains inversions in the form of a 30 Ωmiso-surface along with 1 g/t gold mineralization within the Western conductive ko-matiite unit. Almost no drilling exists in the Eastern conductors, and the focusshifts primarily to the Western anomaly. A rock model of the komatiite unit it-self does not exist, therefore, the 1 g/t gold shape that is known to reside within96the komatiite acts as a proxy for highly conductive zones. A full model view anda closeup of the Western conductor is provided at y = 7333300 m. Figure 4.17ashows the frequency-domain voxel iso-surface, which provides a broad spatial cor-relation with the gold shape but no clear dip information. Figure 4.17b displays thetime-domain voxel inversion iso-surface where a better sense of the dip is achievedbut the anomaly is still much broader compared to the thin nature of the mineraliza-tion. The frequency and time-domain hybrid parametric inversion iso-surfaces aredisplayed in Figures 4.17c and 4.17d respectively. Both of these inversion demon-strate an excellent match with the compact dipping nature of the gold mineraliza-tion within the conductive komatiite, although the time-domain model achieves aslightly tighter spatial correlation.97a) b)x (m) x (m)y (m)y (m)c) d)478210479000479500480000480850310-130z (m)7333490733400073335007333000733295007331640 478210479000479500480000480850310-1307333490733400073335007333000733295007331640z (m)x (m) x (m)y (m)y (m)478210479000479500480000480850310-130z (m)73334907334000733350073330007332950073316407333490310z (m)-130480850480000479500479000478210733164073329500733300073340007333500Figure 4.17: 30Ωm iso-surfaces from voxel and hybrid inversions from WestPlains field data. 1 g/t gold shapes from drilling in the West conductorare shown in yellow and correspond with the conductive komatiite unit.a) Frequency-domain voxel. b) Time-domain voxel. c) Frequency-domain hybrid parametric. d) Time-domain hybrid parametric.984.5 ConclusionsA multiple body parametric inversion framework has been developed for frequencyand time-domain AEM data with unique resistivity values. The merits of the ap-proach are demonstrated through synthetic and field inversions where the anoma-lies of interest are thin dipping conductors. The targets are designed to representa typical orogenic gold exploration setting, such as the West Plains area of theCommittee Bay greenstone belt in the Canadian Arctic. The parametric inversionis coupled with a voxel approach to form a hybrid inversion in order to recoversharp targets in addition to background features. The parametric anomalies arecomposed of skewed Gaussian ellipsoids that can change location, radius and ro-tation angle in all three Cartesian planes. The synthetic example shows that theparametric stage of the hybrid inversion successfully recovers the spatial location,dip and resistivity of all three anomalies. The parametric results make interpretingthe dip and depth extent of the conductors much easier compared to conventionalvoxel-based approaches.In the West Plains field example, three near-vertical conductors are imagedwith the parametric algorithm. The frequency and time-domain field inversionscontain similarities with regards to the spatial extent of the targets in plan viewand the dip of the Western conductor in cross-section, although discrepancies existwhen examining the dip and depth extent of the two Eastern conductors. To achievean improved data fit and to fill in any remaining details missed by the parametricstage, a hybrid parametric technique is employed where the parametric results areused as initial and reference models for voxel-based inversions. Final hybrid resultshave improved agreement with the gold mineralization from drilling compared topure voxel inversions, and contain added detail not present within the parametricstage. Like the Caber example, no physical property information was available atWest Plains for constrained inversion, but would be beneficial for future inversions.As multiple models have been produced over a common spatial area at WestPlains, it is well suited for a cooperative approach to produce one consistent model,and this topic will be addressed in Chapter 5.99Chapter 5Cooperative Multiple AnomalyParametric Inversion of AEMDataIn this chapter, I apply both the cooperative and parametric methods to a scenariowhere spatially overlapping frequency and time-domain AEM data sets exist overmultiple thin dipping conductors. This approach is first tested with synthetic mea-surements before being applied to field data over the West Plains orogenic golddeposit. The cooperative parametric results highlight the benefits of this methodwhen data sets are compatible and provide complementary information, but alsoidentify potential pitfalls when data inconsistencies arise between surveys.5.1 IntroductionThe discovery of a natural resource deposit generally takes years if not decades,and over the course of that time frame, many rounds of geophysical surveys areoften collected. Due to changes in technology, budget, geologic target, or a myriadof other reasons, multiple AEM data sets can be flown over the same area. Eachof these data sets contains valuable information regarding the complex nature ofthe resistivity distribution of the target of interest. Chapter 2 introduced a cooper-ative technique for inverting overlapping surveys with voxel-based codes, but the100approach is equally applicable for parametric inversions where sharp anomaliesare sought with large resistivity contrasts compared to the background. This set-ting exists at the West Plains orogenic gold deposit discussed in Chapter 4, wheretwo overlapping frequency and time-domain AEM surveys were flown in differentyears over a greenstone geologic setting with multiple thin gold-bearing conduc-tors.The two research goals in this chapter are listed below.• To apply the cooperative approach to a parametric inversion framework.• To demonstrate and discuss the pros and cons of a cooperative parametricinversion with synthetic and field data.This chapter first revisits the cooperative algorithm from Chapter 2, adding afew modifications designed for the parametric problem. The cooperative paramet-ric method is then applied to overlapping synthetic frequency and time-domainAEM data introduced in Chapter4. Finally the process is implemented for over-lapping field data at West Plains in an attempt to produce one consistent resistivityinversion model for the area. Results are compared to models from Chapter 4 andto information from drilling, and the observations are discussed.5.2 Inversion MethodologyAs shown in Figure 2.8 the voxel-based cooperative algorithm operates by iteratingthrough each data set with a cooling schedule for the regularization parameter β ,and swapping the data back and forth between inversion codes. Recall the strengthof the cooperative method is its ability to rely on all possible sources of informationinstead of a single data set alone. However, since the parametric technique assignsβ = 0, the methodology is altered slightly to iterate between data sets after a fixednumber of Gauss-Newton iterations. Also, the cooperative inversion now passesparameters relating to the size and shape of the anomalies back and forth insteadof conductivity values as shown in Chapter 2. How many Gauss-Newton iterationsto choose depends on how well the data sets complement one another, with fewerneeded if the data misfit reduces steadily after each cooperative iteration, and moreif the data misfit has trouble decreasing. This research chooses to start with five101Gauss-Newton iterations for the synthetic data while a few choices are tested withfield measurements. Once this number has been selected, the model swappingcontinues until both data sets ideally hit the target misfit or the data misfits startincreasing. The cooperative parametric result can then be passed along as the initialand reference model to a voxel-based cooperative inversion where the process startsover as described in Chapter 2. The model at the end of this process is consideredthe final hybrid cooperative parametric model.5.3 Synthetic Results5.3.1 Cooperative parametric inversionThe cooperative parametric approach is tested on synthetic frequency and time-domain AEM data sets from Chapter 4, where the model consists of three dippingconductors in a resistive background. The parametric inversion, with β set to zero,begins with time-domain data, although results are similar by starting with fre-quency measurements. The process continues until the inversion reaches the endof the eighth cooperative iteration, when both the frequency and time-domain in-versions have a higher data misfit compared to the seventh iteration. The processterminates at this point and the result after the seventh cooperative iteration is con-sidered the final model. The data misfit progression is outlined in Figure 5.1 wheretwo black circles represent final data misfits. For clarity, five Gauss-Newton stepsare first taken for the time-domain data, numbered 1) on Figure 5.1. The param-eters from this result are used as starting guesses for five Gauss-Newton steps forthe frequency-domain data, numbered 2) on Figure 5.1. The parameters after thesecond segment are subsequently incorporated as starting values for the next fivetime-domain Gauss-Newton steps, as numbered 3) on Figure 5.1, and the processcontinues as such.Figure 5.2a displays in plan view and cross-section the true model, Figures 5.2band 5.2c re-shows the frequency and time-domain parametric results from Chap-ter 4, and Figure 5.2d portrays the cooperative parametric model. The cooperativeparametric inversion produces resistivity values of 0.16 Ωm and 1.87 Ωm respec-tively, and the shape and dip of the anomalies represent the closest approximation102Cooperative iterationData mist1)2)3)Figure 5.1: Synthetic cooperative parametric inversion data misfit progres-sion with the final model circled in black. The first three segments arenumbered for explanation purposes.to the true model of any of the synthetic inversions. Moreover, the quantitativeresidual for the cooperative model is the lowest number of any inversion at Rc= 2.15 for the conductive cells only and RAll = 0.18 for all cells. The syntheticinversion results collectively from Chapter 4 and Chapter 5 are summarized in Ta-ble 5.1. It is worth noting a comparable number of Gauss-Newton iterations arerequired per data set in the cooperative approach versus the individual parametricinversions, and final data misfit levels are also similar. Altogether, it is proposedthat the synthetic cooperative parametric inversion is an improvement over bothindividual parametric and voxel inversions.103y(m)x(m)y(m)x(m)x(m)z(m)a)c)ΩmΩmz(m)x(m)x(m)z(m)z(m)x(m)True ModelTime domainy(m)x(m)x(m)z(m)b)Ωmz(m)x(m) Frequency domainy(m)x(m)x(m)z(m)d)Ωmz(m)x(m) CooperativeFigure 5.2: Synthetic parametric inversions in plan view at z = 250 m (leftpanel) and cross-section across solid white lines at y = 7333080 mand y = 7332780 m (right panel). Dashed white lines represent start-ing guess locations for parametric anomalies. a) True synthetic model.b) Frequency-domain parametric inversion. c) Time-domain parametricinversion. d) Cooperative parametric inversion.104Table 5.1: Synthetic data parametric inversion summary.Inversion model ρ1 ρ2 Rc RAll # GN Final φdTrue model 0.20 Ωm 2.00 Ωm 0.0 0.0 - -Parametric (freq) 0.17 Ωm 2.62 Ωm 2.18 0.26 38 19.6Parametric (time) 0.15 Ωm 1.99 Ωm 3.33 0.31 28 17.6Parametric (coop) 0.16 Ωm 1.87 Ωm 2.15 0.18 35 18.7 (freq)35 27.0 (time)5.4 West Plains Case Study Results5.4.1 Cooperative parametric inversionsThe cooperative parametric method is now applied to West Plains AEM data in-troduced in Chapter 4, and once again inversion parameters are kept the same asbefore. However, as shown in Chapter 4, the best fitting half-spaces of 12,609 Ωmand 495 Ωm for the frequency and time-domain observations respectively suggestthat a discrepancy exists between the two data sets. Not surprisingly this createsproblems for the cooperative inversion, because during each cooperative iterationthe inversion drastically changes the background level, which takes away from thefocus of adjusting the shape of the target conductors.From research presented in Fournier et al. (2014); McMillan et al. (2015a) anddiscussed further in a full paper currently in revision, an analogous discrepancywas encountered over the Tli-Kwi-Cho (TKC) diamond exploration project in theNorth-West Territories, Canada. At TKC, two spatially overlapping AEM data setswere flown over a kimberlite complex, a 1992 DIGHEM frequency-domain (Foun-tain, 1998) and a 2004 VTEM time-domain survey. Similarly to West Plains, thetwo AEM data sets over TKC generate best-fitting half-spaces that are an order ofmagnitude apart. At TKC, the frequency-domain data also produces a more resis-tive half-space, which corroborates the observed discrepancy seen at West Plains.Literature states that unweathered igneous or metamorphic rock, representa-tive of background geology at TKC and West Plains, has an average backgroundresistivity of 10,000 Ωm (Palacky, 1988). This suggests that the frequency-domain105half-space value is more trust worthy. One possible explanation stems from thelack of transmitter waveform knowledge from 2003 or 2004 VTEM data, as theexact transmitter current used during the survey was not recorded in these years.Instead, it was presumed that the transmitter produced an ideal waveform with anelectric current level of zero during the off-time. Modern systems monitor the ac-tual transmitted waveform, which decays to zero during the early off-time. Thissmall time-varying current produces an additional primary field that increases theobserved ∂bz∂ t response measured by the receiver. Incorrectly assuming a currentlevel of zero in the off-time generates a spuriously elevated conductivity level inthe inversion model to compensate for the increased observed response, and wouldexplain the background level inconsistency at West Plains and TKC.Figure 5.3 illustrates this waveform effect. Here, ∂bz∂ t responses at the firstfive time channels from West Plains are modeled with both a practical and idealwaveform shown in Figure 5.3a. The data are modeled using a VTEM setup withthe transmitter and receiver 30 m above a 1000 Ωm half-space. 3% Gaussian noiseis added to the responses shown in Figure 5.3b, and it is clear that a small decayingcurrent in the off-time shifts the response curve up significantly. Although everytime-domain AEM system has a different current shut-off signature, this waveformis based off the 2012 VTEM waveform used in Chapter 3.With this data incompatibility in mind, the cooperative parametric algorithm isaltered to start with frequency-domain data, and after five Gauss-Newton iterationsthe best-fitting background level is fixed and the parameters are used as startingvalues for the time-domain data. To de-emphasize the early time channels, whichare the most likely to be contaminated by waveform issues, a 20% uncertainty as-signment is added to the first three channels as opposed to 3% for all others. Ahigher uncertainty is added to the first three time channels only, because in resis-tive regions, later time channels generally fall below the perceived noise thresholdand are discarded. In contrast, the effect of the unknown waveform is minimizedin conductive terrain as the response level is higher, and the data are deemed suit-able for a lower uncertainty level. Figure 5.4a shows the data misfit progressionwith the final model circled in black, which illustrates that after four cooperativeiterations both data misfits increase and the process stops. Even with adjustmentsto the background level and uncertainty assignments, the time-domain data misfit106Time (s) Time (s)Current (A)Response (V / Am2 )a) b)Figure 5.3: Comparison of a practical waveform with a small decaying cur-rent during early off-time ∂bz∂ t measurements compared to an ideal in-stantaneous shut-off waveform. a) Practical waveform (red) and idealwaveform (black) in early off-times. b) ∂bz∂ t responses at early timesfor the practical and ideal waveforms modeled 30 m above a 1000 Ωmhalf-space.has trouble decreasing with each progressive cooperative iteration. Therefore, asecond inversion is run with ten Gauss-Newton steps per cooperative iteration, andthe data misfit curve is displayed in Figure 5.4b. Once again, there is no great re-duction in misfit with subsequent cooperative iterations, and the process terminatesafter the fifth round of data swapping. It appears in this case that multiple cooper-ative iterations adds little value to the overall data fit, so a third trial is completedwith 25 Gauss-Newton steps. This inversion terminates after only two coopera-tive iterations, and the data misfit graph is displayed in Figure 5.4c. Although thetime-domain data misfit noticeably decreases in the first cooperative iteration, theincompatibility between data sets is evident by the extremely high initial misfitvalue at the start of the second frequency-domain cooperative iteration. An alter-native option is to let each stage run as many Gauss-Newton iterations as possiblebefore swapping, but this process not only takes longer to compute, but also pro-duces a similar model to the 25 Gauss-Newton iteration case; consequently, thisresult is not shown.At this point, no cooperative inversion has shown an ability to match both data107Cooperative iterationData mista) Cooperative iterationData mistb)c) Cooperative iteration d) Cooperative iterationData mistData mist1)2)3)4)5)Figure 5.4: West Plains cooperative parametric inversion data misfit progres-sions with varying numbers of Gauss-Newton iterations per cooperativeiteration. Final models circled in black. a) 5 Gauss-Newton iterations.b) 10 Gauss-Newton iterations. c) 25 Gauss-Newton iterations. d) Hy-brid cooperative inversion with 25 Gauss-Newton iterations, with eachsegment numbered for explanation purposes.sets sufficiently, and the relative importance of each survey must be evaluated.Frequency-domain measurements provide the best estimate for the background re-sistivity, but time-domain data have more time channels with a greater depth ofpenetration. Therefore, even with an uncertain waveform, it is deemed that fittingtime-domain data will provide the most information with relation to the shape, dipand resistivity of the target conductors. With this in mind, the cooperative para-metric result using 25 Gauss-Newton steps is chosen as the final parametric modelsince it produces the lowest time-domain data misfit.1085.4.2 Hybrid cooperative parametric inversionA hybrid inversion is now completed using the final cooperative parametric result.This involves computing a further 25 Gauss-Newton steps for each data set withvoxel-based codes. Figure 5.4d shows the data misfit progression, where every seg-ment is numbered for clarity. The first segment represents the frequency-domainparametric stage, and the parameters after this segment are passed along as startingguesses for the time-domain parametric stage, numbered 2) on Figure 5.4d. Thethird segment represents the frequency-domain voxel inversion, the latter half ofthe hybrid method, using the parameters from the second stage. The fourth andfinal segment represents the time-domain voxel inversion using the result from thethird segment as an initial model. After this stage, the time-domain data misfitreaches a low point of 1.56, and this is considered the final hybrid cooperativeparametric model. Like in Chapter 4, the regularization parameter β is set to zerofor the duration of this process and all model cells are kept active in the voxelstage. This model still does not accurately fit the frequency-domain data as shownby the subsequent initial data misfit of 150, numbered 5) on Figure 5.4d, but this isconsidered a consequence of incompatible data sets and the goal of fitting the time-domain measurements is deemed a success. Note that with only one cooperativeiteration, this is a similar style of sequential voxel inversion as used for electromag-netic and magnetotelluric data in Commer and Newman (2009). The difference inthis example is the use of a parametric algorithm to provide the initial guess for thecooperative voxel inversion.Plan view slices and cross-sections from the cooperative parametric model arepresented in Figure 5.5a and through the final hybrid cooperative parametric resultin Figure 5.5b. The similarities between the parametric and hybrid models are ob-vious, but in plan view the hybrid approach outlines subtle concentrated areas ofconductivity highs within the parametric bodies, which could be helpful to targetareas of greatest mineralization. In both parametric and hybrid models, the Westernconductor has a steep vertical dip and a comparable shape to individual parametricmodels of Chapter 4. In contrast, the Eastern conductors are both wider com-pared to parametric results, but the central conductor maintains a similar steep dipto the West. The Eastern conductor changes the most and possesses an elongate109near-vertical shape instead of being a near-surface compact body. This anomalyrepresents the most uncertain target, especially at depth, but as a whole this hybridmodel can be used to make new interpretations regarding the location and dip ofthe Eastern conductors.With such a high final data misfit value, it may not seem like the frequency-domain data contributes much to the cooperative process apart from the back-ground value, but it provides a valuable starting model, making it easier to fittime-domain responses. The frequency measurements also provide near-surfaceinformation not present in the time domain data, as seen in a surficial view of thehybrid parametric model in Figure 5.5c. Once again, some caution must be takeninterpreting these surface anomalies due to the poor fit of the frequency data, butthere is a general correlation with these features to near-surface clay layers in lakes.110y(m)x(m)x(m)a)Ωmx(m) Cooperative parametricy(m)x(m)x(m)b)Ωmx(m)z(m)z(m)z(m)z(m)Hybrid cooperative parametricy(m)x(m)Ωmc)Figure 5.5: Cooperative and hybrid cooperative parametric inversions at WestPlains. Dashed white lines represent starting guess locations for para-metric anomalies. Plan view slices at z = 190 m (left) and cross-sectionacross solid white lines at y = 7333080 m and y = 7332780 m (right).a) Cooperative parametric model. b) Hybrid cooperative parametricmodel. c) Hybrid cooperative parametric model with a plan view sliceat surface.1115.5 Cooperative Inversion DiscussionThe purpose of the cooperative approach is to find one consistent 3D resistivitymodel from two overlapping AEM data sets at West Plains, but this proved prob-lematic as no one model matched this criteria. The incompatibility between datasets led to fixing the background resistivity at the best-fitting level derived fromfrequency-domain data and focusing on fitting time-domain observations there-after. Figure 5.6a portrays the 30 Ωm iso-surface through the final hybrid co-operative parametric model. A closeup view in Figure 5.6b demonstrates a highcorrelation with gold mineralization within the Western komatiite. Figure 5.6c cir-cles an area of widening in the Western conductor at an elevation slice of z = 280m that corresponds to a zone of concentrated gold mineralization. This is highlybeneficial information achieved through the hybrid approach. Finally Figure 5.6dshows the excellent correlation between the dip of the Western conductor and goldmineralization.112b)x (m)y (m)a)c)-130310479000479500480000481010478130733491473340007333500733300073325007331560x (m)d)z (m)ΩmFigure 5.6: Final hybrid cooperative parametric inversion at West Plains with1 g/t gold shapes in yellow. a) 30 Ωm iso-surface. b) Closeup view ofWestern conductor with 30Ωm iso-surface. c) Closeup view of Westernconductor with 280 m elevation slice. Orange circle highlights areaof condensed mineralization that corresponds with a widening of therecovered conductor. d) Cross section at y = 7333230 m.113Even though the recovered model extends to a greater depth compared to thegold shapes, this deeper zone contains higher uncertainty as few drill holes reachedthis depth, and mineralization could extend further. Overall, the accurate corre-lations with gold mineralization in the Western conductor suggests that valuableinterpretations can now be made regarding the Eastern conductors. To this end, a190 m elevation slice from the final hybrid cooperative parametric model is over-laid on the simplified geology map in Figure 5.7. The recovered conductors aremuch thinner compared to the mapped komatiites, and new spatial outlines in planview and cross-section can be applied to update the geology image.KomatiiteSedimentsGranodioriteTonaliteSurvey outlineFigure 5.7: Final hybrid cooperative parametric inversion at West Plains. El-evation slice at z = 190 m overlaid on simplified geology map.This field example is meant to demonstrate a few difficulties that can be en-countered with overlapping field data that do not necessarily agree. Weightingschemes could also be applied to down-weight one data set over the other such asdescribed in Meqbel and Ritter (2014) and Commer and Newman (2009); alter-natively, an ADMM approach (Wahlberg et al., 2012) could be implemented. AtWest Plains, certain elements of each data set are deemed important, such as the114frequency-domain background level and time-domain information from later timechannels. Finding the optimal way to treat overlapping incompatible data sets isstill an open question, but this chapter sheds some light on ways to progress for-ward in this scenario. Fortunately, time-domain AEM data in recent years havestarted recording the exact waveform transmitted, and the background level issuepresented here should not be as prevalent in present-day and future surveys.5.6 ConclusionsA cooperative parametric approach is applied to spatially overlapping frequencyand time-domain AEM data sets where the targets of interest are thin dipping con-ductors. The method is tested on both synthetic and field examples. In the syntheticexample, the cooperative method produces one consistent model that achieves acloser match to the true model compared to individual parametric or voxel inver-sions.In the West Plains field example, the cooperative method fails to find one modelthat accurately fits both data sets, and it is deemed that the overlapping AEM sur-veys contain incompatible measurements. This is possibly due to a small decayingcurrent in the transmitter waveform for time-domain data during off-time mea-surements. The proposed solution is to use frequency-domain data to generate aninitial parametric model with a best-fitting background resistivity, and to fix thisbackground level for the time-domain parametric inversion stage. This process isthen repeated with voxel based algorithms to produce a final hybrid cooperativeparametric model.The corresponding agreement to gold mineralization is excellent, and new in-terpretations can now be made with regards to the shape of the Eastern conductors.Although this is perceived to be the best possible model from this data, it is ac-knowledged that the final data fit on the frequency-domain data is still poor andmore research is needed to fully rectify this discrepancy.115Chapter 6ConclusionsAirborne electromagnetic data remains one of the most commonly used tools forgeophysical exploration as a means to image the subsurface of the earth. Thismethod is sensitive to contrasts in electrical resistivity, which can be used to helplocate and characterize buried natural resources such as metal, water and oil. Whilesurveys have been flown for decades, 3D AEM inversion is a relatively new andrapidly developing field. As such, I designed this thesis to contribute to the scien-tific community by advancing 3D AEM inversion accuracy over a wide range ofgeologic environments. In particular, I present practical strategies that address twoimportant and unresolved questions.1. How to improve 3D AEM inversion with spatially overlapping surveys.2. How to improve 3D AEM inversion to recover thin, high contrast anomalies.The following sections summarize the findings of this thesis followed by adiscussion on future research questions.6.1 Cooperative InversionIn Chapter 2 I investigated the first topic listed above, and looked at how to invertmultiple overlapping electromagnetic data sets to produce one consistent 3D re-sistivity model that satisfied all geophysical observations. Three overlapping fielddata sets, time-domain AEM, CSAMT and DC, were examined over the Antonio116high-sulfidation epithermal deposit in Peru. After comparing inversion results fromeach individual data set, I found that resulting models contained many discrepan-cies in their resistivity distributions. This phenomenon was studied in more detailthrough a synthetic example which was designed to emulate the geology at Antoniowith a large resistor that hosted two buried conductors. The individual syntheticinversions were able to recover similar looking resistive regions, but the detectionof the contained conductive bodies varied greatly. A cooperative approach wassubsequently developed, where the inversion iterated between data sets and con-structed a result that satisfied all geophysical information provided. This methodcreated one consistent 3D model that successfully imaged the resistive region aswell as the conductive bodies to a closer degree to the true model compared to anyindividual inversion.Following the success of the synthetic example, I applied the cooperative tech-nique to the Antonio field data and generated a common physical property modelthat mapped a large gold-bearing silica-rich resistor interspersed with smaller con-ductive zones. The final cooperative model resolved this resistive region with ahigh level of agreement compared to known silica alteration from drilling and ge-ologic mapping, and highlighted small conductive regions that could be indicativeof sulfide and gold mineralization. These inferences represented new potentialtargets of interest at this project, and showed the applicability of the cooperativeapproach to a field setting. This chapter also demonstrated how physical prop-erty information derived from drill holes assisted in constraining both the syntheticand field cooperative inversion models through upper and lower resistivity bounds.This showed that drill hole information should be incorporated into cooperativeinversions when available.6.2 Parametric Inversion - Single AnomalyChapter 3 addressed the second thesis topic of how to improve 3D AEM inver-sion over thin conductive targets. Here I developed a parametric solution wherethe algorithm recovered a single anomaly in the shape of a best-fitting skewedGaussian ellipsoid that represented the target of interest in a smooth background.The method was first tested on a synthetic example composed of a thin conductive117plate-like anomaly in a uniform resistive background, and the parametric methodwas able to recover the target to a much higher degree of accuracy compared to apurely voxel-based solution.The approach was then applied to time-domain AEM field data over the Cabervolcanogenic massive sulfide deposit in Quebec where the parametric inversionmodel was able to replicate the thin dipping nature of the copper-zinc resource asdefined by drilling. The Caber parametric model was subsequently used as a refer-ence model for a voxel-based inversion to fill in missing features such as a conduc-tive overburden layer. The hybrid parametric inversion was designed to combineelements from both parametric and voxel methods, and this proved beneficial tomap out the Caber target in addition to relevant surrounding geology.6.3 Parametric Inversion - Multiple AnomaliesIn Chapter 4 I modified the parametric framework to include multiple skewedGaussian ellipsoids with individual resistivity values. This extension made theparametric approach suitable for a wider range of geologic settings and targetstyles. The multiple anomaly parametric code was first tested on synthetic fre-quency and time-domain AEM data over three narrow dipping conductors in aresistive background. This synthetic model was built to replicate an orogenic goldexploration setting such as the West Plains area of the Committee Bay greenstonebelt in Nunavut. The synthetic results recovered the shape, dip and resistivity valueof the dipping conductors with excellent precision, and made interpreting the mod-els much simpler compared to voxel-based inversions.The technique was subsequently tested on the West Plains frequency and time-domain AEM field data, and recovered three narrow, elongate, near-vertical con-ductors. The Western anomaly from both parametric models displayed a high levelof correlation to the shape and dip of gold mineralization hosted in a conductive ko-matiite unit known from drilling, and managed to depict its thin nature with greaterease compared to conventional voxel inversions. The Eastern conductors were lessconsistent between the parametric models and were not as thoroughly evaluatedsince little ground truth was available in this area. I applied the hybrid approachto the parametric field results which helped fit the observed data to a higher degree118and added near-surface features not possible with a parametric method alone. Thehybrid inversion also generated concentrated areas of conductivity within the targetanomalies, which could be valuable to detect high grade mineralization.6.4 Cooperative Parametric InversionFor Chapter 5 I combined the cooperative and parametric methods for spatiallyoverlapping AEM data sets. This approach saw improved inversion accuracy withsynthetic data, but encountered issues with incompatible field data. With syntheticdata from Chapter 4, the cooperative parametric result imaged three dipping con-ductors with greater precision compared to individual parametric inversions. Thelatter were already a close approximation to the synthetic model, but the coopera-tive method was able to extract additional information and created subtle improve-ments. As these data were computer generated, they represented a test examplewith perfectly calibrated and compatible data.The cooperative parametric algorithm was subsequently tested on overlappingfrequency and time-domain AEM field measurements at West Plains, and it be-came quickly apparent that the data sets were incompatible and did not providefully complementary information. The biggest discrepancy came in the form ofthe best-fitting recovered background resistivity, where the frequency observationsproduced a highly resistive value and the time-domain data suggested a muchmore conductive background. A similar inconsistency was previously identifiedat the TKC kimberlite deposit in the Canadian Arctic with comparable AEM sys-tems. It is postulated to be caused from the time-domain transmitter possessingan unknown non-zero decaying current level during the system off-time. Con-sequently, it was determined that the frequency-domain data produced the moreaccurate background value and this was fixed for time-domain iterations.The final hybrid cooperative parametric model at West Plains imaged the West-ern conductor with excellent agreement to drilling information, and was able toisolate areas of high conductivity within the Western conductor that correspondedwith areas of concentrated gold mineralization. Similar interpretations can nowbe made for the shape and dip of the two Eastern conductors where little drillingknowledge exists. However, as the surveys were deemed incompatible, the final119model was unable to fit both data sets, and in this example the priority was givento time-domain measurements. This was due to more available time channels com-pared to frequencies coupled with a greater depth of penetration. However, thefrequency-domain data were still beneficial for providing the background value andnear-surface features. Collectively this chapter showed that a hybrid cooperativeparametric method can be advantageous when the observations provide comple-mentary information, but caution must be taken when compatibility issues arisebetween data sets.6.5 Future ResearchIn this thesis I have introduced practical cooperative and parametric strategies for3D AEM inversion, both in the frequency and time domain, but further researchis needed to optimize algorithm robustness and efficiency. The ADMM method(Wahlberg et al., 2012) would be worthwhile testing for multiple spatially over-lapping data sets and comparing the results against the cooperative algorithm pre-sented. In terms of the parametric algorithm, future trials should raise the com-plexity level of the target geometry to push the limits of the inversion capability.It will be important to include small-scale compact targets and also large-scalesmoothly varying features to thoroughly test the hybrid parametric method. Mul-tiple data sets can be synthetically generated over these complex models to testthe hybrid cooperative parametric method, and additional field data sets should besought for practical applications. These future trials should incorporate alternativeparameterizations such as multiple spheres and potentially other relevant shapes inorder to evaluate a broad scope of parametric methods. Finally this dissertation hasfocused on airborne electromagnetic data but the cooperative and parametric strate-gies developed are designed to be general, and are equally applicable to ground EMmeasurements, potential field or induced-polarization surveys.120BibliographyAbubakar, A., T. M. Habashy, V. Druskin, L. Knizhnerman, and S. Davydycheva,2006, A 3D parametric inversion algorithm for triaxial induction data:Geophysics, 71, G1–G9. → pages 44Adair, R., 2011, Technical Report on the Resource Calculation for the MatagamiProject , Que´bec: Technical report. → pages xiv, 63, 64Aghasi, A., M. Kilmer, and E. L. Miller, 2011, Parametric Level Set Methods forInverse Problems: SIAM Journal on Imaging Sciences, 4, 618–650. → pages44, 45, 60, 129Albouy, Y., P. Andrieux, G. Rakotondrasoa, M. Ritz, M. Descloitres, J.-L. Join,and E. Rasolomanana, 2001, Mapping coastal aquifers by joint inversion of DCand TEM soundings-three case histories: Groundwater, 39, 87–97. → pages 11Allard, M., 2007, On the Origin of the HTEM Species in B. Milkereit ed.:Proceedings of Exploration 07: Fifth Decennial International Conference onMineral Exploration, 355–374. → pages 19, 46, 56, 65, 78, 81Arribas, A., 1995, Characteristics of high-sulfidation epithermal deposits, andtheir relation to magmatic fluid: Mineralogical Association of Canada ShortCourse Series, 23, 419–454. → pages 16Ascher, U., and F. Roosta-khorasani, 2016, Algorithms that Satisfy a StoppingCriterion , Probably: Vietnam Journal of Mathematics, 44, 49–69. → pages 44Bezanson, J., S. Karpinski, V. B. Shah, and A. Edelman, 2012, Julia: A FastDynamic Language for Technical Computing: arXiv preprint arXiv:1209.5145.→ pages 10Boyko, W., S. Balch, and N. Paterson, 2003, The AeroTEM airborneelectromagnetic system: The Leading Edge, 22, 562–566. → pages 4, 81Carr, P. M., L. M. Cathles III, and C. T. Barrie, 2008, On the size and spacing ofvolcanogenic massive sulfide deposits within a district with application to theMatagami district, Quebec: Economic Geology, 103, 1395–1409. → pages 63Carson, C. J., R. G. Berman, R. A. Stern, M. Sanborn-barrie, T. Skulski, andH. A. I. Sandeman, 2004, Age constraints on the Paleoproterozoictectonometamorphic history of the Committee Bay region , western Churchill121Province , Canada : evidence from zircon and in situ monazite SHRIMPgeochronology: Canadian Journal of Earth Sciences, 41, 1049–1076. → pages87Caudillo-Mata, L. A., E. Haber, and C. Schwarzbach, 2016, An oversamplingtechnique for the multiscale finite volume method to simulate electromagneticresponses in the frequency domain: arXiv preprint arXiv:1610.02105, 1–15. →pages 26Chen, T., G. Hodges, and P. Miles, 2015, MULTIPULSE high resolution and highpower in one TDEM system: Exploration Geophysics, 46, 49–57. → pages 4Commer, M., and G. A. Newman, 2004, A parallel finite-difference approach for3D transient electromagnetic modeling with galvanic sources: Geophysics, 5,1192–1202. → pages 7——–, 2009, Three-dimensional controlled-source electromagnetic andmagnetotelluric joint inversion: Geophysical Journal International, 178,1305–1316. → pages 10, 15, 109, 114Constable, S. C., R. L. Parker, and C. G. Constable, 1987, Occam’s inversion: Apractical algorithm for generating smooth models from electromagneticsounding data: Geophysics, 52, 289–300. → pages 9Cox, L. H., G. A. Wilson, and M. S. Zhdanov, 2010, 3D inversion of airborneelectromagnetic data using a moving footprint: Exploration Geophysics, 41,250–259. → pages 7, 9——–, 2012, 3D inversion of airborne electromagnetic data: Geophysics, 77,WB59–WB69. → pages 9, 10, 43Dorn, O., and D. Lesselier, 2006, Level set methods for inverse scattering: InverseProblems, 22, R67–R131. → pages 44, 60Dorn, O., E. L. Miller, and C. M. Rappaport, 2000, A shape reconstructionmethod for electromagnetic tomography using adjoint fields and level sets:Inverse Problems, 16, 1119–1156. → pages 44Eaton, P., 1998, Application of an improved technique for interpreting transientelectromagnetic data: Exploration Geophysics, 29, 175–183. → pages 9Eaton, P. A., R. G. Anderson, S. V. Queen, B. Y. Nilsson, E. Lauritsen, C. T.Barnett, M. Olm, and S. Mitchell, 2013, Helicopter time-domainelectromagnetics - Newmont and the NEWTEM experience: Geophysics, 6,W45–W56. → pages 4, 19, 81Egbert, G. D., and A. Kelbert, 2012, Computational recipes for electromagneticinverse problems: Geophysical Journal International, 189, 251–267. → pages 7Ekblom, H., 1973, Calculation of linear best Lp approximations: BIT NumericalMathematics, 13, 292–300. → pages 8, 57EMIT, 2005, Maxwell 4.0: Modeling, presentation and visualization of EM andelectrical geophysical data: Technical report. → pages 11, 66122Farquharson, C. G., and D. W. Oldenburg, 1998, Non-linear inversion usinggeneral measures of data misfit and model structure: Geophysical JournalInternational, 134, 213–227. → pages 44Farquharson, G., and D. W. Oldenburg, 1993, Inversion of time-domainelectromagnetic data for a horizontally layered Earth: Geophysical JournalInternational, 114, 433–442. → pages 9Fountain, D., 1998, Airborne electromagnetic systems - 50 years of development:Exploration Geophysics, 29, 1–11. → pages 105Fournier, D., 2015, A Cooperative Magnetic Inversion Method With Lp-normRegularization: PhD thesis. → pages 8, 44Fournier, D., L. Heagy, N. Corcoran, D. Cowan, S. G. R. Devriese, D. Bild-enkin,K. Davis, S. Kang, D. Marchant, M. S. McMillan, M. Mitchell, G. Rosenkjar,and D. Yang, 2014, Multi-EM systems inversion - Towards a commonconductivity model for the Tli Kwi Cho complex: 84th Annual InternationalMeeting, SEG, Expanded Abstracts, 33, 1795–1799. → pages 105Fraser, D. C., 1978, Resistivity mapping with an airborne multicoilelectromagnetic system: Geophysics, 43, 144–172. → pages 9Gallardo, L. A., and M. A. Meju, 2004, Joint two-dimensional DC resistivity andseismic travel time inversion with cross-gradients constraints: Journal ofGeophysical Research, 109, B03311. → pages 11Gill, P. E., W. Murray, and M. H. Wright, 1981, Practical Optimization: AcademicPress Inc. → pages 48Goldie, M. K., 2000, A geophysical case history of the Yanacocha gold district,northern Peru: 70th Annual International Meeting, SEG, Expanded Abstracts.→ pages 17, 18Haber, E., and U. M. Ascher, 2001, Fast finite volume simulation of 3Delectromagnetic problems with highly discontinuous coefficients: Siam Journalof Scientific Computing, 22, 1943–1961. → pages 80Haber, E., U. M. Ascher, and D. W. Oldenburg, 2004, Inversion of 3Delectromagnetic data in frequency and time domain using an inexact all-at-onceapproach: Geophysics, 69, 1216–1228. → pages 23, 46Haber, E., and S. Heldmann, 2007, An octree multigrid method for quasi-staticMaxwell’s equations with highly discontinuous coefficients: Journal ofComputational Physics, 223, 783–796. → pages 46Haber, E., S. Heldmann, and U. Ascher, 2007a, Adaptive finite volume method fordistributed non-smooth parameter identification: Inverse Problems, 23,1659–1676. → pages 22, 46Haber, E., D. W. Oldenburg, and R. Shekhtman, 2007b, Inversion of time domainthree-dimensional electromagnetic data: Geophysical Journal International,171, 550–564. → pages 7, 46123Haber, E., and C. Schwarzbach, 2014, Parallel inversion of large-scale airbornetime-domain electromagnetic data with multiple OcTree meshes: InverseProblems, 30, 1–28. → pages 7, 9, 10, 19, 22, 23, 43, 46, 48, 51, 52, 57Hoschke, T., 2011, Geophysical signatures of copper-gold porphyry andepithermal gold deposits, and implications for exploration: PhD thesis,University of Tasmania. → pages xi, 17Isakov, V., S. Leung, and J. Qian, 2011, A Fast Local Level Set Method forInverse Gravimetry 1 Introduction: Communications in Computational Physics,10, 1044–1070. → pages 44Keating, P. B., and D. J. Crossley, 1990, The inversion of time-domain airborneelectromagnetic data using the plate model: Geophysics, 55, 705–711. →pages 11Keller, G. V., 1988, Rock and mineral properties: Electromagnetic methods inapplied geophysics, 1, 13–52. → pages 2, 43Kerswill, J. A., 1996, Iron-formation-hosted stratabound gold: Geology ofCanadian mineral deposit types. Geology of Canada, 8, 367–382. → pages 77Key, K., and C. Weiss, 2006, Adaptive finite-element modeling using unstructuredgrids: The 2D magnetotelluric example: Geophysics, 71, G291–G299. →pages 7, 10Kunetz, G., 1966, Principles of direct current resistivity prospecting: GebruderBorntraeger, Berlin. → pages 4Last, B. J., and K. Kubik, 1983, Compact gravity inversion: Geophysics, 6,713–721. → pages 8, 44Lelie`vre, P. G., and D. W. Oldenburg, 2009, A comprehensive study of includingstructural orientation information in geophysical inversions: GeophysicalJournal International, 178, 623–637. → pages 11Lelie`vre, P. G., D. W. Oldenburg, and N. C. Williams, 2009, Integratinggeological and geophysical data through advanced constrained inversions *:Exploration Geophysics, 40, 334–341. → pages 11Li, Y., and K. Key, 2007, 2D marine controlled-source electromagnetic modeling :Part 1 An adaptive finite-element algorithm: Geophysics, 72, WA51–WA62.→ pages 7Li, Y., and D. W. Oldenburg, 2000, Incorporating geological dip information intogeophysical inversions: Geophysics, 65, 148–157. → pages 11Lines, L. R., A. K. Schultz, and S. Treitel, 1988, Cooperative inversion ofgeophysical data: Geophysics, 53, 8–20. → pages 10, 15Loayza, C., and J. Reyes, 2010, Geologic overview of the Yanacocha miningdistrict: Presented at the SIMEXMIN10 - IV Brazilian symposium of mineralexploration. → pages 17Lu, W., and J. Qian, 2015, A local level-set method for 3D inversion of124gravity-gradient data: Geophysics, 80, G35–G51. → pages 44Macnae, J., A. King, A. Osmakoff, and A. Blaha, 1998, Fast AEM dataprocessing and inversion: Exploration Geophysics, 29, 163–169. → pages 9Macnae, J., and Y. Lamontagne, 1987, Imaging quasi-layered conductivestructures by simple processing of transient electromagnetic data: Geophysics,52, 545–554. → pages 9Macnae, J. C., R. Smith, B. D. Polzer, Y. Lamontagne, and P. S. Klinkert, 1991,Conductivity-depth imaging of airborne electromagnetic step-response data:Geophysics, 56, 102–114. → pages 3, 9McMillan, M. S., and D. W. Oldenburg, 2012, Three-dimensional electromagneticand electrical inversions over the Antonio gold deposit in Peru: 82nd AnnualInternational Meeting, SEG, Expanded Abstracts, 31, 1316–1320. → pages iv——–, 2014, Cooperative constrained inversion of multiple electromagnetic datasets: Geophysics, 79, B173–B185. → pages iv, xi, 16McMillan, M. S., D. W. Oldenburg, E. Haber, and C. Schwarzbach, 2015a,Parametric 3D inversion of airborne time domain electromagnetics: ASEGExpanded Abstracts, 1–5. → pages 70, 105McMillan, M. S., D. W. Oldenburg, E. Haber, C. Schwarzbach, and E. Holtham,2015b, 3D Multiple Body Parametric Inversion of Time-domain Airborne EMData: First European Airborne Electromagnetics Conference, 1–4. → pages ivMcMillan, M. S., C. Schwarzbach, E. Haber, and D. W. Oldenburg, 2015c, 3Dparametric hybrid inversion of time-domain airborne electromagnetic data:Geophysics, 80, K25–K36. → pages iv——–, 2016, 3D Cooperative Parametric Inversion of Frequency andTime-domain Airborne Electromagnetic Data: Presented at the FirstConference on Geophysics for Mineral Exploration and Mining. → pages vMcMillan, M. S., C. Schwarzbach, D. W. Oldenburg, E. Haber, E. Holtham, andA. Prikhodko, 2014, Recovering a thin dipping conductor with 3Delectromagnetic inversion over the Caber deposit: 84th Annual InternationalMeeting, SEG, Expanded Abstracts, 33, 1720–1724. → pages ivMeqbel, N., and O. Ritter, 2014, New Advances for a Joint 3D Inversion ofMultiple EM Methods: 76th EAGE Conference and Exhibition-Workshops,1–4. → pages 114Modersitzki, J., 2003, Numerical Methods for Image Registration: Oxforduniversity press. → pages 49Mule`, S., R. Miller, H. Carey, and R. Lockwood, 2012, Review of three airborneEM systems: ASEG Expanded Abstracts, 26–29. → pages 4, 81Nabighian, M. N., and J. C. Macnae, 1991, Time domain electromagneticprospecting methods: Electromagnetic methods in applied geophysics, 2,427–479. → pages 9, 20125Nelson, P. H., and G. D. Van Voorhis, 1983, Estimation of sulfide content frominduced polarization data: Geophysics, 48, 62–75. → pages 38Oldenburg, D. W., R. A. Eso, S. Napier, and E. Haber, 2005, Controlled sourceelectromagnetic inversion for resource exploration: First Break, 23, 41–48. →pages 17, 23, 39Oldenburg, D. W., E. Haber, and R. Shekhtman, 2013, Three dimensionalinversion of multisource time domain electromagnetic data: Geophysics, 78,E47–E57. → pages 7, 9, 44, 80Oldenburg, D. W., and Y. Li, 1994, Inversion of induced-polarization data:Geophysics, 59, 1327–1341. → pages 23Oldenburg, D. W., Y. Li, and R. G. Ellis, 1997, Inversion of geophysical data overa copper gold porphyry deposit : A case history for Mt. Milligan: Geophysics,62, 1419–1431. → pages 10, 11, 15Oldenburg, D. W., R. Shehktman, R. A. Eso, C. G. Farquharson, P. Eaton, B.Anderson, and B. Bolin, 2004, Closing the gap between research and practicein EM data interpretation: 74th Annual International Meeting, SEG, ExpandedAbstracts, 23, 1179–1183. → pages 23, 39Osher, S., and R. P. Fedkiw, 2001, Level Set Methods: An Overview and SomeRecent Results: Journal of Computational Physics, 169, 463–502. → pages 44Osher, S., and J. A. Sethian, 1988, Fronts Propagating with Curvature- DependentSpeed : Algorithms Based on Hamilton-Jacobi Formulations: Journal ofComputational Physics, 79, 12–49. → pages 44, 51, 75Palacky, G. J., 1988, Resistivity Characteristics of Geologic Targets:Electromagnetic methods in applied geophysics, 1, 53–129. → pages 66, 105Palacky, G. J., and G. F. West, 1973, Quantitative interpretation of AEMmeasurements: Geophysics, 38, 1145–1158. → pages 3, 9Pidlisecky, A., K. Singha, and F. D. Day-Lewis, 2011, A distribution-basedparametrization for improved tomographic imaging of solute plumes:Geophysical Journal International, 187, 214–224. → pages 45, 129Portniaguine, O., and M. S. Zhdanov, 1999, Focusing geophysical inversionimages: Geophysics, 64, 874–887. → pages 8, 44Prikhodko, A., T. Eadie, N. Fiset, and L. Mathew, 2012, VTEM (35 m) test resultsover the Caber and Caber North Deposits (October 2012): Technical report. →pages xiv, 66, 69, 71Prikhodko, A., E. Morrison, A. Bagrianski, P. Kuzmin, P. Tishin, and J. Legault,2010, Evolution of VTEM technical solutions for effective exploration: ASEGExpanded Abstracts, 1, 1–4. → pages xiv, 63, 64Raiche, A., 1998, Modelling the time-domain response of AEM systems:Exploration Geophysics, 29, 103. → pages 11, 44Raiche, A. P., D. L. B. Juppt, H. Rutter, and K. Vozoff, 1985, The joint use of126coincident loop transient electromagnetic and Schlumberger sounding toresolve layered structures: Geophysics, 50, 1618–1627. → pages 9Ruthotto, L., E. Treister, and E. Haber, 2016, Jinv a flexible julia package for pdeparameter estimation: arXiv preprint arXiv:1606.07399, 1–19. → pages 7, 9,10, 44Schwarzbach, C., and E. Haber, 2013, Finite element based inversion fortime-harmonic electromagnetic problems: Geophysical Journal International,193, 615–634. → pages 7Skulski, T., H. A. I. Sandeman, and M. Sanborn-barrie, 2003, Bedrock geology ofthe Ellice Hills map area and new constraints on the regional geology of theCommittee Bay area: Natural Resources Canada, Geological Survey of Canada,1–11. → pages 87Smith, R., D. Fountain, and M. Allard, 2003, The MEGATEM fixed-wingtransient EM system applied to mineral exploration: A discovery case history.:First Break, 21, 73–77. → pages 81Sørensen, K. I., and E. Auken, 2004, SkyTEM a new high-resolution helicoptertransient electromagnetic system: Exploration Geophysics, 35, 194–202. →pages 81Sosa, A., A. A. Velasco, L. Velazquez, M. Argaez, and R. Romero, 2013,Constrained optimization framework for joint inversion of geophysical datasets: Geophysical Journal International, 195, 1745–1762. → pages 11Steklova, K., and E. Haber, 2016, Joint hydrogeophysical inversion: stateestimation for seawater intrusion models in 3D: Computational Geosciences,1–20. → pages 10Sturler, E. D. E., and M. E. Kilmer, 2011, A regularized GaussNewton trust regionapproach to imaging in optical tomography: Siam Journal of ScientificComputing, 33, 3057–3086. → pages 51Tai, X.-c., and T. F. Chan, 2004, A survey on multiple level set methods withapplications for identifying piecewise constant functions: International Journalof Numerical Analysis and Modeling, 1, 25–47. → pages 51, 75Teal, L., and A. Benavides, 2010, History and geologic overview of theYanacocha mining district, Cajamarca, Peru: Economic Geology, 105,1173–1190. → pages 16, 17van den Doel, K., and U. Ascher, 2006, On level set regularization for highlyill-posed distributed parameter estimation problems: Journal of ComputationalPhysics, 216, 707–723. → pages 44, 60Viezzoli, A., E. Auken, and T. Munday, 2009, Spatially constrained inversion forquasi 3D modelling of airborne electromagnetic data an application forenvironmental assessment in the Lower Murray Region of South Australia:Exploration Geophysics, 40, 173–183. → pages 78127Vozoff, K., and D. L. B. Jupp, 1975, Joint inversion of geophysical data:Geophysical Journal of the Royal Astronomical Society, 42, 977–991. → pages11Wahlberg, B., S. Boyd, M. Annergren, and Y. Wang, 2012, An ADMM Algorithmfor a Class of Total Problems: arXiv preprint arXiv:1203.1828, 1–6. → pages10, 114, 120Ward, S. H., and G. W. Hohmann, 1988, Electromagnetic theory for geophysicalapplications: Electromagnetic methods in applied geophysics, 1, 131–311. →pages 20Weiss, C. J., and S. Constable, 2006, Mapping thin resistors and hydrocarbonswith marine EM methods , Part II Modeling and analysis in 3D: Geophysics,71, G321–G332. → pages 7Williams, N., D. Oldenburg, and P. Lelie`vre, 2009, Constraining gravity andmagnetics inversions for mineral exploration using limited geological data:ASEG Expanded Abstracts. → pages 11Wilson, G. A., A. Raiche, and F. Sugeng, 2006, 2.5D inversion of airborneelectromagnetic data: Exploration Geophysics, 37, 363–371. → pages 9Xiong, Z., and A. C. Tripp, 1995, A block iterative algorithm for 3Delectromagnetic modeling using integral equations with symmetrizedsubstructures: Geophysics, 60, 291–295. → pages 11, 44Yang, D., D. W. Oldenburg, and E. Haber, 2014, 3-D inversion of airborneelectromagnetic data parallelized and accelerated by local mesh and adaptivesoundings.: Geophysical Journal International, 196.3, 1492–1507. → pages 9Zhdanov, M. S., 2010, Electromagnetic geophysics : Notes from the past and theroad ahead: Geophysics, 75, 75A49–75A66. → pages 7Zhdanov, M. S., and L. H. Cox, 2013, Multinary Inversion for Tunnel Detection:IEEE Geoscience and Remote Sensing Letters, 10, 1100–1103. → pages 44Zheglova, P., and C. Farquharson, 2016, Joint level set inversion of gravity andtravel time data : application to mineral exploration: 86th Annual InternationalMeeting, SEG, Expanded Abstracts, 35, 2165–2169. → pages 44, 60Zonge, K. L., and L. J. Hughes, 1991, Controlled source audio-frequencymagnetotellurics: Electromagnetic methods in applied geophysics, 2, 713–809.→ pages 4, 20128Appendix AMultiple Sphere ParametricInversionA.1 IntroductionHere I introduce a multiple sphere parameterization for AEM data as a way tocomplement the Gaussian ellipsoid approach discussed in this dissertation. This ismeant to demonstrate the flexibility of the developed parametric framework and toemphasize how various shape choices are available. The reader is also encouragedto see Aghasi et al. (2011); Pidlisecky et al. (2011) for discussions about othershape reconstructions.A.2 Sphere ParameterizationA sphere is a simplified shape compared to a Gaussian ellipsoid, because althoughthe radius and central location of the sphere can change, it can not stretch or rotate.This means that fewer parameters are needed to describe a sphere compared to anellipsoid, with the obvious drawback of a reduction in shape flexibility.Mathematically, the transition in the parametric algorithm from multiple Gaus-sian ellipsoids to multiple spheres is straight forward. For positions x,y,z in thespatial domain, define a position vector x and a central location vector x0(n) for the129nth anomaly as in Equation 4.2x =xyz x0(n) =x0(n)y0(n)z0(n) . (A.1)The parametric function τn is once again defined asτn = c− (x−x0(n))TMn(x−x0(n)) (A.2)but the former skewing matrix Mn is now a diagonal matrix with the inverse of thesquared sphere radius, r, on the diagonalMn =1r2 0 00 1r2 00 0 1r2 . (A.3)Each sphere can be given its own radius value, but in this example the radii are fixedat r. The resistivity of each sphere, ρ1, is also fixed with the idea that anomaliescan arrange themselves in ways to either increase or decrease the desired resistivityeffect. In this manner, the inversion is designed to be an iterative process that givesthe user the ability to test different scenarios. Finally the resistivity in each meshcell is defined asρ(τn,ρ0,ρ1) = ρ0+12(1+ tanh(aτn))(ρ1−ρ0) (A.4)and the derivatives for the Jacobian matrix are calculated in the same manner asEquation 3.23, without the need for derivatives with respect to skewing parame-ters which no longer exist. To compensate for the rigid form of the sphere, themesh is populated with many spheres, i.e. 100 or more, allowing the inversion theflexibility to fit non-linear target features not well suited for a Gaussian ellipsoidparameterization. With a large number of spheres in the inversion, the regulariza-tion trade-off term β is set to unity and follows a cooling schedule as described inChapter 2 with a cooling factor γ = 0.9. This ensures there exists a penalty term fordeviating from a reference model, which prevents the spheres from scattering eas-130ily throughout the model. A much more conservative cooling schedule is neededcompared to previous chapters, otherwise the effect of the regularization is too mildand spheres begin to spread out too quickly.A.3 Synthetic ExampleTo test the multiple sphere approach, synthetic AEM data from Chapter 4 are re-visited. In this example, frequency-domain data are used, but similar results areexpected for time-domain measurements. The data and inversion parameters arekept the same as Chapter 4, but instead of three Gaussian ellipsoids, the startingguess is composed of 100 spheres of 75 m radius with a fixed resistivity of 1 Ωm.Three cluster centers are spatially chosen from the observed synthetic data, and theinitial spheres are split up and populated randomly around these points with a stan-dard deviation of 200 m, 200 m, 50 m in x,y,z respectively. The depth of the clustercenters is chosen to be 60 m, or z = 250 m. This keeps the spheres near the surfaceand ensures that any depth or dip information from the ensuing inversion comesfrom the data and not from the initial model. The left panel of Figure A.1a showsa plan slice through the true model at z = 250 m, where cluster centers are shownas black crosses and the location of two East-West cross-sections at y = 7332880m and y = 7332640 m are displayed in white. These cross-sections through thesynthetic model are shown in the right panel of Figure A.1a. Figure A.1b presentsa plan view and two cross-sections through the initial guess. The exact nature ofthe initial guess is clearly user-dependent and will vary, but this example is meantto demonstrate a typical guess when little a priori information is available.A.4 Synthetic ResultFigure A.1c portrays the synthetic multiple sphere parametric inversion results inplan view and cross-section. In plan view, the spheres are able to reconstruct thetrue model quite well, and even in cross-section there is some sense of the truedip in all three anomalies. The initial data misfit is 6.3e5, and after 74 Gauss-Newton iterations the final data misfit levels off at 38.1. Recall the data misfitfor the Gaussian ellipsoid example was 19.6, so the multiple spheres method doesnot achieve a similar level of data fit, but it does provide valuable information. If131the anomalies were curved or non-linear the sphere approach would have a muchgreater advantage, and is worthy of future research.y(m)x(m)y(m)x(m)x(m)z(m)a)c)ΩmΩmz(m)x(m)x(m)z(m)z(m)x(m)True ModelSphere modely(m)x(m)x(m)z(m)b)Ωmz(m)x(m) Initial modelxxxFigure A.1: Multiple sphere parametric inversion, frequency-domain AEMdata. Plan view at z = 250 m (left) and cross-section across solid whitelines at y = 7332880 m and y = 7332640 m (right). a) True model withcluster centers marked with black x’s. b) Initial guess. c) Parametricsphere inversion.132Appendix BAirborne Electromagnetic DataSamplingB.1 Total Horizontal Gradient Data SamplingThe survey locations in Chapter 4 are selected through a data reduction samplingmethod based on the total horizontal gradient (THG) magnitude of the data, definedfor time-domain ∂bz∂ t measurements asTHG =√√√√(∂ ( ∂bz∂ t )∂x)2+(∂ ( ∂bz∂ t )∂y)2. (B.1)Figure B.1 visually demonstrates the THG sampling method, where Figure B.1ashows the synthetic time-domain AEM ∂bz∂ t data from chapter 4 at 150 µs with thefull data set of 1172 source locations presented as black dots. Figure B.1b thenillustrates a standard evenly down-sampled ∂bz∂ t data set at 150 µs where everysecond source location along-line is discarded. Figure B.1c depicts a THG down-sampled data set overlaid on the THG of ∂bz∂ t data at 150 µs. Here, all data pointsare kept in regions where the THG value remains above a user-defined threshold,0.06 in this case, whereas outside this core region only every fourth point is se-lected.133a)b)c)x (m)x (m)x (m)y (m)y (m)y (m)log10(V/Am2)log10(V/Am2)log10Figure B.1: Time-domain AEM data sampling techniques. Data locationsshown as black dots. a) Full synthetic ∂bz∂ t data set at 150 µs with 1172total source locations. b) Evenly down-sampled synthetic ∂bz∂ t data setat 150 µs with 533 total source locations. c) THG of ∂bz∂ t data at 150 µswith 522 selected source locations based from the THG.134The THG down-sampling technique focuses on keeping data points where thesignal is rapidly changing to retain the most valuable information from the data setwhile eliminating data where the signal is constant or slowly varying. This methodis clearly dependent on which time channel or frequency is chosen to calculate theTHG. In this example the first time-channel is selected because its THG image isrepresentative of the entire data set, but this may not be the case for every scenario.Synthetic parametric inversion results are displayed in Figure B.2 for each ofthe aforementioned data sampling methods. Inversion parameters and meshes arekept to match those from Chapter 4. The left panel shows a 250 m constant eleva-tion slice through the inversion models and the thick white East-West lines depicttwo cross-section locations at y = 7333080 m and y = 7332780 m that are displayedin the right panel. The true model conductors are outlined in thin solid white andthe location of the starting guesses are in thin dashed white. All methods recoveran accurate reconstruction of the true conductors in both plan and cross-section;however, the full and THG sampled data inversions have a better overall recoveryof the dip of the Western conductor and a closer match to the true resistivities.Recovered resistivities, total number of Gauss-Newton iterations and final datamisfits are displayed in Table B.1. This table shows that the THG sampling methodquantitatively has the best result of the three methods with a reconstruction accu-racy Rc of 1.52 when looking only at conductive cells and RAll = 0.17 using allcells in the true model. Note that the evenly down-sampled recovery has a worsereconstruction accuracy of Rc = 3.27, and also requires 269 Gauss-Newton itera-tions before terminating and achieves a final data misfit of 31.1. This is comparedto only 23 Gauss-Newton iterations and a final misfit of 21.4 for the THG dataset. This suggests that a THG down-sampling approach is preferable over evenlydown-sampling the data, and little information is lost compared to inverting the fulldata set. Although this method is presented with a time-domain example, similarresults occur for frequency measurements, but are not shown.135Table B.1: Synthetic inversions summary: down-sampling methods.Model or method ρ1 ρ2 Rc RAll # GN Final φdTrue model 0.20 Ωm 2.00 Ωm 0.0 0.0 - -Full data 0.14 Ωm 1.94 Ωm 2.01 0.21 22 15.8Even sampling 0.45 Ωm 1.91 Ωm 3.27 0.30 269 31.1THG sampling 0.16 Ωm 1.95 Ωm 1.52 0.17 23 21.4y(m)x(m)y(m)x(m)x(m)z(m)a)c)ΩmΩmz(m)x(m)x(m)z(m)z(m)x(m)Full samplingTHG samplingy(m)x(m)x(m)z(m)b)Ωmz(m)x(m) Even samplingFigure B.2: a) Full sampling. b) Even sampling. c) THG sampling.136
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Cooperative and parametric strategies for 3D electromagnetic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Cooperative and parametric strategies for 3D electromagnetic inversion McMillan, Michael S. G. 2017
pdf
Warning
You are currently on our download blacklist and unable to view media. You will be unbanned within an hour.
To un-ban yourself please visit the following link and solve the reCAPTCHA, we will then redirect you back here.
Page Metadata
Item Metadata
Title | Cooperative and parametric strategies for 3D electromagnetic inversion |
Creator |
McMillan, Michael S. G. |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | Airborne electromagnetic (AEM) data is commonly collected for detecting buried natural resources, and this technique is sensitive to subsurface electrical resistivity distributions. The subsequent process of 3D AEM inversion constructs a physical property model from this data in order to better understand the size and shape of potential hidden resources. This thesis is designed to develop practical strategies to improve 3D AEM inversion accuracy in geologic settings where AEM data sets produce inconsistent or unsatisfactory inversion results. In this research, two overarching problematic scenarios are examined. First, in regions where an AEM survey overlaps with other electromagnetic data sets, a novel cooperative approach is introduced. This method is first tested on synthetic data where instead of producing an inversion model from each data set, the cooperative algorithm finds one consistent 3D resistivity model with improved resolution. The approach is then applied to field data over a high-sulfidation epithermal gold deposit where similar improvements are observed. The second scenario relates to improving 3D AEM inversions over thin conductive anomalies, a common geophysical target for copper and gold deposits. A new parametric inversion is developed using skewed Gaussian ellipsoids to represent target bodies. The approach is general but applied to frequency and time-domain AEM data with one or multiple anomalies. Combined with a voxel algorithm, the parametric inversion forms a hybrid approach capable of resolving thin conductive targets with smooth surrounding features. This hybrid technique is tested on synthetic data over conductive targets in a resistive background, and consistently produces models that are easier to interpret compared to voxel inversions alone. Field examples from a volcanogenic massive sulfide and an orogenic gold deposit highlight the practical nature of this method to image conductive mineralization with increased precision. The thesis concludes by analyzing a setting where multiple spatially overlapping AEM data sets exist over narrow conductive anomalies. Here, parametric, cooperative and voxel inversions are applied together to generate a consistent 3D resistivity model with thin targets and smooth background features. This section includes a discussion about potential pitfalls of such an approach when incompatible field measurements are encountered. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-04-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NoDerivatives 4.0 International |
DOI | 10.14288/1.0343483 |
URI | http://hdl.handle.net/2429/61125 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2017-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_may_mcmillan_michael.pdf [ 49MB ]
- Metadata
- JSON: 24-1.0343483.json
- JSON-LD: 24-1.0343483-ld.json
- RDF/XML (Pretty): 24-1.0343483-rdf.xml
- RDF/JSON: 24-1.0343483-rdf.json
- Turtle: 24-1.0343483-turtle.txt
- N-Triples: 24-1.0343483-rdf-ntriples.txt
- Original Record: 24-1.0343483-source.json
- Full Text
- 24-1.0343483-fulltext.txt
- Citation
- 24-1.0343483.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0343483/manifest