P r a c t i c a l I n v e r s i o n o f 3D T i m e D o m a i n E l e c t r o m a g n e t i c D a t a Applicat ion to the San Nicolas deposit by Scott R. Napier B.Sc. . The University of British Columbia, 2003 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F Master of Science in The Faculty of Graduate Studies (Earth and Ocean Sciences) The University Of British Columbia April 11, 2007 .© Scott R. Napier 2007 Abstrac t The purpose of this thesis is to provide a simple procedure by which time domain electromag-netic (TEM) data, which has been collected in the field, can be inverted to recover a three dimensional distribution of electrical conductivity. An algorithm, EH3DTDinv, which has been developed in order to invert T E M data in 3D is briefly presented. The algorithm has been pre-viously shown to work well on synthetic data. However, its ability to invert real data collected in field surveys was previously unclear. To illustrate the effectiveness of the inversion algorithm when used in conjunction with an inversion procedure, the example of the U T E M survey at the San Nicolas deposit in Zacatecas Mexico is used. The inversion procedure provides a detailed description of the preparations required for 3D T E M inversion. In order to clarify the individual stages in the procedure, the inversion preparations are outlined as they apply to the San Nicolas survey. Problems of non-uniqueness in the inverted models at San Nicolas are encountered and the solution is determined to be a combination of a sensitivity weighting scheme and the use of the data from multiple transmitters simultaneously in a single inversion. The iterative nature of the inversion procedure is demonstrated through the San Nicolas example. The San Nicolas inversion example culminates in the simultaneous inversion of data from 3 transmitters where a priori information from the regional geology has been incorporated into the result. The 3D in-version result is the first large scale 3D T E M inversion for a mineral deposit and agrees well with the geology and conductivities interpreted from drilling. Practical aspects of inversion are also considered and some recommendations on data formats, software utilities, and survey design are made in order to enhance future results. A brief study into inversion efficiency determined that a significant reduction in computation times can be achieved without loss of accuracy or resolution by fine tuning the cooling schedule of the algorithm. Further reductions in compu-tation times are made through parallelization of the forward modelling process. Finally, yet more reductions in computation times are possible through the coarsening of the spatial and temporal discretizations. However, this final method of computational reduction is shown to ii Abstract come at the cost of accuracy and resolution in the inverted model. The thesis concludes with a discussion of the potential for 3D T E M inversion and brief outline of future advancements in the field. iii Acknowledgements I would like to thank Teck Cominco Ltd. for providing the U T E M data, and SJ Geophysics Ltd. for collecting it. I would like to thank Nigel Phillips for his extensive work done on San Nicolas as a part of his Master's thesis. Phillips compiled the geologic sections and physical property data shown here in addition to inverting the numerous previous E M surveys which have been helpful in evaluating the current inversion results. Finally I would like to thank Doug Oldenburg for his guidance and support, without which this project would never have been possible. iv Table of Contents Abstract ii Acknowledgements iv Table of Contents v List of Tables viii List of Figures ix Nomenclature xi 1 Introduction 1 1.1 Maxwell's equations 1 1.2 Introduction to electromagnetic methods 2 1.3 Electrical conductivity 3 1.4 Role of electrical conductivity in geophysical methods 3 1.5 Practical electromagnetic surveying 4 1.6 Time domain interpretation 4 1.7 The time domain survey 5 1.8 Statement of purpose and outline of thesis 7 2 On three-dimensional T E M forward modelling and inversion 10 2.1 Introduction 10 2.2 Forward modeling with EH3DT 11 2.3 Inversion with EH3DTinv 13 2.4 Conclusion 16 v Table of Contents 3 A m e t h o d o l o g y fo r t h r e e - d i m e n s i o n a l T E M i n v e r s i o n 18 3.1 Introduction 18 3.2 A simplified schematic diagram 18 3.3 Background on the San Nicolas deposit 19 3.3.1 Geologic background 19 3.3.2 Geophysical background 20 4 O n r e q u i r e d i n v e r s i o n i n p u t a n d p r e p a r a t i o n s 25 4.1 Introduction 25 4.2 Understanding the data 25 4.2.1 Introduction to the U T E M survey 26 4.2.2 U T E M data reduction and normalization at San Nicolas 27 4.2.3 U T E M orientation conventions 30 4.2.4 Synthetic examples of U T E M data 30 4.2.5 The San Nicolas U T E M survey 34 4.3 Background modelling 35 4.3.1 Background models from a priori information 36 4.3.2 Background models from ID forward modelling 36 4.3.3 Background models from ID inversion 37 4.4 Discretization 39 4.4.1 Discretization in time 40 4.4.2 Discretization in space 43 4.4.3 Reducing the volume 44 4.5 Understanding the data through forward modelling 46 4.5.1 Validation of discretization 46 4.5.2 Forward modelling using a priori information 48 4.6 Assigning standard deviations to the Data 53 4.7 Conclusion 56 5 O n p r a c t i c a l t h r e e - d i m e n s i o n a l T E M i n v e r s i o n 57 5.1 Introduction 57 5.2 Inversion of synthetic data 57 Table of Contents 5.2.1 Inversion for a conductive block 59 5.3 Preliminary inversion of field data 62 5.4 Evaluation of preliminary results 63 5.4.1 Sensitivity weighting 66 5.4.2 Inversion of synthetic data using sensitivity weighting 68 5.4.3 Multiple transmitter synthetic inversion 69 5.5 Further improvements to inversion 71 5.5.1 Data editing 71 5.5.2 Refinement of discretization 72 5.5.3 Multiple transmitter inversion 73 5.6 Inversion of field data incorporating improvements 73 5.7 Evaluation of improved results 74 5.8 Inversion of field data using a priori information 77 5.9 Evaluation of results with a priori information 79 5.10 Conclusion 83 6 On practical and efficient T E M inversion 84 6.1 Introduction 84 6.2 Practical inversion of 3D T E M data 84 6.2.1 Survey design 85 6.2.2 Analysis of measurement errors 85 6.2.3 Data Formats 85 6.2.4 Construction of 3D models 86 6.3 Efficient inversion of 3D T E M data 86 6.3.1 Parallelization of forward modelling 87 6.3.2 Initial models 88 - 6.3.3 Cooling schedule for the regularization parameter 89 6.3.4 Discretization factors 90 6.3.5 Application to the San Nicolas deposit 91 6.4 Conclusion 93 7 Conclusion and Discussion 96 Table of Contents Bibliography 100 vm Lis t of Tables 3.1 Table of physical properties for the San Nicolas deposit 22 4.1 Mean delay times for the U T E M system operating at a base frequency of 30.974 Hz 34 ix Lis t of Figures 1.1 Source terms and their fourier transforms 6 1.2 Transmitter current waveforms for some common T D E M systems 9 3.1 A simplified procedure for T E M inversion 19 3.2 East-West oriented geologic cross section of the San Nicolas Deposit 20 3.3 Survey geometry for the San Nicolas U T E M dataset 21 3.4 U T E M data showing an inductive response to the San Nicolas deposit 23 4.1 A simplified procedure for T E M inversion 26 4.2 Schematic diagram showing locations of U T E M time windows 28 4.3 Transmitter current waveforms for the U T E M system 31 4.4 Comparison of U T E M data and a step-off response 32 4.5 Synthetic example of U T E M data 33 4.6 Forward modelling and inversion for background models 38 4.7 Step response of conductive sphere and halfspace 41 4.8 Steady-state and square wave responses for a 100 Vtm halfspace 42 4.9 U T E M transmitter current waveform 43 4.10 Validation of time and 3D spatial discretization 47 4.11 Forward modelling a conductive block 50 4.12 Models of San Nicolas with and without, deposit present 51 4.13 Forward modelled data for 3D San Nicolas conductivity model with and without deposit present, 52 4.14 Comparison of real and forward modelled data for San Nicolas 54 5.1 A simplified procedure for T E M inversion 58 5.2 Predicted and observed data for inversion of conductive block 59 x List of Figures 5.3 Convergence curve for synthetic block inversion 61 5.4 Inversion results for a conductive block . 61 5.5 Preliminary inversion result for San Nicolas 63 5.6 Predicted and observed data for preliminary inversion 64 5.7 Inversion of synthetic data with and without sensitivity weighting 70 5.8 Inversion of synthetic data with data from multiple transmitters 71 5.9 Inversion of 2 sources with improved parameters 74 5.10 Predicted and observed data for improved inversion 75 5.11 Inversion of 3 sources with a priori information 80 5.12 Inversion of 3 sources with a priori information 81 6.1 Convergence curve with parallelization and efficient starting regularization . . . . 92 6.2 Comparative profiles through efficient inversion results 94 xi Nomenclature Lat in Letters A vector potential y_ m equation (2.3b), 12 A ( m ) forward operator matrix various equation (2.5), 13 o j jth U T E M channel normalized at location Xj , yj, Zj 'unities s equation (4.1), 27 ci i t h U T E M channel at location T s equation (4.1), 27 d o b s observed data various equation (2.7), 13 predicted data various equation (2.6), 13 D electric displacement field C m 2 equation (1.1), 1 E electric field y_ 71) equation (1.1), 1 E° electric field at time to v_ m equation (2.0), 11 E n electric field at time tn y_ 11} equation (2.1d), 11 g gradient of objective function equation (2.12), 14 H approximate Hessian (iogsy°2 equation (2.14), 15 H Hessian (logS)- 2 equation (2.12), 14 H magnetic field A_ in equation (1.1), 1 H ° , magnetic field at time to V_ in. equation (2.0), 11 H n magnetic field at time tn A_ m equation (2.1d). 11 H P ( x j .> •j.Zj) ith U T E M primary magnetic •'•>//./•-./ field at location of A in equation (4.1), 27 J approximate sensitvity unities s equation (5.4), 65 J c current density A ni2 equation (1.1), 1 J s source current density A in'2 equation (1.1), 1 M preconditioning matrix unitless equation (2.16), 15 h normal vector unities s equation (2.0), 11 xii Nomenclature p model perturbation log(S) equation (2.12), 14 Q data sampling matrix various equation (2.6), 13 q sources and BC terms various equation (2.5), 13 R regularization function various equation (2.7), 13 current source A in2 equation (2.Id), 11 m model log(S) equation (2.5), 13 So corrective sources various equation (4.7), 43 u fields various equation (2.5), 13 w d data weights various equation (2.7), 13 w s cell volume based weighting matrix unitless equation (2.10), 14 w t user defined weighting matrix . unitless equation (2.10), 14 w x horizontal finite difference matrix unitless equation (2.10), 14 w y horizontal finite difference matrix unitless equation (2.10), 14 w z vertical finite difference matrix unitless equation (2.10), 14 wt weighting function unitless equation (2.8), 14 Greek Letters a reciprocal of time step 1 equation (2.Id), 11 3 regularization paramter unitless equation (2.7), 13 e dielectric permittivity j ~ 7 equation (1.2), 1 7 / line search parameter unitless equation (2.15), 15 //, magnetic permeability equation (1.4), 1 CT conductivity ^ equation (1.3), 1 Pf free charge density equation (1.1), 1 <f> objective function unitless equation (2.7), 13 d> scalar potential V equation (2.3b), 12 X diagnostic for convergence unitless equation (2.17), 16 xiii Chapter 1 I n t r o d u c t i o n 1.1 Maxwell 's equations Maxwell's equations in their modern form as composed by Ward & Hohmann ( 1 9 8 7 ) are shown in equation 1.1. These equations, along with the constitutive relations, from equations 1.2, 1.3 and 1.4, provide the basis of the theory behind electromagnetic surveying. V x E + — = 0 ot V x H - — = J C + J S 3 D dt V - B = 0 V - D = P / (1.1) Here E is the electric field, B is the magnetic flux density, H is the magnetic field, J c is the current density, J s is the source current density and D is the electric displacement field. The quantity pf is the free charge density. The constitutive relations are D = eE (1.2) where e is the dielectric permittivity J C = C T E (1.3) where a is the electrical conductivity, and H = ^ B (1.4) where p is the magnetic permeability. 1 Chapter 1. Introduction The constitutive relations relate the various field quantities to each other through the key physical properties relevant to electromagnetics, electrical conductivity, magnetic permeability and dielectric permittivity. 1.2 I n t r o d u c t i o n to electromagnetic m e t h o d s It is possible to define geophysical electromagnetic (EM) methods as consisting of any type of survey which uses a time dependent source of electric or magnetic fields and proceeds to make measurements of electric or magnetic fields, or alternatively their time derivatives, in order to make inferences about physical properties within the ground. Maxwell's equations along with the constitutive relations and both initial and boundary conditions contain all the information necessary for precisely predicting a field measured in a receiver given an arbitrary source. In practice, however, these equations are difficult to solve in many realistic situations where arbitrary, and sometimes complicated distributions of electrical conductivity, magnetic permeability and dielectric permittivity are present. The solution of Maxwell's equations for a specified model of the physical properties of the earth is referred to as the forward problem of electromagnetics. Further detailed discussion on the E M forward problem will be undertaken in chapter 2. Generally when electromagnetic surveys are conducted in typical earth materials the phys-ical properties of magnetic permeability and dielectric permittivity are much less important to the observed data than the electrical conductivity. In the case of magnetic permeability this is because, with the notable exceptions of the minerals magnetite and pyrrhotite, the magnetic permeability of common earth materials does not vary widely. In the case of dielectric per-mittivity, most geophysical E M surveys either utilize sources which are limited to a frequency band or take place over a material which has a conductivity such that the contribution of the dielectric permittivity is at best secondary. For most E M surveys the electrical conductivity is the dominant factor in the earth which affects the fields recorded by the receiver. This of course, is stated with exceptions, the most notable of which is one of the most common E M methods, ground penetrating radar (GPR). G P R takes place at frequencies high enough that the dielectric permittivity becomes very important to results. This thesis, however, confines itself to situations where the most important physical property of interest for electromagnetic methods is the electrical conductivity of the ground. 2 Chapter 1. Introduction 1.3 E l e c t r i c a l c o n d u c t i v i t y Electrical conductivity is defined by the Oxford English Dictionary to be the quality of "[GJuidance, conveyance ... [r]elating to or connected with electricity" O E D (1989). In the field of geophysics the electrical conductivity is thought of as the ability of rocks to transmit electrical energy. Rocks themselves can be defined as consisting of an aggregate of mineral grains and the fluids in the pore spaces. The conductivity of such a material can be complicated. Electrical conductivity can be largely dependant on the inter-connectivity of minerals which have a dom-inant conductivity or it may be dependent on the conductivity and inter-connectivity of fluids in the pore space. For these reasons electrical conductivity is often a scale dependant quantity as shown by Bahr (1997) and Jones (2005). Electrical conductivity of rocks is also dependant in a complicated manner on the frequency content of the source as shown by Fuller & Ward (1970). Telford et al. (1990) show that electrical conductivity can vary among earth materials by more than 12 orders of magnitude. This variation presents its own difficulties from the point of view of interpretation, nevertheless, this extreme variation also makes conductivity a useful diagnostic tool. Designing a survey which is sensitive to the electrical conductivity of the earth has proven to be very useful in a variety of geophysical applications. 1.4 R o l e of electrical c o n d u c t i v i t y in geophysical m e t h o d s There are numerous examples of geologic situations where determination of the electrical con-ductivity distribution within the ground can be useful. These situations include, but are not limited to, the location and delineation of mineral deposits as in Macnae (1985), hydrocarbons as in Constable Sz Wiess (2006), aquifers and water resources as in Fitterman & Stewart (1986), pollutants and dangerous or hazardous chemical spills as in Atekwana et al. (2000), and even imaging large scale crustal and mantle structures in the earth as in Unsworth et al. (2005). Any method which is sensitive to variation in electrical conductivity has the potential to be useful in these types of situations. Electromagnetic methods are important because the method has a sensitivity to conductivity structures at a wide variety of depths. Electromagnetic methods are used to define and detect targets, from the centimeter scale at one end of the spectrum, in G P R , to mantle structures at tens or even hundreds of kilometers of depth, in the case of 3 Chapter 1. Introduction Magnetotellurics. The type of survey that is performed depends on the depth to the target, its expected size, its conductivity and its conductivity contrast with the host. 1.5 Practical electromagnetic surveying Electromagnetic surveys require two main components, a source and a receiver. Sources can be either a time-varying magnetic field, commonly generated by a current loop, or a time-varying electric field which is usually generated by a configuration of one or more grounded electric poles. Magnetic sources can be oriented in any direction and can be inside, on, or above the earth. Electric sources are usually grounded on or within the earth. Receivers are generally anything capable of measuring the electric or magnetic field, and are located either above, on, or inside the earth. The magnetic field time derivative is usually measured using coil sensors while the magnetic field itself can be measured using fluxgate magnetometers or SQUIDs. The electric field is generally measured using grounded electric dipoles. Ground based electromagnetic methods in geophysics emerged in the 1920's in Scandinavia, the United States and Canada in the exploration for shallow base-metal deposits. Airborne E M methods emerged some 30 years later, around the time when the method of Magnetotellurics, which uses a natural source field, was being developed as a geophysical method Telford et al. (1990). The modern electromagnetic survey consists of a wide range of possible surveys and are practical for a variety of applications. Survey types can run from simple metal detectors, at one end of the spectrum, to sophisticated 3D arrays. These arrays can be used to measure electric and magnetic fields generated from both natural and man-made sources and have the ability to detect E M anomalies from objects near surface and at great depth. 1.6 Time domain interpretation What is meant by time domain electromagnetic surveys is that the data are interpreted directly, as a time series, as opposed to being interpreted by Fourier analysis. There is a major division in the way in which electromagnetic surveys can be interpreted. It is possible to transform Maxwell's equations from the time domain into the Fourier, or frequency domain. Any method which make use of this type of analysis is referred to as frequency domain electromagnetics ( F E M or F D E M ) . Any method which interprets the data using Maxwell's equations directly as 4 Chapter 1. Introduction shown in equation 1.1 is referred to as time domain electromagnetics ( T E M or T D E M ) . In practice the choice of which type of analysis to use depends largely on the time-frequency dependance of the source. Sources which are band-limited are generally interpreted in the frequency domain, while sources with a wide-band are often interpreted in the time domain. Figure 1.1 provides an illustrative example. It shows two source functions and their respective Fourier transforms. It is clear from the Fourier transforms of the source functions shown that while the sinusoid is band-limited, to essentially a single frequency, the step function has a wide-band Fourier domain representation. The decision as to which method to use is generally one of efficiency. The data collected using a sinusoidal source are most efficiently interpreted in the frequency domain while data collected using the step source may be better interpreted in the time domain. This does not mean that a particular method can not be used in either situation, it is merely a question of which method provides the maximum simplicity and the greatest ease of use. One advantage of frequency domain interpretation is that the time derivative terms in Maxwell's equation are simplified, since time derivatives are merely multiplication in the Fourier domain. The advantage of using a wide-band source such as the step-off shown in Figure 1.1, is that E M signals can be measured in the absence of a primary field. This can have advantages from a signal to noise standpoint. 1.7 The time domain survey As we have discussed, the time domain survey generally uses a source waveform which has a relatively broadband spectrum. Many time domain surveys use a transmitter waveform which consists of both an on-time and an off-time. The off-time, when the transmitter has zero current, is useful since the only fields which are generated at this stage are those due to diffusive currents traveling in the earth and the air. These fields are termed the secondary fields, as opposed to the primary fields, those generated by the transmitter directly. It should be clear that the secondary fields are those that are of real interest to the geophysicists. When trying to characterize the conductivity of the earth the relative strengths of the two types of fields is a major consideration. Reference to Figure 1.1 shows that for the sinusoidal source there is no significant off-time for the transmitter; thus the secondary field must be determined always in the presence of the primary field. This tends to negatively affect signal to noise ratios 5 Chapter 1. Introduction a) b) c CD 3 O 1 0.5 0 -0.5 -1 II V c) c CD 3 o 1 ' 0.8 0.6 0.4 0.2 0 • -0.2 • •:r: ; ii;: t ime [s] 0 t ime [s] : ! : ^ CD 1 ~ O-E CO 120 100 80 60 40 20 O'-0 50 100 150 200 angular f requency [rad/s] d) 140 120 CD 100 -*J —. 80 mpl 60 CO 40 20 0 50 100 150 200 angular f requency [rad/s] Figure 1.1: a) Source with sinusoidal variation in current b) Amplitude of fourier transform of sinusoid c) Source with step variation in current d) Amplitude of fourier transform of step function 6 Chapter 1. Introduction of the recovered fields when compared to a source which has an off-time component. Figure 1.2 shows the transmitter current waveform for several commonly used time domain systems. Note that the U T E M system does not have an off-time component to its source. The reason for this will be discussed in detail in chapter 3. 1.8 Statement of purpose and outline of thesis The E M forward problem predicts the data given a model of the electrical conductivity of the earth. The geophysicist is almost always faced with the situation where the data have been recorded yet what is truly desired is a model of the electrical conductivity in the earth. The goal of the inverse problem for time domain E M is to attempt to find a model of electrical conductivity in the ground in the area where the E M were data recorded. This tends to be a difficult problem to solve. Historically the interpretation of surface and borehole time domain data has been hampered because of the inability to invert data to recover a 3D conductivity model. Past interpretations have often been limited to a direct interpretation of the data or at best the use of simplistic plate modeling and inversion or combining the results of ID inversions. However, for many geologic situations where complicated conductivity structures are present, these approaches have proved to be insufficient. An algorithm for the 3D inversion of time domain E M data has been previously developed at the University of British Columbia's Geophysical Inversion Facility by Haber et al. (2007). This algorithm has been shown to be useful on synthetically generated data. However, a comprehensive demonstration of the utility of the algorithm on a realistic field data set has yet to be completed. The purpose of this thesis is to show that: • Time domain electromagnetic data collected in a field survey can be inverted to recover a 3D distribution of electrical conductivity. • The inversion of T D E M data can yield substantial benefits for interpreters and that the pro-cess of inversion of T D E M data is practical and could become routine for many geophysical applications. • T D E M inversion is useful in situations where the geology is complicated and varies in three dimensions, and where interpretation of the data in a traditional fashion becomes difficult. 7 Chapter 1. Introduction This chapter of the thesis has begun with an introduction to E M methods and discussed the relevant physical properties for E M surveying. The thesis will continue in chapter two with a brief description of the forward modeling and inversion algorithms used in this study. The discussion will attempt to summarize the key elements of the algorithms, making reference to the published literature on the subject for the details. In chapter three the thesis will proceed to outline an iterative procedure for inverting 3D time domain data. To illustrate the effectiveness of the 3D inversion procedure, an archetypal example, the U T E M survey at the San Nicolas deposit in Zacatecas Mexico, will be used throughout. This is an excellent study area since the geology is complex but well understood and numerous other geophysical surveys have been performed and interpreted at the site. The San Nicolas deposit will be introduced and the U T E M survey conducted there will be described in detail. In chapter four the thesis will discuss the necessary steps that must be completed to pre-pare a field data set for inversion. The topics will include developing a background model, discretization of the spatial domain, discretization of the source waveform in time, reducing the spatial domain and assigning standard deviations to the data In chapter five the thesis will present a series of inversion results along with interpretations of those results for the San Nicolas deposit. The iterative nature of the inversion procedure presented here will become clear as successive inversions completed on the San Nicolas U T E M data set are presented. Because of the large scale of the numerical computations involved in 3D inversion the process of inversion can be very time consuming. In chapter six the topic of improving the efficiency of the inversion will be discussed along with recommendations on possible improvements to data acquisition and processing. Finally the thesis will conclude with a summary of the achievements and provide some recommendations on future work in the field of time domain E M inversion. 8 Chapter 1. Introduction a) =t 1 CD 0.5 13 O 0 . . • X „ . n — . — b) 3 time [s] x 10 ' < 1 CD op -1 c) 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 time [s] CD 0.5 3 O 0 -0.03 -0.02 -0.01 0.01 0.02 0.03 time [s] Figure 1.2: a) Transmitter current waveform for the airborne M E G A T E M system b) Transmit-ter current waveform for the ground based U T E M system operating at 30 Hz c) Transmitter current waveform for the ground based E M 5 7 system. 9 Chapter 2 O n three-dimensional T E M forward model l ing and inversion 2.1 I n t r o d u c t i o n In order to pursue the practical application of 3D time domain inversion to realistic data sets it is necessary to discuss the theoretical background for the forward modelling and inversion algorithms that will be used in this study. The forward modelling algorithm used here is E H 3 D T D and the inversion algorithm is EH3DTDinv. They are both codes developed by the University of British Columbia's Geophysical Inversion Facility in 2003 and 2005 respectively. In both the forward modelling and the inversion algorithm the earth is represented as a set of cuboidal cells of constant conductivity and the prediction of the data from an earth model can be calculated by solving Maxwell's equations using a finite volume method. The inversion methodology is similar to that used to invert 3D gravity, magnetic, or D C / I P data, although the scale of computations is considerably larger. Receivers can be in the air, on the ground or in boreholes, and they can measure either electric or magnetic fields in the on-time or off-time of the transmitter. The codes also include the capability to model the time derivative of the magnetic field. Because computation time increases substantially with the number of transmitters, the algorithm is most effectively used in ground-based geophysics where there are only a few source locations. The discussion of forward modelling and inversion within this chapter is by no means com-prehensive and some basic understanding of the processes of geophysical forward modelling and inversion is necessary in order to understand these topics in the format presented here. Non-linear inversion is discussed at length in Nodecal & Wright (1999) and Parker (1994). Ad-ditionally for a much more detailed view of the two algorithms used here readers are referred to the papers by Haber et al. (2004) and Haber et al. (2007). 10 Chapter 2. On three-dimensional TEM forward modelling and inversion 2.2 Forward modeling with E H 3 D T The details of the forward modelling techniques used in EH3DTD are outlined in the paper by Haber et al. (2004) and the references therein. Beginning with Maxwell's equations as shown in 1.1 we include the constitutive relations to eliminate displacement field and magnetic flux density. Initial conditions E°, H°, are introduced along with boundary conditions, n x H = 0, and the initial conditions are propagated through time using a Backward Differentiation Formula (BDF). The algorithm is capable of BDF schemes which are first order and second order accurate in time. The first order version of the BDF scheme is widely known as the backward Euler method. This is the method which will be presented here for the sake of simplicity of appearance; however the implementation of the second order scheme proceeds in an analogous way. Employing the BDF scheme and semi-discretizing Maxwell's equations over the time interval [ i n _ i , i n ] leads to the system of equations V x E n + a n ^ H n = a^ H " " 1 (2.1a) V x H n - (a + a ne)E n = srn - QneE""1 (2.1b) h x H" = 0 (2.1c) where an *n — * n - i A i At this stage a process of substitution and elimination in the equations 2.1a and 2.1b yields the curl-curl system in equation 2.2. V x ^ V x E n + Q n r j E n + a*eE n = o n(V x H n _ 1 + a n e E n _ 1 - sTn) (2.2) In the spatial domain, the earth is discretized into cuboidal cells and a finite volume tech-nique is used to calculate the electromagnetic fields on a staggered grid where E is defined on cell faces and H is defined on cell edges. Equation 2.2 can be difficult to solve due to the existence of a nullspace in the highest order operator, the curl-curl term. For this reason the 11 Chapter 2. On three-dimensional TEM forward modelling and inversion electric field is decomposed into electromagnetic potentials A and <f> defined by E = A + Vc/> (2.3a) V - A = 0 (2.3b) The divergence free nature of the vector potential term A, as shown in equation 2.3b, elim-inates the nullspace difficulties which are associated with the curl-curl system in equation 2.2. The effect is that this transformation helps to stabilize the equations and simplifies solution process. The formulation in terms of potentials has the consequence that in the spatial dis-cretization A is defined on the faces of cells while (j) is defined at each cell center. The problem can be composed in terms of A and (f) as follows where S(o) = a + ae and s is a term containing sources and boundary conditions. The system of linear equations in 2.4 must be solved in order to give the fields everywhere inside the discretized volume at some time tn provided that the fields at time tn-\ are known. The solution to Maxwell's equations can be propagated forward through time in this manner at the discrete times specified by a discretization of the transmitter current waveform. The advantage of the B D F method is that it is unconditionally stable and tends to rapidly attenuate high frequency components of discretization error as discussed in Haber et al. (2004). This is an important quality for stiff equations such as Maxwell's equations. The implicit B D F method requires a solution to a large matrix-vector system at each time step. The solver used here is a conjugate gradient routine. This can be time consuming when compared to explicit methods which require only matrix multiplication at each time step. Despite the computational demands, the advantages of an unconditionally stable and, as a consequence, much more flexible implicit schemes, are thought to outweigh the more restrictive and potentially unstable explicit schemes. V 2 - aiiS{o) -apS(o)V \ (2.4) V - S ( t r ) V - 5 ( C T ) V J 12 Chapter 2. On three-dimensional TEM forward modelling and inversion 2.3 Inversion with E H 3 D T i n v The inversion algorithm is described in accepted, but as yet, unpublished work by Haber et al. (2007). The system of equations for forward modelling shown in equation 2.4 can be abbreviated to a simpler form A(m)u = q (2.5) where m is the conductivity model of the earth. A(m) is the forward modelling matrix which depends upon the model, u contains the potentials, A and 4>, and q contains sources and boundary conditions. A common technique in E M inversion is to define the model as m = log(cr), see Whittall & Oldenburg (1992). This is done because electrical conductivity is a strictly positive quantity and earth materials tend to display a wide dynamic range in values of conductivity. The use of m = log(cr) helps to ameliorate problems associated with both of these issues. The predicted data, d p r e d , are sampled from the forward modeled potentials in the vector u by a transformation matrix Q . d p r e d = Qu (2.6) The inverse problem is defined as the minimization of the functional defined as min $ — — 2 2 0, W d ( Q A ( m ) - ' q - d o b s ) + ^ R ( m - m r e f ) (2.7) where d o b s are the observed data and m r e f is a reference model. Q is the unknown trade-off or regularization parameter which controls the relative importance between the first and second terms in equation 2.7. The two terms are often referred to as the data misfit, (pa, and model objective function, (f>m, respectively. Here is data weighting matrix whose diagonal elements consist of the reciprocals of the estimated standard deviations, ej, for each datum. R is a regularization term, which can in principle be any measure of any property of the model. The discrete regularization function, R(m — m r ef), used in the algorithm EH3DTDinv is the discretization of the continuous integral form 13 Chapter 2. On three-dimensional TEM forward modelling and inversion R(m — mre{) = / [ W f ( m — m r e f ) ] d Q + Jn Jn + Wi d(m - m r e f ) dy 12 dQ + Wtd(m-mTe{)dn,2 ox d(m - m r e f ) dz Wi dQ + (2.8) where Wt is a user defined weighting function and Q is the inversion domain. In the algorithm EH3DTDinv the discrete representation of this function is denned as R ( m - m r e f ) = asHWsWtCm-mref)!^+ a x | | W t W x ( m - m r e f ) | | | + (2.9) Oy | | W t W y ( m - mref)||| + etz | | W t W z ( m - mref)||2 here W t is a user defined model weighting matrix, W s is a diagonal weighting matrix containing cell volumes and W x , W y , W z are first order finite difference matrices of the corresponding orientation multiplied by a cell volume term. The minimization of the objective function, <£, is achieved by the Gauss-Newton approach. The Gauss-Newton method is related to Newton's method in which the objective function is expanded in a Taylor series around the current model. $(m + <5m) = #(m) + [V<£(m)](5m + ^ [ V 2 $(m)]8m2 + ... (2.10) All terms beyond second order are neglected. This is equivalent to approximating the objective function as a quadratic function. It is now possible to derive a model perturbation, <5m, which yields the minimizer of this quadratic function. Differentiating the system in equation 2.10 with respect to 5m and setting the right hand side equal to zero yields [V 2<£](5m = - V # (2.11) Often for the sake of simplicity and appearances the equation shown in 2.11 is displayed in the somewhat simpler notation shown below Hp = -g (2.12) 14 Chapter 2. On three-dimensional TEM forward modelling and inversion where p is the model perturbation <5m, H is the Hessian or V2(3?, and g is the gradient of the objective function or V<&. The gradient of the objective function can be written explicitly as g = J T W j W d ( Q A ( m ) - ' q - d ° b s ) + / R ( m g ~ m r e f ) (2-13) where J is the Jacobian and is defined as d^Ag^—^- For most realistic problems the true Hessian is time consuming to calculate and makes excessive demands on computing memory. The Gauss-Newton method helps to alleviate this problem. The method proceeds by making an approximation of the Hessian matrix often denoted as H . The approximation neglects second order derivative terms in the Hessian itself to give H = J T w T w d J + [3Wl W m (2.14) Given these quantities, namely an approximation of the Hessian and the gradient of the objective function, the resulting Gauss-Newton system can be solved with standard methods for linear systems. In the case of EH3DTDinv a stabilized bi-conjugate gradient (BICGSTAB) algorithm of Saad (1996) is used. The true objective function is not necessarily quadratic, however the solution to this system should give a perturbation to the model which will cause a reduction in the value of the objective function provided the objective function is convex. In order to ensure that a reduction is made in the objective function, the algorithm employs a line search technique where the updated model m n is given by m" = m""1 + np (2.15) The choice of the parameter rj is handled through a line search where n G [0,1]. The conjugate gradient solutions to the Gauss-Newton system in equation 2.12 can be time consuming, however the speed of such a calculation can be improved by the use of a pre-conditioning matrix M . M ~ J H p = - M " J g (2.16) A good pre-conditioner for the Gauss-Newton system in 2.12 is a matrix which is quickly in-vertible while remaining a close approximation to the true Hessian. The Quasi-Newton method 15 Chapter 2. On three-dimensional TEM forward modelling and inversion using BFGS vectors from Dennis Sz Schnabel (1996), allows the calculation an approximate Hessian matrix. This is done by the differencing of the gradient of the objective function at successive stages of the inversion process. The BFGS method utilizes these differences to create a positive definite approximation of the Hessian, which is easily invertible to give M _ i . This approximation can be shown to be close to the true value of the Hessian, given that the steps are not too large, several steps are incorporated, and the current solution is, in some sense, near the true solution. This may seem somewhat restrictive but in practice the Quasi-Newton BFGS method provides an efficient method of finding an effective pre-conditioner. Finally the determination of the unknown regularization parameter, j3, is handled through an iterative cooling process. A target misfit is defined and a successive series of inverse sub-problems are solved with a decreasing value of ft until the resulting solution has a misfit which is sufficiently close to the target value. The target can be defined through the use of the parameter X2 and a diagnostic for convergence is when x 2 = U though this is a parameter which may be adjusted by the user. X2 is defined as and here n is the number of data. This iterative cooling technique is particularly useful for highly non-linear inverse problems as it helps to avoid local minima in the objective function. By starting with a heavily weighted quadratic model objective function and cooling ft gradually, the complicated non-linear data misfit function becomes gradually dominant. Proceeding carefully in this way results in a much lower likelihood of becoming trapped by local minima in the objective function. 2.4 C o n c l u s i o n The forward modelling and inversion algorithms developed at U B C for the 3D inversion of time domain data inversion have been previously shown to work well on synthetically generated data in Haber et al. (2007) but there are numerous practical aspects that need to be dealt with when working with field data. The purpose of this thesis is to outline the practical considerations that must be made in order to apply this algorithm to actual geophysical data recorded in a (2.17) 16 Chapter 2. On three-dimensional TEM forward modelling and inversion field situation. It is desirable to develop a procedure by which the geophysicist may routinely obtain inversion results which are useful and enhance the understanding of the local geology. 17 Chapter 3 A methodology for three-dimensional T E M inversion 3.1 Introduction The primary purpose of this thesis is to outline a procedure that can help others successfully complete 3D T E M inversions and ensure that the results are useful for geophysicists. 3D inversion of time domain data is a complicated procedure with many interrelated factors which affect the quality of the results. Combining an algorithm of this level of sophistication with a field data set, often having many of its own complications and sources of uncertainty, can be difficult. The intention is to provide a simple procedure, which when followed carefully, will allow 3D inversion of time domain electromagnetic data, not only to produce useful results, but also produce results in a timely fashion. Eventually we imagine such a procedure will aid in the process of making 3D time domain inversion routine for use on almost any field data set. 3.2 A simplified schematic diagram Figure 3.1 shows a simplified schematic of a procedure for the 3D inversion of time domain electromagnetic data. Defining a procedure that describes a process that has this level of sophistication could easily result in something that is too complex to easily comprehend or is too cumbersome to use effectively. For this reason the procedure is presented as simply as we believe possible. Each step in the procedure will be discussed in the following chapters within the context of an archetypal inversion example, the inversion of U T E M data recorded at the San Nicolas deposit. For this reason, in order to further discuss this inversion procedure, it is necessary to first introduce both the San Nicolas deposit and the U T E M survey conducted over it. 18 Chapter 3. A methodology for three-dimensional TEM inversion understanding the data - • f \ * background model j i discretize * | j validate forward model j I { error assignment ; invert data i evaluate results j Figure 3.1: A simplified procedure for the inversion of 3D time domain electromagnetic data 3.3 Background on the San Nicolas deposit All of the following geologic sections and results are presented in a local coordinate system. A n approximate transformation to U T M zone 14 coordinates is provided in the appendices of Phillips (2001). 3.3.1 Geologic background San Nicolas is a Copper-Zinc volcanogenic massive sulphide deposit located in central Mexico in the state of Zacatecas. The deposit is a continuous but geometrically complex body of sulphides which is covered by 175-250 m of variable composition overburden. The local geology is somewhat complex and contains numerous sedimentary and volcanic units. Figure 3.2 shows an east-west cross section through the deposit at 400 m south as interpreted from the drilling program. Major geologic units are labeled. In previous work Phillips (2001) has determined that the deposit is hosted in a series of mafic and felsic volcanic rocks of the Upper Chilitos Formation, which lie unconformably over graphitic mudstones of the Zacatecas Formation. The deposit is entirely bounded to the east by a southwest-dipping fault. Mineralization continues to follow the fault at depth in a deeper a priori information 19 Chapter 3. A methodology for three-dimensional TEM inversion -2000 Easting (m) -1100 Figure 3.2: Geologic cross section of the San Nicolas deposit located at 400 m South. Geology-interpreted from drilling by Phillips (2001). part of the deposit known as the keel. The deposit is flanked by a thick package of rhyolites to the west. The volcanic units that envelope the deposit are overlain by a Tertiary-aged breccia overburden, which varies in thickness from 50 to 150m. Finally an overlying thin veneer of quaternary alluvium is present in the vicinity of the deposit. The upper portion of the deposit is a polymetallic sulphide with some gold and silver present. The lower portion of the deposit is copper-rich. The reserve estimates as of 2001 were 72 million tonnes grading 1.4 % copper and 2.27 % zinc. At present the deposit has not been mined due to the depth of the deposit, the lack of other proximal ore bodies, and mineralogical issues which affect the efficiency of the milling process. 3.3.2 Geophysical background The San Nicolas deposit has a long history of being a test case for the application of geophysics to the exploration for massive sulphides. The deposit has been investigated in the past using a wide array of geophysical surveys including gravity, magnetics, DC resistivity, induced polarization, airborne E M , and ground-based C S E M / C S A M T surveys. The deposit has also been thoroughly drilled and the recovered drill core has been logged for physical properties. As a result, both the geology and physical properties of the deposit are reasonably well understood. The sulphide deposit presents a electrical conductivity contrast with most of the geologic 20 Chapter 3. A methodology for three-dimensional TEM inversion c o 400 200 0 -200 -400 -600 -800 -1000 -1200 • 1 0 0 • 4 • I • 11 • ! -2500 -2000 -1500 -1000 -500 Easting [m] Figure 3.3: Survey geometry of the San Nicolas U T E M survey showing the location of the 3 large transmitter loops from the San Nicolas survey and their corresponding colour coded receiver locations. Blue is loop 1, red is loop 2, and green is loop 9. 21 Chapter 3. A methodology for three-dimensional TEM inversion units in the area, however the main overburden unit, the tertiary volcanic breccia, has a con-ductivity close to that found in the sulphide. This, combined with the depth to the target, makes E M interpretation at San Nicolas challenging. Table 3.3.2 contains a summary of typical values for the various physical properties of the deposit and its host rocks taken from Phillips (2001). These values are tabulated as a result of borehole surveys and core sample testing. Unit Name Physical Property Density Susceptibility Chargeabillity Resistivity g/cc SI x l 0 ~ 3 msec Qm Tertiary Brecias .2.3 0-5 10-30 20-30 Mafic Volcanics 2.7 0-5 30-50 80 Massive Sulphide 3.5 10 200 20 Quartz Rhyolite 2.4 0-10 10-20 100 Table 3.1: Table of physical properties for the San Nicolas deposit A U T E M survey was conducted over the San Nicolas deposit in December of 1998. A diagram showing the location of the receivers and transmitters used in this study is shown in Figure 3.3. A further 12 smaller loops with only a few data locations were deployed over the deposit in the original survey, however they are not considered in this study. The reasons for this are twofold. As discussed previously, for reasons of efficiency, the inversion algorithm EH3DTDinv is best used with a limited number of sources, each with as many data as possible. Secondly, as Macnae (1985) explains, smaller loops tend to couple poorly with deep conductors due to their smaller dipole moment. The smaller loops at San Nicolas tend to exhibit extremely small responses to the San Nicolas deposit. The large loops of the survey detected the deposit, as shown by the initial report by Visser & Shledrake (1999). A more complete interpretation of the data was hampered, at least partly, due to the complexity of the local geology. In particular, the conductive overburden tends to mask the deposit response, making the response of the deposit quite subtle for many of the transmitter-receiver configurations used here. The best example of a line of U T E M data which have clearly detected the deposit is shown in Figure 3.4. The transmitter in this case is a large loop to the south of the deposit, referred to as loop 9. Measurements of were made along North-South oriented lines to the north of 22 Chapter 3. A methodology for three-dimensional TEM inversion -706 a) -1300 -2670 8 8 g * ( D N O n Loop 9 -2225 -1779 Soult] §5 100 j ° -100 o , -200 - 300 * cf 5 • ch 6 0 ch 7 V ch 8 A ch 9 O ch 10 -700 -600 -500 -400 300 -200 -100 northing [m] 5 10 X s 0 |lnductive response] |current channeling | EEiiEIEE j • ch 1 [ii I ] D : C C : 3 : C C : y : i i i""""" ;'' l''*'' r!j* J= -700 -600 -500 -400 -300 -200 -100 northing [m] b) Figure 3.4: a) Showing the location of the transmitter loop 9 along with location of associated receivers, b) U T E M data for the line 1600W are plotted in channel 1 reduced continuously normalized format. Different scales for 3 groups of time channels are used to emphasize features in the data. The location of a deep seated inductive response due to the deposit along with the location of a current channeling response are indicated. These features are based on the original interpretation of the data in the technical report by Visser & Shledrake (1999) the transmitter loop. The data are presented employing a method called channel 1 reduction with continuous normalization. We will discuss this type of data normalization and reduction in detail along with many other types of data plotting subsequently in chapter 4. It is sufficient to note at this point that the data are reduced and normalized to emphasize features which are related to the secondary fields generated by conductors in the ground. The location of the inductive response noted in Figure 3.4 corresponds well with the lateral location of the San Nicolas deposit as can be seen from the geologic cross section shown in Figure 3.2. This type of information, which can be visually extracted from the data is very useful, however it gives little detail as to the 3D geometry of the conductor or its true conductivity. Geophysical inversion provides a method of potentially extracting much more information from the data. In the following document we apply our 3D forward modeling and inversion to these data. We detail some of the challenges encountered and attempt to develop a systematic method 23 Chapter 3. A methodology for three-dimensional TEM inversion for inverting field data which can be applied to time domain surveys in general. 24 C h a p t e r 4 On required inversion input and preparations 4.1 Introduction Inversion is a complex, multistage procedure and as a consequence, it requires a significant amount of preparation of the data and many other input parameters in order to produce useful results. In this chapter we will discuss the preparatory steps required to invert a T E M dataset. These steps include all stages prior to actual inversion that are shown in Figure 4.1. They are: understanding the data, developing a background model, developing a discretization, validation of the discretization, forward modeling and the assignment of errors Throughout we will apply the methodology that is developed here to the archetypical example for 3D T E M inversion, the San Nicolas U T E M survey. 4.2 Understanding the data The first and perhaps the most important stage of the inversion process is to understand the data. Without a clear understanding of the data it is impossible to evaluate modelling and inversion results in any meaningful way. Unfortunately this step happens to be often neglected and this may lead to failure of the inversion or poor results. Individual T E M systems use transmitters with a variety of different configurations and pro-vide data in a variety of different units. T E M systems often provide data with some kind of normalization applied. For inversion with the code EH3DTDinv the data must have all normalization procedures carefully removed to provide E , H, or ^ field values in System In-ternationale (SI) units. Not only must the normalizations be removed but the data must be converted into a coordinate system that is consistent with the 3D modelling domain. Addi-tionally the location of the transmitters and receivers must be precisely known along with the 25 Chapter 4. On required inversion input and preparations understanding the data < • background model j y*~ •* » discretize * | j validate forward model j I •{ error assignment j invert data i evaluate results j Figure 4.1: A simplified procedure for the inversion of 3D time domain electromagnetic data. In this chapter the topics discussed will include: understanding the data, developing a background model, discretization, validation, forward modeling and the assignment of errors specifics of the transmitter current waveform. Since the San Nicolas survey used the U T E M system we will discuss these issues as they relate to U T E M data. This may seem like minute detail, however a similarly detailed analysis must be done for any other type of data set encountered by the geophysicist. 4.2.1 Introduction to the U T E M survey The U T E M system uses a continuous on-time loop transmitter with a sawtooth shaped wave-form. Coil receivers measure the time derivative of the magnetic field in 3 possible components. Measurements of the horizontal components of the electric field are also possible though less common. The system is described in detail in the paper by West et al. (1984). This section will undertake to discuss several aspects important to understanding the U T E M data. These include: the nature of U T E M data, the process required to convert it to a suitable form for inversion with EH3DTDinv, and the specific details of the San Nicolas U T E M survey itself. This discussion will include a discussion of the necessary conditions required to convert U T E M data in ^ to step responses in H. 26 a priori information Chapter 4. On required inversion input and preparations U T E M ^ data are averages of the rate of change of the magnetic field recorded in specific time windows and labeled by a mean delay time. The delay times are referenced to the time of the reversal in polarity of the current derivative in the transmitter current waveform. The windows defined by binary spacing such that the latest time window encompasses the data recorded during last half of the transmitter ramp. This time window is referred to as U T E M channel 1. The remaining early time is divided in two equal parts and the adjacent window to channel 1 is defined as channel 2. Subsequently the remaining window comprising the earliest time data is again divided in two equal parts and the procedure continues to be repeated until ten channels are defined which encompass the entire range of a single ramp or half cycle in the U T E M waveform. The U T E M system thus has ten time channels with the earliest time being labeled, in a somewhat counterintuitive manner, channel 10 and the latest time channel labeled channel 1. An image modified from West et al. (1984) illustrating the positions of last three U T E M time windows with respect to the transmitter current is shown in Figure 4.2. The earliest time windows are omitted from the diagram in the interest of clarity of presentation. U T E M data are traditionally interpreted by viewing profiles plotted after a process of data reduction and normalization are applied to each channel. The data reduction and normalization is designed to have the effect of compensating for spatial and temporal variation in the primary field of the transmitter. This allows the fields from many different locations and time channels to be viewed on a similar scales and emphasizes the anomalies in the recorded fields due to subsurface conductors. This is useful for visual interpretation but as a result U T E M data, as it is generally provided by geophysical contractors, is not in a form suitable for 3D inversion. As there are several types of data reduction and normalization possible, it is essential that the data are well documented to ensure that accurate removal of the reduction and normalization is possible. 4.2.2 U T E M data reduction and normalization at San Nicolas The industry standard is to provide U T E M data in a reduced normalized form. There are several types of normalization and reduction possible. Here we will confine ourselves to a discussion of a subset of the most common types of reduction and normalization, including those that are used at San Nicolas. There are two common types of normalization, point normalization and continuous normal-27 Chapter 4. On required inversion input and preparations TRANSMITTER CURRENT Figure 4.2: Schematic diagram showing locations of U T E M time windows corresponding to channels 1 through 3. Labels C\ through C 3 indicate the location of each time windows' mean delay times. The locations of U T E M channels 4 through 10 are not shown here. The figure is modified from West et al. (1984) 28 Chapter 4. On required inversion input and preparations ization. There are also two common reduction possibilities, channel 1 reduction and primary field reduction. The U T E M dataset at San Nicolas used continuous normalization exclusively however the data were reported using both types of reduction. Continuous normalization means that the primary field is used continuously in a spatial sense in order to normalize the data. Point normalization means that the primary field a single point in space is used to normalize the data. Channel 1 reduced, continuous normalization is defined as follows Cf = ( C * ~ °^ (4.1) Hp(%ji Vj-i Zj) where Cf and Cf are the normalized and the un-normalized data respectively, for the iih time channel at location (XJ,VJ,ZJ). H p ( X J , Vj, Zj) is the primary field of the transmitter in free-space at the point (XJ , yj, Zj) Primary field reduced, continuous normalization is defined as follows $ = iCi-Hp{xj,yi,Zj)) Hp{xj J Vj) Zj) Continuous normalization helps to account for exponential variation in field strength as a function of distance from the transmitter. This means that anomalies will tend to plot on similar scales regardless of the distance from the transmitter. In channel 1 reduction, C\ itself is reduced by the primary field. Channel 1 reduction tends to emphasize the secondary fields due to discrete conductors in the ground since the background effects are effectively removed by subtracting the late time channel. However when a discrete conductor has significant late time response this may have the effect of distorting the early time channels. Primary field reduction, by contrast, ignores the effects due to the background as it removes only the primary field of the transmitter calculated in free space. This leaves the early time channels 2 through 10 undistorted by any possible residual response in C\ and can aid in interpretation in cases where significant variations in channel 1 data are present and where that variation is not attributable to the background conductivity. There are other types of normalization possible. Point normalization normalizes all data by the primary field at a single point. It is often used in the interpretation of thin vertical conductors where preservation of scale has advantages to the interpreter Macnae (1985). 29 Chapter 4. On required inversion input and preparations Conversion of reduced and normalized formats to data in SI units is relatively simple pro-vided the normalization type is known. It is always best to get clear confirmation from the contractor as to which format the data is provided before proceeding with the conversion to SI units and the preparations of the data set outlined in the following sections. We find that proper conversion of data formats, though seemingly a trivial detail, is often one of the major bottlenecks in efficient and accurate inversion of nearly every type of geophysical data set. 4.2.3 U T E M orientation conventions Care must also be taken with respect to orientations of data as geophysical contractors often have different definitions for the coordinate axes to those used in electromagnetic modeling. For the E M modeler the x and y coordinates are usually absolute and often correspond to East and North. Traditionally for U T E M surveys, due to the prevalence of 2D interpretation, the data are provided in a form where the appellations x and y are relative to the recording line direction and imply in-line and cross-line respectively. It must B E made abundantly clear precisely which orientation is which and how those orientations correspond to the coordinate system used in the 3D modelling software. 4.2.4 Synthetic examples of U T E M data The intention here is to demonstrate several of the key features of U T E M data through ID forward modelling. The first and most important feature of U T E M data is the fact that the survey can, in certain situations, be viewed as equivalent to measuring the step response of the magnetic field due to the ground. This is because the U T E M survey uses a receiver which measures the time derivative of the magnetic field. Maxwell's equations in the earth are a linear system. The 100 percent duty cycle square wave shown in figure 4.3b is a scaled version of the time derivative of the sawtooth waveform used in the U T E M system. If coils are used to measure the time derivative of the magnetic field, then it can be shown that this is essentially equivalent to measuring the response of H to a 100 % duty cycle square wave. This is because B and H are related by the constitutive relation in equation 1.4 and the response to the square wave differs only by an additive constant and a scale factor from that of a step off due to a unit current in the transmitter. The recorded data can be converted to an equivalent 30 Chapter 4. On required inversion input and preparations Figure 4.3: a) Sawtooth shaped U T E M current waveform at base frequency of 30.974 Hz. b) 100 % duty cycle square wave. Note that the square wave shown here is a scaled version of the time derivative of the sawtooth wave. step off response using the following equation BsTEP = -T— 2m fdB\ (dB\ V dt J U T E M V dt Jo (4.3) where (^)o is the steady state coil response to a constant rate of change of current, m is the rate of change of current in the transmitter and the quantity ( ^ )T J T E M is a U T E M datum. Equivalently this can expressed in term of the magnetic field as follows HSTEP = i ( ^ ) u T E M " H ° ( 4 - 4 ) where Ho is the free space response in the magnetic field due unit current in the transmitter. A numerical example illustrating this conversion is shown in Figure 4.4. Here the ID forward modelling code EMIDTMfwd Farquharson et al. (2005) is used to model the magnetic field step response to a very resistive halfspace. The ID forward modelling code is also used to generate the U T E M ^ response to the same halfspace. The U T E M data are converted using equation 4.4 to produce the equivalent step response in magnetic field. The two step responses are plotted in Figure 4.4. The two responses are virtually identical and the slight differences are attributable to small numerical errors in the calculations of the two types of fields. It is also useful for the purposes of understanding the data to generate a simple conductivity 31 Chapter 4. On required inversion input and preparations :c 1 1 * — l - l M - H - l M t . .1 . t J.I.I.I. i + UTEM i } • step r r r • • • • ! C[ ; J : a': :!:c i ill:::: :::::rr::y.r. iii •i-rr ; £ : : : : -!-r • .1 c c c : : Ui i i -- j . [ £ " " : : : : : : : : : : : : ^ : : : : ; : : -\ "H" i— Hi' j. I--- -N-i-i-i 10"5 10'4 1D'3 10'2 10'' t ime [sec] Figure 4.4: Showing a comparison of the converted response in H z due to the U T E M system and a step off current source. Conversion was done using equation 4.4. Both responses are modelled for a 20, OOOfim halfspace using the ID forward modelling code EMIDTMfwd. The U T E M system was modelled as operating at a base frequency of 30.974 Hz model and forward model the response of the U T E M system. The response can be plotted in reduced normalized form to illustrate a typical profile over a conductor. This has been done done in Figure 4.5 where the response to a conductive block model is shown in the bottom right panel. The synthetic data resulting from this model are plotted in channel 1 reduced continuously normalized form in the upper right panels. The block has a resistivity of 2.0 Vtm and is buried at a depth of 200 m in a background halfspace with resistivity of 200 Clm. The figure illustrates a few key points. The response to the block is largely confined to the latest few time channels. The earliest time channels are highly affected by the presence of the background halfspace and in this survey configuration the response to both the background and block are not completely dissipated by the latest time, channel 1. The U T E M responses to the block appear large when compared to the responses that were interpreted to be due to the San Nicolas deposit in the recorded U T E M data for this same transmitter configuration. The conductive overburden that is present at the deposit masks the response of the deposit. The overburden is missing from this synthetic model thus explaining the comparatively large responses measured in the synthetic block modelling. Now that some synthetically generated U T E M data have been presented and the data formats explained we move on to review the configuration of the San Nicolas U T E M survey. 32 Chapter 4. On required inversion input and preparations •2300 3000 -ia00 -HDD -l«D -1300 -1000 4900 400 | | northing [m] [ Figure 4.5: a ) Survey configuration showing the location of transmitter loop and a receiver line. Also shown is the outline of the surface projection of the buried conductive block, b) U T E M data for early-time channels plotted in primary field reduced continuously normalized form c ) U T E M data for intermediate channels plotted in primary field reduced continuously normalized form d) U T E M data for channel 1 plotted in primary field reduced continuously normalized form, e ) East west cross-section through synthetic conductivity model used to generate data in the upper panels. 33 Chapter 4. On required inversion input and preparations 4.2.5 The San Nicolas U T E M survey The survey geometry has been shown in Figure 3.3. The U T E M transmitters were set to a base frequency of 30.97 Hz. The mean delay times for each channel in the U T E M system operating at this frequency is provided. At San Nicolas the portion of the U T E M survey considered here consists of 3 large loops: a large loop which surrounds the deposit, a loop south of the deposit and a loop to the east of the deposit. A number of smaller loops were recorded over the San Nicolas deposit however the data recorded from them are not considered in this particular study. Data were recorded in vertically oriented coil sensors on north-south and east-west lines over the deposit. The data were primarily z component ^ however some x component data were recorded on lines collected using southern transmitter, loop 9. Additionally some electric field data were recorded at the survey site for large loop transmitter, loop 2, however these data were not provided by Teck Cominco and thus are not considered in this study. The data were provided in both channel 1 reduced and primary field reduced continuously normalized format. The ^ data were converted to units of volts, or T/s, for inversion with EH3DTDinv. Channel Delay Time [ms] 1 12.11 2 6.053 3 3.027 4 1.513 5 0.757 6 0.378 7 0.189 8 0.095 9 0.047 10 0.024 Table 4.1: Mean delay times for the U T E M system operating at a base frequency of 30.974 Hz 34 Chapter 4. On required inversion input and preparations 4.3 B a c k g r o u n d m o d e l l i n g Determining a background model is an important stage in the inversion process of 3D inversion of tiiue domain data. The background model determines the way in which the transmitter waveform and the spatial domain is discretized. It also helps dictate the minimum size of cells in the central portion of the mesh as well as the number of padding cells required to obtain accurate forward modelling results. This will be discussed in more detail in the sections following on the subject of discretization. The background model is also essential to the use of a method of reducing the forward modelling volume referred to as the corrective source technique. The corrective source technique enables inversion for a much smaller volume by calculating a set of artificial sources which reproduce, on a much smaller mesh, the fields forward modelled on the larger volume. The topic of corrective sources will be discussed in detail in section 4.4.3. For the purposes of this discussion it suffices to say that in order to perform the preparatory steps for inversion it is necessary to determine some kind of background model which defines the gross characteristics of the survey area. We believe that the use of the simplest possible background model is advantageous even in situations where the background is known to have complicated ID, 2D or even 3D structure. The reason for this is because, though complicated subsurface structures may be known to exist, it is never completely clear where the boundaries of such objects actually lie, even when detailed a priori information exists. The resulting output inversions can become biased in an unpredictable way if a boundary or object is, for any reason, placed in the incorrect location in the background model. The correct place to use detailed a priori information is in the reference model, not the background model. The exception to this is situations where significant topography or bathymetry is present. We believe that this information should be incorporated into the background model. This means that a topographic halfspace of a single conductivity is generally the best choice for a background model, especially for initial inversions. As such the subject of the following discussion is how to determine the "best" halfspace to use for any particular survey. We envision three distinct methods for determining a background model. First a priori information about the host rock in a targeted area can be used. In the second method one dimensional forward modelling can be used to predict the data from various halfspaces. The predicted data from this can then be compared to the recorded data in order to determine a 35 Chapter 4. On required inversion input and preparations background model that fits the data in an approximate way. The third method is to use one of a number of existing I D T E M inversion codes in order to determine the best fit halfspace or simple layered space. 4.3.1 Background models from a priori information A priori information can often provide all the necessary information for background modelling, removing the necessity for inversion and taking much of the aspect of trial and error out of forward modelling. For example, if previous surveys have been completed or drilling has been done in the area, then an abundance of a priori information may be available to the geophysicist. This type of information often presents its own difficulties. Conductivities obtained from drill cores tend to display a bias associated with the scale dependence of electrical conductivity. Thoughtful consideration of this fact must be made in order to reconcile the two types of data and produce a reasonable background model. In the case of San Nicolas, the conductivity of the host rock was well determined due to the comprehensive work done by Phillips (2001) including, not only extensive interpretation of drill core conductivity testing but also the inversion of numerous E M surveys at the site. Table 3.3.2 illustrates this fact. It is, however, recognized that this abundance of a priori information about lithology and their associated physical property values may not always be present. In this case it is still possible to obtain a background model by the process of I D forward modelling or inversion. The San Nicolas background model at the most extreme level of simplification, based on a priori information, was determined to be a halfspace of roughly 100-200 Qm. This, as it will be shown, is consistent with both forward modelling and inversion results used to determine the background model. Note that the conductive overburden at San Nicolas, which is known to exist at the site, has been neglected intentionally in this model. 4.3.2 Background models from ID forward modelling As discussed previously we suggest that, initially, the simplest background possible should be used. The method of forward modelling is useful for developing background models in situations where little or no a priori information is available. The methodology involves conducting a series of forward modellings with halfspaces of differing conductivity. The method then proceeds to 36 Chapter 4. On required inversion input and preparations make a comparison of the predicted data with the observed data. There are several advantages to using forward modelling as a method of determining a background conductivity. First the process of forward modelling is useful in identifying anomalous data, as halfspaces tend to produce smooth curves, of very specific shape. Deviations from this shape in the observed data can indicate the locations of anomalies in the data due to objects of interest. Secondly the modelling of T E M responses can often be useful in locating errors in the conversion of data formats including ensuring that polarity of the modelled fields matches the recorded data. Additionally the process of forward modelling halfspace responses and comparison with the data is insightful. It can lend a hand in estimating data quality especially in situations where no information about measurement errors is available. Often early time results are affected by complexities in the near surface which we do not wish to include in the background model. For most common situations it is important to match the later times more closely than the early times. Figure 4.6b shows the observed C2 data for loop 2 and 4.6c shows the forward modeled responses to a 100 flm halfspace. The data, particularly at late time, shows a reasonably good agreement with the data recorded at the San Nicolas deposit. A direct comparison showing decay curves for the U T E M data and the forward modelled data is shown in shown in Figure 4.6d . 4.3.3 Background models from I D inversion If no a priori information about the background model is available then inversion can also be used, in concert with forward modelling, to develop a background model. The technique that we envision is the use of averages based on a large group of ID models developed with ID inversion. This large group may have statistical properties which reflect the behavior of the background model as a whole. The method of inversion has some advantages over that of forward modelling as it is an automated procedure and avoids the potentially time consuming trial and error inherent in forward modelling. However the process of inversion is complicated, and as such, it has the potential to produce results which are misleading if used incorrectly. In the example of the San Nicolas deposit we use the ID inversion algorithm E M 1 D T M of Farquharson et al. (2005) to invert 315 soundings from the largest loop shown in Figure 4.6. The results have been collated to produce an average conductivity as a function of depth. It is 37 Chapter 4. On required inversion input and preparations Figure 4.6: a) Showing the average background model as a function of depth based on 315 separate ID inversions of U T E M soundings taken at San Nicolas. Curves also show the location of 1 standard deviation from the average, b) Observed U T E M Ci data from loop 2. c) Forward modelled U T E M C2 data from loop 2 using a 100 fim halfspace. d) Comparison of decay curves showing predicted data for a 100 halfspace and observed U T E M data. Both curves plotted for the location [-1500, -556.633,1]. 38 Chapter 4. On required inversion input and preparations clear from the results that the recovered conductivity at San Nicolas is highly variable but the average conductivity shows some interesting details including the fact that the surface is much more conductive than the deeper material. This seems to confirm the presence of conductive overburden at the site. The conductance, the conductivity-thickness product, is a key parameter for electromag-netics. We use the concept of the conductance to calculate the average halfspace conductivity through the use of the weighted averaging equation - = E U i <nu ( 4 5 ) Z^n=l bi where o is the average halfspace conductivity, <7j is the average conductivity the ith layer and ti is the thickness of the ith layer. For the collection of I D inverted soundings plotted in Figure 4.6 this equation gives an average conductivity of 0.017 S/m. This corresponds to an average resistivity of roughly 59f2m. Additionally a quick visual inspection of Figure 4.6 reveals that a value of 100 Qm lies within the limits of one standard deviation of the average model at every depth. It is clear that the I D inversion analysis, to a large degree, agrees with the previous infor-mation acquired from both forward modelling and a priori information. For the San Nicolas survey, a background model of 100 Qm was used to develop a discretization. This decision was consistent with all three of the techniques discussed here, a priori information, forward modelling and I D inversion. 4.4 D i s c r e t i z a t i o n The importance of discretization on inversion results cannot be underestimated. The process of discretization in T E M inversion is more difficult than in other electrical and E M methods such as D C resistivity and frequency domain E M because the domain must be discretized in both time and space. This leads to complications as these discretizations are intimately related to each other and changes in time discretizaton must be accommodated by changes in spatial discretization and vice versa. Here we attempt to outline a procedure whereby an adequate discretization can be achieved for both space and time. 39 Chapter 4. On required inversion input and preparations AAA Discretization in time A background model assists with development of a discretization in time. The time discretiza-tion should ensure that the accuracy of the forward modelling is such that any errors due to modelling are below the estimated noise level for the survey, or at a minimum, well below the expected response to the target. There are several key issues one must address in order to dis-cretize the transmitter waveform in time. The discretization must contain an adequate density of sampling and the discretizion must begin an appropriate amount of time before the first data are sampled. These issues are discussed at length in Phillips & Oldenburg (2005) as part of UBC-GIF's T I M E project consortium report. The report suggests that in order to maintain a modelling error below 2.0% we must sample fields at least 2 decades before the first modelled channel and maintain a sampling density of at least 4 samples per decade. For the example of the San Nicolas deposit U T E M survey, which uses an on-time system, one must be especially careful. Not only does adequate time sampling have to be maintained but there is the potential for subsequent cycles of the transmitter waveform to interfere with each other. Ideally we wish to convert U T E M responses to step responses in H as this necessitates far fewer points in the discretization than would a continuous on-time waveform. In order to determine whether this is possible we must determine if modelled electromagnetic fields reach an equilibrium within the earth within a single half cycle, or ramp, of the U T E M waveform. Ideally step-off conversion should be undertaken only if it can be shown that there is no significant response in the model which outlasts the duration of a single ramp in the U T E M system. In order to determine whether this may be the case at San Nicolas we propose some simple modelling exercises. In a very approximate way the nature of the secondary fields over the San Nicolas deposit will be the superposition of two types of fields. The first type will be those due to a background layered space or halfspace while the second type of fields will be those generated by the conductive target. The first goal is to determine whether the response in the magnetic field time derivative for the background model reaches a steady state during a single ramp of the U T E M transmitter waveform. The second issue is to determine if the magnetic field response of a target with similar conductivity and size to that of the San Nicolas target will be significant at late times. We begin simply with ID analytic and semi-analytic modelling for both halfspaces and conductive spheres respectively. Figure 4.7 shows the step-off response in H for the San Nicolas 40 Chapter 4. On required inversion input and preparations —ha l f space 0 30.974 Hz UTEM channels time [s] Figure 4.7: Free space time domain response of a sphere buried at 200 m depth with a resistivity of 2.0 Jim and radius of 150 m (dashed lined) for a large circular loop transmitter. Also shown (solid line) is the response of a 100 flm halfspace. The transmitter is a circular loop with a radius of 500 m and responses are measured at the center of the loop directly above the sphere. Diamonds show the location of the 10 U T E M time channels. background halfspace model following Ward &; Hohmann (1987). Also shown is the step-off response in H to a large conductive sphere which approximates the San Nicolas deposit in both size and average estimated conductivity. The sphere modelling is based on the formulation presented by Wait &: Spies (1969) along with a discrete implementation of a frequency-to-time conversion. The survey configuration is loop transmitter which is a circular approximation of the situation at San Nicolas for transmitter loop 2. The figure shows that while the step response for a large spherical conductive object drops off rapidly at late times the response to the background halfspace does not. This indicates that, while the deposit does not have a significant effect on subsequent cycles of the waveform, the conductive host rocks may have a significant impact. This makes the feasibility of a step-off conversion of the U T E M data uncertain. To quantify the interaction of the halfspace background model with the U T E M fields from previous cycles of the U T E M waveform we compare two situations. First we calculate the steady state response in H z due to a unit current in the transmitter over the background model. Second we calculate the value of H z over the halfspace background model at the end of a series of 4.25 full transmitter cycles with a rectangular 100 % duty cycle waveform. Figure 4.8 shows that we find that there is roughly a 4% difference in the calculated fields. While this amount is 41 Chapter 4. On required inversion input and preparations -1678 -1402 -1126 easting [m] 7.91e-004 7.49e-004 7.06e-004 6.63e-004 6.21e-004 5.78e-004 5.35e-004 F's] -1678 -1402 easting [m] -1126 5.35e-004 F's] -1678 -1402 -1126 easting [m] 4.00e+000 3.74e+000 3.49e+000 3.23e+000 2.97e+000 2.71e+000 2.46e+000 F/s] d) 7.40e-004 F/s] 5 15e-004 Figure 4.8: a) H z Steady-state response to San Nicolas transmitter loop 2 over a 100 tlm halfspace b) A map showing percent difference in the responses between steady-state and the square wave c) Map of H Z response after 4.25 cycles of a 100 % duty cycle square wave, d) profile through points [A,B]. Blue is the response due to steady-state initial conditions while green is the response after multiple cycles of the square wave. not overwhelmingly large, and probably represents an acceptable modelling error for many E M targets, consider the magnitude of the anomaly identified in the original technical report by Visser & Shledrake (1999). The report indicates that the magnitude of the San Nicolas anomaly is less than 5% of the primary field. The conclusion here is that introducing a 4% modelling error due to step off conversion is unacceptable. This leaves two avenues of pursuit for time discretization. It is possible to modify the initial conditions of the step off to match those achieved after extended cycles of the U T E M waveform. This however is a model dependent correction and would need to be re-calculated if the model changes in a significant way over the course of the inversion. The alternative is to develop a single waveform using the true sawtooth shaped current in the U T E M system and calculate ^ at data locations for the model in question. This second method exactly reproduces the true U T E M data and does not require a transformation of the data to the analogous situation using the magnetic field. This method turns out to be a far 42 Chapter 4. On required inversion input and preparations Figure 4.9: The time discretization of the U T E M waveform to used for inversion. simpler solution in application and for this reason it was the method chosen for this study. A synthetic modeling exercise with multiple cycles of the true, sawtooth transmitter wave-form was carried out using the background model developed previously. We find, for the error level ascribed to the data, that data after | of a U T E M cycle are sufficient to achieve an approx-imate equilibrium to within less than 1 % error. Note however, that in general, the modelling error is a function of the background conductivity and, of course, the data sampling times used in the particular survey and this issue should be re-investigated for each new survey. The initial inversions of San Nicolas data used a discretization of the U T E M waveform with 38 time steps over lj cycles of the sawtooth wave. The discretization is shown in Figure 4.9 4.4.2 Discretization in space Discretization in space necessitates two considerations. First, mesh boundaries must be suffi-ciently distant from the survey area so that the boundary conditions of the forward modeling problem are approximately satisfied. This means that the fields should be, in some sense, "close to" zero at the boundary. How close these fields actually must be to zero depends entirely on how accurate the solution is required to be. Second, the cells in the center of the mesh should be fine enough to describe the early time data adequately. Some work on the subject of 3D spatial discretization for T E M problems is presented in Napier & Oldenburg (2005). The key elements of this work indicate that ideally we should let the mesh boundary extend away from the survey location to several diffusion 43 Chapter 4. On required inversion input and preparations distances, based on the latest time recorded. The size of the smallest cells in the mesh are required to be a fraction of the diffusion distance of the earliest time that is to be used in the inversion. We define the diffusion distance as being the quantity L in equation 4.6 as defined by Ward & Hohmann (1987). here the conductivity, o is that of the background model and t is the time channel in question. Following these recommendations the forward modelling mesh for San Nicolas used cells with a length of 25 m on each edge in the central region and had a mesh boundary that was roughly 8 km from the transmitter loop. The forward modelling mesh consisted of 331,000 cells. This is a large mesh and computation of an accurate forward model on it can take several hours. Because inversion requires many forward models to be computed for each iteration, simply following these mesh design recommendations to the letter often leads to a mesh which is far too large for an inversion to be completed in a reasonable amount of time. Even provided unlimited time in which to perform the necessary forward models, these large meshes can lead to an inversion which makes R A M memory demands which cannot be satisfied by present day computers. We discuss how to reduce the volume that must be modelled in any inversion in the next section 4.4.3 Reducing the volume It is often necessary to reduce the volume of interest by removing padding cells in the inversion domain. The technique used to reduce the volume for inversion is one whereby we can use an accurate forward solution to the E M problem on the background model. This solution is used in order to assist with forward modelling on the smaller inversion domain. This accurate solution may come from analytic solutions, ID forward modelling or 3D forward modelling on a large and finely discretized mesh. In the corrective source technique the fields due to some simple background are used to calculate corrective sources which reproduce, on the small inversion domain, the fields calculated using the accurate solution. This is closely related to the primary secondary formulation described in the paper by Haber et al. (2007). Recall from chapter 2 that the forward problem can be described by a linear system of equations. Given a background model mo, along with sources and boundary conditions q the (4.6) 44 Chapter 4. On required inversion input and preparations accurate fields u produce a corrective source, s c when modelling on the reduced domain. The inversion can then be solved on the reduced domain using the corrective source to compensate for errors due to the boundary conditions being imperfectly satisfied. The equation that must be solved for forward modelling is now This produces adequate results provided that the new model, m, does not differ from the background model, mo, by too much. The practicalities of using the corrective source technique must be considered here. In order to achieve high accuracy in forward modelling the background model must be a reasonable approximation of the true conductivity and the conductivity of the inversion model, at any stage, must not differ by too much from that of the background. One issue that may be a concern is the presence of conductors in the inversion model which have strong secondary fields associated with them. The boundary of the reduced volume must be far enough from these conductors such that the secondary fields are essentially zero on the boundary of the reduced volume. Generally this does not present a problem but in the case of extreme conductivity contrast, such as is common in nickel exploration, this may become a serious issue. This can be investigated on a case by case basis by the use of the same kinds of methods employed in the previous section on time discretization. Where conductors are discrete and either plate-like or roughly spherical a number of forward modeling codes exist which are capable of calculating the secondary fields of these types of objects. These fields can then be interrogated to determine their magnitude on the boundary of the inversion domain and these magnitudes can be compared against that of the features of interest that are being fit in the observed data. The methodology is exactly that used to develop mesh design criteria in Napier & Oldenburg (2005) but extended to secondary rather than primary fields. For the San Nicolas U T E M survey the conductivity contrast is not large, thus secondary fields do not present a significant problem in the volume reduction. For the San Nicolas inversion the reduction of volume technique allowed the inversion domain to be reduced from a forward modelling domain of 313,000 cells to roughly 138, 000 cells. This enabled substantial time A ( m 0 ) u - q = s c (4.7) A ( m ) u = q + s c (4.8) 45 Chapter 4. On required inversion input and preparations savings in the inversion process and kept memory demands for the computation within limits dictated by the presently existing hardware. 4.5 U n d e r s t a n d i n g the d a t a t h r o u g h forward m o d e l l i n g Forward modelling serves the important purpose of giving the geophysicist familiarity with, and hopefully intuition about, the data. This can be invaluable in the interpretation of both the data itself and the inversion results. Equally important is the use of forward modelling as a check on the accuracy of the discretization. Forward modelling also has important uses at this stage of inversion preparation, especially if a priori information about the survey area or suspected target exists. It can help locate anomalies in the data due to the target and this can help determine a level of both modeling and measurement error that is tolerable. This means that forward modelling is important not only in survey design but also in the evaluation of previously acquired data. This step is often neglected in the process of inverting geophysical data, and in doing so, much diagnostic information about the ability and degree to which inversion is able to succeed is lost. 4.5.1 Validation of discretization Validation of discretization is a key stage in the preparation for inversion. It is necessary to determine if the space and time discretizations are accurate enough to model the fields to an acceptable level of error. What is meant by an acceptable level of error will be discussed in detail in the next section. It suffices to say, for the purposes of this discussion, that we must ensure that the discretization in both time and space reproduces the results acquired through an independent solution for the E M fields, given some pre-defined level of error. Halfspaces make excellent test cases for validation as analytic solutions for the electric and magnetic fields exist, at least for many common transmitter current waveforms and geometrical configurations. Another method which is useful to perform validation is to use numerical codes which are computationally distinct from the method used in E H 3 D T D . Integral equation codes or ID forward modelling codes which use a frequency-to-time conversion may be useful for these purposes. For the San Nicolas deposit the discretization was checked using a numerical solution from a ID time domain forward modelling code EMIDTMfwd written by Farquharson et al. (2005). 46 Chapter 4. On required inversion input and preparations Figure 4.10: A validation of the time and space discretization developed for use for the inversion of the San Nicolas U T E M data set. are calculated over a 200 Qm halfspace using both E H 3 D T D and the independent ID forward modeling code EMIDTMfwd. 47 Chapter 4. On required inversion input and preparations The code employs propagator matrices and Shelknoff potentials and as a result does not suffer from spatial discretization error. It does, however, use a frequency-to-time conversion which is dependant on a discretization of frequency space. Since the ID forward computation is fast, enough frequencies can be used so as to make the solution exact, for all intents and purposes. Figure 4.10 shows a comparison between the value of calculated using E H 3 D T D and the value calculated using EMIDTMfwd using the time and space discretization outlined above. Note that the values are visually indistinguishable, and for no U T E M channel do they differ by more than 1.5% from each other. 4.5 .2 Forward modelling using a priori information If a priori information exists then forward modelling can provide substantial insight into the data which have been recorded. Forward modelling can also provide useful input for error assignment as we shall see in the next section. Here we undertake to describe forward modelling in two situations: the first where very little a priori information is available and the other where a substantial amount of information is available about the geometry and conductivity of the target. Both of these examples relate to the San Nicolas deposit however the methodology is generic and can be used for any T E M inversion problem. The situation where little a priori information is available is the conductive block scenario discussed in the next section. The situation where detailed a priori information is available is discussed in the following 2 examples of forward modelling. A conduct ive block We wish to forward model the response of a simple synthetic model based upon the San Nicolas survey geometry and some very simple a priori information we have about the deposit. Later we will invert this forward modelled data in order to determine the resolving capability and depth of investigation of the 3D inversion given the survey parameters dictated by the San Nicolas U T E M survey. The synthetic model consists of a 2.0Qm block buried at 200m depth emplaced in a 200 Vtm halfspace background. This model represents San Nicolas at the most extreme level of simplification, the model does differ in a key way in that the block extends to a depth much deeper than the true deposit. By doing this we are able to extract information about the depth of investigation in the inversion of the San Nicolas U T E M survey. This will 48 Chapter 4. On required inversion input and preparations become apparent in the next chapter. The data for this model were generated by forward modelling for the loop 2 transmitter configuration shown in Figure 4.11. The model, along with a selection of the data generated along a profile passing directly above the conductor, are also shown. The anomaly due to the block in the halfspace is less than 4 % of the halfspace field for all time channels with only one exception. The maximum anomaly in the data due to the conductive block is on the order of 20% of the value of the field due to background halfspace. This relatively large anomaly exists for a single time channel where the background halfspace field is small. These responses correspond to the approximate strength of the response attributed to the San Nicolas deposit in the original technical report by Visser & Shledrake (1999). Forward modelling of the San Nicolas deposit Much work has been done by Phillips (2001) on the San Nicolas deposit to develop a 3D conductivity model using drilling information. Here we use this 3D conductivity model to generate some synthetic responses to the San Nicolas deposit. A number of different forward modelling exercises may be undertaken here. The first exercise that may be of interest is to compare the response of the forward modelled data using the 3D model with the both deposit and overburden present and compare it to the response generated with only the overburden present. This should help provide information about the relative strengths of the response to the two features. The two models used in the forward modelling exercise along with the configuration of the transmitter loop 9 relative to the deposit and the location of a line profile at 1600 m west are shown in Figure 4.12. The forward modelled data are shown in Figure 4.13 The original report by Visser & Shledrake (1999) indicates that there is a roughly 5% re-sponse in channel 1 data to the deposit for the actual U T E M data recorded in this configuration. Figure 4.13b shows a decay curve directly above the deposit while figure 4.13d shows a decay curve that is slightly south of the deposit. Both of these figures show that the difference in late time channels due to the 3D conductivity model is exceedingly small. In fact the channel 1 response of the deposit is much less than 1 % of the field due to the conductive overburden and the resistive basement layer. There is less than a 5 % response in channel 2 and channel 3 for the 3D deposit model. From this lack of late time response we infer that the conductivity of the target is probably underestimated in the drill core conductivity based model. This situation 49 Chapter 4. On required inversion input and preparations a) i -363 -1490 -544 -2131 South Figure 4.11: a) 3D model showing location of conductive block, transmitter loop and location of data profile [A,B]. The conductivity of the background halfspace is 200 fim b) Map of U T E M channel 4 data for model with block c) Profile of data through locations [A,B] with and without the conductive block d) Map of U T E M channel 4 data for 200 fim halfspace e) Decay curve for location at [-1659,-500,1] showing forward modelled data with and without the conductive block 50 Chapter 4. On required inversion input and preparations Figure 4.12: a) 3D conductivity model showing the San Nicolas deposit and the overburden from analysis of drill core, b) 3D conductivity model showing the San Nicolas overburden from analysis of drill core. The deposit has been removed from this model and replaced with a background halfspace of 100 Qm. In both figures the transmitter loop is shown along with the location of a profile where data was collected is as expected since the drill core conductivities suffer from a downward bias due to the scale dependance of electrical conductivity. This effect was discussed in chapter 1. Despite the fact that the conductivities in the 3D model are likely to be underestimated it may be interesting to compare the real recorded U T E M data at San Nicolas with the fields modelled from the 3D conductivity model. The 3D conductivity model was used to forward model the values of for the locations where data were recorded at San Nicolas. Once again the forward modelling was done for transmitter loop 9 since this transmitter shows the clearest response to the deposit. The forward modelled data are compared directly with the data recorded at San Nicolas and are presented in Figure 4.14. Decay curves are shown in 4.14b and 4.14d for locations directly above the deposit and to the north of the deposit respectively. The figure indicates that while the San Nicolas U T E M data recorded over the deposit are quite similar at early times to the forward modelled data, there are large differences between the modelled data and the observed data in later channels, including channel 2, in locations proximal to the deposit. It has not been possible to show here due to the large quantity of data considered but significant differences exist over all time ranges in other locations away from the deposit, particularly to the north-west, south and north-east of the deposit. The similarities in the early time data in the vicinity of the deposit imply that the 3D conductivity model is a reasonable approximation of the true conductivity in the near surface, at 51 Chapter 4. On required inversion input and preparations • 1800 .1725 fo^ 1 18e-008. 6.99e-009. 5 87e-008. 3 91e-008. 1 95e-008. [T/s] -1.26e-010. -1.97e-0q8. -5 deposit and overburden overburden only 4 -3 log t [s] deposit and overburden overburden only -4 -3 -2 log t [s] deposit and overburden overburden only 9 80.008 rj/s] . S M . O O S B Figure 4.13: a) Map of forward modeled U T E M channel 2 data with overburden and deposit present, b) Decay curves for the point [-1600,-209,1] c) Map of forward modeled U T E M channel 2 data with overburden only and deposit replaced by a background conductivity, d) Decay curves for the point [-1600,-500,1] e) Profile spanning the line [A,B] showing response with and without the deposit 52 Chapter 4. On required inversion input and preparations least in the vicinity of the deposit. Because of the differences between the forward modelled data and the U T E M data at late times, the 3D conductivity model is thought to be an underestimate of the true conductivity of deposit. The source of the differences in the data at all time scales in areas away from the deposit is thought to be largely due to variations in the overburden conductivity. This conductivity variation is not reflected in the 3D model which is based more on physical property averages than actual conductivity data. Also the physical property model outside the immediate deposit area is not reliable as it is based on sparsely sampled data with most drill holes being spotted only in the immediate vicinity of the deposit. The differences between the two datasets outside the immediate vicinity of the deposit is really highlighting a spatial variability in the overburden thickness or conductivity. The comparison between the data generated from the 3D model and the true data indicates a few important points about the overburden and deposit and these are details that should help to evaluate the results of an inversion. Interpreting the various forward modelling results presented here it was possible to determine the following The overburden is variable in either thickness or conductivity The overburden is thick or more conductive over the center of the deposit and continues to thicken or become more conductive to the north-west The overburden thins or becomes more resistive to the south, south-west and possibly the north-east of the deposit The true deposit response appears to be much stronger than the 3D conductivity model indi-cates. This final point is important since it indicates that the true deposit is probably more con-ductive than previously thought. 4.6 Assigning standard deviations to the Data The misfit functional used in the inversion algorithm's objective function, as shown in equation 2.7, requires that each datum be normalized by an estimated standard deviation. This is important in order to ensure that no datum contributes unduly to the solution, especially if the datum is thought to be suspect, or if there is variation in the measurement error associated 53 Chapter 4. On required inversion input and preparations 308*006 •457 • -1800 -1725 -1650 -(ps -1500 7 74e-009 1 SfcQQS b) rj/s] 4.62e-008 decay curve [-1600,-500,1] UTEM data 3D o model 4 -3 -2 log t [s] 3 G M O S 3 oe*ooe -457' -1800 -1725 -1650 .4575 -1500 7 74e-009 1 97e-009 d) rr/s] decay curve [-1600,-57,1] 2 50e-008 5 786-009". 1 93e-O08 3 19e-0O9. rj/s] 5 92e-010. 1.35e-00e 200e4XI9. 460e-0Q9 -5 UTEM data 3D a model -4 -3 -2 log t [s] Figure 4.14: a) Map of observed U T E M channel 2 data recorded at San Nicolas for transmitter loop 9. b) Decay curve at location [-1600,-57,1] c) Map of predicted channel 2 data forward modelled from a 3D conductivity model for transmitter loop 9. d) Decay curve at location [-1600,-500,1] e) A profile through the points [A,B] showing San Nicolas U T E M data and forwaM modelled data for U T E M channel 2 Chapter 4. On required inversion input and preparations with different locations or different time channels. This is another key stage in the inversion process. Unrealistic error estimates can degrade the solution sometimes causing the inversion program to fail long before it reaches a geologically interpretable result. The factors contributing to the total error in a geophysical survey include instrumental bias errors, additive noise sources, mis-locations or mis-orientations of receivers, and numerical errors associated with discretization and modelling. These aspects arise in all inverse problems but are particularly difficult for time domain E M problems where we have many orders of magnitude change in the amplitude of the data as a function of both distance from the source and over differing time channels. It is important not to exclusively use the stacking error which is calculated and stored by some T E M instruments. While the stacking error can represent the background noise level for a particular datum, it does not account for errors due to mis-location, mis-orientation or discretization in the forward model. In practice the stacking error should be used as a minimum floor value of error along with reasonable estimates of the other types of error. It is important to consider the fact that modeling errors become large if a datum is close to a source. This is because the forward modelling code is incapable of a completely accurate representation of a singularity. The forward modelling code uses distributed sources over cell faces and thus data must be far enough from the source that the distributed source appears as an accurate representation of a point or line source. In practice it is desirable to remove any data which are within 5-6 cells of a source. If this is not done it invariably leads to poor results and slow convergence of the algorithm due to the modelling problems associated with singularities. For San Nicolas, the assignment of data error was complicated by the presence of many zero crossings in the data, in both time and space. Data standard deviations were assigned by separating time channels, and giving the data in each time channel an error estimate based on a percentage of the datum and a constant floor value. In consideration of the modelling error caused by having large cell sizes relative to the diffusion distance at earliest times, the percentage error was larger for early time channels and smaller for late time channels. The floor value for each time channel was a constant and based on a fraction of the median of the absolute value of all the data in each time channel. This procedure results in a fairly effective and robust method of assigning standard deviations to field T E M data. Additionally the data 55 Chapter 4. On required inversion input and preparations which were located within 6 cells of the U T E M transmitter loop were removed from the inverted dataset. This ensures that the spatial distribution of modelling error caused by an inaccurate representation of a point or line source does not unduly affect results. 4.7 Conclusion It is useful before moving on to the subject of inversion of T E M data to briefly summarize the accomplishments made in this chapter. The U T E M data type has been introduced and explained in detail. Various methods of determining a background model representative of an entire survey area have been discussed. A procedure employing the background model in order to develop a time discretization of the transmitter current waveform has been presented along with a method to develop an accurate forward modelling mesh. The background conductivity has been used in a technique to reduce the volume of interest for the purposes of inversion. The dicretizations developed have been validated by independent means of calculating T E M fields. The data from a priori models of the survey target have been forward modelled and the implications of the forward modelling results have been considered. It is only at this stage that we are finally ready to invert the data. The following chapter will discuss the topic of 3D T E M inversion in general terms and present the result of the inversion of the San Nicolas U T E M survey. 56 Chapter 5 On practical three-dimensional T E M inversion 5.1 Introduction The previous chapter exclusively dealt with the preparations necessary for inversion. This chap-ter will discuss the process of inversion itself; however, constant reference to topics in chapter 4 will be made. It is worthwhile at this point, to once again quickly review the inversion proce-dure reproduced here in Figure 5.1. This schematic outline of the 3D T E M inversion procedure shows that the process of inversion produces results which feed back into an understanding of the data. In fact the results attained from an inversion, and particularly initial inversions, pro-vide additional information about every other stage in the inversion process. What this means in practice is that the process of completing an inversion and evaluating the results will have an impact on subsequent choices of background model, discretization and error assignment. Fur-thermore any changes in these quantities may require further changes to be made and require more forward modelling to be completed. In this view inversion is an iterative process and each successive result inspires further refinements until, eventually, a result is attained which satisfies the interpreter. The concept of inversion as a process dependant on iterative feedback should become clear to the reader as we proceed with the inversion of the archetypical T E M field survey example, the San Nicolas U T E M survey. This example required several iterations of the inversion procedure outlined in Figure 5.1 in order to produce useful results. 5.2 Inversion of synthetic data It is often useful to invert synthetic data before proceeding to invert the field data. This gives a clear understanding of the resolution of the inversion process for a given noise level, and survey 57 Chapter 5. On practical three-dimensional TEM inversion understanding the data background model a prion information discretize forward model error assignment invert data evaluate results validate Figure 5.1: A simplified procedure for the inversion of 3D time domain electromagnetic data. In this chapter the topics discussed will include inversion, evaluation of results, and the inclusion of a priori information. The iterative nature of the inversion process will become apparent as successive inversions results are presented. 58 Chapter 5. On practical three-dimensional TEM inversion Figure 5.2: a) Observed UTEM channel 2 data for the conductive block model with 2 % gaussian random noise added b) Predicted UTEM channel 2 data for inversion of conductive block data configuration. The synthetic data should be generated by creating a model which typifies the gross characteristics of the exploration target. The data for this model should be generated to high degree of accuracy using a finely discretized mesh with adequate padding such as the mesh one would use to generate primary fields for the reduction of volume technique. It is not a major concern if forward modelling takes several hours to complete since only a single forward solution is required for these purposes. Once the data are generated they should be corrupted with random gaussian noise of a magnitude of some small percentage of the datum in question. The actual magnitude of the synthetic noise may vary from case to case and should depend on an estimate of the measurement error present in the actual data. Standard deviations assigned in the inversion of the synthetic data should use this percentage noise as a guide but also include a small floor value to help deal with situations where zeros crossings occur in the data. Inversion should take place on a reduced volume in order to achieve results quickly. The following section deals with the inversion of synthetic data generated for the San Nicolas deposit and should help illustrate the application of the ideas presented here. 5.2.1 Inversion for a conductive block In the previous chapter a very simple model of a conductive block in a background halfspace was developed based on the structures and associated conductivities known to exist at the San Nicolas deposit. The data for a simple rectangular conductive block were forward modeled using the survey parameters for the large transmitter loop of the San Nicolas survey. The 59 Chapter 5. On practical three-dimensional TEM inversion results of this forward model l ing have been presented previously in a reduced normalized form in Figure 4.5 and have also been compared to a halfspace response in Figure 4.11. These da ta were corrupted wi th 2.0 % Gaussian random noise. T h e data were inverted using an assigned standard deviation of 2 .0% along wi th a small floor value. The reference model and in i t i a l model used were a 200 Qm halfspace. T h e inversion took place on a reduced mesh consisting of 142,002 cells and wi th waveform consisting of 38 t ime steps. T h e objective function min imized was min = — 2 2 W d ( Q A ( m ) - J q - d o b s ) + ^R(m - m r e f) (5.1) 2 2 where the regularization function was defined as R ( m - m r e f ) = a s | | W s W t ( m - m r e f ) | | 2 + a x | | W t W x ( m - m r e f ) | | 2 + (5.2) ay | | W t W y ( m - m r e f ) | | 2 + az | |W t W z (m - m r e f ) | | 2 with the value of as of 1 x 10~ 5 and ax,ay,az of 1. T h e matrices W s , W x , W y , W z are as defined in Chapter 2 and the user defined weighting mat r ix W t was left as the identity. These are the default conditions for the code E H 3 D T D i n v and have been found to produce reasonable results in previous inversions of synthetically generated E M data. Convergence required 39 total iterations over 10 values of 6. T h e cooling of the objective function proceeded from an in i t ia l value of 6 of 10 and the value of Q was reduced by a factor of 5 unt i l reaching a final value of 5.12 x 10~ 5 , at which point the target misfit criteria, x 2 = lj was met. T h e value of the data misfit, ^ as a function of iteration is plotted in Figure 5.3. T h e predicted and observed data are presented in Figure 5.2. The inversion converged to a final model wi th a value of the model objective of 8, 300. T h e true value of the model objective for the actual block model used to generate the data was slightly over 386,000 confirming that, as is necessary for this algori thm, the inversion is producing a model w i th less structure than the true model . T h e final result of this inversion is shown in Figure 5.4. A profile through the true conduc-t iv i ty model along wi th the same profile through the inverted result are provided here, along wi th a map view of the surface of the inverted result. It is worth noting is that the block is blurred and has a lower conduct ivi ty than the true model however this is as expected given the model norm-used here. T h i s inversion of the 60 Chapter 5. On practical three-dimensional TEM inversion T3 -e-15 20 25 30 iteration 35 40 Figure 5.3: Data misfit as a function of iteration in the inversion of the synthetic block model. The total number of data was 2640. The target misfit was 2640. The inversion achieved the target misfit in 39 iterations. -ZOB1 EBUtfl - u w 400 243 147 69 4 54 3 ' 20 -1490 -899 Figure 5.4: a) East-west cross section though conductive block model at 450 m south. The view is from the south, b) Map view of recovered model for conductive block inversion at depth of Om showing the locations of data (red points) and transmitter loop (red lines), c) East-west cross section though recovered model at 450 m south. The view is from the south. The outline of the conductive block is superposed. 61 Chapter 5. On practical three-dimensional TEM inversion synthetic data functions correctly however the result does contain some undesirable features. First the recovered conductivity model shows an object which is too small, and too close to the surface as compared to the true model. Second there are a number of high amplitude artifacts apparent at the surface and shown in Figure 5.4b which appear to be strongly related to the location of the receivers and are located directly above the target block in the key area of interest for the survey. These issues and others will be dealt with in detail in the subsequent sections by the incorporation of a sensitivity based weighting scheme on the model objective function. Additionally the recovered results indicate that the inversion, even at this very low noise level, appears to be insensitive to the existence of the portions of the conductive block below roughly 400 m depth. This indicates that inversion of the San Nicolas U T E M survey will likely be only be capable of defining the top of the deposit since the actual deposit extends down to a depth of more than 550 m. Before addressing the issues with this inversion of synthetic data in detail we will perform a preliminary inversion of the San Nicolas field data set in order to see what a typical inversion using generic inversion parameters produces given the San Nicolas U T E M data set. 5.3 Preliminary inversion of field data Preliminary inversions should take place on a relatively coarse mesh so that they can be com-pleted in a short amount of time. Often preliminary inversions highlight some inadequacy with the data, discretization or inversion methodology and in many cases the first inversions of a new field data set may fail to produce results because of an error made at a preparatory stage in the inversion process. For these reasons the inversion should be coarse with an emphasis on speed and merely provide enough detail to determine whether the inversion process is producing results that are, to first order, reasonable. For the San Nicolas dataset the preliminary inversion involved a spatial discretization of 138,240 cells in a reduced volume approach. The U T E M waveform was discretized into 38 time steps over 5/4 of a cycle with the data sampled on the final ramp. 9 of 10 possible U T E M channels were used. Standard deviations were assigned based on a percentage ranging from 7% for channel 9 data to 3 % for channel 1. The floor value of the standard deviation was set to 5 % of the median of the absolute value of the data from each channel. In the absence of detailed 62 Chapter 5. On practical three-dimensional TEM inversion a) 128. South Figure 5.5: a) East-west profile through preliminary inversion result for transmitter loop 2 at 450 m south. The view is from the south b) Map of preliminary inversion result at 0 m depth showing location of transmitter loop and receivers information about a system noise floor for the U T E M survey this formulation has been found to be robust. The data were fit to value of x 2 = 1- The inversion took 30 iterations over 8 values of the regularization parameter 0. The results of the first inversion of the actual San Nicolas data are shown in Figure 5.5 while the predicted and observed data are shown in Figure 5.6. The inversion result suffers from a number of inadequacies which will be discussed in detail in the next section. 5.4 Evaluation of preliminary results Upon achieving some sort of inversion result, the result itself must be evaluated. There are a number of key questions that can provide insight into the results at this point • Did the inversion converge? 63 Chapter 5. On practical three-dimensional TEM inversion -1753 -1402 -1051 easting -700 -6 66e-008 -8 24e-008 -9 82e-008 -1 14e-007 -1 30e-007 -145e-007 1.-1 61e-007 T/s -1753 -1402 easting d) 7 796-008 3 85e-008 14e-007h7s -8 83e-010 •4 03e-008. -7 9Se-008 -5 -4 log?) -3 observed predicted i°g(t) Figure 5.6: a) Observed U T E M channel 2 data map b) Decay curve at the location [-850,-310,1] c) Predicted U T E M channel 2 data map d) Decay curve at the location [-1700,-350,1] Does the data fit look reasonable? Does the model look reasonable? Not only is it important that these questions are asked, it is important that they are asked in the order presented here. For example, the usefulness of examining an inversion result in which the algorithm has failed to converge is questionable, as is the usefulness of an inversion result where the data fit is poor. Generally if the inversion fails to converge, or if the data fit is poor, it is necessary to revisit the preparatory stages of the inversion in order to determine the source of the problem. Poor convergence and data fit are most commonly caused by errors in data formatting, improper assignment of standard deviations, or excessive erroneous and outlying data points. Less commonly the discretization, either in space or in time, may be the source of the problem. If the algorithm converges and the data fit is good then it is possible to proceed with an evaluation of the inversion model itself. The resultant inversion model should be something that is geologically reasonable; exactly what this means will vary from case to case. This is a difficult statement to quantify. If the model appears to be "reasonable" then it is possible to 64 Chapter 5. On practical three-dimensional TEM inversion proceed to an inversion on a much finer spatial discretization. If the model is "not reasonable" then further investigation is necessary including perhaps inversion of additional synthetic data. A specific example of how to evaluate a preliminary inversion result, including making reference to a priori information will be discussed for the San Nicolas deposit. This should aid in the clarification of this somewhat subjective issue. For the San Nicolas preliminary inversion the algorithm converged as expected and the data fit was, in general, quite good. Figure 5.6 shows the predicted and observed data and it is apparent that the data are fit to higher degree at late time than early time. Recall that an analysis of modelling errors showed a distinct time dependance, with a tendency for error to decrease with increasing time, so this type of fit is exactly that which was intended. In general the data fit is excellent with the exception of a few data points which are considered to be outliers. This leads to an evaluation of the inversion model. A profile through the inversion model is presented in Figure 5.5a. Inspection reveals that the inversion model does not show the type of features that a priori information has indicated should be present at San Nicolas. Most notably a thick and largely continuous conductive overburden is absent from the result. In its place a number of conductive bodies have been generated at locations to the north-west and to the east of the deposit. There is a conductive body situated in a similar lateral location as the San Nicolas deposit in this inversion result, however its depth is far too shallow. Additionally in Figure 5.5b, a plot of the surface conductivity of the recovered model shows that a number of extremely shallow conductive artifacts are present which appear to be closely related to the position of the receivers. This is a situation which has been seen previously in the results for a synthetic inversion of the conductive block. In a number of ways it is a very similar result to that of the synthetic inversion shown previously. This analysis leads to the implication that there may be a fundamental problem with in-verting data using this particular configuration of source and receivers. It turns out that the problem is one of non-uniqueness. The inversion is able to find a minimum norm solution by, primarily, changing the conductivity of the cells near the transmitter and receivers. It is the high sensitivity of the data to the conductivity of these model cells which is causing the problem. What is required is a method by which the effect of the variation in sensitivity between cells near the transmitters and receivers and those farther away can be compensated for. Ideally we 65 Chapter 5. On practical three-dimensional TEM inversion would like the situation to be such that each cell is equally likely to change from the value of the reference model. EH3DTDinv has a mechanism for weighting cells in the model objective function by the use of the weighting matrix W t - The method that we will use here is commonly referred to as sensitivity weighting and has been previously used in several other types of geophysical inversion. An example of its use is provided in Magnetometric Resistivity (MMR) inversion by Chen et al. (2002). 5.4.1 Sensitivity weighting Due to the nature of the San Nicolas U T E M survey, specifically the fact that the E M loop transmitters are very close to the area where data were collected and only a few source ge-ometries were used, we find that the magnitude of sensitivity in the model cells varies in an extreme manner over the volume of investigation. Small changes in the conductivity of the cells closest to the transmitters, or more importantly, close to the data recording locations, tends to produce large changes in the predicted data. We find that without the information provided by many source configurations, the inversion tends to focus on changing the conductivity of these few highly sensitive cells, and to a large degree neglects changing the conductivity of other cells in the volume of interest. This produces undesirable artifacts even in the highly controlled and low noise situation of synthetic data, as seen previously in section 5.2.1. These artifacts are also apparent in the inversion of the San Nicolas U T E M data as illustrated by the similar conductivity anomalies shown in the synthetic inversion in Figure 5.5. In the field data situation we note that the problem tends to be complicated by the fact that the data contain unknown distributions of measurement error. The consequence is that it is necessary to develop a method of attenuating these artifacts in order to improve the inversion results for both synthetic and real data . This situation has been noted before in the literature of the geophysical inverse problem by others such as Chen et al. (2002) and it is a common situation in downhole electrical resistivity tomography (ERT) and M M R surveys. In these cases a weighting matrix, whose individual elements are scaled integrations of the magnitude of the sensitivity of each cell over the entire set of data, is often used to balance these undesirable effects. The program E H 3 D T D does not explicitly calculate the sensitivity matrix J . We follow the 66 Chapter 5. On practical three-dimensional TEM inversion approximate sensitivity calculations as outlined by Chen et al. (2002). This method is based on the current density in a representative background halfspace. The calculations are somewhat simplified in the case of the San Nicolas survey, as the survey used only measurements. In order to proceed with the sensitivity weight calculation, we require the current density J C for each cell in the model at each time when measurements were made. Given that the U T E M survey at San Nicolas is approximately the same as measuring the step off response, recall that we saw in a previous section that the difference between the two is roughly 4%, then the approximate sensitivity for can be calculated by using the sensitivities from H Z due to a step off. This saves considerable time in the forward calculation. The current densities are calculated using E H 3 D T D . The basis of the method is that the existence of a current in each cell can be thought of as a source of a magnetic dipole field. Under these assumptions the individual elements of the sensitivity matrix, [ ^ f ] can approximated analytically. = ^[(y-y')Jcx-(x-x>)jy} (5.3) where AV is the volume of the cell in question, is the current density in the x direction and J | is the current density in the y direction in the cell. The primed coordinates are the location of the cell's center while the unprimed coordinates define the location of the datum in question. See Chen et al. (2002) for details on the derivation of this equation. Given the approximate sensitivities calculated in this manner the approximate sensitivity matrix, J , can be assembled. We define the sensitivity weight as follows ^ = 7 ( K ) 1 / 4 (5-4) Here is an element of the approximate sensitivity matrix and we define 7 as a constant to be chosen by the user. TV is the number of data. For application to the San Nicolas survey we chose 7 such that the largest element in w is roughly 10. This results in a large range in sensitivity weight values and in order to ensure that the weights do not unduly affect cells which are largely undetermined by the data, for example in padding region of the mesh, we threshold w at a value of 1. This is the procedure used to weight the models' objective functions in the following inversions of both the synthetic and San Nicolas U T E M data. dHz dm 67 Chapter 5. On practical three-dimensional TEM inversion 5.4.2 Inversion of synthetic data using sensitivity weighting After making a significant change to the inversion methodology such as using a sensitivity weighting matrix the effects should be evaluated on a synthetic example before using it on real field data. We use the forward modelled data for the conductive block corrupted with 2.0 % Gaussian random noise as before. The data were inverted using an assigned standard deviation of 2.0 % along with a small floor value. The reference model and initial model used was a 200 Clm halfspace. The inversion took place on a reduced mesh consisting of 142,002 cells and with waveform consisting of 38 time steps. The objective function minimized was exactly as before with the exception that the user defined weighting matrix, W t contained the sensitivity weights as described in the previous section. Convergence required 22 total iterations over 6 values of (3. The cooling of the objective function proceeded from an initial value of 0 of 10 and the value of P was reduced by a factor of 5 until reaching a final value of 3.2 x 1 0 - 5 , at which point the target misfit criteria, x 2 = 1, was met. The inversion converged to a final model which had a value of the model objective of 11,181. The true value of the model objective for the block model used to generate the data was slightly over 386,000. The results of the inversion are shown in figure 5.7 and it is clear that the sensitivity weighted result is vastly improved from the previous result. In the sensitivity weighted inversion the block is recovered at the correct location, at the correct depth and with a conductivity that is close to the true value. Considering the fact that the data were corrupted with noise and considering the inherent blurring effect of the model objective function this result is quite good. It is also apparent from the two map views of the surface of the recovered models in Figures 5.7b and 5.7d that the near surface artifacts associated with the position of the receivers have been highly attenuated by the sensitivity weighting. To conclude, the sensitivity weighted inversion converged faster, to a better result, with more structure and produced a model whose structure resembled the true model much more closely than the previous inversion of this data set. Sensitivity weighting provides a first order improvement in the recovered result for the block in a halfspace model given the survey geometry used at San Nicolas. The result still shows an unexpected asymmetry with the block being resolved slightly better to the west than the east. This is primarily due to the existence of much higher current densities on the west side of block than the east side as the block is positioned off-center in the transmitter loop. This 68 Chapter 5. On practical three-dimensional TEM inversion indicates that better results may be achieved by simultaneously inverting data from multiple transmitter loops. This should improve resolution of the recovered model. Not only will this help to resolve some of the non-uniqueness issues with the data but this should produce a more even distribution of current density, and as a consequence, sensitivity in the inverse problem. 5.4.3 Inversion of synthetic data using multiple transmitters In this section the data, forward modelled from 3 large loop transmitters over the synthetic model of the conductive block shown in Figure 5.4, will be simultaneously inverted and the results will be presented. We predict that results will illustrate the effectiveness of the multiple transmitter inversion and the improvement in results that can be achieved by using data from several transmitters to reduce the non-uniqueness of the inversion problem. We also predict that the notable asymmetry of the inverted result shown in figure 5.7 will be attenuated by the simultaneous inversion of multiple transmitters due to a more even distribution of current density throughout the model. The transmitter and data locations are those of the San Nicolas U T E M survey. The data from all 3 loops had 2.0% gaussian random noise added. The combined forward modelled data for 3 loop transmitters were inverted using an assigned standard deviation of 2.0% along with a small floor value. Both the reference model and initial model were a 200 flm halfspace. The inversion took place on a reduced volume consisting of 142,002 cells and with waveform consisting of 38 time steps. The objective function to be minimized was that of equation 5.1. Convergence required 26 total iterations over 8 values of 0. The cooling of the objective function proceeded using an initial value of 0 of 10 and the value of 0 was reduced by a factor of 5 until reaching a final value of 6.4 x 10~4, at which point the target misfit criteria, x 2 = 1, was met. The inversion converged to a final model which had a value of the model objective of 18, 669. The true value of the model objective for the actual block model used to generate the data was slightly over 386,000 confirming that, as is necessary for this algorithm, the inversion is producing a model with less structure than the true model. Figure 5.8d shows the result of an inversion of the now familiar block model and provides a direct comparison with previous inversions of the same model. The figure clearly illustrates the ascending quality of synthetic inversion results that have been achieved through iteratively inverting data and evaluating the results. In the following sections we wish to apply the tech-69 Chapter 5. On practical three-dimensional TEM inversion Figure 5.7: a) East-west cross section though conductive block model at 450 m south. The view is from the south, b) Map view of surface conductivity for inversion of conductive block without sensitivity weighting at depth of Om. Showing the locations of data (black points) and transmitter loop (black lines), c) East-west cross section though recovered model without sensitivity weighting at 450 m south. The view is from the south, d) Map view of surface conductivity for inversion of conductive block inversion with sensitivity weighting at depth of Om. Showing the locations of data (black points) and transmitter loop (black lines), e) East-west cross section though recovered model with sensitivity weighting at 450 m south. The view is from the south. 70 Chapter 5. On practical three-dimensional TEM inversion Figure 5.8: a) East-west profile through block model used to generate synthetic data. The view is from the south, b) East-west profile through inverted result using and data from transmitter loop 2 and employing a sensitivity weighting. The view is from the south, c) East-west profile through inverted result using data from transmitter loop 2 without sensitivity weighting. The view is from the south, d) East-west profile through inverted result using data from 3 transmitter loops and employing a sensitivity weighting. The view is from the south. niques developed here, using synthetic data, to enhance the results for the inversion of the San Nicolas U T E M data set, 5.5 Further improvements to inversion As the process of inversion continues improvements can be made to refine the data set itself and improve the parameters in some of the preparatory stages. These improvements include editing the data and refinements to the discretizations. 5.5.1 Data editing This is an important stage in the inversion process. Because of the statistical framework of the data misfit function, that is Gaussian statistics, the existence of severe data outliers can cause the inversion to fail, proceed slowly, or result in models with undesirable artifacts. For this reason all data outliers should be removed from the dataset. The key question here is to determine how to identify a data outlier given a field data set. It is very difficult to identify 7 1 Chapter 5. On practical three-dimensional TEM inversion an outlying data point with complete confidence, however it is possible to identify a set of data points which may be outliers. The basic concept is that since the inverse problem already under-determined and the data recorded tend to be linearly dependent on each other we do not suffer in a significant way by mistakenly eliminating a few good data points, at least when compared the damage which is done from keeping a few bad data points. The procedure here is to first visually inspect the complete data set. A visual inspection will help identify any severe outliers, that is, any point which defies the smoothness properties that diffusive E M fields are known to have. Note that while the magnetic field is necessarily spatially smooth the electric field is not. Care must be taken when editing electric field data for this reason. The second stage is to compare the predicted and observed data from the preliminary inversion. In many cases the majority of the total misfit from a preliminary inversion will be attributable to a very small number of data points. These data should also be treated as outliers and be removed from the dataset. This process should be performed iteratively and outlier removal may require several successive inversion results to complete. For the San Nicolas U T E M survey the outlier removal was done over the course of two successive preliminary inversions with the result that each individual transmitter's data set contained approximately 5-15% fewer data points than originally provided. This reduction included data removed due to the data's proximity to the source wires. 5.5.2 Refinement of discretization Discretizations, especially in the spatial domain, should be improved throughout the course of an inversion. As mentioned previously, the preliminary discretization should be optimized for quick results. Once the inversion can be shown to be working correctly then the spatial discretization can be redefined to allow an increase in the accuracy of the forward modelling and achieve higher resolution of geologic formations. The preliminary inversion at San Nicolas took place, as previously stated, on a mesh consisting of 138,240 cells. Several other finer discretizations were developed for the San Nicolas U T E M survey. The results presented for the following inversions in this section will employ a very highly refined mesh of 241,920 cells. This mesh is optimized to have cell aspect ratios as close to 1 as possible. This is important because by keeping aspect ratios close to 1 we improve the speed at which conjugate gradient solutions 72 Chapter 5. On practical three-dimensional TEM inversion of the forward modelling systems proceed. 5.5.3 Mult iple transmitter inversion Because computation times increase rapidly with the number of sources the algorithm EH3DTDinv is best used with only a few sources and where there are many receivers for each source. The use of multiple sources, however, can lead to an improved inversion result. Adding multiple sources tends to reduce the non-uniqueness of the resultant solutions as we have shown in the synthetic data examples. This has the effect of further improving inversion results. The San Nicolas U T E M survey has 3 large loop sources which have many hundreds, and in some cases thou-sands, of data recorded for a single transmitter configuration. Subsequent inversions presented here will use some or all of these large loops in conjunction to provide increased resolution of geologic structures. 5.6 Inversion of field data incorporating improvements Once the inversion is perceived to be proceeding smoothly it is possible to include all of the previously mentioned improvements in the inversion procedure. Note however, that this process of parameter improvement is iterative and one must proceed cautiously in order to completely understand the impact of each change in the inversion process. At a minimum to get to this stage, one must have completed at least 5, and possibly many more inversions, depending on the number of transmitters available in the survey area. For the San Nicolas data set three transmitters were available, multiple background models were used and several iterations of data editing were done. This intermediate stage in the inversion process was reached only after more than fifty inversions of the individual transmitters. This is many more than would be strictly necessary. However since the intention is that the San Nicolas U T E M data set will be used as archetypal case history for the inversion software EH3DTDinv, exhaustive testing on the data set was undertaken. In this section we will present the results of an inversion incorporating data from 2 transmitters, loop 1 and loop 2. The inversion employed a sensitivity weighting based on the arithmetic average of the sensitivity weight calculated for each of the two loops. The spatial discretization consisted of 241,920 cells with aspect ratios of 1 and the time discretization was, the now standard, 38 time steps over 5/4 of a cycle of the U T E M waveform. The 2 individual 73 Chapter 5. On practical three-dimensional TEM inversion hi! • b) Q m Figure 5.9: a) East-West profile through the recovered model at 450 m south from the im-proved or intermediate inversion. An outline of the geologic boundaries as interpreted from drilling is overlain on the result b) 20 Qm iso-surface through inverted result. Locations of both transmitter loops are shown, c) East-West profile at 450 m south through a 3D conductivity model generated from drill core conductivity testing. An outline of the geologic boundaries as interpreted from drilling is overlain on the conductivity model data sets from each transmitter underwent a two-stage editing process for outlying data points as discussed in the previous section. Both the initial and reference model for this inversion were 100 Qm halfspaces. 9 of 10 possible U T E M channels were used. Standard deviations were assigned based on a percentage ranging from 7% for channel 9 data to 3% for channel 1. The floor value of the standard deviation was set to 5 % of the median of the absolute value for the data in each channel. The inversion proceeded over 34 total iterations and 9 values of the regularization parameter 0. The data were fit to value of x2 = 0.5 as a higher degree of fit was desired for this inversion than in previous inversions. The result is shown in Figure 5.9. The predicted and observed data are shown in Figure 5.10. The data fit is, in general very good, and the model is consistent with both data sets. 5.7 Evaluat ion of improved results Once these types of improvements have been implemented to the satisfaction of the inversion operator it is necessary to evaluate the results. The question to ask is whether the changes to the 71 Chapter 5. On practical three-dimensional TEA1 inversion -1350 5.95e-008 I 4.90e-008 . 3.85e-008 . 2.79e-008 . 1.74e-008 . 6.§5e-009 . -3.64e-009 T/s 2216 -1927 -1639 Easting [m] -1350 5.95e-008 4.90e-008 3.85e-008 2.79e-008 1.74e-008 6.^ 8e-009 -3.64e-009 T/s Figure 5.10: a) Observed U T E M channel 3 data map for transmitter loop 2 b) Observed U T E M channel 3 data map for transmitter loop 1 c) Predicted U T E M channel 3 data map for transmitter loop 2 d) Predicted U T E M channel 3 data map for transmitter loop 1 75 Chapter 5. On practical three-dimensional TEM inversion inversion parameters improve the quality of the result. The answer hinges on how one evaluates the result of an inversion. Any sort of evaluation of results requires some kind of knowledge about the geology in the area. Additionally some knowledge is necessary to determine which features of the model are artifacts of the inversion process and which features of the model are demanded by the data. The only way to really get an understanding of these issues is to have access to geologic knowledge and physical property estimates in the area and to have completed multiple inversions with many differing parameters. Additionally the experience of the user in the field of inversion tends to be invaluable when interpreting and evaluating results. For the San Nicolas U T E M survey we are at a significant advantage when it comes to evaluating the results of inversion. The geology is well understood as are the variations in physical properties in the area. In fact a 3D conductivity model of the immediate area around the deposit has been assembled by Phillips (2001) and has been previously presented in chapter 4. This 3D model will be used to help evaluate the inversion results achieved in this and subsequent sections. The result of the improved, 2 transmitter, inversion is shown in Figure 5.9. Here an outline showing geologic boundaries is overlain on a profile through inversion result at 450 m south. Also shown is the same profile through a 3D conductivity model which was generated from conductivity measurements made on the drill core. From this figure it is apparent that this result is vastly improved over the preliminary results. A conductive overburden layer is present. Additionally there is a conductive anomaly in a similar lateral location as the deposit. This is encouraging. However there is no clear separation between the overburden and the conductive anomaly and the depth to the large conductive anomaly is underestimated in the inverted results. It is apparent that, though a conductive overburden layer is present, the layer's structure does precisely not match the situation interpreted from drilling. It is interesting to note that the surface cells in this inversion have very similar values to those of the reference model while those directly below differ considerably. Considering the magnitude of the sensitivity weighting at the surface this should not be surprising as there is an extremely large penalty to paid in the objective function for any deviation from the reference model for these surface cells. It is well known that the surface resistivity in the San Nicolas area, ignoring the thin veneer of the quaternary alluvium which in places is less than l m thick, is actually much less than 100Vim. It is likely that by enforcing the conductivity of the near surface cells to be that of the 100 fim 76 Chapter 5. On practical three-dimensional TEM inversion reference model we are biasing the inversion. This is the penalty to be paid for using the sensitivity weighting. In the areas of high weighting the model is largely pre-determined by the reference model. This is why it is essential that we use the correct conductivity for the surface cells. This is a subject that will be discussed at length in the next section on the subject of including a priori information in 3D inversions. 5.8 Inversion of field data using a priori information Once the inversion is working correctly and reasonable results are being produced it is advan-tageous to include some of the a priori information that is available and which may improve the results of the inversion. It is often the case that a small amount of additional a priori information can produce a substantial improvement in inversion results. The type of infor-mation used may come from a study of the surface geology, boreholes or even some inferred knowledge about the dip or probable shape of structures. Sometimes incorporating this type of information into the inversion can be difficult or time consuming and as a result practitioners of inversion often do not make the effort to proceed to this stage. This is a key problem for the field of geophysical inversion and it degrades the usefulness of inversion in general. Recall that inversion is an under-determined problem and whether or not the user is actively aware of it, the algorithm is using some estimate of the true conductivity to regularize the solution. Often substantial improvements can be made by spending time to create an estimate which is as close as possible to the true conductivity. The San Nicolas U T E M survey provides a prime example of a situation where a small amount of a priori information can provide an important improvement in the inversion result. Recall from the previous section that the sensitivity weighting, in combination with the inaccurate reference model for the near the surface, was tending to create an undesirable resistive layer at the surface. The solution to this problem is to incorporate a priori information about the overburden into the inversion. The problem is that, although we know an overburden exists, we do not know exactly its thickness nor how the thickness varies laterally. Recall that in the previous chapter some of the forward modelling analysis done on the San Nicolas U T E M data indicated that the overburden thickness and possibly its conductivity varied significantly over the survey area. We do not want to bias the inversion result by the inclusion of an incorrect thickness 77 Chapter 5. On practical three-dimensional TEM inversion of overburden in the reference model. The method we use to solve this problem involves a relaxation of the smoothness part of the objective function. Recall that the regularization function R(m — mref) used in the previous inversions shown here was given by R ( m - m r e f ) = a s | | W s W t ( m - m r e f ) | | 2 + a I | | W t W x ( m - m r e f ) | | 2 + (5.5) ay | | W t W y ( m - mref)!!2 + « z | |W t W z (m - m r e f)||f This means if a horizontal conductive layer is included in the reference model, mref, and the thickness is not correct then the smoothness term, az | |WtW z (m — m r ef)|| 2, will tend to require an interface at the location specified in the reference model regardless of whether one exists or not. The solution is to remove the reference model from the derivative terms of the regularization function. The new regularization function will be redefined to be R ( m - m r e f ) = a s | | W s W t ( m - m r e f ) | | 2 + a x | | W t W x m | | 2 + (5.6) ay | | W t W y m | | 2 + az | |W t W z m| | f This allows the algorithm to preferentially choose models with overburden layers but be flexible about locating the interface between overburden and host. The reference model used for this case was a 100 Qm basement with a 90 m thick layer of 20 Qm overburden. These values reflect the interpreted average overburden thickness and conductivity given the a priori information available at San Nicolas. It is not out of the question to imagine that this type of information would be available even in the most "greenfields" sort of exploration areas. Even if the detailed drilling information did not exist at the San Nicolas survey site the overburden conductivity could have been acquired inexpensively through a number of independent means, from laboratory measurements of surface samples to the ,use of the reconnaissance airborne E M surveys that have been conducted in the area. The average thickness could have been estimated from an analysis of the local surface geology or from well logs available from previous water well, mineral, oil exploration drilling in the area, or analysis of airborne E M data. In the interest of correctly recovering the overburden it is possible to use part of the reg-ularization function to give a selective advantage to models which contain horizontal layered structures. This can be done through the use of the a parameters by reducing the value of the az parameter such that the ratio of ^ and ^ is larger than 1. This will decrease the penalty 78 Chapter 5. On practical three-dimensional TEM inversion on vertical structures relative to horizontal structures allowing horizonal layered models to be selectively chosen from the suite of those that fit the data . In this section we will present the results of an inversion incorporating ^ data from all 3 large loop transmitters. The inversion employed a sensitivity weighting based on the arithmetic average of the sensitivity weight calculated for each loop. The spatial discretization used 241,920 cells with aspect ratios close to a value of 1 and the time discretization was, the now standard, 38 time steps over 5/4 of the U T E M waveform. The 3 individual data sets from each transmitter underwent a two-stage editing process for outlying data points as discussed previously. 9 of 10 possible U T E M channels were used. Standard deviations were assigned based on a percentage ranging from 7% for channel 9 data to 3% for channel 1. The floor value of the standard deviation was set to 5 % of the median of the absolute value for the data in each channel. The value of — and ^ was set to 5. The reference model for this inversion was a 2 Cfz OLz layer earth with a 90 m thick 20f2m layer at the surface overtop of a 100 flm basement. The regularization function was modified as discussed previously in this section. In order to save time, the inversion was initiated with a result from the previous inversion shown in Figure 5.9. The starting value of 0 was modified to be 100 times higher than the final 0 from the previous inversion. This process of restarting an inversion with an initial model from a previous result will be discussed at length in the next section. The inversion proceeded over 13 total iterations and 4 values of the regularization parameter 0. The data were fit to value of x 2 = 0.5. The result is shown in Figure 5.11. The predicted and observed data are shown in Figure 5.12. The data fit is, in general, very good and the model produced is consistent with all data sets. 5.9 Evaluation of results with a priori information Adding a priori information is an iterative process. Simple information should be used first and more complex information added in later inversions. By simple information what is meant is information about the generalities of the survey area. Complex information includes detailed information acquired as as result of an independent geophysical survey or some kind of drilling program. For example simple information may be something as basic as recognizing the exis-tence of an overburden layer. Complex information may consist of a detailed description of an overburden conductivity and thickness. By proceeding cautiously it is possible to differentiate 79 Chapter 5. On practical three-dimensional TEM inversion Figure 5.11: a) East-west profile at 450 m south through the recovered model from an inversion with a priori information added. An outline of the geology, as interpreted from drilling, is overlayed on the result. The view is from the south, b) 20 Qm iso-surface through inverted result. Locations of all 3 transmitter loops are shown c) East-west profile at 450 m south through a 3D conductivity model generated from drill core conductivity testing. An outline of the geology, as interpreted from drilling, is overlayed on the conductivity model. The view is from the south. SO Chapter 5. On practical three-dimensional TEM inversion Figure 5.12: a) Observed U T E M channel 4 data map for loop 2 b) Observed U T E M channel 4 data map for loop 1 c) Observed U T E M channel 4 data map for loop 9 d) Predicted U T E M channel 4 data map for loop 2 e) Predicted U T E M channel 4 data map for loop 1 f) predicted U T E M channel 4 data for loop 9. Units here are volts s i Chapter 5. On practical three-dimensional TEM inversion aspects of the model which are determined by the reference model and aspects of the model which are determined by the data itself. This helps to avoid the major hazard inherent in adding a priori information, letting preconceived notions of the geology, rather than the actual data, drive the interpretation. The inversion presented in Figure 5.11 for the San Nicolas U T E M survey represents our best estimate for the conductivity of the San Nicolas area and is perhaps the most accurate inversion result for any E M survey conducted at the site. The result clearly defines a conductive overburden whose structure agrees very closely with the drilling results. The inversion also locates the upper portion of the San Nicolas deposit accurately, both laterally and at depth. It is clear that, given this recovered conductivity model, a very efficient drilling program could be implemented to define the deposit and quickly determine its value. The inversion result does not define the bottom portion of the deposit referred to as the deposit's "keel". This is not a failure of the inversion methodology. This is simply a result of the fact that the bottom of the deposit is below the depth of investigation of the U T E M survey given the base frequency of the U T E M system used by the contractor. This is essentially a survey design problem and if information about the bottom of the deposit is required then modifications to survey type or design would be able provide more information. Certainly by conducting the U T E M survey at a base frequency of 2 Hz, a better sensitivity to the deeper structure would be achieved. At 2 Hz the diffusion depth of U T E M Channel 1 is roughly a factor of four deeper than that of the 30.974 Hz U T E M survey. However if the purpose of the survey was to define the deeper portions of the deposit the question of whether U T E M would be as effective as other deep-seeing E M surveys such as Audio Magnetotellurics (AMT) is worth considering. A n A M T survey would certainly be sensitive to conductors at, and below, 600 m depth. The effectiveness of A M T surveys is highly dependant on the strength of the natural source fields used in the method and the strengths of these fields are dependant on the geographic location of the survey and the amount of solar activity at the time of acquisition. Under certain conditions a U T E M survey conducted at a low base frequency may be more effective than A M T at locating the deep portions of the deposit. In either case this would be a difficult task and very high accuracy data would need to be acquired as the large, conductive upper regions of the deposit tend to mask the response to smaller lower portions. Borehole E M surveys may be necessary to properly detect and define the deposit "keel". 82 Chapter 5. On practical three-dimensional TEM inversion 5.10 Conclusion To summarize the achievements in this chapter, we have generated synthetic data from a simple conductivity model and inverted them to produce a result which has minimal model structure and fits the data to the required degree. The results of this first inversion of synthetic data suffered from some inadequacies and we have shown that the preliminary inversion of the San Nicolas U T E M data suffered from very similar sorts of problems. This led to the inclusion of a sensitivity weighting in subsequent inversions which improved results. Despite the improve-ments gained through the sensitivity based weighting scheme, the synthetic inversions showed some residual deficiencies in resolving the conductive block. Further inversions were carried out on the synthetic data using data from multiple transmitter configurations. The inclusion of data from multiple transmitters led to good results for the inversion of the synthetic dataset. This methodology was applied to the U T E M data set where a sensitivity weighted, multiple transmitter inversion was conducted on the San Nicolas U T E M data set with much improved results when compared to the preliminary inversions. The next stage was to incorporate some simple a priori information about overburden conductivity and thickness into the inversion. Inversions using this a priori information produced the results shown in the previous section which are very similar to the conductivity model generated from drill core measurements and contain boundaries which corresponded well with the interpreted geology. The final result shown in the previous section is an example of a highly successful application of a 3D inversion algorithm for T E M data to a field data set. 3D inversion of time domain E M data has previously been completed by Newman & Commer (2005) and by Commer et al. (2003) however in the first instance only synthetic data are inverted while in the second a very limited type of model is inverted for. The result presented here is thought to be the first example of a large scale inversion of actual time domain E M data that have been recorded in a field survey. Having achieved an inversion result of this quality it is possible to bring the discussion back to some practical matters. In the next section we will discuss the practical aspects of 3D T E M inversion with emphasis on using the various parameters of the inversion algorithm to achieve results quickly and efficiently. 83 C h a p t e r 6 On practical and efficient T E M inversion 6.1 Introduction In this chapter the discission will be centered around the process of making inversion of 3D T E M data both practical and efficient. As these topics are intimately related they will be discussed together here. The topic of practical inversion involves many considerations including survey design, analysis of measurement errors, data formatting and 3D model building. Efficient inversion means ensuring that inversion results are achieved in a timely manner. This also includes consideration of survey design along with the choice of initial models, the cooling schedule, and level of discretization. Once these issues of efficiency have been discussed in a general sense the procedures developed here will be employed on the example of the San Nicolas U T E M data set to illustrate the effectiveness of the strategies and demonstrate the substantial time savings that can be achieved by using these techniques. 6.2 Practical inversion of 3D T E M data The process of making inversion of 3D T E M data truly practical requires a change in the way geophysical data are collected, processed, stored and interpreted. This will require efforts by geophysical contractors, clients and possibly the creation of a new field of specialized geophysical processors and "inversionists". The process of changing the way T E M data are used will necessarily take some time and effort to implement. However, the benefits that can be achieved should make the effort worthwhile. 8 4 Chapter 6. On practical and efficient TEM inversion 6.2.1 Survey design There is an improvement that can be made in the interpretation of geophysical data by designing geophysical surveys with the concept of inversion in mind. Surveys should be designed to have receivers positioned such that the fields can be accurately modeled. Recall that in order for numerical modelling to achieve a sufficiently accurate result receivers must be located far enough from the transmitter for numerical errors associated with the discretization of the source to be minimal. The results from San Nicolas indicate that inversions employing data from several transmit-ters produce substantially better results, however, each additional transmitter adds considerable computational cost to the inversion process. In contrast the addition of more data locations for each transmitter adds very little computational cost. This means that the ideal survey should contain a few transmitters, each with very different geometries, in order to illuminate differing portions of the target area with currents. Each transmitter in this survey should provide fields at many data locations and if possible at many times. From a logistics point of view this is desirable since setting up the transmitter often represents a large portion of the time spent acquiring data. This configuration is, in effect, a description of the survey design for the 3 large loops used in the preceding inversion of the San Nicolas U T E M survey and it explains, at least in part, the reason for the success that was achieved with the data set. 6.2.2 Analysis of measurement errors If inversion is to become truly practical then serious attempts must be made to analyze system measurement errors. This includes simple measurement and analysis of background noise levels due to natural and anthropogenic sources. It also includes analysis of errors due to mis-location and mis-orientation of transmitters and receivers. These are topics which are most often ne-glected by both geophysical contractors and clients. The provision of this kind of error analysis would greatly enhance the process of assigning standard deviations to the data in the inversion process and as a result improve the quality of the inversion result. 6.2.3 Data Formats Data formats, including the way in which data are stored, should be modified to enhance the utility of inversion. Data are often normalized into a format which attempts to be easily 85 Chapter 6. On practical and efficient TEM inversion interpretable by the human eye. While this process of data inspection and visual interpretation is invaluable to the geophysicist, all too commonly the data are converted to a normalized format and then normalization factors are difficult or time consuming to locate and remove after the fact. The data should to be recorded and output in at least two formats, in SI units and in a processed and normalized format. Furthermore the cost of data storage has dropped to levels where stacking time domain signals before recording makes little sense. We suggest all time domain systems should now record full waveform data. This would allow many beneficial processing steps to be done to the data before stacking and allow a more sophisticated analysis of the measurement error due to background noise levels. As 3D interpretation of time domain data becomes more popular, orientations of data need to become clearly and uniformly labeled. While the directions x and y are sufficient if 2D data are to be acquired and interpreted, when data are acquired with the intention of 3D processing then consistent coordinate systems and orientations need to be used for all data regardless of the line on which they were recorded. 6.2.4 Construction of 3D models Data acquisition is not the only area related to T E M inversion where improvements can be made. In order that the geophysical inversion can proceed as outlined here it is required that the geophysicist have access to utilities whereby three dimensional models can be created quickly and efficiently. There are a number of reasons for this, including the purposes of hypothesis testing through forward modelling, inclusion of a priori information, and the creation of weighting functions. Without 3D model building utilities such activities can become time consuming to the point where they are frequently neglected and the substantial benefits to interpretation that they produce are lost to the user. 6.3 Efficient inversion of 3D T E M data The inversion result achieved in the previous chapter is exceptional, yet, there has been little mention as to the amount of time that it takes to achieve such a result. In this section, the factors that affect the efficiency of the 3D inversion process will be identified and various techniques will be discussed to improve the speed at which 3D T E M inversion proceeds. Working towards 86 Chapter 6. On practical and efficient TEM inversion an efficient inversion is a critical step in the process of making 3D T E M inversion practical and routine. Once again the San Nicolas U T E M survey will be used as the archetypal example, in this case to illustrate the effectiveness of the various techniques for improving the efficiency of inversion of 3D T E M data. In order to provide some additional motivation for this section consider that the final inversion result for San Nicolas presented in section 5.8 would have taken 28 days to complete, had the initial model not been set to an inversion model achieved from a previous result. Even with the starting model set to a previous inversion result this inversion took nearly 16 days to converge. All inversions here were performed on Opteron 244 A M D 64 bit processors, the clock speeds were set to 1.8 Ghz. In an exploration setting, where time is often of the essence, and any particular inversion is merely one of many that must be done to fully investigate a data set, this represents an unacceptable length of computation time. In the following chapter an effort will be made to reduce this computational time while maintaining, as much as possible, the quality of the result. The major algorithmically determined factors which affect the length of time it takes to complete a 3D T E M inversion with the code EH3DTDinv are listed below. They will be discussed in general in the following sections and then presented together as they apply to the inversion of the U T E M data set at San Nicolas. Parallelization of forward modelling Choice of initial model Cooling of the regularization parameter Level of discretization 6.3.1 Parallelization of forward modelling The bulk of the computation for 3D E M inversion is the process of forward modelling. When multiple transmitters are simultaneously inverted, the forward modelling operations that must be performed in order to invert the data are independent of each other. In the original code, which was designed to run on a single processor, each forward modelling had to be completed in series. This meant that the length of time that a particular inversion took was roughly proportional to the number of transmitters in the inverted data set. 87 Chapter 6. On practical and efficient TEM inversion To solve this problem the algorithm was parallelized such that for each transmitter a single processor was used to calculate all forward models necessary for an inversion. Since forward modelling proceeds independently for each transmitter in this scheme, the time savings resulting from this implementation are nearly proportional to the number of transmitters used in the inversion. 6.3.2 Initial models The choice of an initial model can make a significant impact on the efficiency of the inversion process. A commonly used method to improve inversion efficiency is to choose an initial model which is thought to be a good guess at the final model. The difficulty is to still allow the algorithm the flexibility to change that model if, for any reason, a part of the initial model is incorrect or is inconsistent with the data. It is because of this need for flexibility that the choice of an initial model is intimately related to both the choice of reference model and the value of the regularization parameter. When nothing is known about the nature of the true model then it is often best to initiate the inversion with a simple model such as a halfspace. However as the process of inversion proceeds, or when a significant amount of a priori information exists, it is often possible to choose an initial model which is closer to the final acceptable model. This can have the effect of greatly diminishing the amount of time that it takes to complete an inversion since the algorithm has to step through much less distance in model space to reach convergence. In order to initiate an inversion with an estimate of the final model one must carefully con-sider the choice of the reference model and regularization parameter. This may require having completed previous inversions in order to give a general idea of the value of the regulariza-tion parameter which provides convergence of the algorithm for a similar situation. Given this knowledge it is possible to choose a starting value of the regularization parameter and cooling schedule that provides the required characteristics. The topic of the cooling schedule will be discussed in detail in the next section. We have found that good results are produced, when initiating an inversion with a previous inversion's result by using a value of 0 of 100 times that of the convergence value for the previous inversion. 88 Chapter 6. On practical and efficient TEM inversion 6.3.3 Cooling schedule for the regularization parameter The starting value of the regularization parameter, and the cooling schedule by which the regularization parameter is reduced, are important areas of the inversion process where gains in efficiency can be made. However, modifying parameters in this area is more complicated to implement in a practical setting than some of the previous subjects mentioned here. Recall that the inversion algorithm proceeds by solving a sequence of inverse sub-problems. Each sub-problem has a progressively smaller value of the regularization parameter until the algorithm proceeds to find a model which has predicted data which are acceptably close to the observed data. The key questions to be considered for development of a cooling schedule are What initial value of 3 should be used At what rate should the value of the regularization parameter be reduced How well to solve the inverse sub-problem at a particular value of p. Here we will attempt to address each of these issues in turn. Choosing a starting value of the regularization parameter that is too large will result in a series of initial inversion steps which do not change the model or the data fit in any significant way. The algorithm, when used on a default setting often chooses a value which is too conser-vative. If previous inversions of the data set have been done then this value of starting 3 can be reduced significantly by making a informed guess about what the final value will tend to be. The methodology that was used for the San Nicolas survey was to inspect the convergence curves of preliminary inversions of similar data sets. 3 values that correspond to regions of the convergence curves where the inversion begins to make significant reduction in the data misfit appears to be an acceptable compromise between efficiency and a conservative estimate at where to initiate 3 m subsequent inversions. Once an initial value is chosen the cooling schedule must be devised. The cooling schedule is made up of two components. They are the accuracy to which the inverse sub-problem is solved and the rate of reduction of the value of regularization. The discussion will focus, first, on the degree of accuracy with which to solve the inverse sub-problem for each value of 8. This is controlled through two parameters in the inversion algorithm, the number of inner iterations and the inversion tolerance based on the smallest step size that the inversion will take. While a preliminary inversion my use 2 iteration per 3 89 Chapter 6. On practical and efficient TEM inversion to achieve increased efficiency, an inversion at the final stages may use 4 or 5 iterations per (3 or perhaps even more. We do not recommend changing the minimum step size from its default value of 1 x 1 0 _ 3 5 / m at this time. The rate of reduction of the regularization parameter is controlled by a parameter called the (3 factor, f@ , which is defined as follows Aiew = Po\d X fp (6.1) The choice of this parameter should be optimized for the problem at hand. A preliminary inversion might use a factor of 0.2 or smaller in order to quickly reduce the regularization factor, achieve a high efficiency, and quick convergence. Later inversions may use higher values to proceed more slowly and produce more accurate solutions of the particular inverse problem that is being solved. 6.3.4 Discretization factors The' temporal and spatial disretization can also be changed in order to achieve a higher ef-ficiency. Note however that any reduction in the number of cells in the inverted model will generally have the effect of increasing the discretization error in the forward calculation. The magnitude of an acceptable discretization error is, as previously discussed, an issue which is largely problem dependant. Often meshes are designed in an overcautious manner and this can have a negative impact on the efficiency of the inversion process. The most common area where meshes contain excessive cells is in the padding area. For an inversion which employs the re-duced volume technique very few padding cells are required in the reduced mesh provided that large conductivity contrasts are not present. Note that the process of refining discretization for the purposes of efficiency can be difficult and may result in poor inverted models if too many cells are removed. It is best to proceed cautiously in this, starting from a discretization which produces good results and coarsening the discretization incrementally. The aspect ratios of cells are also important in improving speed at which the forward problem can be solved. Ideally the reduced volume should have aspect ratios of 1 everywhere. This may not be feasible for all applications and can be relaxed in situations where domains are required to have unusual shapes or are extremely large in a particular dimension. In the next section we will see how this applies to the San Nicolas inversion and the example will demonstrate the 90 Chapter 6. On practical and efficient TEM inversion reduction of computation time that can be achieved by applying all the efficiency considerations to a large inversion problem. 6.3.5 Application to the San Nicolas deposit In this section the recommendations for practical and efficient inversion will be applied to the San Nicolas U T E M data set. Unfortunately the data formats, orientations and recording techniques are fixed once the data have been recorded and cannot be changed after the fact. The survey design also cannot be changed at this point though it is encouraging to note that the San Nicolas survey design is near optimal for inversion with the code EH3DTDinv. When it comes to the subject of efficient inversion however, there is a considerable improvement that can be made. Before any of the efficiency considerations were addressed in the San Nicolas U T E M inversions achieving results for a large scale inversion such as the inversion shown in section 5.8 would have taken 28 days to complete. The following discussion will illustrate the impressive reduction in computation time using the recommendations of the previous sections. First the parallelization of the forward modelling portions of the inversion code was im-plemented and the new code was employed to invert the San Nicolas data. The reduction in the time required to reach convergence was roughly proportional to the number of transmitters included in the inversion. For the final San Nicolas inversion result of the previous chapter this meant a reduction in computation time of almost a factor of 3. Second the starting value of the regularization parameter was changed from the default value by using knowledge based on typical convergence curves obtained through preliminary inversions. Using the methodology outlined in the previous section a value was chosen based on values of 3 in previous inversions where large reductions in the misfit were achieved. Using an appropriate value for a starting 3, along with fewer iterations per 3 resulted in a reduction in the total number of iterations for the inversion to reach convergence. The convergence curves for two inversions of the San Nicolas U T E M data sets are shown in Figure 6.1. The efficient choice of 3 results in a 21 % percent reduction in number of iterations and as a result, computation time. The resulting models from these two inversion are virtually indistinguishable from each other. All of these efficiency enhancements, with the exception of any discretization changes, were applied to the final large scale inversion in Section 5.8. The spatial discretization used 241,920 91 Chapter 6. On practical and efficient TEM inversion x default cooling • fasl cooling Iteration number Figure 6.1: Comparison of the data mistfit function, (fid, as a function of iteration for an inversion of a 3 transmitter data set using an efficient cooling schedule and using the default cooling parameters. cells with aspect ratios close to a value of 1 and the time discretization used 38 time steps over 5/4 of the U T E M waveform. The 3 individual data sets from each transmitter underwent a two-stage editing process for outlying data points and 9 of 10 possible U T E M channels were used. Standard deviations were assigned based on a percentage ranging from 7% for channel 9 data to 3 % for channel 1. The floor value of the standard deviation was set to 5 % of the median of the absolute value for the data in each channel. The reference model for this inversion was a 2 layer earth with a 90 m thick 20fim layer at the surface overlying a 100 fim basement. The initial model was a result from a previous inversion of a single transmitter data set. The inversion proceeded over 8 iterations and 4 values of 0 from an initial value of 0.001 to a final value of 8 x 10 - 6 . The data were fit to a value of x = 0.5. A profile through the inverted result is shown in Figure 6.2a. This inversion reduced computation time to roughly 4.5 days and the resultant model is indistinguishable from, a default parameter inversion without parallelization which took slightly more than 28 days. This is a substantial reduction in computation and it has has not degraded the quality of the result. If further efficiency is required then some compromises must be made with respect to the quality of the result. The discretization can be coarsened to achieve very quick results. However, 92 Chapter 6. On practical and efficient TEM inversion this will necessarily reduce the accuracy of the forward modelling and as a consequence of the loss of accuracy, the quality of the result will be reduced. To pursue further efficiency gains, a coarse discretization of the U T E M waveform was developed employing 28 points, a 27% reduction from previous discretizations. A spatial discretization was developed, in order to reduce the total number of cells, which coarsened all padding cells on the reduced domain and further reduced the size of the domain. Interior cells remained at 30 m on a side but cells with larger aspect ratios were present near the boundary of the domain. The total number of cells in the coarsened inversion domain was 135,090, a 46 % reduction from the previous domain. The inversion employing these discretizations took only slightly more than 52 hours to complete. This means that by considering all the recommended matters of efficiency in conjunction it is possible reduce the total computation time by a factor of roughly 13. A profile through the inverted result in shown in Figure 6.2b. The quality of the result is somewhat poorer than that shown in section 5.8. The inversion shows a conductor slightly offset from the true location of the San Nicolas deposit. The over-burden layer is present and corresponds well with the 3D conductivity from drilling however the resolution is somewhat poorer than previous results. The high accuracy in the location of the deposit, the value of its recovered conductivity and much of the distinctive shape of the deposit were lost as a result of the coarse discretization. Also the continuity and resolution of the overburden layer was poorer in this result. Nevertheless this result may be satisfactory as a preliminary inversion and would certainly provide enough accuracy to spot initial drill holes that would intersect the massive sulphide. We consider this result to be the coarsest discretiza-tion possible where the inversion produces useful results. Any further significant coarsening in the inversion mesh would likely result in instabilities in the inversion due to the lack of accuracy in forward modelling the data. 6.4 C o n c l u s i o n A few of the more critical matters that must be addressed in order to make 3D T E M inver-sion truly practical have been discussed in this chapter. They include efforts in the stages of geophysical data acquisition, T E M system design, survey design and the creation of software utilities for the purposes of T E M processing and inversion. Part of making the process of 3D T E M inversion practical is to make the process efficient and allow inversion results to be pro-93 Chapter 6. On practical and efficient TEM inversion Figure 6.2: a) East-west cross section from original inversion without parallelization and using default parameters for cooling. The view is from south, b) East-west cross section though inverted model with efficient cooling schedule. The view is from south, c) East-west cross section though inverted model with efficient cooling schedule and coarse discretization. The view is from south, d) East-west cross section though 3D conductivity model from drill core testing. The view is from south. 94 Chapter 6. On practical and efficient TEM inversion duced in a reasonable amount of time. The factors affecting the efficiency of the algorithm Eh3DTDinv have been investigated and using the San Nicolas U T E M survey as a example we have shown that it is possible to reduce computation times by a factor of 6 while still achieving high quality results. This factor can be increased to nearly 13 if concessions are made with regards to the resolution and quality of the desired result. This reduction in computation times represents a significant advancement in the process of making T E M inversion practical. 95 Chapter 7 Conclusion and Discussion The purpose of this thesis was 1. To show that T D E M data collected in a field survey can be inverted to recover a 3D distribution of electrical conductivity. 2. To provide a simple, easy to follow, yet highly effective procedure, that can aid in the inversion of T E M data. 3. To show that the inversion of T D E M data can yield substantial benefits for interpreters. 4. To show that the process of inversion of T D E M data is practical and could become routine for many applications. In order to determine the extent to which we have succeeded in these goals, the achievements made will be briefly summarized here. The thesis began by introducing Maxwell's equations, the relevant physical properties and the time domain method of electromagnetic surveying. The theoretical background for forward modelling and inversion of 3D T E M data was discussed. For reasons of brevity the discussion of these subjects was intended to be more of an introduction, rather than a comprehensive treatment. The necessary references were made to literature should further detail be required. A n outline was presented for a simple procedure which can be used to help invert T E M data recorded in the field. In order to illustrate the effectiveness of this procedure an archetypal example was used, that of the U T E M survey over the San Nicolas deposit in Zacatecas, Mexico. The numerous preparations required for inversion of a field data set were addressed in the specific order required by the procedure developed for T E M inversion. It was necessary to develop a background model model based on the gross characteristics of the survey domain. This background model was found to assist with discretization of the transmitter waveform in time. The discretization of the spatial domain was based on several guidelines developed, both here, 96 Chapter 7. Conclusion and Discussion and elsewhere in the literature. The process of reducing the volume using the background model was explained and the utility of the process was demonstrated. Further inversion preparations were considered, such as validation of discretization, and the use of forward modelling exercises to help to determine a tolerable level of modelling and measurement error. The thesis has clearly demonstrated the utility of forward modelling in providing an understanding of the data, and developing some simple interpretations of the data. Interpretations which are often useful in evaluating inversion results. The final stage of preparation before inversion should be to assign estimates of the error in the data based on consideration of both, the modelling error, and an estimate of the measurement error in the data set. It was determined that, before proceeding with inversion of field data, forward modelling and inversion of data generated from synthetic models can be useful. For the San Nicolas deposit a synthetic model was developed which contained a simple conductive block in a more resistive halfspace. The forward modelled data was inverted to aid in the determination of the resolu-tion ability of the San Nicolas U T E M survey. The thesis provided methods by which one could evaluate preliminary inversion results. As a result of evaluations of both synthetic inversions and preliminary inversions of the San Nicolas data set it became clear that the non-uniqueness of the T E M inversion problem was responsible for undesirable artifacts in the recovered mod-els. The solution to the non-uniqueness problem was determined to be a combination of a sensitivity weighting scheme and the use of data from multiple transmitters simultaneously in a single inversion. Inversions of this type provided excellent results from synthetic models and a substantial improvement to the inversion of the San Nicolas U T E M data. Another important stage in an inversion process should generally be the inclusion of a priori information about the local geology in order to further constrain the results. For the inversion of San Nicolas data it was possible to incorporate a priori information about overburden thickness and conductivity. The final inversion of the San Nicolas U T E M data produced a model which fit the data from all three transmitters to a high degree of accuracy and agreed well with the interpreted geology at the site. The results are, to date, the most accurate inversion model produced from any of the numerous E M surveys that have been conducted and interpreted at the site. The result of this inversion is, to our knowledge, the first completely successful application of a 3D T E M inversion scheme to a large scale field data set. 97 Chapter 7. Conclusion and Discussion The thesis continued to discuss the practical matters facing the geophysicist who intends to use 3 D T E M inversion in routine analysis of data. The discussion focused primarily on increasing the efficiency of the inversion algorithm. However, other aspects are necessary to consider including recommendations for survey and instrument design, data format changes, and software development. When one considers efficiency of inversion, one must keep in mind that the goal at a research level is to to be able to complete an inversion in a matter of a day or perhaps a few days. By concentrating on increasing the efficiency of the T E M inversion code, large scale 3 D T E M inversion problems were solved with a degree of high accuracy in less than 5 days. This is a significant achievement considering that, prior to this work, similar problems took nearly 1 month to achieve an acceptable solution. It is clear that we are very close to achieving the research goals for 3 D time domain inversion. For the technology to be useful to industry, however, 3 D inversion results need to be pro-duced in a few hours or perhaps overnight at most. It has been noted that computers tend to increase in speed rapidly and in an amazingly consistent manner. Moore's "Law" states that affordable, readily available computers will double in speed every eighteen to twenty four months. Given that this trend continues to hold we can envision that in the space of four to five years it may be possible to complete high quality inversions, such as those shown here for the San Nicolas survey, in less than twelve hours. This means that the process of 3 D inversion could become a routine processing step for many time domain data sets within the space of just a few years. To conclude, the work done here indicates the utility and practicality of the 3 D T E M inversion, particularly in situations where geology is complex and where conductivity structures tend to vary in three dimensions. There has also been considerable success in our effort to develop a procedure for the inversion of time domain E M data. This procedure should assist in applying the process of 3 D inversion to a wide variety of T E M data types and geophysical applications. We anticipate that the ability for 3 D inversion of T E M data will have an beneficial impact on E M surveying as a whole. This means helping to create new and more efficient survey designs and expanding the range of targets that E M surveys are capable of locating and describing. The results presented here are the first of their kind and represent a significant milestone in the road to making 3 D inversion routinely applicable to T E M data. Eventually we imagine that the procedure developed here will aid in the process of making 3 D time domain 98 Chapter 7. Conclusion and Discussion inversion a routine processing procedure. There are a number of crucial advancements that remain to made in the field of 3D T E M inversion. Before concluding this document it is possible to define some critical future research objectives. First, as mentioned previously it will prove necessary to link the 3D T E M inversion with geological modelling software to facilitate the incorporation of a priori information. Second there are many situations where T E M is used to detect and define highly conductive targets. Theses targets often have conductivity contrasts with their hosts of six or more orders of magnitude. This means that they have extremely long time constants and need to be very finely discretized for inversion to proceed correctly with EH3DTDinv. Much work has been completed on this subject and some good results have been produced from field data however some questions remain. Third, many T E M surveys, such as airborne geophysics, tend to record relatively few data for each transmitter location. This is for reasons of logistic limitations which would be very difficult to surmount. Because the code EH3DTDinv is best used with only a few sources each with many data, airborne T E M surveys are difficult or impossible to invert using today's generation of computers. It may be possible to modify the forward modelling portions of the inversion code to allow highly efficient calculation of the electromagnetic fields for situations where many transmitters are present. Fourth, there are considerable advances that could be made in the area of spatial and temporal discretization. It may be desirable to use unstructured spatial grids to calculate forward models. Though the use of unstructured grids requires considerable programming effort to implement, it could have a profound effect towards reducing the total number of cells in the inversion model. Unstructured grids could also improve problems associated with large aspect ratios of the cells in the padding region of a mesh. These large aspect ratios tend to slow the forward modelling calculations. Additionally it may be possible to divorce the inversion grid from the forward modelling grid. This can result in inverse problems with significantly fewer cells than presently required for E H 3 D T D . Furthermore it may be possible to develop 3D T E M inversions that perform on adaptive grids which refine themselves as the inversion progresses in order to more accurate represent interfaces and objects of interest in the inverted model. All of these potential advancements are either current or future research topics at U B C - G I F and together should help make 3D T E M inversion a widely applicable, routine processing procedure, that yields high quality results and helps to advance the understanding of all time domain electromagnetic data. 99 B i b l i o g r a p h y Atekwana, Estella A, Sauck, William A, & Jr., Douglas D. Werkema. 2000. Investigations of geoelectrical signatures at a hydrocarbon contaminated site. Journal of Applied Geophysics, 44(2-3), 167-180. Bahr, K . 1997. Electrical anisotropy and conductivity distribution functions of fractal random networks and of the crust: the scale effect of connectivity. Geophysical Journal International, 130(3), 649-660. Chen, J , Oldenburg, D W, Haber, Eldad, & Bishop, J . 2002. The down-hole magnetometric resistivty method: A 3-D inversion scheme, unpublished. Commer, Micheal, Hordt, Andreas, Helwig, Stefan L, & Scholl, Carsten. 2003 (September). Three-dimensional inversion of time-domain EM data with highly constrained model complex-ities. Kolloquium Elektromagnetishe Tiefenforshung. Constable, S, & Wiess, C J . 2006. Mapping thin resistors and hydrocarbons with marine E M methods: insights from ID modeling. Geophysics, 71(2), G43-G51. Dennis, J E , & Schnabel, R B. 1996. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SI A M . Farquharson, Colin G, Oldenburg, Douglas W, & Shekhtman, Roman. 2005. Manual for Program EM1DTM. Tech. rept. University of British Columbia Geophysical. Fitterman, David V , & Stewart, Mark T . 1986. Transient Electromagnetic Sounding for Groundwater. Geophysics, 51(4), 995-1005. Fuller, B D, & Ward, S H . 1970. Linear system description of the electrical parameters of rocks. IEEE Transactions on Geoscience Electronics,, (1), 7-18. 100 Chapter 7. Conclusion and Discussion Haber, E , Oldenburg, D W, & Shekhtman, R. 2007. Inversion of Time Domain 3D electro-magnetic data. Geophysical Journal international. Accepted for Publication Haber, Eldad, Ascher, Uri, & Oldenburg, Douglas W. 2004. Inversion of 3D electromagnetic data in frequency and time using an inexact all-at-once approach. Geophysics, 69, 1216-1228. Jones, J W. 2005. Scale-dependent resistivity measurements of Oracle granite. Geophysical Research Letters, 22(11), 1453-1456. Macnae, James. 1985. Manual of Large Loop Interpretation. Lamontagne Geophysics. Napier, Scott, & Oldenburg, Douglas W. 2005 (August). Workflow Elements: Discretization in Space. Tech. rept. University of British Columbia Geophysical Inversion Facility. Newman, Gregory A, & Commer, Micheal. 2005. New advances in three dimensional transient electromagnetic inversion. Geophysical Journal International, 160, 5-32. Nodecal, J , & Wright, S. 1999. Numerical Optimization. Springer. O E D . 1989. electrical conductivity, The Oxford English Dictionary. 2nd edn. Oxford University Press. Parker, Robert L . 1994. Geophysical Inverse Theory. Princeton University Press. Phillips, Nigel. 2001. Geophysical Inversion in an Integrated Exploration Environment: Ex-amples from the San Nicolas deposit. Ph.D. thesis, University of British Columbia. Phillips, Nigel, & Oldenburg, Douglas W. 2005 (August). Workflow Elements: Discretization in Time. Tech. rept. University of British Columbia Geophysical Inversion Facility. Saad, Y . 1996. Iterative Methods for Sparse Linear Systems. PWS publishing company. Telford, W E , Geldart, L P, & Sheriff, R E . 1990. Applied Geophysics. 2 edn. Cambridge University Press. Unsworth, M J, Jones, A G, Wei, W, Marquis, G, Gokarn, S G, Spratt, J E , Booker, P Bedrosian J , Leshou, C, Clarke, G, Li , S, Chanhong, L , Ming, D, Sheng, J , Solon, K , Handong, T , Ledo, J , &: Roberts, B. 2005. Crustal rheology of the Himalaya and Southern Tibet inferred from magnetotelluric data. Nature, 438(7064), 78-81. 101 Chapter 7. Conclusion and Discussion Visser, S, & Shledrake, R F. 1999. UTEM-3 electromagnetic survey: San Nicolas massive sulphide deposit Zacatecas State, Mexico. Tech. rept. S J geophysics Ltd. Wait, J R, & Spies, K P. 1969. Quasi-Static Transient Response of a Conducting Permeable Sphere. Geophysics, 34(5), 789-792. Ward, S, & Hohmann, G. 1987. Electromagnetic theory for geophysical applications, in Elec-tromagnetic methods in applied geophysics, vol. 3, no. 1. Society of Exploration Geophysicists Investigations in Geophysics 3, 131-308. West, G F , Macnae, J C, & Lamontagne, Y . 1984. A Time Domain E M System Measuring the Step Response of the Ground. Geophysics, 26, 1010-1026. Whittall, Kenneth P, & Oldenburg, Douglas W. 1992. Inversion of Magnetotelluric Data for a One-Dimensional Conductivity. Geophysical Monograph Series, no. 5. Society of Exploration Geophysicists. 102
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Practical Inversion of 3D time domain electromagnetic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Practical Inversion of 3D time domain electromagnetic data Napier, Scott 2007
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Practical Inversion of 3D time domain electromagnetic data |
Creator |
Napier, Scott |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | The purpose of this thesis is to provide a simple procedure by which time domain electromagnetic (TEM) data, which has been collected in the field, can be inverted to recover a three dimensional distribution of electrical conductivity. An algorithm, EH3DTDinv, which has been developed in order to invert TEM data in 3D is briefly presented. The algorithm has been previously shown to work well on synthetic data. However, its ability to invert real data collected in field surveys was previously unclear. To illustrate the effectiveness of the inversion algorithm when used in conjunction with an inversion procedure, the example of the UTEM survey at the San Nicolas deposit in Zacatecas Mexico is used. The inversion procedure provides a detailed description of the preparations required for 3D TEM inversion. In order to clarify the individual stages in the procedure, the inversion preparations are outlined as they apply to the San Nicolas survey. Problems of non-uniqueness in the inverted models at San Nicolas are encountered and the solution is determined to be a combination of a sensitivity weighting scheme and the use of the data from multiple transmitters simultaneously in a single inversion. The iterative nature of the inversion procedure is demonstrated through the San Nicolas example. The San Nicolas inversion example culminates in the simultaneous inversion of data from 3 transmitters where a priori information from the regional geology has been incorporated into the result. The 3D inversion result is the first large scale 3D TEM inversion for a mineral deposit and agrees well with the geology and conductivities interpreted from drilling. Practical aspects of inversion are also considered and some recommendations on data formats, software utilities, and survey design are made in order to enhance future results. A brief study into inversion efficiency determined that a significant reduction in computation times can be achieved without loss of accuracy or resolution by fine tuning the cooling schedule of the algorithm. Further reductions in computation times are made through parallelization of the forward modelling process. Finally, yet more reductions in computation times are possible through the coarsening of the spatial and temporal discretizations. However, this final method of computational reduction is shown to come at the cost of accuracy and resolution in the inverted model. The thesis concludes with a discussion of the potential for 3D TEM inversion and brief outline of future advancements in the field. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-03-07 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0053309 |
URI | http://hdl.handle.net/2429/32113 |
Degree |
Master of Science - MSc |
Program |
Oceanography |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2007-0199a.pdf [ 11.72MB ]
- Metadata
- JSON: 831-1.0053309.json
- JSON-LD: 831-1.0053309-ld.json
- RDF/XML (Pretty): 831-1.0053309-rdf.xml
- RDF/JSON: 831-1.0053309-rdf.json
- Turtle: 831-1.0053309-turtle.txt
- N-Triples: 831-1.0053309-rdf-ntriples.txt
- Original Record: 831-1.0053309-source.json
- Full Text
- 831-1.0053309-fulltext.txt
- Citation
- 831-1.0053309.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0053309/manifest