"Science, Faculty of"@en . "Earth, Ocean and Atmospheric Sciences, Department of"@en . "DSpace"@en . "UBCV"@en . "Holtham, Elliot Mark"@en . "2012-06-29T18:03:57Z"@en . "2012"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "Z-Axis Tipper Electromagnetic (ZTEM) data are airborne electromagnetic data which relate the vertical magnetic field measured by a helicopter to the orthogonal horizontal magnetic field measured at a ground-based reference station. The 30-720 Hz frequency range makes it possible to image structures at a kilometer depth for many geologic settings. This depth of penetration, coupled with rapid spatial acquisition of an airborne system, means that ZTEM data can be used to map large-scale\nstructures that are difficult to survey with ground-based surveys. I present some fundamentals of the ZTEM signatures and then develop a Gauss-Newton inversion algorithm which I test on a synthetic and a field example from the Mt.\nMilligan porphyry deposit.\nFor large airborne surveys, the high number of cells required to discretize the area at a reasonable resolution can make the computational cost of inverting the dataset prohibitively expensive. I present an iterative methodology that can be used to invert large natural source surveys by using a combination of coarse and fine meshes and a domain decomposition that splits the full model into smaller subproblems which can be run in parallel. After each round of tiled inversions, the\ntiles are merged together to form an updated model. I demonstrate this procedure first by inverting the data computed from a large synthetic model, before working with a survey over the Pebble porphyry deposit. The inverted ZTEM results are consistent with other electromagnetic datasets.\nSince a 1D conductivity structure produces no vertical magnetic fields, the technique cannot recover information about a purely 1D model. In order to increase sensitivity to the background structure, and greatly improve the depth of investigation over just ZTEM data, magnetotelluric (MT) and ZTEM data can both be collected. The combination of sparse MT data, with the economical and rapid spatial acquisition of airborne ZTEM data, creates a cost effective exploration technique that can map structures at depths. I develop an algorithm to jointly invert ZTEM and MT data, and apply this technique to a synthetic model and a field example from the Reese River geothermal property."@en . "https://circle.library.ubc.ca/rest/handle/2429/42582?expand=metadata"@en . "3D Inversion of Natural Source Electromagnetic Data by Elliot Mark Holtham BASc, The University of British Columbia, 2007 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Geophysics) The University Of British Columbia (Vancouver) June 2012 c\u00C2\u00A9 Elliot Mark Holtham, 2012 Abstract Z-Axis Tipper Electromagnetic (ZTEM) data are airborne electromagnetic data which relate the vertical magnetic field measured by a helicopter to the orthogonal horizontal magnetic field measured at a ground-based reference station. The 30- 720 Hz frequency range makes it possible to image structures at a kilometer depth for many geologic settings. This depth of penetration, coupled with rapid spatial acquisition of an airborne system, means that ZTEM data can be used to map large- scale structures that are difficult to survey with ground-based surveys. I present some fundamentals of the ZTEM signatures and then develop a Gauss-Newton inversion algorithm which I test on a synthetic and a field example from the Mt. Milligan porphyry deposit. For large airborne surveys, the high number of cells required to discretize the area at a reasonable resolution can make the computational cost of inverting the dataset prohibitively expensive. I present an iterative methodology that can be used to invert large natural source surveys by using a combination of coarse and fine meshes and a domain decomposition that splits the full model into smaller subproblems which can be run in parallel. After each round of tiled inversions, the tiles are merged together to form an updated model. I demonstrate this procedure ii first by inverting the data computed from a large synthetic model, before working with a survey over the Pebble porphyry deposit. The inverted ZTEM results are consistent with other electromagnetic datasets. Since a 1D conductivity structure produces no vertical magnetic fields, the technique cannot recover information about a purely 1D model. In order to in- crease sensitivity to the background structure, and greatly improve the depth of investigation over just ZTEM data, magnetotelluric (MT) and ZTEM data can both be collected. The combination of sparse MT data, with the economical and rapid spatial acquisition of airborne ZTEM data, creates a cost effective exploration tech- nique that can map structures at depths. I develop an algorithm to jointly invert ZTEM and MT data, and apply this technique to a synthetic model and a field example from the Reese River geothermal property. iii Preface Chapters 2,3 and 4 of this thesis are versions of research papers published or sub- mitted to peer-reviewed scientific journals. The publication details of each chapter, as well as the contributions of each co-author are outlined below Chapter 2 - 3D Inversion of ZTEM Data Authors: Elliot Holtham and Doug Oldenburg A version of Chapter 2 has been published. Holtham, E. and D. Oldenburg, 2010, Three-dimensional inversion of ZTEM data: Geophysical Journal International, 182, 168-182. I completed the research and prepared the manuscript. Doug Old- enburg provided guidance on the direction of the research, and helped with editing and revising the manuscript. Chapter 3 - Large-scale inversion of ZTEM data Authors: Elliot Holtham and Doug Oldenburg This chapter has been accepted for publication in Holtham, E., and D. W. Olden- burg, 2012, Large-scale inversion of ztem data: Geophysics, 77,1-9. I completed iv the research and prepared the manuscript. Doug Oldenburg provided guidance on the direction of the research, and helped with editing and revising the manuscript. Chapter 4 - Joint Inversion of ZTEM and MT Data Authors: Elliot Holtham and Doug Oldenburg A version of this chapter has been submitted in 2012 for publication and this chapter has also been published in Holtham, E. and D. Oldenburg, 2010, Three- dimensional inversion of MT and ZTEM data, SEG Extended Abstracts, Denver. I completed the research and prepared the manuscript. Doug Oldenburg provided guidance on the direction of the research, and helped with editing and revising the manuscript. v Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Water Resources . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.3 Mineral Exploration . . . . . . . . . . . . . . . . . . . . 6 1.2.4 Hydrocarbon Exploration . . . . . . . . . . . . . . . . . . 9 1.3 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.1 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . 10 2 3D Inversion of ZTEM Data . . . . . . . . . . . . . . . . . . . . . . 13 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.1 ZTEM System and Data Acquisition . . . . . . . . . . . . 16 2.1.2 Magnetic Transfer Functions . . . . . . . . . . . . . . . . 17 2.1.3 Forward Modelling . . . . . . . . . . . . . . . . . . . . . 20 2.1.4 Currents and Transfer Functions . . . . . . . . . . . . . . 24 2.2 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3 Synthetic Inversion . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.1 Effect of the Reference Station . . . . . . . . . . . . . . . 33 2.4 3D Inversion of ZTEM Field Data . . . . . . . . . . . . . . . . . 38 2.4.1 Discretizing the Earth . . . . . . . . . . . . . . . . . . . 43 vi 2.4.2 Determining a Background Model . . . . . . . . . . . . . 45 2.4.3 Assigning Uncertainties . . . . . . . . . . . . . . . . . . 51 2.4.4 Mt. Milligan Inversion Results . . . . . . . . . . . . . . . 54 3 Large-scale Inversion of ZTEM Data . . . . . . . . . . . . . . . . . 59 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2 Inversion Methodology . . . . . . . . . . . . . . . . . . . . . . . 62 3.3 Synthetic Example . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.3.1 Preliminary Coarse Inversion . . . . . . . . . . . . . . . . 70 3.3.2 Interpolating Models and Computing Primary Fields . . . 73 3.3.3 Setting Up the Subdomains . . . . . . . . . . . . . . . . . 74 3.3.4 Perform Subdomain Inversions . . . . . . . . . . . . . . . 76 3.3.5 Reference Station Considerations . . . . . . . . . . . . . 78 3.3.6 Model Updates and Merging Tiled Inversions . . . . . . . 79 3.4 Pebble Field Data Example . . . . . . . . . . . . . . . . . . . . . 81 4 Joint Inversion of ZTEM and MT Data . . . . . . . . . . . . . . . . 96 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.2.1 Determining a Background Model . . . . . . . . . . . . . 103 4.2.2 Synthetic Example . . . . . . . . . . . . . . . . . . . . . 103 4.3 Field Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.3.1 Reese River Background . . . . . . . . . . . . . . . . . . 108 4.3.2 ZTEM and MT Inversion Results . . . . . . . . . . . . . 111 4.3.3 Interpretations . . . . . . . . . . . . . . . . . . . . . . . 116 5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.1 Research Summary . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.1.1 Inversion of ZTEM Data . . . . . . . . . . . . . . . . . . 122 5.1.2 Inversion of Large-scale Datasets . . . . . . . . . . . . . 123 5.1.3 Joint ZTEM and MT Inversion . . . . . . . . . . . . . . . 123 5.1.4 Future Research Opportunities . . . . . . . . . . . . . . . 125 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 vii List of Figures Figure 1.1 General correlations between common rock types and their ex- pected conductivities. ( c\u00C2\u00A9Society of Exploration Geophysi- cists, 1987, adapted by permission from Palacky (1987)) . . . 4 Figure 2.1 ZTEM system components. ( c\u00C2\u00A9Geophysical Journal Interna- tional, 2010, adapted by permission from Holtham and Olden- burg (2010)) . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Figure 2.2 Electromagnetic skin depths for 3 ZTEM frequencies com- puted for typical geologic conductivity values. . . . . . . . . 19 Figure 2.3 Regular rectangular mesh structure used throughout the thesis. 21 Figure 2.4 a) 1 \u00E2\u0084\u00A6m conducting block buried in a 200 \u00E2\u0084\u00A6m background. b) 10000 \u00E2\u0084\u00A6m resistive block buried in a 200 \u00E2\u0084\u00A6m background. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by per- mission from Holtham and Oldenburg (2010)) . . . . . . . . . 25 Figure 2.5 The real anomalous current densities (Log10( Am2 )) in the z-y plane (x= -25 m) for a source whose electric field is polarized in the (y\u00CC\u0082) direction. The white rectangles are the outlines of the original blocks. Figure a) is the current for the conduc- tive block, while figure b) is the current for the resistive block. The anomalous current is mostly localized inside the block. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by per- mission from Holtham and Oldenburg (2010)) . . . . . . . . . 26 Figure 2.6 Resulting transfer functions for the two block example of fig- ure 2.4. Figure a) corresponds to a conductive block while fig- ure b) corresponds to a resistive block. The transfer functions closely resemble the magnetic fields due to a finite line current inside the blocks. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . 27 viii Figure 2.7 Conductivity model: a resistor surrounded by a conductor, all buried in a 100 \u00E2\u0084\u00A6m background halfspace. (a) Cross section through y = 0 m. The bottom table summarizes the block di- mensions and resistivities. (b) Depth slice at 1000 m. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . . . . . . . . . . . . . . . . . 32 Figure 2.8 Conductivity model (S/m) from the inversion of the synthetic model data. Figure a) is a cross section at the center of the model and a planview section at a depth of 800 m. b-d) are results of inverting data at 1 Hz. b) is the inversion of the real parts, c) is the imaginary parts, and d) is the real and imaginary components. e) is the simultaneous inversion of 6 frequencies. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by per- mission from Holtham and Oldenburg (2010)) . . . . . . . . . 34 Figure 2.9 Convergence curve for the inversion of six frequencies for the synthetic inversion model. The algorithm reduced the misfit to 1.38\u00C3\u0097 106 (target misfit 1.21\u00C3\u0097 106). ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . . . . . . . . . . . . . . . . . . . 35 Figure 2.10 Figure a) is the observed data and figure b) is the predicted data for the synthetic model at 1 Hz. There is excellent agreement between the observed and predicted data. . . . . . . . . . . . 36 Figure 2.11 Conductivity model (S/m) from the inversion of the synthetic model using 6 frequencies. Slices of the true model (left hand figures) and recovered model (right hand figures) are shown at 500, 1000, 1500m. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Figure 2.12 Percent difference in the x\u00CC\u0082 component of the magnetic field between the reference location and the measurement location (source electric field in the y\u00CC\u0082 direction). ( c\u00C2\u00A9Geophysical Jour- nal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . . . . . . . . . . . . . . . . . 38 ix Figure 2.13 Conductivity model (S/m) from the inversion of the synthetic model data using different reference station locations. For each inversion all of the parameters except the location of the ref- erence station remain unchanged. Figure a) corresponds to a reference station located at x = -3000 m, y = -3000 m. Figure b) corresponds to a reference station located at x = -900 m, y = 0 m. Figure c) corresponds to data computed using a reference station at x = -3000 m, y = -3000 m, but during the inversion process the reference station is mislocated at x = -900 m, y = 0 m. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . . . 39 Figure 2.14 Location of Mt. Milligan after Devine (2011b) in north cen- tral British Columbia. The deposit is located within the Ques- nel Terrane, a well known porphyry belt that extends through much of British Columbia. . . . . . . . . . . . . . . . . . . . 42 Figure 2.15 Regional geology around Mt. Milligan after Devine (2011b) . 43 Figure 2.16 Geologic cartoon of a typical alkalic porphyry deposit from Devine (2011a) . . . . . . . . . . . . . . . . . . . . . . . . . 44 Figure 2.17 Planview section after Mitchinson and Enkin (2011) detailing the main geologic features in the area. The main features of interest are the monzonite intrusives as well as the tertiary sed- iments to the east of the deposit. . . . . . . . . . . . . . . . . 45 Figure 2.18 Topography over the Mt. Milligan test ZTEM block. The to- pography over the survey area ranged from 934 m to 1475 m for a total of 541 m of relief. The flight lines are shown in black and the outlines of the main monzonite intrusives are outlined in red to get a sense of scale from the deposit scale geologic map in figure 2.17. . . . . . . . . . . . . . . . . . . . . . . . 46 Figure 2.19 Work flow process summarizing the critical step in a ZTEM in- version. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . 47 Figure 2.20 a) Conductivity (S/m) of the synthetic true model at a depth of 800 m. b) Misfit curve. c) Inverted model at 800 m, using a 10000 \u00E2\u0084\u00A6m starting model, taken at the sharp change of slope on the misfit curve. d) Final inversion model starting with a 10000 \u00E2\u0084\u00A6m initial model. The recovered model structure and conductivity values are not reasonable. Here \u00CE\u00B1s = 0, \u00CE\u00B1x =\u00CE\u00B1y = \u00CE\u00B1z = 1. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . 48 x Figure 2.21 Data misfit for various conductivity scaling factors. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . . . . . . . . . . . . . . . . . . . . . 50 Figure 2.22 a) Inverted model(S/m) at a depth of 800 m starting with the correct initial background conductivity. b) Inverted model at a depth of 800 m using the scaling method to determine the ini- tial conductivity. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) . 50 Figure 2.23 Synthetic test example of the Mt. Milligan topography to show that the method of fitting a best fitting halfspace to the data may be used if the topography response is much greater than the ge- ologic response. The synthetic model has a constant resistivity of 250 \u00E2\u0084\u00A6m. . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Figure 2.24 Misfit between the true data and the data produced from dif- ferent conductivity background models. As the conductivity values get further from the true conductivity, the data misfit increases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Figure 2.25 Observed data over Mt. Milligan at 180 Hz. Real part of Tzx, imaginary part of Tzx, real part of Tzy, and the imaginary part of Tzy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Figure 2.26 Predicted data over Mt. Milligan at 180 Hz. Real part of Tzx, imaginary part of Tzx, real part of Tzy, and the imaginary part of Tzy. The inversion fit the data well. . . . . . . . . . . . . . 55 Figure 2.27 Normalized misfits (normalized by the standard deviations of the data) over Mt. Milligan at 180 Hz. Real part of Tzx, imag- inary part of Tzx, real part of Tzy, and the imaginary part of Tzy. The misfits generally have little correlation with the ob- served data. . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Figure 2.28 Depth slice of the VTEM inversion results over the Mt. Mil- ligan deposit from Yang and Oldenburg (2012). (a) Slice at an elevation of 1030 m overlain by geology and (b) 420 \u00E2\u0084\u00A6m contour of DC resistivity model. The MBX stock is moder- ately resistive in both the 3D VTEM and DC resistivity inver- sions. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Yang and Oldenburg (2012)) . . . . . . . 57 xi Figure 2.29 Depth slice of the ZTEM inversion results over the Mt. Milli- gan deposit at an elevation of 1028m, the closest depth section to the 1030 m slice in figure 2.28. Both the 3D VTEM and ZTEM inversions recover a moderately resistive MBX stock. Also the conductive region north of the MBX stock is well mapped in both inversions. . . . . . . . . . . . . . . . . . . . 58 Figure 2.30 Depth slice of the ZTEM inversion results over the Mt. Milli- gan deposit at a depth of 728m. The conductive Great Eastern Fault zone and Tertiary Sediments are still well mapped at this depth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 Figure 3.1 Example of an L-shaped domain split up into 3 smaller do- mains with two shared boundaries after Saad (2003). . . . . . 64 Figure 3.2 Depth slice at -275 m of the true model. The topography for the synthetic model ranges from -25 m to -250 m. The conduc- tivity of the ore was chosen to be 0.2 S/m; however, the color scale on the model has been clipped at 0.01 S/m to improve the visualization of the other geologic units. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . 71 Figure 3.3 Depth slice at -475 m of the true model. The conductivity of the ore was chosen to be 0.2 S/m; however, the color scale on the model has been clipped at 0.01 S/m to improve the visual- ization of the other geologic units. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . . . . . . . . 72 Figure 3.4 Coarse inversion result at -475 m shown on the same color scale as the true model. Even with coarse 500 m x 500 m cells in the x and y directions, the main geologic features of the model have been recovered. ( c\u00C2\u00A9Society of Exploration Geo- physicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . . . . . . . . 73 Figure 3.5 Mesh layouts for the example using three overlapping domains. Each subdomain mesh contained approximately 3.5 million cells, with 353 x 140 x 70 cells in the easting, northing and vertical directions. The central core region of each domain overlapped by 20 cells, each with dimensions of 50 m, in both the x and y directions. ( c\u00C2\u00A9Society of Exploration Geophysi- cists, 2012, adapted by permission from Holtham and Olden- burg (2012)) . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 xii Figure 3.6 Mesh layouts for the example using nine overlapping domains. Each mesh contained approximately 1.4 million cells with 139 x 140 x 70 cells in the easting, northing and vertical direc- tions. The central core region of each domain overlapped by 20 cells, each with dimensions of 50 m, in both the easting and northing directions. ( c\u00C2\u00A9Society of Exploration Geophysi- cists, 2012, adapted by permission from Holtham and Olden- burg (2012)) . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Figure 3.7 Depth slice at -275 m of the inversion result after merging the three tiles together, shown on the same color scale as the true model. The resolution in this model is greatly improved com- pared to the coarse result. Several of the larger known deposits have now been recovered; they were not visible from the initial coarse inversion. The recovered mineralized body in the inver- sion between easting -4250 m and 0 m is a mineralized body in the synthetic model that only extends to -250 m. The ex- tend of this body that is not apparent in the -275 m depth slice of the true model, can be seen in figure 3.11. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . 82 Figure 3.8 Depth slice at -275 m of the inversion result after merging the nine tiles together, shown on the same color scale as the true model. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . 83 Figure 3.9 Depth slice at -475 m of the inversion result after merging the three tiles together, shown on the same color scale as the true model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Figure 3.10 Depth slice at -475 m of the inversion result after merging the nine tiles together, shown on the same color scale as the true model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Figure 3.11 Comparison of the larger southern mineralized regions. a) Re- covered bodies more conductive than 0.005 S/m that corre- spond to the mineralized bodies in the southern section of the inverted model. b) Actual mineralized bodies. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . 86 Figure 3.12 Location of the Pebble deposit in the Bristol Bay region of Alaska. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . 87 xiii Figure 3.13 Regional geology after Lang et al. (2009). The Pebble deposit sits in the Lake Clark Fault Zone. . . . . . . . . . . . . . . . 88 Figure 3.14 Map of the local geology and timeframe of the major geologic events near the Pebble deposit after Lang et al. (2009) . . . . 88 Figure 3.15 Plan view of the Pebble Deposit showing the deposit outline and grade distribution across the deposit. ( c\u00C2\u00A9Society of Ex- ploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . 89 Figure 3.16 Cross section of the Pebble Deposit showing the grade dis- tributions and the intrusive stocks. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . . . . . . . . 90 Figure 3.17 6504 line-km ZTEM survey block flown in the summer of 2010. The 30 km x 30 km block inverted in this section is outlined by the dashed square. The full survey block is shown as the solid black region. ( c\u00C2\u00A9Society of Exploration Geophysi- cists, 2012, adapted by permission from Holtham and Olden- burg (2012)) . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Figure 3.18 Coarse ZTEM inversion result at a depth of 100 m. The edge of the data coverage is shown by the dashed lines and the de- posit is outlined by the solid lines. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . . . . . . . . 92 Figure 3.19 Fine ZTEM inversion result at an elevation of 100 m formed by merging the three subproblem tiles together. The edge of the data coverage is shown by the dashed lines and the de- posit is outlined by the solid lines. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . . . . . . . . . . . . . . . . . . . 93 Figure 3.20 Spectrem Z component time constant plot over the same region as the inverted ZTEM data. The outline of the deposit is out- lined in black. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . 94 Figure 3.21 Inversion results shown on the same color scale and elevation of 100 m. The deposit is outlined in black. a) 3D DC Inver- sion b) ZTEM inversion shown over the same area as the DC resistivity data coverage. To first order the results seem consis- tent. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) . . . . . 95 xiv Figure 4.1 Typical ground MT site from Smirnov et al. (2008). Two sets of electrodes measure the electric field in two orthogonal direc- tions. The two horizontal magnetic fields are measured (possi- bly the vertical component as well). Because the MT site may be left for a long period of time, additional devices such as power sources, data loggers etc. may also be required. . . . . 99 Figure 4.2 Seafloor MT site. Each site measures the electric and magnetic fields, and after sufficient data has been recorded the MT site is released from the bottom of the ocean and collected by the control ship. ( c\u00C2\u00A9Society of Exploration Geophysicists, 1998, adapted by permission from Constable et al. (1998)) . . . . . 100 Figure 4.3 Cross section through northing = 0 m of the synthetic conduc- tivity structure of a resistor surrounded, except on the top and bottom, by a conductor. Both blocks are embedded in a 0.01 S/m background half-space. . . . . . . . . . . . . . . . . . . 104 Figure 4.4 Survey grids for the synthetic inversions. The three grids are the coarse MT (large squares), fine MT (small squares), and ZTEM (solid lines). . . . . . . . . . . . . . . . . . . . . . . . 106 Figure 4.5 Inversion of the synthetic model in figure 4.3, all shown on the true color scale. Adding the coarse MT data to the ZTEM data greatly improved the inversion result in the 2-3 km depth range. 107 Figure 4.6 Reese River geothermal site located in Nevada. On the east- ern portion of the survey area is the Range Front Mountains. The 7 MT lines are shown in red. The map of the location of the Reese River site in Nevada is provided curtesy of Sierra Geothermal Power Corp. . . . . . . . . . . . . . . . . . . . . 109 Figure 4.7 Reese River surface geology provided from Sierra Geothermal Power Corp. The 7 MT lines are shown in red along with the location of the 5 drill holes. A possible interpretation for the fluid flow upwelling and outflow is also outlined. The location of the Range Front Fault is outlined by the thick black line. . . 110 Figure 4.8 Reese River well logs. From Palacky (1987), one would expect a conductive structure in the near surface due to the conductive sediments, and a more resistive structure at depth because of the limestone and granite. . . . . . . . . . . . . . . . . . . . 111 Figure 4.9 ZTEM inversion going through drill holes 13-4 and 56-4. The inversion recovers a conductive structure on the west, due to the near surface sediments. The resistive structure east of the Range Front Fault is a result of the limestone unit shown in the geologic map in figure 4.7. . . . . . . . . . . . . . . . . . . . 113 xv Figure 4.10 MT inversion going through drill holes 13-4 and 56-4. The inversion recovers a conductive structure on the west, due to the near surface sediments. The resistive structure east of the Range Front Fault results from the limestone unit shown in the geologic map in figure 4.7. The MT data shows that the resistive structure extends to depths of greater than 1.5 km . . 114 Figure 4.11 Joint ZTEM and MT inversion going through drill holes 13- 4 and 56-4. The inversion recovers a conductive structure on the west, due to the near surface sediments. The complex re- sistive structure east of the Range Front Fault results from the limestone unit shown in the geologic map in figure 4.7. In the complex geologic region near the limestone unit, the joint in- version result recovers features that are a combination of the features in the two individual inversions. . . . . . . . . . . . . 115 Figure 4.12 Conceptual model reproduced from Cumming (2009) of a 150 \u00E2\u0097\u00A6C to 200 \u00E2\u0097\u00A6C geothermal reservoir with isotherms, alteration zones, and structures. . . . . . . . . . . . . . . . . . . . . . . . . . 119 Figure 4.13 Joint ZTEM/MT inversion showing the inverted model with a conductivity cutoff (only showing model cells less conductive than 0.05 S/m). The large resistive structure comes up from depth under the potential accommodation zone. The dip on the western side of the resistive structure seems to agree fairly well with the interpreted dip of the Range Front Fault. . . . . . . . 120 Figure 4.14 Inverted ZTEM and MT model across the prospective geother- mal upflow zone shown in figure 4.7. Superimposed on the model is a potential geothermal interpretation. . . . . . . . . . 120 xvi Acknowledgments Firstly I would like to thank the parties that have provided me with financial support throughout my graduate studies. In particular, the government of Canada NSERC and Vanier Canada Graduate Scholarship programs have provided me with first class funding for which I am extremely grateful. I would also like to thank the family and employees at Zonge Engineering, and the family of Stanley Ward for providing funding through the Society of Exploration Geophysics scholarship pro- grams. Also, I would like to thank Geoscience BC for their Geoscience BC schol- arship program. Next I would like to thank those individuals without which my thesis would not have been possible. Firstly, the person that had to put up with my endless requests, Roman Shektman, without whom nobody would graduate from UBC- GIF. I am extremely grateful to my supervisor Doug Oldenburg for mentoring me and providing countless opportunities over the years. Thank you for accepting and supporting my athletic pursuits. I don\u00E2\u0080\u0099t think anyone has more confidence in our abilities than you. To all the members of UBC-GIF thank you for your support and help over the years. To my family thank you for instilling in me the value of education, and to Megan thank you for supporting me through the ups and downs. Finally, I would like to thank the companies and individuals who have been kind enough to provide data and support for the various examples used in my the- sis. In the second chapter I would like to acknowledge Terrane Metals Corp. for providing support on the UBC-GIF field trip to Mt. Milligan. In the third chapter I would like to thank Sierra Geothermal Power and in particular Jeff Witter, for their help on the Reese River Project. In the final chapter I would like to thank Mike Allard from Xstrata Zinc, Louis Martin from Xstrata Copper and Gervais Peron from Mira Geosciences, for providing the geologic shape model for the Noranda mining camp that has formed the basis our synthetic example. Finally, I would like to thank Pascal Pare\u00CC\u0081 of Anglo American and The Pebble Partnership for providing the data and support on the pebble project. xvii Chapter 1 Introduction 1.1 Research Motivation Many of the earth\u00E2\u0080\u0099s remaining natural resources, such as aquifers, fossil fuels, and mineral deposits are buried deep in the subsurface of the earth and hidden from traditional surface investigation techniques. Demand for commodities such as oil, copper, iron ore and other manufacturing minerals remains high due to the rapidly developing economies of countries such as India and China (Bosworth and Collins, 2007; Nel and Cooper, 2008; Yellishettya et al., 2010). While the demand for these commodities remains high, as the remaining reserves are depleted, it becomes increasingly difficult to find replacements. In fact for most commodities, the easy to find deposits have already been discovered. Because of this it will be become more difficult and more expensive to find new resources (Blain, 2000). 1 Exploration programs on established targets are generally quite effective. For these targets, particularly when surface outcroppings exist, geologic mapping, geo- chemistry and drilling often play more prominent roles in the exploration process than geophysical methods. While following up on existing targets is relatively straightforward, finding new exploration targets and resources is accomplished through grassroot exploration on district and regional scales. On these scales, par- ticularly when the underlying geology is deeply hidden under cover, deep probing geophysical methods must be used to gather information. In fact, investigation of the earth\u00E2\u0080\u0099s interior, without direct contact such as digging or drilling, requires geophysical methods. While ground based geophysical methods can provide high quality datasets, when exploring over large areas, airborne methods (Palacky and West, 1987) must be employed. One available geophysical exploration tool is electromagnetic (EM) surveys that can be used to recover a 3D conductivity model of the earth. While the electrical conductivity itself is unlikely to be the final property of interest, by un- derstanding the link between the physical properties and lithology, alteration and other petrophysical properties of interest (see for example figure 1.1 that lists the expected conductivity of several common rock types), conductivity information can be used to map geologic and economic features of interest (Pellerin et al., 1996; Meju, 2002; Sheard et al., 2004). Not only can electromagnetic methods 2 be used to find resources, but they can also be used in countless other ways, such as finding unexploded munitions (Pasion et al., 2007), looking for water on Mars (Grimm, 2002), and environmental examples such as monitoring carbon sequestra- tion (Kirkendall and Roberts, 2004), and managing waste sites (Pellerin, 2002). Commonly EM investigations are carried out using local man-made controlled sources. The most significant drawback to this approach is the limited penetra- tion depth due to the geometric decay of the source fields. These drawbacks can be overcome by using planewave electromagnetic source fields such as lightning and solar events (Garcia and Jones, 2002; Simpson and Bahr, 2005). Using these plane wave natural source fields, which do not decay geometrically as they penetrate the earth, allows for exploration to greater depths, and for targets buried under conductive cover. While the traditional natural source method, the magnetotel- luric (MT) method (Tikhonov, 1950; Cagniard, 1953; Simpson and Bahr, 2005; Berdichevsky and Vladimir, 2008), has played a significant role in crustal studies (Eisel and Haak, 1999; Clark, 2000; Wei et al., 2001; Gokarn et al., 2001; Jones and Ferguson, 2001; Wannamaker et al., 2009) as well as in resource exploration (Constable et al., 1998; Hoversten et al., 1998; Mitsuhata et al., 1999; Farquharson and Craven, 2009), a practical limitation of this technique and other deep prob- ing controlled source electromagnetic methods is that surveys are costly and time consuming because they are extremely labor intensive. It would be preferable to 3 Figure 1.1: General correlations between common rock types and their ex- pected conductivities. ( c\u00C2\u00A9Society of Exploration Geophysicists, 1987, adapted by permission from Palacky (1987)) collect MT data in an aircraft but this goal has not yet been achieved because of the difficulty in measuring the electric fields in the air. While the first airborne nat- ural source system the AFMAG (Audio Frequency Magnetics) (Ward, 1959) was proposed and developed several decades ago, the concept had many problems and limitations outlined in Ward et al. (1966). Recently Geotech Inc., a geophysical 4 airborne contractor, has started flying a new natural source airborne electromag- netic (AEM) system, the ZTEM system (Lo and Zang, 2008). The current ZTEM system tows a coil from a helicopter that is used to measure the vertical component of the magnetic field over the survey area. The data depend on the relationship between this vertical magnetic field component at the location of the aircraft and the horizontal magnetic field recorded at a ground-based base station. In this thesis I examine this new type of data, and how it and MT data can be used to effectively image the earth at depth over large areas. 1.2 Applications Because of the deep penetration of natural source electromagnetic waves, the ZTEM and MT methods can be used for a variety of applications. The electrical conduc- tivity models obtained by inverting the geophysical data can be linked to properties of interest for the particular application. In this section a brief overview of some of the potential applications for this research is given while at the same time, I introduce the synthetic models and field data examples used in this thesis. 1.2.1 Water Resources Water content, particularly saline water, can affect the electrical properties of rocks (Herrick and Kennedy, 1994; Carcione et al., 2003; Lesmes and Friedman, 2005; Friedman, 2005). Underground aquifers that have contrasting electrical properties 5 from the surrounding materials can be mapped with electrical methods. Because underground water sources can extend over large regions, airborne electromagnetic methods, particularly deep probing methods such as the ZTEM method, can be used to map underground aquifers. 1.2.2 Energy Geothermal The joint inversion of ZTEM and MT data is an excellent technique to search for geothermal targets because it has the capability to penetrate the highly conductive environments often associated with geothermal reservoirs (Pellerin, 2002). One such geothermal resource that is examined in detail in this thesis is the Reese River property located in Nevada. When exploring for geothermal resources, electrical methods may be used to map the geologic structures with the hope that these struc- tures will give insight into structural controls for the fluid pathways. 1.2.3 Mineral Exploration Electrical conductivity can be used to map regional geology as well as localized lithology and alteration for mineral exploration projects. On the deposit scale, the inverted conductivity models can be used to extend known deposits or extrapolate between sparsely spaced drill holes. Known physical properties from drill holes can be used to constrain subsequent inversions and improve the inverted models. 6 Porphyry Systems Porphyry deposits are desirable exploration targets because of the potential large size and the polymetallic nature of the deposits. The possible gold, copper and molybdenum can provide economic diversity for the deposit during differing eco- nomic climates. In this thesis data over two porphyry deposits is examined. The Mt. Milligan deposit is an alkalic Cu-Au porphyry deposit located in north central British Columbia. In 2008 a ZTEM survey was flown for Geoscience BC and Terrane Metals Corp. over the deposit and surrounding areas. The test survey was flown with the hopes that the results from this survey could help when using airborne electromagnetic data to find additional Mt. Milligan style deposits in British Columbia. ZTEM data over this deposit will be examined in this thesis. The Pebble deposit is a world-class calc-alkalic copper-gold molybdenum por- phyry deposit located in the Bristol Bay region of southwest Alaska. Because of the size of the deposit, and the quality and diversity of geophysical data collected over the area, this deposit is an ideal test site for geophysical techniques and inver- sion methodologies. The Pebble Deposit consists of two zones, the Pebble West and East Zones. Because the Pebble West Zone extends near the surface, this part of the deposit was discovered in 1986, far before its deeper and higher grade neigh- bor, Pebble East, which was not discovered until 2005. Pebble East is completely hidden by up to 600 m of volcano-sedimentary cover which explains why it took al- 7 most two decades to discover the zone and highlights the importance in developing geophysical techniques that can image deep targets. Unconformity Uranium Canada hosts many unconformity type uranium deposits (Jefferson et al., 2007) lo- cated in the Athabasca Basin in Saskatchewan. The basic geologic structure is a re- sistive sandstone layer that extends to over 1 km in depth depending on the specific location in the Basin. Below the unconformity may lie resistive granitic gneiss or more conductive metasediments. These metasediments, especially those with high concentrations of graphite, can form strong conductors. Because these metasedi- ments are often spatially located near unconformity uranium deposits, these units are important imaging targets. Often these linear features can extend for several kilometers and can easily be mapped with EM data (in particular ZTEM data since the long conductors will create a good ZTEM response). Although the unconfor- mity may be buried deeper then 1 km, because of the resistive nature of the sand- stone, the ZTEM frequencies should easily penetrate beneath the unconformity. District Scale Exploration Natural source data, particularly large airborne ZTEM surveys, are excellent ex- ploration tools for promising mining districts. One Canadian mining district exam- ined in this thesis is the Noranda district which hosts 20 economic volcanogenic 8 massive sulphide deposits (VMS), 19 orogenic Au deposits, and several intrusion hosted Cu-Mo deposits (Gibson and Galley, 2007). 1.2.4 Hydrocarbon Exploration Land-based and marine-based magnetotelleric data have been used for hydrocar- bon exploration (Constable et al., 1998; Hoversten et al., 1998; Mitsuhata et al., 1999; Key et al., 2006; Zhdanov et al., 2011; Key, 2011). Because hydrocarbon resources are often deeply buried, natural source methods with their deep depth of penetration are excellent techniques to image the reservoirs. 1.3 Thesis Objectives In order for geophysical methods to be useful in answering questions about the subsurface of the earth, several critical conditions must be met. Firstly, the chosen geophysical method must be capable of imaging the desired physical properties at the depths of interest. Secondly, the collected geophysical measurements must be repeatable and of sufficiently high quality that the necessary information can be ex- tracted. Finally, robust interpretation schemes must exist to create 3D conductivity models of the subsurface such that they can be used for interpretation. In this thesis I try and answer the question of how to determine the electrical conductivity in the top few kilometers of the earth. I examine this problem by using natural source electromagnetic methods and in particular I address and attempt to answer: 9 \u00E2\u0080\u00A2 Can ZTEM data be inverted to recover useful information about the electrical properties of the earth ? \u00E2\u0080\u00A2 Do the models produced by inverting ZTEM data seem consistent with those from other electrical geophysical methods ? \u00E2\u0080\u00A2 Do the models produced by inverting ZTEM data seem consistent with the known geology of the region ? \u00E2\u0080\u00A2 Can I develop methodologies to handle large natural source datasets ? \u00E2\u0080\u00A2 Can adding MT data provide additional information for ZTEM datasets ? In this thesis I attempt to answer these questions by developing inversion al- gorithms and methodologies to invert both ZTEM and MT data. I apply these methods to a series of synthetic and field data examples. 1.3.1 Thesis Outline While the thesis consists of a series of chapters that have been published, or submit- ted into peer reviewed journals, the thesis flows naturally under the common theme of developing natural source inversion algorithms and methodologies. While the newly developed ZTEM system is a promising geophysical technique, in order that meaningful information be extracted from the data, interpretation methodolo- gies must exist. In chapter 3 of this thesis I explore this new data type in detail 10 and develop a 3D inversion algorithm for ZTEM data. I test the algorithm on some synthetic examples and from a field dataset from the Mt. Milligan porphyry deposit in British Columbia. The real strength of natural source datesets lies in the ability to image over large areas at depth. Both government and industry have come to recognize the importance of electromagnetic techniques in imaging regional and district scale structure (for instance the QUEST and QUEST WEST projects (Phillips et al., 2009) flown over much of central British Columbia, and a multitude of datasets flown by Geoscience Australia). While collecting large high quality datasets is a step in the right direction, in order to extract the most information from this data, one must have methodologies and workflow procedures to invert extremely large datasets. In chapter 3 I develop an efficient workflow methodology to invert large- scale natural source datasets. I test this methodology on a geologic model of the Noranda mining camp, as well as the Pebble porphyry deposit in Alaska. Since a 1D conductivity structure produces no vertical magnetic field, the ZTEM technique can have limited sensitivity to models that are nearly 1D in nature. In order to increase the sensitivity to layered structures, reduce the non-uniqueness of the inverse problem and greatly improve the depth of investigation, ZTEM and MT data can both be collected. In chapter 4 I demonstrate how the combination of sparse MT data, with the economical and rapid spatial acquisition of airborne 11 ZTEM data, creates a cost effective geophysical technique that can map large- scale structures at depths that are difficult to image with other techniques. I test this technique on both a synthetic example and a field geothermal example from Reese River in Nevada. Chapter 5, the concluding chapter of this thesis, summarizes the results from the previous chapters and lays out future paths of research. 12 Chapter 2 3D Inversion of ZTEM Data 2.1 Introduction Natural source electromagnetics play an important role in understanding the elec- trical conductivity of the upper regions of the earth. Their primary advantage, compared to controlled source methods, is the large depth of penetration that is a consequence of the plane wave excitation. The MT method, relates the electric and magnetic fields as data and it has played a significant role in crustal studies as well as in mining and hydrocarbon exploration. Its recognized importance, combined with increased computational abilities, has been the impetus for significant work in three-dimensional MT inversion (among others Mackie and Madden (1993), New- man and Alumbaugh (2000), Zhdanov et al. (2000), Farquharson et al. (2002), Sasaki (2004), Siripunvaraporn et al. (2005)). A practical limitation of the MT 13 technique and other deep probing controlled source electromagnetic methods, is that surveys are costly and time consuming because they are extremely labor inten- sive. It would be preferable to collect MT data in an aircraft but this goal has not yet been achieved because of the difficulty in measuring the electric fields in the air. In an effort to continue to use the penetration advantage of natural sources, it has long been recognized that tipper data, a measurement relating the local ver- tical magnetic field to the horizontal magnetic field, provide information about 3D electrical conductivity structures. The underlying reason is that the inducing electromagnetic fields are vertically propagating plane waves and if the earth is 1D then the vertical component of the magnetic field is zero. Non-zero values of the tipper data are thus directly related to anomalous currents due to non 1D conductivity structures. It was this understanding that prompted the development of AFMAG (Audio Frequency Magnetics) (Ward (1959)). The original airborne AFMAG technique used the amplitude outputs from two orthogonal coils towed behind an aircraft to determine the tilt of the plane of polarization of the natural magnetic field. The tilt angle is zero over a 1D earth and hence the technique was particularly effective when crossing conductors. However, because the direction of the inducing field varied with time, AFMAG results were not always repeatable. Limitations of the AFMAG technique were outlined in Ward et al. (1966). 14 Some of the AFMAG problems can be removed by using improved signal pro- cessing and data acquisition techniques. In particular Labson et al. (1985) devel- oped a technique that used ground-based horizontal and vertical coils to measure the magnetic fields. They used MT processing techniques to show how tipper data, a measurement relating the vertical to horizontal magnetic fields, could be obtained from the measured magnetic fields. A further improvement to the original AFMAG technique combines improved instrumentation and MT data processing techniques. This has resulted in the ZTEM technique (Lo and Zang (2008)). In ZTEM, the vertical component of the mag- netic field is recorded above the entire survey area, while the horizontal fields are recorded at a ground-based reference station. MT processing techniques yield fre- quency domain transfer functions that relate the vertical fields over the survey area to the horizontal fields at the reference station. Since new instrumentation exists to measure the vertical magnetic fields by helicopter, data over large survey areas can quickly be collected. The result is a cost effective procedure for collecting natural source EM data that provides information about the 3D conductivity structure of the earth. Over the last several years industry has recognized the potential in this technique and the need to be able to invert these data. In this chapter we present a three dimensional inversion algorithm for ZTEM data. Forward modelling is fundamental to any inversion algorithm and we first 15 show how to model the magnetic transfer functions. A simple example is used to illustrate the principles associated with ZTEM data. Next, we present our algorithm for inverting ZTEM data and apply it to a synthetic test problem. Finally, the technique is tested on a field dataset from the Mt. Milligan porphyry deposit in British Columbia. We discuss practical issues of implementing the inversion and develop a workflow methodology to make the procedure more efficient. 2.1.1 ZTEM System and Data Acquisition ZTEM instrumentation, developed by Geotech Inc, uses a helicopter to tow a coil which measures the vertical component of the magnetic field. Orthogonal coils from a single reference station measure the horizontal components of the magnetic field. The essential components of the system can be seen in figures 2.1a- 2.1c. Time series of the magnetic fields are recorded with a fixed sampling rate and data are binned and processed to generate transfer functions in the frequency domain. The lowest frequency of the transfer function depends upon the speed of the heli- copter, and the highest depends upon the sampling rate, signal strength and noise. Although frequencies up to 2800 Hz are potentially available, in practise the use- able bandwidth is about 22 - 720 Hz. The lower frequency has an electromagnetic skin depth (the depth at which the amplitude of the plane wave attenuates to 1e of the amplitude at the surface of the earth) of 3300 m in a 1000 \u00E2\u0084\u00A6m background and 330 m for a 10\u00E2\u0084\u00A6m background (skin depths for 3 ZTEM frequencies are computed 16 over a range of conductivity values and shown in figure 2.2). These are depth lim- its that are extremely difficult to achieve with controlled source ground or airborne systems. In addition to the collected ZTEM data, the system also collects aeromag- netic data that can provide additional information for geologic interpretation 2.1.2 Magnetic Transfer Functions ZTEM data are transfer functions that relate the vertical magnetic fields computed above the earth to the horizontal magnetic field at some fixed reference station. This relation is given by Hz(r) = Tzx(r,r0)Hx(r0)+Tzy(r,r0)Hy(r0), (2.1) where r is the location for the vertical field, r0 is the location of the ground based reference station, Tzx and Tzy are the vertical field transfer functions, Hz is the ver- tical magnetic field and Hx and Hy are the horizontal magnetic fields in the x and y directions. Solving for the transfer functions during the numerical modelling requires that the vertical fields are known for two independent polarizations. The fields for each polarization are given by H(1)z = TzxH (1) x +TzyH (1) y , (2.2a) 17 (a) ZTEM helicopter system components including the coil used to measure the vertical magnetic field. (b) Ground based reference station (c) Coil used to measure the ver- tical fields Figure 2.1: ZTEM system components. ( c\u00C2\u00A9Geophysical Journal Interna- tional, 2010, adapted by permission from Holtham and Oldenburg (2010)) 18 Figure 2.2: Electromagnetic skin depths for 3 ZTEM frequencies computed for typical geologic conductivity values. H(2)z = TzxH (2) x +TzyH (2) y , (2.2b) or \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD H (1) z (r) H(2)z (r) \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8= \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD H (1) x (r0) H (1) y (r0) H(2)x (r0) H (2) y (r0) \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD Tzx Tzy \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 , (2.3) where the superscripts (1) and (2) refer to two independent source field polariza- tions. Solving equation 2.3 yields the following expressions for the transfer func- 19 tions, Tzx = H(2)y H (1) z \u00E2\u0088\u0092H(1)y H(2)z H(1)x H (2) y \u00E2\u0088\u0092H(2)x H(1)y , (2.4a) Tzy = \u00E2\u0088\u0092H(2)x H(1)z +H(1)x H(2)z H(1)x H (2) y \u00E2\u0088\u0092H(2)x H(1)y . (2.4b) In order to forward model the transfer functions, one must have the capability to solve the natural source field problem for different polarizations of the source field. A procedure to do this is presented next. 2.1.3 Forward Modelling While the final goal is to invert a given dataset in order to recover a physical prop- erty model, before this can be accomplished the forward problem must understood and solved. The geophysical forward problem is given as, F(m)+ \u00CE\u00B5 = d (2.5) where F is the forward modelling operator, m the physical property model, \u00CE\u00B5 is a vector of measurement errors, and d is the data observations. In order to compute the forward response, F(m) for a model m to which simple analytic solutions do not exist, the earth must be discretized and the necessary equations solved. In this thesis the earth is discretized using regular rectangular meshes such as the one shown in figure 2.3 (although using less structured meshes is also possible, see for 20 Figure 2.3: Regular rectangular mesh structure used throughout the thesis. example Haber and Heldmann (2007) and Schwarzbach et al. (2011)). Our ZTEM forward modelling procedure is essentially that outlined in Far- quharson et al. (2002), where the solution of Maxwell\u00E2\u0080\u0099s equations is that of Haber et al. (2000b). Briefly, Maxwell\u00E2\u0080\u0099s equations in the quasi-static regime, when com- bined with the constitutive relations of charge conservation and Ohm\u00E2\u0080\u0099s law, form the necessary equations for ZTEM modelling. These equations read, \u00E2\u0088\u0087\u00C3\u0097E\u00E2\u0088\u0092 i\u00CF\u0089\u00C2\u00B5H= 0, (2.6a) \u00E2\u0088\u0087\u00C3\u0097H\u00E2\u0088\u0092\u00CF\u0083E= 0, (2.6b) \u00E2\u0088\u0087 \u00C2\u00B7J= 0, (2.6c) 21 J\u00E2\u0088\u0092\u00CF\u0083E= 0, (2.6d) where E is the electric field, H is the magnetic field, \u00CF\u0089 is the angular frequency, \u00CF\u0083 is the conductivity, and J is the current density. We have assumed an e\u00E2\u0088\u0092i\u00CF\u0089t time- dependence and that the magnetic permeability \u00C2\u00B5 is constant and equal to that of free space, \u00C2\u00B50. To solve for the fields, the electric field is decomposed into vector and scalar potentials E = A+\u00E2\u0088\u0087\u00CF\u0086 and the Coulomb gauge condition \u00E2\u0088\u0087 \u00C2\u00B7A = 0 is imposed for uniqueness. Using the decomposition and eliminating H yields the potential equations \u00E2\u0088\u00872A+ i\u00CF\u0089\u00C2\u00B50\u00CF\u0083(A+\u00E2\u0088\u0087\u00CF\u0086) = 0, (2.7) \u00E2\u0088\u0087 \u00C2\u00B7J= \u00E2\u0088\u0087 \u00C2\u00B7 (\u00CF\u0083(A+\u00E2\u0088\u0087\u00CF\u0086)) = 0. (2.8) For our modelling we discretize the earth into rectangular cells and, after applying a finite volume technique (Haber et al. (2000b)) to equations 2.7 and 2.8, we obtain the system of equations to be solved. In this staggered discretization, the scalar potential is defined on the cell centers while the vector potential is defined by the normal components at the center of the cell faces. To solve for the total fields we use a primary-secondary decomposition. The primary fields equations, which must 22 only be solved once on the background conductivity model, are \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD L+ i\u00CF\u0089\u00C2\u00B50S i\u00CF\u0089\u00C2\u00B50SG DS DSG \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD Ap \u00CF\u0086p \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8= 0. (2.9) Ap are the primary vector potentials on the mesh and \u00CF\u0086p are the primary scalar potentials. L represents the discretization of the Laplacian operator, S represents the harmonically averaged background cell conductivities, and G and D are the discretizations of the gradient and divergence operators. From Farquharson et al. (2002), the boundary conditions for the primary fields are that the normal compo- nent of J and the tangential component of \u00E2\u0088\u0087\u00C3\u0097A (that is the H-field) are specified. These values are computed by solving the one or two-dimensional magnetotelluric problem. Once the background fields have been calculated, the source term is known for the secondary field calculation. The equations for the secondary field are \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD L+ i\u00CF\u0089\u00C2\u00B50S i\u00CF\u0089\u00C2\u00B50SG DS DSG \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD As \u00CF\u0086s \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8= \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD \u00E2\u0088\u0092iw\u00C2\u00B50SdEp \u00E2\u0088\u0092DSdEp \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 , (2.10) where Sd = S\u00E2\u0088\u0092 Sp is the difference between the averaged conductivities of the actual model (S) and those of the background (Sp), and Ep is the primary electric field. From Farquharson et al. (2002) and (Haber et al. (2000b)), the boundary 23 conditions over the boundary of the domain for the secondary fields are that (\u00E2\u0088\u0087\u00C3\u0097As)\u00C3\u0097n= 0, (2.11a) As \u00C2\u00B7n= 0, (2.11b) \u00E2\u0088\u0082\u00CF\u0086s \u00E2\u0088\u0082n = 0, (2.11c) where n is a unit vector normal to the domain boundary. For each source polar- ization the resulting H fields are computed by summing the computed primary and secondary field components. 2.1.4 Currents and Transfer Functions Since a 1D earth produces no vertical magnetic fields, it is deviations from a 1D conductivity structure that cause a ZTEM response. With the addition of two and three dimensional structure, new anomalous currents will flow in the earth. It is these anomalous currents, Janomalous = Jtotal\u00E2\u0088\u0092Jprimary (2.12) the give rise to vertical magnetic fields. We investigate the currents and associated vertical magnetic fields by considering the models in figure 2.4 of conductive and resistive prisms in a moderate conductivity host. Figure 2.5 shows the anomalous 24 (a) (b) X Y Z X Y Z !\"# $%&'#()*'+)%,#- .%&'#()*'+)%,#- /%&'#()*'+)%,#- !\"#$%&'%\"#$ !#$(%&'%#$( !$(%&'%!\"$( \" 0(*'*1'2'13%,4#- ! \" # $%&'#()*'+)%,#- .%&'#()*'+)%,#- /%&'#()*'+)%,#- !\"#$%&'%\"#$ !#$(%&'%#$( !$(%&'%!\"$( \"(((( 0(*'*1'2'13%,4#- !\"# $%&'#()*'+)%,#- .%&'#()*'+)%,#- /%&'#()*'+)%,#- !\"#$%&'%\"#$ !#$(%&'%#$( !$(%&'%!\"$( \" 0(*'*1'2'13%,4#- ! \" # $%&'#()*'+)%,#- .%&'#()*'+)%,#- /%&'#()*'+)%,#- !\"#$%&'%\"#$ !#$(%&'%#$( !$(%&'%!\"$( \"(((( 0(*'*1'2'13%,4#- Figure 2.4: a) 1 \u00E2\u0084\u00A6m conducting block buried in a 200 \u00E2\u0084\u00A6m background. b) 10000 \u00E2\u0084\u00A6m resistive block buried in a 200 \u00E2\u0084\u00A6m background. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) current density, computed using the 3D forward modelling code described in the previous section, for a vertical slice through the prisms. The anomalous current is fairly uniform and strongest inside the block. The magnetic field and transfer functions (figure 2.6) look like those obtained from an extended dipole current. The anomalous currents for a resistive block flow in the opposite direction compared to the conductive block. The direction of the anomalous current determines the sign of the response while the distance between positive and negative maxima is related to the size of the block. 25 (a) (b) Figure 2.5: The real anomalous current densities (Log10( Am2 )) in the z-y plane (x= -25 m) for a source whose electric field is polarized in the (y\u00CC\u0082) direc- tion. The white rectangles are the outlines of the original blocks. Figure a) is the current for the conductive block, while figure b) is the current for the resistive block. The anomalous current is mostly localized in- side the block. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) 2.2 Inverse Problem For simple geometries the concepts of anomalous currents can help interpret ZTEM data. For instance, a linear filter can be applied to the vertical magnetic field (Karous and Hjelt, 1983) to generate current pseudosections. For this simple block in a halfspace example, the block geometry may be adequately interpreted from us- ing these current pseudosections. Nonetheless, in the case of complex 3D models, 26 (a) X (m)X (m) Y (m ) X (m) Y (m ) Y (m ) X (m) Y (m ) (b) X (m)X (m) Y (m ) X (m) Y (m ) Y (m ) X (m) Y (m ) Figure 2.6: Resulting transfer functions for the two block example of fig- ure 2.4. Figure a) corresponds to a conductive block while figure b) corresponds to a resistive block. The transfer functions closely resem- ble the magnetic fields due to a finite line current inside the blocks. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) 27 3D inversions are necessary. The geophysical inversion problem is stated as follows: given a dataset, mea- surement errors, and forward modelling operator, find a model that generates the data measurements. Because airborne frequency domain electromagnetic data are only collected at a few frequencies on a surface above the earth, the inverse prob- lem is highly non-unique because data is only being collected at a few points over the surface of the earth and then being inverted to recover a 3D conductivity model. In order to obtain a solution to the inverse problem we must regularize the solution by attempting to incorporate as much additional information into the inversion as possible. Our inversion algorithm is based on the MT inversion code of Farquhar- son et al. (2002). By minimizing the objective function \u00CE\u00A6= \u00CF\u0086d +\u00CE\u00B2\u00CF\u0086m, (2.13) we obtain our solution to the inverse problem. \u00CF\u0086d is the data-misfit, \u00CF\u0086m is the amount of structure in the model, and \u00CE\u00B2 is the trade-off or regularization parameter. We use the sum of squares as the measure of misfit \u00CF\u0086d = ||Wd(dobs\u00E2\u0088\u0092dprd)||22, (2.14) where dobs is the observation vector, dprd is the vector of predicted data and Wd is 28 a diagonal matrix whose elements are the reciprocals of the standard deviations of the data errors. The measure of model structure is \u00CF\u0086m = \u00CE\u00B1s||Ws(m\u00E2\u0088\u0092mre f )||22 +\u00CE\u00B1x||Wx(m\u00E2\u0088\u0092mre f )||22 +\u00CE\u00B1y||Wy(m\u00E2\u0088\u0092mre f )||22 +\u00CE\u00B1z||Wz(m\u00E2\u0088\u0092mre f )||22, whereWs is a diagonal matrix andWx,Wy andWz are the first-order finite-difference matrices in the x, y and z directions, and mre f is a reference model. The \u00CE\u00B1 \u00E2\u0080\u00B2s are adjustable parameters. \u00CE\u00B1s controls the closeness of the recovered model to the ref- erence, while \u00CE\u00B1x,y,z determine the smoothing in the x, y and z directions. For the minimization of the objective function at the (n+1)th iteration, the Gauss-Newton method requires the solution of (JTWTdWdJ+\u00CE\u00B2W TW)\u00CE\u00B4m= JTWTdWd(d obs\u00E2\u0088\u0092dn)\u00E2\u0088\u0092\u00CE\u00B2WTW(mn\u00E2\u0088\u0092mre f ), (2.15) where mn is the vector of model parameters from the previous iteration, J is the Jacobian matrix of sensitivities, W is such that WTW =\u00E2\u0088\u0091\u00CE\u00B1kWTkWk, and \u00CE\u00B4m is the model perturbation. Here the entries of Wk are normalized against the volume of each cell. Since ZTEM data are different from MT data, the Jacobian matrix must be modified accordingly. The implementation however follows that outlined 29 in Farquharson et al. (2002). Equation 2.15 is solved using an inexact precondi- tioned conjugate gradient (IPCG) (Saad (2003)) algorithm preconditioned using an incomplete LU decomposition of the matrix (WTW+ 0.1I), where I is the iden- tity matrix. Typically, 5 IPCG iterations are performed to solve equation 2.15 for each value of \u00CE\u00B2 . Since J is prohibitively large to compute, this solution method only requires the sparse matrix operations WTW,J,JT on a vector. \u00CE\u00B2 is initially chosen to be large such that \u00CE\u00B2WTW dominates the approximate Hessian in equa- tion 2.15. It is then reduced in a cooling schedule such that at each subsequent iteration \u00CE\u00B2k+1 = \u00CE\u00B3\u00CE\u00B2k, where \u00CE\u00B3 is a constant (the default value of \u00CE\u00B3 is 17 ). Efficient implementation of this inversion in application to field datasets requires a workflow procedure but we postpone details about that until section 4. 2.3 Synthetic Inversion We illustrate the inversion of ZTEM data using a synthetic example that was cho- sen to emulate a porphyry deposit. The model has a central resistive core, an outer region of high conductivity, and is buried in a moderate conductivity background. The conductivity structure is shown in figure 2.7. The Earth was discretized into a mesh containing 42 x 42 x 57 cells in the x, y, and z directions respectively. The cell dimensions in the central core region were 250 x 250 x 100 m in the x, y, and z directions respectively. The vertical fields at 1, 3.2, 5.6 , 10 , 18 , and 32 Hz were computed at a fixed height of 100 m over the surface. This frequency range 30 is shifted somewhat from the traditional ZTEM survey but this is inconsequential for the purpose of testing the code. The data were computed at 10 m intervals and 50 m line-spacing over a 5 km by 5 km region. The horizontal field components were computed at a ground based reference station at a location r0 of x = - 3000 m, y = - 3000 m. For each component of each frequency Gaussian noise, equal to 3.5 percent of the difference between the 10 th percentile and 90 th percentile data, was added. The initial model and reference models were homogeneous halfspaces of 100 \u00E2\u0084\u00A6m. We first inverted different subsets of the data to determine their informa- tion content. The 1Hz data were inverted with the real and imaginary components separately and then the 1Hz data were inverted with both the real and imaginary components. From the two sections of the model inversions in figure 2.8, it is clear that both the real and imaginary components contain significant information and that both components of the data should be obtained in a ZTEM survey. Since 1Hz is the lowest frequency used, it also resolves the deepest structure. Next all six fre- quencies were inverted simultaneously. The inversion was executed in parallel over 6 dual processors. The convergence curve for the inversion is shown in figure 2.9. The recovered model fits the data well and select data are shown in figure 2.10. The final model in figure 2.11 recovers both the conductivity and block dimensions in the x and y directions very well. There is significant improvement in the block res- olution in the x and y directions over the single 1Hz inversion. The inverted result 31 (a) X Dimension (m) Y Dimension (m) Z Dimension (m) -500 to 500 -1000 to 1000 -500 to -3500 1000 -1500 to 1500 -2000 to 2000 -500 to -3500 10 Resistivity (\u00CE\u00A9m) X Y Z S/m (b) X Y Figure 2.7: Conductivity model: a resistor surrounded by a conductor, all buried in a 100 \u00E2\u0084\u00A6m background halfspace. (a) Cross section through y = 0 m. The bottom table summarizes the block dimensions and resistiv- ities. (b) Depth slice at 1000 m. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) 32 shows that the top interface is well resolved, however, the base of the blocks are not. This is consistent with skin depth arguments of the lowest frequency. 2.3.1 Effect of the Reference Station One component of the ZTEM system is the reference station that is used to compen- sate for the unknown source amplitude. The transfer functions related the vertical magnetic field to the horizontal magnetic field at the reference station. Ideally the fields at the reference station are those due to a 1D conductivity structure. Gen- erally, the variation between the horizontal field at the reference station and the observation point are small even if there are lateral changes in conductivity. To illustrate that the choice of reference station does not drastically change the in- version result, we compute the relative difference, ||(Hx(r)\u00E2\u0088\u0092Hx(rre f )||||Hx(rre f )|| for one source polarization on the synthetic inversion model using a reference station at x = -3000 m, y = -3000 m. As can be seen in figure 2.12 the variations in the horizontal field are less than 9.5% over the entire survey area. The location of the maximum deviations occurs at x = -900 m, y = 0 m. We now compute a new dataset using a reference station located at this point of max- imum deviation. We invert this new dataset using the same inversion parameters as an inversion with a reference station at x = -3000 m, y = -3000 m. The com- parison between the two inversion results can be seen in figure 2.13. As we would expect, when the reference location is correctly modelled, the inversion result is 33 (a) X Y Z X Y (b) (c) (d) (e) Figure 2.8: Conductivity model (S/m) from the inversion of the synthetic model data. Figure a) is a cross section at the center of the model and a planview section at a depth of 800 m. b-d) are results of inverting data at 1 Hz. b) is the inversion of the real parts, c) is the imaginary parts, and d) is the real and imaginary components. e) is the simultaneous in- version of 6 frequencies. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) 34 Figure 2.9: Convergence curve for the inversion of six frequencies for the synthetic inversion model. The algorithm reduced the misfit to 1.38\u00C3\u0097 106 (target misfit 1.21\u00C3\u0097 106). ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) not strongly dependent upon the location of the reference station. A more interesting example concerns the result of working with an incorrect location of the reference station. In some surveys, the model domain can be exces- sively large if the reference station is located far from the survey area. It would be desirable to have an artificial location for the reference station inside the modelling domain. We simulate this by computing the data using the reference location out- side the survey area (x = -3000 m, y = -3000 m), and then positioning the reference location at x = -900 m, y = 0 m during the inversion. The results from this inver- sion in figure 2.13 show that even if the reference station position is moved in the inversion, we still recover a reasonable model. 35 (a) X (m)X (m) Y (m ) X (m) Y (m ) Y (m ) X (m) Y (m ) (b) X (m)X (m) Y (m ) X (m) Y (m ) Y (m ) X (m) Y (m ) Figure 2.10: Figure a) is the observed data and figure b) is the predicted data for the synthetic model at 1 Hz. There is excellent agreement between the observed and predicted data. 36 (a) (b) (c) Figure 2.11: Conductivity model (S/m) from the inversion of the synthetic model using 6 frequencies. Slices of the true model (left hand fig- ures) and recovered model (right hand figures) are shown at 500, 1000, 1500m. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by per- mission from Holtham and Oldenburg (2010)) 37 Figure 2.12: Percent difference in the x\u00CC\u0082 component of the magnetic field be- tween the reference location and the measurement location (source electric field in the y\u00CC\u0082 direction). ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) 2.4 3D Inversion of ZTEM Field Data In the fall of 2008, a ZTEM survey was flown for Geoscience British Columbia and Terrane Metals Corp. over the Mt. Milligan test block located in north central British Columbia. The Mt. Milligan deposit which was owned by Terrane Met- als Corp. in 2008 and is currently owned by Thompson Creek Metals Company Inc., is an alkalic Cu-Au porphyry deposit located within the Quesnel Terrane (see figures 2.14 and 2.15). The initial deposit outcrop was first discovered by a local prospector in 1983 and since then the deposit and surroundings areas have seen a wealth of exploration activity including geochemistry, geologic mapping (Nelson et al., 1990), ground and airborne geophysics (Oldenburg et al., 1997; Yang and 38 (a) X Y Z X Y (b) (c) Figure 2.13: Conductivity model (S/m) from the inversion of the synthetic model data using different reference station locations. For each inver- sion all of the parameters except the location of the reference station remain unchanged. Figure a) corresponds to a reference station located at x = -3000 m, y = -3000 m. Figure b) corresponds to a reference sta- tion located at x = -900 m, y = 0 m. Figure c) corresponds to data computed using a reference station at x = -3000 m, y = -3000 m, but during the inversion process the reference station is mislocated at x = -900 m, y = 0 m. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) 39 Oldenburg, 2012) as well as extensive drill programs. While the geologic details and alteration events of every deposit are unique, being classified as an alkalic porphyry deposit, Mt. Milligan can be assumed to have the general structure and alteration events as shown in figure 2.16. From the Terrane Metals Corp. 2009 NI 43-101 report (Mills et al., 2009), Mt. Milligan is a tabular, near-surface, alkalic copper-gold porphyry deposit that mea- sures some 2,500 m north-south, 1,500 m east-west and is 400 m thick. The impor- tant geologic features of the deposit and the surrounding areas are shown in figure 2.17. From Mills et al. (2009), the Main Zone is spatially associated with the MBX monzonite stock (see figure 2.17) and Rainbow Dyke. The mineralization and asso- ciated alteration are primarily hosted in volcanic rocks. Mineralization consists of pyrite, chalcopyrite and magnetite with bornite localized along intrusive-volcanic contacts. Copper-gold mineralization is primarily associated with potassic alter- ation which decreases in intensity outwards from the monzonite stocks. Pyrite content increases significantly outward from the stocks where it occurs in asso- ciation with propylitic alteration, which forms a halo around the potassic-altered rocks. The Mt. Milligan copper-gold porphyry deposits contain Proven and Prob- able reserves of (Mills et al., 2009) 482 million tonnes (Mt) averaging 0.20 % Cu and 0.39 grams per tonne (g/t) Au totaling 2.1 billion pounds copper and 6.0 mil- lion ounces gold. 40 While the deposit has been extensively drilled over the years, little information is available on the expected conductivity of the different rock units and alteration events. Recently a Geoscience BC initiative has been undertaken to try and im- prove this understanding by performing laboratory studies on core samples taken from several British Columbia porphyry systems (Mitchinson and Enkin, 2011). While more research and down-hole conductivity measurements are needed to gain a more complete understanding how lithology, alteration, and mineral grade con- trol the conductivity on the deposit scale, it is nonetheless possible to understand the relative conductivity of several geologic units in the area. Firstly, when un- altered, the intrusive monozonite stocks are likely the most resistive structures in the deposit area. The various alteration events distal to the stocks are expected to be more conductive than the unaltered intrusives. Finally, the faults and tertiary sediments are likely the most conductive structures in the area. The 2008 ZTEM test survey was flown over the Mt. Milligan deposit with the hope of showing that electromagnetic methods can be used to find additional deposits that share similar features to the Mt. Milligan deposit. The test survey consisted of 25 lines over an approximate area of 6 km x 8 km for a total of 211 line-km of data (see figure 2.18). The data were collected as time series which were then stacked and processed into transfer functions at five ZTEM frequencies (30, 45, 90, 180 and 360 Hz) using proprietary Geotech software. The results were 41 Figure 2.14: Location of Mt. Milligan after Devine (2011b) in north cen- tral British Columbia. The deposit is located within the Quesnel Ter- rane, a well known porphyry belt that extends through much of British Columbia. filtered and attitude corrected to compensate for the pitch and roll of the coil. The terrain in the survey block is moderately severe and varies between 934 m to 1475 m (see figure 2.18) over a relatively small area which increases the difficulty of interpreting and inverting the data. The 3D inversion of any geophysical field dataset is a complex and time con- suming process that requires a logical work flow in order to be completed accu- rately and efficiently. We present our workflow in (figure 2.19). In this section of 42 Figure 2.15: Regional geology around Mt. Milligan after Devine (2011b) the chapter we invert the field data and at the same time illustrate some specific parts of this workflow. 2.4.1 Discretizing the Earth In order to solve the necessary differential equations, the earth must be discretized into a mesh and the discrete system for Maxwell\u00E2\u0080\u0099s equations formed. In this dis- cretization step a compromise must be achieved between the numerical accuracy in the forward modelling and the computational resources required. As with MT discretization problems, the usual guidelines that cell sizes must not violate skin depth rules apply. Also, since ZTEM data are sensitive to topography, it is im- portant to finely discretize the topography. These requirements, when combined 43 Figure 2.16: Geologic cartoon of a typical alkalic porphyry deposit from Devine (2011a) with the high frequency and large geographic areas covered in ZTEM surveys, can require large meshes and significant computational overhead. The topography for the inversions was obtained from U.S. Geological Survey (USGS) Global Multi-resolution Terrain Elevation Data (GMTED). The Earth was discretized into a mesh containing 188 x 136 x 85 cells in the easting, northing and vertical directions. The central region had 172 x 120 x 60 cells that were 50 m x 50 m x 25 m in the easting, northing, and vertical directions. This 25 m vertical cell dimension allowed for the topography to be reasonably modelled. 44 Figure 2.17: Planview section after Mitchinson and Enkin (2011) detailing the main geologic features in the area. The main features of interest are the monzonite intrusives as well as the tertiary sediments to the east of the deposit. 2.4.2 Determining a Background Model In the MT problem, where electric fields are measured, the background can be cho- sen either using the apparent resistivities or by choosing the halfspace conductivity which minimizes the data misfit between the halfspace data and the true data. These methods cannot be used for ZTEM since all 1D conductivity models produce zero data. Nevertheless, the anomalous currents and data are affected by the background conductivity for 3D models. Our experience shows that it is important to initialize the inversion with a good estimate of the background conductivity, especially if the anomalous bodies are deeply buried. If this is not done, then the inversion may recover a reasonable estimate of the relative conductivities, but fail to adequately 45 Figure 2.18: Topography over the Mt. Milligan test ZTEM block. The to- pography over the survey area ranged from 934 m to 1475 m for a total of 541 m of relief. The flight lines are shown in black and the outlines of the main monzonite intrusives are outlined in red to get a sense of scale from the deposit scale geologic map in figure 2.17. reproduce the data and also be a poor approximation to the true model. We il- lustrate these principles by starting the 1Hz inversion of the synthetic model with an initial background of 10000 \u00E2\u0084\u00A6m instead of the 100 \u00E2\u0084\u00A6m true background. The model smallness parameter, \u00CE\u00B1s, was set to be 0 so that the reference model does not affect the inversion. From the misfit curve in figure 2.20 we see that the inversion was initially able to significantly reduce the misfit. During the first part of the in- 46 Obtain the data View and understand the data Obtain topography Discretize, validate forward modelling Determine background model for inversion Assign uncertainties Invert each frequency individually Evaluate Inversion results Simultaneously invert all data Evaluate results Final Result Incorporate additional information Figure 2.19: Work flow process summarizing the critical step in a ZTEM in- version. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by per- mission from Holtham and Oldenburg (2010)) version, some estimate of the relative conductivity has been established, however, it becomes increasingly difficult to fit the data since the background conductivity is drastically incorrect. Eventually, the inversion algorithm terminates with the final model in figure 2.20. Clearly this is not a satisfactory result. This inversion result motivated the following procedure. We used the inversion shown in figure 2.20 but selected a model just prior to the sharp change in slope on the misfit curve. We assumed that the relative conductivities on that model were 47 (a) (b) (c) (d) Figure 2.20: a) Conductivity (S/m) of the synthetic true model at a depth of 800 m. b) Misfit curve. c) Inverted model at 800 m, using a 10000 \u00E2\u0084\u00A6m starting model, taken at the sharp change of slope on the misfit curve. d) Final inversion model starting with a 10000 \u00E2\u0084\u00A6m initial model. The recovered model structure and conductivity values are not reasonable. Here \u00CE\u00B1s = 0, \u00CE\u00B1x = \u00CE\u00B1y = \u00CE\u00B1z = 1. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permission from Holtham and Oldenburg (2010)) 48 correct. That conductivity was scaled by various constants and the data from the models were computed. The misfit between the true data and data produced from the scaled models is shown in figure 2.21. The smallest data misfit occurred when the conductivity of the preliminary model was scaled by a factor of 10. The in- verted conductivity model underneath the 5 km by 5 km survey area of the scaled model was then averaged to obtain a background conductivity. This new conductiv- ity was determined to be 169 \u00E2\u0084\u00A6m instead of the actual 100 \u00E2\u0084\u00A6m. This procedure of determining relative conductivities, scaling, and averaging, provides a method for determining appropriate starting models even if the starting background conductiv- ity is incorrect by a few orders of magnitude. The final inversion result using the new 169 \u00E2\u0084\u00A6m background is shown in figure 2.22. It is a significant improvement over the preliminary inversion in figure 2.20. With significant topography, even a uniform conductivity model will generate charge build up and anomalous currents. For these geometries it may be possible to find the conductivity that best fits the given data. This method is illustrated using a synthetic model based on the Mt. Miligan ZTEM survey (survey layout shown in figure 2.18). Using the Mt. Milligan topography, a synthetic model was created with a constant 250 \u00E2\u0084\u00A6m conductivity. The response from just the topography was then forward modelled at 90 Hz. Uniform background models (10, 50, 100, 250, 500, 1000 \u00E2\u0084\u00A6m) with the correct topography were forward modelled to obtain data 49 Figure 2.21: Data misfit for various conductivity scaling factors. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permis- sion from Holtham and Oldenburg (2010)) (a) (b) Figure 2.22: a) Inverted model(S/m) at a depth of 800 m starting with the correct initial background conductivity. b) Inverted model at a depth of 800 m using the scaling method to determine the initial conductiv- ity. ( c\u00C2\u00A9Geophysical Journal International, 2010, adapted by permis- sion from Holtham and Oldenburg (2010)) 50 Easting (m) Elevation (m ) Figure 2.23: Synthetic test example of the Mt. Milligan topography to show that the method of fitting a best fitting halfspace to the data may be used if the topography response is much greater than the geologic response. The synthetic model has a constant resistivity of 250 \u00E2\u0084\u00A6m. at the same locations. The data misfits between the various conductivity models and the true 250 \u00E2\u0084\u00A6m model were computed and are plotted in figure 2.24. As expected, the background that produces the smallest data misfit matches the true conductivity. As the conductivity gets further from the true value, the data misfits increase. This method may allow for background models to be created when the topographic response is strong relative to the geologic response. 2.4.3 Assigning Uncertainties Assigning uncertainties to data is a critical step in any geophysical inversion. Over- estimating the noise will result in a loss of information which could have been extracted from the data. Similarly, underestimating the noise may create artifacts 51 Figure 2.24: Misfit between the true data and the data produced from differ- ent conductivity background models. As the conductivity values get further from the true conductivity, the data misfit increases. in the inverted model that are not supported by the data. Assigning uncertainties for any geophysical field dataset is difficult because there are many errors sources, each of which is difficult, if not impossible to quantify. Survey errors such as in- strument repeatability, human error, external noise, and mis-orientation of receivers are often thought to be the only sources of errors that need to be considered in an inversion. These sources of errors are important, however, there are many additional sources of errors. Another source of error is discretizing the earth into rectan- gular cells. Since both the topography and conductivity vary on fine scales, by discretizing the earth into relatively large cells, the earth is only being approxi- mately represented. Errors are also associated with modelling difficulties. In nat- 52 ural source techniques, approximations such as the quasi-static approximation and plane wave sources are made. Furthermore, the forward modelling equations are iteratively solved which can lead to convergence issues and solution errors if the solution tolerance is not small enough. Ultimately, it is difficult to quantify the individual error sources, let alone their sum. Nevertheless it is essential that there is a robust and repeatable method to assign uncertainties to the data. Generally, uncertainties for geophysical data are assigned as some constant value or a floor across a data component, plus some percentage of the individual datum. Using this approach helps represent various sources of noise in the data, modelling, and inversion. While ZTEM data generally has a small dynamic range, the data contains many zero crossings. Because of this, the constant floor value assigned to the data can be quite important, and make up a significant portion of the assigned errors. Due to the remote location of the Mt. Milligan deposit (the deposit was not being put into production at the time of the ZTEM survey), the ZTEM data was of high quality and devoid of significant cultural noise such as powerlines and other made structures. The errors on the data were assigned as a constant floor of 0.015 on each datum plus a percentage error which was based on the frequency of the data. 53 2.4.4 Mt. Milligan Inversion Results Both the real and imaginary components of 5 frequencies were inverted for the Mt. Milligan ZTEM survey. For the inversion, the earth was discretized on a mesh that contained 188, 136 and 85 cells in the easting, northing and vertical directions respectively. The cell dimensions in the central core region were 50 m x 50 m x 25 m. The 180 Hz observed and predicted data can be seen in figures 2.25 and 2.26. Figure 2.27 shows the normalized data misfits. The inversion was completed using no prior knowledge of the geology, so no additional information was incorporated into the inversion, and no bound constraints were applied to the inversion. The initial and reference models had a constant resistivity of 250 \u00E2\u0084\u00A6m, and the model smallness parameter, \u00CE\u00B1s was set to be 2.5 \u00E2\u0088\u0097 10\u00E2\u0088\u00926. The inversion took 33 hours to run on two Xeon X5660 dual hex core processors. After the inversion was completed, the geologic model for the deposit area (fig- ure 2.17) was examined. The ZTEM inversion recovers a very conductive structure east of the deposit, due to the Great Eastern Fault Zone and the tertiary sediments. The largest MBX stock sits in a region of relatively low conductivity, and surround- ing the stock on the north is a region of higher conductivity, likely due to alteration. As a comparison to the ZTEM inversion we can examine the 3D inversion of a ver- satile time-domain electromagnetic (VTEM) survey, also flown by Geotech Inc., and inverted in Yang and Oldenburg (2012). The VTEM inversion result in figure 54 Easting (m) N or th in g (m ) N or th in g (m ) N or th in g (m ) N or th in g (m ) Easting (m) Figure 2.25: Observed data over Mt. Milligan at 180 Hz. Real part of Tzx, imaginary part of Tzx, real part of Tzy, and the imaginary part of Tzy. Easting (m) N or th in g (m ) N or th in g (m ) N or th in g (m ) N or th in g (m ) Easting (m) Figure 2.26: Predicted data over Mt. Milligan at 180 Hz. Real part of Tzx, imaginary part of Tzx, real part of Tzy, and the imaginary part of Tzy. The inversion fit the data well. 55 Easting (m) N or th in g (m ) N or th in g (m ) N or th in g (m ) N or th in g (m ) Easting (m) Figure 2.27: Normalized misfits (normalized by the standard deviations of the data) over Mt. Milligan at 180 Hz. Real part of Tzx, imaginary part of Tzx, real part of Tzy, and the imaginary part of Tzy. The misfits generally have little correlation with the observed data. 2.28 shows a very similar conductivity structure to the ZTEM inversion result at the same depth in figure 2.30. 56 6110400 6109950 6109500 6109050 6108600 N o rt h in g ( m ) \u00E2\u0080\u00933.64 \u00E2\u0080\u00933.24 \u00E2\u0080\u00932.84 \u00E2\u0080\u00932.44 \u00E2\u0080\u00932.04 \u00E2\u0080\u00931.63 \u00E2\u0080\u00931.23 \u00CF\u0083 (log S/m) 43 35 00 43 39 50 43 44 00 43 48 50 43 53 00 43 35 00 43 39 50 43 44 00 43 48 50 43 53 00 Easting (m)Easting (m) (a) (b) MBX Stock DWBX 66 MBX Zone T e r tia r y s e d im e n ts 4 2 0 \u00CE\u00A9 m 4 2 0 \u00CE\u00A9 m < 4 2 0 \u00CE\u00A9 m < 4 2 0 \u00CE\u00A9 m >420\u00CE\u00A9m Figure 2.28: Depth slice of the VTEM inversion results over the Mt. Milli- gan deposit from Yang and Oldenburg (2012). (a) Slice at an elevation of 1030 m overlain by geology and (b) 420 \u00E2\u0084\u00A6m contour of DC resis- tivity model. The MBX stock is moderately resistive in both the 3D VTEM and DC resistivity inversions. ( c\u00C2\u00A9Society of Exploration Geo- physicists, 2012, adapted by permission from Yang and Oldenburg (2012)) 57 Easting (m) N or th in g (m ) S/m Figure 2.29: Depth slice of the ZTEM inversion results over the Mt. Milligan deposit at an elevation of 1028m, the closest depth section to the 1030 m slice in figure 2.28. Both the 3D VTEM and ZTEM inversions recover a moderately resistive MBX stock. Also the conductive region north of the MBX stock is well mapped in both inversions. Easting (m) N or th in g (m ) S/m Figure 2.30: Depth slice of the ZTEM inversion results over the Mt. Milligan deposit at a depth of 728m. The conductive Great Eastern Fault zone and Tertiary Sediments are still well mapped at this depth. 58 Chapter 3 Large-scale Inversion of ZTEM Data 3.1 Introduction Since most shallow and outcropping deposits have already been discovered, many of the earth\u00E2\u0080\u0099s remaining natural resources are buried deep and under cover mak- ing them difficult to find. In order to satisfy the strong global resource demand, techniques to explore for these big targets must be developed. Electromagnetic methods can be used to map electrical conductivity which can then be linked to ge- ologic features of interest. When exploring for large buried targets, natural source electromagnetic methods can be advantageous over controlled source methods be- 59 cause of the deeper penetration of plane wave sources. While the traditional natural source method, the magnetotelluric (MT) method, has been effectively applied to both mining and hydrocarbon exploration, a practi- cal limitation of the MT technique is that surveys are costly and time consuming, since many expensive stations must be installed to measure the needed electromag- netic field components on the earth\u00E2\u0080\u0099s surface. It was the desire to collect airborne natural source data that prompted the development of AFMAG (Audio Frequency Magnetics) (Ward, 1959). Unfortunately, because the direction of the inducing source fields varies with time and the derived data were scalar instead of a vec- tor quantities, AFMAG results were not always repeatable (limitations of AFMAG are outlined in Ward et al. (1966)). More recently Geotech Ltd. has modified the AFMAG concept and developed the ZTEM (Lo and Zang, 2008) technique. In the ZTEM technique the vertical component of the magnetic field is recorded above the entire survey area, while the horizontal magnetic fields are recorded at a ground- based reference station. MT processing techniques yield frequency domain transfer functions typically between 30-720 Hz that relate the vertical fields over the survey area to the horizontal fields at the reference station. While the ZTEM method has no sensitivity to purely 1D conductivity structures (purely 1D structures produce zero vertical magnetic field component), the method is highly sensitive to lateral conductivity contrasts. If ZTEM inversions are started with good estimates for the 60 background conductivity, then the inversions can recover reasonable conductivity estimates. Large ZTEM datasets can be effective exploration tools, particularly when ex- ploring on the district and regional scale where subsurface geology is hidden under cover. However, in order to justify the financial expense of collecting the data, viable interpretation methods must exist. This requires that inversions are fairly fast, especially for field data where multiple inversions are often performed using different parameters. Unfortunately, for surveys that cover large areas, and where even moderate inverted resolution is desired, the discrete 3D Maxwell systems for solving the electromagnetic (EM) problem quickly become very large. Although computer technology continues to improve, clock speed and instruction level par- allelization have started to flatten since 2003 (Shalf, 2007). This means that the majority of future improvements in the size of inverse problems that can be tackled must come from improved methodologies and increased parallelization. One viable method to solve large inverse problems is to reduce the modelling domain by using footprint methods, where computations are simplified by only considering model cells that influence the data above some threshold criterion. These methods have been used to invert large airborne controlled source elec- tromagnetic datasets (Cox et al., 2010) and MT datasets (Gribenko et al., 2010). Another solution is to use domain decomposition methods (DDM) and split the 61 computational domain into smaller manageable subproblems which can be solved quickly in parallel before being merged together to form the final solution. The concept of using domain decomposition methods and splitting the compu- tation domain into smaller subproblems has long be used in applied math to solve boundary value problems (for example see Saad (2003)). One of the simplest and earliest domain decomposition procedures proposed by Schwartz (1870) was to split the computation domain into two overlapping domains, and alternate between the subproblems, while updating the parameters on the next domain based on the recent solution from the other domain (Saad, 2003). Another more parallel varia- tion on this algorithm, the additive Schwartz procedure, differs from the previously described procedure in that the components in each subdomain are not updated un- til a whole cycle of updates through all domains are completed (Saad, 2003). In this chapter we present a highly parallel and effective procedure, based loosely on Schwartz alternating procedures, to invert large natural source surveys. 3.2 Inversion Methodology While this chapter focuses on inverting ZTEM data, the approach and methodolo- gies are equally applicable to MT as well as a combination of MT and ZTEM data (to be discussed in chapter 4). Our original natural source inversion algorithm is that of Farquharson et al. (2002), and our ZTEM inversion algorithm (discussed in chapter 2) was created by modifying this code. We wish to solve the inverse 62 problem by minimizing the objective function \u00CE\u00A6= ||Wd(F(m)\u00E2\u0088\u0092dobs)||22+\u00CE\u00B2 (\u00CE\u00B1s,x,y,z||Ws,x,y,z(m\u00E2\u0088\u0092mref)||22), (3.1) where, m is the full fine cell model, mre f is a reference model, F is the forward modelling operator, and dobs is the observation vector. Wd is a diagonal matrix whose elements are the reciprocals of the standard deviations of the data errors, Ws is a diagonal matrix, Wx, Wy, and Wz are the first-order finite-difference ma- trices in the x, y and z directions. The \u00CE\u00B1 \u00E2\u0080\u00B2s are adjustable parameters. \u00CE\u00B1s controls the closeness of the recovered model to the reference, while \u00CE\u00B1x,y,z determine the smoothing in the x, y and z directions. \u00CE\u00B2 is the regularization parameter which is reduced throughout the inversion process until the desired data misfit, \u00CF\u0086 \u00E2\u0088\u0097d , has been achieved. Although the goal is to solve this problem, it is difficult to tackle directly because of the size of the problem and the cost of evaluating F(m). In order to avoid inverting the full domain directly, we use a algorithm based on the Schwartz alternating procedure (Saad, 2003). The basic concept of Schwartz methods is to split the full domain into s smaller domains\u00E2\u0084\u00A6i, i= 1...s. Using the notation consistent with Saad (2003), the boundary of subdomain \u00E2\u0084\u00A6i that is included in subdomain j is given as \u00CE\u0093i j (see for example figure 3.1). The basic Schwartz alternating procedure is then given as, 63 Figure 3.1: Example of an L-shaped domain split up into 3 smaller domains with two shared boundaries after Saad (2003). Algorithm 1 Schwartz Alternating Procedure Choose an initial guess u to the solution while Not converged do for i= 1, ...,s do Solve problem in \u00E2\u0084\u00A6i with current boundary conditions on \u00CE\u0093i j end for Update values on all the boundaries \u00CE\u0093i j end while The key elements in the algorithm are choosing an initial guess to the solu- tion, solving the problem on each subdomain, and updating the solution on the boundaries. In the next section these three basic concepts are incorporated into our algorithm (algorithm 2) to solve large-scale natural source problems. The first step is to determine an initial guess to the solution. This is accom- plished by first inverting the full dataset on a coarse mesh to determine an initial 64 coarse inversion result, mc. This coarse mesh result can be interpolated onto the desired fine mesh to form an initial solution model, m= Lfcmc, (3.2) to the full inverse problem. Here Lfc is an interpolation operator that goes from the coarse mesh to the fine mesh. Next, we can decompose the full fine-scale model domain m, and the data set into T smaller tiles or domains, m\u00CC\u0083t = Pmt m, t = 1, ..,T, (3.3) and d\u00CC\u0083t = Pdt d obs, t = 1, ..,T, (3.4) where Pmt and Pdt are projection operators which form the model and data subdo- mains respectively. Here Pmt and Pdt are constructed such that the m\u00CC\u0083t\u00E2\u0080\u0099s and d\u00CC\u0083t\u00E2\u0080\u0099s are overlapping. We also define the projection operator, u\u00CC\u0083t = Put u, t = 1, ..,T, (3.5) which maps the fields, u, from the full conductivity structure to each subdomain. We would like to invert these T overlapping subproblems separately and in paral- 65 lel, however, in doing so directly we would neglect interactions between subdo- mains. In order to mitigate this problem we can incorporate the domain interac- tions through source terms obtained on a full domain model. To accomplish this, the coarse model mc is first interpolated onto the full fine mesh to form an initial solution, m1. This model is then forward modelled to compute initial fields u1 at each mesh location. Both m1 and u1 are then projected onto the T domains to form m\u00CC\u00831t and u\u00CC\u00831t . Although these new fields and models are only defined on an individ- ual subdomain, they were still obtained by solving the inverse problem on the full model domain and hence contain information about the conductivity outside their respective domains. To understand how the fields computed from a full conductivity structure can be exploited in the subdomain inversions, we examine our forward modelling pro- cedure of Farquharson et al. (2002) (the solution of Maxwell\u00E2\u0080\u0099s equations is that of Haber et al. (2000b)) where the total fields are computed by summing com- puted primary and secondary fields. The electric field is decomposed into vector and scalar potentials E = A+\u00E2\u0088\u0087\u00CF\u0086 and the Coulomb gauge condition \u00E2\u0088\u0087 \u00C2\u00B7A = 0 is imposed for uniqueness. The secondary fields, As and \u00CF\u0086s, can then be computed 66 as, \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD L+ i\u00CF\u0089\u00C2\u00B50S i\u00CF\u0089\u00C2\u00B50SG DS DSG \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD As \u00CF\u0086s \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8= \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD \u00E2\u0088\u0092iw\u00C2\u00B50(S\u00E2\u0088\u0092Sp)Ep \u00E2\u0088\u0092D(S\u00E2\u0088\u0092Sp)Ep \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 . (3.6) Here L represents the discretization of the Laplacian operator, S represents the harmonically averaged cell conductivities, G and D are the discretizations of the gradient and divergence operators, and Sp contains the averaged conductivities of the primary field model, and Ep is the electric field computed on the primary con- ductivity structure. We can write the block matrix system in equation 3.6 more compactly as, Bus = q(S\u00E2\u0088\u0092Sp)Qup (3.7) where Q is some operator that maps the A-\u00CF\u0086 solution up to Ep. In this approach, the primary conductivity structure and fields are obtained from a full conductivity structure and used as source terms on the right hand side of equation 3.6. Using full model domain source terms ensures that the large-scale physics of the problem will be modelled more accurately. For each subdomain mesh we minimize the objective functions, \u00CE\u00A6\u00CC\u0083t = \u00CF\u0086d(m\u00CC\u0083t)+\u00CE\u00B2\u00CF\u0086m(m\u00CC\u0083t) (3.8) to solve the tth inverse problem. Once the m\u00CC\u0083t\u00E2\u0080\u0099s have been determined, the kth 67 update to the full model mk, is given by merging each subdomain model together, mk =M(m\u00CC\u0083kt=1:T ), (3.9) where M is an operator that merges the subdomain conductivity structures. Ad- ditional details on how the tiles are merged are discussed later in the subsection, model updates and merging tiled inversions. We can determine if the current model mk is a satisfactory solution to our full inverse problem by forward modelling this new conductivity structure and computing the data misfit \u00CF\u0086d . If the data misfit falls within that allowed by the initial inverse problem in equation 3.1, then a solution to the desired inverse problem has been obtained. If the data have not been suffi- ciently fit, then we can perform another iteration of this procedure where we now use the fields and models from mk as the source terms in the kth + 1 subdomain inversions. This procedure outlined in algorithm 2, can be repeated until the data has been sufficiently fit. 3.3 Synthetic Example In this section a synthetic model is used to develop a workflow type procedure by examining in detail a few elements of the algorithm described in the previous section. The synthetic model is from the Noranda District in Canada, home to 20 economic volcanogenic massive sulphide deposits (VMS), nineteen orogenic Au 68 Algorithm 2 Large-Scale Inversion mc\u00E2\u0086\u0090 argminmc(\u00CF\u0086d(mc)+\u00CE\u00B2\u00CF\u0086m(mc)) d\u00CC\u0083t=1:T = Pdt dobs m1 = Lfcmc u1\u00E2\u0086\u0090 F(m1) Initialize: \u00CE\u00B2\u00CC\u0083 for k = 1,2, ... do for t = 1...T (in parallel) do Set primary model: m\u00CC\u0083kt p = Pmt mk u\u00CC\u0083kt p = Put uk Invert each tile: Initialize: m\u00CC\u0083kt = Pmt mk Solve: m\u00CC\u0083kt \u00E2\u0086\u0090 argminm\u00CC\u0083kt (\u00CF\u0086d(m\u00CC\u0083kt ) + \u00CE\u00B2\u00CC\u0083 \u00CF\u0086m(m\u00CC\u0083kt )) where: B\u00CC\u0083(m\u00CC\u0083kt )u\u00CC\u0083kts = q(S\u00CC\u0083kt \u00E2\u0088\u0092 S\u00CC\u0083ktp)Qu\u00CC\u0083ktp from equation 3.7 end for Merge models: mk+1 =M(m\u00CC\u0083kt=1:T) Compute fields: uk+1\u00E2\u0086\u0090 F(mk+1) Check misfit: if \u00CF\u0086d(mk+1)<= \u00CF\u0086 \u00E2\u0088\u0097d then break end if \u00CE\u00B2\u00CC\u0083 k+1 = \u00CE\u00B3\u00CE\u00B2\u00CC\u0083 k, where 0 < \u00CE\u00B3 < 1 end for 69 deposits, and several intrusion-hosted Cu-Mo deposits (Gibson and Galley, 2007). The original model, provided courtesy of the Xstrata mining group, contained 12.7 million cells covering an area of almost 20 km x 20 km. The 38 geologic units in the model were converted into expected conductivities and the entire 12.7 million cell model was forward modelled at 30, 45, 90, 180, 360 and 720 Hz. The data were then corrupted with noise equal to 2 percent of the data plus a floor of 0.005 to form the observed synthetic data. Two slices at -275 m and -475 m of the synthetic model can be seen in figures 3.2 and 3.3. 3.3.1 Preliminary Coarse Inversion The first step to find an approximate solution to the full inverse problem is to per- form a coarse inversion to determine the large-scale conductivity structure that will be used to compute the primary fields and starting models for subsequent finer cell size inversions. The earth is discretized into relatively large cells such that the total number of cells in the mesh is small and the inverse problem can be solved quickly. In this example the coarse mesh contained 42 x 42 x 73 cells. The cell lengths in the x and y directions were 500 m. Working initially on a coarse mesh, which for this synthetic example can be inverted in approximately 30 minutes, allows multiple inversions to be run simultaneous with different parameters and starting models. The coarse inversion result can be seen in figure 3.4. Because of the large cell dimensions, some geologic structures such as accurate body boundaries 70 Easting (m) N or th in g (m ) S/m -8500 -4250 0 4250 8500 - 8500 - 4250 0 4250 8500 Figure 3.2: Depth slice at -275 m of the true model. The topography for the synthetic model ranges from -25 m to -250 m. The conductivity of the ore was chosen to be 0.2 S/m; however, the color scale on the model has been clipped at 0.01 S/m to improve the visualization of the other geologic units. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) and fine scale features may not be recovered by the discretization. Therefore it is important not to overfit the data and risk adding discretization artifacts into the inversion result. This is particularly true for higher frequency data which will have smaller skin depths and contain more information about fine scale features. In fact, because the initial goal of the coarse scale inversion is to quickly determine 71 Easting (m) N or th in g (m ) S/m -8500 -4250 0 4250 8500 - 8500 - 4250 0 4250 8500 Figure 3.3: Depth slice at -475 m of the true model. The conductivity of the ore was chosen to be 0.2 S/m; however, the color scale on the model has been clipped at 0.01 S/m to improve the visualization of the other geologic units. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) the large-scale conductivity features, some data, particularly the higher frequency data, may be omitted to reduce the computational cost, and prevent the large cell sizes of the coarse mesh from violating the shorter skin depths at the higher fre- quencies. In this example the 360 Hz and 720 Hz data were omitted for the coarse mesh inversion. 72 Easting (m) N or th in g (m ) S/m -8500 -4250 0 4250 8500 - 8500 - 4250 0 4250 8500 Figure 3.4: Coarse inversion result at -475 m shown on the same color scale as the true model. Even with coarse 500 m x 500 m cells in the x and y directions, the main geologic features of the model have been recovered. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) 3.3.2 Interpolating Models and Computing Primary Fields In order to compute the secondary fields on each subdomain tile, the primary fields are needed according to equation 3.6. These computed primary fields contain in- formation about the large-scale structure, and are first obtained by linearly inter- polating the coarse inversion model onto the fine model domain. Linear averaging when interpolating from a coarse mesh to a fine mesh will create a smoother fine 73 model which is consistent with our smooth model regularization in equation 3.1. The new interpolated fine model is then forward modelled to compute the primary electric fields which are used (equation 3.6) to calculate the total fields on each subdomain mesh. While forward modelling operations on the full fine mesh can be time consuming and resource intensive, this operation only needs to be performed once for each frequency and outer tile iteration. Our experience shows that working with primary fields computed using the interpolated fine model is superior to the fields directly from the coarse mesh, and the additional cost of forward modelling on the fine mesh is worthwhile. 3.3.3 Setting Up the Subdomains The coarse inversion result can be used to generate an initial starting model and to compute the primary fields for subsequent finer discretized inversions. Once the initial models and primary fields have been computed, the subdomain meshes must be designed. When decomposing the computational domain and designing the sub- domain meshes there are two critical elements that must be considered. Firstly, the decomposition should satisfy the physics between domain interactions. Secondly, in the final inverted model there should be continuity of conductivity structures across the various domain boundaries. Incorporating domain interactions is done by using primary fields in the tiled inversions which have been computed from full domain conductivity structures. Continuity of the conductivity structure across 74 domain boundaries is accomplished by using both data and mesh overlap. Gener- ally the core regions of the meshes are chosen to overlap by slightly less than the skin depth of the lower frequencies (additional padding cells are added outside this overlap region). Here a tradeoff must be achieved between increased overlap and increased computational requirements. Another tradeoff that must be considered is the number of tiles used to decompose the domain. As the number of tiles in- creases, while the individual inversion time for each subproblem might decrease, the percentage of overlap cells vs total number of cells in each mesh increases and the decomposition may become less efficient. Here it is up to the user to choose a reasonable number of tiles such that each tile runs efficiently on their available computer hardware. For this synthetic example, two different sets of meshes, one with three tiles and the other with nine tiles are used to demonstrate the scalability and robustness of the methodology. For the three tile example, each subdomain mesh contained approximately 3.5 million cells with 353 x 140 x 70 cells in the easting, northing and vertical directions. Each mesh in the nine tile subdomain example contained approximately 1.4 million cells with 139 x 140 x 70 cells in the easting, northing and vertical directions. Both the three and nine tile examples contained the same sized cells (50 m x 50 m cells in the easting and northing directions) in the core and overlapping regions. In the vertical direction the cells dimensions started at 25 75 21 3 Figure 3.5: Mesh layouts for the example using three overlapping domains. Each subdomain mesh contained approximately 3.5 million cells, with 353 x 140 x 70 cells in the easting, northing and vertical directions. The central core region of each domain overlapped by 20 cells, each with dimensions of 50 m, in both the x and y directions. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) m, and then expanded with depth. Both sets of meshes used 20 overlapping cells on each mesh boundary. The mesh layouts and overlaps for the three and nine tile examples can be seen in figures 3.5 and 3.6. 3.3.4 Perform Subdomain Inversions Inverting large field datasets without some sort of decomposition method is not feasible. Selecting the optimal tile size is a complex process which depends on 76 45 6 7 8 9 1 2 3 Figure 3.6: Mesh layouts for the example using nine overlapping domains. Each mesh contained approximately 1.4 million cells with 139 x 140 x 70 cells in the easting, northing and vertical directions. The central core region of each domain overlapped by 20 cells, each with dimensions of 50 m, in both the easting and northing directions. ( c\u00C2\u00A9Society of Explo- ration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) many factors such as system CPU, RAM, and inversion parameters. As previously stated, as the tile size becomes smaller, while the individual inversion time for each subproblem might decrease, the percentage of overlap cells vs total number of cells in each mesh increases and the decomposition may become less efficient. The user must decide between working with smaller tiles with quicker inversion times and increased parallelism (since all of the tiles can be run simultaneously in parallel) and larger tiles that may take longer to run, but have less overlap cells vs total num- 77 ber of cells, and the full problem has less overlapping tile boundaries. Generally our methodology is to select the largest tile size that can be inverted in an accept- able time. For the synthetic example, both the three tile and nine tile examples were inverted on Intel Xeon X5660 processors each with 64 GB of random access memory (RAM). For the nine tile example there was sufficient memory on each computer to run the six frequencies in parallel. Each beta iteration took around 26 hours. For the three tile example, only three frequencies could be run in parallel on each computer with the 64 GB of available memory. Here each beta iteration took around 104 hours. During the subdomain inversions the initial beta is chosen to be the final beta from the initial coarse inversion result, and the same cooling scheme is used on the tiled inversions as in the initial coarse inversion. 3.3.5 Reference Station Considerations An additional complication when tiling ZTEM inversions, is the reference station that is used to normalize the data. While the location of the reference station may be included in full dataset inversions (in the full dataset inversion the reference sta- tion is positioned near the edge of the data grid), when the full modelling domain is subdivided, the true reference station location cannot practically be included in each domain. To circumvent this problem, the horizontal field values are fixed to those obtained from the initial coarse inversion, which should be a good approxi- mation to the true fields. 78 3.3.6 Model Updates and Merging Tiled Inversions Once the subdomain inversions have been completed, a full model update can be obtained by merging the inverted tiles. For updating model cells that correspond to only one inversion tile, it is clear what conductivity value should be assigned to the update model. The situation becomes more complicated for update cells which correspond to multiple subdomain cells (overlapping tile regions). A simple merging technique that will mitigate discontinuities across tile boundaries is to use an averaging scheme that assigns the conductivity of the update cell to be the weighted average of the overlapping conductivity tiles. The weights for each cell can be based on the distance from that cell to the respective tile boundary. From our practical experience, this method works quite well; however, we can improve on this method by trying to incorporate some estimate of each cell\u00E2\u0080\u0099s sensitivity to the individual tiles. Firstly, it is assumed that the sensitivity Ji j of the jth cell to the ith datum has a geometric decay that drops off as 1r2 according to Biot-Savart\u00E2\u0080\u0099s Law. Secondly, it is assumed that there is an attenuation component that drops off as exp(\u00E2\u0088\u0092r/\u00CE\u00B4 ) where \u00CE\u00B4 is some approximation to the skin depth. Although clearly the skin depth will vary with frequency and model conductivities, to first order we can gain some insight into the relative sensitivity of a particular cell to each data by using a reasonable fixed skin depth. That is we may approximate the sensitivity 79 to be proportional to, Ji j = K r2i j exp( \u00E2\u0088\u0092ri j \u00CE\u00B4 ), (3.10) where K is the constant of proportionality. For cells that fall in an overlapping region, we can determine weighting schemes based on the relative sensitivity of a particular cell to each tile of interest. For some cell j, the normalized sensitivities for each tile (using a simple example of two overlapping tiles) can be given as a sum over the N1 and N2 data points in each tile (data points that overlap data points in other tiles are omitted in the calculation). W 1j = 1 N1 N1 \u00E2\u0088\u0091 i=1 J(1)i j (3.11) and, W 2j = 1 N2 N2 \u00E2\u0088\u0091 i=1 J(2)i j . (3.12) Here we have assumed that each data should be weighted equally. The conductivity for the jth cell, \u00CF\u0083\u00CC\u0084 j, can then be given as a weighted average of the overlapping inverted tile conductivities \u00CF\u00831 and \u00CF\u00832. \u00CF\u0083\u00CC\u0084 j = 1 W 1 +W 2 (W 1\u00CF\u00831j +W 2\u00CF\u00832j ) (3.13) 80 This method is easily extended to more complex tile overlaps. The model updates for the synthetic and field data examples were computed using this update scheme. The final inversion models from the three and nine tile examples in figures 3.7 - 3.10, recover the large-scale geological features as well as several of the larger known mineralized regions which were not evident from the coarse inversion. The inversion results of the three and nine tile examples are very similar, illustrating that the inversion results are not highly dependent on the choice of subdomain meshes. Figure 3.11 compares the recovered high conductivity regions (mineralized zones) with the true mineralized zones. While one cannot expect to recover the exact complex geometry of these relatively small mineralized zones, the inverted model does a good job of recovering the general geometry of the ore bodies. 3.4 Pebble Field Data Example The Pebble deposit is a world-class calc-alkalic copper-gold molybdenum por- phyry deposit located in the Bristol Bay region of southwest Alaska (geographic location shown in figure 3.12). The deposit is managed by the Pebble Partnership, a joint venture between Northern Dynasty Mines and Anglo American plc. Be- cause of the size of the deposit and the quality and diversity of geophysical data collected over the area, this deposit is an ideal test site for geophysical techniques and inversion methodologies. Over the years, many groups have looked at the geo- logic and geophysical data from the area, but in particular Pare and Legault (2010) 81 Easting (m) N or th in g (m ) S/m -8500 -4250 0 4250 8500 - 8500 - 4250 0 4250 8500 Figure 3.7: Depth slice at -275 m of the inversion result after merging the three tiles together, shown on the same color scale as the true model. The resolution in this model is greatly improved compared to the coarse result. Several of the larger known deposits have now been recovered; they were not visible from the initial coarse inversion. The recovered mineralized body in the inversion between easting -4250 m and 0 m is a mineralized body in the synthetic model that only extends to -250 m. The extend of this body that is not apparent in the -275 m depth slice of the true model, can be seen in figure 3.11. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Olden- burg (2012)) 82 Easting (m) N or th in g (m ) S/m -8500 -4250 0 4250 8500 - 8500 - 4250 0 4250 8500 Figure 3.8: Depth slice at -275 m of the inversion result after merging the nine tiles together, shown on the same color scale as the true model. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) provide a fairly comprehensive overview of the general geology and the available geophysical datasets. A regional geologic map for the area can be seen in fig- ure 3.13 and a more local scale geologic map of the area can be seen in figure 3.14 In summary, from Pare and Legault (2010), the Pebble Deposit consists of two contiguous zones, the Pebble West and East Zones. The Pebble West Zone oc- curs around small granodioritic stocks that intrude the country rock, and the Pebble 83 Easting (m) N or th in g (m ) S/m -8500 -4250 0 4250 8500 - 8500 - 4250 0 4250 8500 Figure 3.9: Depth slice at -475 m of the inversion result after merging the three tiles together, shown on the same color scale as the true model. East Zone occurs within a granodioritic stock and in sills that cut the country rocks. Because the Pebble West Zone extends to surface, the deposit was discovered in 1986, far before its deeper and higher grade neighbor which was not discovered un- til 2005. Pebble East is completely hidden by up to 600 m of volcano-sedimentary cover which explains why it took almost two decades to discover, and highlights the importance to develop deep probing geophysical techniques. More recently, the Far East Zone has been discovered at a depth of greater than 1.5 km. The deposit 84 Easting (m) N or th in g (m ) S/m -8500 -4250 0 4250 8500 - 8500 - 4250 0 4250 8500 Figure 3.10: Depth slice at -475 m of the inversion result after merging the nine tiles together, shown on the same color scale as the true model. hosts K-silicate alteration associated with quartz-sulphide veins, overprinted by phyllo-silicate alteration. The higher grade mineralization at Pebble East is associ- ated with advanced argillic alteration. From Rebagliati et al. (2009), the measured and indicated resource of both zones is given at a 0.3% equivalent copper cut-off grade of 5.096 billon tonnes at 0.43% Cu, 0.35g/t Au, and 256 ppm Mo. A plan view and cross section of the Pebble West and Pebble East zones, can be seen in figures 3.15 and 3.16. 85 a) b) -6600 -8400 -7500 -6600 -8400 -7500 Easting (m ) -2325 3825 Easting (m ) -2325 Northin g (m)3825 -900 -125 -125 -900 Northin g (m) Z (m) Z (m) Figure 3.11: Comparison of the larger southern mineralized regions. a) Re- covered bodies more conductive than 0.005 S/m that correspond to the mineralized bodies in the southern section of the inverted model. b) Actual mineralized bodies. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) 86 Figure 3.12: Location of the Pebble deposit in the Bristol Bay region of Alaska. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) Over the years a wealth of geophysical data has been collected over the Pebble deposits. Along with potential field data, multiple DC-IP surveys have been col- lected along with two small MT surveys. More recently in the summer of 2009, a 3840 line-km Spectrem survey (Leggatt et al., 2000) was performed. In the same summer, a small 250 line-km ZTEM survey was flown over the Pebble deposit. In the summer of 2010, a much larger 6504 line-km ZTEM survey was flown over 87 Figure 3.13: Regional geology after Lang et al. (2009). The Pebble deposit sits in the Lake Clark Fault Zone. Figure 3.14: Map of the local geology and timeframe of the major geologic events near the Pebble deposit after Lang et al. (2009) 88 Pebble West Pebble East N Jura-Cretaceous to Eocene Igneous and Sedimentary Rocks Late- Cretaceous to Early Tertiary Volcanic and Sedimentary Rocks Figure 3.15: Plan view of the Pebble Deposit showing the deposit outline and grade distribution across the deposit. ( c\u00C2\u00A9Society of Exploration Geo- physicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) the deposit and surrounding areas. The 2010 survey used 250 m line spacing and collected 30,45, 90, 180, 360, 720 Hz data. In this section of the chapter, we apply our large-scale inversion methodology outlined in the previous sections to invert a large 30 km x 30 km section of the newer 6504 line-km ZTEM survey. The results are compared to deposit scale 2D DC resistivity inversions and the regional scale time constant analysis from the Spectrem AEM data. A central (30 km x 30 km) block of the 2010 ZTEM survey was chosen to be inverted in 3D. The ZTEM survey flight path geometry and the inverted subset of 89 Pebble West Pebble East Stock A East Zone Stock Jura-Cretaceous to Eocene Igneous and Sedimentary Rocks Figure 3.16: Cross section of the Pebble Deposit showing the grade dis- tributions and the intrusive stocks. ( c\u00C2\u00A9Society of Exploration Geo- physicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) the full survey is shown in figure 3.17. Initially the lower frequency (30, 45, 90 and 180 Hz) data were inverted on a coarse 83 x 83 x 52 mesh. The cell dimensions in the central core region were 400 m x 400 m in the x and y directions and started at 25 m in the vertical direction before expanding with depth. While the 400 m cell sizes in the horizontal dimensions may not be adequate to recover the fine deposit scale features, the coarse discretization allowed the full 30 km x 30 km block to be inverted in around 2 hours on a single dual Xeon X5660 processor. The coarse inversion result at an elevation of 100 m can be seen in figure 3.18. Next, the modelling domain was split into three tiles which were then run in parallel. Each tile contained 314, 112, and 52 cells in the x,y and z directions. Each cell in the coarse inversion was subdivided into sixteen cells by using 100 m x 100 m 90 \u0001\u0002\u0003\u0004\u0005\u0006\u0004 \u0007\u0008 \u0006 \u0004 \u0001 \u0008\u000E Figure 3.17: 6504 line-km ZTEM survey block flown in the summer of 2010. The 30 km x 30 km block inverted in this section is outlined by the dashed square. The full survey block is shown as the solid black re- gion. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by per- mission from Holtham and Oldenburg (2012)) horizontal cells in the tiled inversion. The tiles took around 36 hours to run on dual Xeon X5660 processors. The tiles were merged using the approximate sensitivity procedure outlined in the synthetic inversion section of this chapter. Figure 3.19 shows the fine scale model inverted by merging the three tiles together. The regional scale inversion result shows good agreement with the Spec- trem AEM Z component time constant tau plot shown in figure 3.20. On the deposit scale, figure 3.21 shows a comparison between the ZTEM inversion re- sult and a 3D DC inversion performed using a UBC OcTree DC/IP inversion code. To first order there seems to be good agreement on the regional scale between the Spectrem and ZTEM and on the deposit scale between the ZTEM data and the DC 91 NS/m 30 km 0.1 0.013 0.00168 2.17e-004 Figure 3.18: Coarse ZTEM inversion result at a depth of 100 m. The edge of the data coverage is shown by the dashed lines and the deposit is outlined by the solid lines. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) data. 92 NS/m 30 km 0.25 0.015 8.5e-004 5e-005 Figure 3.19: Fine ZTEM inversion result at an elevation of 100 m formed by merging the three subproblem tiles together. The edge of the data coverage is shown by the dashed lines and the deposit is outlined by the solid lines. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) 93 \u00C2\u00B5s 30 km N Figure 3.20: Spectrem Z component time constant plot over the same region as the inverted ZTEM data. The outline of the deposit is outlined in black. ( c\u00C2\u00A9Society of Exploration Geophysicists, 2012, adapted by per- mission from Holtham and Oldenburg (2012)) 94 23 km S/m N S/m 23 km N a) b) 0.25 0.02 0.0016 1.25e-004 0.25 0.02 0.0016 1.25e-004 N Figure 3.21: Inversion results shown on the same color scale and elevation of 100 m. The deposit is outlined in black. a) 3D DC Inversion b) ZTEM inversion shown over the same area as the DC resistivity data coverage. To first order the results seem consistent. ( c\u00C2\u00A9Society of Ex- ploration Geophysicists, 2012, adapted by permission from Holtham and Oldenburg (2012)) 95 Chapter 4 Joint Inversion of ZTEM and MT Data 4.1 Introduction Natural source electromagnetics play an important role in understanding the elec- trical conductivity of the upper regions of the earth. Their deeper depth of in- vestigation compared to controlled source electromagnetic (CSEM) methods, is a consequence of the plane wave excitation. Furthermore, computing the data at each frequency only requires the solution of a matrix system with two right hand sides which can result in significant computational advantages over multi- transmitter controlled source electromagnetics. Unfortunately, because of the ex- 96 pense of ground based surveys, the traditional natural source MT method can be impractical as a stand alone method for large surveys where high spatial resolution is required. While an airborne natural source system that measures both the electric and magnetic fields does not currently exist, the airborne AFMAG (Audio Frequency Magnetics) (Ward, 1959) and recently the ZTEM technique (Lo and Zang, 2008) have been developed that measure magnetic fields due to natural sources. For the ZTEM system, the data (Tzx and Tzy) are frequency domain relations that relate the vertical magnetic fields over the survey area to the horizontal magnetic fields at the ground based reference station. The typical 30-720 Hz collected data makes the method effective for finding targets at moderate depths of up to a few kilometers for resistive geologies; however, the depth of investigation is significantly shallower in more conductive regions. Although the ZTEM system may have a greater depth of investigation than many other airborne electromagnetic methods, the depth of investigation is still far below that possible using a ground MT survey which can collect much lower frequency data. For deep investigation, particularly in conduc- tive environments, the ZTEM depth of penetration may not be adequate. Also, since 1D layered structures produce no vertical magnetic fields, ZTEM data alone cannot recover purely 1D structures and it is desirable to obtain robust conductiv- ity information from other sources and reduce the non-uniqueness of the inverse 97 problem. This information could come from a priori geologic and petrophysical information or from additional geophysical data such as MT data. The magnetotelluric technique (Tikhonov, 1950; Cagniard, 1953) measures electric and magnetic fields at the surface of the earth using MT sites (see fig- ures 4.1 and 4.2). The data are impedance tensors Z, formed by taking the ratios of the electric fields Ex and Ey and the magnetic fields Hx and Hy. This relation is given by, \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD Ex Ey \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8= \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD Zxx Zxy Zyx Zyy \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 \u00EF\u00A3\u00AB\u00EF\u00A3\u00AC\u00EF\u00A3\u00AC\u00EF\u00A3\u00AD Hx Hy. \u00EF\u00A3\u00B6\u00EF\u00A3\u00B7\u00EF\u00A3\u00B7\u00EF\u00A3\u00B8 , (4.1) where Zxx, Zxy, Zyx and Zyy are the four components of the MT impedance ten- sor. Because the MT technique is often used to image at depth over large areas, significant station spacing is often employed to keep the survey costs affordable. Since each geophysical exploration technique has different practical strengths and limitations, combining multiple geophysical surveys can provide complementary information that will reduce the non-uniqueness of the inverse problem and recover an improved inversion result. An effective method that has sensitivity at depth and high resolution at moderate depths is to jointly invert CSEM and MT data as done for example by Mackie et al. (2007) and Commer and Newman (2009). Although it is possible to invert multi-transmitter electromagnetic data (for instance Olden- burg et al., 2008), it is time-consuming and computationally intensive. With the 98 Figure 4.1: Typical ground MT site from Smirnov et al. (2008). Two sets of electrodes measure the electric field in two orthogonal directions. The two horizontal magnetic fields are measured (possibly the vertical component as well). Because the MT site may be left for a long period of time, additional devices such as power sources, data loggers etc. may also be required. advent of the ZTEM system, it is now possible to collect and invert ZTEM and MT data simultaneously. The ZTEM data are used to gather high spacial informa- tion at moderate depths, while sparse MT data are used to recover deep structures. Since the forward modelling operation required to compute the ZTEM and MT responses is the same, this procedure is computationally efficient and much easier than inverting CSEM and MT data. In this paper we show how ZTEM and MT data can be jointly inverted in three dimensions to provide an effective method to image the earth. 99 Figure 4.2: Seafloor MT site. Each site measures the electric and magnetic fields, and after sufficient data has been recorded the MT site is re- leased from the bottom of the ocean and collected by the control ship. ( c\u00C2\u00A9Society of Exploration Geophysicists, 1998, adapted by permission from Constable et al. (1998)) 4.2 Methods Our MT inversion algorithm is that of Farquharson et al. (2002) and our ZTEM inversion algorithm is that outlined in chapter 2. Our algorithm to invert ZTEM and MT data simultaneous was obtained by modifying these codes. The solution is obtained by an iterative Gauss-Newton procedure based on Haber et al. (2000a). By minimizing the objective function \u00CE\u00A6= \u00CF\u0086d +\u00CE\u00B2 (\u00CE\u00B1s||Ws(m\u00E2\u0088\u0092mre f )||22 +\u00CE\u00B1x||Wx(m\u00E2\u0088\u0092mre f )||22 +\u00CE\u00B1y||Wy(m\u00E2\u0088\u0092mre f )||22 +\u00CE\u00B1z||Wz(m\u00E2\u0088\u0092mre f )||22) we obtain our solution to the inverse problem. \u00CF\u0086d is the data-misfit, m is the in- verted model, mre f is a reference model, Ws is a diagonal matrix, Wx, Wy, and 100 Wz are the first-order finite-difference matrices in the x, y and z directions. The \u00CE\u00B1 \u00E2\u0080\u00B2s are adjustable parameters. \u00CE\u00B1s controls the closeness of the recovered model to the reference, while \u00CE\u00B1x,y,z determine the smoothing in the x, y and z directions. \u00CE\u00B2 is the regularization parameter which is reduced throughout the inversion process until the desired data misfit, \u00CF\u0086 \u00E2\u0088\u0097d , has been achieved. Although the forward problems for ZTEM and MT are extremely similar, when considering the data there are dif- ferences that can have important implications for the inversion process. Because ZTEM data are collected from a helicopter, while MT data are collected on the ground, the data errors can be different. When inverting multiple datasets simul- taneously, a tradeoff must be achieved between the different data misfit terms and the model regularization in the objective function, \u00CE\u00A6= \u00CF\u0086ZTEMd +\u00CF\u0086 MT d +\u00CE\u00B2\u00CF\u0086m, (4.2) where \u00CF\u0086ZTEMd is the ZTEM data misfit, \u00CF\u0086 MT d is the MT data misfit and \u00CF\u0086m is the model regularization. Because typical surveys collect many times more ZTEM data than MT data, if no additional data re-weighting is imposed, then the inversion can easily add shallower structure to the model that fits the ZTEM data, and not fit the few low frequency MT data points without greatly effecting the total data misfit. While selectively fitting the ZTEM data over the MT data may achieve the 101 final target misfit, this method may not extract the desired information from both datasets since the ZTEM data may be selectively fit over the lower frequency MT data. Our approach is that ZTEM and MT data should play approximately equal roles in the inversion and hence the two data misfit terms, \u00CF\u0086ZTEMd and \u00CF\u0086 MT d should be approximately equal. To accomplish this, we have incorporated an additional weighting term in our data misfit. The data misfit is then computed via, \u00CF\u0086d = ||Wd1(dobsZTEM\u00E2\u0088\u0092dprdZTEM)||22 + \u00CE\u00B3||Wd2(dobsMT \u00E2\u0088\u0092dprdMT )||22, (4.3) where \u00CE\u00B3 is a constant that must be determined. dobsZTEM and dobsMT are the observation vectors, dprdZTEM and d prd MT are the predicted data, and Wd1 and Wd2 are diagonal matrices whose elements are the reciprocals of the standard deviations of the datum errors. Assuming a gaussian noise model for each datum, the target misfits for the ZTEM and MT datasets would respectively be the number of ZTEM data, N1, and the number of MT data, N2. In order to make the contributions to the final misfit equal for each datatype, the constant \u00CE\u00B3 is then equal to N1N2 . This approach to determine the relative weighting is the same as the joint CSEM and MT inversion problem of Mackie et al. (2007) and Commer and Newman (2009) where the re- weighting takes into account the number of data from each survey. Initially the ZTEM and MT data are inverted individually to determine appropriate standard 102 deviations for each datum, before being inverted together with the appropriate re- weight \u00CE\u00B3 . 4.2.1 Determining a Background Model Determining a reasonable starting model is critical for any inversion. Because of the highly non-unique nature of the geophysical inverse problem, inversions started with different initial models can result in variations in the final model even if the data are sufficiently fit. Since 1D layered structures produce no vertical mag- netic field, traditional methods that determine a best fitting half-space conductivity structure cannot be used for ZTEM data with a flat earth. In order to obtain robust information about the background conductivity structure, additional information is needed. If only ZTEM data are being collected, even a few electric field measure- ments can provide additional information about the 1D background structure. In this thesis where both ZTEM and MT surveys are considered, our solution is to use the MT data to determine the background and initial models. 4.2.2 Synthetic Example We demonstrate the effectiveness of jointly working with ZTEM and MT data by examining the synthetic model in figure 4.3. The model has a central resistive core, an outer region of high conductivity, all intruded in a background host. The earth was discretized into a mesh containing 62, 62, and 45 cells in the x, y, and z 103 Easting (m) Northing (m) Z (m) Conductivity (S/m) -500 to 500 -1000 to 1000 -250 to -3500 0.001 -1500 to 1500 -2000 to 2000 -250 to -3500 0.1 Easting (m) E le va tio n (m ) S/m Figure 4.3: Cross section through northing = 0 m of the synthetic conduc- tivity structure of a resistor surrounded, except on the top and bottom, by a conductor. Both blocks are embedded in a 0.01 S/m background half-space. directions respectively. ZTEM data (30, 45, 90, 180, 360 Hz) and MT data (0.04, 0.1, 0.4, 1, 3.5, 7, 30, 45, 90, 180, 360 Hz) were computed from the synthetic model using three different survey configurations over a 9 km x 10 km area whose details are shown in figure 4.4. The first survey is an extremely coarse MT survey using 9 stations in total with 1.5 km station spacing and 2 km line spacing. This station spacing will not be adequate to achieve the high resolution required for many detailed exploration projects. The second survey configuration is a finer MT survey using 1196 stations in total with 200 m station spacing and 400 m line spacing. Having almost 1200 stations over such a small area will be too expensive 104 for most exploration programs considering a typical MT station can cost greater than a thousand dollars. The ZTEM survey has 10 m data spacing and 300 m line spacing. Gaussian noise was added to all the data before being inverted. For the ZTEM data, the standard deviations of the noise were assigned based on a floor of 0.005 plus 2 percent of the data. The for the MT data, the standard deviations on the diagonal elements were assigned as 3.5 percent of the data plus a floor of 0.005. For the off diagonal elements the standard deviations were assigned as 1 percent of the data plus a floor of 0.0001. To demonstrate the advantages of working with ZTEM and sparse MT, the data from the ZTEM and coarse MT grids were jointly inverted. For all of the inversions, the initial and reference models were set to be 0.01 S/m, and the predicted data fit the observed synthetic data very well. The final conductivity models can be seen in figures 4.5b - 4.5e. The coarse MT inversion result (figure 4.5b) recovers a crude estimate for the outline of the conductive block. Because of the limited data coverage from the 1.5 km x 2 km MT station spacing, the boundaries of the blocks are not very well resolved. Nevertheless, the 9 MT stations are enough to resolve the approximate outline of the conductive block. The single MT station over the resistive struc- ture does manage to recover some estimate of the restive core. As expected, by increasing the number of MT stations, the boundaries of the blocks become bet- ter resolved (figure 4.5c) at the expense of drastically increasing the survey costs. 105 N or th in g (m ) Easting (m) Figure 4.4: Survey grids for the synthetic inversions. The three grids are the coarse MT (large squares), fine MT (small squares), and ZTEM (solid lines). The ZTEM inversion (figure 4.5d) shows good resolution down to around 1 km which is impressive for an airborne method. By jointly inverting the ZTEM and coarse MT grid (figure 4.5e), the inversion recovers both the deep structure from the low frequency MT data and the fine structure from the ZTEM data. The 1.5 km x 2 km station spacing, which was inadequate to recover the near-surface details, greatly enhances the ZTEM inversion and adds significant information about the large-scale features at the 2-3 km depth range. 106 Easting (m) El ev at io n (m ) S/m (a) True cross section of the model Easting (m) E le va tio n (m ) S/m (b) Inversion of the coarse MT grid S/m E le va tio n (m ) Easting (m) (c) Inversion of the fine MT grid Easting (m) E le va tio n (m ) S/m (d) Inversion of the ZTEM grid E le va tio n (m ) Easting (m) S/m (e) ZTEM and coarse MT grid Figure 4.5: Inversion of the synthetic model in figure 4.3, all shown on the true color scale. Adding the coarse MT data to the ZTEM data greatly improved the inversion result in the 2-3 km depth range. 107 4.3 Field Data The joint inversion of ZTEM and MT data is a viable technique to search for geothermal exploration targets (among other potential applications) because it may have the capability to penetrate the high conductivity environments often associ- ated with geothermal reservoirs (Pellerin, 2002). We demonstrate the joint inver- sion of the two data types by inverting overlapping ZTEM and MT data from the Reese River geothermal site in Nevada. 4.3.1 Reese River Background The Reese River property is a blind geothermal resource located in north central Nevada (27 miles north of Austin and 200 miles northeast of Reno, see figure 4.6). The concept of a blind geothermal resource implies that there is no surface expres- sion of the resource such as hot springs, and hence the property must be found through geophysical methods and drilling. The property lies within the Humboldt structural zone, a region of NE-trending faults that may be structural controls for geothermal reservoirs (Faulds et al., 2003). The main feature in the property is the NE-trending Range Front Fault that runs through the property. The surface geology can be seen in figure 4.7. Over the years, significant geophysical and geologic mapping has been per- formed on the property, the details of which are outlined in Witter et al. (2009). 108 Figure 4.6: Reese River geothermal site located in Nevada. On the eastern portion of the survey area is the Range Front Mountains. The 7 MT lines are shown in red. The map of the location of the Reese River site in Nevada is provided curtesy of Sierra Geothermal Power Corp. Briefly in summary, the geothermal reservoir was first mapped in the 1970\u00E2\u0080\u0099s to 1980\u00E2\u0080\u0099s using shallow temperature gradient wells after uranium prospectors discov- ered hot water in their drill holes. In 2007, a 2D reflection and refraction seismic survey was performed, and in 2008, ground based gravity data was collected. In the same year, a MT survey was performed then finally a ZTEM survey was flown in 2009. In addition to the geophysical data collected over the prospect area, ge- ologic mapping, geological logging of drill cuttings from 5 deep wells (between 487 m to 1580m in depth), chemical geothermometry of well waters, and shallow 109 North 88-5 56-4 13-4 18-33 62-4 Upflow Zone?Outflow? 2 km Figure 4.7: Reese River surface geology provided from Sierra Geothermal Power Corp. The 7 MT lines are shown in red along with the location of the 5 drill holes. A possible interpretation for the fluid flow upwelling and outflow is also outlined. The location of the Range Front Fault is outlined by the thick black line. (2 m) temperature studies have also been performed. As outlined in Witter et al. (2009), Sierra Geothermal Power (SGP) drilled four temperature gradient wells and one slim well at the Reese River prospect. The slim well had a total drilled depth of 1200 m and the temperature gradient wells had depths of 487 m, 495 m, 610 m and 1580 m. Drill cuttings from all of the wells were logged. All the well logs shown in figure 4.8 show intersections with tertiary sedimentary rocks in the top few hundred meters which one would expect to be quite conductive. At greater depths, the wells intersected a combination of shales, siltstones and other generally assumed conductive rocks. The two deepest wells, 110 13-4 56-4 18-33 88-5 62-4 Volcanic Tuff Paleozoic siltstone TD = 1600 feet TD = 1620 feet TD = 2000 feet TD = 3930 feet TD = 5200 feet Granite Tertiary Sedi- mentary Rocks Volcanic Tuff Dolomite (Vinini Fmtn) Siltstone, Sandstone, & Shale (Vinini Fmtn) Si lty , Pa le o zo ic Li m es to n e & Sh al e Tertiary Sedimentary Rocks Volcanic Tuff Tertiary Sedimentary Rocks Shale (Vinini Fmtn) Pa le o zo ic Li m es to n e Tertiary Sedimentary Rocks Quartzite, siltite, argillite (Vinini Fmtn) T er tia ry Se di m en ta ry R o ck s Paleozoic siltstone (Vinini Fmtn) Figure 4.8: Reese River well logs. From Palacky (1987), one would expect a conductive structure in the near surface due to the conductive sediments, and a more resistive structure at depth because of the limestone and granite. 13-4 (1580 m depth) and 56-4 (1200 m depth), intersect limestone at depth, which when unaltered, should be more resistive than the shallower rocks. The deepest well 13-4, intersected granite at depths greater than 1500 m. 4.3.2 ZTEM and MT Inversion Results In the winter of 2008, 7 lines of ground based MT data was collected by Quantec Geoscience. The 30 line-km of Quantec Spartan MT data was collected with the lines oriented perpendicular to the known NE trending structure. The survey used 111 400 m site separation, and the electric fields were measured using a dipole spacing of 200 m in the x direction and 100 m in the y direction. Following up on the ground MT data, an airborne ZTEM survey was flown by Geotech in August 2009. The survey collected data at 30, 45, 90, 180, 360, and 720 Hz, and was flown over a larger area than the MT survey. The survey consisted of 7 lines each of which was approximately 10 km in length. For the inversion, the earth was discretized on a mesh that contained 121, 88 and 70 cells in the x, y and z directions respectively. The topography was obtained from the U.S. Geological Survey (USGS) National Elevation Dataset (NED). The cell dimensions in the central core region were 100 m x 100 m x 30 m. Small cells were used in the vertical dimension to limit the modelling errors due to inaccuracies in discretizing the topography. The ZTEM data (30, 45, 90, 180, 360 Hz) and the MT data (27 frequencies between 0.1953 and 238.3 Hz) were inverted individually before being jointly inverted. The standard deviations assigned to the ZTEM data were 2.5 percent of the data plus a floor of 0.02. The standard deviations assigned to the MT data were 4 percent of the data plus a floor which was equal to 1 percent of the maximum absolute value of all the elements of the impedance tensor at the frequency of interest. The initial background conductivity was determined to be 0.125 S/m, by forward modelling different halfspace conductivity models and selecting the the halfspace conductivity which best fit the MT data. 112 El ev at io n (m ) 13- 4 56- 4 Range Front Fault S/m Easting (m) 1 0.1 0.01 0.001 Figure 4.9: ZTEM inversion going through drill holes 13-4 and 56-4. The inversion recovers a conductive structure on the west, due to the near surface sediments. The resistive structure east of the Range Front Fault is a result of the limestone unit shown in the geologic map in figure 4.7. The ZTEM inversion took around 12 hours to run on Intel Xeon X5660 dual hex core processors. The MT and joint ZTEM/MT inversions took longer to run because of the additional frequencies (these inversions were not run fully paral- lelized in order to limit the computer resources used). The final inversions fit the large-scale features of the data well; however, the joint ZTEM/MT inversion did not fit the data quite as well as the separate MT and ZTEM inversions. Cross- sections of the ZTEM, MT and joint ZTEM/MT inversion models between the two 113 El ev at io n (m ) 13- 4 56- 4 Range Front Fault S/m Easting (m) 1 0.1 0.01 0.001 Figure 4.10: MT inversion going through drill holes 13-4 and 56-4. The in- version recovers a conductive structure on the west, due to the near sur- face sediments. The resistive structure east of the Range Front Fault results from the limestone unit shown in the geologic map in figure 4.7. The MT data shows that the resistive structure extends to depths of greater than 1.5 km deep drill holes (13-4 and 56-4) can be seen in figures 4.9 - 4.11. On the western portion of the inverted slices, all three inversions recover a similar conductivity structure. The inversions recover a conductive feature likely related to tertiary sed- imentary rocks (Palacky, 1987) which were found at similar depths in the adjacent drill hole 13-4. In such a conductive environment, since the ZTEM system can typ- ically only measure data down to near 30 Hz, the method does not see very deeply, 114 Easting (m) El ev at io n (m ) 13- 4 56- 4 Range Front Fault S/m 1 0.1 0.01 Figure 4.11: Joint ZTEM and MT inversion going through drill holes 13-4 and 56-4. The inversion recovers a conductive structure on the west, due to the near surface sediments. The complex resistive structure east of the Range Front Fault results from the limestone unit shown in the geologic map in figure 4.7. In the complex geologic region near the limestone unit, the joint inversion result recovers features that are a combination of the features in the two individual inversions. and consequently the ZTEM inversion only recovers the shallower features. All of the deep information must come from the MT data. To the east of the Range Front Fault, the structure becomes much more complex. In the top few hundred meters, all three models show a resistive structure in agreement with the limestone unit from the geologic map in figure 4.7. While the general structure to the east 115 of the fault is resistive, the ZTEM and MT inversions in figures 4.9 and 4.10, each have some smaller conductive features which are not consistent between the two inversion models. This is understandable in such a complex 3D environment with significant topography. Relative to the survey station spacing and line spacing, the surface limestone structure and surrounding conductive rocks are small features. If it were desired to accurately resolve these small shallower discrepancies between the two inverted models, a more detailed ground electromagnetic survey with a smaller line and station spacing would be required. 4.3.3 Interpretations Determining large-scale structure is critical for understanding structural controls for fluid flow in geothermal exploration (Faulds et al., 2006). For a typical low temperature geothermal play (for example, see figure 4.12), there are many fea- tures that may help delineate potential prospect zones. For the Reese River Prop- erty there are three pieces of evidence that point towards a geothermal prospect zone on the western edge of the the Range Front Fault as seen in figure 4.7. Firstly, geologic mapping at the surface identified an accommodation zone where faults are interpreted to overlap and potentially create a zone of increased permeability at depth. Secondly, a shallow temperature gradient drilling program showed that the strongest thermal anomaly coincides with this zone. Finally, the large-scale struc- ture obtained by inverting the ZTEM and MT data shows a large resistive structure 116 coming up from over 1.5 km in depth under the geothermal prospect zone (as can be seen in figure 4.13 showing inversion model cells less conductive than 0.05 S/m). The dip on the western side of this deep structure is in agreement with the interpreted dip of the Range Front Fault. If hot fluid is up-welling in this region, then this resistive structure and the Range Front Fault are the main structures in the area and likely control the fluid flow in the region. Hot fluids could either travel in fractured rock along the fault boundaries, or within the rocks themselves if there is sufficient permeability. Because these rocks are less conductive than the surround- ing rocks, this structure could be interpreted to be low in clay content. With low clay content, the limestone and fault zones likely have increased permeability over the surrounding more conductive rocks (Cumming, 2009; Cumming and Mackie, 2010). Regardless of the fluid flow mechanism, the structure terminates below the interpreted accommodation zone, and the structure seems to follow the path of the fault fairly closely over the northern section of the dataset. Another well recovered feature in figure 4.13 is the limestone body to the east of the Range Front Fault, outlined in figure 4.7. Figure 4.14 shows a potential interpretation for the inverted conductivity model across the interpreted upflow zone from figure 4.7. In the top few hundred meters on the western portion of the cross section, are tertiary sedimentary rocks as logged by all 5 drill holes in the area. At deeper depths on the western portion of the cross 117 section, are more conductive rocks, which, from the drill logs, could be interpreted as Vinini Fmtn. Under the interpreted accommodation zone is the large resistive structure which coincides with the interpreted dip of the Range Front Fault. One interesting feature in the inverted model is the conductive feature buried a few hundred meters below the top of the Range Front mountain range. One possible explanation of this would be a clay cap associated with hot up-welling fluid along the eastern edge of the resistive structure. The hypothesis that the fluid is actually coming from the eastern side of the resistive structure, before cooling and trav- elling to the west into the sedimentary basin, may explain why many moderate temperature wells have been drilled to the west of the fault, yet the main source of the hot fluid has not been found. 118 Figure 4.12: Conceptual model reproduced from Cumming (2009) of a 150 \u00E2\u0097\u00A6C to 200 \u00E2\u0097\u00A6C geothermal reservoir with isotherms, alteration zones, and structures. 119 1924 -296 814 Easting (m) 485953 489053 487503 Range Front Fault 13-4 56-4 N S/m 1 0.1 0.01 Elevation (m ) Figure 4.13: Joint ZTEM/MT inversion showing the inverted model with a conductivity cutoff (only showing model cells less conductive than 0.05 S/m). The large resistive structure comes up from depth under the potential accommodation zone. The dip on the western side of the resistive structure seems to agree fairly well with the interpreted dip of the Range Front Fault. Possible Up-flow Paths Outflow? Range Front Fault Clay alteration? Limestone ? Tertiary Sediments ? Tertiary Sediments Siltstone & Shale (Vinini Fmtn) ? 62-4 Easting (m) El ev at io n (m ) 1 0.1 0.01 S/m Figure 4.14: Inverted ZTEM and MT model across the prospective geother- mal upflow zone shown in figure 4.7. Superimposed on the model is a potential geothermal interpretation. 120 Chapter 5 Conclusion 5.1 Research Summary Developing methods to image the subsurface of the earth is required for resource exploration, environmental monitoring and a variety of other applications. In this thesis I have examined how to image the electrical properties in the top few kilo- meters of the earth (the depths of investigation needed to find large deeply buried resource targets). To address this problem, I have examined natural source elec- tromagnetic methods and developed 3D inversion algorithms for these important data types. In particular, I have developed and demonstrated the first 3D inversion algorithm for the newly developed airborne ZTEM system. I have shown the ef- fectiveness of the data and our inversion methodology on both synthetic and field data examples. I am also the first to propose and implement the idea of jointly 121 using sparse ground MT data and airborne ZTEM data as a cost effective explo- ration method. Finally, I have developed a workflow type procedure that allows large natural source datasets to be inverted by working with multiple meshes and performing a domain decomposition of the full problem. The method has been demonstrated on both a synthetic and a field data example from the Pebble deposit in Alaska, where the large-scale ZTEM inversion seems consistent with both the DC resistivity and airborne time domain data results. 5.1.1 Inversion of ZTEM Data The ZTEM technique is a promising method to explore for large-scale targets at depth because of easy data acquisition and deep penetration of natural source fields. The data are transfer functions that relate the vertical magnetic field at the obser- vation point to the horizontal fields at a ground based reference station. In sum- mary, in the second chapter of this thesis, I have developed the first 3D inversion algorithm for ZTEM data. The inversion algorithm shows promising results on a synthetic test problem and illustrates the depth of penetration that can be ob- tained by using natural source fields. Because the 3D inversion of electromagnetic field data is a complex task, I present the essential inversion elements as a work- flow process. Using this approach, field data were inverted from the Mt. Milligan alkalic porphyry deposit in British Columbia. The final inverted model showed correspondences with some known geologic features from the area. The ZTEM 122 inversion model was also consistent in showing the same conductivity structures as the inversion results from a 3D airborne time domain dataset. 5.1.2 Inversion of Large-scale Datasets In the third chapter of this thesis, I have shown how data from a large-scale ZTEM survey can be inverted by using a strategy involving coarse and fine meshes, as well as a domain decomposition that splits the computational domain into smaller manageable subproblems which can be solved in parallel. I have presented an al- gorithm and practical workflow procedure to carry this out. Interactions between tiles are incorporated using primary fields in each tiled inversion that have been computed from a full domain conductivity model. To ensure continuity between overlapping tiles, I have used both data and model cell overlaps. The methodology shows promising results on a synthetic example from the Noranda mining camp. The procedure was also used to invert a 30 km x 30 km block of the 6504 line-km ZTEM survey collected over the Pebble Deposit in the Bristol Bay region of south- west Alaska. The results over this world class porphyry deposit are encouraging, and seem consistent with the other electromagnetic datasets. 5.1.3 Joint ZTEM and MT Inversion The ZTEM method is a promising method to explore for large scale targets at depth because of the rapid and cost effective airborne data acquisition, and the 123 deep penetration of natural source fields. While the depths of investigation are lim- ited somewhat by the lowest frequency signal that can be recorded in the air, the technique can still recover structures at depths that are difficult to image with other controlled source airborne and ground based systems. MT, the traditional natural source electromagnetic method, has the greatest depth of investigation of all elec- tromagnetic methods. Despite this, the costly nature of each MT station makes it generally impractical for large surveys that require high spacial resolution. A clear solution is to collect both ZTEM and MT data and exploit the strengths of each technique. High spacial resolution at moderate depths can then be achieved by fly- ing multiple ZTEM lines. Together with the ZTEM data, MT data can be used to estimate the background conductivity and resolve the deep structure. Based on the expected geologic model, synthetic modelling and survey design experiments can be performed to determine an appropriate combination of ZTEM and MT data. The result is a relatively easy to invert, cost effective solution that can map deep struc- tures over large areas. I have developed a Gauss-Newton style inversion algorithm to jointly invert ZTEM and MT data. The algorithm shows promising results on both a synthetic test problem and a field example from the Reese River geothermal project. 124 5.1.4 Future Research Opportunities Although the existing ZTEM method may be a good system to image deep tar- gets, there are possible developments with respect to the reference station that may improve ZTEM inversion results. One option would be to use multiple reference locations. A remote reference station (Clarke et al., 1983; Gamble et al., 1979b,a; Goubau et al., 1978) can be used to remove noise from measured magnetotelluric signals, and this could be applied to further remove unwanted noise from ZTEM signals. Since fixed wing airborne systems are much cheaper to operate than heli- copter systems, it would be desirable to develop a fixed wing version of the ZTEM system. Apart from the resulting differences in signal and noise characteristics between a helicopter and fixed-wing platform, there would be differences in flight heights and other survey parameters. These differences could be examined by com- paring inversion result differences between synthetic fixed-wing data and standard helicopter ZTEM data. While in this thesis I have shown that the ZTEM and MT techniques are effec- tive at imaging the earth using a number of synthetic and field examples, I have not performed a comprehensive comparison between these natural source methods and other deep probing electromagnetic methods. A logical future research direction is to perform a comprehensive synthetic modelling study to compare leading industry EM systems. During such a study, one could compare both the resolution and depth 125 of investigation of various EM systems. An extension of this study would be to in- corporate a survey design element, and determine the cheapest survey type together with line and station configurations, that would recover the desired resolution and depth information. While performing synthetic modelling to compare inversion results from var- ious airborne geophysical methods is a first and very important step in airborne systems and inversion codes validation, in order to make real world conclusions, system performances should be compared over a known test area. In recent years, this sort of rigorous validation approach has been gaining popularity, especially by government agencies such as Geoscience Australia and major mineral exploration companies that want to ensure they have a good understanding of each systems per- formance. For example, Geoscience Australia has established the Kauring airborne gravity test site (Dransfield et al., 2010) which provides a location for benchmark- ing the capabilities of airborne gravity (AG), airborne gravity gradiometer (AGG) and other airborne sensing systems. Many electromagnetic test ranges such as the Reid-Mahaffy (Reford and Fyon, 2000), Elura (Spies, 1980), Forrestania (10 air- borne and ground EM datasets) and Nepean (6 airborne and ground EM datasets operated by Southern Geoscience Consultants), now have multiple publicly avail- able datasets, which could and should be used to perform a comprehensive com- parison between the leading geophysical systems over known targets. 126 Another extremely important tool that has not been examined in this thesis is the use of geophysical constraints to improve inversion results. Adding geo- logic constraints as a priori information can significantly improve the inversion re- sult (for example Williams (2008) and Lelievre (2009)). Adding constraints could prove to be particularly valuable for the ZTEM method, which has limited sensitiv- ity to 1D layered structures and the background conductivity. It is often the case, as in the Pebble field data example in chapter 4, that regional surveys are flown over and around well established deposits, in the hopes of finding additional reserves in the surrounding area. In such a situation, there may be numerous drill holes with downhole conductivity information which can be used to constrain the inversion. Adding upper and lower bounds to the allowed conductivity values near the drill holes, even if only over a small region of the model, could affect the conductivity structure at other regions in the model and the final inversion result. 127 Bibliography Berdichevsky, M. N., and D. I. Vladimir, 2008, Models and methods of magne- totellurics: Springer. Blain, C., 2000, Fifty-year trends in minerals discovery \u00E2\u0080\u0093 commodity and ore-type targets: Exploration Mining Geology, 9, 1\u00E2\u0080\u009311. Bosworth, B., and S. M. Collins, 2007, Accounting for growth: Comparing China and India: Technical report, National Bureau of Economic Research. Cagniard, L., 1953, Basic theory of the magneto-telluric method of geophysical prospecting: Geophysics, 18, 605\u00E2\u0080\u0093635. Carcione, J. M., G. Seriani, and D. Gei, 2003, Acoustic and electromagnetic prop- erties of soils saturated with salt water and NAPL: Journal of Applied Geo- physics, 52, 177\u00E2\u0080\u0093191. Clark, M. K., 2000, Topographic ooze: Building the eastern margin of Tibet by lower crustal flow: Geology, 28, 703\u00E2\u0080\u0093706. Clarke, J., T. D. Gamble, W. M. Goubau, R. H. Koch, and R. F. Miracky, 1983, Remote-reference magnetotellurics: Equipment and procedures: Geophysical Prospecting, 31, 149\u00E2\u0080\u0093170. Commer, M., and G. A. Newman, 2009, Three-dimensional controlled-source elec- tromagnetic and magnetotelluric joint inversion: Geophysical Journal Interna- tional, 178, 1305\u00E2\u0080\u00931316. 128 Constable, S. C., A. S. Orange, G. M. Hoversten, and H. F. Morrison, 1998, Ma- rine magnetotellurics for petroleum exploration, part i: A sea-floor equipment system: Geophysics, 63, 816\u00E2\u0080\u0093825. Cox, L. H., G. A. Wilson, and M. S. Zhdanov, 2010, 3D inversion of airborne elec- tromagnetic data using a moving footprint: Exploration Geophysics, 41, 250\u00E2\u0080\u0093 259. Cumming, W., 2009, Geothermal resource conceptual models using surface ex- ploration data: Thirty-Fourth Workshop on Geothermal Reservoir Engineering Stanford University, 6. Cumming, W., and R. Mackie, 2010, Resistivity imaging of geothermal resources using 1D, 2D and 3D MT inversion and TDEM static shift correction illustrated by a Glass Mountain case history: Proceedings World Geothermal Congress. Devine, F., 2011b, Porphyry integration project: Bringing together geoscience and exploration datasets for British Columbias porphyry districts: Technical report, Geoscience BC. \u00E2\u0080\u0094\u00E2\u0080\u0094\u00E2\u0080\u0093, October 12, 2011a, Alkalic porphyry deposits in British Columbia deposit- to district-scale characteristics: Presented at the Geoscience BC: Exploration Undercover; a practical example using the QUEST study area. Dransfield, M., L. Roux, and T. Burrows, 2010, Airborne gravimetry and gravity gradiometry at Fugro Airborne Surveys, in Airborne Gravity 2010 abstracts from the ASEG-PESA Airborne Gravity 2010 Workshop: Geoscience Australia and the Geological Survey of New South Wales. Eisel, M., and V. Haak, 1999, Macro-anisotropy of the electrical conductivity of the crust: a magnetotelluric study of the German continental deep drilling site (ktb): Geophysical Journal International, 136, 109\u00E2\u0080\u0093122. Farquharson, C. G., and J. A. Craven, 2009, Three-dimensional inversion of mag- netotelluric data for mineral exploration: An example from the McArthur River 129 uranium deposit, Saskatchewan, Canada: Journal of Applied Geophysics, 68, 450\u00E2\u0080\u0093458. Farquharson, C. G., D. W. Oldenburg, E. Haber, and R. Shekhtman, 2002, An algo- rithm for the three-dimensional inversion of magnetotelluric data: 72nd Annual International Meeting, SEG, Expanded Abstracts, 649\u00E2\u0080\u0093652. Faulds, J. E., M. F. Coolbaugh, G. S. Vice, and M. L. Edwards, 2006, Charac- terizing structural controls of geothermal fields in northwestern Great Basin: A progress report: GRC Transactions, 30, 69\u00E2\u0080\u009376. Faulds, J. E., L. J. Garside, and G. L. Oppliger, 2003, Structural analysis of the Desert Peak-Brady geothermal fields, northwestern Nevada: Implications for understanding linkages between northeast-trending structures and geother- mal reservoirs in the Humboldt structural zone: Geothermal Resources Council Transactions, 27, 859\u00E2\u0080\u0093864. Friedman, S. P., 2005, Soil properties influencing apparent electrical conductivity: a review: Computers and Electronics in Agriculture, 46, 45\u00E2\u0080\u009370. Gamble, T. D., W. M. Goubau, and J. Clarke, 1979a, Error analysis for remote reference magnetotellurics: Geophysics, 44, 959\u00E2\u0080\u0093968. \u00E2\u0080\u0094\u00E2\u0080\u0094\u00E2\u0080\u0093, 1979b, Magnetotellurics with a remote magnetic reference: Geophysics, 44, 53\u00E2\u0080\u009368. Garcia, X., and A. G. Jones, 2002, Atmospheric sources for audio-magnetotelluric (AMT) sounding: Geophysics, 67, 448\u00E2\u0080\u0093458. Gibson, H., and A. G. Galley, 2007, in Mineral Deposits of Canada: A Synthe- sis of Major Deposit-Types, District Metallogeny, the Evolution of Geological Provinces: Mineral Deposits Division, 5, 533\u00E2\u0080\u0093552. Gokarn, S. G., C. K. Rao, G. Gupta, B. P. Singh, and M. Yamashita, 2001, Deep crustal structure in central India using magnetotelluric studies: Geophysical Journal International, 144, 685\u00E2\u0080\u0093694. 130 Goubau, W. M., T. D. Gamble, and J. Clarke, 1978, Magnetotelluric data analysis: Removal of bias: Geophysics, 43, 1157\u00E2\u0080\u00931166. Gribenko, A., A. M. Green, M. Cuma, and M. S. Zhdanov, 2010, Efficient 3D inversion of MT data using integral equations method and the receiver footprint approach: application to the large-scale inversion of the EarthScope MT data: SEG Expanded Abstracts Denver. Grimm, R. E., 2002, Low-frequency electromagnetic exploration for groundwater on Mars: Journal of Geophysical Research, 107. Haber, E., U. Ascher, and D. Oldenburg, 2000a, On optimization techniques for solving nonlinear inverse problems: Inverse Problems, 16, 1263\u00E2\u0080\u00931280. Haber, E., U. M. Ascher, D. A. Aruliah, and D. W. Oldenburg, 2000b, Fast simula- tion of 3D electromagnetic problems using potentials: Journal of Computational Physics, 163, 150\u00E2\u0080\u0093171. Haber, E., and S. Heldmann, 2007, An OcTree multigrid method for quasi-static Maxwell\u00E2\u0080\u0099s equations with highly discontinuous coefficients: Journal of Compu- tational Physics, 223, 783\u00E2\u0080\u0093796. Herrick, D. C., and W. D. Kennedy, 1994, Electrical efficiency - A pore geometric theory for interpreting the electrical properties of reservoir rock: Geophysics, 59, 918\u00E2\u0080\u0093927. Holtham, E., and D. W. Oldenburg, 2010, Three-dimensional inversion of ZTEM data: Geophysical Journal International, 182, 168\u00E2\u0080\u0093182. \u00E2\u0080\u0094\u00E2\u0080\u0094\u00E2\u0080\u0093, 2012, Large-scale inversion of ZTEM data: Geophysics, 77, 1\u00E2\u0080\u00939. Hoversten, G. M., H. F. Morrison, and S. C. Constable, 1998, Marine magnetotel- lurics for petroleum exploration, part ii: Numerical analysis of subsalt resolu- tion: Geophysics, 3, 826\u00E2\u0080\u0093840. Jefferson, C. W., D. J. Thomas, S. S. Gandhi, P. Ramaekers, G. Delaney, D. Brisbin, C. Cutts, D. Quirt, P. Portella, and R. A. Olson, 2007, Unconformity-associated 131 uranium deposits of the Athabasca Basin, Saskatchewan and Alberta: Geologi- cal Association of Canada, Mineral Deposits Division, Special Publication, St. Johns, 5, 273305. Jones, A., and I. Ferguson, 2001, The electric moho: Nature, 409, 331\u00E2\u0080\u0093333. Karous, M., and S. E. Hjelt, 1983, Linear filtering of vlp dip-angle measurements: Geophysical Prospecting, 31, 782\u00E2\u0080\u0093794. Key, K., 2011, Marine electromagnetic studies of seafloor resources and tectonics: Surveys in Geophysics, 33, 135\u00E2\u0080\u0093167. Key, K. W., S. C. Constable, and C. J. Weiss, 2006, Mapping 3D salt using the 2D marine magnetotelluric method: Case study from Gemini Prospect, Gulf of Mexico: Geophysics, 71, B17\u00E2\u0080\u0093B27. Kirkendall, B., and J. Roberts, 2004, Electromagnetic imaging of CO2 seques- tration at an enhanced oil recovery site: Techn. Rep. No. UCRL-TR-204708, Lawrence Livermore Nat. Lab. Labson, V. F., A. Becker, H. F. Morrison, and U. Conti, 1985, Geophysical explo- ration with audiofrequency natural magnetic fields: Geophysics, 50, 656\u00E2\u0080\u0093664. Lang, J., K. Roberts, M. Rebagliati, M. Gregory, C. Harradan, B. McNulty, R. Smith, and G. Wober, 2009, Geological relationships & controls on mineraliza- tion in the Pebble Cu-Au-Mo porphyry deposit, Alaska: Technical report, The Pebble Partnership. Leggatt, P. B., P. S. Klinkert, and T. B. Hage, 2000, The Spectrem airborne elec- tromagnetic system\u00E2\u0080\u0094Further developments: Geophysics, 65, 1976\u00E2\u0080\u00931982. Lelievre, P. G., 2009, Integrating geologic and geophysical data through advanced constrained inversions: PhD thesis, UBC. Lesmes, D. P., and S. P. Friedman, 2005, Relationships between the electrical and hydrogeological properties of rocks and soils: Hydrogeophysics, 50, 87\u00E2\u0080\u0093128. Lo, B., and M. Zang, 2008, Numerical modeling of Z-TEM (airborne AFMAG) 132 responses to guide exploration strategies: SEG Technical Program Expanded Abstracts, 27, 1098\u00E2\u0080\u00931102. Mackie, R. L., and T. R. Madden, 1993, Three-dimensional magnetotelluric in- version using conjugate gradients: Geophysical Journal International, 115, 215\u00E2\u0080\u0093 229. Mackie, R. M., M. D. Watts, and W. Rodi, 2007, Joint 3D inversion of marine CSEM and MT data: SEG Technical Program Expanded Abstracts, 26, 574\u00E2\u0080\u0093 578. Meju, M. A., 2002, Geoelectromagnetic exploration for natural resources: models, case studies and challenges: Surveys in Geophysics, 23, 133\u00E2\u0080\u0093205. Mills, K., J. Huang, G. Bosworth, S. Cowie, B. Borntraeger, H. Welhener, J. Collins, T. Bekhuys, and D. Labrenz, 2009, Technical Report - Feasibility Up- date Mt. Milligan Property - Northern BC: Technical report, Report to: Terrane Metals Corp. prepared by Wardrop. Mitchinson, D. E., and R. J. Enkin, 2011, Continued investigations of physical property - geology relationships in porphyry-deposit settings in the QUEST and QUEST-West Project Areas, central British Columbia (NTS 093E, K, L, M, N): Technical report, Geoscience BC. Mitsuhata, Y., K. Matsuo, and M. Minegish, 1999, Magnetotelluric survey for ex- ploration of a volcanic-rock reservoir in the Yurihara oil and gas field, Japan: Geophysical Prospecting, 47, 195\u00E2\u0080\u0093218. Nel, W. P., and C. J. Cooper, 2008, A critical review of IEA\u00E2\u0080\u0099s oil demand forecast for China: Energy Policy, 36, 1096\u00E2\u0080\u00931106. Nelson, J., K. Bellefontaine, K. Green, and M. MacLean, 1990, Regional geolog- ical mapping near the Mount Milligan copper-gold deposit: British Columbia Geological Survey Geological Fieldwork, 89\u00E2\u0080\u0093110. Newman, G. A., and D. L. Alumbaugh, 2000, Three-dimensional magnetotelluric 133 inversion using non-linear conjugate gradients: Geophysical Journal Interna- tional, 140, 410\u00E2\u0080\u0093424. Oldenburg, D., Y. Li, and R. Ellis, 1997, Inversion of geophysical data over a copper gold porphyry deposit: a case history for Mt. Milligan: Geophysics, 62, 1419\u00E2\u0080\u00931431. Oldenburg, D. W., E. Haber, and R. Shekhtman, 2008, Forward modelling and inversion of multi-source TEM data: SEG Annual Meeting Las Vegas. Palacky, G. J., 1987, Resistivity characteristics of geologic targets: in Electromag- netic Methods in Applied Geophysics, Vol 1, Theory, 53\u00E2\u0080\u0093129. Palacky, G. J., and G. F. West, 1987, Airborne electromagnetic methods: in Elec- tromagnetic Methods in Applied Geophysics, Vol 1, Applications, 811\u00E2\u0080\u0093877. Pare, P., and J. M. Legault, 2010, Ground IP-Resistivity, and airborne Spectrem and helicopter ZTEM survey results over Pebble copper-moly-gold porphyry deposit, Alaska: SEG Annual Meeting Meeting expanded abstracts, 1734\u00E2\u0080\u00931738. Pasion, L. R., S. D. Billings, D. W. Oldenburg, and S. E. Walker, 2007, Application of a library based method to time domain electromagnetic data for the identifi- cation of unexploded ordnance: Journal of Applied Geophysics, 61, 279\u00E2\u0080\u0093291. Pellerin, L., 2002, Applications of electrical and electromagnetic methods for envi- ronmental and geotechnical investigations: Surveys in Geophysics, 23, 101\u00E2\u0080\u0093132. Pellerin, L., J. M. Johnston, and G. W. Hohmann, 1996, A numerical evaluation of electromagnetic methods in geothermal exploration: Geophysics, 61, 121\u00E2\u0080\u0093137. Phillips, N., T. N. Nguyen, and V. Thompson, 2009, QUEST Project: 3D inversion modeling, integration, and visualization of airborne gravity, magnetic, and elec- tromagnetic data, BC, Canada: Technical report, Mira Geoscience Advanced Geophysical Interpretation Centre. Rebagliati, C. M., J. Lang, E. Titley, D. Rennie, L. Melis, D. Barratt, and S. Hodg- son, 2009, Technical report on the 2008 program and update on mineral re- 134 sources and metallurgy, Pebble copper-goldmolybdenum project, Iliamna Lake area, southwestern Alaska, U.S.A.: Technical report, Northern Dynasty Miner- als Ltd. Reford, S., and A. Fyon, 2000, Airborne magnetic and electromagnetic surveys Reid-Mahaffy airborne geophysical test site survey miscellaneous release data (MRD) 55: Technical report, Ontario Geological Survey. Saad, Y., 2003, Iterative methods for sparse linear systems: Society for Industrial and Applied Mathematics. Sasaki, Y., 2004, Three-dimensional inversion of static-shifted magnetotelluric data: Earth, Planets Space, 56, 239\u00E2\u0080\u0093248. Schwartz, H. A., 1870, Gesammelte mathematische abhandlungen: Vierteljahrss- chrift der Naturforscheden Gesellschaft, 15, 272\u00E2\u0080\u0093286. Schwarzbach, C., R.-U. Bo\u00CC\u0088rner, and K. Spitzer, 2011, Three-dimensional adap- tive higher order finite element simulation for geo-electromagnetics\u00E2\u0080\u0094a marine CSEM example: Geophysical Journal International, 187, 63\u00E2\u0080\u009374. Shalf, J., 2007, The new landscape of parallel computer architecture: Journal of Physics: Conference Series, 78. Sheard, S. N., T. J. Ritchie, K. R. Christopherson, and E. Brand, 2004, Mining, environmental, petroleum, and engineering industry applications of electromag- netic techniques in environmental and geotechnical investigations: Surveys in Geophysics, 23, 101\u00E2\u0080\u0093132. Simpson, F., and K. Bahr, 2005, Practical magnetotellurics: Cambridge University Press. Siripunvaraporn, W., G. Egbert, Y. Lenbury, and M. Uyeshima, 2005, Three- dimensional magnetotelluric inversion: data-space method: Physics of the Earth and Planetary Interiors, 150, 3 \u00E2\u0080\u0093 14. (Electromagnetic Induction in the Earth). Smirnov, M., T. Korja, L. Dynesius, L. B. Pedersen, and E. Laukkanen, 2008, 135 Broadband magnetotelluric instruments for near-surface and lithospheric studies of electrical conductivity: A fennoscandian pool of magnetotelluric instruments: Geophysica, 44, 31\u00E2\u0080\u009344. Spies, B., 1980, Results of experimental and test TEM surveys, Elura deposit, Cobar, NSW: Bulletin of the Australian Society of Exploration Geophysicists, 11, 147\u00E2\u0080\u0093152. Tikhonov, A. N., 1950, The determination of the electrical properties of deep layers of the earth\u00E2\u0080\u0099s crust: Dokl. Acad. Nauk. SSR, 73, 295\u00E2\u0080\u0093297. Wannamaker, P. E., T. G. Caldwell, G. R. Jiracek, V. Maris, G. J. Hill, Y. Ogawa, H. M. Bibby, S. L. Bennie, and W. Heise, 2009, Fluid and deformation regime of an advancing subduction system at Marlborough, New Zealand: Nature, 460, 733\u00E2\u0080\u0093736. Ward, S. H., 1959, AFMAG \u00E2\u0080\u0094airborne and ground: Geophysics, 24, 761\u00E2\u0080\u0093787. Ward, S. H., J. O\u00E2\u0080\u0099Donnell, R. Rivera, G. H. Ware, and D. C. Fraser, 1966, AFMAG - applications and limitations: Geophysics, 31, 576\u00E2\u0080\u0093605. Wei, W., M. Unsworth, A. Jones, J. Booker, H. Tan, D. Nelson, L. Chen, S. Li, K. Solon, P. Bedrosian, S. Jin, M. Deng, J. Ledo, D. Kay, and B. Roberts, 2001, Detection of widespread fluids in the Tibetan crust by magnetotelluric studies: Science, 292, 716\u00E2\u0080\u0093719. Williams, N. C., 2008, Geologically-constrained UBC-GIF gravity and magnetic inversions with examples from the Agnew-Wiluna Greenstone Belt, Western Australia: PhD thesis, UBC. Witter, J. B., J. A. Ronne, and G. R. Thompson, 2009, Exploration at the Reese River geothermal prospect in Nevada: Geothermal Resources Council Transac- tions, 33, 549\u00E2\u0080\u0093556. Yang, D., and D. W. Oldenburg, 2012, Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit: Geo- 136 physics, 77, B23\u00E2\u0080\u0093B34. Yellishettya, M., P. G. Ranjitha, and A. Tharumarajah, 2010, Iron ore and steel production trends and material flows in the world: Is this really sustainable?: Resources, Conservation and Recycling, 54, 1084\u00E2\u0080\u00931094. Zhdanov, M. S., S. Fang, and G. Hursan, 2000, Electromagnetic inversion using quasi-linear approximation: Geophysics, 65, 1501\u00E2\u0080\u00931513. Zhdanov, M. S., L. Wan, A. Gribenko, M. Cuma, K. Key, and S. Constable, 2011, Large-scale 3D inversion of marine magnetotelluric data: Case study from the Gemini prospect, Gulf of Mexico: Geophysics, 76, F77\u00E2\u0080\u0093F87. 137"@en . "Thesis/Dissertation"@en . "2012-11"@en . "10.14288/1.0053566"@en . "eng"@en . "Geophysics"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "Attribution-NonCommercial-NoDerivatives 4.0 International"@en . "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en . "Graduate"@en . "3D inversion of natural source electromagnetic data"@en . "Text"@en . "http://hdl.handle.net/2429/42582"@en .