UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Geophysical inversion in an integrated exploration program : examples from the San Nicolas deposit Phillips, Nigel David 2002-08-14

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
831-ubc_2002-0206.pdf [ 12.54MB ]
Metadata
JSON: 831-1.0052682.json
JSON-LD: 831-1.0052682-ld.json
RDF/XML (Pretty): 831-1.0052682-rdf.xml
RDF/JSON: 831-1.0052682-rdf.json
Turtle: 831-1.0052682-turtle.txt
N-Triples: 831-1.0052682-rdf-ntriples.txt
Original Record: 831-1.0052682-source.json
Full Text
831-1.0052682-fulltext.txt
Citation
831-1.0052682.ris

Full Text

GEOPHYSICAL INVERSION IN AN INTEGRATED EXPLORATION PROGRAM: EXAMPLES FROM THE SAN NICOLAS DEPOSIT By Nigel David Phillips B. Sc. Colorado School of Mines, 1996 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 2002 © Nigel David Phillips, 2002 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth and Ocean Sciences The University of British Columbia 129-2219 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract The recent ability to produce three-dimensional physical property models of the sub surface from surface geophysical data, coupled with an increasing need to explore for minerals in concealed terranes, results in geophysical inversions providing more signifi cant information to the exploration team. This thesis examines the role that geophysical inversion can play in an integrated mineral exploration program, and the impact it can have on the results. As an example, geophysical data from the San Nicolas copper-zinc massive sulphide deposit in Mexico are inverted. Such deposits are distinguished by high density, magnetic susceptibility, conductivity, and chargeability values. Within the framework of an integrated exploration program (which is communicated through the use of a flowchart) many different, and hard-to-interpret, geophysical data sets are given geologic context. This is achieved by interpreting physical property models that have been generated by inversion modeling. The aim of generating these models, and interpreting geological information from them, is to: 1) assist mineral exploration in the area around the deposit, and roughly define areas of mineralization that can be used as a starting point for more detailed modeling ; 2) define the size and location of the San Nicolas deposit; 3) further improve the models and delineate ore more accurately by combining the detailed modeling with core physical property measurements; and 4) determine how best to use the information that has been acquired through modeling to find additional sulphide ore-bodies, including those that may be deeper than the existing deposit. Density and magnetic susceptibility distribution models, inverted from regional grav ity and magnetic data respectively, define large-scale structures that reflect the tectonic ii setting of the region. Several distinct anomalies that exhibit high density and magnetic susceptibility values are identified. Since massive sulphides are often dense and mag netic, a correlation method is employed that determines volumes that have high density and magnetic susceptibility. The correlation procedure isolates five anomalies. Two of these are easily dismissed as having poor exploration potential, and two of the remaining anomalies, one of which is the San Nicolas deposit, are associated with mineralization. At a more detailed scale, the deposit is well defined by gravity, magnetic, CSAMT, and IP methods individually. However, a drill hole that is targeted on the intersection of these favorable physical property distributions would have intersected the heart of the deposit. This demonstrates the advantages of using these methods in concert. Other sources of information, such as core physical property measurements and geologic constraints, are also used to improve modeling results. The inclusion of data from a single drill-hole is shown to significantly enhance detailed physical property distributions, and produces models that correlate better which known mineralization. •. Finally, forward modeling of the physical property models is used to demonstrate how deep ore can be detected by implementing different survey designs that increase signal-to-noise of the data. This newly acquired information can then be used in the next step of the exploration program as the search for mineralization continues in the surrounding area. In addition to the main focus of the thesis, it was found that polarizable material is needed in order to fit time-domain, airborne electromagnetic data collected over the deposit. This demonstrates the potential for detecting chargeable bodies from the air, and has significant implications for future mineral exploration. iii Table of Contents Abstract ii List of Tables xList of Figures xii Acknowledgments xxviiDedication xxix 1 Introduction 1 1.1 Exploration methodology 3 1.2 San Nicolas 5 1.2.1 Geology 6 1.3 Physical Properties 10 1.4 Geophysical Exploration Model 12 1.5 Data 11.6 Modeling 3 1.6.1 Forward modeling 11.6.2 Inverse modeling1.7 Outline of the thesis 4 1.7.1 Inversion methodology 11.7.2 Regional-scale modeling 4 1.7.3 Deposit-scale modeling 1iv 1.7.4 Detailed modeling using drill-hole information 16 1.7.5 Detection of deep ore: survey design 11.8 Summary 17 2 Geophysical Inversion 18 2.1 Introduction2.2 Inversion Basics2.2.1 Data 18 2.2.2 Discretization and forward modeling 19 2.2.3 Non-uniqueness and the Model Objective Function 20 2.2.4 Signal-to-noise and data misfit 22 2.2.5 Optimization 23 2.2.6 Appraisal 4 2.3 A synthetic example 5 2.3.1 The data 26 2.3.2 Presentation and interpretation of inversion models 27 2.3.3 Recovered models 22.4 Conclusions 32 3 Regional Scale Modeling 34 3.1 Introduction3.2 Regional geology 6 3.3 Introduction to the inversion results 40 3.4 Density-contrast models 43.4.1 Regional gravity inversion results 42 3.4.2 Summary of regional gravity results 9 3.5 Magnetic susceptibility models 4v 3.5.1 Data 50 3.5.2 Composite susceptibility model 52 3.5.3 Regional susceptibility model analysis 53.5.4 Geologic interpretation 56 3.6 Regional conductivity models 7 3.6.1 Frequency-domain data 53.6.2 Time-domain data 60 3.6.3 Forward modeling the Questem data 67 3.6.4 Summary of regional conductivity methods 73 3.7 Chargeability models 74 3.8 Physical property correlation 8 3.9 Conclusions 81 4 Deposit Scale Modeling 84 4.1 Introduction4.1.1 Local geology 6 4.2 Density 89 4.2.1 Regional Removal 84.2.2 Inversion of Residual Gravity Data 89 4.3 Magnetic Susceptibility • • • 92 4.4 Conductivity Data 96 4.4.1 Testing the 3D resistivity model 97 4.5 Chargeability 8 4.5.1 Chargeability data 101 4.5.2 The resistivity model4.5.3 Removal of electrode effects 103 vi 4.5.4 The chargeability model 103 4.6 Conclusions 105 5 Detailed Modeling: Incorporating a priori information 108 5.1 Introduction 105.1.1 A priori information at San Nicolas 109 5.2 Inversion of gravity observations Ill 5.2.1 Setting bounds on the density-contrast model 113 5.2.2 A priori density values from drill-hole SAL-31 117 5.2.3 A priori density values from all drill-holes 119 5.3 Magnetic Susceptibility 123 5.3.1 Core Measurements5.3.2 Bounding the magnetic susceptibility model 125 5.3.3 Magnetic inversion using a priori information from SAL-69 .... 125 5.3.4 Magnetic inversion using a priori information from all drill-holes . 126 5.4 Discussion on the discontinuous nature of the models 129 5.5 Conclusion 126 Forward Modeling and Survey Design: Detecting the keel 132 6.1 Introduction 136.2 Gravity Observations 135 6.2.1 Simulation of gravity field observations 136 6.2.2 Synthetic, downhole, gravity observations 138 6.2.3 Summary of gravity forward modeling 142 6.3 Magnetic Observations 146.3.1 Simulating the ground magnetic field data 143 6.3.2 Synthetic downhole magnetic observations 14vii 6.3.3 Summary of magnetic forward modeling 147 6.4 CSAMT/CSEM observations 146.4.1 Simulated CSAMT observations .... . . 148 6.4.2 CSAMT/CSEM survey design 152 6.4.3 Summary of CSAMT/CSEM forward modeling 158 6.5 DC resistivity 159 6.5.1 Simulation of Realsection resistivity observations 159 6.5.2 Forward modeling of dipole-dipole resistivity observations 160 6.5.3 Summary of DC resistivity forward modeling ;. ... . 162 6.6 Chargeability 164 6.6.1 Simulating Realsection IP field observations 166 6.6.2 Synthetic dipole-dipole chargeability observations 167 6.6.3 Summary of chargeability forward modeling 167 6.7 Conclusions 169 7 Conclusions and Suggestions for Further Work 171 7.1 Summary of results 173 7.2 Further work 5 References 178 Appendices 184 A Coordinate System 18B Physical Property Analysis 185 B.l Introduction 18B.2 Densityviii B.3 Magnetic Susceptibility 185 B.3.1 Magnetic Susceptibility core measurements 186 B. 3.2 Composite models for bounding inversions 188 B.4 Resistivity 189 B.5 ChargeabilityB. 6 Physical property rock-models 191 C Geophysical methods 193 Cl The gravity methodC. 2 The magnetic method 195 C.3 Inductive electromagnetic methods 197 C. 3.1 Frequency-domain electromagnetic airborne surveys. . 197 C.3.2 Time-domain electromagnetic airborne surveys 199 C.4 CSAMT surveys 200 C.5 DC resistivity methods 203 C.6 IP methods 204 D Data 206 E Inversion 220 E.l IntroductionE.2 GRAV3D 22E.3 MAG3D 1 E.4 CSAMT1D 22E.5 EM1DFM 2 E.6 TEM ID 22E.6.1 Complex forward modeling 222 ix E.7 DCIP3D 222 E.8 EHSD 3 E. 9 EH3DTD 22F Gradient IP data modeled using MAG3D 225 F. l Four synthetic examples 226 F.2 Modeling the electric field 9 F.3 Summary 233 G IP electrode weighting 234 x List of Tables 1.1 Estimated physical properties for the five major rock types encountered at San Nicolas. The massive sulphide has contrasting values in all physical properties apart from resistivity. A low resistivity value of 20 Ohm-m for sulphides, similar to the Tertiary breccia, makes detection of the deposit more difficult when using some geophysical methods 11 3.1 Questem time-decay (^) data for one station located over the deposit. The center time of each time channel is given (time after turn-off of primary field) along with the z-component data in micro-volts and in ppm, as normalized by the peak z-component of the primary field 65 B.l Range of densities for rock-types encountered at San Nicolas 186 B.2 Range of resistivities for rock-types encountered at San Nicolas. ..... 190 B.3 Range of chargeabilities for rock-types encountered at San Nicolas 191 D. l Geophysical survey specifications 219 E. l Inversion parameters 224 xi List of Figures 1.1 Geophysical Exploration Methodology Flowchart: Defining the geophysi cal exploration model 2 1.2 Map of Mexico (upper panel) and Zacatecas State (lower panel) showing the location of San Nicolas 7 1.3 Data coverage over the San Nicolas deposit and the surrounding area. Top: the regional-scale area with data coverage and outcrop geology (leg end shown). Center: Intermediate-scale data coverage shown on top of outcrop geology. Bottom: deposit-scale data coverage. The size of the rock-hammer symbol approximately corresponds to the extents of the sur face projection of the San Nicolas deposit 15 2.1 The earth is discretized and represented by an orthogonal 3D mesh. ... 20 2.2 Perspective view of synthetic density-contrast model. The model has been volume rendered so that cells with density contrasts less than lg/cm3 have been removed 25 2.3 Surficial synthetic gravity observations with 3% + 0.0005 mGal random noise added 6 2.4 Perspective view of recovered density-contrast model. The model is viewed from the southeast (local coordinate system), and is volume-rendered with a cut off at 0.3 g/cm3 28 2.5 North-facing cross-section through recovered density-contrast model with location of true model shown 28 xii 2.6 North-facing cross-sections through true model and seven recovered models for seven different inversions. All of the inversion models are shown on the same color-scale with dark blue representing 0 g/cm3 and below, and magenta representing 0.5 g/cm3 and above. The original model, an outline of which is shown on all the recovered models, is shown at a different scale. 30 3.1 Geophysical Exploration Methodology Flowchart: Regional Scale Modeling. 35 3.2 Plan-view of regional geologic outcrop at San Nicolas. Modified from John son et al., 2000. Most of the region is overlain by Quarternary alluvium. Inset shows a key to the stratigraphic units 37 3.3 Plan view of the regional density contrast model with cells below a value of 0.1 g/cm3 invisible. Outcrop geology is overlaid, and the boundary of the gravity survey is represented by a purple line. All coordinates are in the local grid system 43 3.4 West-facing cross section of the regional density contrast model with cells below a value of 0.1 g/cc invisible. This cross-section is located at -1600m east and shows an interpreted horst 44 3.5 Surface projection of the regional density-contrast model with interpreted faults overlaid. Northwest-southeast trending faults are shown in red, while north-northeast trending faults are shown in black. The mapped fault is shown in yellow. Material with a density-contrast of less than 0.1 g/cm3 has been removed 45 3.6 Horizontal slice through the regional density-contrast model at an eleva tion of 1700m (400m below the surface at the location of San Nicolas). An outline of geology outcrops is overlaid, and the extent of data coverage is shown by a purple line 47 xiii 3.7 North facing cross-section through a regional density contrast model at a northing of -400m (coincident with the deposit) 48 3.8 Comparison of individual magnetic susceptibility inversion models with composite model of mean values. The models are shown in plan-view looking at a depth of 300m (within the depth of the deposit) 51 3.9 Plan-view of horizontal slices through the composite magnetic susceptibil ity model. The upper panel shows a slice at 300m depth, while the lower panel shows a slice through the model at 1200m depth 53 3.10 North- and west-facing cross-sections (upper and lower sections respec tively) through the magnetic susceptibility model generated by combin ing inversion results from Dighem, Fugro, and BHP airborne magnetic data. The north-facing cross-section is positioned at -400m south, while the west-facing cross-section is positioned at -1600m east 55 3.11 Horizontal slice 60m below the surface of the 3D resistivity model gener ated by interpolating between lines of ID models inverted from Dighem frequency-domain airborne electromagnetic data 59 3.12 Plan view of a slice through the 3-D resistivity model generated by in terpolating between lines of ID models from the inversion of Questem time-domain airborne electromagnetic data. The slice (at a depth of 200m below the surface) shows, in red, the thicker areas of overburden that con tinue to this depth. The underlying, more resistive, volcanic rocks are represented by cooler colors 62 xiv 3.13 Vertical component, time-channel data, flight height, total field magnet ics, and resistivity inversion model for line 20110 of the Questem data. The data are presented as parts-per-million of the peak z-component of the waveform, with channels 1 and 15 corresponding to the earliest and latest times respectively. The inversion model is produced by "stitching" together ID models. The location of the deposit corresponds to a coherent pattern of negative data values in time-channels 12 through 15 (1), and a magnetic high. Other regions of negative data that could be above the noise level are also observed (2, 3 and 4) 64 3.14 Logarithmic plot of Questem vertical-component time-decay data given in Table 3.1. Positive data are shown by squares, and negative data (the last four) are shown with triangles. The data are from over the deposit and show a typical decay pattern for stations that contains negative data. . . 66 3.15 Observed and predicted time-decay data, and corresponding ID conduc tivity models for Questem station located over the deposit. Frequency-independent conductivity layers do not reproduce the observed negative data . . . . . 69 3.16 Observed and predicted data, and ID conductivity models for Questem station located over the deposit. A layer with complex (frequency-dependent) conductivity (based on the Cole-Cole dispersion model), positioned at the approximate depth of the deposit, can produce data that has the same character as the observations. The conductivity of the overburden was also decreased in an attempt to more closely fit the early-time data. ... 71 xv 3.17 Plan-view of a slice at 300m below the surface through the chargeability model generated by inverting gradient-array IP data as magnetic ..data using MAG3D. The deposit is the largest chargeable body detected by the survey. The extents of the survey are shown by the red outline and outcrop geology is shown in black 75 3.18 North- and west-facing cross-sections (upper and lower sections respec tively) through the chargeability model. The chargeable anomaly is well defined yet extends below the elevation of the deposit, as defined by drilling. 76 3.19 Plan-views of density contrast, magnetic susceptibility, and cut-off corre lation models. The upper images show the projection-to-surface of cells, within the first 1000m from the surface, that exhibit high density contrast and magnetic susceptibility values. The lower image shows regions that are high in both density and susceptibility. The image is a projection to surface of the volume intersection of the upper models. A blue border shows the extents of the correlation analysis 79 3.20 Normalized density-susceptibility regional correlation values. The values are the product of normalized regional density-contrast and magnetic sus ceptibility models 82 4.1 Geophysical Exploration Methodology Flowchart: Deposit scale modeling. 85 4.2 Geologic cross-sections through the San Nicolas deposit, modified from sections provided by Teck Corp. The main geologic units are labeled. . . 87 4.3 Simplified stratigraphic section in the vicinity of the San Nicolas deposit, modified from [Danielson, 2000] 88 xvi 4.4 Perspective view, facing north-west, of regional density-contrast model used for calculation of the regional signal. The volume that has been removed contains the deposit and two smaller high density anomalies, and is the volume on which the more detailed inversion is focused 90 4.5 South-east perspective view of the detailed, density-contrast model shown as a volume-rendered isosurface plot using a cut-off of 0.18 g/cc. Cells with density contrasts less than this value are invisible 91 4.6 Cross-sections through deposit-scale density-contrast model at the San Nicolas deposit 93 4.7 Perspective view of the detailed, magnetic susceptibility model 94 4.8 Cross-sections through deposit-scale magnetic susceptibility model at the San Nicolas deposit 95 4.9 Perspective view of resistivity model created by interpolating between lines of "stitched" ID models. The inversion has successfully discriminated between the overburden and the deposit 96 4.10 Cross section of resistivity model with geology. The resistivity model was produced by "stitching" together 60 ID inversion models 97 4.11 Comparison of CSAMT field data with forward modeled data from the stitched 1-D inversion resistivity model(left plots) and a physical property resistivity model (right plots). The field observations are displayed in red while the modeled data is shown in blue. At 8Hz (upper plots), the physical property model fits the field data well while the 1-D inversion model fits the field data better at a frequency of 256Hz (lower plots). . . 99 xvii 4.12 North-facing, cross-sectional diagram of the geology at San Nicolas show ing the skin-depths associated with propagating plane-waves at frequen cies of 256Hz and 8Hz. The skin-depths were calculated for a 25 Ohm-m whole space. The higher frequency signal samples geology that is much more "one-dimensional" than that sampled by the low frequency signal. . 100 4.13 Plan-views of the recovered chargeability models at the surface. The high chargeability values that correspond with electrode positions have been reduced through the use of a weighting function 104 4.14 Perspective view of the volume-rendered chargeability model with a cut-off value of 40 msec 105 4.15 Cross-sections through deposit-scale chargeability model at the San Nicolas deposit 106 5.1 Geophysical Exploration Methodology Flowchart: Using drill-hole infor mation 110 5.2 Plan-view of drill-hole locations. Density measurements on core from.drill-hole SAL-31 are used to set bounds on the model values in Section 5.2.2 (shown by a red square), while measurements made on core from all of the drill-holes shown (apart from SAL-90 and SAL-91) are used to pro vide a priori density contrast information in Section 5.2.3. Susceptibility measurements made on drill-cores from SAL-69 are used to set bounds in Section 5.3.3 (shown by a green triangle), while bounds are assigned in Section 5.3.4 based on magnetic susceptibility values made on cores from thirteen drill-holes that are denoted by a blue diamond 112 xviii 5.3 Drill-hole SAL-31. Displayed as a function of depth from left to right: 1) Geology, 2) density (g/cm3), 3) percent copper, and 4) percent zinc. The average density of the sulphide is 4.2 g/cm3 and the density of the rocks below the deposit is approximately 2.8 g/cm3 114 5.4 North- and west-facing cross-sections through recovered deposit-scale density-contrast models. The upper cross-sections show the results from the orig inal deposit-scale gravity inversion, described in Section 4.2.2. The lower cross-sections show the results of inverting the gravity data with the model bounded using the density measurements made on core from drill-hole SAL-31. Geologic contacts have been overlaid 118 5.5 North-facing cross-section of interpreted geology with drill-hole locations. Density measurements made on core from the drill-holes were used to bound model cells during the inversion of gravity data. Only the drill holes along line -400m S are displayed, although many more were used to bound the model 120 5.6 Cross-sections through recovered deposit-scale density-contrast models. The upper cross-sections show the results from the original deposit-scale gravity inversion, described in Section 4.2.2. The lower cross-sections show the results of inverting the gravity data with the model bounded by den sity measurements made on core from all drill-holes. Interpreted geologic contacts are shown on each cross-section. Note: the lower cross-sections are shown with a range of density-contrast that is different from those previously presented 122 xix 5.7 Magnetic susceptibility core measurements and geology for drill-holes SAL-25, SAL-34, SAL-47, and SAL-70. Note: the upper-most geologic unit is the Tertiary breccia, although it may appear a similar colour to the sulphide ore, which is shown in red. Reference the geologic section and recognize it as an interpretation 124 5.8 North- and west-facing cross-sections through recovered deposit-scale mag netic susceptibility models. The upper cross-sections show the results from the original deposit-scale magnetic susceptibility inversion, described in Section 4.3. The lower cross-sections show the results of inverting the magnetic data with the model bounded by magnetic susceptibility mea surements made on core from drill-hole SAL-69. An outline of geology interpreted from drill-holes is overlaid on each section 127 5.9 North and west-facing cross-sections through recovered deposit-scale mag netic susceptibility models. The upper cross-sections show the results from the original deposit-scale magnetic inversion, described in Section 4.3. The lower cross-sections show the results of inverting the magnetic data with the model bounded by magnetic susceptibility measurements made on core from twelve drill-holes. Geologic contacts as interpreted from drill-holes are overlaid on each section 128 6.1 Geophysical Exploration Methodology: Survey design 134 6.2 Size and position of the keel relative to the deposit 135 6.3 The two density rock-models that were used to generate synthetic gravity measurements 137 xx 6.4 Synthetic gravity data produced from a density rock-model of San Nicolas. Data are produced from models with and without the keel. The difference is also plotted 139 6.5 West-facing cross-section of density rock-model. The main ore-zone is represented by the higher density values in the top half of the model while the keel is represented by the anomaly in the lower half of the model. The locations of six drill-holes, within which synthetic gravity measurements are made, are also shown. Drill-hole A is 25m north of the keel and each subsequent drill-hole is positioned 25m further north than the last one. . 140 6.6 Down-hole plots of synthetic gravity data made in six drill-holes. The data are the difference between the response of a density rock-model with the keel, and a model with the keel removed. Note the change in horizontal scale for plots C, E, and F 141 6.7 Synthetic magnetic data produced from a susceptibility rock-model of San Nicolas. Data are produced from models with and without the keel. The difference is also plotted 144 6.8 West-facing cross-section of the magnetic susceptibility rock-model' with the keel shown to the lower left. Four drill-holes are also shown, within which synthetic total magnetic field measurements are performed 145 6.9 Synthetic total-field, down-hole, magnetic data produced from a magnetic susceptibility rock-model of San Nicolas. Data are produced from models with and without the keel. The difference is also plotted. Drill-hole A is just north of the keel and drill-holes B,C, and D are each 50m further north that the previous one so that drill-hole D is 150m north of the keel. 146 xxi 6.10 A perspective view of conductivity rock-model (including keel), cut to show the conductivity structure at -400m south. The approximate resistivity values of the model are shown, as is the location of the transmitter that was used in the field CSAMT survey 149 6.11 Color plot of synthetic real and imaginary Ex and Hy fields, at the surface over the San Nicolas deposit, that emanate from a far, 16Hz, grounded transmitter. The first and second rows show the field components that would be produced over conductivity models with, and without the keel, respectively (these are plotted on the same scales). The third row of plots show the difference between fields produced with and without the keel for each component. Each plot covers the area over the deposit, as described previously, although the northing and easting values are not shown. . . . 150 6.12 Perspective views of conductivity rock-model showing locations of the transmitters used for forward modeling of electromagnetic data 153 6.13 North facing cross-sections through a rock-model, which includes the keel, showing the magnitude of current density in each cell for two, (a) far and (b) near, transmitter locations. The cross-section is located at -400m north and shows an increase in current density within the keel (the small anomaly to the lower left of each section) from about 5xl0-8 to 6xl0~7 A/m2 155 6.14 Real part of Ex surface measurements made for four different transmitter source locations over conductivity rock-models of the San Nicolas deposit with and without the keel 156 6.15 Imaginary part of Ex surface measurements made for four different trans mitter source locations over conductivity rock-models of the San Nicolas deposit with and without the keel 157 xxii 6.16 Synthetic, DC, apparent resistivity, Realsection data produced from con ductivity rock-models of San Nicolas. The x-locations refer to the mid point of the receiver dipoles, while the y-axis refers to different transmitter dipole positions; the bottom line corresponds to data collected with the longest transmitter, and the top line corresponds to data collected with the smallest transmitter length. The section is just a representation of the data and should not be taken as the true structure of the earth. The percent difference between data produced from models with, and models without the keel is also shown. Note: the apparent resistivity color-bar is on a log scale with low resistivities (high conductivities) displayed in red, while the percent difference is plotted on a linear scale 161 6.17 Schematic diagram of dipole-dipole array for one transmitter location, showing the pseudo-section plotting points relative to the electrode geometry. 162 6.18 Synthetic apparent resistivity dipole-dipole pseudo-sections produced from conductivity rock-models of San Nicolas. Pseudo-sections of data produced over conductivity rock models with and without the keel are plotted on a log scale along with the percent difference between the two, which is plotted on a linear scale. The plotting convention use is shown in Figure 6.17.163 6.19 Synthetic apparent chargeability Realsection data produced from charge-ability rock-models of San Nicolas. The x-locations refer to the midpoint of the receiving dipoles, while the y-axis refers to different transmitter di pole positions; the bottom line of data being collected with the longest transmitter, and the top line of data being collected with the smallest transmitter length. The difference between the two data sets is also plotted. 166 xxiii 6.20 Synthetic apparent chargeability dipole-dipole data produced from charge-ability rock-models of San Nicolas. The plotting convention used is shown in Figure 6.17. The difference between the two data sets is also plotted. . 168 7.1 Geophysical Exploration Methodology Flowchart 172 B.l The range of susceptibility values for different rock types at San Nicolas. 187 B.2 Coarse and fine magnetic susceptibility measurements on a section of core from drill-hole SAL 70 188 B.3 Coarse and fine magnetic susceptibility measurements on a section of core from drill-hole SAL 34 189 B.4 Magnetic susceptibility measurements made on core from drill-hole SAL 25. A composite model was produced that was used to constrain inversion of ground magnetic data 190 B. 5 Perspective view of the physical property conductivity model used for test forward modeling in Section 4.4.1. The model has been sliced at -400m N to expose a north-facing cross-section 192 Cl The force exerted on mass mi by mass m2 is in the direction of r 194 C. 2 Orientation of transmitter-receiver (Tx-Rx) pairs in the Dighem frequency-domain airborne electromagnetic system 198 C.3 Diagrammatic sketch of the Questem time-domain electromagnetic sys tem. The transmitter (Tx) is placed around the aeroplane while the three-component receiver (Rx) is towed in a bird behind the plarie. GPS positioning is used to measure the position of the receiver relative to the bird ... 200 XXIV C.4 The current in the transmitter as a function of time for 50% of one cycle. Receiver measurements are made over 15 time-windows during the 16 msec following the current on-time 201 C.5 Diagrammatic sketch of the CSAMT survey layout (modified from Zonge and Hughes) 202 C.6 "Real-section" array geometry used for collection of DC resistivity and IP data in Section 4.5 203 C.7 The voltage response as a function of time, V(t), in the presence of polar-izable material, when a current is suddenly applied and then interrupted. VT = VP + VS 205 D.l Topographic elevations at San Nicolas. Apart from a hill at 6000m east, -1000m north, most of the area exhibits small topographic elevation changes. Outlines of outcropping geology are also shown, along with a hammer and pick symbol that locates San Nicolas at about -1700m east, -400m north. 207 D.2 Ground gravity observations. The data are the original data supplied by Quantec but with observations on the hill removed. The topographic map shown in Figure D.l indicates a hill with over 100m of elevation gain. Initial inversions produced a large negative anomaly coincident with the hill. Observations made on and close to the hill were subsequently discarded, although the affects of the hill probably persist in some of the remaining observations in this area. A constant was also removed from the data in order to remove a first order regional field. An outline of outcrop geology is overlain 208 D.3 Dighem airborne total-field magnetic data. The deposit is identified by a subtle anomaly at -1600m east, -400m north 209 xxv D.4 Subset of Dighem airborne total-field magnetic data used for inversion. . 210 D.5 Falcon airborne total field magnetic data 21D.6 Questem airborne total field magnetic data 211 D.7 In-phase and quadrature frequency-domain electromagnetic data from the Dighem system: 466Hz, 1046Hz, and 6495Hz 212 D.8 In-phase and quadrature frequency-domain electromagnetic data from the Dighem system: 7206Hz and 56kHz 213 D.9 Profiles of Z-component dB/dT data for the fifteen time-channels (0 to 18 msec after the transmitter current turn-off) along line 10. The upper panel shows data for the first five time channels, the center plot shows data for the next five time channels, and the lower panel shows data for the last five time channels. The deposit corresponds to the strong negative response in the later time channels. Thirty-one lines of data were collected and different orientations but this is the only line that has been presented. 214 D.10 Gradient Array IP data. The boundary shows the extension of the survey measurements and an outline of outcropping geology is also shown. . . . 215 D.ll Reduced gravity data. This gravity data has had a regional signal removed. The resulting gravity data were used for producing the density-contrast models in Chapters four and five 216 D.12 Total-field ground magnetic data. The observation locations are displayed and the gaps where the data have been discarded are clearly visible. ... 216 D.13 Observed and predicted CSAMT apparent resistivity data. The data was collected along three lines, -200m N, -400m N, and -600m N. Each data plot is plotted with frequency in the vertical axis (ranging from 8192Hz at the top to 0.5Hz at the bottom), and easting along the horizontal axis (-2200m E to -800m E) 217 xxvi D.14 Observed and predicted CSAMT phase data. Each data plot is plotted with frequency in the vertical axis (ranging from 8192Hz at the top to 0.5Hz at the bottom), and easting along the horizontal axis (-2200m E to -800m E) 217 D.15 Observed and predicted "Realsection" IP data. Five rows of data are plotted in each panel. Each row corresponds to a different transmitting dipole length 218 F.l Synthetic chargeable body 227 F.2 North-facing cross-sections through the conductivity models: A) 0.01 S/m half-space, B) 0.1 S/m block in a 0.01 S/m background, C) 1 S/m conduc tive block in a 0.01 S/m background, D) 0.1 S/m block and a 40m thick 0.1 S/m overburden in a 0.01 S/m background 228 F.3 North-facing cross-section through recovered chargeability models 229 F.4 North-facing cross-section through vector model of real component of the electric field 231 F. 5 North-facing cross-section through vector model of real component of the current density 232 G. l Removal of electrode-effect via weighting 236 xxvii Acknowledgments I would like to thank my advisor, Dr. Doug Oldenburg, for his guidance, motivation, and enthusiasm for this project. His mentorship has made this thesis a very enjoyable and interesting study to be involved with. My sincere appreciation to former and present Geophysical Inversion Facility mem bers without whom this thesis wouldn't be what it is, specifically: Christophe for his endless patience and help with time-domain electromagnetic modeling; Roman for his amazing technical support; and Colin and Eldad for their ability to answer my questions at a simplified level of understanding. Other thanks go to Partha, Chad, Jiuping, Yuji, Francis, Yaoguo, Peter, Len, and Steve for their help, input, and friendship. My thanks to other members and former members of the Department of Earth and Ocean Sciences at the University of British Columbia who have made my time here fun and fulfilling, including: Christina, Chris, Nicolas, Kevin, Phil, Charlie, Tim, Dave, James, Daniel, Kim, Kate, Jeff, Jeff, Andrew, John, Fern, Andy, Stephan, Parisa, John, Steve, Matt, and of course Jane. Included in this acknowledgment are the members of the Mineral Deposit Research Unit. Also, my thanks go to the following companies for their support in this venture: New-mont Gold Company, for their encouragement and financial support (specifically: Jim, Rick, Perry, and Colin); Teck-Cominco, for providing an excellent data set with which to work; Mira Geoscience; Quantec Geoscience; Fugro; BHP-Billiton; SJ Geophysics; and Kennecott, all for helping with issues relating to the data sets. Finally, my thanks to Ursula and Ellen for helping English with his English. xxviii Dedication I would like to dedicate this thesis to my family who support me in exploring the riches of the world. xxix Chapter 1 Introduction Geophysical techniques can be used to detect and delineate mineral deposits. This has become more evident in recent years for the following two reasons: 1. The recently developed ability to invert geophysical data to recover three-dimensional physical property models. 2. The increasing need to explore in concealed terranes. As a result, geophysics is providing a larger amount of information than in the past, which can impact an exploration program. At this time, it is important to re-evaluate how geophysical inversion can be employed in an integrated exploration program that includes geological, geophysical, geochemical, and remote sensing components. The possible role of geophysics in an integrated exploration program is outlined by a flowchart in Figure 1.1. The aim of this thesis is to demonstrate the use of geophysical inversion techniques at various stages of the exploration program. As an example, the results of inverting geophysical data from the San Nicolas copper and zinc deposit, in Mexico, are presented. An exploration methodology is presented that outlines the steps that might be taken in an exploration program. Within this larger framework, the role of geophysics is defined. The geophysical role is an integrated part of the program with geologic and physical property information being essential for the optimal use of geophysics. This geophysical exploration methodology is followed throughout the thesis from defining an exploration model to delineating the deposit. 1 Chapter 1. Introduction 2 Exploration Goal: •find economic mineralization Geologic Constraints; •existing mapping , 'historical workings •tectonic setting Geologic 'descriptive model [ | Models: 'genetic model j •> Exploration Model I J ™1 \ Regional Geologic I j Mapping: j I 'Outcrop geology j j •Structure " j I 'Alteration \ Geochemistry •Stream sediments Remote Sensing •Multispectral 'Hyperspectral Geologic Interpretation No Targets 3* Target Generation Terminate Exploration Fro gram Update Geologic L Exploration Model j Geologic Logging 4 * Detailed Geologic Mapping: •Stratigraphy •Structure >, | •Alteration j Trenching Rock samples •Assays •Alteration mineralogy (X-ray/PIMA) Geochemistry •soil surveys .... Geologic 1 Interpretation Mineralization Not Intersected Locate Drill-hole Terminate Exploration Program Mineralization Intersected Goal Realized: Mineral Resource Defining the Geophysical Exploration Model Geophysical .cxpcctoJ sja/depth of target Exploration .^i^ propertfe, contrasts Model: Regional/Reconnaissance Scale Geophysics: Physical Property Measurements Existing Data Survey Design Collection I Gravity AMAO (+/- ARAD) Other VLF, MT Modeling . Individual Forward Correlation Inverse Joint/Simultaneous Forward Inverse Conclusions about physical property distributions of the the earth Update Geophysical Exploration Model Local/Deposit Scale Geophysics: Additional Physical Property Measurements Survey Design 3 1 • i • Gravity GMAG EM: CSAMT, UTEM DC/IP Other: VLF Modeling Individual Forward Inverse Joint/Simultaneous Forward Correlation Conclusions about physical property distributions of the the earth Physical property core measurements Down-hole methods z Delineation If Figure 1.1: Geophysical Exploration Methodology Flowchart: Defining the geophysical exploration model. Chapter 1. Introduction 3 In this thesis, sophisticated inversion techniques are used to generate digital, physical property models of the subsurface from geophysical data. These techniques are imple mented at the regional scale to identify additional areas that might host mineralization, and at the deposit scale in order to define the deposit and delineate mineralization. Physical property measurements made on drill-core have been integrated in the inver sion modeling at the delineation stage. The models are more easily interpreted than the data, and can be related to geology using physical property measurements of local rocks. Furthermore, different inversion models are combined to establish a more comprehensive understanding of the subsurface. Interpretations of the inverted models are proposed based on the geologic information presented in the thesis. In addition, the San Nicolas deposit is almost entirely bounded to the east by a southwest-dipping fault, which could have been a feeder structure. Mineralization con tinues to follow the fault at depth in an unconstrained part of the deposit called the "keel". Modeling has been performed to demonstrate the sensitivity of different survey designs to this deep high-grade ore-zone. The remainder of this chapter introduces: (1) the exploration methodology that is followed throughout this thesis, (2) the geologic background of San Nicolas, which is used to produce a geologic model, (3) the physical properties of rock units at San Nicolas, (4) the geophysical exploration model, (5) the data that have been collected, and (6) the tools that are used to model the geophysical data. 1.1 Exploration methodology Mineral exploration programs commonly employ a staged approach when assessing large amounts of land or data. The aim of this approach is to systematically discern the most Chapter 1. Introduction 4 prospective targets from a large amount of land without allowing any economic miner alization to be overlooked. These targets can then be evaluated with a local exploration program and drill-testing. The flowchart in Figure 1.1 starts with an exploration goal and proceeds with defining a geologic exploration model on which a regional scale exploration program is based. The geologic exploration model provides exploration criteria, which are an estimate of the essential controls on mineralization (favorable stratigraphy, association with faults, etc.). These criteria are based on descriptive models (including mineralogy) and possibly genetic models that provide information about how the deposit might have formed. Also the model can be defined from knowledge about existing deposits in the area, and their tectonic setting. Geophysical exploration methodology requires a geophysical exploration model. This builds on the geologic model and, with the inclusion of physical property information, defines geophysical exploration criteria. The criteria provide the expected physical prop erties of the mineralization itself, or physical properties of critical geologic features asso ciated with mineralization. At the regional scale the aim is to identify smaller, more focused, areas of interest. Different geophysical techniques can be used with survey designs that are tailored to the geophysical exploration model. Modeling of the geophysical data can produce inferences about the physical property distribution of the earth that, along with physical property information, will aid in regional geologic interpretations. The regional program may result in generation of targets that warrant more detailed investigation. As shown in the flowchart, target generation is followed by the local- or deposit-scale exploration program. The local geophysical program will require different survey techniques and designs to those used at the regional scale, and this will be based on an updated geophysical exploration model. At this stage, geophysical inversion can play a Chapter 1. Introduction 5 role in defining more detailed aspects of the subsurface, such as the depth and size of a body of favorable physical properties, or the geometry of deposit-critical geologic features. Along with physical property information of the local rocks, geophysical inversion models contribute to local-scale geologic interpretation. This geologic interpretation can result in the spotting of drill-holes. The iterative nature of an integrated exploration program can continue during drilling. Information from drill-holes or drill-cores can be used to further refine geologic and geo physical exploration models. This will result in better inversion modeling of the physical property distributions. Subsequently, these models can better guide drilling, which could result in defining mineralization with fewer drill-holes. Furthermore, surveys can be de signed to delineate specific areas of mineralization that may have been intersected, but not completely defined by drilling. This exploration model, while only outlining two stages of the exploration process, allows for several iterations and further focusing of a program through the use of feedback paths. Within the flowchart, emphasis is placed on how geophysics and inversion is integrated into such a program, but this is not meant to reflect a reduced importance of the other methods. It is also acknowledged that the exploration model may be over-simplified; however, it does demonstrate how geophysics has been used at San Nicolas. 1.2 San Nicolas San Nicolas, owned by Teck Corporation and Western Copper Holdings Ltd., is an un-mined, volcanic-hosted, massive sulphide deposit containing ore-grade copper and zinc with associated gold and silver. It is located in Zacatecas State, Mexico (Figure 1.2), a region that, until the discovery of San Nicolas and the nearby El Salvador deposit, Chapter 1. Introduction 6 was unrecognized as having potential to host such deposits [Johnson et al., 2000]1. The reserve estimates are 72 million tonnes grading 1.4% copper and 2.27% zinc, making it the largest massive sulphide deposit yet discovered in Mexico, and a world class deposit. 1.2.1 Geology Discovery of economic ore is of primary importance. Although we hope to detect mas sive sulphide mineralization using geophysics directly, the geologic environment must be considered in order to effectively use geophysical information. An understanding of the re gional geology can focus exploration efforts in more favorable structural and stratigraphic settings, thus expediting the whole exploration process. The association of mineraliza tion with geologic features that could be identified using geophysics, such as: favorable horizons, plutons, faults, shallow intrusions, rhyolite domes, and alteration zoning, aid in producing informed exploration decisions. Discrimination of geologic features that could be mistaken for sulphide ore is also an important part of the exploration process. Fur thermore, knowledge of the local geology, and the associated physical property values, allow for expected geophysical contrasts to be determined, and appropriate geophysical methods to be applied. The geophysical data that have been collected is not only used for direct detection of massive sulphide ore, but for providing information that can help form a geologic model of the area. A background to the geologic setting and deposit genesis is given below, while more detailed regional- and local-scale geologic information is given at the beginning of Chap ters 3 and 4 respectively. 1The area has been known mostly for the historically worked, epithermal silver and gold deposits. Chapter 1. Introduction K .'I w Pacific MEXICO UNITED STATES Gulf of Mi Zacatecas State 1 ; Mexico City*'. ' 'ft BELIZE \ fUATKMAlf & S9ft8 Natie*jai Geograpiwc Society Zacatecas 0 50 m! 8i ' tt SO km OAlitlllA . • SOOI National Caographir Socsty DURANCO T AM4UI ;?AS IANAJUATO •" •' Figure 1.2: Map of Mexico (upper panel) and Zacatecas State (lower panel) showing location of San Nicolas. Chapter 1. Introduction 8 Geologic setting The paleotectonic regime in which the San Nicolas deposit formed is a transition from a back-arc basin to an arc, and possibly a fore-arc setting2, as documented by Daniel-son [Danielson, 2000]. The deposit is hosted in marine volcanic and sedimentary rocks, locally known as the Chilitos Formation, of Upper Jurassic to Lower Cretaceous age, [Johnson et al., 2000]. Mineralization is predominately found along a time-line that oc curs between the formation of rhyolitic lava domes and the deposition of mafic extrusives above. These rocks include: mafic, intermediate, and rhyolite flows and shallow intru sions, felsic tuffs, volcanoclastic sediments, mudstone, chert, and limestone. This package of rocks was subsequently accreted to North America; possibly in the mid-Cretaceous. Northeast-directed thrusting of the Late Cretaceous to Early Tertiary resulted in the juxtaposition of the Chilitos formation against carbonates and siliciclas-tic sedimentary rocks. Throughout the region, the rock assemblages are intruded by Late Cretaceous to Early Tertiary granitic plutons, and are unconformably overlain by Tertiary aged felsic volcanic flows and tuffs, [Johnson et al., 2000]. North-, northwest-, and west-striking faults cut the assemblage; the age and displacement of these faults are unknown, although they could be synvolcanic extensional or postvolcanic strike-slip faults, [Johnson et al, 2000]. A section of rocks, which includes the Chilitos Forma tion, has also been uplifted to form a horst, within which the San Nicolas deposit if found. Tertiary volcanic rocks have intruded along the northwest-striking faults, and along north-striking faults at the west of the horst. The deposit, itself, is bounded to the southwest by the northeast flank of a rhyolite 2 A back-arc to fore-arc setting refers to relative locations within an island arc system. Island arcs are a chain of islands, such as the Aleutian Islands, that form at the convergent margin of two oceanic plates. They are often associated with deep oceanic trenches, and volcanic activity. A fore-arc setting would exhibit compressional tectonic forces, while a back-arc setting exhibits tensional tectonic forces. Chapter 1. Introduction 9 dome, and by a southwest dipping fault to the northeast. The deposit is copper-rich towards the bottom with a high-grade polymetallic cap at the top. Also, the sulphide is granular and fragmented in the lower zone and laminated in the upper zone, typical of VMS deposits. The two zones could have different physical properties that geophysics could help differentiate. Silicification and barite are present along the fault that bounds the deposit to the northeast. Extensive hydrothermal alteration of the host rocks is also observed; this could affect physical properties of the rocks. Deposit Genesis Massive sulphide deposits are often formed in the proximity to plate margins [Lydon, 1988]. Mineral-rich hydrothermal fluids, either originating from a magmatic source or circulat ing in the sediments and being driven by the heat from a magmatic source, react with sea water and deposit sulphide minerals. These deposit types are usually quite small and the duration of mineral deposition short-lived. Several deposits can be found in the same area often in association with submarine volcanic flows and sediments [Lydon, 1988]. San Nicolas, although probably formed by a similar process, is larger than usual and was possibly formed in a tectonic environment not usually associated with massive sulphide deposits. This makes San Nicolas an unusual massive sulphide deposit, and provides a new and exciting type of exploration target for this region. Consequently, the regional potential for discovering other deposits is high because exploration programs in this area have not previously focused on this type of deposit. Geologic Exploration Model The main features of the deposit that could be critical to similar deposits being found comprises the geologic exploration model. At San Nicolas this could include the following: Chapter 1. Introduction 10 1. Favorable host stratigraphy between the rhyolitic and mafic volcanics that corre sponds to a time-line 2. Association of the deposit with a quartz-rhyolite dome complex. 3. The shallow depth of the Chilitos formation due to the presence of a horst. 4. Faults that may have provided conduits for hydrothermal fluids. 1.3 Physical Properties Geophysics can help define the physical property distribution of the subsurface. This is useful when objects can not easily be investigated firsthand, and have physical prop erties that distinguish them from their surroundings. Estimates of physical property measurements are also needed to relate inversion results to geology. Massive sulphides are commonly known to have contrasting values in density, conductivity, chargeability, and possibly magnetic susceptibility. These physical properties stem from the minerals which constitute the ore. In the case of San Nicolas, the ore is primarily made up of pyrite, chalcopyrite, and sphalerite, with some magnetite (not a sulphide mineral, al though frequently associated with sulphide assemblages), and possibly pyrrhotite. The total sulphide content is greater than fifty percent by volume within the main ore body. From the mineralogy described above, the deposit is expected to have the following physical attributes: 1. High density (a minimum of 3.8 g/cc [Ward, 1966]) 2. High magnetic susceptibility (due to the presence of magnetite and pyrrhotite) 3. High conductivity (this variable property often depends on the connectivity of the Chapter 1. Introduction 11 Rock Type Density Magnetic Susceptibility Resistivity Chargeability (g/cm3) (S.I. xlO"3) (Ohm-m) (mV/V) Tertiary Breccia 2.3 0-5 20 10-30 Mafic Volcanics 2.7 2 80 30-80 Sulphide 3.5 10 20-30 200+ Quartz rhyolite 2.4 0-1 100 10-20 Graphitic Mudstone 2.4 0 100+ 30-70 Table 1.1: Estimated physical properties for the five major rock types encountered at San Nicolas. The massive sulphide has contrasting values in all physical properties apart from resistivity. A low resistivity value of 20 Ohm-m for sulphides, similar to the Ter tiary breccia, makes detection of the deposit more difficult when using some geophysical methods. conducting paths3) 4. High chargeability due to the presence of metallic minerals [Telford, Geldart, and Sheriff, 1990c] Apart from the sulphide ore, the geological units that are of significant interest are: the Tertiary overburden, which is conductive4 and contains some susceptible units; a mafic volcanic unit that is magnetically susceptible; and dense mafic intrusions. Ta ble 1.1 provides physical property information for the major rock types found at San Nicolas5, while Appendix B contains a complete review of physical property values used to constrain the models and interpret the inversion results. 3many sulphides, excluding sphalerite, are semiconductors with conductivities often above 1000 S/m in their pure form [Sherman, 1997] Conductivity (CT [Siemens per meter]) is the inverse of resistivity (p [Ohm-meters]). 5The density of the sulphide has been reduced due to the presence of semi-massive sulphide material. Chapter 1. Introduction 12 1.4 Geophysical Exploration Model Based on the physical property values and the geologic exploration model we can define a geophysical exploration model (Figure 1.1). For San Nicolas this can be defined as follows: 1. Bodies with contrasting density, chargeability, conductivity and magnetic suscep tibility will be good candidates for mineralization. 2. In combination, the physical property contrasts can be used to differentiate between mineralization and shallow intrusives that may be just dense and magnetic. 3. The density contrast between the host volcanic rocks and the overburden will help define the horst, which is important in bringing the favorable horizon near to the surface. 1.5 Data Many geophysical data have been collected over San Nicolas, and the surrounding area. Data were collected at San Nicolas during the last five years, and provided to U.B.C. by Teck Corporation. The data used in this study include: airborne electromagnetic and magnetic data, gravity, ground magnetic, CSAMT (controlled source audio-frequency magnetotelluric), IP (induced polarization) data. The locations of these data sets are shown in Figure 1.3. Plots of the data are presented in Appendix D along with a summary of survey specifications for each data set (Table D.l). In order to limit the number of plots in the main text I present only: (1) physical property inversion models, (2) geologic plan maps and sections, and (3) interpretative figures. Geological information, as interpreted from drill holes, is also used to demonstrate the success of inversions performed, and imply the benefits of inversion methods in exploration. Unless stated otherwise, all coordinates Chapter 1. Introduction 13 are given in meters and are denned on a local San Nicolas grid system. The conversion between the local grid system and the UTM Zone 13, NAD 27 system that was employed, is given in Appendix A. 1.6 Modeling Geophysical modeling involves relating, or mapping, data (a numerical representation of a physical signal that has been measured or observed) to a numerical representation of a physical property distribution within the earth, known as a model. The mapping procedure can occur in both directions. 1.6.1 Forward modeling Forward modeling is a procedure in which a unique set of synthetic geophysical data are calculated based on a known physical property distribution. As well as being an inherent part of the inversion process, forward modeling plays an important role in designing geophysical surveys, and testing geologic scenarios. 1.6.2 Inverse modeling Inversion is the process of constructing physical property models of the earth based on information contained in geophysical data. Unlike forward modeling, the inversion prob lem is often ill-posed. An infinite number of models could have produced the data. We are able to choose one model by only considering those with simple structural charac teristics. Inversion models are generally much easier to interpret than the original data and provide a superior understanding of the subsurface. A description of the inversion procedure, as it applies to this thesis, is given in Chapter 2 and the inversion parameters used to generate the physical property distributions presented are given in Table E.l. Chapter 1. Introduction 14 1.7 Outline of the thesis The following provides a synopsis of each chapter. Having already defined the geophysical exploration model, the following chapters (apart from Chapter two, which introduces geo physical inversion) proceed through the geophysical exploration steps while emphasizing the role of inversion. The work is presented within the context of San Nicolas. 1.7.1 Inversion methodology Chapter two provides an overview of the inversion process. This is needed in order for the inversion models to be given their correct context regarding the information they supply to the geologist. 1.7.2 Regional-scale modeling The regional scale for this thesis is defined by the coverage of gravity and airborne magnetic and electromagnetic data. This large-scale region is shown by the upper plot of geology with data boundaries in Figure 1.3. Not all of the airborne magnetic and electromagnetic, and reconnaissance IP data sets cover such a large area as the regional-scale coverage, so an intermediate-scale is introduced and shown in the center plot in Figure 1.3. Inversions presented in Chapter three are performed at either the regional or intermediate scales. In the context of this thesis, only inversion models of geophysical data are used to answer relevant regional geologic questions and focus exploration efforts. 1.7.3 Deposit-scale modeling This is the scale at which anomalies are explored and deposits are found, or grossly delineated. More detailed surveys are carried out on the best areas and targets that were identified at the previous stage of the exploration process. It may be, that work carried Chapter 1. Introduction 15 15000 10000 5000 e o -5000' -10000. -15000 ^ \ \\ \ \ \ \ * \ .A f—- ' M \ \ \ V Local Grid North UTM Grid North -1(1000 -5000 0 5000 ; Easting (m) 10000 15000' 5000 E cn c. E o -5000 -10000 Regional Scale Gravity Data Dighem Electromagnetic and Magnetic Data Intermediate Scale Questem Electromagnetic and Magnetic Data Falcon Magnetic Data Subset of Dighem Magnetic Data Gradient Array IP Data -5000 Easting (m) 5000 Deposit Scale Ground Magnetic Data Realsection IP\DC resistivity and CSAMT Data X San Nicolas -2500 -1500 Easting (m) Figure 1.3: Data coverage over the San Nicolas deposit and the surrounding area. Top: the regional-scale area with data coverage and outcrop geology (legend shown). Center: Intermediate-scale data coverage shown on top of outcrop geology. Bottom: deposit-scale data coverage. The size of the rock-hammer symbol approximately corresponds to the extents of the surface projection of the San Nicolas deposit. Chapter 1. Introduction 16 out on a more regional scale is sufficient to identify a discrete target, but more detailed information about the deposit is needed to determine whether it is worth drilling with one or five drill-holes, or at all. Alternatively, drilling may already have commenced and the geologist may require information about the geometry of the deposit before drilling the next hole. These scenarios provide the motivation behind the deposit-scale inversions of gravity, magnetic, electromagnetic, and induced polarization data that are presented in Chapter four. This work is performed at the scale depicted by the lower plot in Figure 1.3. 1.7.4 Detailed modeling using drill-hole information Chapter five addresses the use of a priori information in refining the physical property models constructed at the local-scale. Density and magnetic susceptibility measurements provide information about the earth that is independent of the data. Including this information in the inversion process forces only models that are consistent with known physical property distributions to be constructed. Such a priori information could be used to help spot a second drill-hole based on the results of the first, at the early stages of drilling, or could help model what is below a drilled out portion of a deposit at the end of the first phase of drilling. Furthermore, inversions with a priori information could be used to incrementally update physical property models with new information as each drill-hole is completed and logged. The evolving physical property models could help guide a drill program that would potentially result in less drilling needed to delineate the mineralization or conclusively test the anomaly. 1.7.5 Detection of deep ore: survey design Forward modeling is an important tool that generates synthetic geophysical data based on physical property distribution models. Chapter six employs forward modeling to examine Chapter 1. Introduction 17 the success (based on signal-to noise estimates) of different survey designs in detecting the keel and any extensions of mineralization that may exist. Forward modeling can also be used to test the geophysical response of different geological scenarios against field data. 1.8 Summary The flowchart of Figure 1.1 provides the framework around which this thesis is based. The flowchart is presented at the beginning of each results chapter in order to set the scene and provide an exploration context for the results. Equipped with knowledge of the geologic environment, physical properties, and in version methods, geophysical inversions are performed on data from San Nicolas. The inversion results demonstrate the assistance that geophysical inversion can provide to geologists regarding a variety of issues that range from regional exploration to detection and delineation of deep ore. Understanding the geology, and providing physical property models with a geologic context, are essential at all stages of the exploration process and are important to this study. Geologic information is also used to determine the success of using geophysics alone to define the deposit. Chapter 2 Geophysical Inversion 2.1 Introduction This chapter aims to construct a framework of understanding about the inversion process so that important issues such as non-uniqueness and data misfit are appreciated, and the inversion models are given their correct context regarding the information they supply to the geologist. Geophysical inversion produces physical property models from geophysical data, whereas forward modeling produces data from a physical property model of the earth. Several factors influence the inversion procedure and, therefore, the outcome of the model. These factors include: (1) estimating data errors, (2) using different reference models, (3) choosing smoother models or smaller model types, and (4) choosing a different balance between honoring the data and introducing features in the model. This chapter will briefly address the aforementioned issues, along with providing a synthetic example, in order to demonstrate the construction of the physical property models presented in this thesis. 2.2 Inversion Basics 2.2.1 Data What are the data? This may seem like a basic question, but it is one that needs to be explicity answered in order for an inversion to be successful. Data include not only the actual datum values, such as the strength of a magnetic field or the voltage between 18 Chapter 2. Geophysical Inversion 19 two electrodes, but, an error associated with each datum (where possible), locations of transmitters and receivers, orientations of instruments, data normalizations or changes in units, and detailed knowledge of any processing that is applied to the data. These are all crucial elements that define what the datum is and how it is connected to the spatial distribution of physical property values within the earth. Other information that can be thought of as a form of data could be physical property measurements made directly on the earth or any geologic knowledge that is pertinent to the region of the earth that produced the geophysical data. In the context of this thesis I deem such information to be a priori, that is information known prior to inverting the geophysical data. We can use these data to constrain and guide inversions. 2.2.2 Discretization and forward modeling The earth is a continuous medium having physical property values everywhere. We can describe the physical properties of the earth as a continuous function of position, m(x, y, z). When modeling, we usually simplify the earth through discretization, that is, dividing the earth into many cells whose physical property values are constant and whose size is smaller than the resolving power of the experiment or survey. For the inversions presented in this thesis, the cells are cuboidal and aligned with an orthogonal coordinate system. A 3D grid system or mesh defines the discretization of the model (Figure 2.1). The earth is now represented by a column vector m = (mi, m%,... , mM)T of length M (the number of cells), where each element of the vector corresponds to a the physical property value of a different cell in the model. Forward modeling then consists of solving a system of equations to predict the responses at each observation location based on the physical property distribution of the discretized earth. For a linear system of equations, the data are calculated by performing a matrix-vector multiplication. The matrix notation of the discretized forward problem is stated Chapter 2. Geophysical Inversion 20 d = Gm, (2.1) where G is the forward operating matrix (of size N x M), m is a column vector com prised of M model parameters, and d is a column vector of length N, containing the observations, or data. For the inversion problems addressed in this thesis, the number of data N is much less than the number of model parameters M. 2.2.3 Non-uniqueness and the Model Objective Function In the inverse problem we have the data (d) and a forward operator (G), and we want to find the model (m). However, as with any data set, we are limited to dealing with a finite number of inaccurate observations. This means there are many more model parameters than data, which makes the system under deter mined. Also, we don't want to produce a model that will exactly reproduce the observed data, as this will also reproduce noise. Therefore the data don't constrain all of the model parameters uniquely, and there are an infinite number of models that, when forward modeled, will reproduce, or fit, the data to a specified degree. Chapter 2. Geophysical Inversion 21 In order to chose one model from the infinite many that fit the data, we select the type of model we desire using a model objective function (sometimes called a model norm). The inverse problem is posed as an optimization problem by minimizing the model objective function subject to the requirement that the recovered model produces data that adequately resembles the original observations. When a model objective function is suitably designed the recovered model has simple structural characteristics that should emulate geology. Flexibility should be built into the model objective function so most simple geologic situations can be accommodated. The model objective function for a continuous 3D model rn(x, y, z) can be written as: = as J wa(m(x,y,z) - m(x,y, z)ref)2dv + ax J wx ^ dm(x,y,z)\2 dx 2 f fdm(x,: 2 dy ) dv + a*Jv w*{—z with the equivalent matrix notation written as: + oyj^ ^^1) dv + az[Wz() dv, (2.2) $m= ||Wm(m-rnre/)||2. (2.3) The first term in equation 2.2 measures the closeness of the model to a reference model (m(x, y, z)ref), while the following terms measure the smoothness of the model in the x, y, and z directions. Coefficients as, ax, ay, and az are used to weight the different parts of the model objective function relative to each other. This enables the model objective function to recover models that are close to the reference model, or preferentially smooth in one direction, by only changing one or two parameters. Spatially dependent weighting functions ws, wx, wy, and wz are used to weight the importance of one model parameter versus another. The reference and weighting functions allow for a priori information to be incorporated into the inversion that may come from other geological or geophysical sources. Bounds can be implemented to limit the range of values a cell can exhibit and Chapter 2. Geophysical Inversion 22 dip information can also be included [Li and Oldenburg, 2000b]. Equation 2.3 is the matrix notation of the model objective function, where Wm is a weighting matrix of size M x M that includes all the weighting and derivative functions of the model objective function, and || • || is the l<i norm. 2.2.4 Signal-to-noise and data misfit Signal-to-noise plays an important part in inversion because, no matter how good the field work and instrumentation, the data always contain noise. It follows that we do not wish to produce a model that can exactly reproduce the data as it would also reproduce the noise (known as over-fitting the data). Unfortunately noise levels are often unknown or only approximately known, but some estimate as to either absolute or relative noise levels is usually required. In order to measure how well the model response, or predicted data, fit the observations, we design a data misfit function as stated below. where N is the number of data, d°bs is the ith observed datum, dfre is the ith predicted datum generated by forward modeling the model, and represents the standard deviation of the ith datum, and is included in when the expression is written using matrix notation (Equation 2.5). We assume the noise is Gaussian in nature, has zero mean, and is non-correlative. This being said, if noise levels are well estimated, the data misfit should be approximately equal to the number of observations « N). When data have been over fit « N), the misfit is too small, and models often display erroneous structural complexity. This complexity is needed to fit noise in the signal and has no reflection on the structure of (2.4) Wd(d - dpre) (2.5) Chapter 2. Geophysical Inversion 23 the earth. If the data misfit is too large (<3>d » N), the model often lacks structure signifying that not all of the information that the data contain about the earth has been extracted. In practice, several inversions are generally needed in order to produce a model that has some structure, as defined by signal in the data, but not too much structure or complexity, as originating from noise in the data. Estimating the correct noise levels often comes down to the experience of the user. Obviously, the better the signal-to-noise ratio and the noise estimates of the data, the better the inversion model. 2.2.5 Optimization The inversion procedure is stated as the following optimization problem: find a model that minimizes the model objective function subject to fitting the data to a specified tolerance. The objective function is comprised of the data misfit and the model objective function (Equation 2.6), with (3 representing a trade-off, or regularization, parameter. Minimization of the objective function will generate the model. Choosing the regularization parameter is key to producing a model that has the max imum amount of information about the earth while containing the minimum amount of erroneous structure that relates to noise in the data. When noise levels are known, the reg ularization parameter can be determined using the discrepancy principle [Hansen, 1998]. Methods that can be applied when noise levels are unknown include: (1) generalized cross-validation (GCV), which is based on the principal that a good model would be able to predict a missing datum1 [Haber, 1997], (2) the L-curve criteria that searches 1For the generalized cross validation method to be most effective, a sense of the relative size of data errors is necessary, although the actual magnitude of the error is not needed. (2.6) Chapter 2. Geophysical Inversion 24 for the corner of the curve when data misfit and model objective function values are plotted for different values of the trade-off parameter2, both of which are described by Hansen [Hansen, 1998], and (3) manually performing inversions with different noise es timates. In practice it is often useful to apply GCV or L-curve methods to get an idea of the approximate magnitude of the errors, and then use the discrepancy principle to slightly increase the misfit. 2.2.6 Appraisal When a model is generated via the inversion procedure it is often desirable to determine how good the model is, that is, how well does it represent the earth? As geophysics is usually a minimally invasive technique we don't tend to dig up the earth, but revert to the data and model to answer this question. What is really being asked is: what can be concluded about the earth, given the data? As it is already known that the constructed model is not the only one that will fit the data, a range of models that fit the data could be generated. Construction of many models using inversion could be quite time intensive, although a rudimentary exploration of this model space is a good idea. This can be accomplished by changing the model objective function. The user soon gets an idea of features in the model that don't change (which are therefore probably required in order for the model to fit the data) and features that do change (which are not constrained by the data). It is important to be aware that the model may not define the earth well in all places, and over-interpretation of features that are not constrained by the data is 2On a log-log scale, a plot of the data misfit (<£<,) verses model norm (4>M) for different values of the trade-off parameter (j3) produces a curve which resembles the letter L. In the upper part of the curve, a decreasing value of /3 will produce large decreases in data misfit for small increases in model norm. Conversely, along the lower part of the curve, as (3 decreases further a small decrease in data misfit will correspond to a large increase in model complexity, as measured by the model norm. A model that corresponds to the point of maximum curvature will have a reasonable amount of complexity but should not over fit the data. Chapter 2. Geophysical Inversion 25 possible. As an alternative to defining how good the model is, hypothesis testing can be performed to confirm, or eliminate, geological scenarios of interest. Other, more sophisticated methods of appraisal, which are not used in this thesis, might include analysis of averaging functions [Backus and Gilbert, 1968] and point spread functions of the model resolution matrix [Menke, 1984]. 2.3 A synthetic example In order to demonstrate the inversion process, and the level of understanding that is needed in order to gain the most from an inversion result, a synthetic example has been devised. A density-contrast model is used to forward model synthetic, vertical-component, gravity data. The data are contaminated with noise and then inverted using GRAV3D (Appendix E.2) to recover a density-contrast model. g/cm3 Figure 2.2: Perspective view of synthetic density-contrast model. The model has been volume rendered so that cells with density contrasts less than lg/cm3 have been removed. Chapter 2. Geophysical Inversion 26 2.3.1 The data To start, the data are generated by making measurements over a three-dimensional den sity contrast distribution (Figure 2.2), much like we would in the field. This is performed using a forward modeling procedure that is described in Appendix Cl. Figure 2.2 shows a volume-rendered, perspective view of the synthetic density contrast model. The model consists of 7200 cells of size 10m x 10m x 5m. The feature we wish to model is a 40m x 40m x 40m cube located 20m below the origin with a density-contrast of 1 g/cm3. The cube is hosted in a background of 0 g/cm3 density-contrast. Observed Gravity Data 400 data Easting (m) Figure 2.3: Surficial synthetic gravity observations with 3% + 0.0005 mGal random noise added. Four hundred synthetic gravity measurements were made at the surface on a regular 10m grid. The data were contaminated with 3% + 0.0005 mGal random noise (Figure 2.3) and standard deviations were assigned equal to 3% of the magnitude of each datum value. The minimum standard deviation was limited to 0.0005 mGal. Chapter 2. Geophysical Inversion 27 2.3.2 Presentation and interpretation of inversion models. Mineral deposits and the geologic structures can exhibit very three-dimensional geome tries. The geology of San Nicolas and the surrounding area is no exception, and much confidence is gained when we are able to model the earth in three-dimensions. Models often contain complicated 3D structures of varying physical property values that are diffi cult to present in a document. In this thesis 3D models are presented using a combination of perspective view images, plan views at different elevations, and cross-sections. Visu alization tools such as: coloration based on physical property values, shading, volume rendering, and view angle are factors that can change what is perceived when viewing a model. As a result, an attempt is made to objectively convey the salient features of each model with one or two images. For comparison between models, north-facing cross-sections at -400m north, and west-facing cross-sections at -1700m east are frequently presented. Furthermore, in order to address non-uniqueness issues, information may be distilled from the analysis of several models generated using different objective functions. Care is needed when making inferences about geology from an inversion model, and when the models are discussed in the text the reader may assume the features have been appraised using the methods outlined in Section 2.2.6. 2.3.3 Recovered models The inversion problem is solved by minimizing Equation 2.6. The solution depends on the definition of the model objective function ($m), the data misfit and the regularization parameter (B) that controls the trade-off between misfit and structure. The following sections show how changing the desired misfit or alpha coefficients in the model objective function effects the recovered model. Examples of changing the reference model and assigning bounds on some of the model cells are also shown. Chapter 2. Geophysical Inversion 28 g/cm3 Figure 2.4: Perspective view of recovered density-contrast model. The model is viewed from the southeast (local coordinate system), and is volume-rendered with a cut off' at 0.3 g/cm3. Figure 2.5: North-facing cross-section through recovered density-contrast model with location of true model shown. Chapter 2. Geophysical Inversion 29 A generic model For this first example the errors are assigned appropriately and so we have an estimate of the desired misfit. A model is recovered that fits the data to a value equal to the number of data = TV). Figure 2.4 shows a perspective view of the model recovered. The model has been volume-rendered using a cut-off value of 0.3 g/cm3. The view angle is the same as that used in Figure 2.3, but the coloration has been changed. A north-facing cross-section through the model is shown in Figure 2.5 with the outline of the true model also displayed. As can be seen from these two images, the recovered model has a concentration of anomalous density coincident with the location of the true model. The recovered model seems smeared out relative to the true model, especially at depth, producing an anomaly of larger volume, and, as a result, lower density contrast values. This is because the model objective function was designed to generate a model that is close to a reference density-contrast model of 0 g/cm3, and have smooth features. The density values of the recovered model are smaller than the true model; however, the total mass within the model is the same. Effects of over- or under-fitting the data Figure 2.6 shows north-facing cross-sections through the true model and seven recovered density contrast models from inversions of the synthetic gravity data. Each model is overlain with an outline of the true model, and all of the models are plotted on the same colour scale (0 to 0.5 g/cc). For comparison, Figures 2.6(a) and 2.6(b) show cross-sections of the true and recovered model with $d = TV respectively. Figure 2.6(c) shows a model chosen with the data misfit one tenth of the number of data. Because the correct noise levels were assigned to the data, the model over-fits the data. In comparison to Fig ure 2.6(b), where the data misfit was appropriate, the density-contrast values are slightly Chapter 2. Geophysical Inversion 30 (a) True density-contrast model (b) Model fits the data: $d = N (c) Model over fits the data: $d = O.liV (d) Model under fits the data: $d = 5N (e) Model smoother in z-direction:az > ctx,cty (f) Model smoother in horizontal direction: ax,ay > az (g) Reference model is half of true model (h) Half of model bounded Figure 2.6: North-facing cross-sections through true model and seven recovered mod els for seven different inversions. All of the inversion models are shown on the same color-scale with dark blue representing 0 g/cm3 and below, and magenta representing 0.5 g/cm3 and above. The original model, an outline of which is shown on all the recovered models, is shown at a different scale. Chapter 2. Geophysical Inversion 31 closer to the value of the true model. However, erroneous artifacts, represented by more complex structures, are produced in the recovered model as the noise in the data mani fests itself. Conversely, the density contrasts are lower when the data misfit is too large, as shown by the model in Figure 2.6(d). This model was produced with a misfit of five times the number of data. The under-fitting arises because the model objective function has more control over the model produced when the data misfit is larger, resulting in a smaller and smoother model with more simple structure. If noise levels were truly unknown then both of these models would be an equally viable a representation of the earth. In such cases determining how well the data have been fit becomes more difficult. GCV and L-curve methods can be useful in such cases as the models are chosen based on methods other than data misfit. Effects of changing the model objective function Figures 2.6(e) and 2.6(f) show the results of altering the model objective function to pro duce models smoother in the vertical direction and smoother in the horizontal directions, respectively. This is done simple by changing the relative values of the a coefficients in equation 2.2. As can be seen, the model shown in figure 2.6(e) is elongated in the vertical direction, which extends the anomaly at depth, while the model shown in fig ure 2.6(f) is truncated at depth but has broader horizontal extents. Both models still have a concentration of mass in the vicinity of the true model. The models also both fit the data to the same degree ($d = N) so, in the absence of any external biases, both are an equally reasonable representation of the earth. This is a demonstration of non-uniqueness. Together, the models produced represent a very brief exploration of model space. Chapter 2. Geophysical Inversion 32 Including a priori information As mentioned in section 2.2.1, other forms of information (such as that from a drill hole) can be incorporated into the inversion. If we assume that half of the cube is known to exist prior to inverting the data, we can use this information directly in the inversion. In the first case a reference model is constructed that is 1 g/cm3 where the cube is known to exist, and 0 g/cm3 everywhere else. The result of using the reference model with the same parameters (data misfit and a coefficients) as those used in the original inversion (figure 2.5) is given in figure 2.6(g). The recovered model is slightly improved compared to the one generated without the reference model in that the anomaly doesn't have values extending to such depths. This result could be changed further by altering the relative weighting of each cell to the reference model using the weighting matrix in equation 2.3. This would allow regions of the model that are known prior to inverting the data to be weighted more strongly with respect to the reference model than other cells. An alternative method is to restrict the range of values a cell can assume. In figure 2.6(h) a model was constructed using this method. For the part of the cube known to exist, the model was restricted to ranges between 0.9 and 1.1 g/cm3, for the rest of the model space the bounds were set from 0 to 1 g/cm3. Most of the mass is concentrated in the correct place. It seems evident that this technique has produced a very good model that is the most representative of the true model. This demonstrates how constraining the model in one area can improve the model in the other areas. 2.4 Conclusions Inverting geophysical data to construct physical property models can allow for greater understanding of the subsurface. However, in order to get the most from the inversion process the following issues must be appreciated: (1) the data have to be fully understood Chapter 2. Geophysical Inversion 33 and, where possible, errors should be assigned; (2) fitting the data to an appropriate de gree can be difficult when errors are unknown; (3) users should be aware of non-uniqueness and the role that the model objective function plays; and (4) significant understanding of the main model features can be obtained through rudimentary exploration of model space. The issues addressed above, and in the synthetic example, underlie the importance of accurate data and the need for well-designed surveys that constrain models leaving as little room for ambiguity, and non-uniqueness, as possible. Chapter 3 Regional Scale Modeling 3.1 Introduction This chapter demonstrates the use of regional geophysical inversion in an exploration program. The goals of the regional geophysical program have been determined by the geologic and geophysical exploration models of Chapter one. This chapter presents the results of inverting regional geophysical data provided to U.B.C. by Teck Corporation, and includes a preliminary geologic interpretation based on physical property values and geologic information. The steps that are followed in this chapter are within the context of the exploration model, and are shown by shaded boxes in Figure 3.1. Extensive Tertiary and Quarternary cover makes exposure of the bedrock limited. As well as obscuring the sulphide-hosting Chilitos Formation, this has resulted in a limited understanding of the regional geology and hampers exploration for VMS deposits. Geophysical methods have been applied in order to shed light on the regional geology and search for massive sulphides underlying the overburden. Regional-scale inversion of gravity, magnetic, electromagnetic, and IP data is per formed in order to generate areas for detailed follow-up exploration. The: generation of target areas is carried out by identifying physical property distributions that fit the expected response from massive sulphide deposits and by assigning regional geologic con text to physical property features observed in the models. Both of these approaches can help focus exploration efforts. 34 Chapter 3. Regional Scale Modeling 35 Exploration Goal: •find economic mineralization Geologic Constraints: •existing mapping •historical workings •tectonic setting Geologic 'descriptive model Models: 'genetic model *H => Exploration Model Regional Geologic Mapping: •Outcrop geology •Structure •Alteration Geochemistry 'Stream sediments Remote Sensing •Multi spectral 'Hyperspectral j Geologic j -|«J Interpretation -4 No Targets Target Generation Terminate Exploratioo Program Update Geologic Exploration Model Geologic Logging A 4 I Detailed Geologic Mapping: 'Stratigraphy •Structure •Aheration Trenching Rock samples •Assays •Alteration mineralogy (X-ray/PlMA) Geochemistry •soil surveys Geologic Interpretation Mineralization Not Intersected Locate Drill-hole Terminate Exploration Program Mineralization Intersected Geophysical Exploration: Regional Scale Modeling Geophysical .expectodltf2e,(fcpUloftal8(!l Exploration .phy^ properties contrasts Model: Regional/Reconnaissance Scale Geophysics: Physical Property; Measurements Existing Data Survey Design Data Collection Gravity AMAQ (+/- ARAD) Other. VLF. MT Modeling' Individual •Forward Correlation Inverse , Joint/Simultaneous Forward Inverse ^1 of the the earth Update Geophysical Exploration Model Additional Physical Property Measurements uTK. Local/Deposit Scale Geophysics: Data Collection Survey Design |r1\ Gravity EM: CSAMT, UTEM Modeling Individual Forward Joint/Simu Itaneous Forward Inverse Correlation Conclusions about physical property distributions of the the earth Physical property core measurements Down-hole methods Delineation Goal Realized: Mineral Resource Figure 3.1: Geophysical Exploration Methodology Flowchart: Regional Scale Modeling. Chapter 3. Regional Scale Modeling 36 In the first instance, correlation methods are applied to three-dimensional models to determine regions of favorable combinations of physical property values. This involves a rudimentary spatial correlation of cells in inversion models that exhibit physical prop erty values comparable to the physical properties of massive sulphide, as determined in section 1.3. In the second instance, geologic features are interpreted from models, and the rela tion of features to mineralization is determined by using conceptualized geologic models. Interpretation of the models is performed while keeping in mind the notions of non-uniqueness and fitting the data. Ideally, the two methods will complement each other by generating fewer, but better, exploration targets than when one method is used in exclusion. This approach will hopefully result in a more efficient use of exploration resources. 3.2 Regional geology A regional map showing mapped surficial geology is presented in Figure 3.2. Apart from the geological setting that was described in Chapter 1, and the direct detection of massive sulphides, features that are of regional interest in this chapter include: (1) the varying thickness and physical properties of the overburden, (2) the location of shallow mafic intrusions, (3) the extent of a horst block, and (4) widespread fault systems and how they might relate to mineralization. Some of these features will undoubtedly be related to each other because faults define the boundaries of the horst, and overburden thickens away from the horst. The scenarios mentioned above represent a subset of the many possible geologic attributes of the area that maybe related to mineralization. However, the geologic conceptual model considered in this thesis assumes that the Chilitos formation is where we expect to discover deposits like San Nicolas, and faults play an Chapter 3. Regional Scale Modeling 37 Dlas -5000 Qal Tert Breccia Mst/Lst Mafic Vol MaficAnt Vol Qtz Rhyolite Graphitic Mst -10000 -15000 I 1 ' 1 1 1 -10000 -5000 0 5000 10000 15000 Easting (m) Figure 3.2: Plan-view of regional geologic outcrop at San Nicolas. Modified from Johnson et al., 2000. Most of the region is overlain by Quarternary alluvium. Inset shows a key to the stratigraphic units. Chapter 3. Regional Scale Modeling 38 important role in delivering mineral-laden hydrothermal fluids to areas where deposition occurs. Proximity to heat sources could also be significant to the creation of deposits. (1) Overburden: The overburden varies in thickness throughout the region, and can change even on a local scale. An increased thickness in overburden could produce a gravity variation that could mask the response of high density sulphides, while false anomalies could be caused by a thinning overburden. Away from the deposit the overburden also contains large volumes of Tertiary volcanic tuff that has high magnetic susceptibility1. Identifying the overburden, even with varying properties, plays an important part of geophysical regional exploration by allowing false anomalies, and features that could potentially mask a deposit, to be determined. In addition, the possible propagation through the overburden of reactivated, deeper faults can give clues about the underlying structure. (2) Mafic Intrusions: A small, intrusive, mafic body outcrops to the northeast of the deposit (shown on Fig ure 3.2). The mafic material has a high density and, due to the high iron content of mafic rocks, could also have high magnetic susceptibility values. These intrusives have the potential for producing physical distributions similar to massive sulphide concentra tions when using gravity and magnetic methods. Other methods are needed in order to discriminate between false anomalies and possible mineralization. Large, magnetically susceptible, intrusive bodies could also be present at depth. It is 1 Another consideration to the physical characterization of the overburden involves the interaction of water. Seasonal changes in the level of the water table may change the conductivity, and IP responses have been observed coinciding with the boundary between areas of contrasting hydraulic.permeability within the overburden. Both of these phenomena, while not being entirely explained or understood in this thesis, may have the potential for producing anomalous responses that could be mistaken for responses from mineralized bodies. Chapter 3. Regional Scale Modeling 39 expected that the deposit is proximal to a former heat source that would have driven the circulation of hydrothermal fluids. The relationship between a heat source and the deposit maybe very complicated and could depend, among other things, on the availability of fluid conduits. Nevertheless, it is hoped the identification of large intrusive bodies will make the regional picture more complete and allow geologists to make their own hypotheses about the importance of such a feature. (3) Horst: The oldest exposed rocks in the region are Upper Triassic phyllite, slate, quartzite and marble of the Zacatecas Formation ([Danielson, 2000]). The host rocks of the Chilitos formation lay unconformably on these metamorphic rocks. These rocks outcrop to the north and south of the deposit on a large horst-block of upthrust geologic section. The formation of the horst is assumed to be more recent than the deposition of the volcanic rocks of the Chilitos Formation, and more recent than the genesis of the deposit. The proximity of favorable host rocks to the surface throughout the horst makes defining the limits of the horst intrinsic to the exploration process. The gravity data show high values where these rocks outcrop and it is expected that these data will be instrumental in defining the horst. (4) Faults: Finally, the data imply the presence of several similarly orientated fault sets; faults that may have played a role in the circulation of hydrothermal fluids. If a fault can be shown to be related to the formation of the deposit, then identification of similar faults becomes an important part of exploring the surrounding area for more deposits. Regional inversions can help define these faults and determine their association with the deposit. It is hoped that airborne EM data will give insight into fault orientations nearer the surface, while Chapter 3. Regional Scale Modeling 40 gravity and magnetic methods could give information about faults at depth 2. While the existence of faults will define the extents of the horst, faults internal to the horst are also of interest. 3.3 Introduction to the inversion results. The following sections describe the results of inverting the regional geophysical data sets. Theory and methodology behind each data type is contained in Appendix C, and descriptions about the inversion programs used are in (Appendix E). The main features of each inversion model are described and, where appropriate, geologic information is inferred in an attempt to address the issues mentioned above. Information pertaining to each data set is provided in Table D.l, while Table E.l shows the inversion parameters used for each inversion. Every attempt is made to provide a clear image of the feature or features being discussed in the text although it is realized that this may not be achieved when only a few images of each model are provided. Interpretations are made from the model displayed and also from the numerous other models that were generated with different, yet just as reasonable, model objective functions and data misfit values. 3.4 Density-contrast models Quantec Geofisica de Mexico and Geociencias Consultores collected gravity data at San Nicolas in 1998. Traditional corrections, including instrument drift, tidal, and Bouguer, were applied to the data by Quantec; however, no terrain corrections were included 2The inversions of the gravity and magnetic data are designed, through implementation of the model objective function, to produce small or smooth models. Ideally a blocky model is sought/when the aim is the definition of faults. However, physical property transitions in the models may still be indicative of fault contacts and their orientation. Chapter 3. Regional Scale Modeling 41 because of the generally low topographic relief in the majority of the survey area (Fig ure D.l). Neglecting the terrain correction is probably valid except in the vicinity of a hill located at 7000m East and -1000 North on the local coordinate system. After an initial inversion it was realized that data in and around the hill were erroneous. The decision was made to discard them. The resultant data, further altered by removing their mean value, are shown in Figure D.2. Surface gravity data contain relatively strong information about lateral variations in density but they don't have any inherent depth resolution. Therefore a depth distrib ution comes from incorporating an additional depth-weighting into the inversion. The weighting procedure has been developed through the use of synthetic forward modeling and inversion [Li and Oldenburg, 1998b]. GRAV3D, a program developed by the UBC-GIF, was used for inverting the gravity data. A brief overview of this program is given in Appendix E.2. For the inversion, each datum was supplied with a location and assigned a standard deviation of 0.01 mGal. Topographic relief was included in the inversion pro cedure. Table E.l provides more details about the inversion parameters, as it does for all of the inversion results presented in this thesis. The minimum data misfit that was achieved, without producing too much complexity in the model, was about twenty five times the number of data. This number seems larger than expected even if the data errors were assigned inappropriately3. The reason for this may reside with the discretization of the model. In order to model such a large volume, within a reasonable amount of time, a cell size of 250m by 250m by 100m was used for the majority of the model. The data were collected on intervals ranging from 100m to 50m. As a result, shallow density variations that may have been detected by the data, would 3Errors of 0.01 mGal magnitude correspond to the normal operating sensitivity of most gravime-ters. Applying gravity corrections (Instrument, Bouguer, etc.) to the data does not usually introduce significant errors that would degrade the accuracy of the data. Chapter 3. Regional Scale Modeling 42 not be modeled by a cell of uniform density-contrast and the model will not be able to fit the data. However, it is expected that large, regional-scale features are acceptably modeled. 3.4.1 Regional gravity inversion results A plan-view of the regional density-contrast model is shown in Figure 3.3. All the cells that have a density-contrast of 0.1 g/cm3 or lower have been removed. Outcrop geology is overlaid and a boundary shows the extent of data coverage. A cut-off of 0.1 g/cm3 is somewhat arbritrary because the transition from areas of low density-contrast to areas of high density-contrast is smooth. However, this does show the distribution of the denser materials, and the western edge of the dominant north-south feature moves only slightly when the cut-off value is varied. Outcropping rocks that correspond to the Chilitos and Zacatecas Formations roughly coincide with the surface projection of the northern half of a large dense north-south body, which approximately extends from -5,000m to 2,000m east, and -10,000m to 10,000m north. The continuation of this high density feature to the south, together with knowledge of the regional geology, suggests the rocks persist under the Tertiary and Quarternary overburden to about -10,000m north. This large feature probably represents uplifted geologic section and is interpreted as a horst. A north-south trending string of high density anomalies seems to mark the western edge of the horst. For the southern half of the horst, the east and west edges show a north-northeast orientation. This orientation is consistent with a mapped fault located between -750m east, 0m north, and 1,500m east, 7,000m north, which is also shown in Figure 3.2). Four northwest-southeast linear anomalies are also present east of the main horst feature. Westerly projections of these linear features coincide with breaks and inflections along the east and west edges of the horst, suggesting a relationship between the features. The relationship described above can be explained with northwest-southeast faults that Chapter 3. Regional Scale Modeling 43 15000 10000 5000 £ o sz o -5000 -10000 -15000 -10000 t *> fc«* s IfogJi a « t \ \l ^ HI |K i —*\ /•- •r' • \ r\ » y \ - ^ / -5000 0 5000 Easting (m) 10000 15000 0.05 0.175 g/cm3 0.3 Figure 3.3: Plan view of the regional density contrast model with cells below a value of 0.1 g/cm3 invisible. Outcrop geology is overlaid, and the boundary of the gravity survey is represented by a purple line. All coordinates are in the local grid system. Chapter 3. Regional Scale Modeling 11 Figure 3.4: West-facing cross section of the regional density contrast model with cells below a value of 0.1 g/cc invisible. This cross-section is located at -1600m east and shows an interpreted horst. laterally off-set sections of the horst. The linear features themselves, could be faults containing higher density materials. A west-facing cross-section located at -1600m east shows a side view of the horst (Figure 3.4). The image shows an east projection of material between -1,600m and -3,000m east that has a density-contrast of greater than 0.1 g/cm3. The image shows that the horst, as defined by the density cut-off, is not continuous. In fact three sections of the horst could be interpreted, with the two most southerly sections being joined at depth and the northern section being larger and deeper. This agrees with the idea of northwest-southeast faults being present, as the faults could have been responsible for changing the elevation of the rocks as well as offsetting them laterally. Figure 3.5 shows the inferred locations of faults in the area based on evaluation of this model. The straight trace of the northwest-southeast faults, and the cross-cutting of the east and west flanks of the horst, imply an age relationship between faults in the Chapter 3. Regional Scale Modeling 45 15000 -10000 -15000 -10000 -5000 5000 10000 15000 Easting (m) 0.05 0.175 g/cm3 0.3 Figure 3.5: Surface projection of the regional density-contrast model with inter preted faults overlaid. Northwest-southeast trending faults are shown in red, while north-northeast trending faults are shown in black. The mapped fault is shown in yellow. Material with a density-contrast of less than 0.1 g/cm3 has been removed. Chapter 3. Regional Scale Modeling 46 area, with the north-northeast fault being older than the northwest-southeast faults. Figure 3.6 shows a horizontal slice through the regional density contrast model at an elevation of 1700m (400m below the surface). The main features are the linear trends of high density that correspond with the east and west edges of the horst and northwest-southeast trending linear features. It is interesting that high densities would be associated with faults. An explanation for this could be the presence of higher density material within faults. This could be representative of igneous rocks that were intruded along the fault. An anomaly located at -1,600m east and -400m north represents the deposit. This is a small anomaly and selecting this feature over the others would be difficult. Based on the interpreted faults that are shown in Figure 3.5, the deposit lies within a coherent block of bedrock. This suggests that any faults associated with the genesis of the deposit are not detected and are probably older than the formation of the horst. Searching for density anomalies that are not coincident with younger, interpreted faults seems like a reasonable way to search for massive sulphides. However, other interpretations might place faults directly through the deposit. A north-facing cross-section through the regional density-contrast model is shown in Figure 3.7. The negative density-contrast located near the surface at 6,000m east is in the vicinity of the hill where data were removed prior to inverting, and is thought to be an artifact created by not terrain correcting the data around the hill. The deposit is shown as a discrete density high near the surface. No significant depth extent to the deposit has been modeled. Another isolated density high is shown at about 4,000m east, this is one of the northwest-southeast linear features. A lateral change in density contrast at -6,000m east represents a fault, which is probably representative of the near vertical western flank of the horst block. The eastern flank of the horst block maybe dipping to the east. Chapter 3. Regional Scale Modeling 47 Figure 3.6: Horizontal slice through the regional density-contrast model at an elevation of 1700m (400m below the surface at the location of San Nicolas). An outline of geology outcrops is overlaid, and the extent of data coverage is shown by a purple line. Chapter 3. Regional Scale Modeling 48 5000 Q. Q 10000 -10000 -0.4 San Nicolas -5000 0 5000 Easting [m] 0.05 10000 15000 g/cm3 0.3 Figure 3.7: North facing cross-section through a regional density contrast model at a northing of -400m (coincident with the deposit). Chapter 3. Regional Scale Modeling 49 At depths greater than 4 kilometers, a large amount of higher density material is present in the model, and it is coincident with the bounds of the horst block as defined from the shallower parts of the model. The deep large high density features, needed to explain the high gravity anomalies shown in the data, could be uplifted limestone, mudstone, or metamorphic rocks of the Zacatecas Formation. All of these rocks could have higher densities than the overlying volcanic and sedimentary units. In addition, the deep-rooted, large bodies at depth could represent large intrusive bodies rich in mafic material and high in density. 3.4.2 Summary of regional gravity results The density contrast model is consistent with the previously interpreted geology and has contributed to regional geological understanding with the following interpretations: (1) the delineation of a horst block, (2) the interpretation of faults below the overburden that cross-cut the horst (3) a possible discrimination of discrete high density anomalies similar in nature to the deposit, and (4) the identification of deep-rooted bodies that could represent large amounts of sedimentary or intrusive igneous rocks, or both. 3.5 Magnetic susceptibility models This section addresses the the inversion of regional magnetic data. Total-field magnetic data is inverted with the aim of producing a three-dimensional magnetic susceptibility distribution of the earth. Background information on the mag netic method is given in Appendix C.2. Airborne methods are employed at San Nicolas to provide generalized geologic information on large areas in a short period of time as well as identify anomalous regions of high susceptibility that might indicate mineralization. The data are displayed in Appendix D. Three data sets are inverted and subsequently Chapter 3. Regional Scale Modeling 50 combined to form a single susceptibility distribution model of the subsurface. Surface magnetic data also have no inherent depth resolution and so, as with the gravity inversion, depth weighting is needed. 3.5.1 Data The three airborne magnetic data sets include: 1. Dighem collected airborne magnetic data while also collecting frequency-domain electromagnetic and radiometric data. The Dighem data has the most extensive coverage of the three data sets. Northwest-southwest trending features that can be interpreted from the regional density-contrast models are also apparent in Fig ure D.3. This texture betrays the regional structural setting. A subset of observa tions from the Dighem data set was chosen for inversion (Figure D.4). 2. The airborne magnetic data, shown in Figure D.5 were acquired by BHP at the same time as airborne gravity-gradient data were being recorded with their Falcon system. 3. Fugro collected airborne magnetic data (Figure D.6) while collecting time-domain electromagnetic data with their Questem system. The coverage and survey specifications for each data set are slightly different although similar features are observed. The deposit is represented in the data by a subtle magnetic high at about -1600m east and -400m north, and the response from highly magnetic Tertiary tuffs are seen to the north west of the deposit. All of the data sets had regional fields removed prior to inversion and standard deviations of lOnT were assigned to each datum. Chapter 3. Regional Scale Modeling (a) Dighem susceptibility model. (b) BHP susceptibility model. 5000 -5000 0 Easting (m) 5000 -5000 m 1 -5000 0 Easting (m) xirrs.1. xKTS.I. (c) Fugro susceptibility model. (d) Composite susceptibil ity model. Figure 3.8: Comparison of individual magnetic susceptibility inversion models with com posite model of mean values. The models are shown in plan-view looking at a depth of 300m (within the depth of the deposit). Chapter 3. Regional Scale Modeling 52 3.5.2 Composite susceptibility model Initially an attempted was made to invert all three data sets simultaneously. The results from this inversion were not successful; this is possibly due to the different regional fields that were removed prior to inverting. Consequently each data set was inverted separately. A single susceptibility model was produced from the three individual inversion models by placing them all on the same mesh and calculating mean value of each cell. Volumes of the model that do not lie directly below the data coverage have been removed. Figure 3.5.1 shows a comparison between the three individual models and the resulting model. 3.5.3 Regional susceptibility model analysis Figure 3.9(a) shows a plan view of the regional magnetic susceptibility model at a depth of 300m4. At this depth the deposit is modeled as a discrete anomaly located at about -1700m east and -500m south (labeled). However, features with the highest susceptibility include two elongated regions to the west and northwest of the deposit (between -4000m and -2000m east, and -500 and 750m north), and a discrete anomaly located at -1500m east and 2000m north. The two elongate features correspond with out-cropping Tertiary breccia. The other anomaly to the north (located at -1500m east and 2000m north) corresponds with out-cropping mafic volcanics. Except for a small anomaly located at -2000m east and -2000m north, other mafic volcanic outcrops have low susceptibilities. Several other low magnitude anomalies of unknown origin are also observed at this depth. At a depth of 1200m, as shown in Figure 3.9(b), four large features are observed. An anomaly at -1500m east and 2000m north is probably an extension of the anomaly seen at 300m depth, and could represent an intrusive body. A large anomaly is observed at this depth, located between -4000 and -2000m, and -1000m and -3000m. It is interesting 4The deposit is approximately located between depths of 200m and 500m. Chapter 3. Regional Scale Modeling 53 -10000 -5000 _ „ , , 0 Easting (m) 5000 xlO"3S.I. 8 (a) Composite regional susceptibility model: 300m below the surface. 5000 cn c 'sz tl o -5000 -10000 1 1 pi'' i 1 P -5000 Easting (m) 0 5000 xlO"3S.I. (b) Composite regional susceptibility model: 1200m below the surface. Figure 3.9: Plan-view of horizontal slices through the composite magnetic susceptibility model. The upper panel shows a slice at 300m depth, while the lower panel shows a slice through the model at 1200m depth. Chapter 3. Regional Scale Modeling 54 to note that a shallow ring of small, moderate ( 2xl0~3 S.I.) susceptibility concentra tions coincides with the edge of this body (Figure 3.9(a)). This volume of anomalous susceptibility could be a deep rooted intrusion, and the shallower anomalies may have originated from this deep source. Alternatively, this susceptible body might represent magnetic units below the Chilitos Formation. The anomaly located between -4000 and -2000m, and -1000m and -3000m also has a distinct northern edge that is orientated in a northwest-southeast direction. This orientation corresponds to those interpreted from the regional density-contrast model. This suggests that it was in place prior to a faulting event or its emplacement was controlled by an existing fault. The western edge of this susceptible body is also distinct but is likely an artifact of combining the three separate susceptibility models. Anomalous susceptibility distributions have also been modeled at -6000m east, 0m north, and -2500m east, 3000m north. Both of these are constructed to fit incomplete anomalies at the edge of the data coverage. They could also represent deep intrusive bodies or magnetic sequences within the Zacatecas Formation. Some of the features mentioned above are observed in the north- and west-facing cross-sections of Figure 3.10. The north-facing cross-section is located at -400m south and corresponds to the location of the deposit. The deposit is clearly represented by a discrete susceptibility anomaly (labeled). The anomaly near -4000m east is the southern boundary of the Tertiary tuffaceous unit. The large anomaly in the far west of the cross-section, between depths of 100m and 2500m, corresponds to the most westerly anomaly seen in figure 3.9(b). This anomaly does not extend to great depths although it might extend beyond the modeled volume. The west-facing cross section shows the deposit which is connected to a smaller mag netic feature to the north. The high susceptibility concentration at -1500m east and 2000m north is clearly seen on this cross-section, along with the large deep body to the Chapter 3. Regional Scale Modeling 55 San Nicolas Q-"a 5000 -10000 5000 Easting (m) xl03S.I. N Northing (m) 5000 Figure 3.10: North- and west-facing cross-sections (upper and lower sections respectively) through the magnetic susceptibility model generated by combining inversion results from Dighem, Fugro, and BHP airborne magnetic data. The north-facing cross-section is positioned at -400m south, while the west-facing cross-section is positioned at -1600m east. Chapter 3. Regional Scale Modeling 56 south. This large, most southerly anomaly is the only anomaly that persists below a depth of 2500m and continues to the bottom of the mesh at a depth greater than 5000m. 3.5.4 Geologic interpretation Mafic rocks make up part of the volcanic section within which the deposit was formed. The anomalies observed in the mafic volcanics could represent concentrations of mag netite or pyrrhotite associated with massive sulphides, as seen at San Nicolas. Alterna tively, the mafic rocks may naturally exhibit a higher magnetite content5. The Tertiary volcanic breccia, while predominately being comprised of subordinate volcanic debris, contain some ash-flow tuffs that are magnetically susceptible6. These tuffs are responsible for the elongate anomalies seen in figure 3.9(a). The extention of the susceptibility values to the north of the deposit may be indicative of an extension of mineralization and warrants further investigation. The large susceptibility features seen at depth could indicate large intrusive bodies. This is especially true when the relationship between some of the deeper and shallower anomalies are considered. It seems unlikely that thick sections of susceptible, unaltered volcanic rocks reside at depth. Evidence for the alteration or metamorphism of rocks on a regional scale comes from the rocks of the Zacatecas Formation, which underly most of the region. Some linear features can be seen in the model, such as an approximate north-northwest line of susceptibility anomalies, which includes the deposit, (seen in figure 3.9(a)), and well defined boundaries to some features that share the same orientation. These linear 5Volcanic rocks in the San Nicolas region generally exhibit subdued magnetic susceptibility values, possibly due to the effects of wide spread alteration. Less altered volcanic rocks could posses higher magnetic susceptibility values. 6Drill hole Sal-133 intersected Tertiary breccia to a depth of more than 400m [Jansen, 2000]. Sus ceptibility measurement made on the core exceed 2 xl0~3 S.I. between a depth of 200m and 300m (Appendix B.3.1). Chapter 3. Regional Scale Modeling 57 features might represent faults although collaborative information would be necessary to confirm this. 3.6 Regional conductivity models Airborne electromagnetic data are used to determine the regional conductivity structure of the subsurface and are usually acquired contemporaneously with airborne magnetic data. The development of airborne electromagnetic systems has been hand-in-hand with the exploration of massive sulphides as the targets are usually very conductive and can be hosted in quite resistive rocks. Dighem's frequency-domain system and Fugro's Questem time-domain system were employed at San Nicolas and modeling of both these data sets was performed. 3.6.1 Frequency-domain data An overview of the frequency-domain, airborne electromagnetic method is given in Ap pendix C.3.1. The overview describes the nature of the data and Appendix E.5 gives a description of EM1DFM, the algorithm that is used to invert frequency-domain EM data to recover a one-dimensional conductivity distribution at each station. The data cover a very large area (Appendix D) and it was hoped that a regional conductivity model of the earth could be produced to complement the 3D density contrast and magnetic susceptibility information that was obtained by inverting the gravity and magnetic data. Unfortunately the inversions have shown that the data are not sensitive to conductive features below the overburden. This lack of sensitivity was determined by performing several inversions with different reference models and observing the depth at which the model reverts to the value of the reference model. Both GCV and L-curve methods of choosing the appropriate misfit were used and generate similar outcomes. This suggests Chapter 3. Regional Scale Modeling 58 that inappropriate error assignments do not contribute to the shallow penetration. The depth at which the data were deemed insensitive to the model was about 150m. The reason for the observed sensitivity of this system is probably because of the high overburden conductivity, the frequencies employed, and the geometry of the system. Insight into the observed depth of penetration can be explained by considering the skin depth for a plane wave in a homogenous medium [Ward and Hohmann, 1987]7. For a 30 Ohm-m medium at the lowest frequency of the system (466 Hz), the skin depth is approximately 127m. This implies that the signal has dramatically deteriorated by the time it passes through the conductive over burden and reaches the volcanic sequence. At higher frequencies the signal will be even more deteriorated by passing through the overburden. Aside from the lack of penetration, the survey is useful in that it covers a large area and provides detailed information about the thickness of the overburden. Additionally features can be interpreted in the overburden that might be the manifestation of deeper structures. After the ID inversion models were produced an interpolation routine was used to construct a 3D model8. The 3D models are presented below in the same manner that the 3D density and magnetic susceptibility models are presented, although it should be kept in mind that they have been produced by stitching together many ID conductivity models. 7The skin depth (<5) is the depth at which a sinusoidal electric field (from a planar source in a whole-space) diminishes to about 37% of it's original strength, and is approximately stated as 5 = 503 [m], where p is the resistivity, and / is the frequency [Ward and Hohmann, 1987]. In the Dighem system we are not dealing with plane waves but a dipolar source. Also the geometry between source and receiver influences depth of penetration. However, this approximation can give an idea of the strength of different frequency signals in different mediums. 8The interpolation routine uses triangulation to assign values to cells at each depth. The routine was applied to the Logio of the conductivity values in order to reduce homogeneity artifacts. Chapter 3. Regional Scale Modeling 59 15000 10000 5000 .£ 0 o z -5000 -10000 -15000 -10000 \ o <*•#•:» 0 TV O \ * « -TO) '"' 1 s il ) -5000 0 5000 Easting (m) 10000 15000 100 Ohm-m 2000 Figure 3.11: Horizontal slice 60m below the surface of the 3D resistivity model generated by interpolating between lines of ID models inverted from Dighem frequency-domain airborne electromagnetic data. Chapter 3. Regional Scale Modeling 60 A horizontal slice through the stitched conductivity9 model is shown in Figure 3.11. The slice is at a depth of 60m below the surface and it clearly maps out the outcropping geology with high resistivity values, and the overburden with low resistivity values. Out cropping volcanic units of the Tertiary overburden exhibit low resistivity values. This is especially significant in the south, where a tuffaceous unit outcrops. No anomaly is detected in the vicinity of the San Nicolas deposit. Although this model doesn't seem that useful in the direct search for massive sul phides, the method can be used as a reconnaissance mapping tool for determining the extent of cover and inferring regional structures. Linear features seen in the model have the same orientation as structures observed in the regional density contrast model. Specif ically, three discrete resistive regions are defined that correspond well with the limits of the horst (as defined by the density contrast model). Northwest-southeast trending linear resistivity contrasts are seen that may also relate to the structural nature of the region. 3.6.2 Time-domain data Airborne time-domain electromagnetic (TEM) surveys are used at the regional or recon naissance stage of exploration programs. They often enjoy more success than frequency-domain systems at resolving the conductivity structure of the earth below conductive overburdens, and can be used to detect conductors as well as delineate structures. This section details the inversion of such a data set collected by Fugro in 1998. A cursory investigation is also made into the possible cause of negative data that is observed at late time-channels. Time-domain electromagnetic data were collected using Fugro's Questem system. The nature of the data and how they are collected is discussed in Appendix C.3.2, while 9The conductivity (a) models are presented as resistivity (where p = 1/<T) for consistency with models presented later in the thesis. Chapter 3. Regional Scale Modeling 61 the specifications particular to this survey are given in Table D.l. The x-, y-, and z-components of the rate of change of the magnetic field (^) are measured as a function of time, but only the vertical component is considered here. One stacked time-decay series (dB^ measured in fifteen time-windows) is collected for each station. Each time-decay series was inverted using TEM1D to recover a layered-earth conductivity model. The resulting ID models were combined, using an interpolation routine, to produce a 3D conductivity structure of the earth. An overview of TEM1D is given in Appendix E.6. A horizontal slice through the "stitched" conductivity model at a depth of 200m is presented in Figure 3.12. The 3D conductivity model shows the distribution of conductive overburden and maps out the outcrops quite effectively, similar to the frequency-domain results presented in Section 3.6.1. However, not much structure is modeled below the overburden, and the absence of a conductive body in the vicinity of the deposit is con spicuous. It was expected that data collected with this system would be sensitive to conductivity structure below the overburden10, and would probably result in a more deeply resolved inversion model than that recovered from the inversion of the Dighem frequency-domain data. The reason for this inconsistency might lie in how standard deviations were assigned to the data to account for the "noise". 10The depth of the peak amplitude of a decaying and propagating, impulsive, electromagnetic plane-wave, in a whole-space, at a given time, can be stated as: z = f ^) [m] [Ward and Hohmann, 1987]. This is sometimes referred to as the diffusion depth. Although this depth would be different for a layered earth, and a magnetic dipole source, as used in the Questem survey, this gives a general idea of the penetration of time-domain electromagnetic systems. For an overburden with a conductivity of 20 Ohm-m, the diffusion depth is 280m for late times, which is well below the top of the deposit. Chapter 3. Regional Scale Modeling 62 Figure 3.12: Plan view of a slice through the 3-D resistivity model generated by interpo lating between lines of ID models from the inversion of Questem time-domain airborne electromagnetic data. The slice (at a depth of 200m below the surface) shows, in red, the thicker areas of overburden that continue to this depth. The underlying, more resistive, volcanic rocks are represented by cooler colors. Chapter 3. Regional Scale Modeling 63 Negative data A significant amount of negative data in late time-channels, initially thought to be noise, was present throughout the survey. Large errors were assigned to these data. Assigning these errors had the effect of making the inverted models insensitive to the negative data, but also to structures at depth. As a result, the resulting model only detects the shallow i feature, namely the overburden. Sixty percent of the stations contained negative data, usually at late time-windows. However, sixteen percent of the stations contained negative data from time-windows that correspond to earlier times (time-channels 8 to 11); channels that usually contain little noise. Furthermore, these negative values persisted locally from station to station and appear to be above the noise levels, particularly over the deposit. A line of Questem, z-component data is presented in Figure 3.13, along with a total magnetic field profile, and the results of inverting the data using TEM1D. A significant negative z-component response is observed over the deposit in time-channels 10 to 15 (annotated with a number 1). Smaller depressions that are possibly above the noise level, are also observed at late times (labeled with 2, 3, and 3). Line-by-line inspection of the whole data set showed that no other late-time response was as compelling as that seen over the deposit. An example of time-decay data that exhibit negative values for one station over the deposit is provided in Table 3.1 and plotted in Figure 3.14. Although this time-decay series is from over the deposit, it's characteristics are typical of other time-decay series in that it exhibits a faster decay at late time and then changes sign. The negative data are part of a coherent decay pattern observed in the time-decay series. It is strongly suspected that the negative data are a real part of the signal and are not noise. Such negative data are commonly referred to as IP effects. Negative data occur when the rate of change of the magnetic field, as measure by a receiver coil, opposes Chapter 3. Regional Scale Modeling til D Channel 10 Channel 11 Channel 12 Channel 13 -Channel 16 • • 43830 43840 law fieis MaQnafict 43820 !nl 43810 43800 43790 E. 500--2000 -3000 Resistivity Model Approximate location of San Nicolas deposit Ohm-m -1250 Easting [m] 2300 Northing |m] Figure 3.13: Vertical component, time-channel data, flight height, total field magnetics, and resistivity inversion model for line 20110 of the Questem data. The data are pre sented as parts-per-million of the peak z-component of the waveform, with channels 1 and 15 corresponding to the earliest and latest times respectively. The inversion model is produced by "stitching" together ID models. The location of the deposit corresponds to a coherent pattern of negative data values in time-channels 12 through 15 (1), and a magnetic high. Other regions of negative data that could be above the noise level are also observed (2, 3 and 4). Chapter 3. Regional Scale Modeling 65 Time Window center time z-component z-component Channel (msec) (MV) (ppm) 1 0.375 0.080425279 15834.6 2 0.531 0.038734626 7626.3 3 0.688 0.022706049 4470.5 4 0.922 0.012765264 2513.3 5 1.234 0.006915682 1361.6 6 1.625 0.003740238 736.4 7 2.094 0.002087504 411.0 8 2.641 0.001168190 230.0 9 3.344 0.000640473 126.1 10 4.281 0.000316427 62.3 11 5.531 0.000087868 17.3 12 7.094 -0.000008127 -1.6 13 9.047 -0.000083297 -16.4 14 11.469 -0.000070091 -13.8 15 14.359 -0.000042156 -8.3 Table 3.1: Questem time-decay (^) data for one station located over the deposit. The center time of each time channel is given (time after turn-off of primary field) along with the z-component data in micro-volts and in ppm, as normalized by the peak z-component of the primary field. Chapter 3. Regional Scale Modeling 66 1.E-06 -I I I I I II II | I I I I I I I I | I I I I Mil) 0.1 1 10 100 Time [msec] a Observed data: positive A Observed data: negative Figure 3.14: Logarithmic plot of Questem vertical-component time-decay data given in Table 3.1. Positive data are shown by squares, and negative data (the last four) are shown with triangles. The data are from over the deposit and show a typical decay pattern for stations that contains negative data. Chapter 3. Regional Scale Modeling 67 that of the primary field. This could be caused by geometric effects distorting the field, or by the presence of polarizable material that reverses the field. In addition, for ID conductivity structures with realistic ranges of magnetic susceptibility, negative data can only be produced when polarizable is present material [Weidelt, 1982]. Flis, Newman, and Hohmann [Flis et al., 1989] have also shown IP effects in TEM data to be caused by polarizable, or chargeable, material at the surface. As the massive sulphides at San Nicolas are expected to be chargeable, further investigation is performed on the negative data, through the use of ID forward modeling. 3.6.3 Forward modeling the Questem data Forward modeling is performed in an attempt to reproduce negative data that are ob served at late times in the Questem data set. Two investigations are undertaken; one to determine if the geometry of the earth can produce the negative data, and one to determine the role of polarizable material in producing these data. Forward modeling a 3D conductivity structure The affect of geometry is investigated by forward modeling synthetic Questem data over a three-dimensional conductivity structure. The conductivity model is generated by dig itizing geologic cross-sections and assigning conductivity values to different rock types based on physical property measurements (Appendix B.4). The conductivity sections are interpolated to produce a three-dimensional model (Figure B.5 shows a north-facing cross-section through the conductivity model at -400m south). EH3DTD (described in Appendix E.9) was used to collect synthetic ^f- data at twenty four stations, spaced 100m apart along one east-west line, over the conductivity structure. The same sur vey parameters (transmitter-receiver geometry, wave-form, time-channels) were used to collect the synthetic data as those employed by the Questem system. Chapter 3. Regional Scale Modeling 68 All of the synthetic data were positive. Moreover, the synthetic time-decay data over the deposit exhibit a slower decay at late times, which is in contrast to that observed in the field data. The decay-rate is also slower than that observed over a half-space. This result suggests that the negative data were not due to the geometry of the subsurface. ID forward modeling of polarizable material The role polarizable material might play in producing negative data is investigated by forward modeling one-dimensional, layered-earth structures that contain polarizable layers11. A ID forward modeling program was available that uses the Cole-Cole [Pelton et al., 1978] model to represent complex conductivities12. Models were produced in an attempt to match time-decay data observed at one station over the deposit. The time-decay cho sen for forward modeling analysis was recorded with the receiver located at -1600N and -500N and is representative of many of the decay curves over the deposit. This data is presented in Table 3.1 Figure 3.14. The observed data show a reversal in polarity at the twelfth time-channel (about 7 msec). Figure 3.15 shows observed and predicted time-decay data and the corresponding real 11In this case, polarization refers to a relaxation mechanism at the pore-scale that occurs when al ternating, or time-dependent currents are present. Polarizable material has frequency-dependent, or complex, conductivity. 12 Complex conductivity is introduced in the model by using the empirical Cole-Cole model. This is a simple relaxation model that allows conductivity to change with frequency based on a few parameters. An equation for the Cole-Cole model is given below: where Z(u>) is the frequency-dependent impedance, to is the frequency, Ro is the zero-frequency resis tivity, c controls the rate at which the relaxation occurs, r is a time constant that determines when the relaxation occurs, and m is the chargeability, and relates the zero-frequency resistivity to the resistivity at infinite frequency. Chapter 3. Regional Scale Modeling 69 1200 (a) Observed and predicted time-decay data. (b) TEM1D inverted conductivity model. The predicted data are from the TEM1D inver sion. (c) Observed and predicted time-decay data for (d) Conductivity model with additional conduc-conductivity model with additional conductor at tor,at depth. depth. Figure 3.15: Observed and predicted time-decay data, and corresponding ID conductivity models for Questem station located over the deposit. Frequency-independent conductiv ity layers do not reproduce the observed negative data. Chapter 3. Regional Scale Modeling 70 (frequency independent) one-dimensional conductivity models. Figure 3.15(a) shows the observed field data and predicted data from the inversion using TEM1D. The predicted data fit the early time data very well, negative signals are not produced at late times. As a result the overburden is well represented in the model, but most depth information is not realized. Based on the expected conductivity structure of the deposit from physical property measurements, a conductive layer of 200 Ohm-m was introduced into the model between 200m and 400m depth. This is shown in Figure 3.15(d), and the data resulting from this model is shown in Figure 3.15(c). The introduction of a conductor into the model increases the magnitude of the data, and a worse fit to the observations is produced. When the conductive layer at depth is assigned a complex (frequency-dependent) conductivity, the data values are reduced, negative data are produced, and the character of the observed time-decay is emulated. The results from two complex conductivity models are shown in Figure 3.16. Model "A" (Figure 3.16(b)) consists of the one-dimensional real conductivity distrib ution that was recovered from the inversion with a layer of complex conductivity between depths of 200m and 400m. The conductivity of the overburden was also decreased to 100 Ohm-m in order to more closely fit the early-time data. This model produces data (Figure 3.16(a)) that emulate the observed time-decay data. The predicted late-time data are negative and have the same decay pattern as the observed data. The early-time data are also a good fit. Model "B" (Figure 3.16(d)) consists of the one-dimensional real conductivity distrib ution that was recovered from the inversion with a layer of complex conductivity between depths of 350m and 480m. This model produces predicted data (Figure 3.16(c)) that also emulate the observed time-decay data. The predicted late-time data are negative and have the same decay pattern as the observed data. The late time data also have a better Chapter 3. Regional Scale Modeling 71 1 Time [msec] 10 • Observed data: positive Observed data: negative • Predicted data: positive s Predicted data: negative (a) Observed field data and predicted time-decay data for complex conductivity model A. 10000 1000 ? E 100 O •> "35 10 V c=0.4 m=0.7 T=0.005 200 400 600 Depth [m] 800 1000 1200 (b) Complex conductivity model A. A layer of complex conductivity (Cole-Cole parameters are shown) is present between depths of 200m and 400m. 1 Time [msec] 10 Observed data: positive ; Observed data: negative |4 Predicted data: positive « Predicted data: negative (c) Observed field data and predicted time-decay data for complex conductivity model B. 10000 1000 E £ 100 .c O * 10 > "95 c=0.785 m=0.8 r=0.005 200 400 600 Depth [m] 800 1000 1200 (d) Complex conductivity model B. A layer of complex conductivity (Cole-Cole parameters are shown) is present between depths of 350m and 480m. Figure 3.16: Observed and predicted data, and ID conductivity models for Questem station located over the deposit. A layer with complex (frequency-dependent) conduc tivity (based on the Cole-Cole dispersion model), positioned at the approximate depth of the deposit, can produce data that has the same character as the observations. The conductivity of the overburden was also decreased in an attempt to more closely fit the early-time data. Chapter 3. Regional Scale Modeling 72 fit to the observations than those in Figure 3.16(a). However, the early-time data do not fit the observed data so well. The Cole-Cole parameters were manipulated to produce data that resemble the ob servations. Although the data have not been fit exactly, the results seem acceptable for only having introduced one conductive layer. The Cole-Cole parameters used were Ro = 200 Ohm-m, c = 0.4, r = 5 msec, and m = 0.7, for model "A", and R0 = 100 Ohm-m, c = 0.785, r = 5 msec, and m = 0.8, for model "B". The parameters correspond quite well with numbers determined by Pelton [Pelton et al., 1978] for massive sulphides, although values of c less than 0.55 are expected. It is not unreasonable to expect that, by slightly changing the conductivities and Cole-Cole parameters at different depths, a model can be found that fits the data within the expected noise levels and has chargeable material at depth. Discussion on inverting complex conductivity data Since a forward modeling is available that produces time-decay data from complex con ductivities, it might seem reasonable to use this to invert the data. However, there are issues that should be considered before this is undertaken. First the appropriateness of the Cole-Cole model should be demonstrated. Pelton [Pelton et al, 1978] has shown that many sulphide minerals cannot be described by one set of Cole-Cole parameters alone. The use of two or more Cole-Cole models is necessary to fit laboratory observations. Secondly, inverting for more unknowns could introduce more non-uniqueness into the problem. That is, each layer of complex conductivity will have its own frequency dependence, and will be characterized by four parameters, a, m, c, and r. In contrast, each frequency-independent layer needs only one parameter (a) for full characterization. This means that the number of unknowns will increase by a factor of four if only one Cole-Cole model is employed. Chapter 3. Regional Scale Modeling 73 Alternatively, so as not to be restricted by any parameterized frequency-dependent models, real and imaginary parts of the complex conductivity could be inverted for at each frequency. Routh (1999) investigated this for multiple grounded sources in order to remove coupling effects form IP data. However, data from one source may not contain enough information to constrain both depth and frequency-dependent information. Negative data are in the noisiest part of the time-decay (late times). These data may not be sufficient to resolve between different combinations of the Cole-Cole parameters. In such cases, additional information would be needed to constrain or eliminate some of the parameters from the inversion. Better constraints could come from the data collected with TEM systems with increased sensitivity and better sampling of the decay series. Alternatively, rock samples in the area may reveal only small variations in some of the Cole-Cole parameters, thus allowing them to be constrained or removed from the inversion altogether. The discussion so far has only dealt with one-dimensional complex conductivity struc tures. It is possible that no combination of complex conductivity layers will fit the data. The combined effect of three-dimensional geometries and complex conductivities is un known, but may be necessary in order predict some decay-series. 3.6.4 Summary of regional conductivity methods. Both frequency-domain data and early-time, time-domain data are able the model the overburden quite well, and regional structures that manifest themselves in the overburden structure can be identified. While it is not expected that the frequency-domain data can provide much more information than that which has already been modeled, the time-domain data has significant potential for providing more information about the subsurface. Forward modeling has shown that the decay response produced over the deposit could be explained with a buried, polarizable conductor. While other models Chapter 3. Regional Scale Modeling 74 might also be able to explain the observed data, and cannot be ruled out, the concept of detecting massive sulphides with this data becomes conceivable. A benefit of accounting for polarizable material would be in modeling non-polarizable, conductive bodies more accurately because this would effectively remove IP effects that come from geologic noise. Until an inversion routine that models complex conductivity is available, it is likely the majority of the data can only be inverted to define the overburden, and complex conductivity forward modeling analysis could be applied to select parts of the data. 3.7 Chargeability models At San Nicolas, gradient-array and "Real-section" array configuration IP data were ac quired by Quantec. An overview of how data are generated using the IP method is given in Appendix C.6. The gradient-array data is used for regional modeling as it covers a rel atively large area (Figure D.10). The large chargeability anomaly, indicated by a strong positive anomaly at -1600m east, -400m north, is associated with the deposit. There is no doubt that this is valuable information regarding the existence of an ore-body. In fact, this was the data that prompted the discovery of San Nicolas. In the gradient array configuration, the current electrodes are located outside of a rectangular area to be sur veyed and the IP potential data are collected using a roving-receiver dipole within the rectangular area. The data are generally used only as a mapping tool because detailed depth information cannot be obtained without multiple transmitter locations. However, the similarity of gradient array IP data to magnetic field data allows for the data to be inverted using MAG3D through which depth information is introduced with sensitivity weighting. A description of how gradient IP data can be inverted using MAG3D is given in Appendix F, along with an example. Inversion parameters such as noise estimates and cell sizes are provided in Table E.l. Chapter 3. Regional Scale Modeling 75 Figure 3.17: Plan-view of a slice at 300m below the surface through the chargeability model generated by inverting gradient-array IP data as magnetic data using MAG3D. The deposit is the largest chargeable body detected by the survey. The extents of the survey are shown by the red outline and outcrop geology is shown in black. Chapter 3. Regional Scale Modeling 76 Q. OJ TD 5000 -10000 5000 Easting (m) 0 50001 L_ _J -5000 0 5000 Northing (m) Figure 3.18: North- and west-facing cross-sections (upper and lower sections respectively) through the chargeability model. The chargeable anomaly is well defined yet extends below the elevation of the deposit, as defined by drilling. Chapter 3. Regional Scale Modeling 77 Figures 3.17 and 3.18 show the modeled distribution of chargeable material on slices through the deposit in plan and cross-sectional views. As expected from the data, the anomaly corresponding to the deposit is the largest feature seen in the model. In plan view, the anomaly seems to be concentrated slightly south of the deposit and has exten sions of lesser chargeable material to the east and the southwest (the eastern extension of chargeability anomaly is also shown in the north-facing cross-section). The deposit does have mineralized extensions to the east and southwest, and although the elevation of this mineralization does not quite correspond with the chargeability model, the model maybe sensitive to these areas. Of course, the distribution of chargeable material may not be one-to-one with ore. Chargeable minerals could have been deposited within the surrounding host-rocks that might result in anomalies peripheral to the deposit. The modeled concentration of chargeable material also extends deeper than the extent of the deposit, as defined by drilling. This is probably a result of having no depth constraints within the data themselves. However, the overall results suggest that the assumptions made in order to invert the data using MAG3D, were appropriate in this case. A cluster of smaller chargeability anomalies are located to the northeast of San Nicolas centered at about 500m east and -200m north. This is in the vicinity of El Salvador, a smaller lense of massive sulphide hosted within the same volcanic rocks13. The only other anomaly of significant magnitude is a small, shallow discrete body located at about -1200m east and -2800m north. This anomaly has a magnitude of about 20 msec and is located on the edge of the survey. This anomalous region of chargeable material coincides with regions of anomalous density and susceptibility and does contain mineralization. An analysis of this anomaly is presented in the next section — physical property correlation. 13E1 Salvador was discovered prior to the discovery of San Nicolas. Chapter 3. Regional Scale Modeling 78 The only anomalies that have been modeled correspond to mineralization. However, away from this area, chargeability anomalies have been associated with some of the Tertiary rocks. In a more widespread survey, other methods such as gravity would be needed to distinguish between mineralization and Tertiary rocks. 3.8 Physical property correlation The regional physical property models that have been generated contain three-dimensional structures that are, hopefully, easier to interpret and relate to geology than the original data. However, combining multiple sources of information, in order to target areas of favorable physical properties, becomes more involved. As a result, models are evaluated in three-dimensions using two simple spatial correlation methods. For each correlation method the physical property models are first imported on to a common mesh. This allows for direct correlation between cell values of different property types. For the first method, cut-off values for the physical properties are determined based on the expected physical property values and the recovered values over San Nicolas. For the density-contrast model a cut-off of 0.13 g/cm3 is used and a cut-off value of 0.0026 S.I. is used for the regional susceptibility model. The models have cells removed that correspond to values below these cut-off levels. As a result only the cells with values above the cut-off remain. These relate to volumes of high density-contrast and magnetic susceptibility. A model is then generated, on the same mesh, that represents the intersection of regions of high density-contrast and high magnetic susceptibility. The resulting cut-off correlation model should be useful in locating more massive sulphides because it represents regions with density and magnetic susceptibility values similar to that of the deposit. Figure 3.19 shows pro jection-to-surf ace plan-views of density contrast, magnetic sus ceptibility, and correlation models of the upper 1000m. Only the upper 1000m of the Chapter 3. Regional Scale Modeling 79 -5000 •sir -10000 5000 -5000 0 Easting (m) 5000 -5000 -10000 -5000 0 Easting (m) 5000 (a) Cells with a density contrast greater than (b) Cells with magnetic susceptibility greater 0.13 g/cm3 than 0.0026 S.I. 5000 I 0 t! o -5000 II \J 1 / Q I r . . N ~r\ ( -10000 -5000 0 Easting (m) 5000 (c) Cut-off correlation model: regions high in density and magnetic susceptibility. Figure 3.19: Plan-views of density contrast, magnetic susceptibility, and cut-off correla tion models. The upper images show the projection-to-surface of cells, within the first 1000m from the surface, that exhibit high density contrast and magnetic susceptibility values. The lower image shows regions that are high in both density and susceptibility. The image is a projection to surface of the volume intersection of the upper models. A blue border shows the extents of the correlation analysis. Chapter 3. Regional Scale Modeling 80 earth is considered, because finding ore below that depth is not within the objective of most exploration programs. Cells below a value of 0.13 g/cm3 in the density contrast model are invisible, as are cells below a value of 0.0026 S.I. in the magnetic susceptibility model. The cut-off correlation model (Figure 3.19(c)) shows five separate anomalies,one of which coincides with the deposit at -1600m east, -400m north. The most northern anomaly lies coincident with out cropping mafic volcanics. The small (one model cell), shallow anomaly about 750m to the west-southwest of San Nicolas, upon field inspection, proved to be a small unmapped outcrop of unaltered intrusive mafic material. The large anomaly to the west is the only dense-susceptible body that extends beyond a depth of 1000m. This is the intersection of the large deep intrusive features seen in the regional density-contrast and magnetic susceptibility models. The remaining anomaly is located at -1300m east, -2700m north within the first 100m from the surface, and is coincident with the chargeability anomaly mentioned in section 3.7. Three drill-holes have been completed each to a depth of 450m around, but not directly into, this anomaly. Of the three drill-holes only one encountered mineral ization similar to that found at San Nicolas14. If the chargeability model was correlated with both the susceptibility and the density models, the resulting three-way correlation model would contain only two anomalies; one being the deposit and the other being this small occurrence of mineralization. This is an important result and demonstrates how the inversion models can be used together to find mineralization. However, the known mineralization at El Salvador has not been found with this method. The method presented above has been successful in reducing the number of targets, however, it is dependent on the cut-off values and doesn't indicate whether one area is more dense or susceptible than another. Consequently, a second correlation method 14Hole SAL-61 intersected 0.7% copper at 430m depth and 0.2% zinc at 207m. Chapter 3. Regional Scale Modeling 81 was employed. For the second method, the density-contrast and magnetic susceptibility models are normalized to each of their range of values. This means that high density and high susceptibility regions will have a value of one, while low density and susceptibility values will be represented by zero. The product of these models is then calculated. The resulting model (Figure 3.20) will have high values that corresponds to high density and susceptibility values. Figure 3.20(a) shows a perspective view of the product of normalized density-contrast and magnetic susceptibility values. Again, the number of anomalies has been significantly reduced with this method. The model has similar features to that shown with the pre vious correlation method, however one important difference is present. A volume has been produced to the northwest of the deposit that corresponds to very high magnetic susceptibility but only moderate density-contrast values. This is a limitation of this method, however most of the other anomalies produced are from a combination of the two physical properties. The deposit along with a dense mafic intrusion to the north and a deep intrusion to the south have all been modeled. In addition, the small anomaly, located at -1000m E -2700m N, which also corresponds to the chargeability anomaly, is also present. Both of the correlation methods work well at San Nicolas because all of the physical properties considered give high values over the deposit. If low or intermediate physical property values are of interest then further manipulation of the models would be required. 3.9 Conclusions Inversion of regional geophysical data sets has produced large-scale physical property distributions of the earth. Together with physical property values and an understanding of the geologic setting, these models have been interpreted and given a geologic context. Chapter 3. Regional Scale Modeling 82 North^f 3500 (a) Northwest-facing perspective view of density-susceptibility correlation values. The model has been volume-rendered at a value of 0.2. (b) Projection-to-surface of density-susceptibility correlation values within the first 1000m. Cells with a value of less than 0.2 have been removed. A blue border shows the extents of the correlation analysis. Figure 3.20: Normalized density-susceptibility regional correlation values. The values are the product of normalized regional density-contrast and magnetic susceptibility models. Chapter 3. Regional Scale Modeling 83 This interpretation would be much harder, if not impossible, without considering three-dimensional physical property distributions. Geologic features such as the extents of a horst, the detection and definition of a deep rooted intrusion, the identification of faults of two different ages, and the detection of shallow bodies that could be mafic intrusions have all contributed to the regional geologic understanding of the area. In addition to the regional geologic interpretation, the models have been correlated in order to directly detect mineralization. The correlation of density and susceptibility models results in a reduced number of regions that might warrant further investigation. When combined with results from the inversion of gradient array chargeability data only two anomalous regions would remain, both of which contain mineralization. This method of selecting targets would hopefully be as effective in areas further from the deposit as it is here. False anomalies can occur from all of the methods used, as a result there isn't one combination of physical properties that will detect mineralization without missing some that is already known. As a result different priorities of targets could be determined. These could be based on the different intersections of dense, magnetic, and susceptible material. If a comprehensive conductivity model could be constructed for the same region the correlation results could change and a better discrimination method might result. Within the context of our exploration program, regional inversion modeling of geo physical data has been successful in targeting areas for detailed work based on geologic interpretations and direct detection methods. An integrated geophysical approach has been shown to be most efficient at selecting targets. When combined with results from other exploration methods such as geological mapping, and geochemical surveys, the target selection process could be even more refined. In addition, the knowledge obtained from this work can be used to update the explo ration models and thus refine exploration directives. Chapter 4 Deposit Scale Modeling 4.1 Introduction The work presented in chapter three was successful in directly detecting mineralization and contributing to the regional geologic picture. From the correlation methods presented possible targets could be identified. Ground surveys are now conducted in order to provide more information about these targets. The flowchart in Figure 4.1 shows how the work presented in this chapter is integrated into the second part of an exploration program, and builds on the knowledge already obtained. This chapter deals with detailed geophysical modeling of the deposit, using finely dis cretized meshes, for further physical property characterization and enhanced definition. Results from the regional modeling are used for removal of regional trends from the grav ity data prior to inversion, and for constructing reasonable starting and reference models when inverting the other data. Along with the gravity and gradient-array IP data, which was used for the regional modeling, more localized data sets are considered, including ground magnetics, CSAMT, and realsection-array IP. The goal is to invert the data to see how much information about the deposit can be recovered. Unlike our example exploration program, we are in the position to check the inversion results of this chapter with geologic sections that have been interpreted from drill-hole information. However, the models are constructed through inversion of the geophysical data alone. 84 Chapter 4. Deposit Scale Modeling 85 Exploration Goal: •find economic mineralization \ Geologic Constraints: j existing mapping | •historical workings •tectonic setting Geologic •descriptive model Models: •genetic model -> Exploration Model Regional Geologic Mapping: •Outcrop geology •Structure •Alteration Geochemistry •Stream sediments Remote Sensing •Multi spectral •Hyperspectral No Targets Geologic Interpretation Target Generation Terminate Exploration Program Update Geologic Exploration Model Geologic Logging I Detailed Geologic Mapping: •Stratigraphy •Structure •Alteration Trenching Rock samples •Assays •Alteration mineralogy (X-ray/PlMA) Geochemistry •soil surveys j Geologic Interpretation Mineralization Not Intersected Locate Drill-hole Terminate Eiplo ration Program Mineralization Intersected Goal Realized: Mineral Resource Geophysical Exploration: Deposit Scale Modeling Geophysical .expect gije/depth of target Exploration .physical properties contrasts Model: Regional/Reconnaissance Scale Geophysics: Physical Property Measurements Existing Data Survey Design Collection I Gravity AMAG (+/- ARAD) AEM II IP || Other: VLF, MT Modeling Individual Forward Correlation Inverse Joint/Simultaneous Forward Inverse Conclusions about physical property distributions of the the earth Update Geophysical Exploration Model Local/Deposit Scale Geophysics: 5 Additional Physical Property Measurements Survey Design Data r I Collection Gravity GMAG EM: CSAMT, UTEM DC/IP Other VLF Modeling Individual Inverse Correlation Joint/Simultaneous Forward Inverse Conclusions about physical property distributions of the the earth Physical property core measurements Down-hole methods Delineation Figure 4.1: Geophysical Exploration Methodology Flowchart: Deposit scale modeling. Chapter 4. Deposit Scale Modeling 86 4.1.1 Local geology Figure 4.2 shows the local geology in the vicinity of the deposit as interpreted from drilling. The deposit is hosted in a series of interwoven mafic and felsic volcanic rocks of the Upper Jurassic to Lower Cretaceous Chilitos Formation, which lie unconformably over graphitic mudstones of the Zacatecas Formation. The deposit is almost entirely bounded to the east by a southwest-dipping fault, which could have been a feeder structure. Mineralization continues to follow the fault at depth in an unconstrained part of the deposit referred to as the keel. The deposit is flanked by a thick package of rhyolites to the west that contains a steeply dipping north-south fault (not shown on this section). The volcanic succession that envelops the deposit is overlain by a Tertiary-aged breccia overburden, which varies in thickness from 50 to 150m. The breccia includes tuffs and clasts derived from the underlying volcanics. Outcrops of the breccia have been mapped to the northwest, but an overlying thin veneer of Quaternary alluvium is present in the vicinity of the deposit. Hydrothermal alteration, which may have caused a change in physical property values, is prevalent throughout the deposit and surrounding geology. Figure 4.3 shows a simplified stratigraphic section as interpreted in the vicinity of the deposit. Introduction to inversion results The models constructed by inverting local-scale geophysical data are presented in the following sections. Only the models are presented; the data are displayed in Appendix D. The models are presented in plan, cross-sectional, and perspective views in much the same way as Chapter 3. The geologic sections shown in Figure 4.2 are shown with each inversion result in order to evaluate the success of the inversion. All of the inversion parameters used are provided in Table E.l. Chapter 4. Deposit Scale Modeling 87 -2000 -1600 -1200 Easting (m) (a) North-facing geologic cross-section located at -400m north. -800 -400 0 Northing (m) (b) West-facing geologic cross-section located at -1700m east. Figure 4.2: Geologic cross-sections through the San Nicolas deposit, modified from sec tions provided by Teck Corp. The main geologic units are labeled. Chapter 4. Deposit Scale Modeling 88 Stratigraphic Section u/c u/c • Alluvium Mafic Volcanics Mafic to Intermediate Flows & Breccia Massive Sulphide Rhyolite Flows and Breccias Graphitic Mudstones Quaternary | Tuff, Volcanoclastic Breccia Tertiary S Sedimentary - Mudstone, Limestone Cretaceous U. Jurassic -» L. Cretaceous =3 s w ,5 Figure 4.3: Simplified stratigraphic section in the vicinity of the San Nicolas deposit, modified from [Danielson, 2000]. Chapter 4. Deposit Scale Modeling 89 4.2 Density Inverting the gravity data on a local-scale requires a two-pass procedure. First, a "re gional" field is removed from a subset of the original data, and secondly, the residual data are inverted on a finely discretized mesh, centered on the deposit. The removal of regional signal from the local data is undertaken in order to remove signal in the data that doesn't originate from the localized region of interest, thus allowing the inversion to be concerned with only the deposit and the immediately surrounding volume. 4.2.1 Regional Removal The regional density contrast model created in Chapter 3 is used to generate a regional gravitational field over a 1.8 km x 2.4 km area centered on the deposit. This is done by calculating the response of the regional density contrast model with a volume removed (Figure 4.4). The removed volume encompasses the deposit and two smaller, dense bodies. The resulting regional signal is subtracted from the original observations and residual data are produced that can be used for the detailed inversion (Figure D.ll). Li and Oldenburg [1998a] describe removal of regional magnetic fields using this technique. 4.2.2 Inversion of Residual Gravity Data Data errors of 0.01 mGal were assigned to the residual gravity data and the 3D inversion was carried out with GRAV3D using the parameters listed in Table E.l1. The observed data are reproduced very well by responses from the detailed density contrast model (Figure 4.5) and the errors that were assigned to the data seem to have been appropriate. xThe regional density contrast model used to create the regional response is based on a mesh with large cells (100m x 100m x 50m). Consequently, this is a coarse representation of the earth, and the regional gravity field contains small inaccuracies that will propagate through to the residual gravity data. In order to account for these inaccuracies, padding cells are used in the inversion of the residual gravity. Chapter 4. Deposit Scale Modeling 90 Figure 4.4: Perspective view, facing north-west, of regional density-contrast model used for calculation of the regional signal. The volume that has been removed contains the deposit and two smaller high density anomalies, and is the volume on which the more detailed inversion is focused. The central density contrast anomaly represents the San Nicolas deposit and correlates well with the dense massive sulphides. This is an improvement over the regional model probably because the deposit is defined with small cells, and discretization is not affecting the results. Figure 4.6 shows cross-sections of the density-contrast model overlain with geologic boundaries. The centroid of the anomalous density coincides with the San Nicolas deposit and corresponds well with the dense massive sulphides. The lateral dimensions of the main body of the deposit are reasonably well defined, and the density anomaly doesn't continue to the surface, as it does in the regional inversion. Large density-contrast values are present above the expected elevation of the deposit, and neither the keel nor the thin extension of sulphides to the east have been modeled in this inversion. However, this is a useful result when one considers the lack of depth information contained in the data, and locating a drill-hole based on this model would Chapter 4. Deposit Scale Modeling 91 035 Figure 4.5: South-east perspective view of the detailed, density-contrast model shown as a volume-rendered isosurface plot using a cut-off of 0.18 g/cc. Cells with density contrasts less than this value are invisible. Chapter 4. Deposit Scale Modeling 92 lead to a successful intersection of the ore-zone. Two other density anomalies are also well defined by the model. The north-eastern anomaly was drilled and intersected 400m of intrusive, mafic rock. Field inspection of the western anomaly revealed similar mafic intrusives outcropping at the surface. 4.3 Magnetic Susceptibility Finely sampled, magnetic field data were available for detailed modeling in the form of ground magnetic data. A base-level of 44,000 nano-Tesla was removed from the 614 diurnally corrected data as a processing step prior to inverting. In addition, several data in the center of the survey that were contaminated with cultural noise, such as a fence or steel-cased drill hole, were discarded. The 3D magnetic inversion was performed using MAG3D and the predicted data were in good agreement with the observations (Figure D.12). A volume-rendered representation of the 3D susceptibility structure is shown in Fig ure 4.7. The cut-off value is 5 xlO-3 S. I. A distinct body of higher susceptibility is modeled. The majority of the susceptibility high is coincident with the deposit, however, high values continue to the north (as seen in the data and the regional model). A good correlation between magnetic susceptibility values and geology in the vicinity of the de posit can be seen in the north-facing cross section presented in Figure 4.8. The high magnetic susceptibilities associated with the deposit align well with the boundaries of the sulphide body. A concentration of high susceptibility is present in the western half of the deposit. As with the regional magnetic susceptibility model (Figure 3.10), there appears to be a continuation of moderately susceptible material to the east and below the deposit. Mafic rocks below the ore-bounding fault zone are less altered, and may be the source of Chapter 4. Deposit Scale Modeling 93 (a) North-facing cross-section through density-contrast model located at -400m north. Geology from Figure 4.2(a) is overlain. 700 I 1 -800 0 Northing [m] (b) West-facing cross-section through density-contrast model located at -1700m east. Geology from Figure 4.2(b) is overlain. Figure 4.6: Cross-sections through deposit-scale density-contrast model at the San Nicolas deposit. Chapter 4. Deposit Scale Modeling 91 1 ,0.008 0 00667 0 00533 0 004 0.00267 0.00133 0 S.I. Figure 4.7: Perspective view of the detailed, magnetic susceptibility model. this susceptible material. There is no indication of the keel or the eastern extension of sulphide in this model. The result is a smoothed version of the true susceptibility values, and it does not contain highly detailed structural information, particularly at depth. A reasonable magnetic susceptibility model has been recovered from the ground mag netic data, however, it appears that obtaining information about the location of sulphides from the model is not a straight forward exercise. Evidence from physical property mea surements (Appendix B.3.1) show that magnetic minerals, while initially expected to be associated with sulphides, are not distributed throughout the deposit. Moreover, they can also be present in host rocks and might have been deposited by hydrothermal events not associated with the deposit itself. This might explain the high susceptibility values that persist to the north, and those that persist to the east and below the deposit, and the lower values that are modeled in the eastern half of the deposit. Chapter 4. Deposit Scale Modeling 95 (a) North-facing cross-section through susceptibility model located at -400m north. Geology from Figure 4.2(a) is overlain. Figure 4.8: Cross-sections through deposit-scale magnetic susceptibility model at the San Nicolas deposit. Chapter 4. Deposit Scale Modeling 96 4.4 Conductivity Data Three lines of controlled source audio-frequency magneto-telluric (CSAMT) data were collected over the San Nicolas deposit. The CSAMT method is an electromagnetic tech nique that uses a grounded dipole current source2. At San Nicolas, CSAMT data were collected with the receiver lines orientated perpendicular to perceived geologic strike. It is not obvious from the observed apparent resistivity or phase data (Figures D.13 and D.14) that a conductive ore deposit is present. Ohm-m Figure 4.9: Perspective view of resistivity model created by interpolating between lines of "stitched" ID models. The inversion has successfully discriminated between the over burden and the deposit. Apparent resistivity and phase data are inverted to recover a one-dimensional, resistiv ity model beneath each station. The one hundred and eighty ID models were "stitched" together to form a 2D model along each of the three lines of data. An interpolation rou tine was used to fill in the volumes between the lines and a three-dimensional resistivity 2An overview of the CSAMT method is given in Appendix C.4. Chapter 4. Deposit Scale Modeling 97 model was produce. A volume-rendered image of the model, with an isosurface cutoff of 30 Ohm-m, is shown in Figure 4.9. The large conductive feature at depth is the deposit. It is separated from a conductive overburden by an intervening layer of higher resistivity. The cross-section of the resistivity model, shown in Figure 4.10, shows that the deposit boundaries interpreted from the inversion agree reasonably well with those inferred from drilling. Easting [m] Fi gure 4.10: Cross section of resistivity model with geology. The resistivity model was produced by "stitching" together 60 ID inversion models. 4.4.1 Testing the 3D resistivity model In order to test the validity of the one-dimensional inversion, the "stitched" inversion models and a 3D physical property model3 were forward modeled using EH3D (see Ap pendix E.8). Apparent resistivity and phase data were calculated from electric and mag netic field measurements and are plotted in Figure 4.11. The forward-modeled data from 3The physical property model was generated by digitizing the geologic sections, assigning conductivity values based on the physical property measurements, and interpolating between the sections to create a 3D model (Appendix B.5). Chapter 4. Deposit Scale Modeling 98 the inversion models replicated the original data quite well at higher frequencies (Fig ure 4.11(b) - left plots). At low frequencies, the 3D conductivity model produced data that came closer to the observations than the "stitched" inversion models (Figure 4.11(a) - right plots). This result is plausible when one considers the geometry of the deposit and the volume of the subsurface that low frequency measurements will sample. At low frequencies the measurements will sample a larger volume, including the 3D deposit, and this information is hard to represent in only one dimension. Conversely, at high frequencies the part of the subsurface that the measurement has sampled will be more one-dimensional. Since the high frequency data are sensitive to shallow features, this suggests that structure in the upper region of the "stitched" inversion models is valid, and the structure of the physical property model is more representative of the deposit at depth. Discrepancies between the observed data and the forward modeled physical property data possibly arise due to inaccurate assignment of conductivity values. The forward modeling mentioned above suggests that this data set might produce a model with greater resolution at depth if it were to be modeled in three dimensions. With the inversion codes available, the raw electric and magnetic field data, from which the apparent resistivity and phase data were derived, are needed in order to invert this data in 3D. Unfortunately, at this time, these raw data have not been recovered. Of specific interest in a 3D inversion would be detection of the deep high-grade keel. A sensitivity analysis on the detectability of the keel using different CSAMT and CSEM (Controlled Source Electromagnetic) survey designs is presented in Section 6.4. 4.5 Chargeability In this section "Real-section" array IP data is simultaneously inverted along with a subset of the gradient array data in order to produce a chargeability model that has more detail, Chapter 4. Deposit Scale Modeling 99 1000 ^ a; 1 -D Inversion Model Physical Property Model 30 600 200 1000 is S ^3 a* 30 600 -2500 Easting -500 -2500 Easting -500 (a) 8Hz data 1-D Inversion Model Physical Property Model 1 1 r ; 5S g 200 -2500 Easting -500 -2500 Easting -500 (b) 256Hz data Figure 4.11: Comparison of CSAMT field data with forward modeled data from the stitched 1-D inversion resistivity model(left plots) and a physical property resistivity model (right plots). The field observations are displayed in red while the modeled data is shown in blue. At 8Hz (upper plots), the physical property model fits the field data well while the 1-D inversion model fits the field data better at a frequency of 256Hz (lower plots). Chapter 4. Deposit Scale Modeling 100 256Hz 8Hz -2000 -1600 -1200 Easting (m) Figure 4.12: North-facing, cross-sectional diagram of the geology at San Nicolas showing the skin-depths associated with propagating plane-waves at frequencies of 256Hz and 8Hz. The skin-depths were calculated for a 25 Ohm-m whole space. The higher frequency signal samples geology that is much more "one-dimensional" than that sampled by the low frequency signal. Chapter 4. Deposit Scale Modeling 101 and better depth information, than the regional model that was produced using MAG3D in Section 3.7. 4.5.1 Chargeability data In real-section surveys, the voltages from consecutive dipoles are measured and plotted in a pseudosection format. The current transmitter electrodes straddle the potential electrode array and the transmitter electrode spacing is continually reduced to change current flow in the subsurface. Data are recorded only at those potential electrodes lying interior to the transmitters. The configuration is a modified version of that employed in gradient-array surveys but should have better depth information due to the increase in number, and varied spacing of the transmitter electrodes. The calculated chargeability is plotted beneath the receiver dipole and each row of the resultant pseudosection-type plot corresponds to a different location of the current elec trodes. The final plot has an inverted appearance compared to pseudosections obtained with more traditional pole-dipole or dipole-dipole plots. The Realsection plots for San Nicolas are shown in Figure D.15. It is apparent from the data that a chargeable body is present, but these pseudosections provide no tangible information about the depth of the target. That can only be obtained through inversion. The inversion procedure for chargeability is a three-step process. First the resistivity of the model volume needs to be calculated. Next, the sensitivities needed for the IP inversion are calculated using the resistivity model. Lastly, the IP data are inverted using IPINV3D to recover a 3-D chargeability model. The details for the IP inversion are provided in Table E.l. 4.5.2 The resistivity model DC resistivity data are collected at the same time as IP data and, as codes are available to invert these data in 3D, they are usually used to construct the resistivity model needed Chapter 4. Deposit Scale Modeling 102 for the sensitivity calculation when inverting IP data. At San Nicolas, these data do not produce a resistivity model that corresponds to the expected resistivity distribution. A surficial conductive layer was modeled that resembles the overburden, but no conductive body is recovered in the vicinity of the deposit. One possible explanation could be that the currents emitted from the transmitters do not penetrate into the volcanic rocks and are only concentrated within the more conductive overburden; in other words, taking the path of least electrical resistance. This would seem to make sense, although some currents would have to penetrate to the deposit in order for a chargeability response to be recorded (assuming the deposit is the only chargeable source)4. Only an approximate conductivity structure is needed for the sensitivity calculations, so chargeability inversions were performed using sensitivities calculated from the inversion of DC resistivity data. This produced a chargeability model that corresponded well with the deposit, but seemed to have layered artifacts corresponding to the Real-section transmitter locations. However, other, more accurate, conductivity models were available; the conductivity model produced from the inversion of C.S.A.M.T. data, and a model constructed by digitizing the geologic sections and assigning conductivities based on physical property measurements. Both of these conductivity models were used to construct sensitivities for the inversion of the IP data. The resulting models did not have the artifacts observed before. The model produced using the CSAMT conductivity model is presented in this section as it contains no apparent artifacts and also is derived only from geophysical data collected at the surface. 4The phenomenon! of a disconnected conductive and chargeable body, lying below a conductive overburden, producing a negligible DC resistivity yet observable IP response has been observed before and is documented by Brescianini et al. [Brescianini et al., 1992]. Chapter 4. Deposit Scale Modeling 103 4.5.3 Removal of electrode effects When inverting IP data using IPINV3D, especially data that have limited inherent depth information, such as gradient-array data, accumulations of chargeable material are often observed around the transmitter electrodes. This "electrode-effect" often leads to the recovery of bodies with smaller chargeability values than expected. For this reason, an empirical weighting function is applied in the inversion of the San Nicolas data that penal izes against chargeability values near the transmitter electrodes. The weighting function is included in the spatially dependent weighting matrix (ws) introduced in Equation 2.2. It is designed to strongly weight the model close to a reference model of zero at the transmitters. Appendix G shows how this weighting function has been derived using synthetic models. The result of implementing this weighting function on the model is shown in Fig ure 4.13. Most of the surficial chargeability has been reduced. The magnitude of charge-ability material at the anomaly that corresponds to the deposit has doubled from about 25 msec to 50 msec. An alternative method to using this weighting function would be to weight the sensi tivity matrix used in the inversion. This could make cells near the transmitter locations less sensitive to the data, and cells further away more sensitive to the data. 4.5.4 The chargeability model Figure 4.14 shows a perspective view of the recovered chargeability model. The spatial distribution of the high chargeability values agrees quite well with the location of the deposit. The high chargeability values are off-set slightly to the south in relation to the deposit (this is also seen in the observed real-section data in Figure D. 15). The high chargeability values present at the surface are probably remaining artifacts of the Chapter 4. Deposit Scale Modeling 104 -3000 -2500 -2000 -1500 -1000 -500 0 Easting [m] (a) Plan-view of recovered chargeability model at the surface. The black circles show the location of the transmitters; these locations correspond with most of the high chargeability material that has been recovered at the surface. -3000 -2500 -2000 -1500 -1000 -500 0 Easting [m] (b) Plan-view of recovered chargeability model at the surface. The black circles show the location of the transmitters. Much of the surficial chargeable material that was modeled has been subdued with a weighting function. Figure 4.13: Plan-views of the recovered chargeability models at the surface. The high chargeability values that correspond with electrode positions have been reduced through the use of a weighting function. Chapter 4. Deposit Scale Modeling 105 Figure 4.14: Perspective view of the volume-rendered chargeability model with a cut-off value of 40 msec. electrode-effect mentioned above, although they could be real. By viewing a north-facing cross-section (Figure 4.15), one can see a concentration of chargeable material centered on the south-west dipping fault. An exploration of model-space was conducted by changing the inversion parameters and the chargeable material around the fault was consistent for all models produced. The fault zone is known to contain semi-massive sulphides and stringer zones, which contain high-grade veins mixed in with host rock. It might be possible that this fault zone could be more chargeable than other parts of the deposit. 4.6 Conclusions All of the models produced at the deposit-scale build on the understanding obtained at the regional scale. The density model has been refined, through removal of a regional signal and finer discretization, in order to produce an anomaly that better represents the deposit. A much more detailed magnetic susceptibility model has also been produced. Chapter 4. Deposit Scale Modeling 106 -2200 _800 Easting [m] (a) North-facing cross-section through chargeability model located at -400m north. Geology from Figure 4.2(a) is overlain. (b) West-facing cross-section through chargeability model located at -1700m east. Geology from Figure 4.2(b) is overlain. Figure 4.15: Cross-sections through deposit-scale chargeability model at the San Nicolas deposit. Chapter 4. Deposit Scale Modeling 107 However, the new model raises questions about the relationship between the distribu tion of magnetic minerals and massive sulphides at San Nicolas. The distribution of magnetic material is offset from the distribution of ore in what can be described as an over-printing effect. This over-printing effect suggests that the magnetic content of San Nicolas was produced after the deposition of the sulphides in a later hydrothermal event. Layered-earth inversions of the CSAMT data have produced a conductivity model that differentiates between the deposit and the overburden. When the models are "stitched" together to produce a 3D model, they have been shown to fit the data well at high fre quencies. The simultaneous inversion of gradient-array data with "Real-section" IP data has produced a chargeability model that defines the deposit well, and has high values co-located with the main south-east dipping fault zone. These models were all constructed by inverting the geophysical data alone, and information in the form of physical prop erty measurements and geologic interpretations have only been used to demonstrate the success of these models, and the validity of the techniques that produced them. This work clearly demonstrates how geophysics can be used as part of a more focused exploration effort to spot drill-holes. All of these methods have proved successful in defining the deposit at the correct depths and lateral extents. Furthermore," the models start to produce more detailed information regarding the geometry of the deposit. The gravity inversion seems to be sensitive to the keel and the chargeability model produces increased values around the fault-zone. The next step in the exploration program would be to use information that might result from drilling these well defined targets. The next chapter addresses this in an attempt to further refine the models and delineate more ore by including physical property information directly in the inversion procedure. Chapter 5 Detailed Modeling: Incorporating a priori information 5.1 Introduction Inversion of surface geophysical data to produce deposit-scale physical property models has been successful in recovering the general features of the San Nicolas deposit (as shown in Chapter 4). These models could have been used to successfully direct a drilling program that would have resulted in the discovery of San Nicolas. The next stage in the exploration process is to further delineate any mineralization that may have been found. Geophysical inversion can play a role at this delineation scale, as outlined by the shaded elements of the flowchart presented in Figure 5.1. If a target is drill tested, and whether or not mineralization is intersected, each drill-hole can provide additional information that can be used to further refine the inversion models. The inclusion of this data-independent information will help eliminate models inconsistent with known geology that fit the data [Scales and Tenorio, 2001] The additional information can be in the form of analysis performed on the drill-core itself, down-hole geophysical measurements, or even interpreted geology that can be assigned physical property values (this last option would obviously not provide as much information unless the physical properties were very constrained before hand). These updated models can, in turn, guide further drilling. This iterative process can expedite 1In a similar maimer additional information would be an important addition to inversions performed on in-mine data sets where mine geologists or mining engineers might require detailed information about parts of the earth where drilling has not yet reached. Again, the inclusion of known physical property information would constrain the result and focus the inversion on areas of interest. 108 Chapter 5. Detailed Modeling: Incorporating a priori information 109 a drilling program and more effectively define a resource. The extensive drilling program that was completed at San Nicolas (Figure 5.2 shows the location of drill-holes that make up over 25,000m of drilling2) presents an opportunity to improve upon the deposit-scale models by including additional information about the subsurface in the inversion process. The idea of including drill-hole information to constrain an inversion seems entirely acceptable, however, some major issues are raised when this put into practice. The issues include how to incorporate detailed point measurements in a large-scale model, and how to perform inversions that are consistent with the drill-hole information. This chapter presents a first look at these issues. The results demonstrate the potential benefits that may arise from including information from one drill-hole, and those that may arise from including information from many drill-holes. The process of generating these models identifies potential research areas. 5.1.1 A priori information at San Nicolas As well as providing geologic information about the subsurface, core samples were mea sured for density and magnetic susceptibility physical property values. Density measure ments were made on sulphides from most drill-holes in order to estimate mine tonnages. The susceptibility measurements were made by myself during a visit to San Nicolas in February of 2001. Over one thousand magnetic susceptibility measurements were made on core samples from thirteen drill-holes. The measurements were made on all rock types in the local geologic section. These physical property measurements are considered to be a source of explicit a priori information, that is, information that can provide direct knowledge about the physical 2The drill-hole notation of SAL represents Salvador—the San Nicolas discovery came from a drilling program that was initiated on the El Salvador project. Chapter 5. Detailed Modeling: Incorporating a priori information 110 Exploration Goal: *flnd economic mineralization i Geologic Constraints: > 'existing mapping | •historical workings | 'tectonic setting Geologic -descriptive model Models: 'genetic model => Exploration Model j Regional Geologic j Mapping: |. "Outcrop geology •Structure •Alteration Geochemistry •Stream sediments Remote Sensing •Multispectral •Hyperspectral No Targets Geologic Interpretation Target Generation Terminate Exploration Program i ^ i Update Geologic j-- j Exploration Model j Geologic Logging 4 • Detailed Geologic Mapping; •Stratigraphy •Structure •Alteration Trenching Rock samples •Assays •Alteration mineralogy (X-ray/PlMA) Geochemistry •soil surveys Mineralization Not Intersected Geologic Interpretation Locate Drill-hole Terminate Exploration Program Mineralization Intersected Goal Realized: Mineral Resource Geophysical Exploration: Using Drill-hole Information Geophysical .expected5ize/depthoftelget Exploration .physical properties contrasts Model: Regional/Reconnaissance Scale Geophysics: Physical Property Measurements Existing Data Survey Design Collection Onyity AMAG (+/• ARAD) Other VLF, MT Modeling Individual Inverse Joint/Simultaneous Correlation Conclusions about physical property distributions of the the earth Update Geophysical Exploration Model Additional Physical Property Measurements in. Local/Deposit Scale Geophysics: Survey Design Data 1 1 1 1 1 Collection Gravity GMAQ EM: CSAMT. UTEM DC/IP Other: VLF Modeling ^Individual Forward Inverse Correlation Joint/Simultaneous Forward : Conclusions about physical property distributions of the the earth : Physical : property core measurements: Down-hole methods x Delineation Figure 5.1: Geophysical Exploration Methodology Flowchart: Using drill-hole informa tion. Chapter 5. Detailed Modeling: Incorporating a priori information 111 property distribution prior to inverting the data. For this study, a priori information is incorporated in the inversion of gravity and magnetic data by directly constraining the range of values model cells can assume. For the inversion of gravity and magnetic data, the a priori information is used in two stages. First, information from one drill-hole that intersects the deposit is used to constrain the inversion. This is done to determine if constraining part of the deposit can improve definition to other parts of the deposit. Secondly, all the available drill-hole information is included. This is performed in an attempt to identify any sulphide ore that has not been already drilled. If the physical property distribution of the deposit is well accounted for, it is hoped other parts of the model might be more accurately represented, thus improving the chances of finding more ore. Including a priori information in these two scenarios demonstrates the approaches that might be used at different stages of a drilling program. 5.2 Inversion of gravity observations Density measurements were made on core samples of the sulphide ore, and some of the host rock material, that came from all of the drill-holes that intersected sulphide ore, as well as a few other holes. This is performed for reserve estimates, however, it provides an excellent source of information for constraining inversions. Measurements were made about every 1.5m along the core and a plot of density verses depth, with geology and assay results, for drill-hole SAL-31 is shown in Figure 5.3. The higher density values clearly correspond to the massive sulphides3, which have an average value of about 4.3 g/cm3, while the density of the mafic rocks below the deposit are about 2.8 g/cm3. 3 The copper and zinc assay results show a change in ore-grade with depth, having larger concentra tions of zinc at the top of the deposit and more copper present towards the base. The copper assay values seem to display a negative correlation with density; this probably reflects the relative concentrations of pyrite and chalcopyrite. Chapter 5. Detailed Modeling: Incorporating a priori information 112 O) O 0 -100 -200 -300 -400 -500 --600 --700 -800 -2200 f f C^rf^W^ € €' 4 t 4 f <4 ^ y * _c£ .1* -2000 -1800 -1600 -1400 -1200 -1000 Easting [m] Figure 5.2: Plan-view of drill-hole locations. Density measurements on core from drill-hole SAL-31 are used to set bounds on the model values in Section 5.2.2 (shown by a red square), while measurements made on core from all of the drill-holes shown (apart from SAL-90 and SAL-91) are used to provide a priori density contrast infor mation in Section 5.2.3. Susceptibility measurements made on drill-cores from SAL-69 are used to set bounds in Section 5.3.3 (shown by a green triangle), while bounds are assigned in Section 5.3.4 based on magnetic susceptibility values made on cores from thirteen drill-holes that are denoted by a blue diamond. Chapter 5. Detailed Modeling: Incorporating a priori information 113 Discussion of bounds The density measurements are used to constrain the density inversion model through the use of bounds that limit the range of values a cell can assume. The thought of bounding a 25m x 25m x 25m cell from values measured along drill-core that is five or six centimeters in diameter raises questions about the reliability of the process. There is obviously a discrepancy between the information we have and the information we require in order to constrain the inversion. This discrepancy is both a scale issue, and an issue relating to lateral inhomogeneities. The information recorded in fine detail by the core measurements cannot be ade quately represented by one cell value. In order to alleviate this problem we have two options. Either we can reduce the size of the cells or we can try to scale-up the drill hole information. Wherever possible it would be better to reduce the size of the cells so there is less inconsistency in scale between drill-hole information and the model. A smaller cell will probably contain a smaller range of values as sampled by a drill hole, and will probably be less likely to have large lateral inhomogeneities. Unfortunately, the size of the model being used at San Nicolas was already large (close to half a million cells) so the decision was made to scale-up the drill-hole information. In doing so we lose valuable information. 5.2.1 Setting bounds on the density-contrast model To bound a model, we require upper and lower bounds for each model cell. Drill-core den sity measurements provide very detailed and localized information (values every 1.5m) whereas cells within the density-contrast model exhibit constant values throughout 25m x 25m x 25m volumes. Therefore, processing of these core measurements is necessary in order for this information to be consistent with the density-contrast model, and to Chapter 5. Detailed Modeling: Incorporating a priori information 114 a Q Massive Sulphide Quartz Rhyolite and Mafic Volcanics Density [g/cc] Copper [%] 500 p 1 2 : 1 4 5 1 Zinc [%] 10 20 3.0 100 200 300 400 500 Figure 5.3: Drill-hole SAL-31. Displayed as a function of depth from left to right: 1) Geology, 2) density (g/cm3), 3) percent copper, and 4) percent zinc. The average density of the sulphide is 4.2 g/cm3 and the density of the rocks below the deposit is approximately 2.8 g/cm3. Chapter 5. Detailed Modeling: Incorporating a priori information 115 obtain the minimum and maximum density-contrast bounds for each model cell. Fur thermore, the core measurements are absolute density values, while the model values are density contrasts relative to an unknown background value. To that end, inversions were conducted with a background value removed from the core density measurements. Considering the volcanic host rocks, it is expected that a background value of between 2.35 g/cm3 and 2.5 g/cm3 is reasonable. One method of setting bounds on a model cell could be to assign the minimum and maximum core density values that are located within the model cell. It doesn't seem realistic, however, to expect the cell, which has a consistent value throughout, to exhibit these extreme values; a cell value closer to the mean core-density might be more probable. The core-density measurements were sampled regularly throughout the vertical extent of each cell, and from these measurements the mean density value from the core was determined. The upper and lower bounds were calculated by adding and subtracting a small density variation from the mean. This variation could be based on the standard deviation of the measurements. However, the uncertainty in the result would decrease as more samples are used. For independent random variables, the uncertainly associated with the mean would be the standard deviation divided by the square-root of the number of samples. Although, as more samples are used the likelihood of them being independent will probably reduce. In addition, other factors have to be taken into account. We are only sampling the cell at one position in x and y. The drill-hole would sample about three parts-per-million of the model cell (a five centimeter diameter drill-hole intersecting a 25m x 25m x 25m model cell) and, although we can consider vertical density variations within the cell at one horizontal location, we know nothing about lateral density variations within the cell. We might expect horizontal inhomogeneities to exist that could increase the range of values that the cell might have to represent. As a result, the difference between the upper and Chapter 5. Detailed Modeling: Incorporating a priori information 116 lower bounds would probably increase. Apart from a knowledge of the geological setting, we have know way of knowing the possible homogeneities that might exist close to a drill hole. At San Nicolas, the variation from the mean value ranged from 0.05 g/cm3 to 0.2 g/cm3 for different inversion runs. In addition, if the core-measurements only populated half of the cell height, the variation from the mean was doubled relative to cells with core-measurements throughout the cell; and cells with no core-measurements were as signed upper and lower density-contrast bounds of 5 g/cm3 and 2 g/cm3 respectively. These values might be a bit small although, based on the discussion above, it is hard to accurately determine these bounds at such large scales. This is one area where further investigation is needed. Implementing the bounds Once values for the upper and lower bounds have been determined, they are used as inputs in the inversion process. GRAV3D uses a logarithmic barrier method, directly in corporated in the objective function, to implement bounds on the model during inversion. An equation for the objective function with the barrier term is given in Equation 5.1. This function is minimized for subsequently decreasing values of the A while keeping the regularization parameter ((3) constant. When A reaches zero the resulting model honors the a priori information provided. M d?(A) = $d + /JSTO - 2A Y$*iPi - Pfn) + ln(P?aX - Pi)]. C5-1) where: M is the number of model cells, pj is the density of the jth cell, and p™ax and p-rain are tfiQ upper and lower bounds the of the jth cell respectively. Chapter 5. Detailed Modeling: Incorporating a priori information 117 For the subsequent inversions presented, unless otherwise stated, the same inversion parameters where used as those presented in Section 4.2.2 and Table E.l. That is, the same noise levels, reference model, and choice of regularization parameter are used. 5.2.2 A priori density values from drill-hole SAL-31 Using the method described above, density values from drill-hole SAL-31 were used to constrain the deposit scale density-contrast model. The bounds implemented assume a background density value of 2.5 g/cm3, and constrain the model along a vertical column of cells, from 150m to 475m depth. The bounded cells correspond mostly with massive sulphides, and about twenty-five meters of interlayered, quartz rhyolite and mafic vol-canics. The inversion result, displayed as east-west, and north-south cross-sections in Figure 5.4, shows an improvement over the results of inverting the gravity data.without a priori information. The anomaly that corresponds with the sulphide ore has a higher density-contrast throughout most of the deposit, not just around the bounded cells (the density contrast increases from 0.3 g/cm3 to 0.4 g/cm3 75m away from bounded cells). These higher density-contrast values make the anomaly less spread out, and are con sistent with expected values and a priori information. The model has also been well constrained along the fault that marks the lower extent of the deposit, and there is much less mass contained in the model below the deposit than previously modeled. Further more, the small density-contrast anomaly that was thought to correspond with the keel, in the inversion results presented in Section 4.2.2, is again present. The model is more affected by the introduction of a priori information from drill-hole SAL-31. Changes over the original inversion are observed vertically above and below, and laterally adjacent to the bounded cells. This modeling suggests that high density material is not present below the fault, which is in contrast to the previous gravity inversion result that had a more significant continuation of high density-contrast values at depth. Chapter 5. Detailed Modeling: Incorporating a priori information 118 Easting [m] 700 -800 Northing [m] (a) Initial density model: north-facing cross- (b) Initial density model: west-facing cross-section with geology overlaid. section with geology overlaid. o rs-Easting [m] g/cc 0.35 700 -800 Northing [m] (c) Model bounded using density values from (d) Model bounded using density values from drill-hole SAL-31: north-facing cross-section drill-hole SAL-31: west-facing cross-section with with geology overlaid. geology overlaid. Figure 5.4: North- and west-facing cross-sections through recovered deposit-scale den sity-contrast models. The upper cross-sections show the results from the original de posit-scale gravity inversion, described in Section 4.2.2. The lower cross-sections show the results of inverting the gravity data with the model bounded using the density mea surements made on core from drill-hole SAL-31. Geologic contacts have been overlaid. Chapter 5. Detailed Modeling: Incorporating a priori information 119 Away from drill-hole, the model has not been significantly changed, and the next section addresses the results of incorporating all the available drill-hole information in an attempt to further improve the model, and determine the density structure below, and to the side of where drill-holes have sampled the earth. Summary of using a single borehole to constrain a gravity inversion The models presented above show the results of inverting gravity data with additional information obtained from a single drill-core. The results have improved over the previ ously inverted model and now we investigate the use of multiple drill-hole information. 5.2.3 A priori density values from all drill-holes Figure 5.5 shows the north-facing interpreted geology section at -400S with the location of drill-holes. These drill-holes, along with many others to the north and south, are used to bound the density-contrast model during the inversion of gravity data, by providing explicit a priori density information. This is performed in order to account for the density distribution of the deposit and, hopefully, enable other parts of the model to be more accurately defined, particularly any extension to the keel. The same method of bounding the model cells as described previously, is employed. Different inversions were performed using background densities of 2.35 g/cm2, and 2.5 g/cm3, and bounds calculated with deviations of between 0.05 g/cm3 and 0.2 g/cm3 above and below the mean value for each cell. All the results were very comparable, showing no major differences, and the results of using a background density of 2.35 g/cm3 with a 0.1 g/cm3 deviation from the mean are presented in Figure 5.6. As expected, the deposit is now very well defined by the density-contrast model, and the density structure matches the interpreted geology. Also, the density anomaly that corresponds with the deposit exhibits even greater density-contrast values, which is Chapter 5. Detailed Modeling: Incorporating a priori information 120 Figure 5.5: North-facing cross-section of interpreted geology with drill-hole locations. Density measurements made on core from the drill-holes were used to bound model cells during the inversion of gravity data. Only the drill-holes along line -400m S are displayed, although many more were used to bound the model. Chapter 5. Detailed Modeling: Incorporating a priori information 121 consistent with the a priori information. Below the penetration depth of the drill holes the model is very uniform; giving more credence to the absence of high-density sulphides at depth. Apart from where a priori information has been explicitly incorporated, the keel has not been well defined as a distinct density-contrast structure. The smooth, deep, south-west extension to the density anomaly seen in the previous inversions of this data has been replaced by a few cells of high density-contrast. In addition there is no extension of the keel present in this model. This does not mean the keel terminates at the limit of the drilling. The result merely indicates a model can exist that doesn't exhibit an extension to the keel, and meets the following requirements: (1) it fits the data to a reasonable degree, (2) it is consistent with the a priori density information provided, and the way in which it has been implemented. Summary of using multiple borehole to constrain a gravity inversion Density contrast models have been produced by inverting gravity data while constraining the range of values some of the models can exhibit. This has been useful in defining the deposit, but more importantly accounting for the mass of the deposit will increase the resolution away from the drill-holes, especially below the deposit. In this density anomaly that was observed when drill-hole information was not included virtually disappeared when the drill-hole information was included. The recovered model does seem to be discontinuous in some areas. This is mainly due to the fact that in cells next to those that are bounded, the bounds are opened right up to accept almost any reasonable density value. Constraining cells next to those that contain drill-hole information is a philosophical one and relates back to the issue of lateral inhomogeneities that are expected in the geology. A discussion is presented at the end of this chapter that deals with some of the issues regarding the discontinuous nature of the recovered models. Chapter 5. Detailed Modeling: Incorporating a priori information 122 g/cc I 0.35 800--2200 Easting [m] (a) Initial density model: north-facing cross-section with geology overlaid. 700 -800 Northing [m] (b) Initial density model: west-facing cross-section with geology overlaid. 800 Easting [m] (c) Density-contrast model bounded by all drill holes: north-facing cross-section with geology overlaid. 700 -800 Northing [m] g/cc • 0.76 (d) Density contrast model bounded by all drill holes: west-facing cross-section with geology overlaid. Figure 5.6: Cross-sections through recovered deposit-scale density-contrast models. The upper cross-sections show the results from the original deposit-scale gravity inversion, described in Section 4.2.2. The lower cross-sections show the results of inverting the gravity data with the model bounded by density measurements made on core from all drill-holes. Interpreted geologic contacts are shown on each cross-section. Note: the lower cross-sections are shown with a range of density-contrast that is different from those previously presented. Chapter 5. Detailed Modeling: Incorporating a priori information 123 5.3 Magnetic Susceptibility Inversion of ground magnetic data has produced a model that shows part of the deposit is magnetically susceptible, parts of the deposit are void of susceptible material, and some rocks outside of the deposit are also likely to be magnetically susceptible. Over one thousand magnetic susceptibility measurements, made on core samples from thirteen drill-holes, confirm the complexity of the susceptibility distribution (Appendix B.3.1). Figure 5.7 shows magnetic susceptibility values along with geology for core from the four drill-holes: SAL-25, SAL-34, SAL-47, and SAL-70. This drill-hole information, and that from other drill-holes, is used to constrain the inversion of ground magnetic data in an attempt to improve upon the existing models. Firstly, measurements from drill-hole SAL-69 are included, and secondly, all the available measurements from the drill-holes are included. 5.3.1 Core Measurements A KT-9 magnetic susceptibility meter was used to make measurements on drill-core samples at San Nicolas. The instrument gets 90% of its signal from the first 20mm of the sample, can be calibrated for core measurements, and has a maximum sensitivity of lxl0~5 S.I. units. The instrument was also calibrated in free air before each reading. Several test measurements showed that the susceptibility distribution within the de posit could dramatically change from less than lxlO-4 to over 0.1 S.I. within a few centimeters. As time constraints did not permit very fine sampling of the core, I took readings approximately every 3 meters in rocks that were expected to exhibit low suscep tibilities, every meter in rocks with high or variable susceptibilities, and several readings were collected either side of faults and geologic contacts. Tests were conducted over dif ferent rock types to see how effectively this sampling characterized the finer susceptibility Chapter 5. Detailed Modeling: Incorporating a priori information 124 SAL-25 SAL-34 SAL-47 SAL-70 Figure 5.7: Magnetic susceptibility core measurements and geology for drill-holes SAL-25, SAL-34, SAL-47, and SAL-70. Note: the upper-most geologic unit is the Tertiary breccia, although it may appear a similar colour to the sulphide ore, which is shown in red. Reference the geologic section and recognize it as an interpretation. Chapter 5. Detailed Modeling: Incorporating a priori information 125 distribution, the results of which are presented in Appendix B.3.1. 5.3.2 Bounding the magnetic susceptibility model The core susceptibility measurements are used to set upper and lower bounds on model cells, which are implemented using the logarithmic barrier method. As with the density measurements, many core susceptibility readings were used to make up one value for each model cell around which bounds were centered. Unlike the density measurements, how ever, no background level need be removed prior to including this data in the inversion. Also the samples were not evenly collected along the core. The uneven sampling had to be taken into account when calculating the composite susceptibility value to use for each cell (Appendix B.3.2). Upper and lower bounds were centered around the composite susceptibility value. As with the density measurements, the magnitude that the upper and lower bounds vary from the composite value is debatable. Most of the core susceptibility values were around 0.1 x 10~3 S.I. with some values as high as 1 x 10-3 S.I. and a few as high as 100 x 10~3 S.I.. The value that was chosen to vary the bounds by was 5xl0~4 S.I., however, more information might be contained in the inversion by using a bound that was related to the core measurements themselves. 5.3.3 Magnetic inversion using a priori information from SAL-69 Constraining the inversion of ground magnetic data by using bounds on the susceptibility model has produced a model of the San Nicolas deposit that has significant differences with respect to that which was produced without a priori information. Figure 5.8 shows a comparison between the magnetic susceptibility models recovered with and without a priori information from SAL-69. Chapter 5. Detailed Modeling: Incorporating a priori information 126 When the model is bounded, the concentration of susceptible material in the west part of the deposit exhibits higher susceptibility values, occupies less volume, and is tightly constrained between the depths of 275m and 350m. Moderate (2-3 xl0~3 S.I.) susceptibility values are still present below the deposit and the north extension of the susceptibility anomaly is also unaffected by including the a priori information. 5.3.4 Magnetic inversion using a priori information from all drill-holes A priori susceptibility information was used from all twelve drill-holes measured in the inversion of ground magnetic data. The results show a susceptibility distribution that is concentrated in the southwest part of the main ore-zone, the keel, and to the northwest of the delineated deposit. Figure 5.9. The additional a priori information helps constrain the susceptibility distribution within the main ore-zone to a small volume that is located against the rhyolite dome to the west and the southwest dipping fault to the northeast. Also, high susceptibility values have been introduced in the keel that have not been previously modeled. These high susceptibility values in the keel are introduced through the presence of drill-hole information but persist in cells next to the drill-holes where the model was not constrained. This suggests that drill-hole information is improving the model in cells that are not bounded by drill-hole information. A continuation of high susceptibility values persist to the north and moderately high values are present below the deposit. Less-altered, mafic volcanic rocks could be respon sible for the susceptible concentrations below the deposit. Different inversion models were constructed using different alpha coefficients to smooth the model in different directions. The results show that the effect the bounded cells have on laterally surrounding cells is dependent on the relative amount of horizontal (x-component and y-component) weighting that is introduced. This could be used to construct more coherent models when geologic stratigraphy can be assumed. Chapter 5. Detailed Modeling: Incorporating a priori information 127 xl O3 S.I. Easting [m] -800 xlO'S.I. Northing [m] (a) Initial magnetic susceptibility model: north-facing cross-section (b) Initial magnetic susceptibility model: west-facing cross-section SAL-69 Easting [m] -800 700 Northing [m] x!0> S.I. (c) Magnetic susceptibility model bounded by (d) Model bounded by SAL-69: west-facing SAL-69: north-facing cross-section cross-section Figure 5.8: North- and west-facing cross-sections through recovered deposit-scale mag netic susceptibility models. The upper cross-sections show the results from the origi nal deposit-scale magnetic susceptibility inversion, described in Section 4.3. The lower cross-sections show the results of inverting the magnetic data with the model bounded by magnetic susceptibility measurements made on core from drill-hole SAL-69. An outline of geology interpreted from drill-holes is overlaid on each section. Chapter 5. Detailed Modeling: Incorporating a priori information 128 xl O3 S.I. Easting [m] 800 700 -800 Northing [m] (a) Initial magnetic susceptibility model: north-facing cross-section (b) Initial magnetic susceptibility model: west-facing cross-section Easting [m] -800 700 -800 xlO'S.I. Northing [m] (c) Magnetic susceptibility model bounded by thirteen drill-holes: north-facing cross-section (d) Model bounded by thirteen drill-holes: west-facing cross-section Figure 5.9: North and west-facing cross-sections through recovered deposit-scale mag netic susceptibility models. The upper cross-sections show the results from the original deposit-scale magnetic inversion, described in Section 4.3. The lower cross-sections show the results of inverting the magnetic data with the model bounded by magnetic sus ceptibility measurements made on core from twelve drill-holes. Geologic contacts as interpreted from drill-holes are overlaid on each section. Chapter 5. Detailed Modeling: Incorporating a priori information 129 5.4 Discussion on the discontinuous nature of the models When examining Figures 5.4(c), 5.4(d), 5.6(c), 5.6(d), 5.8(c), 5.8(c), 5.9(c), and 5.9(d), it seems apparent that models produced with a priori information are somewhat blocky and discontinuous. This is consistent with the a priori information that reflects discrete geologic boundaries, but it is not consistent with the smooth requirements of the model objective function. In fact, the blockish nature of the bounded model is counteracting the efforts of the model objective function to produce a smooth model. A more consistent way of producing a model would be to have both the model bounds and the model objective function acting in concert; both seeking a model of the same type. This could be achieved by changing the bounds to be smooth, or by changing the model objective function to allow blocky models. In this case, changing the bounds would be dishonoring the a priori information, so the appropriate alternative would be to favor blocky models. This can be accomplished in the following two ways: (1) the Wx, Wy, and Wz weighting matrices could be changed so the model doesn't have to be smooth in the vicinity of sharp discontinuities, as determined from the a priori data, and (2) blocky, or piecewise-continuous, models could be sought through the use of the L\ norm in the model objective function, as discussed in [Farquharson and Oldenburg, 1998]. 5.5 Conclusion The introduction of a priori information into an inversion can improve the resulting model. For the potential field examples presented here, the results have improved around where drill-hole information has been included, and below the known parts of the deposit. All of the results produced fit the data to the same degree to that used in Chapter 4. With the incorporation of information from one drill-hole, expected physical property values were increased over much of the deposit. The density contrast was reduced below Chapter 5. Detailed Modeling: Incorporating a priori information 130 the deposit, while the susceptibility distribution was constrained to a range of depths within the deposit. This demonstrates how a priori core measurements could be used in the early stages of drilling to help direct subsequent holes. Incorporating all a priori information has produced models that define the deposit, in the case of the density-contrast model, and define a localized concentration of suscep tibility in the western part of the deposit, in the case of the susceptibility model. The a priori information is useful is defining the deposit so that parts of the signal represented by the data can be accounted for, and areas that have not been drilled may be more accurately modeled. The gravity inversion has shown a decrease in density values below the deposit that suggests there is no high density sulphide ore below the level of drilling. None of the models presented exhibit an extension to the keel but the data may not be sensitive to such a small, and deep, body (the sensitivity of the data to the keel is addressed in Chapter 6. The issues involved with incorporating drill-hole information have been presented and only a first pass solution has been formulated. It is understood that the work presented here represents a preliminary attempt to resolve these issues but also demonstrates the potential benefits than come about if drill-hole information can be included. Issues that still require investigation include: 1. A priori information can originate from finely sampled measurements made on drill-core. This information is processed in order for it to be consistent with the large model cell volumes that are used to constrain the inversion using bounds. This processing is not trivial. Errors are introduced that relate to the sampling size of the cells, variations of physical property within a cell along a drill hole and lateral variations that can only be assumed. In addition, an estimate of this error is needed to determine upper and lower limits of the bounds. It should also be realized that Chapter 5. Detailed Modeling: Incorporating a priori information 131 the finely sample drill-hole information is itself a coarse representation of the actual susceptibility distribution within the earth (Appendix B.3.1). 2. Removal of regional gravity signals and assignment of appropriate background den sity values are issues that need to be addressed when using density information in gravity inversions. 3. The recovered models that have been bounded have a blocky appearance. This could be due to the inappropriate application of a smooth component model objec tive function when a block model is known to exist. Methods of recovering models that are consistent with known geology has a greater scope than just incorporating bounds in a model. We have to address how the model is constructed. The difficulties in implementing drill-hole information accurately in inversions sug gests that actual down-hole measurements might be better suited to producing inversion models. In such cases all the data will be consistent with each other. In summary, this chapter has addressed the methods that might be employed to incor porate drill-hole information in an inversion. This study reveals important aspects that require further investigation. However, removing ambiguities from an inversion model by constraining a model could have major benefits. Superior models can be produced through the use of just one drill-hole. This can lead to more accurate spotting of drill holes in an iterative exploration program. Chapter 6 Forward Modeling and Survey Design: Detecting the keel 6.1 Introduction Chapter five began looking at the delineation of mineralization. That work is continued here. In the context of an exploration program (Figure 6.1) we address the issue of delineation by trying to conduct surveys that are adequately sensitive to the target. In the case of San Nicolas, detecting a deep ore zone involves modification of survey designs. The main objective of delineation is to define the major features of the deposit. At San Nicolas a major feature is the keel. The keel makes up about one sixth of the volume of the deposit, is located below and to the west of the main ore-zone, and was found through drilling. Despite the success of inverting surficial geophysical data at San Nicolas to define the deposit, the keel was not readily identified (although possibly by the gravity inversion). The primary reason for this is that the signal (anomalous field) from the keel is generally small. We quantify these signals through forward modeling. As will be shown, some of the data collected at San Nicolas were probably not ad equate to detect the keel. This work leads naturally to the question: "What data are required to find the keel?" As a result, the focus of this chapter is in two parts. Firstly, the signal from the keel is quantitatively evaluated. This is performed through forward modeling of synthetic physical property models with and without the keel. The difference between the two responses is compared with assumed noise levels in order to 132 Chapter 6. Forward Modeling and Survey Design: Detecting the keel 133 determine if the response from the keel is observable. The second focus of this chapter is to obtain data that are substantially influenced by the keel. For gravity and magnetic data synthetic measurements are made down-hole. The simulations provide insight about the increase in signal strength as readings are made closer to the keel. For electromagnetic and DC/IP data (active source methods), there is much more flexibility for designing an informative survey. This is because both sources and receivers can be placed close to the keel in order to increase signal strength. The methodology for designing optimum surveys is a separate research project; one that is not addressed in this thesis. However, computing the response of the keel for simple survey design changes is performed. This may provide a foundation on which a more comprehensive analysis of the survey design problem cab be addressed. Finding these small ore-zones1 is more difficult in the absence of a large amount of known sulphide, as it becomes harder to instigate follow-up, detailed exploration. As a result, a method of detecting these from the surface becomes important. Information about the keel is obtained from geology sections that were interpreted from drill-holes. This information, along with linear interpolation methods, is used to construct physical property models, or rock-models (digitized on 25m x 25m x25m cells), on which forward modeling can be conducted. The main ore-zone, containing massive sulphide, semi-massive sulphide, and sulphide stringer zones, is modeled as having a volume of about seventeen million meters cubic meters, and is located between a depth of 175m and 425m. The keel, on the other hand, is modeled as having a volume of about three million cubic meters (one sixth the size of the main zone), and is located between depths of 425m and 550m. This helps put in perspective the different signals that might come from the keel and the main-ore zone, and suggests that different surveys maybe 1 Smaller amounts of ore, that are not economically feasible to mine on their own, could be economi cally viable because of the mining infrastructure that might be in place due to San Nicolas. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 134 s Geologic Constraints: \ 'existing mapping | 'historical workings I 'tectonic setting Exploration Goal: •find economic mineralization 1 Geologic 'descriptive model Models: •gene* mode' -> Exploration Model I } Regional Geologic | Mapping: | 'Outcrop geology | •Structure | 'Alteration Geochemistry •Stream sediments Remote Sensing •Multi spectral •Hyperspectral Geologic Interpretation No Targets 1- Target Generation Terminate Exploration Program Update Geologic Exploration Model Geologic Logging A k I Detailed Geologic Mapping: •Stratigraphy •Structure * •Alteration Trenching Rock samples •Assays •Alteration mineralogy (X-ray/PlMA) Geochemistry •soil surveys Mineralization Not Intersected Geologic Interpretation Locate Drill-hole Terminate Exploration Program Mineralization Intersected Goal Realized: Mineral Resource Geophysical Exploration: Survey design for the delineation of deep ore Geophysical ,tXpK^ size/depth of target Exploration .physical properties contrasts Model: Regional/Reconnaissance Scale Geophysics: Physical Property Measurements Existing Data Survey Design Data Collection Gravity A MAO {H- ARAD) AEM IP Other VLF. MT Modeling Forward Correlation Joint/Simultaneous Forward Inverse Conclusions about physical property distributions of the the earth Update Geophysical Exploration Model Additional Physical Property Measurements Local/Deposit Scale Geophysics: iSurvey Designs Data Collection j Gravity 1 1 1 EMrCSAMT.UTHM DOU> Other VLF Modeling Individual Forward Joint/Simultaneous Conclusions about physical property distributions of the the earth Physical property core measurements Down-hole methods Delineation Y Figure 6.1: Geophysical Exploration Methodology: Survey design Chapter 6. Forward Modeling and Survey Design: Detecting the keel 135 -1100 Figure 6.2: Size and position of the keel relative to the deposit, needed in order to detect the keel. 6.2 Gravity Observations Inversion of gravity observations in Section 4.2 have produced a model that, along with approximately defining the main ore zone, had a distribution of higher density-contrast values in the vicinity of the keel. This suggests that the data contain information about the keel above the noise level. In this section the magnitude of the gravitational response of the keel is investigated by performing synthetic gravity surveys on digitized density rock-models. The results are compared to the estimated instrument sensitivity and noise levels in order to see if gravity readings could be used to constrain an inversion model of the keel. A synthetic survey is conducted that is similar to that conducted in the field, and several down-hole gravity surveys are also performed. As we are not able the change the mass or location of the keel, the only way to Chapter 6. Forward Modeling and Survey Design: Detecting the keel 136 increase the signal from the keel is to decrease the distance from the receiver to the keel. To that end, synthetic down-hole gravity measurements are made in order to determine how near to the keel a drill-hole need pass for a discernable signal to be observed. 6.2.1 Simulation of gravity field observations Two sets of synthetic gravity data are made on different density rock-models. One model contained density distributions of the main ore-zone and the keel, while the other model had density distributions of only the main ore-zone. By comparing gravity observations made on each model one can determine the magnitude of the response of the keel. This response is compared with estimated noise levels of between 0.02 and 0.05 mGal. In order to confidently distinguish the response of the keel in the data, signal-to-noise levels of greater than 3 are assumed to be needed, therefore, a gravitational response from the keel of greater than 0.06 mGal is required. The 3D density rock-model was created by digitizing geologic cross-sections provided by Teck, assigning density values to the rock units, and interpolating between the 2D sections. The density value assigned to each geologic unit is determined from the mea surements made on the core, but small density variations within each unit are not taken into account (a constant value of 4 g/cm3 is assigned to the massive sulphide ore). The density rock-model, which covered only a limited volume, was then embedded within a background of 2.5 g/cm3 material. Two density rock-models were constructed (Fig ure 6.3); one with the keel, as digitized from the geologic sections, and one with all density values below a depth of 425m replaced by with values of 2.5 g/cm3. Replacing the density values below 425m effectively removes the keel from model and, although some minor density variations are also changed during this process, it is assumed that the keel represents most of the anomalous mass at these depths. Five hundred and eighty-five synthetic gravity observations were made on the surface Chapter 6. Forward Modeling and Survey Design: Detecting the keel 137 g/cm3 3.15 (a) Perspective view of density rock-model showing an east-west section at -400m north. I 14 I g/cm3 (b) Perspective view of density rock-model with the keel removed. The model is cut at -400m north to reveal an east-west cross-section. Figure 6.3: The two density rock-models that were used to generate synthetic gravity measurements. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 138 of each density rock-model. The observations, shown in Figure 6.4, have a station spacing of 50m and a line spacing of 100m, which is similar to some of the more closely spaced readings taken at San Nicolas. The difference between the observations made on the model with the keel and those made on the model with the keel removed is shown in Figure 6.4(c). The difference ranges in magnitude from 0.004 to 0.12 mGal and, in comparison to estimated noise levels (between 0.01 and 0.05 mGal), the response from the keel gives a discernable response. Assuming the densities were appropriately assigned to the rock-model, this has the important implication that information specific to the keel could be extracted from surface gravity data, and is consistent with the inversion model, which exhibits an extension to the deposit where drilling indicates the presence of the keel. 6.2.2 Synthetic, downhole, gravity observations Accurate determination of the keel's density distribution would likely come from inverting gravity observations that contain information about the keel significantly above noise levels. If we assume that noise levels are constant, the only other option would be to increase the gravitational response of the keel. Because we are dealing with a potential field, this can only be done by making observations nearer to the source; we do this by performing down-hole measurements. The use of drill-holes in geophysical surveys is a logical progression of survey design at San Nicolas because of the extensive drilling that has occured, and although down-hole gravity measurements are not often carried out in the minerals industry2, they have the potential to produce data with low noise levels because of the quiet environment in which they are made. Six synthetic drill holes were used to collect down-hole gravity measurements, at 10m 2Down-hole gravity measurements are not often performed due to the expense of the instrumentation and the possibility of losing it down a hole. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 139 II i i i i mGal (.III •21S3 -1867 -1550 -1233 -917 -600 -2600 -2183 -1867 1550 1233 -917 EastiriE(m) Esstine(m) (a) Synthetic gravity data produced from density (b) Synthetic gravity data produced from density rock-model. rock-model without keel. Eattlng (m) (c) Difference between synthetic gravity data with keel and synthetic gravity data without keel. Figure 6.4: Synthetic gravity data produced from a density rock-model of San Nicolas. Data are produced from models with and without the keel. The difference is also plotted. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 140 intervals, to a depth of 650m on the rock-model. Figure 6.5 shows the locations of these drill-holes plotted on a west-facing cross-section through the density rock-model. ABCDEF Northing [m] Figure 6.5: West-facing cross-section of density rock-model. The main ore-zone is rep resented by the higher density values in the top half of the model while the keel is rep resented by the anomaly in the lower half of the model. The locations of six drill-holes, within which synthetic gravity measurements are made, are also shown. Drill-hole A is 25m north of the keel and each subsequent drill-hole is positioned 25m further north than the last one. The difference between the gravitational response of the model with the keel and the model with the keel removed is shown for measurements made down the six drill-holes in Figure 6.6. Drill-hole A is located on the north edge of the keel and, at a depth of between 450m and 500m, which corresponds to the depth of the keel, the vertical gravitational response changes from 0.7mGal to about -0.5 mGal. This response is very much above the expected noise levels, and would definitely indicate an anomalous mass. As measurements are made down drill-holes further to the north, the response from the Chapter 6. Forward Modeling and Survey Design: Detecting the keel 141 G, [mGal] A B C D E F Figure 6.6: Down-hole plots of synthetic gravity data made in six drill-holes. The data are the difference between the response of a density rock-model with the keel, and a model with the keel removed. Note the change in horizontal scale for plots C, E, and F. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 142 keel is seen to smooth out and decrease in value to +/- 0.05 mGal in drill F, 250m north of the keel. The magnitude of response from the keel in drill-hole F is close to the expected noise levels. It is likely that gravity data collected in drill-holes further than 150m from the keel would have difficulty in detecting the keel from the noise. 6.2.3 Summary of gravity forward modeling. A simulation of the observations made in the field show that surface gravity data could contain sufficient information for the keel to be determined. This result adds credence to the inversion results obtained in Section 4.2. With assumed noise levels, down-hole measurements should be made within about 150m of the location of the keel. 6.3 Magnetic Observations Inversion of the ground magnetic data was not successful in modeling the keel, even though magnetic susceptibility measurements made on drill-core show the keel to have high susceptible values. Forward modeling of magnetic data is performed to determine if the field ground magnetic data that were inverted in Section 4.3 could detect the keel. In addition, down-hole synthetic measurements are made in an attempt to increase signal levels of the magnetic response of the keel, and determine how near such down-hole measurements need be in order to detect the keel. Forward modeling is conducted on a magnetic susceptibility rock-model. As outlined in Appendix B, the susceptibility rock-model is constructed using interpreted geologic sections and assigned magnetic susceptibility values based on core measurements3. A second susceptibility model is constructed which is the same as the one just described, 3Even though core susceptibility measurements show that the distribution susceptible material.isn't very predictable within the deposit, it does seem to be consistently high within the keel, especially along the fault. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 143 except all susceptibility values below 425m are replaced with a value of zero. The differ ence between the data collected on each of these models will be used to determine the magnitude of the response of the keel. For the field ground magnetic data a noise level of 2nT was used during the inversions. This resulted in a model that corresponded with the deposit yet was devoid of overly complex structures. Based on this information, a signal level of about 5nT is expected from sources in order for them to be adequately discerned from the noise. 6.3.1 Simulating the ground magnetic field data Five hundred and eighty-five surface measurements of the total magnetic field are calcu lated over both susceptibility rock-models. The inducing field strength and orientation was approximately that in which the original field data were collected. The synthetic data are shown in Figure 6.7 along with a plot of the difference between the two. The response from the keel (as depicted by the plot of the difference between the two data sets) ranges from about -1 to 3 nT in magnitude. This signal level is above a noise level of 2 nT, but only exhibits a signal-to-noise ratio of 1.5; smaller than that required for confident detection. 6.3.2 Synthetic downhole magnetic observations So how can magnetic field observations be used to detect and model the keel? The observations need to contain information about the keel above the expected noise levels. As the source of the magnetic field cannot be moved nearer the receiver, the receiver is moved nearer the source; below the surface. Synthetic magnetic field observations are made below the surface in four drill-holes that extend to a depth of 650m. The drill-holes, shown on a cross-section of the susceptibility rock-model in Figure 6.8, are spaced 50m apart and are located at -1800m west, and from -350m north (against the northern edge Chapter 6. Forward Modeling and Survey Design: Detecting the keel 144 585 data, I • 50.6, D = -13.4 III tiii HT •2500 -2183 -1867 .1550 -1233 -937 -600 Easting (m) (a) Synthetic magnetic data produced from sus ceptibility rock-model. 585 data, I • 50.6, D • -13.4 750 .................x^-r................... ==:±:= | lioo. . . . ... t i „ , -t.. ..•«*....... •_ i i t i i i i nT •2500 -2183 -1867 • 1560 • 1233 -917 -600 East Ins (m) (b) Synthetic magnetic data produced from sus ceptibility rock-model without keel. (c) Difference between synthetic magnetic data with keel and synthetic magnetic data without keel. Figure 6.7: Synthetic magnetic data produced from a susceptibility rock-model of San Nicolas. Data are produced from models with and without the keel. The difference is also plotted. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 145 of the keel) to -200m north. Sixty-six measurements are down the hole to a depth of 650m. The measurements are made every 10m. A B C D Figure 6.8: West-facing cross-section of the magnetic susceptibility rock-model with the keel shown to the lower left. Four drill-holes are also shown, within which synthetic total magnetic field measurements are performed. The total field magnetic response for observations made down-hole on both suscepti bility models are shown in Figure 6.9(a). The response from the model with the keel is shown as a dashed line, while the response from the model without the keel is shown as a solid line. The difference between the two responses are plotted in Figure 6.9(b). The signal that emanates from the keel is easily observed in data from drill-hole A at a depth of between 450 and 500m. The observations made down drill-holes at more northern locations show the signal from the keel decreasing to a magnitude of about +/- 15nT at a distance of 150m from the northern edge of the keel. According to the noise and signal-to-noise levels prescribed above, this response should be sufficient to detect the keel. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 146 Total magnetic field (nTJ I I I...I-1 1 1 I I Lil I I I I Lil I i I I Ul I i A B C D (a) Down-hole, synthetic total magnetic field data. For each hole the response from the model containing the keel is shown by a dashed line, and the response from the model without the keel is shown by a solid line. Total magnetic field [nT] f U i » i i i S ! s . i J S is = a 8 i i . • • A B C D (b) Difference between down-hole data collected on a model with the keel and those collected on a model without the keel. This is equal to the response from the keel alone. Note the change of scale for plots from drill-holes C and D. Figure 6.9: Synthetic total-field, down-hole, magnetic data produced from a magnetic susceptibility rock-model of San Nicolas. Data are produced from models with and without the keel. The difference is also plotted. Drill-hole A is just north of the keel and drill-holes B,C, and D are each 50m further north that the previous one so that drill-hole D is 150m north of the keel. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 147 6.3.3 Summary of magnetic forward modeling Forward modeling of total field magnetic data over a susceptibility rock-model of San Nicolas has produced results that suggest noise levels are too high for surface data to detect the keel. This result is consistent with the inversion results that did not model the keel. Also, drill-hole data collected at a distance of 200m from the keel will still contain information about the keel that is sufficiently above assumed noise levels. 6.4 CSAMT/CSEM observations One-dimensional inverse modeling of CSAMT data collected over San Nicolas has, when the models are 'stitched' together, successfully produced a conductivity structure of the subsurface that defines the ore-body and the Tertiary overburden. A conductivity distri bution that models the keel, however, was not recovered. Through the forward modeling presented in Section 4.4.1 it has already been shown that 3D modeling is probably nec essary in order to recover all the information contained in the CSAMT data, especially information about structure at depth encoded in lower frequency data (0.5 Hz to 32 Hz). But, even if this data were to be modeled in three dimensions, would it contain information about the keel that is above expected noise levels? In order to address this question, and examine ways of increasing signal-to-noise levels of the response from the keel, forward modeling of CSAMT data with different transmitter-receiver layouts is performed over a conductivity rock-model of San Nicolas. Although the field CSAMT data were collected at fifteen frequencies ranging from 0.5 Hz to 8192Hz, forward modeling is only presented for a frequency of 16Hz. This frequency is chosen because the plane-wave skin-depth for this frequency, within a 50 Ohm-m half-space (a reasonable approximation for the conductive environment at San Nicolas), is about 900m. Therefore, signals at this frequency would probably contain sufficient Chapter 6. Forward Modeling and Survey Design: Detecting the keel 148 information about the keel. Of course, as the geometry of the survey changes, the plane-wave approximation will become less valid, but an investigation into the fields generated at 16Hz may give some insight into the problem of detecting the keel using the CSAMT method. When inverting the field data, standard deviations of between 2 and 5 percent were assigned to the apparent resistivity data, and between 1 and 2 degrees were assigned to the phase. For this sensitivity study of the keel, apparent resistivity and phases are not calculated, rather the data remain in terms of magnitudes of electric and magnetic field components. For the purpose of drawing conclusions here, it is assumed that noise levels of greater than 5 percent are necessary for a coherent signal from the keel to be useful in detecting the keel. 6.4.1 Simulated CSAMT observations In a similar manner as before, rock-models were made for conductivity structures with and without the keel. The rock-model was embedded in a half-space of 150 Ohm-m material and a perspective view of the conductivity rock-model (with keel), cut at -400m south to reveal the deposit, is shown in Figure 6.10. The location of the transmitter (grounded at each end) as used in the field CSAMT survey is also shown. Receiver locations were centered over the deposit along lines -600m, -400m, and -200m north. Using EH3D (Appendix E.8), electric and magnetic fields were modeled for the conductivity models due to a current source of 1 Amp, oscillating at 16Hz. As with the field survey, real and imaginary parts of the x-component electric field (.E^and y-component magnetic field (Hy) were measured (an orthogonal coordinate system is orientated with the local coordinate system). These measurements were made at 264 stations located at the surface on a regular, 50m x 50m grid between -1000m and -2150m east, and -150m and -650m north. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 149 Figure 6.10: A perspective view of conductivity rock-model (including keel), cut to show the conductivity structure at -400m south. The approximate resistivity values of the model are shown, as is the location of the transmitter that was used in the field CSAMT survey. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 150 ox e With keel | © Without keel Ex Re [V/m] Ex Im [V/m] Hy Re [A/m] Hy Im [A/m] Easting -2 6 -2.5 -2 .7 -1.5 -2 -1.5 .5 -7 -6 6 x1° x10 x10 x10 Difference 4.9 "2 "1 8 0 -1 -0.5 0 x10 x10 -8 X 10 Figure 6.11: Color plot of synthetic real and imaginary EX and HY fields, at the surface over the San Nicolas deposit, that emanate from a far, 16Hz, grounded transmitter. The first and second rows show the field components that would be produced over conductivity models with, and without the keel, respectively (these are plotted on the same scales). The third row of plots show the difference between fields produced with and without the keel for each component. Each plot covers the area over the deposit, as described previously, although the northing and easting values are not shown. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 151 Plan-view, color plots of real (Re) and imaginary (Im) parts of Ex and Hy, for all stations, are shown in Figure 6.11. The first row of plots shows the fields that would be produced over a conductivity structure of the deposit that includes the keel, the second row of plots shows those that would be produced if the keel were absent, and the third row of plots displays the difference. The difference the keel makes to the real part of the Ex field ranges from -5xl0-9 to 10xl0~9. When compared to the actual measured values, we can see this is less than one percent of the total signal. This signal level will be most likely saturated by noise, and therefore not useful for keel detection. Similar observations are made for the remaining field component measurements, where the effect of the keel is about two percent of the imaginary part of the Ex field, and less than one percent of both the real and imaginary parts of the Hy component. This forward modeling suggests that, for a grounded current source oscillating at 16Hz, noise levels much smaller than five percent would be needed for the keel to be modeled even if the information contained in the Ex and Hy field components can be extracted through 3D inverse modeling. The analysis for other frequencies, or for other E and H components has not been carried out, so it is difficult to draw general conclusions. Furthermore, only relative signal magnitudes and noise levels have been considered, however, absolute signal-to-noise levels are also important — if more current were to enter the ground the magnitudes of the electric and magnetic fields would be larger and, assuming that noise levels remain the same, and a coherent response from the keel above the noise levels could be measured. It will probably be difficult to detect the keel using a plane wave source, that is a transmitter a large distance away from the target. Greater signal is needed for detection and this can come about through a change in survey design. The next section looks at bringing the source in closer to the target in order to increase the signal. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 152 6.4.2 CSAMT/CSEM survey design Another way to increase the signal from the keel would be to change the design of the survey. In contrast to potential field methods, whose sources are fixed, active-source methods, such as CSAMT, allow for the repositioning of sources as well as receivers4. In addition, many more components of the resulting fields can be measured, including real and imaginary parts of x-, y-, and z-components of electric and magnetic fields, and different source frequencies can also be employed. Accordingly, there are many more variables to consider when determining survey design, and optimization methods can be used to design the appropriate surveys given a specific target. Glenn and Ward (1976), Boerner and West (1989), and Maurer and Boerner (1998) address optimization of elec tromagnetic survey designs. For this study, however, only simple survey-design changes are considered, for a single frequency of 16Hz, in an attempt to increase the response of the keel by increasing the current-density within the keel5. The forward modeling is performed in order to demonstrate the types of surveys that could be conducted in order to detect the keel, rather than conclusively show the best survey design for the job. By bringing the current source, or transmitter, closer to the deposit, it is expected that more current will enter the deposit, and therefore larger signals can be measured. Consequently, the signal from the keel should also be larger and the absolute response from the keel may become higher than the noise levels. Figure 6.12 shows the four transmitter locations that are used for the forward modeling. The characteristics or each transmitter setup is described below.: 4When sources are used that generate primary electric and magnetic fields that are not planar in nature at the receiver locations, the survey is considered to be a controlled source electromagnetic survey (CSEM), and no longer a controlled source audio magneto-telluric survey (CSAMT). 5 An increase in current density within the keel is expected to result in larger magnitude secondary fields emanating from the keel. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 153 (a) Far location of transmitter; as used in the field survey. (c) Transmitter with one end below the deposit down a drill-hole. (b) Near location of transmitter. (d) Transmitter located with one end within the deposit down a drill-hole. Figure 6.12: Perspective views of conductivity rock-model showing locations of the trans mitters used for forward modeling of electromagnetic data. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 154 1. Figure 6.12(a) shows the far grounded transmitter location that was used in the field survey. This is 1.6 km long, parallel to the perceived geologic strike, and is 2.5 km north of the center of the deposit (-400m N). 2. For the near transmitter location, the source remains grounded on the surface and is brought closer to the deposit (Figure 6.12(b)). The transmitter is now located along line -200m N (very close to the deposit), and is still 1.6 km long. 3. Figure 6.12(c) shows a very different transmitter location. One end of the trans mitter is located at the surface just north of the deposit (-1600m E, -200m N) and one end is located in a drill-hole below the deposit at a depth of 525m at -1600m E, -400m N. 4. The final transmitter location is shown in 6.12(d). This is also located with one electrode at the surface and one down a drill-hole, but now the electrode is grounded within the deposit at a depth of 325m. A plot of current density magnitudes for models cells on a north-facing cross-section through the deposit is shown in Figure 6.13. An increase of current densities in the keel of more than one order of magnitude results from decreasing the distance to the transmitter. The current density within the keel increases in magnitude to over lxlO-5 A/m2 for the down-hole source locations (not shown), which is an increase of more than two orders of magnitude over the original field survey. The new transmitter locations have increased the currents within the keel as well as the rest of the deposit. Synthetic measurements of Ex and Hy are now made at the surface for all three transmitter locations on conductivity structures with, and without, the keel. The following figures show the measurements made with these transmitter locations and compare them with the results of the simulated field observations (Section 6.4.1). Chapter 6. Forward Modeling and Survey Design: Detecting the keel 155 Eosttng [m] Keel (a) Magnitude of current density through a north-facing cross-section of the rock-model at -400m north for the far transmitter location (2500m north of the deposit) used in the field survey. Easting [m] Keel (b) Magnitude of current density through a north-facing cross-section of the rock-model at -400m north for the near transmitter (at the northern edge of the deposit) location. Figure 6.13: North facing cross-sections through a rock-model, which includes the keel, showing the magnitude of current density in each cell for two, (a) far and (b) near, trans mitter locations. The cross-section is located at -400m north and shows an increase in current density within the keel (the small anomaly to the lower left of each section) from about 5xl0"8 to 6xl0~7 A/m2. Four figures are presented; one for each part (real and imaginary) of the Ex and Hy components. Within each figure there are four columns of plots; each column represents one transmitter location. Each column of plots contain synthetic readings calculated for: (1) a conductivity structure with the keel, (2) a conductivity structure without the keel, and (3) the difference between the two. Each plot is for receivers at the surface and covers an area over the deposit which extends from -2150m to -1000m east, and -150m to -650m north (the northing and eastings are not shown). The plots shown for observations with all of the transmitter locations, with and without the keel, are presented at the same scale (given by a color bar to the left), while the difference plots, for all transmitter locations, are presented at a different color-scale. As shown in Figure 6.14, the real part of the Ex component of the electric field measured at the surface increases by two orders of magnitude when the transmitter Chapter 6. Forward Modeling and Survey Design: Detecting the keel 156 ox c -o Far source Near source DH: source DH: source below deposit within deposit Difference of surface Real E, data lxltr* mm Easting Surface Real E, data with keel "WW IBB Surface Real E, data without keel lxlO9 in 1 xlu"4 -1 V/m I1 xlO6 1-1 Figure 6.14: Real part of Ex surface measurements made for four different transmit ter source locations over conductivity rock-models of the San Nicolas deposit with and without the keel. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 157 location is brought close to the deposit. The difference, between measurements made on conductivity structures with and without the keel, increases by just under two orders of magnitude when the transmitter is brought near to the deposit, whereas it increases by nearly three orders of magnitude when one end of the transmitter is placed near the keel, down a drill-hole. Also, the response from the keel seems larger when the transmitter is placed below the deposit, as opposed to within the deposit. A similar result is shown for the imaginary part of the x-component of the electric field at surface receiver locations. While the observed signal increases by about two orders of magnitude when the distance to the transmitter is decreased, the difference between measurements made with and measurements made without the keel is only significant when one end of the transmitter is located below the surface. Far source DH: source DH: source below deposit within deposit Near source Easting Surface Imaginary E„ data with keel Surface Imaginary E, data without keel Difference of surface Imaginary Ex data 1 xlO1 1-1 V/m l xlO7 ftl.5 Figure 6.15: Imaginary part of EX surface measurements made for four different trans mitter source locations over conductivity rock-models of the San Nicolas deposit with and without the keel. Very similar results are also observed for surface measurements of the magnetic field, Chapter 6. Forward Modeling and Survey Design: Detecting the keel 158 and are summarized below. For the three transmitter locations close to the deposit, the real part of the Hy field increases by about two orders of magnitude, while the imaginary part increases by just under two orders of magnitude. The difference between Hy measurements made on conductivity rock-models with and without the keel increases by two orders of magnitude only when one end of the transmitter is grounded close to the keel, down a drill-hole. 6.4.3 Summary of CSAMT/CSEM forward modeling Forward modeling has shown that current densities, and the resulting electromagnetic fields, all increase when the transmitter is brought closer to the deposit. The response of the keel relative to the total signal observed also increases, but does not exceed one percent of the total signal. However, with these elevated field strengths, the noise levels may remain close to their original level, and the signal-to-noise could be improved. The results presented are for two field components (Ex and Hy) at one frequency (16 Hz). These components are measured in CSAMT surveys because they have maximum coupling with the propagating plane-wave source with the transmitter geometry used, and therefore will produce the strongest signals. As other survey geometries are considered, the source does not approximate a plane wave, and the orientation of maximum signal is harder to determine. Therefore, other field components should be measured so all information is recorded. In addition, the aspect of placing the receiver locations down a drill-hole has not been presented. This would probably result in an increase in signal-to-noise response of the keel depending on the proximity of the drill-hole to the keel. Finally, the response of the keel at different frequencies might also show different results. The preliminary study made on the design of CSAMT surveys, targeted at detecting the keel, clearly demonstrates the scope of this problem. As such, definite conclusions cannot be made regarding the detectability and definition of the keel with CSAMT Chapter 6. Forward Modeling and Survey Design: Detecting the keel 159 techniques, apart from stating that the range of survey designs available probably includes one that could provide data to sufficiently model the keel. 6.5 DC resistivity Inversion of DC (direct-current) resistivity measurements, made using the 'Realsection' array configuration, has not been successful in identifying the San Nicolas deposit. This may be due to current channeling in the conductive overburden, or the small number of transmitter locations (more transmitter locations seem to have a greater constraint on conductivity distributions at depth). Even so, forward modeling is performed in order to quantify how sensitive the survey is to the keel, and to determine if other survey designs might provide more information about the keel. In order to determine if the signal from the keel is above noise levels, the noise levels must be known, however, only estimates of the noise are available. Noise levels of five percent were assigned to the DC data for inverting, and, although the recovered models didn't seem very successful, noise levels between two and eight percent typical for such surveys . For a better understanding of the noise levels, instrument sensitivity and specific field conditions would have to be taken into consideration. 6.5.1 Simulation of Realsection resistivity observations Two conductivity models are used for the forward modeling; one contains the structure of the keel and one does not. A synthetic Realsection6 survey was conducted over both of these models, the results of which are presented as apparent resistivity values in Fig ure 6.19. The survey, which has a layout as described in Appendix C.6, was conducted ^Realsection is the name given to a type of modified gradient array survey that was devised by Quantec Geoscience, and employs several transmitter dipoles of different lengths. When plotted on a section, the resulting picture merely represents the data and not the true electrical structure of the earth. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 160 over the deposit on line -400m north and involves measuring the potential difference between two receiver electrodes. This potential difference is the actual measurement made in the field, but apparent resistivities are commonly calculated using the following equation, [Telford, Geldart, and Sheriff, 1990b]. where A; is a geometric factor that is determined by the position of the transmitter and receiver electrodes, V is the potential difference measured with the receiver dipoles, and / is the transmitted current. The synthetic Realsection data are plotted in Figure 6.16 and show a response from the keel that is about one percent of the total signal. This amplitude is below the expected noise levels (5 %). Another survey design is employed to see if the response from the keel can be in creased. The dipole-dipole survey design is chosen. Dipole-dipole arrays are usually more sensitive to lateral variations in conductivity, and, although the depth penetration may not be as great as the Realsection survey, it is hoped that a greater sensitivity to the keel is observed. Dipole-dipole arrays also employ many more transmitter locations and therefore can provide better constraint of inversion models at depth. 6.5.2 Forward modeling of dipole-dipole resistivity observations. A dipole-dipole array was used to collect synthetic resistivity data over two rock-models of the San Nicolas deposit; one that includes the keel, and one that does not. The dipole-dipole array employed has thirteen 100m transmitter dipole locations, each with five 50m, and five 100m, associated receiver dipoles. For display purposes, resistivity values are plotted at the intersection of two forty-five degree lines that lead down from the center Chapter 6. Forward Modeling and Survey Design: Detecting the keel 161 | 40 Ohm-m I 77.5 I 150 •2275 -2033 •1792 -1550 -1308 Easting (m) -2275 -2033 -1792 -1550 -1308 -1067 Easting (m) (a) Synthetic line of Realsection apparent resis- (b) Synthetic line of Realsection apparent resis tivity data collected over a conductivity rock- tivity data collected over a conductivity rock-model with the keel. model without the keel. -2275 -2033 -1792 -1550 -1308 -1067 -825 Easting (m) (c) Percent difference between synthetic Realsection apparent resistivity data with keel and synthetic apparent resistivity data calculated without keel. Figure 6.16: Synthetic, DC, apparent resistivity, Realsection data produced from conduc tivity rock-models of San Nicolas. The x-locations refer to the midpoint of the receiver dipoles, while the y-axis refers to different transmitter dipole positions; the bottom line corresponds to data collected with the longest transmitter, and the top line corresponds to data collected with the smallest transmitter length. The section is just a representa tion of the data and should not be taken as the true structure of the earth. The percent difference between data produced from models with, and models without the keel is also shown. Note: the apparent resistivity color-bar is on a log scale with low resistivities (high conductivities) displayed in red, while the percent difference is plotted on a linear scale. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 162 of the transmitter (I) and receiver (vx) dipoles, as shown in Figure 6.17. This pseudo-section is a convenient way of plotting the data but does not represent the structure of the earth. The values are interpolated and displayed in Figure 6.18. surface Figure 6.17: Schematic diagram of dipole-dipole array for one transmitter location, show ing the pseudo-section plotting points relative to the electrode geometry. The synthetic dipole-dipole data collected using the two conductivity rock-models show apparent resistivity values that range from 25 to 75 Ohm-m. The signal from the keel, as determined by calculating the difference between the data collected over the two models, has resistivity contrast values ranging from -2 to 0.3 Ohm-m. This response is up to 2.75 % of the magnitude of the total signal, which suggests that this dipole-dipole array configuration is more sensitive to the keel that the Realsection array, however, it is still below the assumed noise level of 5%. 6.5.3 Summary of DC resistivity forward modeling While inverse modeling of DC Realsection resistivity data hasn't resulted in an appro priate conductivity model of San Nicolas, the data have been shown to be sensitive to the presence of a conductive keel. Unfortunately, the response from the keel is below expected noise values and so other survey designs are sought in order to increase the signal-to-noise. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 163 25 I Ohm m •2394 2136 1879 1622 1365 1107 -850 -2394 -2136 1879 1622 1365 1107 -850 Easting (m) Easting (m) (a) Pseudo-section of synthetic, dipole-dipole, (b) Pseudo-section of synthetic, dipole-dipole, resistivity data collected over a rock-model that resistivity data collected over a rock-model with-includes the keel. out the keel. 13 I 2.75 % 1.25 -0.25 •2394 -2136 1879 1622 1365 1107 850 Easting (m) (c) Percent difference between synthetic dipole-dipole resistivity data with keel and synthetic dipole-dipole data without keel. Figure 6.18: Synthetic apparent resistivity dipole-dipole pseudo-sections produced from conductivity rock-models of San Nicolas. Pseudo-sections of data produced over con ductivity rock models with and without the keel are plotted on a log scale along with the percent difference between the two, which is plotted on a linear scale. The plotting convention use is shown in Figure 6.17. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 164 A traditional dipole-dipole survey design was chosen as an improvement over the Realsection array in detecting the keel. In changing from a Realsection array to a dipole-dipole array, the response from the keel has been shown to increase by almost three-fold (from 1% of the total signal to 2.75%). Furthermore, if this data were to be inverted it is expected that the dipole-dipole data would constrain the model better at depth than the Realsection data. Unfortunately, the response from the keel is still below estimated noise levels so other survey designs are required, and down-hole methods may prove necessary in order to receive a coherent signal from the keel. The increase in response from the keel demonstrates an improvement in survey design, but comes at the expense of more field work. Five transmitter and eighty-five receiver dipoles were used in the Realsection survey, verses thirteen transmitter and one hundred and thirty receiver dipoles for the dipole-dipole survey. Optimization of DC resistivity surveys for a given target is a very broad subject due to the many different transmitter-receiver geometries that can be employed. The work presented in this section introduces the subject of DC resistivity survey design with two-dimensional surveys over a three-dimensional object, but investigations into three-dimensional designs would provide more complete results. 6.6 Chargeability Realsection IP data collected at San Nicolas shown in Figure D.15, show apparent charge-ability values up to 24 msec. Physical property measurements suggest the deposit as a likely source of the signal, and inversion of this data, along with gradient array IP data, has shown a distribution of anomalously chargeable material that corresponds well with the deposit (Section 4.5). Although the keel was not modeled through inversion of the IP data, a concentration of chargeable material is observed along the fault-zone that bounds Chapter 6. Forward Modeling and Survey Design: Detecting the keel 165 the deposit to the east. It is possible that these chargeability values persist within the keel and IP surveys could be used to detect the keel. As a result, forward modeling of IP data is used to determine the sensitivity of IP data to the presence of the keel. Synthetic IP data are calculated using DCIPF3D — a forward modeling routine that first calculates the primary voltages over a conductivity distribution, and then calculates the secondary voltages based on the chargeability of the material. An overview of the forward modeling method is given in Appendix C.6. A chargeability rock-model is created from digitizing geologic cross-sections and as signing values to the rock units based on physical property values. A maximum value of seventy-five milli-Volts per Volt7 is assigned to the massive-sulphide. As before, syn thetic observations are made over two chargeability models; one with the keel and one without. The difference between the two sets of observations will indicate the sensitivity of the survey technique to the presence of the keel. Two survey designs are examined for sensitivity to the keel: the Realsection survey design that was used in the field, and the dipole-dipole survey design that was used in Section 6.5.2. The success of a survey in providing information about the keel depends on the strength of the signal from the keel and the noise levels encountered during the survey. Standard deviations of about 3 percent of the apparent chargeability plus a base level of 1 mV/V were assumed for the inversion of apparent chargeability data, and resulted in a reasonable model with some structure, but without too much noise. Therefore, a signal from the keel would have to be above this level in order for it to add information to the inversion model. 7Here, the apparent chargeability is calculated by taking the ratio of secondary voltages to primary voltages. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 166 6.6.1 Simulating Realsection IP field observations. Synthetic apparent chargeability data are first calculated for the survey layout that was used in the field. I mV/V 6 mWV 2.5 -2275 -2033 -1792 -1550 -1308 -1067 -825 -2275 -2033 -1792 -1550 -1308 -1067 -825 Eastmg (m) Easting (m) (a) Synthetic line of Realsection apparent (b) Synthetic line of Realsection apparent chargeability data collected over chargeability chargeability data collected over chargeability rock-model with the keel. rock-model without the keel. (c) Difference between synthetic realsection apparent chargeability data from a chargeability rock-model including the keel and synthetic apparent chargeability data from a chargeability rock-model without the keel. Figure 6.19: Synthetic apparent chargeability Realsection data produced from charge-ability rock-models of San Nicolas. The x-locations refer to the midpoint of the receiving dipoles, while the y-axis refers to different transmitter dipole positions; the bottom line of data being collected with the longest transmitter, and the top line of data being col lected with the smallest transmitter length. The difference between the two data sets is also plotted. The synthetic Real-section array apparent chargeability data are displayed in Fig ure 6.19. The difference between data collected on models with and without the keel range between -0.15 and 0.3 mV/V. The difference is below the prescribed noise level and so, if the noise levels are appropriate, this survey design would not be useful in identifying the keel. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 167 In an attempt to increase the chargeability signal from the keel, a dipole-dipole array configuration is examined using forward modeling methods. The results from this study are presented below. 6.6.2 Synthetic dipole-dipole chargeability observations. Synthetic dipole-dipole IP surveys, conducted on two chargeability models (one with the keel and one without), employ the same array configuration used in the DC resistivity survey described in Section 6.5.2. The data, displayed as apparent chargeabilities in Figure 6.20, show values that range from -0.5 to 6 mV/V. The response of the keel is shown in Figure 6.20(c) as the difference between data calculated using a chargeability model that includes the keel, and one that does not. This difference has values that range from -0.5 to 1.1 mV/V. These data do not show a response from the keel that is above the noise level, however, they are larger than the response from the keel that was calculated for the Realsection IP data. 6.6.3 Summary of chargeability forward modeling Forward modeling of apparent chargeability data over a chargeability rock-model, which is thought to represent that present at San Nicolas, shows that it is unlikely Realsection IP data contain information about the keel that is above the expected noise levels. Therefore, insight into the presence of the keel is not expected from inversion of this data. Forward modeling of dipole-dipole chargeability data shows that an increase in re sponse from the keel can be achieved, however, the amplitude is still below the expected noise levels. If the chargeability values have been appropriately assigned to the rock-model and the expected noise levels are also approximately correct, then detecting the keel using IP methods may have to involve reducing the distance from the keel to either Chapter 6. Forward Modeling and Survey Design: Detecting the keel 168 •2394 -2136 1879 1622 1365 1107 850 Easting (m) (a) Pseudo-section of synthetic dipole-dipole ap parent chargeability data collected over a charge-ability rock-model with the keel. •2394 -2136 1879 1622 -1365 -1107 850 Easting (m) (b) Pseudo-section of synthetic dipole-dipole ap parent chargeability data collected over a charge-ability rock-model without the keel. °°. °\ **, \ "'. "", \ °\ *•, \ °% *', I mV/V 0.3 1-0.5 2394 2136 1879 1622 1365 1107 -850 Easting (m) (c) Difference between synthetic dipole-dipole apparent chargeability data from a chargeability rock-model including the keel and synthetic apparent chargeability data from a chargeability rock-model without the keel. Figure 6.20: Synthetic apparent chargeability dipole-dipole data produced from charge-ability rock-models of San Nicolas. The plotting convention used is shown in Figure 6.17. The difference between the two data sets is also plotted. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 169 the transmitter, the receivers, or both. This could be achieved by placing electrodes (transmitter or receiver electrodes) down a drill-hole. It should also be noted that the secondary voltages measured in the dipole-dipole survey are about one tenth the magnitude of those measure in the Realsection survey. This could have implications regarding the signal-to-noise of the response of the keel as lower magnitude observations may be closer to ambient noise and instrument sensitivity levels. The examination of two survey layouts hardly constitutes a comprehensive analysis of how to detect the keel using IP methods, however, it has been shown that different types of surveys have different sensitivity to the keel and, should the rock-model values be correct, down-hole methods may be necessary. 6.7 Conclusions Forward modeling, using physical property rock-models at the deposit scale, has shown that all of the data collected at San Nicolas could be sensitive to the presence of the keel, yet only the gravity data would show a response that is significantly above the noise level. This is consistent with the inversion results presented in Chapter 4, in which the gravity data produced a density distribution of the deposit that had an extension in the vicinity of the keel, whereas the inversion of other data did not. The surveys employed in the field were designed to detect anomalous structure and bodies of unknown location and geometry, and they were successful in finding San Nicolas. However, changes to the survey design are needed when deep, small ore-bodies such as the keel are the target. To this end, different survey designs have been examined in an attempt to increase the observable response of the keel. Insight has been provided about what data should be measured in order to see the keel. Chapter 6. Forward Modeling and Survey Design: Detecting the keel 170 The target of the investigation in a potential-field survey is also the source of the signal. For this reason there is a limit to the signal from the keel that can be measured at the surface. Consequently, down-hole methods may be needed. On the other hand, changing transmitter-receiver geometries on the surface can improve the magnitude of the response from the keel, when active-source methods are used. However, down-hole methods also prove useful for active source measurements. The use of drill-holes in geophysical surveys at San Nicolas is a logical progression in the search for deeper ore as many drill-hole already exist. This ties in with the over-all detailed exploration work that has been the motivation for this chapter. Information from drill-holes can be used to constrain inversion models in the form of a priori information (Chapter 5), or geophysical data can be collected down-hole. I have shown how both of these methods can increase the level of understanding of the subsurface, and therefore, increase the efficiency of a deposit-scale exploration program. The issue of increasing the signal-to-noise of the response from the keel has been addressed by assuming noise levels and trying to increase the signal based on survey design. For a complete study in this area the noise should also be investigated. Different sources of noise, and how absolute noise levels, and noise levels relative to the signal, change in response to different survey designs could also be addressed. As well as considering different survey designs of the methods used to collect data at San Nicolas, other methods may be more successful in defining the keel, and investigation into such methods as MMR (magneto-metric resistivity) could prove fruitful. Chapter 7 Conclusions and Suggestions for Further Work The work completed at San Nicolas has successfully demonstrated how geophysical in version can be used at various stages of an integrated exploration program. This work has been presented within the context of an exploration program that is summarized by Figure 7.1. • Geologic and geophysical exploration models were formulated and provided direction for the subsequent work at different stages. Regional-scale inverse modeling allowed conclusions to be made about physical property distributions of the earth. These models were given a geologic context through interpret ion using physical property measurements. The results provide important information for finding mineralization, and could have been used to direct follow-up work. Local-scale, forward and inverse modeling, along with physical property analysis, has shown how geophysics can be used to define mineralization well enough to successfully spot drill-holes. This leads to further iterations at the local scale that include using drill-hole information directly in an inversion, and modifying surveys in order to help delineation. The work completed has demonstrated how geophysical inversion can be used within general exploration programs by providing information specific to San Nicolas. 171 Chapter 7. Conclusions and Suggestions for Further Work 172 Exploration Goal: •find economic mineralization Geologic Constraints; •existing mapping •historical workings •tectonic setting J Geologic •descriptive model | Models: 'genetic model *1 *=> Exploration Model Regional Geologic Mapping: •Outcrop geology •Structure -Alteration Geochemistry •Stream sediments Remote Sensing •Mu hi spectral •Hyperspcctral Geologic Interpretation Target Generation Terminate Exploration Program Update Geologic Exploration Model Geologic Logging i 4 I Detailed Geologic Mapping: •Stratigraphy •Structure •Alteration Trenching Rock samples •Assays •Alteration mineralogy (X-ray/PIMA) Geochemistry •soil surveys Geologic Interpretation Mineralization Not Intersected Locate Drill-hole Terminate Exploration Program Mineralization Intersected Goal Realized: Mineral Resource Geophysical Exploration Flowchart Geophysical .expcctcd 0ftarget Exploration .physical properties contrasts Model: Regional/Reconnaissance Scale Geophysics: Physical Property Measurements Existing Data Survey Design Data Collection I Gravity AMAG (+/- ARAD) AEM I IP Other VLF, MT Modeling Individual Inverse Joint/Simullaneous Forward Inverse Correlation Conclusions abont physical property distributions of the the earth Update Geophysical Exploration Model Additional Physical Property Measurements Local/Deposit Scale Geophysics: Survey Design Data 1 1 • 1 Collection Gravity OMAQ EM: CSAMT, UTEM DDIP Other: VLF Modeling Individual' Forward Joint/Simultaneous Forward Inverse Correlation Conclusions about physical property distributions of the the earth Physical property core measurements Down-hole methods Delineation Figure 7.1: Geophysical Exploration Methodology Flowchart. Chapter 7. Conclusions and Suggestions for Further Work 173 7.1 Summary of results Regional-scale geophysical programs can cover a large amount of land in a small amount of time. The information that they contain can often be more easily interpreted using geophysical inversion. Individually, the inversions can produce a number of targets for selection. In combination, the number of follow-up targets can be greatly reduced, thus focusing an exploration program faster than before. At San Nicolas, regional-scale modeling has provided information about the extents of a horst, within which the deposit is contained. The location and geometry of deep intrusives has also been modeled. These intrusives might have been the source of mineral-rich hydrothermal fluids and are therefore important for future exploration efforts. In addition, faults have been defined that cross-cut the horst and contain tertiary volcanics. Discriminating between these faults and older faults within the horst may be useful in finding structures that control mineralization. Shallow dense and susceptible mafic intrusives have also been identified. Combining the results from several physical property distribution models has signif icantly reduced the number of targets. Apart from the deposit, other areas of mineral ization were also detected. The potential of airborne EM methods to detect the deposit has also been inves tigated. With both the frequency-domain and time-domain methods only information regarding the upper 200m of the earth has been recovered. Forward modeling of the time-domain data, however, suggest that deeper conductivity information might be re coverable. Local-scale inversions can provide detailed information about the subsurface. This can be used to guide drilling programs. Furthermore, feedback from drilling programs (in the form of physical property measurements on the core, or down-hole surveys) can be used Chapter 7. Conclusions and Suggestions for Further Work 174 to improve physical property inversion models and thus increase the success of drilling. The use of drill-core information in such an iterative manner has been demonstrated. Deposit-scale modeling at San Nicolas has produced physical property models of the subsurface that correspond well with the extents of the deposit. The results confirm the ability of geophysics to find the concealed deposit. The models produced also have geometric features that emulate those of the deposit. Details such as the extension of susceptibility material to the north and the concentration of chargeable material along the fault are also seen at this scale. Including drill-hole information within inversions results in improved models being constructed. The physical property range and the extent of the deposit are much better constrained with the inclusion of one or more drill-holes. Physical property distributions away from the drill-holes also change, especially below the deposit. This allows the possibility of more sulphide below the deposit to be all but discounted. The inclusion of drill-hole information in an inversion addresses the discrepancy be tween physical property measurements made on small drill-core samples and the physical property values recovered in inversion models. This discrepancy comes from inability of a large model cell to adequately represent fine-scale structure. Furthermore, the types of models being constructed have to be consistent with known information, that is smooth models should not be expected when sharp discontinuities have already been found. These issues have been addressed and possible solutions have been put forward, however much more understanding is required. The use of forward modeling to evaluate the sensitivity of survey designs has also been demonstrated. The availability of forward modeling routines allows for quantitative evaluation of survey designs that can lead to more appropriate surveys being conducted in the future. While not addressing the topic of optimizing survey designs, ideas with which survey design can proceed have been introduced. Chapter 7. Conclusions and Suggestions for Further Work 175 Forward modeling has shown that most of the geophysical survey designs employed at San Nicolas were probably inadequate to detect the keel. Simple survey design changes are suggested that, when evaluated using forward modeling techniques, show an increased sensitivity to the keel. 7.2 Further work The work completed at San Nicolas has successfully demonstrated how geophysical in version can be used at various stages of an integrated exploration program. In doing so, many issues have been raised that have not been conclusively resolved. Some of these issues form the basis of suggestions for further work. These are presented below: 1. The chargeable nature of the San Nicolas deposit and the forward modeling analysis that was performed, both suggest that information about the chargeability of the earth could be contained in the Questem, time-domain, airborne, electromagnetic data. If such information could be recovered from the data this would be prove valuable for two reasons. Firstly, the ability to search for chargeable bodies from the air could revolutionize regional exploration programs. Secondly, accounting for the chargeable material will allow for more accurate, and possibly more deeply resolved, conductivity models to be generated. So how should these data be inverted? Cole-Cole parameters along with a conductivity value could be recovered for each layer of a one-dimensional model. However, this would dramatically increase the number of unknowns in the problem and non-unique combinations of the Cole-Cole parameters may be able to produce the same data. Additional information may be required to recover unique solutions. The data could provide additional information if all three components of the magnetic field were to be inverted simultaneously. Furthermore, three-dimensional inversion of these data could provide more information and result Chapter 7. Conclusions and Suggestions for Further Work 176 in much more easily interpreted models. Some of these issues are already being addressed by the UBC-GIF group. 2. More accurate modeling of the deposit by inverting surface CSAMT data might have resulted in less-drilling needed for the mineralization at San Nicolas to be defined. In order to define the deposit more accurately, three-dimensional inversion of CSAMT might prove useful. Forward modeling of CSAMT data suggested that only half of the recovered model accurately reflects the information encoded in the data, so a full 3D inversion might better define the deposit. 3. For external information to be used to maximum benefit within an inversion, meth ods of representing fine-scale physical property features in a large model must be developed. But, how can this information be used to constrain an inversion? Is bounding the range of values a model cell can assume the best method, and how much detail should be accounted for? These are questions that certainly need much more investigation before we can produce models that are truly consistent with the earth. 4. The subject of survey design has only been briefly addressed in this thesis. The advent of widespread inversion modeling allows for data to be collected that don't have to provide information at face-value. We can start to depend in inversions to decode the desired information about the subsurface. Optimization techniques could be used to design surveys specific to the task-at-hand that will be able to adequately resolve the earth where needed. This will hopefully reduce ambiguities in the resulting models. Consequently, survey design is an important topic that could receive increased attention in the future. What does it mean for a survey to be well designed? Does it mean that the data are Chapter 7. Conclusions and Suggestions for Further Work 177 sensitive to the target, or does it mean that the data can adequately constrain an inversion that will model the target? These are questions that need to be addressed at the survey design stage of each exploration program, and the answers will depend on the type of subsurface information being sought. Other obvious questions that come to mind with respect to survey design are the following: (1) How efficient are current survey designs at defining the subsurface? and, (2) How might current survey designs be improved? The answers to these questions will depend on the application of the method, and may be contained within the context of an exploration program, even within this thesis. As better survey design methods are developed, surveys could become much more tailored to the exploration stage at which they are implemented. This would result'in an even more integrated and efficient exploration program. References [Backus and Gilbert, 1968] Backus, G., and Gilbert, F., The resolving power of gross earth data, Geophysical Journal of the Royal Astronomical Society, 16, 169-205, 1968. [Brescianini et al., 1992] Brescianini, R. F., Asten, M. W., and .McLean, N., Geophys ical characteristics of Eloise Cu-Au deposit North-West Queensland, Exploration Geophysics, 33-42, 1992. [Boerner and West, 1989] Boerner, D. E., and West, G. F., A spatial analysis of the electromagnetic sensitivity in a layered earth, Geophysical Journal International, 98, 11-21, 1989. [Danielson, 2000] Danielson, T. J., Age, Paleotectonic setting, and common Pb isotope signature of the San Nicolas volcanogenic massive sulfide deposit, southeastern Za-catecas state, central Mexico, University of British Columbia Masters thesis, 2000. [Farquharson and Oldenburg, 1998] Farquharson, C. G., and Oldenburg, D. W., Nonlin ear inversion using general measures of data misfit and model structure, Geophysical Journal International, 134, 213-227, 1998. [Farquharson and Oldenburg, 2000] Farquharson, C. G. and Oldenburg, D. W., Simulta neous one-dimensional inversion of electromagnetic loop-loop data for both magnetic susceptibility and electrical conductivity, presented at GeoCanada 2000, Calgary, Al berta, 29 May - 2 June 2000. 178 References 179 [Flis et al., 1989] Flis, M. F., Newman, G. A., and Hohmann, G. W., Mineral discrimi nation and removal of inductive coupling with multifrequency IP, Geophysics, 54, p 514-523, 1978. [Glenn and Ward, 1976] Glenn, W. E., and Ward, S. H., Statistical evaluation of electri cal sounding methods. Part 1: Experiment design, Geophysics, 41, 1207-1221, 1976. [Grant and West, 1965] Grant, F. S. and West, G. F., Interpretation Theory in Applied Geophysics, McGraw-Hill, 1965. [Haber, 1997] Haber, E., Numerical strategies for the solution of inverse problems, Uni versity of British Columbia Ph. D. Thesis, 1997. [Hansen, 1998] Hansen, P. C, Rank-deficient and discrete ill-posed problems, Society for Industrial and Applied Mathematics, 1998. [Jansen, 2000] Jansen, J. Diamond drill log: Sal-133, provided by Minera Teck S.A. DE C.V., Zacatecas, Mexico, 2000. [Johnson et al., 2000] Johnson, B. J., Montante-Martinez, J. A., Canela-Barboza, M., and Danielson, T. J., Geology of the San Nicolas deposit, Zacatecas, Mexico in VMS deposits of Latin America, eds. Sherlock, R., and Logan, M. A. V., G.A.C. Mineral Deposits Division special publication, No. 2, 71-86, 2000. [Li and Oldenburg, 1995] Li, Y. and Oldenburg, D. W., 3D inversion of gravity-data, in JACI Annual report, produced by the Geophysical Inversion Facility at UBC, 1995. [Li and Oldenburg, 1998a] Li, Y. and Oldenburg, D. W., Separation of regional and resid ual magnetic field data, Geophysics, 63, 431-439, 1998. References 180 [Li and Oldenburg, 2000b] Li, Y. and Oldenburg, D. W., Incorporating geologic dip in formation into geophysical inversions, Geophysics, 65, vol 1, 148-157, 2000. [Li and Oldenburg, 2000a] Li, Y. and Oldenburg, D. W., 3D inversion of induced polar ization data, Geophysics, 65, 2000. [Li and Oldenburg, 2001] Li, Y. and Oldenburg, D. W., Fast inversion of large-scale mag netic data using wavelet transforms, Submitted to Geophysics awaiting publication, 2001. [Li and Oldenburg, 1998b] Li, Y. and Oldenburg, D. W., 3D inversion of gravity data, Geophysics, 63, 109-119, 1998. [Li and Oldenburg, 1997] Li, Y. and Oldenburg, D. W., 3-D inversion of magnetic data, Geophysics, 61, p 394-408, 1997. [Lydon, 1988] Lydon, J. W., Volcanogenic massive sulfide deposits, part 1: A descriptive model in Ore deposit models, eds. Roberts, R. G., and Sheahan, P. A., Geoscience Canada, Reprint Series 3, 145-154, 1988. [Maurer and Boerner, 1998] Maurer, H., and Boerner, D. E., Optimized and robust ex perimental design: a non-linear application to EM sounding, Geophysical Journal International, 132, 458-468, 1998. [Menke, 1984] Menke, W., Geophysical data analysis: Discrete inverse theory, Academic Press, 1984. [Nabighian, 1979] Nabighian, M. N., Quasi-static transient response of a conducting half-space - An approximate representation, Geophysics, 44, p 1700-1705, 1979. References 181 [Nabighian and Macnae, 1987] Nabighian, M. N. and Macnae, J., C, Time domain elec tromagnetic propsecting methods in Electromagnetic Methods in applied geophysics, volume 2, ed. Nabighian, M., Society of Exploration Geophysics, Tulsa, Oklahoma, p 427-520, 1987. [Nagy, 1966] Nagy, D., The gravitational attraction of a right rectangular prism, Geo physics, 31, 362-371, 1966. [Pelton et al., 1978] Pelton, W. H., Ward, S. H., Hallof, P. G., Sill, W. R., and Nel son, P. H., Mineral discrimination and removal of inductive coupling with multifre-quency IP, Geophysics, 43, p 588-609, 1978. [Quinn et al., 1995] Quinn, J. M., Coleman, R. J., Shiel, D. L., and Nigro, J. M., The joint US/UK 1995 epoch world magnetic model, Technical report No. 314 (TR-314), Naval Oceanographic office, Stennis Space Center, MS, 1995. [Routh and Oldenburg, 1999] Routh, P. S. and Oldenburg, D. W., Inversion of controlled source audio magnetotelluric data for a horizontally layered earth, Geophysics, 64, 1689-1697, 1999. [Routh, 1999] Routh, P. S., Electromagnetic coupling in frequency domian induced po larization data, University of British Columbia Ph. D. Thesis, 1999. [Scales and Tenorio, 2001] Scales, J. A., and Tenorio, L., Prior information and uncer tainty in inverse problems, Geophysics, 66, 389-397, 2001. [Seigel, 1959a] Seigel, H. O., The theory of induced polarization effects (for step-function excitation), in Overvoltage research and geophysical applications, Pergamon Press, 1959. References 182 [Seigel, 1959b] Seigel, H. 0., Mathematical formulation and type curves for induced po larization, Geophysics, 24, 547-565, 1959. [Sherman, 1997] http://mineral.gly.bris.ac.uk/Mineralogy/sulphides/main.html Copyrightl997 David M. Sherman [Telford, Geldart, and Sheriff, 1990a] Telford, W. M., Geldart, L. P., and Sherrif, R. E., Applied Geophysics, Chapter 3, second edition, Cambridge University Press, 1990. [Telford, Geldart, and Sheriff, 1990b] Telford, W. M., Geldart, L. P., and Sherrif, R. E., Applied Geophysics, Chapter 8, second edition, Cambridge University Press, 1990. [Telford, Geldart, and Sheriff, 1990c] Telford, W. M., Geldart, L. P., and Sherrif, R. E., Applied Geophysics, Chapter 9, second edition, Cambridge University Press, 1990. [Walker, 1999] Walker, S. E., Inversion of EM data to recover 1-D conductivity and a geometric survey parameter, UBC masters thesis, 1999. [Ward, 1966] Ward, S. H., Introduction to the search for massive sulphides, Mining Geo physics - Case histories, Society of Exploration Geophysics, 1966. [Ward and Hohmann, 1987] Ward, S. H. and Hohmann, G., W., Electromagnetic theory for geophysical applications in Electromagnetic Methods in applied geophysics, vol ume 1, ed. Nabighian, M., Society of Exploration Geophysics, Tulsa, Oklahoma, p 131-311, 1987. [Weidelt, 1982] Weidelt, P., Response characteristics of coincident loop transient electro magnetic systems, Geophysics, 47, 1325-1330, 1982. References 183 [Zonge and Hughes, 1987] Zonge, K. L. and Hughes, L., J., Controlled source audio frequency magnetotellurics in Electromagnetic Methods in applied geophysics, vol ume 2, ed. Nabighian, M., Society of Exploration Geophysics, Tulsa, Oklahoma, p 713-809, 1987. Appendix A Coordinate System Many of the data at San Nicolas were collected using different coordinate systems. In order to compare and analyze different data sets and models with one another, a standard coordinate system was chosen to which all other data were transformed prior to modeling. The system chosen was a local coordinate system that has a north orientated about 23 degrees to the east of true north. This coordinate system was chosen because the geologic sections and much of the data provided by Teck were already in this coordinate system. Drill-hole fences were also orientated with the local grid. The airborne data was provided in the UTM (Universal Transverse Mercator) Zone 13, NAD (North American Datum) 27 coordinate system. As a result a transformation is presented that was used for converting from UTM (meters) to the local grid system, also in meters. This is accurate to within a few meters at the deposit. XLOC = (cos(-0.36824) * (XVTM - 192432)) + (sin(-0.36824) * (YUTM - 2500845)) - 2000 (A.1) YLOC = -(sin(-0.36824) * (XUTM - 192432)) + (cos(-0.36824) * (YUTM - 2500845)) - 400 (A.2) In the above equation, XLOC is the local x-coordinate, XUTM is the UTM x-coordinate, Yhoc is the local y-coordinate, and YUTM is the UTM y-coordinate. 184 Appendix B Physical Property Analysis B.l Introduction Physical property measurements are essential in order to choose which survey to use and to interpret the results of a survey. This section presents the physical property values that were used in the work carried out at San Nicolas. The physical property rock-models are constructed using these physical property values. These rock-models were used for forward-modeling performed in the thesis. B.2 Density The density values that were used for in this thesis come from density measurements performed by Quantec on a variety of rock-types at San Nicolas, from density measure ments that were made mainly on the ore for mine tonnage estimates, and, where specific rock-type information was insufficient, from physical property values found in Applied Geophysics, [Telford et al., 1990]. Density measurements are summarized below. B.3 Magnetic Susceptibility The magnetic susceptibility (MS) of a volume of rock is a function of the amount of magnetic minerals, (mainly magnetite and pyrrhotite), contained within the rock. These measurements can be interpreted to reflect lithological changes and the presence of al teration zones. 185 Appendix B. Physical Property Analysis 186 Rock-Type Density (g/cm3) Quarternary Alluvium. 1.8 - 2.2 Tertiary Breccia 2.1 - 2.5 Volcanoclastic sediments, tuffs, and flows 2.2 - 2.6 Mafic flows and breccias 2.5 - 2.9 Mafic flows and sediments 2.5 - 3.0 Massive sulphide 3.5 - 4.3 Semi-massive sulphide 3.0 - 3.5 Sulphide stringer zone 2.8 - 3.3 Rhyolite flows and breccias 2.2 - 2.6 Quartz feldspar rhyolite 2.3 - 2.6 Mafic volcanics and siliceous sediments 2.6 - 3.0 Graphitic Mudstone 2.2 - 2.6 Table B.l: Range of densities for rock-types encountered at San Nicolas. B.3.1 Magnetic Susceptibility core measurements Over one thousand magnetic susceptibility measurements were collected on core samples from thirteen drill-holes. These measurements are summerize in Figure B.l for each rock type. Reliability of measurements Susceptibility measurements were collected mostly on three to four meter intervals along the drill-core (except in areas of high susceptibility variation). Tests were conducted on four different sections of core to see how well these coarse measurements sampled the detailed susceptibility features. Results for two test sections of core are presented in the following plots (Figures B.2 and B.3). The graphs show the original coarsely-sampled measurements and measurements made at finer intervals. The results show that the coarse measurements are probably adequately representative of the true susceptibility Appendix B. Physical Property Analysis 187 Tertiary volcaniclastic breccia Upper mafic volcanics Massive sulphides Quartz rhyolite Lower mafic volcanics Mudstone O O o mean •(). 19 mean=l .42 mean=23 „ mcan=(). 12 mean=0.21 mean 0.09 j o o log10scale [xlO3 S.L] 464 o © o o o Figure B.l: The range of susceptibility values for different rock types at San Nicolas. Appendix B. Physical Property Analysis 188 distribution. Figure B.2: Coarse and fine magnetic susceptibility measurements on a section of core from drill-hole SAL 70. B.3.2 Composite models for bounding inversions For magnetic susceptibility measurements that were not evenly sampled, the mean would not be a representative value of several readings as it would be bias towards the areas where more readings were taken. When calculating values to assign large cells from twenty or so unevenly sampled measurements, each measurement was weighted based on the amount or core that it might represent. A depth interval was determined for each measurement that corresponded to the distance between two midpoints. The midpoints were half way between the location of the measurement in question and the locations of the two measurements immediately above and below. Each measurement was weighted based on the length of this interval. It is hoped that this is an improvement over the mean but it is realized that this method is still somewhat ad hoc. The composite model Appendix B. Physical Property Analysis 189 Original measurements GO • r*"i t o '£ * r-l <U o CO 1000 100 10 1 0.1 0.01 o uo CN IT) CO w CO in in oo in m o CO Tf in in n •<r Depth [m] Figure B.3: Coarse and fine magnetic susceptibility measurements on a section of core from drill-hole SAL 34. produced for SAL 25 using this method is shown in Figure B.4 along with the actual measurements. B.4 Resistivity Only twenty-one resistivity measurements were available for different rock-types at San Nicolas. The measurements are summerize below but it should be noted that they may not be representative of the values for each rock-type. B.5 Chargeability As with the resistivity values only twenty-one chargeability values were available. These did, however, sample all of the rock-types in the area. Appendix B. Physical Property Analysis 190 5 •B * <D '—' c 100 10 1 0.1 0.01 o_ % • | o pr81 » « « 0 e a o » • «r— e «a < o 8 100 200 Depth [m] 300 400 500 • Core Measurements—Core Composite Figure B.4: Magnetic susceptibility measurements made on core from drill-hole SAL 25. A composite model was produced that was used to constrain inversion of ground magnetic data. Rock-Type Resistivity (Ohm-m) Quarternary Alluvium 30? Tertiary Breccia 15-30 Volcanoclastic sediments, tuffs, and flows 80 Mafic flows and breccias 80 Mafic flows and sediments 80 Massive sulphide 20 Semi-massive sulphide 20 Sulphide stringer zone 30 Rhyolite flows and breccias 100 Quartz feldspar rhyolite 100 Mafic volcanics and siliceous sediments 80 Graphitic Mudstone 100-200 Table B.2: Range of resistivities for rock-types encountered at San Nicolas. Appendix B. Physical Property Analysis 191 Rock-Type Chargeability (msec) Quarternary Alluvium 5 - 10 Tertiary Breccia 10 - 30 Volcanoclastic sediments, tuffs, and flows 10 - 100 Mafic flows and breccias 30 - 80 Mafic flows and sediments 30 - 50 Massive sulphide 500 Semi-massive sulphide 400 Sulphide stringer zone 500 Rhyolite flows and breccias 10 • - 20 Quartz feldspar rhyolite 100 Mafic volcanics and siliceous sediments 50 -- 80 Graphitic Mudstone 30 -70 Table B.3: Range of chargeabilities for rock-types encountered at San Nicolas. B.6 Physical property rock-models A geologic rock-model was constructed by digitizing five geologic sections (four east-west sections and one north-south section). The sections were digitized on a 25m x25m grid and each grid cell was assigned a rock-type based on the predominant content of each cell. In order to construct the physical property rock-models the digitized geology was assigned values based of the physical property data base presented above. Then an routine was used to interpolate between sections and construct a three-dimensional model. The interpolation routine used triangulation between cells at the same depth. For density and chargeability values a linear interpolation was used, for magnetic susceptibility and resistivity the Logio of the value was interpolated. The resulting models were used for the keel sensitivity and survey design analysis presented in Chapter 6. Figure B.5 shows the conductivity physical property rock-model that was used to forward model Questem, CSAMT, and DC resistivity data. Appendix B. Physical Property Analysis 192 Figure B.5: Perspective view of the physical property conductivity model used for test forward modeling in Section 4.4.1. The model has been sliced at -400m N to expose a north-facing cross-section. Appendix C Geophysical methods Many different geophysical methods were used in the collection of data at San Nicolas. This section gives a brief overview of each method and the underlying physical principles, and also defines the data collected in each case. Cl The gravity method Gravity surveys involve measuring local irregularities in the earth's gravitational field with the aim of using these measurements to determine subsurface density variations. Newton's law of gravitation describes the force of attraction that acts at a distance between two masses. In the example shown in Figure Cl, the force exerted by mass m2 on m\ [Fix] is expressed as: Pn = (Cl) \r2 - riY where G is the universal gravitational constant (6.672 x 10~uNm2/kg2), mi to m2 are masses, f is the unit vector from mi to m2, and fi and f2 are displacement vectors from the origin of mi to m2 respectively. Furthermore the acceleration, or force per unit mass, due to m2 at a point, defined by fi can be stated as: m) = (c.2) 193 Appendix C. Geophysical methods 194 Figure Cl: The force exerted on mass m\ by mass m2 is in the direction of f. By rewriting f as j^f^p-, equation C.2 becomes: 7m2 (r2 - rp |f2 - fJ3 (C.3) Due to superposition, for a number of point masses the gravitational acceleration can be written as: |fi-fi|3 9(ri) = -7 (C.4) and, for a distributed density function pfa), equation C.4 can be stated as the fol lowing volume integral: = f p(F2)(r2 -Wo Jv \r2-n\s Gravimeters use a mass on a zero-length spring to measure the vertical component of the gravitational field as stated below. The usual unit of measurement of the gravitational field is mGal, and most ground gravimeters are accurate to within about 0.01 mGal. Appendix C. Geophysical methods 195 (fi) = -7 / Jv p(r2)(z2 - zi)dv \f2 - ri|3 (C.6) Through the removal of regional fields (the field due to a background density distrib ution) we can focus on the gravitational field due to density variations within our depth of exploration. These data are typically collected at the surface of the earth on a 2D grid, and GPS is often used to accurately locate each gravity measurement to within a few centimeters, thus providing information about ri and z\. The forward modeling operator used in the inversion procedure comes from the evaluation of Equation C.7 through discretization of the density distribution. C.2 The magnetic method Local variations in the earth's magnetic field are measured at the surface in order to gain information about subsurface magnetic susceptibility distributions. When a magnetizable body is placed in an external magnetic field it is magnetized by induc tion [Telford, Geldart, and Sheriff, 1990a]. The amount that the body is magnetized is measured by the magnetization vector J, which is the dipole moment per unit volume. For this work we neglect the possible effects of remanent magnetization and demagneti zation. We assume relatively small inducing magnetic fields, in which case the magne tization vector is proportional to, and is in the same direction as, the inducing field H, thus: (C.7) Appendix C. Geophysical methods 196 J = KH, (C.8) where H = B/u-o (no being the magnetic permeability of free space), and K is the degree to which a body is magnetized by the inducing field or magnetic susceptibility. Most of the minerals which exhibit this type of magnetism are described as ferrimagnetic. These minerals have positive and relatively large susceptibilities. In the magnetic experiment the inducing field is the earth's main magnetic field (Ho), which is produced mainly by processes occuring in the outer core. The intensity of the total magnetic field resulting from the superposition of the main field and the secondary field is measured during magnetic surveys (Equation C.2). B = B0 + Bc (C.9) Through processing, the main field is removed and the secondary, or anomalous field is derived. The anomalous magnetic field due to a distribution of magnetizable materials can be stated as [Telford, Geldart, and Sheriff, 1990a]: Ss(r0) = ^ f (CIO) 4TT JV \r0-r\ It is this field (measured in Teslas, [T]) that is inverted in order to obtain a subsurface susceptibility distribution, K(T). Data are usually presented in nano-Tesla (nT). The strength and orientation of the inducing field (Ho) must be known for inversion of this data. This is obtained by using the program GEOMAG, written by John M. Appendix C. Geophysical methods 197 Quinn of the USGS. GEOMAG used numerical models of the earth's magnetic field to calculate inclination, declination, total intensity, and their associated rates of change given a location and a date [Quinn et al., 1995]. C.3 Inductive electromagnetic methods Inductive electromagnetic methods actively energize the ground using a time varying current source. The current source produces a magnetic field that interacts with con ductive material in the earth, causing electric fields and inducing currents in the ground. The currents in the earth also produce a magnetic field. The magnetic fields emanating from the current in the transmitter are termed primary fields, while the magnetic fields emanating from the current in the ground are termed secondary fields. Both of these magnetic fields induce a voltage in the conductive receiver loop but it is the secondary fields that provide information about the conductivity of the subsurface. C.3.1 Frequency-domain electromagnetic airborne surveys. The Dighem frequency-domain airborne system produces an alternating current wave form in the transmitter loop and operates at several different frequencies. The sources and receivers are small loops that are towed behind the aeroplane (away from electro magnetic interference). In this system the loops are orientated in two configurations: horizontal co-planar, and co-axial, as shown in Figure C.2. The position of the receiver relative to the transmitter is kept constant throughout the survey. Ward and Hohmann [Ward and Hohmann, 1987], derive an expression for the voltage generated in the receiving loop when a sinusoidal current waveform is produced in the transmitter loop and the system is in the presence of a conductive body: Appendix C. Geophysical methods 198 Figure C.2: Orientation of transmitter-receiver (Tx-Rx) pairs in the Dighem fre quency-domain airborne electromagnetic system. V = -iunQNA(H?{s) + #f (w, a, h, s)) (C.ll) where cu is the frequency, no is the magnetic permeability of free space, iV is the num ber of turns in the receiver coil, A is the area of the receiver coil, H^(s) is the vertical component of the primary magnetic field at the center of the receiver coil (and is a func tion of the separation, s, between the transmitter and receiver coils), and H^(u>, cr, h, s) is the vertical component of the secondary magnetic field at the center of the receiver coil, which is a function of frequency, conductivity (tv), height above the ground (h), and transmitter-receiver separation. Normalization of this voltage by the transmitter current, and the representation of the signal as a ratio of the voltage generated in the presence of a conducting body to that generated in free space, allows the measured value to be written as: #f(LJ, LT, h, s) (C.12) Appendix C. Geophysical methods 199 In order to remove the effect of the primary field a bucking coil is used. This is an additional coil that measures the primary field at the receiver and allows the secondary field to be calculated. Because the secondary magnetic field is not in phase with the primary magnetic field, equation C.12 can be represented as inphase and quadrature components, and are usually given in parts per million of the primary field. C.3.2 Time-domain electromagnetic airborne surveys In the time-domain method, a current is passed through the transmitter loop and abruptly interrupted. According to Faraday's law, this change in primary field will induce cur rents in nearby conductors. As time progresses the currents diffuse from the surface of the conductors inwards [Nabighian and Macnae, 1987] and, for a conducting half space, downward and outward in a smoke-ring effect [Nabighian, 1979]. Currents will persist and propagate slower in conductors, and diminish and propagate faster in resistors. The secondary magnetic fields emanating from the induced currents generate an electric field and induce a voltage in a conductive receiver loop. These voltages are measured as a function of time. As a result, a slower decay the voltage will indicate the presence of a conductor. This is the classic use of time-domain electromagnetic systems. Measurements are made by integrating the voltage over time-windows following the turn-off of the current source. The absence of a primary field during recording time (current off-time) allows for small signals to be measured, although some systems also measure during the on-time. The data can be presented as: dB/dt for different times, the secondary, magnetic field as a function of time normalized by the transmitter current, or as voltages induced in the receiver loops. The transmitter is six turns of wire making a 186 square-meter loop, which is slung around an aeroplane (Figure C.3). When the peak current of about 500 Amps is reached the transmitter has a magnetic moment of half a million Amps per square-meter. The Appendix C. Geophysical methods 200 Figure C.3: Diagrammatic sketch of the Questem time-domain electromagnetic system. The transmitter (Tx) is placed around the aeroplane while the three-component receiver (Rx) is towed in a bird behind the plane. GPS positioning is used to measure the position of the receiver relative to the bird. current waveform refers to the behavior of the current in the transmitter as a function of time. The Questem system uses a periodic half-sine waveform alternating in polarity to produce the primary magnetic field (shown in Figure C.4). The current is turned on for 4 msec and recordings are made by the three-component receiver over 15 subsequent time-windows during the 16 msec following the termination of the current pulse. The time-windows increase in width over the recording period so as to improve signal-to-noise in later times. C.4 CSAMT surveys The controlled source audio-frequency magnetotelluric (CSAMT) method is an electro magnetic technique that uses a grounded horizontal electric dipole source to directly inject currents into the ground. These currents produce magnetic fields that induce elec tric fields in the ground as the signal propagates from the source. It is assumed that the fields are measured far from the source (far-field1) in which case the incident fields should :The far-field zone is defined as the region in which the source-receiver separation is much greater than the skin depth. The skin depth is approximately stated as 6 = 503 [m], where p is the resistivity, and / is the frequency. Appendix C. Geophysical methods 201 Figure C.4: The current in the transmitter as a function of time for 50% of one cycle. Receiver measurements are made over 15 time-windows during the 16 msec following the current on-time. Appendix C. Geophysical methods 202 be planar in nature. Components of the electric and magnetic field are measured at a number of frequencies in the audio range (0.1 Hz to 10kHz) (Figure C.5). Perpendicular, horizontal electric and magnetic field values are used to calculate apparent resistivity and phase change at different frequencies. Zonge and Hughes, [Zonge and Hughes, 1987], pro vide a detailed description of the CSAMT method and derive an expression for far-field apparent resistivity (Equation C.13). Pa 2nu.f Ex Hi, [Sim] (C.13) The phase difference between the electric and magnetic fields in Equation C.13 is also calculated and is measured in milli-radians (mrads). Figure C.5: Diagrammatic sketch of the CSAMT survey layout (modified from Zonge and Hughes). Appendix C. Geophysical methods 203 C.5 DC resistivity methods Resistivity methods involve injecting a current into the ground using grounded electrodes, and measuring the voltage between electrodes in the current flow (in Volts). The current, the voltage, and the geometry of the survey are measured so apparent resistivity values can be calculated. A B v, v2 v3 v4 v5 v6 v7 v8 M w f 1 1 \ -X-Figure C.6: "Real-section" array geometry used for collection of DC resistivity and IP data in Section 4.5. For the array configuration show in Figure C.6 the apparent resistivity calculated from measurement V2 would be: Pa = 2nVo I [(rl ~~ r?) ~ (r3 ~~ ri)] where V2 is the measured potential, / is the current applied using electrodes A and B, and ri through r4 are the distances shown on the diagram. In formulating the forward modeling, we start with Ohm's law: J = 1E P (CU) The electric field, E, is the gradient of the scalar potential V so we can write: Appendix C. Geophysical methods 204 E = -W (C.15) With the above equation, and assuming that no charge accumulates (V • J = 0), we can write an equation that relates the resistivity, the potential and the current sources for a homogeneous medium: This equation is the basis for the DC resistivity forward modeling that is used in the inverse modeling —DCINV3D (Section E.7). DCINV3D constructs a resistivity distrib ution that will fit the potential data to a specified degree, given the source locations. C.6 IP methods Induced polarization occurs when a current is applied to the earth and there is an ac cumulation of charged particles in the pore fluid due to either the presence of metallic minerals, clay minerals or graphite, or restrictions in the pore itself. The ability for a material to accumulate these charges is summarized by its chargeability. As a current is applied, the voltage measured across a pair of electrodes increases very fast initially (Vp in Figure C.7) and then more slowly (logarithmically) until a maximum voltage (VT) is reached at a finite time [Telford, Geldart, and Sheriff, 1990c]2. This slow increase in overvoltage corresponds to a build up of charges either side of a pore restriction, much the same way a capacitor charges up. This is known as electrode polarization and occurs in the presence of metallic minerals. 2When the applied current is interrupted, a slow decay is also observed after an initial rapid voltage reduction. V • (-W) = / (C.16) P Appendix C. Geophysical methods 205 A VP t Figure C.7: The voltage response as a function of time, V(t), in the presence of polarizable material, when a current is suddenly applied and then interrupted. Vr = Vp + Vs-Measurements can be made in the time-domain or frequency domain. Time-domain measurements made at San Nicolas involved calculating the area under the decay curve between two times (ti and t2 in Figure C.7) and normalizing it by the total voltage (Vr). In this case chargeability is defined in Equation C.17, and is measured in units of time (usually milli-seconds). The data are converted to secondary voltages (Vs) and forward modeling is carried out using the DC resistivity forward modeling (Equation C.16), but with the conductivity (a) replaced with CT(1 — 77) [Seigel, 1959b]. Where 77 is the chargeability3. This produces large secondary voltage for high chargeabilities and small secondary voltages when 77 is small. 3The chargeability, 77 for this formulation is dimensionless and less than one. The apparent charge-ability data rja can be the ratio of Vs to Vp (also dimensionless and usually small). Alternatively M in Equation C.17 could be used (in units of seconds), in which case the modeled values of 77 are also shown in units of time. ti V(t)dt (C.17) Appendix D Data Data used for the inversions are presented in this section. The survey specifications associated with each data set are given in Table D.l. When inversions are performed the predicted data does resemble the observed data to the specified requirement, and often look identical upon first inspection. Data misfits are presented in Table E.l, however, the predicted data are not displayed themselves. 206 Appendix D. Data 207 10000 -15000 -10000 -5000 0 5000 Easting (m) 2220 10000 15000 m 2389 Figure D.l: Topographic elevations at San Nicolas. Apart from a hill at 6000m east, -1000m north, most of the area exhibits small topographic elevation changes. Outlines of outcropping geology are also shown, along with a hammer and pick symbol that locates San Nicolas at about -1700m east, -400m north. Appendix D. Data 208 15000 10000 -10000 -15000 -10000 -5000 0 5000 Easting (m) 10000 15000 -u mGal ll Figure D.2: Ground gravity observations. The data are the original data supplied by Quantec but with observations on the hill removed. The topographic map shown in Figure D.l indicates a hill with over 100m of elevation gain. Initial inversions produced a large negative anomaly coincident with the hill. Observations made on and close to the hill were subsequently discarded, although the affects of the hill probably persist in some of the remaining observations in this area. A constant was also removed from the data in order to remove a first order regional field. An outline of outcrop geology is overlain. Appendix D. Data 209 -15000 -10000 -5000 0 5000 Easting (m) 10000 15000 -50 nT 50 Figure D.3: Dighem airborne total-field magnetic data, subtle anomaly at -1600m east, -400m north. The deposit is identified by a Appendix D. Data 210 Figure D.4: Subset of Dighem airborne total-field magnetic data used for inversion. Figure D.5: Falcon airborne total field magnetic data. Appendix D. Data 211 Figure D.6: Questem airborne total field magnetic data. Appendix D. Data 212 10000 -5000 -10000 -5000 0 5000 10000 15000 Easting (m) ppm -15000 -10000 -5000 0 5000 10000 15000 Easting (m) ppm (a) Inphase coplanar 466Hz (b) Inphase coaxial 1046Hz -10000 -5000 0 5000 10000 15000 Easting (m) ppm (c) Quadrature coplanar 466Hz -10000 -5000 0 5000 Easting (m) 10000 15000 PPm (d) Quadrature coaxial 1046Hz Figure D.7: In-phase and quadrature frequency-domain electromagnetic data from the Dighem system: 466Hz, 1046Hz, and 6495Hz Appendix D. Data 213 -10000 -5000 0 5000 10000 15000 Easting (m) -10000 -5000 0 5000 10000 15000 Easting (m) ppm ppm (a) Inphase coplanar 7206Hz (b) Inphase coplanar 56kHz 15000 -5000 0 5000 Easting (m) 10000 15000 -10000 -5000 0 5000 Easting (m) 10000 15000 ppm ppm (c) Quadrature coplanar 7206Hz (d) Quadrature coplanar 56kHz Figure D.8: In-phase and quadrature frequency-domain electromagnetic data from the Dighem system: 7206Hz and 56kHz Appendix D. Data 214 g 20000 Cin 15000 Questem Z-component data along flight line over the deposit. — Ch. 1 -Ch. 2 - Ch. 3 - Ch.4 • - Ch. 5 105 118 131 144 157 170 183 196 209 222 235 Ch. 6 - • - Ch. 7 Ch. 8 Ch. 9 Ch. 10 1 14 27 40 53 66 79 92 105 118 131 144 157 170 183 196 209 222 235 Ch. 11 - -Ch. 12 Ch 13 Ch. 14 Ch. 15 1 14 27 40 53 66 79 92 - 105 118 131 144 157 170 183 196 209 222 235 NE Station # SW Figure D.9: Profiles of Z-component dB/dT data for the fifteen time-channels (0 to 18 msec after the transmitter current turn-off) along line 10. The upper panel shows data for the first five time channels, the center plot shows data for the next five time channels, and the lower panel shows data for the last five time channels. The deposit corresponds to the strong negative response in the later time channels. Thirty-one lines of data were collected and different orientations but this is the only line that has been presented. Appendix D. Data 215 Figure D.10: Gradient Array IP data. The boundary shows the extension of the survey measurements and an outline of outcropping geology is also shown. Appendix D. Data 216 Figure D.ll: Reduced gravity data. This gravity data has had a regional signal re moved. The resulting gravity data were used for producing the density-contrast models in Chapters four and five. Observed Magnetic Data 614 data. I = 50.6, D = -13.4 . . i nT •2206 -2005 -1804 1603 1402 1201 -1000 Easting (m) Figure D.12: Total-field ground magnetic data. The observation locations are displayed and the gaps where the data have been discarded are clearly visible. Appendix D. Data 217 Line Observed Predicted ..li mati • 3200 Ohm-m 300 Figure D.13: Observed and predicted CSAMT apparent resistivity data. The data was collected along three lines, -200m N, -400m N, and -600m N. Each data plot is plotted with frequency in the vertical axis (ranging from 8192Hz at the top to 0.5Hz at the bottom), and easting along the horizontal axis (-2200m E to -800m E). Figure D.14: Observed and predicted CSAMT phase data. Each data plot is plotted with frequency in the vertical axis (ranging from 8192Hz at the top to 0.5Hz at the bottom), and easting along the horizontal axis (-2200m E to -800m E). Appendix D. Data 218 Figure D.15: Observed and predicted "Realsection" IP data. Five rows of data are plotted in each panel. Each row corresponds to a different transmitting dipole length. ix D. Data 219 O* CD L 0 O 5 cn O "t* - CM ' K <tf ID lo g BUS >H E g <a E S So » Q.E ,.br ID I 9 ^ * 3 «3 co o £ o H ° 'x — £ D 3 « S g 3 m O tf tf tf tf 3 6 » 3 Q E 'SEE %g £ ^ *t 0> 18"; E » o .5 E B i3 E o o Z oD CB E 0 O D o « 0 «j K S E QJ U kafe ,B «3 k E £ S o <u Appendix E Inversion In the body of the thesis the models have been presented without much discussion. Here we provide a brief description of the inversion codes used and the parameters that were used to invert each data set. E.l Introduction This section gives a brief overview of the algorithms used for the inversions performed in this thesis. All of the algorithms were developed by the UBC Geophysical Inversion Facility under the auspices of the JACI (Joint and Cooperative Inversion of Geophysical and Geologic Data) research consortium. E.2 GRAV3D GRAV3D is a collection of algorithms that are used for forward modeling, inversion, and sensitivity calculations of gravity data. The developers of this code give a description of the method in [Li and Oldenburg, 1998b]. As described in Chapter 2 the earth model is discretized into a large number of cuboidal cells of constant density contrast. The user can design a range of model objective functions to be minimized, subject to fitting the data, in order to recover a density contrast distribution. Functionality allows for incorporating a priori information through the weighting of reference models, and bounding the range of values that cells can assume. A depth weighting scheme that counteracts the natural 220 Appendix E. Inversion 221 decay of the gravitational source is needed to invert data that has no inherent depth information. This allows all parts of the model to play a role in the solution, and therefore introduces depth information into the model. E.3 MAG 3D MAG3D is a suite of algorithms used to invert surface and borehole magnetic responses due to a three-dimensional susceptibility distribution ([Li and Oldenburg, 1997]). When susceptible material is magnetized by an inducing field, an anomalous magnetic field is produced that is superimposed on the inducing field. The anomalous, or induced, mag netic field is inverted to recover the susceptibility distribution. The algorithms assume that susceptible values are small (a reasonable assumption for most geologic scenarios with paramagnetic materials), and that no demagnetized or remnantly magnetized ma terial is present. Residual magnetic field observations, locations and error estimates, the inducing field orientation and strength is required for inversions to proceed. The objec tive function minimized in the inversion allows various types of a priori information to be included. In addition weighting is incorporated to distribute the susceptibility in depth. Wavelet transforms are used for compression of the sensitivity matrix, allowing for faster run-times. The logarithmic barrier method ensures that positive magnetic susceptibility values can be enforced and bounds can be applied. E.4 CSAMT1D CSAMT ID is a near- and far-field 1-D code for recovering conductivity structures from E and H field components or the derived quantities of apparent resistivity and phase. The source electric field is a grounded loop and measurements are made in the frequency domain. Appendix E. Inversion 222 E.5 EM1DFM EM1DFM inverts airborne or ground EM frequency-domain data to recover simultane ously, or separately, 1-D conductivity and magnetic susceptibility models. The source and receivers are represented by magnetic dipoles. This is a fast inversion routine and can be applied to regional data sets ([Farquharson and Oldenburg, 2000]). E.6 TEM1D TEM1D is used to invert the vertical component of the time-domain electromagnetic data to recover a layered earth model. The forward modeling is done in the frequency domain and a frequency-to-time transform performed at each iteration. E.6.1 Complex forward modeling A modified version of the TEM1D forward modeling code was used to produce the re sponse over a polarizable layered earth. The Cole-Cole parameters are used to define the frequency-dependence of the model. E.7 DCIP3D DCIP3D is a suite of algorithms that are used to invert DC resistivity and IP data. The algorithms can be used in sequence. First, DCINV3D is used to invert DC re sistivity data (primary potentials) and produce a resistivity distribution of the earth. Secondly, the resistivity model can be used to calculate sensitivities for the IP inversion (IPSEN3D). Thirdly IP data (secondary potentials) can be inverted using IPINV3D. The model objective function is designed to allow for a priori information to be included. Appendix E. Inversion 223 E.8 EH3D EH3D is a forward modeling algorithm that computes frequency-domain electromagnetic responses from 3D conductivity and/or susceptibility distributions. The conductivity and susceptibility distributions can exhibit large discontinuities. The source can be either a grounded electric source, a loop source, or a plane-wave. Real and imaginary parts of electric and magnetic fields due to the source are computed over the whole model for different frequencies. E.9 EH3DTD EH3DTD is a forward modeling algorithm that computes time-domain electromagnetic responses from a given 3D conductivity structure of the earth. As with the frequency-domain equivalent a number of different sources can be used and the electric and magnetic fields are calculated throughout the whole model region. Different time-dependent wave forms can be used as sources. ix E. Inversion 224 Other parameters depth weighting, topogra phy included depth weighting, topogra phy, inclina-tion=50.638, declination=-13.47 (local grid), inducing field=44000nT depth weighting, topogra phy, inclina-tion=50.638, declination=-13.47 (local grid), inducing field=44000nT depth weighting, topogra phy, inclina-tion=50.638, declination=-13.47 (local grid), inducing field=44000nT Inducing current orienta tion: inclination 0°, decli nation 90° depth weighting, topography depth weighting, topogra phy, inclina-tion=50.638, declination=-13.47 (local grid), inducing field=44000nT Electrode weighting Inversion Time 33.2hr 2hr 11.5hr 38.5hr 17.75hr 3.2hr 6min/stn 9hr Achieved Misfit 482646 CO m CN cn tN CO OJ 4642 2654 in o CO OJ Reference Model E U bO o 0 S.I. xl0~3 0 S.I. xl0"3 CO eo 1 o O X 0 msec CO s u \ bO O co" co 1 o O X 100 Ohm-m 0 msec Cell Size 250 x 250 x 100m 100 x 100 x 50m 100 x 100 x 50m 100 x 100 x 50m 50 x 50 x 50m 25 x 25 x 25m 25 x 50 x 12.5m 1 - 10,000m X m OJ x a m m OJ OJ Number of cells 707625 167040 862580 338640 426400 446684 222880 o 117600 Noise estimation 0.01 mGal 4nT 4nT c 1 msec 1 mGal H a Ol 2-5% ap parent re sistivity and 1-2 de grees phase 5% + 1 msec Recovered model 3D density-contrast 3D magnetic susceptibility 3D magnetic susceptibility 3D magnetic susceptibility 3D chargeability 3D density-contrast 3D magnetic susceptibility ID layered-earth resistivity 3D chargeability Number of Data 17671 2566 3318 10115 6631 CO OO CO 60/stn 1182 Inversion Program GRAV3D MAG3D MAG3D MAG3D MAG3D GRAV3D MAG3D Q H 2 < CO O DCIP3D j Inversion Description Regional Gravity Dighem airborne magnetics BHP airborne magnetics Fugro airborne magnetics Gradient-array IP Local scale gravity Ground magnetics CSAMT | "Real-section" and gradient-array IP 'to Appendix F Gradient IP data modeled using MAG3D Gradient-array data has been modeled using MAG3D, an inverse modeling package for inverting magnetic field data. Similarities between the nature of data allow this. The advantage of inverting gradient-array IP data in this way is that depth information can be introduced in the model. As with potential field data, depth information is not contained in these data. Here we show why this method is valid, and test it with some synthetic examples. MAG3D generates 3-D susceptibility distributions based on the following forward modeling of the anomalous magnetic field from susceptibility distributions, K, in an inducing field [Li and Oldenburg, 1997]: where r is the point of observation of the magnetic field. This conservative field is derivable from the scalar magnetic potential [Grant and West, 1965], which is written as: Gradient-array IP data can be interpreted and modeled as magnetic potential field data because of the dipole source of both fields. Seigel [Seigel, 1959a] describes the over-voltage effect associated with a metallic conducting sphere in an electrolyte as equivalent (F.l) (F.2) 225 Appendix F. Gradient IP data modeled using MAG3D 226 to a dipolar source whose moment is proportional to the current density vector (j). The —> —* dipole source has a volume current moment of M = —mj, where m is the chargeability. The potential due to a distribution of current dipoles can be written as: ATTCT 1 r — 1 (F.3) where o is the conductivity of the medium. Expressions of potentials in equations F.2 and F.3 have the same form if the conductivity and electric field are uniform. The negative of the data is used, as the dipoles would be aligned in antiparallel with the inducing current, so the inducing field is set antiparallel with respect to the position of the current electrodes. F.l Four synthetic examples A synthetic example was used to test the use of MAG3D in this manner. It is expected that the conductivity of the deposit is anomalous with respect to it's surroundings, so several synthetic experiments were carried out with different conductivity distributions in order to determine if this method is valid for use at San Nicolas. A 100 msec chargeable cuboid was generated from 20m x 20m x 10m cells positioned between -60m and 60m east, -60m and 60m north, and 80m and 120m depth (Figure F.l). The body is hosted in a background of 0 msec chargeability. A current of 1 Amp was transmitted through grounded electrodes at two locations straddling the target (-600m east, 0m north and 600m east, 0m north). Eighty-four gra dient IP measurements were made using 12, 40m dipole locations on seven lines centered around the chargeable body. DCIPF3D, which is a forward part of the DCIP3D inversion package, was used to generate these measurements. Four different conductivity models were used to model this data (Figure F.2): Appendix F. Gradient IP data modeled using MA G3D 227 Figure F.l: Synthetic chargeable body. 1. a half-space of 0.01 S/m 2. a 0.1 S/m conductive body of the same size and location as the chargeable body in a background of 0.01 S/m 3. a 1 S/m conductive body of the same size and location as the chargeable body in a background of 0.01 S/m 4. a 0.1 S/m conductive body of the same size and location as the chargeable body in a background of 0.01 S/m with a 40m thick conductive overburden of 0.1 S/m (This model is thought to most closely resemble the conductivity distribution at San Nicolas.) The four synthetic data sets (one for each conductivity model) were contaminated with random noise with a mean magnitude equal to 3% of value of the data. Each datum was assigned an appropriate standard deviation (3% of each datum) and the data sets were inverted using MAG3D. An incident field with 0 degrees inclination, 90 degrees declination, and 1 nT field strength was used. Appendix F. Gradient IP data modeled using MAG3D 228 Figure F.2: North-facing cross-sections through the conductivity models: A) 0.01 S/m half-space, B) 0.1 S/m block in a 0.01 S/m background, C) 1 S/m conductive block in a 0.01 S/m background, D) 0.1 S/m block and a 40m thick 0.1 S/m overburden in a 0.01 S/m background. Appendix F. Gradient IP data modeled using MAG3D 229 01 50 100 mV/V Figure F.3: North-facing cross-section through recovered chargeability models. Figure F.3 shows the recovered chargeability models for the four models used. As expected, the half-space conductivity model gives the best results; the recovered values are closest to the true model value of 100 mV/V. As the conductivity increases, the recovered model deteriorates. The model with the conductive overburden is the worst. This is possibly because only a small amount of current was able to penetrate to the depth of the target. (A smaller current density will produce a smaller dipole source). In all cases there is good localization of the anomaly. F.2 Modeling the electric field As well as being affected by the conductivity structure, the electric field also must also be uniform in order for the treatment of gradient array IP data as potential field data. EH3D was used to forward model the electric field generated by the transmitter for each Appendix F. Gradient IP data modeled using MAG3D 230 of the conductivity models. The results are presented in Figure F.4, which shows the vector representation of the real part of the electric field for the north-facing cross-section at -400m north. The colors represent the magnitude of the field, and the line is orientated in the direction of the field vector with the white tip being forward. The electric field is distorted by interactions with the conductivity structure. It is not surprising to see the best recovered chargeability results coming from the uniform conductivity and electric field distribution, while the less successful results are obtained from conductivity structures that distort the electric fields. Although the electric field and conductivity structure isn't too anomalous in the case where a conductive overburden is added, the current is concentrated near the surface (Figure F.5). This current channeling reduces the signal relative to the noise, which is reflected in smaller amplitude data, and results in a recovered model that is reduced in magnitude. Appendix F. Gradient IP data modeled using MAG3D 231 Appendix F. Gradient IP data modeled using MAG3D 232 Appendix F. Gradient IP data modeled using MA G3D 233 F.3 Summary Gradient array IP data can be regarded as a potential field, and can be treated as pseudo-magnetic data if the conductivity structure over which the data are collected isn't anomalous enough to distort the electric field. Conductivity distributions that have contrast of two orders of magnitude still recover a useful chargeability model. Inversion of this data using MAG3D can produce good chargeability models as the data is weighted in order to provide depth information. In the case of San Nicolas the conductivity isn't too anomalous although the overburden has a detrimental effect on the data. Appendix G IP electrode weighting When inverting IP data using DCIP3D, accumulations of high chargeability values are sometimes observed coincident with the location of the transmitting electrodes. This electrode-effect is even more significant when the data have a small number of trans mitter electrodes. The data collected at San Nicolas used the gradient-array, and "Real-section" array electrode configurations; both of which have a small number of transmitter electrodes. As a result, the electrode-effect is observed in the inversion models. This section introduces an ad-hoc weighting procedure that can be used to reduce the electrode-effect. The weighting is invoked using the ws weighting matrix in Equation 2.2. The value of these weights are simply increased in the vicinity of the electrode. The cells around the electrode will be preferentially weighted closer to the reference model, which is set as zero. A weighting distribution is sought that smoothly decreases from a high value at the electrodes to a background weighting value of one away from the electrodes. The electrode weighting must effectively suppress charge accumulations around the electrodes but must also decay fast enough so as not to influence the true chargeability distribution. The following set of conditions was used to empirically develop a weighting function. when wm > c, when wm < c. (G.l) 234 Appendix G. IP electrode weighting 235 where, wm is the weighting value assigned to the mth component of the diagonals of the matrix ws (Equation 2.2), n refers to the nth transmitter electrode, N is the number of transmitter electrodes, rmn is the distance from the nth electrode to the center of the mth cell, and a, b, and c are all empirically derived constants, a increases the weighting value, b determines the rate at which the weight decays away from the electrode location, and c determines the value the weighting function can obtain before the default weights are used. The weighting around each electrode is superimposed to form a weighting function that exhibits hemispheres of high values below each transmitter electrode, a's of one, two, and three were applied to the synthetic model. A value of 100,000 was used for a, and a background value of one was assigned once the weights were less than 5 (the values assigned to c) (this cut-off was used to reduce the size of the volume that is effected by the electrode weights). Weighting functions were tested for a number of simple synthetic chargeability models with different values of b. A value of 2 produced results that seemed to have the least artifacts. The results of using this weighting function is shown in Figure G. The effect of implementing the weights not only reduces the high values around the electrodes but also increases the magnitude of the anomaly. This increase helps the model achieve values closer to that of the true model. Limitations of this weighting function are obvious. This weighting method might not be appropriate if the electrode were to be placed in chargeable material. Also when the electrode is very close to the center of the cell the method might produce undesired effects. An alternative would be to construct a sensitivity based weighting function. This would make parts of the model less sensitive to the data and therefore allow them to be closer to the reference model. Appendix G. IP electrode weighting 236 1000 Easting [m] (a) Synthetic chargeability model. 1000 Easting [m] (b) Inversion model 1000 Figure G.l: Removal of electrode-effect via weighting. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052682/manifest

Comment

Related Items