APPROXIMATE SENSITIVITIES FOR THE MULTI-DIMENSIONALELECTROMAGNETIC INVERSE PROBLEMbyColin Glennie FarquharsonB.Sc. (Geophysics), University of Edinburgh, 1990A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF GEOPHYSICS AND ASTRONOMYWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAAugust 1995© Colin Glennie Farquharson, 1995In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)aDepartment of o4’I-tYSuS T)2OiOryThe University of British ColumbiaVancouver, CanadaDate i4LL-r ‘t’1’1DE-6 (2)88)AbstractInversion of measurements from a geophysical electromagnetic survey to produce atwo- or three-dimensional conductivity model of the Earth is computationally demanding.The classical approach is to linearise the inverse problem and iterate towards the solution.A typical, modern version of this approach is described in the initial part of this thesis,and applied to the one-dimensional inversion of time-domain electromagnetic data.One of the most time-consuming parts of a linearised, iterative inversion procedure isthe generation of the Jacobian matrix of sensitivities. In this thesis, I have developed anapproximate method for generating these sensitivities. The approximation is based on theadjoint-equation method in which the sensitivities are obtained by integrating, over anindividual cell, the scalar product of an adjoint field (the adjoint Green’s function) withthe electric field produced by the sources for the geophysical survey. Instead of calculatingthe adjoint field in the multi-dimensional conductivity model, an approximate adjoint fieldis computed in a homogeneous or layered halfspace, or using the Born approximation.This approximate adjoint field is significantly quicker to compute than the true adjointfield and leads to considerable reductions in the time required to generate the Jacobianmatrix of sensitivities. The time-differences were found to be of one or two orders ofmagnitude for the small examples considered in this thesis, and this saving will furtherincrease with the size of the problem.The approximate sensitivities were compared to the accurate sensitivities for twoand three-dimensional conductivity models, and for both artificial sources and the planewave source of magnetotellurics. In all the examples considered, the approximate sensitivities appeared to be sufficiently accurate to allow an iterative inversion algorithm toconverge to an acceptable model. This was emphasised by the successful inversion, usingapproximate sensitivities, of two sets of magnetotelluric data: a synthetic data-set and asub-set of the COPROD2 data.11Table of ContentsAbstract.iiTable of Contents iiiList of Tables viList of Figures viiAcknowledgments ixDedication xChapter 1 Introduction 1Chapter 2 One-Dimensional Inversion of Time-DomainElectromagnetic Data 72.1 Introduction 72.2 The forward problem 82.3 Calculation of the sensitivities 102.4 The inversion algorithm 132.5 Examples 172.5.1 Synthetic data 172.5.2 Field data 222.6 Conclusions 27Chapter 3 Introduction to Approximate Sensitivities 303.1 Multi-dimensional problems 303.2 The need for approximate sensitivities 323.3 Previous forms of approximate sensitivities 343.3.1 The one-dimensional problem 343.3.2 The two-dimensional problem 351113.3.3 The three-dimensional problem.373.4 Conclusions 37Chapter 4 An Approximate Form of Sensitivity for the GeneralElectromagnetic Inverse Problem 394.1 Introduction 394.2 The adjoint-equation method for accurately calculating the sensitivities ... 394.3 An approximate form of sensitivity 43Chapter 5 Approximate Adjoint Fields for the Two-DimensionalProblem 455.1 Introduction 455.2 True adjoint field 455.3 Approximate adjoint field: homogeneous halfspace 475.4 Approximate adjoint field: layered halfspace 505.5 Approximate adjoint field: Born approximation 515.6 Computation times 525.7 Conclusions 53Chapter 6 Approximate Sensitivities for the Two-DimensionalMagnetotelluric Inverse Problem 576.1 Introduction 576.2 Calculation of the sensitivities 576.2.1 Sensitivities for apparent resistivity and phase 576.2.2 Area integration 586.3 Comparison of approximate and accurate sensitivities 606.3.1 Homogeneous halfspace 606.3.2 Two-dimensional model 636.4 Computation times 666.5 Two-dimensional inversion of magnetotelluric data 686.5.1 Synthetic data 68iv6.5.2 COPROD2 data.736.6 Conclusions 79Chapter 7 Approximate Sensitivities for the 25d andThree-Dimensional Inverse Problems 807.1 Introduction 807.2 Calculation of the sensitivities for the 2.5d problem 817.3 Comparison of approximate and accurate 2.5d sensitivities 837.3.1 Homogeneous halfspace 837.3.2 Two-dimensional model 847.4 Computation times for the 2.5d problem 867.5 Calculation of the sensitivities for the three-dimensional problem 887.6 Comparison of approximate and accurate three-dimensional sensitivities 897.6.1 Homogeneous halfspace 897.6.2 Three-dimensional model 927.7 Computation times for the three-dimensional problem 927.8 Conclusions 95Chapter 8 Summary 97References 100Appendices 105A Computation of the Electric and Magnetic Fields in a Layered Halfspacedue to a Dipole Source 105B Adjoint-equation method for the purely two-dimensional inverse problem .. 111C Electric fields due to two-dimensional dipole sources on ahomogeneous halfspace 112D Sensitivities for the magnetotelluric apparent resistivity and phase 118vList of Tables5.1 Example computation times for the various types of two-dimensionaladjoint fields 536.1 Example computation times for the various types of sensitivities for thetwo-dimensional magnetotelluric inverse problem 687.1 Example computation times for the accurate and approximate sensitivitiesfor the 2.5d inverse problem 887.2 Example computation times for accurate and approximate sensitivitiesfor the three-dimensional inverse problem 94viList of Figures2.1 Linear transformations for the id TEM forward problem 102.2 Models produced by the id inversions of synthetic data 182.3 The synthetic data inverted to give the models in Fig. 2.2 192.4 Convergence curves for the id inversion of synthetic TEM data 212.5 A TEM sounding acquired during an environmental survey 232.6 Final models from the inversion of the data in Fig. 2.5 262.7 A second TEM sounding from the environmental survey 282.8 The model obtained from the inversion of the data in Fig. 2.8 295.1 The two-dimensional conductivity model used throughout this thesis 465.2 The mesh corresponding to the model in Fig. 5.1 475.3 The amplitudes of the true and approximate adjoint fields 485.4 The phases of the true and approximate adjoint fields 495.5 The differences between the logarithm of the amplitudes of the adjoint fields .. 555.6 The differences between the phases of the adjoint fields 566.1 Notation used for the bilinear interpolation 596.2 E-polarisation MT sensitivities for the homogeneous halfspace 626.3 11-polarisation MT sensitivities for the homogeneous halfspace 646.4 E-polarisation MT sensitivities for the model in Fig. 5.1 656.5 11-polarisation MT sensitivities for the model in Fig. 5.1 676.6 Synthetic MT data inverted using approximate sensitivities 706.7 The final model from the inversion of the data in Fig. 6.6 726.8 Convergence curves for the inversion of the synthetic MT data 736.9 The sub-set of the COPROD2 data inverted using approximate sensitivities 756.10 The mesh used in the inversion of the COPROD2 data 776.11 The final model produced by the inversion of the COPROD2 data 79vu7.1 2.5d sensitivities for the homogeneous halfspace 857.2 2.5d sensitivities for the model in Fig. 5.1 877.3 Geometry for the three-dimensional example in section 7.6 907.4 3d sensitivities for the homogeneous halfspace 917.5 3d sensitivities for the 3d conductivity model 93A.1 Notation for the horizontally layered conductivity model 106C.1 The coordinate system and geometry for purely two-dimensional problems ... 114viiiAcknowledgementsI would like to thank my supervisor Dr. Doug Oldenburg for his advice and encouragement, and for his numerous shots of enthusiasm when my own was flagging. I am alsograteful to Drs. Rob Ellis and Yaoguo Li for their advice and assistance on all mattersconcerning computers, inverse theory and electromagnetic geophysics. I would also liketo thank the members of my committee, Drs. Garry Clarke, Richard Froese and MattYedlin, for their advice and support.I am grateful to Cliff Candy of Frontier Geosciences, Inc., for providing the fielddata, the well-log measurements and the output of the TEMEX-GL program that appearin Chapter 2. I would also like to thank Drs. Martyn Unsworth and Art Raiche for theuse of their 2.5d and 3d modelling programs.I would like to acknowledge receipt of a Commonwealth Scholarship which enabledme to carry out the work for this thesis, and to do so in such a wonderful part of theworld.Finally, thanks to everyone I have known in the Department of Geophysics andAstronomy for making my time here so very enjoyable.ixTo my parents.xChapter 1IntroductionThere are many geophysical techniques which use the principles of electromagnetismto investigate the electrical conductivity of the Earth. The techniques involve a source, either natural or artificial, that generates electric currents within the Earth. Measurementsare then made of the electromagnetic fields associated with these currents. The behaviourof the currents is dependent upon the spatial variation of the Earth’s conductivity, and soinformation about the conductivity can be inferred from the measurements of the electromagnetic fields. A compendium of present-day techniques is given by Nabighian (1988).Geophysical electromagnetic techniques were pioneered by the mining industry andunderwent considerable development in the period immediately following World War II(Brant 1966). Exploration for mineral deposits remains the main use of these techniquestoday. However, electromagnetic methods are also used to study both the deep andshallow structure of the Earth’s crust. Techniques such as the magnetotelluric methodare used to investigate the large-scale features of the lithosphere (e.g., Jones 1993), andmethods using small, controlled sources play a role in engineering and environmentalstudies of the shallow subsurface (e.g., Heiland 1940; Ward 1990).Before the general accessibility of computers, quantitative interpretation of the measurements obtained from geophysical electromagnetic techniques was limited to the useof “type” or “master” curves (Grant and West 1965; Telford, Geldart & Sheriff 1990).Such curves showed the electromagnetic fields that would be observed over simple Earthmodels for certain sources and were obtained either from analytic formulae if such existedor from scale-model experiments. If a type curve could be found that matched to someextent the observations, then the corresponding model could be considered to resemblethe region of the Earth under investigation.1CHAPTER 1: INTRODUCTION 2With the coming of accessible computer resources, it became possible to computethe electromagnetic fields over more complex models than those used to produce typecurves. Cagniard’s (1953) procedure for calculating the magnetotelluric response of ahalfspace made up of uniform, horizontal layers was adapted by Wu (1968) for the numerical computation of the response for an arbitrary number of layers. Algorithms werealso developed for computing the electromagnetic fields generated in a layered halfspaceby controlled sources (e.g., Ryu, Morrison & Ward 1970; Scriba 1974). As the speedand memory of computers increased, so did the complexity of the models and sources forwhich the electromagnetic fields could be computed. Programs were written to computethe magnetotelluric response over two-dimensional conductivity models (e.g., Madden &Swift 1969; Brewitt-Taylor & Weaver 1976; Rodi 1976; Jupp & Vozoff 1977) and tocompute the electromagnetic fields generated in such models by controlled sources (e.g.,Everett & Edwards 1992; Unsworth, Travis & Chave 1993). These models were describedby large numbers of parameters, usually the conductivities of rectangular cells into whichthe Earth model was divided, and were general enough to represent arbitrary conductivitydistributions. Computer programs have also been written to calculate the electromagneticfields generated in arbitrary three-dimensional models, both for the plane wave source ofmagnetotellurics (e.g., Jones 1974; Mackie, Madden & Wannamaker 1993) and for controlled sources (e.g., Gupta, Raiche & Sugeng 1989; Livelybrooks 1993). MDonald &Agarwal (1994) and Newman & Alumbaugh (1995) have developed three-dimensionalforward-modelling programs specifically designed for the new generation of massively parallel computers. However, it is not yet possible, given generally available computer technology, to carry out forward modelling for the size and complexity of three-dimensionalmodels that one would like.The development of automated interpretation schemes has followed closely behind ateach stage in the evolution of forward-modelling programs. With the ability to computewhat would be measured over an arbitrary conductivity structure, and to quickly computethe solution to matrix equations, it became possible to apply mathematical techniquesto invert the nonlinear relationship between the parameters describing a model and theCHAPTER 1: INTRODUCTION 3measured values of the electromagnetic fields. The first inversion techniques to be applied,and the ones which have proven to be the most enduring, are variations of the Gauss-Newton technique (Gill, Murray & Wright 1981). Wu (1968) and Jupp & Vozoff (1975)applied such techniques to the inversion of magnetotelluric data to obtain horizontally-layered conductivity models. The basis of this collection of techniques is to linearise therelationship between the model parameters, m, and the observations, d(0I):d(0= d[mj + i[m] Sm + R. (1.1)The data that would be measured over the Earth if it were equivalent to the model definedby m are denoted by d[m]. The matrix J is the Jacobian matrix of “sensitivities”, orpartial derivatives, of the data with respect to the model parameters:[mJ = i = 1,...,M, j = 1,...,N, (1.2)where M is the number of observations and N is the number of model parameters.Assuming the remainder R can be neglected, eq. (1.1) can be inverted, using some regularised procedure, to give a perturbation Sm to the model that will (hopefully) bringthe predicted data closer to the observations. Successive iterations will then convergeto a model that minimises the misfit between the observations and the predicted data,usually measured in terms of the 12 norm. This was the goal of the initial applications ofGauss-Newton techniques.Linearised, iterative inversion procedures for geophysical electromagnetic data havesince undergone many refinements. The increased speed of computers now means forward modeffing no longer constrains possible inversion techniques for the one-dimensionalmagnetotelluric problem. It is in the context of this problem that Gauss-Newton inversion procedures for geophysical electromagnetic data have received much of their development (Constable, Parker & Constable 1987; Smith & Booker 1988; Whittall & Oldenburg1992), and are now concerned not only with convergence and robustness, but also with thenonuniqueness of the inverse problem. These sophisticated forms of the Gauss-NewtonCHAPTER 1: INTRODUCTION 4method involve many more model parameters than there are observations, and attemptto find the model that has the minimum value of some norm subject to the constraintthat the observations are reproduced to an appropriate level of misfit.Sophisticated Gauss-Newton procedures have now been used to invert magnetotelluric data to give two- and three-dimensional conductivity models (e.g., deGroot-Hedlin &Constable 1990; Mackie & Madden 1993), and to invert data obtained in controlled-sourcesurveys for two-dimensional models (Unsworth & Oldenburg 1995). However, computerresources restrict the size and complexity of the multi-dimensional problem that canpresently be handled, both because of the time and memory requirements of the forwardmodeffing and because of the time required at each iteration to generate the Jacobianmatrix of sensitivities and to invert the linear system of equations.Many other techniques have been applied to the inversion of geophysical electromagnetic data. Because the one-dimensional magnetotelluric problem is the simplestelectromagnetic inverse problem it has received the most attention. An excellent discussion and comparison of the various techniques applied to this problem is given byWhittall & Oldenburg (1992). The majority of techniques consider the inverse problemas an optimisation problem, as do the Gauss-Newton methods, in which an objectivefunction is to be minimised. The objective function is usually some combination of adata misfit and a model norm. Two such techniques are simulated annealing (Dosso &Oldenburg 1991) and genetic algorithms (Schultz et al. 1993). A considerable amount ofwork has also been done on the theoretical aspects of the one-dimensional magnetotelluric inverse problem, both to investigate the existence and uniqueness of solutions, andto devise means of constructing solutions (e.g., Bailey 1970; Weidelt 1972; Parker 1980).The alternatives to the linearised, iterative Gauss-Newton procedures are generallynot as efficient and so have not yet been considered for the two- and three-dimensionalelectromagnetic inverse problems. And since even the Gauss-Newton algorithms are restricted by current computing power, the only feasible approaches to multi-dimensionalinversion of geophysical electromagnetic data are those involving some form of approximation. Most of the methods that have been developed are based on a linearised procedureCHAPTER 1: INTRODUCTION 5but with one or more of the component calculations carried out using an approximation.For example, Torres-VerdIn & Habashy (1993) make use of an approximate forwardmodeffing algorithm, both for the forward modeffing within their inversion routine andto generate the Jacobian matrix of sensitivities, in their procedure for inverting cross-borehole measurements for a two-dimensional conductivity model. The generalised sub-space technique developed by Oldenburg, McGiffivray & Ellis (1993) enables the matrixequation at each iteration of a Gauss-Newton procedure to be inverted in much less timethan the corresponding exact calculations. Mackie & Madden (1993) and Ellis & Oldenburg (1994) use conjugate gradient techniques to solve the matrix equation at eachiteration in their procedures for the multi-dimensional inversion of magnetotelluric data.These techniques have the additional advantage that the Jacobian matrix of sensitivities need not be explicitly generated, but rather that additional forward modellings arecarried out to compute the product of the Jacobian matrix with the necessary vectors.Smith & Booker (1991) use a Jacobian matrix of approximate sensitivities in their inversion procedure for the two-dimensional magnetotelluric problem. Until it is possible tocarry out the necessary computations for multi-dimensional inverse problems in the sameamount of time as is currently possible for the one-dimensional magnetotelluric problem, Gauss-Newton algorithms using accelerated, approximate ways of both generatingthe Jacobian matrix of sensitivities and solving the linear system of equations at eachiteration offer the most practicable means of constructing solutions to multi-dimensionalproblems.This thesis is concerned with the development of a rapid, approximate method forgenerating the Jacobian matrix of sensitivities needed by any iterative, linearised inversion procedure. In particular, an approximate form of sensitivity is sought that willbe suitable for the multi-dimensional inversion of geophysical electromagnetic measurements, and that will be sufficiently general to apply to any survey configuration. I beginin Chapter 2 by describing a Gauss-Newton procedure for inverting time-domain electromagnetic measurements to recover a one-dimensional conductivity model of the Earth.The purpose of this chapter is two-fold: (1) to illustrate the inversion philosophy andCHAPTER 1: INTRODUCTION 6major components of a Gauss-Newton algorithm unrestricted by computing power, and(2) to develop an inversion program, albeit one-dimensional, for a common controlled-source survey configuration. The contents of Chapter 2 have been published (Farquharson & Oldenburg 1993) and the inversion program is currently being used by a numberof mineral exploration companies.Chapter 3 is an introduction to the subject of approximate sensitivities. In this chapter, I discuss the need for a rapid method of approximating the sensitivities and describethe previous work on this subject. A statement of the approximate form of sensitivitythat I present in this thesis is given in Chapter 4. The approximation is obtained fromthe adjoint-equation method (a derivation of which is included in Chapter 4) by replacingthe true adjoint field with an approximate adjoint field that is considerably quicker tocompute. In Chapter 5, I discuss the possible forms of this approximate adjoint field,and compare the main candidates with the true adjoint field for a purely two-dimensionalexample. The approximate sensitivities for the two-dimensional magnetotelluric problemobtained using the approximate adjoint fields described in Chapter 5 are investigated inChapter 6. The approximate sensitivities are compared to accurate sensitivities for asimple conductivity model. A comparison of their respective computation times is alsocarried out. To conclude Chapter 6, the approximate sensitivities are tested in a Gauss-Newton inversion procedure analogous to that described in Chapter 2. Magnetotelluricdata from a synthetic data-set and a field data-set (a sub-set of the COPROD2 data) areinverted to recover two-dimensional conductivity models.In Chapter 7, I compute approximate sensitivities for the 2.5d and three-dimensionalcontrolled-source problems. Comparisons are made of the approximate and accurate sensitivities for simple conductivity models. The comparisons for the three-dimensionalproblem are limited by computing power, but they are sufficient to illustrate the major differences in the respective computation times for the approximate and accuratesensitivities. Finally, the major conclusions of this thesis are summarised in Chapter 8.Chapter 2One-Dimensional Inversion ofTime-Domain Electromagnetic Data2.1 IntroductionGauss-Newton inversion procedures for geophysical electromagnetic data have undergone much development in the context of the one-dimensional magnetotelluric problem. These procedures now concentrate on obtaining the model that has the minimumvalue of some norm subject to the constraint that the observations are reproduced to anacceptable level. Such a procedure has not yet been applied to the controlled-source survey configurations commonly used in the exploration industry, even for one-dimensionalmodels. Fullagar & Oldenburg (1984) did construct an iterative, linearised procedurefor the one-dimensional inversion of frequency-domain data collected using the horizontal loop source-receiver configuration. However, Fullagar & Oldenburg minimised thenorm of the model perturbation at each iteration rather than the norm of the modelitself. In this chapter, a Gauss-Newton procedure is developed for inverting data from acontrolled-source sounding. The procedure contains aid the features present in the modern approaches to solving the one-dimensional magnetotelluric inverse problem. As such,this chapter also serves to illustrate the ultimate capabilities of inversion procedures formulti-dimensional problems.The geophysical method for which the inversion procedure is developed in this chapter is the time-domain electromagnetic (TEM) method. The TEM method is used extensively in the exploration, geotechnical and environmental applications of geophysics. Areview of the method and its uses is given by Nabighian & Macnae (1991). Typically, astep or ramp turn-off in the current flowing in a rectangular transmitter loop is used toinduce currents in the Earth, and the variation with time of the vertical component of the7CHAPTER 2: 1 d INVERSION OF TEM DATA 8magnetic field, or its time derivative, resulting from these induced currents is measured.These measurements can be at any point on the surface of the Earth, either inside oroutside the transmitter ioop.The conductivity model for the inversion procedure is composed of a large numberof horizontal layers of fixed thickness, and is terminated by a half-space. In general thereare many more layers than observations, so the inverse problem is under-determined.This increases the nonuniqueness of the mathematical solution but allows the model thatminimises a specific norm to be found from the infinity of models that adequately reproduce the data. Suitable choice of the model norm enables the inversion procedure toproduce the model that is the most plausible representation of the region of the Earthunder investigation. Ideally, the model that is produced should exhibit the right “character” (that is smooth or blocky in accordance with the assumed geology), be as close aspossible to the conductivity section obtained from a well-log or a neighbouring sounding(if such is available) and have a minimum amount of structure. This last point is particularly important since arbitrarily complicated structures, which would seem unlikely toresemble the Earth, can suffice as mathematical solutions to the inverse problem. Theaim is to generate a model that contains just enough structure to fit the observations, butno more. The flexibility to generate models of a particular character, and hence producethe most plausible model for a particular geological setting, is a fundamental feature ofthe inversion algorithm developed in this chapter for the TEM method.2.2 The forward problemConsider firstly the forward problem of calculating either the vertical componentof the magnetic field, h(t), or its time derivative, 8h(t)/ot, induced at some point onthe surface of a layered conductivity structure by a step turn-off in the current in arectangular transmitter loop. To exploit the work that has been done on electromagneticmethods in the frequency domain, the analysis is carried out in the frequency domainand only at the very end are the results transformed to the time domain. The notation ofCHAPTER 2: id INVERSION OF TEM DATA 9Ward & Hohmann (1988) is used in which lower-case letters represent fields in the timedomain and upper-case letters represent fields in the frequency domain.The magnetic field, H, due to a rectangular transmitter loop can be evaluated by integrating H due to a horizontal electric dipole around the transmitter ioop (Poddar 1982).The major portion of the forward problem therefore involves calculating H due to a horizontal electric dipole. The computation of the magnetic and electric fields due to a dipolesource on the surface of a homogeneous or layered halfspace is required throughout thisthesis. The method used to compute these fields is described in Appendix A.The method given in Appendix A for computing the magnetic field in a layeredhalfspace due to a horizontal electric dipole source makes use of the Schelkunoff potentials A = Ae and F = Fe (section 4, Ward & Hohmann 1988) where ê, is the unitvector in the z-direction. When computing just the vertical component of the magneticfield, only the potential F is required. In the th layer of the conductivity model (seeFig. A.1), F, the two-dimensional Fourier transform of F, satisfies the ordinary differential equation( — = 0, (2.1)where u = k + k — ew2 + iwp.o-. and e are the magnetic permeability and electricalpermittivity, respectively, of the layer: they are assumed fixed in the inversion processbut can have different values in different layers. cr is the conductivity of the th layer.The tilde indicates the (k, ky,,) domain. The steps for solving this differential equationare given in Appendix A.Once F (kr, k, z,w) has been found for all j, the time-domain values of either h(t)or Oh(t)/Ot for the rectangular transmitter loop are obtained by the sequence of lineartransformations shown in Fig. 2.1. The inverse Hankel transform of F that gives H as afunction of space and frequency for the dipole is numerically evaluated using the digitalfiltering technique of Anderson (1979b). Integration of the dipole response around therectangular transmitter ioop is done using a Romberg integration routine (e.g., Press etal. 1992). Finally, the time-domain values of h are computed from the frequency-domainCHAPTER 2: 1 d INVERSION OF TEM DATA 10(km, k, z, w)dipoleinverse Hankel transformH (x, y, z, Li))di poleIntegration around transmitter loopH (x, y, z, w)loopInverse Fourier transformh (x, y, z, t) or8h (x, y, z, t)/8tloopFigure 2.1 The sequence of integrations used to transform the values of thevector potential due to a horizontal electric dipole as a function of wavenumber and frequency to values of the vertical component of the h field, or itstime derivative, for the rectangular transmitter ioop as a function of space andtime. This sequence of linear transformations is also used to obtain the sensitivities 8h/9cr (z, y, z, t) for the rectangular transmitter loop from the valuesof OF/9o (k,k,z,w) for the dipole.values of H using the routine of Newman, Hohmann & Anderson (1986). This sequenceof linear transformations is also used in the next section to obtain the sensitivity of hfrom the sensitivity of F.2.3 Calculation of the sensitivitiesIn order to construct an iterative, linearised algorithm for inverting values of h(t)or 8h(t)/8t for a horizontally layered conductivity model, a procedure is required forcalculating the Jacobian matrix of sensitivities where the sensitivities are the partialCHAPTER 2: 1 d INVERSION OF TEM DATA 11derivatives of the data with respect to the model parameters. For the problem in thischapter,=, (2.2)where h is the jth observation and o is the conductivity of the j’ layer. Here, and forthe remainder of this chapter, h is used to represent either h or Oh/Ot, depending onwhich is being considered as the data.The method used here to calculate the sensitivities for the one-dimensional inverseproblem is a special case of the general adjoint-equation method for calculating the sensitivities for the multi-dimensional electromagnetic inverse problem. The general methodwill be described in Chapter 4, and will form the basis for the approximate form ofsensitivities developed in this thesis.The layered conductivity structure that is the model for the one-dimensional inverseproblem (see Fig. A.1) can be described in terms of a linear combination:(z) = (z), (2.3)where the basis function is equal to unity within the th layer and zero everywhereelse. The coefficient is equal to the conductivity of the th layer. Using this definitionfor the layered conductivity structure, eq. (2.1) can be rewritten in a form that is validfor all zE(—oo,00):/2 N— u0 — iw crjb) F = 5, (2.4)where u = k + k — jt ew2 and S represents the dipole source in the appropriate layer.Differentiating with respect to 0k (McGifflvray & Oldenburg 1990), making use of thechain rule and realising that S is independent of o, eq. (2.4) becomes=iLL)/Tk.F, (2.5)CHAPTER 2: 1 d INVERSION OF TEM DATA 12which is an inhomogeneous ordinary differential equation for the sensitivity OF/Ook.The boundary conditions are —* 0 as z — ±00. This boundary value problemcan be solved using the adjoint Green’s function method (Lanczos 1961). Hence,= f iwk(z)F(z)Gt(z; dz, (2.6)-where the adjoint Green’s function Gt(z; c) satisfiesN(_u+iwaj)Gt(z;) S(z-) (2.7)andG(z;) —* 0 as z —* +00. (2.8)To construct the adjoint Green’s function, two linearly independent solutions of thehomogeneous form of eq. (2.7) are needed, one of which satisfies the boundary conditionas z — —cc, and the other which satisfies the boundary condition as z —* +co. For theproblem under discussion here, it is only the value of the sensitivity at the surface of theEarth (c = 0 in eq. 2.6) that is of interest. In the region —cc < z = 0, where theconductivity is zero, the adjoint Green’s function has the form exp(u0z). In the region= 0 z < cc, eq. (2.7) is the complex conjugate of eq. (2.4). And, since P —* 0 asz —* +oo, the adjoint Green’s function is proportional to F* in this region. Hence,Gt(z. — fce_, z_0, 29“ ‘‘— lc+F*(z), z<=0.At z = = 0, Gt(z; ) must be continuous, and its derivative with respect to z must bediscontinuous by an amount equal to 1 (Roach 1982). This determines the two coefficients:*()c_ = I*(o)— uo*(0)1(2.10)c+ = - -FI*(O) — uoF*(0) (2.11)CHAPTER 2: 1 d INVERSION OF TEM DATA 13where the prime denotes the derivative.So, using the explicit form of the adjoint Green’s function given above, and remembering that is unity in the th layer and zero everywhere else, eq. (2.6) becomesUi)= j [E(z)j2 dz. (2.12)F (0) —u0F(0) z—zThis expression for the sensitivity links the change in F(k, k, z = 0,w) to the changein the conductivity of the j’ layer. F is known throughout the layers from the forwardmodeffing mentioned in the previous section. The desired sensitivity, Oh/Oo, can beobtained from 8F/Uu by the series of linear transformations shown in Fig. 2.1.2.4 The inversion algorithm(obs) .Consider a set of M observations, h , i = 1, . . . , M, which can either be of thevertical component of the h field, or of its time derivative. The goal of the inversionalgorithm is to find a set of layer conductivities that adequately reproduces these observations. Because the conductivities found in the Earth can vary by orders of magnitude,it is convenient to work with the logarithm of conductivity in the inverse problem. Also,working with logarithms ensures that o is positive. Let m = ln o, and let the vectorm = (m1,. . . , m)” define the model for the inverse problem.The inverse problem is nonunique: if there is one model that adequately reproducesthe finite set of observations, then there is an infinite number of such models. To find aspecific model, a norm of the model is minimised which has the form4’m = II Lm (m — m) j2 (2.13)where Wm is a weighting matrix, is a reference model and . represents the12 norm. The model, m, that minimises 4m will have a character that depends on theparticular choices for the weighting matrix and the reference model. The reference modelCHAPTER 2: id INVERSION OF TEM DATA 14can be used to include any a priori information that may be available about the possibleconductivity structure.The appropriate numerical values of Wm for generating a model of a specific character can be obtained by considering a functional analogous to m for models that arecontinuous functions of depth, for example,= jw(z) T[m(z) _m(z)1I2 dz. (2.14)The operator T can be the identity operator, or the first- or second-order derivative withrespect to z. The function w(z) is an additional weighting function which can be used toenhance or suppress structure over certain depth ranges. The weighting matrix Wm canbe obtained by making m the discrete equivalent of mTo determine whether or not the data produced by the model conductivity structureare sufficiently close to the observations, the following measure of misfit is used:= (h(0b— h[m]) 2 (2.15)where h[m] represents the data computed for the model m. W is a weighting matrixwhich is usually a diagonal matrix whose elements are the reciprocals of the error estimateof each datum. The objective in the inversion is to find a model which gives a misfit,d’ equal to a target misfit For the examples considered in this chapter (andthroughout this thesis), it will be assumed that the measurement errors of the data areunbiased, independent and Gaussian. In such cases, q is equal to the x2 random variable.From the properties of the x2 random variable (Menke 1984; Parker 1994), the expectedvalue of cbd is equal to M, the number of observations. Hence, the final target misfit inthe inverse problem is usually taken as = M.The relationship between the observations and the model parameters is nonlinearand so an iterative method is required to solve the inverse problem. The data misfit afterCHAPTER 2: id INVERSION OF TEM DATA 15the (n + l)th iteration is given by(n+i)= 1 (h(0bs) — h[m’)]) 2 (2.16)To avoid excessive structure developing in the model during early iterations, the targetth • (n) (tar)misfit for the (n + 1) iteration is chosen to be max (,Bq ) where typically 0.1 <3 < 0.5. Also, it is the norm of the model produced by this iteration that is to beminimised:(n+i)= II Wm (m(’ — m) 2 (2.17)The inverse problem to be solved at the (n + l)th iteration is therefore:Minimise (2.18)subject to (n+1) max(,t) c, (2.19)The general approach to solving this optimisation problem is to construct an objective function(n+1)+((fl+1)-q), (2.20)where -y is a Lagrange multiplier. This objective function is then linearised about themodel using a Taylor series expansion, differentiated with respect to the modelperturbation Sm = m(1) — and the resulting derivatives equated to zero. Theset of simultaneous equations is then solved for the model perturbation. (See, for example, Gill, Murray & Wright 1981; Oldenburg, McGiffivray & Ellis 1993.) It is thismethod of solving the linear inverse problem at each iteration that will be used for multidimensional electromagnetic inverse problems elsewhere in this thesis. However, for theproblem considered in this chapter, the method of singular value decomposition (SVD)was used. This method is particularly efficient when the matrix equation to be solvedat each iteration is not too large (less than, say, 200x200), which is the case for theone-dimensional problem treated here.CHAPTER 2: 1 d INVERSION OF TEM DATA 16Consider again the data misfit at the (ri+ i)tI iteration given in eq. (2.16). Considera Taylor series expansion of h[m(’)] about the modelh[m2419= h[m] + J[m] Sm + R, (2.21)where = Oh/8m = uOh/ao (since the model is in terms of the logarithm of thelayer conductivities) is the Jacobian matrix of sensitivities. Assuming the remainder termR can be neglected, substituting eq. (2.21) into eq. (2.16) gives a linearised estimate, çj’,of the true data misfit:= W4 (h(0b — h[m] — [Sm) 2 (n+1) (2.22)By writing Sm explicitly as the difference between the model parameters at two successiveiterations (e.g., Oldenburg 1983; Constable, Parker & Constable 1987), and introducingthe reference model, eq. (2.22) becomes= II W (h(0) — h[m] — 3 + — .j m’’ + 3 m) 2d — (2.23)whered = W (h(0b — h[m] — Jm + Jm()) (2.24)= JW’, (2.25)— m9. (2.26)The linearised inverse problem to be solved at the (m + i)tI iteration can now be statedas:Minimise cbm = Ilth4’W2 (2.27)subject to q = — jth(n+1)II2 = (2.28)CHAPTER 2: id INVERSION OF TEM DATA 17The SVD solution of an under-determined system of equations is the one with thesmallest 12 norm (Wiggins 1972; Parker 1977; Menke 1984; Golub & Van Loan 1989).Writing the SVD of as J = UAVT the solution to the linear inverse problem ineqs. (2.27) and (2.28) is given by= VTA_1UT (2.29)where=s S/(-y s + 1), s is the th singular value of and is the Kroneckerdelta (Wiggins 1972). The Lagrange multiplier y (the same as that in eq. 2.20) couldbe chosen using a line search so that the constraint q q in eq. (2.28) is satisfied.However, it is a solution to the full nonlinear inverse problem that is required. A line(n+1) *search is therefore used to choose -y such that bd = cbd. The value of -y from theprevious iteration is used as an initial estimate for y for the current iteration. Thecorresponding model is computed using eq. (2.29), and a forward modeffing carried outto determine the data misfit, bd, for this model. 7 is then halved, the correspondingmodel computed using eq. (2.29) and a forward modelling carried out to calculate thedata misfit for this new model. 7 is then repeatedly altered, assuming a locally linearrelationship between in cd and in , until either q5, or the smallest possible data misfit,is obtained.2.5 Examples2.5.1 Synthetic dataSynthetic data were generated from the three-layer conductivity structure shown inFig. 2.2. The transmitter loop was a square of side 50 m and the vertical component ofthe h field due to a step turn-off in a 1 A current was calculated 50 m from the centreof the loop. Both the transmitter and receiver were on the surface of the conductivitystructure, and 20 values of h were calculated over the range of delay times shown inFig. 2.3. Gaussian random noise with a standard deviation of 2.5% was added to theCHAPTER 2: id INVERSION OF TEM DATA 1810—1Scfb 102Figure 2.2 The three-layer model (“true”) used to generate the syntheticdata, and the three models produced by the inversion routine. For clarity, theindividual layers in the smallest, flattest and smoothest models are not shown.values of h. The data, with error bars equal to the standard deviation of the addednoise, are shown in Fig. 2.3.The inversion routine was used to obtain three different conductivity models, eachwith a distinct character, that reproduced the synthetic data to the same level of misfit. These three models minimise, in turn, the difference between the model and areference halfspace of 5 x 10—2 S m1, the gradient of the model, and the curvature.For convenience, these models are referred to respectively as the smallest, flattest andsmoothest. The controlling factor for each model is the weighting matrix YLm To obtain the appropriate form of Wm for the smallest model consider eq. (2.14) with w = 1and T = I, where I is the identity operator. Consider also a description of the layeredmodel analagous to that in eq. (2.3). Eq. (2.14) then becomes= J [m — m] (z)2dz (2.30)z0 j=1100 101 102Depth (m)SNI I I I ICHAPTER 2: id INVERSION OF TEM DATA 191 O1 oio ioTime (s)Figure 2.3 The synthetic data, and associated error bars, generated from thethree-layer model in Fig. 2.2. The values of h were calculated 50 m from thecentre of a 50 m x 50 m transmitter ioop. A step turn-off in a 1 A currentwas used as the transmitter current waveform and Gaussian random noise ofstandard deviation 2.5% was incorporated into the data. The continuous curvesrepresent the data predicted by each of the three models shown in Fig. 2.2.N NJ [m_m] [mk—mj(z)k(z)dz (2.31)z=ON= (m — m )2 f dz (2.32)j=1 z=zi_1= (m — )2 , (2.33)because of the particular nature of the basis functions, t is the thickness of theth layer. Comparing eq. (2.33) with eq. (2.13), the weighting matrix required for thesmallest model is therefore= diag (‘T (2.34)CHAPTER 2: 1 d INVERSION OF TEM DATA 20The thickness of the basement halfspace is infinite, and so the final element of Wm isassigned the same value as the previous element. The value of the model in the basementhalfspace will not therefore dominate the model norm.To obtain the weighting matrix for the flattest model consider eq. (2.14) with w = 1and T = d/dz:d 2= J -- [m(z) — m(z)] dz. (2.35)z=O uZAn analogous form of this expression for the discrete case that results in an appropriatemodel norm to minimise isN—i / (ref) (ref) 2i[m. -m ]-[m.--m. 1’=3Lz.‘ ) (2.36)j=i 3Letting Lz equal the vertical distance between the centre of the th layer and the centreof the layer above, comparison of eqs. (2.36) and (2.13) implies that the weighting matrixfor the flattest model is given by—_____1/2 (tl±t2) 1/2 0 . . . 00—(t+t3)/2 (t2+t3)_h/2 00 _(tN_i)uh!2 (tN1)h/20 ... 0 -C(2.37)where the constant C was chosen as 103(tNi)/2. This value is effectively negligiblecompared with the other elements in the matrix, and yet is sufficiently large so that W’can be computed. For the smoothest model, the weighting matrix is the square of thematrix used for the flattest model.The resulting three models are shown in Fig. 2.2. All have 100 layers (i.e., N=100),the thicknesses of which increase exponentially (to the base 1.05 with the thickness ofthe first layer equal to 0.5 m). During the inversion process, 3 (the desired reduction inmisfit at each iteration) was kept fixed at 0.5. This gave a slow but steady convergencetowards the final model. Fig. 2.4 shows the variation of q and m during the courseCHAPTER 2: 1 d INVERSION OF TEM DATA 211010010-110-2 / bio—3 I0 2 4 6 8 10 12ItcrationFigure 2.4 The variation of (a) the data misfit, c1d (b) the model norm, m’and (c) the Lagrange multiplier, 7, during the inversion of the synthetic data inFig. 2.3 to produce the flattest model shown in Fig. 2.2.of the inversion to produce the flattest model. The thirteen iterations shown in Fig. 2.4took 90 minutes on a Sun SparclO workstation.The data calculated from the three models produced by the inversion routine arerepresented by the continuous lines in Fig. 2.3. The values of c5d for the smallest, flattestand smoothest models were all equal to 20 (M = 20). Although these three modelsreproduce the observations to the same level of misfit, they are noticeably different incharacter. These different characters are a direct consequence of the particular modelnorm that was minimised in the inversion. The most obvious differences appear in thedepth ranges that are poorly constrained by the data: below approximately 200 m thesmallest model returns to its reference halfspace of 5 x 10—2 5m’, the flattest modellevels out to achieve zero gradient, and the smoothest model degenerates to a straightline when no longer influenced by the data. Similar behaviour can be seen at shallowdepths above about 10 m. Even in the depth range to which the data are most sensitiveCHAPTER 2: id INVERSION OF TEM DATA 22(10— 200 m, approximately) there are differences in character between the three models.The smallest model manages to follow the block nature of the true model quite closely,whereas the smoothest model smears out the high conductivity layer as much as it canin order to produce the model with the minimum amount of curvature.From Fig. 2.2 it is obvious that the smallest model is the closest to the “true” Earth.However, the remarkable agreement above 10 m depth is somewhat fortuitous, since themodel at these depths is very poorly constrained by the data. And the exact agreementbelow 200 m is only to be expected since the reference model was chosen to have thesame conductivity as the basement halfspace in the true model. For many of the otherdata-sets that the inversion routine was tested upon, it was found that the smallest modeloften contained an unreasonable amount of structure. A model with the least amount ofvariability required to reproduce the observations is generally preferred.2.5.2 Field dataTo test the inversion routine with more realistic data, two sets of field data acquiredduring an environmental study were inverted. Values of the time derivative of the verticalcomponent of the h field were measured at the centre of a square (60m x 60m) transmitterloop using the Geonics Protem system (see Nabighian & Macnae 1991). The data, andtheir assumed error estimates, obtained from soundings at two different locations arepresented in Figs. 2.5 and 2.7.The Protem instrument uses a linear ramp turn-off of length r instead of a pure stepturn-off, which is impossible to generate in practice. This modification of the sourcecurrent waveform can be taken into account when calculating the resultant time decayof the fields by convolving the time-derivative of the linear ramp turn-off with the valuesof Oh/Ot calculated for a pure step-off current source (eq. 1, Asten 1987):ahr 00 Oh(t) = j (u) r’(t - u) du. (2.38)(1r() = —t/rt < —T,—T t 0,t>0.Its time derivative, ‘1’ is therefore a boxcar of length T and height 1/T. The response,Oh”/Ot, for the ramp turn-off was obtained from the response, Oh/Ot, for the step turn-offby numerically evaluating eq. (2.38) using a Romberg integration routine (e.g., Press etModel fi: 191.Model f2: ç = 179.CHAPTER 2: id INVERSION OF TEM DATA 231 31 21010—1>1 —31 —4101Time (s)Figure 2.5 A TEM sounding acquired during an environmental survey. Thetransmitter ioop was 60 m x 60 m, and the receiver was at its centre. The datawere acquired in three overlapping sweeps: 7 its —0.7 ms, 0.1 —2.8 ms and 0.8—70 ms. The observations and assigned error estimates are represented by theerror bars. The continuous curve indicates the data predicted from both modelsfi and f2 in Fig. 2.6.The linear ramp turn-off is represented by(2.39)i05I 1111111 I 1111111 I I I II III1al. 1992).CHAPTER 2: id INVERSION OF TEM DATA 24The observations shown in both Figs. 2.5 and 2.7 were acquired in three overlappingsweeps, the first sweep from 7 s to 0.7 ms, the second from 0.1 to 2.8 ms, and thethird from 0.8 ms onwards. The length of the ramp, r, for these sweeps was 3.5 jis, 35 ,usand 45 ,us respectively. These different ramp times give the data a somewhat sawtoothappearance in the time-range for which the sweeps overlap.The linear ramp turn-off in the transmitter current must also be taken into accountin the inversion process. Because the inversion is being carried out in the time domain,the sensitivities required for a ramp turn-off in the transmitter current can be obtainedby applying the convolution in eq. (2.38) to the sensitivities calculated for the pure stepturn-off.Estimates of the errors were not available for the data-sets shown in Figs. 2.5 and 2.7.Values that seemed reasonable were therefore assigned to the data. For the data inFig. 2.5, errors of 1% were assigned to the measurements before 50 ,us, and errors of 5%were assigned to all other measurements, except over the interval 90 ,us to 3 ms. Inthe centre of this interval (0.3 ms to 0.8 ms) there is an obvious discrepancy betweenmeasurements from the first and second sweeps, suggesting larger errors in these data.Errors of 10% were therefore assigned to the whole interval (90 ,us to 3 ms) over whichthere was overlapping of the three sweeps, although the discrepancy between the firstand second sweeps in the first part of the interval was subsequently found to be due tothe different ramp times for the two sweeps. For the data in Fig. 2.7, errors of 1% wereassigned to the measurements before 80is, errors of 5% to the measurements between 80isand 5 ms, and errors of 10% to those after 5 ms.Fig. 2.6 shows two inversion results for the data in Fig. 2.5. Model fl, represented bythe solid line, was obtained using the same weighting matrix as for the flattest model inSection 5.2.1. Model f2, represented by the dotted line, was obtained using the weightingCHAPTER 2: id INVERSION OF TEM DATA 25matrix—1 1 0 ... 00—11 0=... (2.40)—1 10 0-Cwhere C = iO. Both models have 200 layers increasing exponentially in thickness (tothe base 1.03, with the first layer of thickness 0.1 m). Because of this increasing layerthickness, the model norm constructed using Wm above is essentially a discretised versionof= x: [:]2donz = jz [d(lnu)]2 (2.41)(Remember that m = ln u.) The model norm constructed using Wm in eq. (2.31) is adiscretised version of= f° [d(ln)]2 (2.42)The dashed line in Fig. 2.6 corresponds to the result obtained from a parameter estimationprogram (TEMEX—GL), based on Anderson’s (1979a) routine, which was restricted tocontain seven layers.The values of the misfit, q, for models fi and f2 were 191 and 179, respectively.The value of d for the seven-layer model was 983. Models fl and f2 produce nearlyidentical fits to the observations, as illustrated in Fig. 2.5. However, the two valuesof misfit quoted above are still significantly larger that the expected value of 56 (thenumber of observations). No model could be found which gave a smaller misfit thanthat for model f2. This suggests that the error estimates that were assigned to the dataare smaller than the true errors in the measurements. This seems entirely plausible,especially for the late-time measurements in each data-set, where errors of 30 or 40% aremore likely to be realistic.The similarity between the two models in Fig. 2.6, despite the different nature ofthe weighting matrices used to produce them, emphasises the robustness of an inversionFigure 2.6 fi and f2 are the two versions of the flattest model obtained fromthe inversion of the data in Fig. 2.5. The best-fitting seven-layer model producedby a parameter estimation routine (TEMEX—GL) is also shown.routine which looks for a minimum-structure model. Since models fi and f2 both containjust enough structure to fit the data, it is not surprising that they agree on the featuresthat are required by the data. This provides confidence that the particular sequence ofconductive and resistive layers to a depth of about 300 m is present in the real Earth.Below 300 m the models are less strongly constrained by the observations, and sothe particular form of the weighting matrix begins to dominates their behaviour. Theweighting term z that appears in eq. (2.41) but not in eq. (2.42) causes structure inmodel f2 to be suppressed at these depths compared with model fi. However, the factthat both models fi and f2 do show increases in conductivity around 500 m depth, ratherthan levelling off, might indicate the presence of a good conductor on the extreme limitof penetration of this sounding. When this increase in conductivity was removed fromCHAPTER 2: id INVERSION OF TEM DATA 26Solid — model fi,Dotted — model f2,Dashed — parameter estimation.1011 O10—11 —2b101 100liii I I 111111 I 1111111 I [____I__I_IIIIII101 1 2Depth (m)1 3CHAPTER 2: id INVERSION OF TEM DATA 27model if, the misfit did increase to a value of 219, and the predicted data from thismodified model did not fit the last two data points in Fig. 2.5 as well as those formodel fi. Similarly, removing the final increase in conductivity from model f2 resultedin an increase in the misfit to 194, suggesting that this feature might be required bythe data. At shallow depths, the different weighting matrices lead to structure beingenhanced in model f2 relative to model fi.As a final example, the data from the second sounding (see Fig. 2.7) were inverted.The weighting matrix was the same as the one used to produce model f2 in Fig. 2.6.The resulting model is represented by the dotted line in Fig. 2.8. It is again made up of200 layers whose thicknesses increase exponentially with depth. The solid line in Fig. 2.7represents the data predicted from this model. The corresponding value of the misfitwas cd = 1216. The assigned error estimates again seem to be too small and it was notpossible to obtain a misfit close to the expected value of 40. The solid line in Fig. 2.8represents well-log measurements of conductivity obtained from a borehole 70 m awayfrom the location of the sounding. There is good agreement between the model and thewell-log measurements in the depth range 10— 150 m, especially on the location of theupper regions of the two conductive zones around 35 and 90 m. The character of thewell-log, comprising transition zones rather than a few layers of constant conductivity, isalso reproduced by the model in Fig. 2.8. This type of model, constructed by minimisingthe 12 norm of its gradient, therefore seems particularly suited to this geological setting.2.6 ConclusionsIn this chapter, I have presented an inversion algorithm that generates a layeredconductivity structure from measurements, at any point on the surface of the Earth,of the time-variation of h, or Oh/8t, induced by a rectangular transmitter loop. Theinversion algorithm is a Gauss-Newton algorithm that minimises a norm of the modelsubject to the constraint that the observations are reproduced to an appropriate level.The form of the model norm to be minimised is flexible, and can be constructed toenable the algorithm to find the most plausible model from the infinity of models thatModel:‘kd = 1216.I 1111111 I 1111111 I IioTimeICHAPTER 2: 1 d INVERSION OF TEM DATA 281 31 210110010—11 —21 310-2(s)Figure 2.7 A second sounding from an environmental survey. The surveygeometry and the three measurement sweeps are the same as for Fig. 2.5. Theerror bars indicate the observations and assigned error estimates, and the continuous curve represents the predicted data from the flattest model producedby the inversion and shown in Fig. 2.8.adequately reproduce the observations. Which model is the most plausible will dependon the geological setting and the amount of prior knowledge, and may be the one withthe least amount of structure, or the one that is the closest to some reference model thatrepresents the preconceived image of the region under investigation.The algorithm described in this chapter is robust: it does not become unstable inthe presence of noisy data. In addition, the combination of (1) finding a model thathas a minimum amount of structure, and (2) requiring that the desired misfit at eachiteration is only gradually decreased, leads to a procedure that consistently finds a modelFigure 2.8 The flattest model (dotted line) obtained from the inversion ofthe data in Fig. 2.7. The solid line represents the values of the conductivitymeasured in a borehole 70 m from the observation location.CHAPTER 2: id INVERSION OF TEM DATA 29Dotted — model,Solid — well—log.10°10_i1 —2(I)b10—1 100NH I 1111111 I 11111111 I I I I I III101Depth (m)1 2that reproduces the observations and is a plausible representation of the Earth. The finalmodel produced by the algorithm is also insensitive to the particular starting model used.Finally, the solution to the one-dimensional controlled-source inverse problem described in this chapter is not restricted by computing power unlike the multi-dimensionalproblems considered later in this thesis. All the calculations are therefore exact, includingthe generation of the Jacobian matrix of sensitivities.Chapter 3Introduction to Approximate Sensitivities3.1 Multi-dimensional problemsBefore I start to discuss approximate sensitivities, I shall summarise some conceptsrelevent to multi-dimensional geophysical electromagnetic inverse problems.The source for the magnetotelluric method is assumed to be a plane wave impingingon the surface of the Earth. Because of this simple source, the magnetotelluric method isthe one for which the forward and inverse problems can be most easily carried out. Themeasured quantities in a magnetotelluric experiment are the horizontal components ofthe electric and magnetic fields. Since the details of the plane-wave source are unknown,the ratio of orthogonal components of the electric and magnetic fields is calculated (inthe frequency domain). An apparent resistivity and phase are obtained from this ratio(see eq. D.1), and it is these quantities that are considered as the data acquired in theexperiment. To investigate the two- or three-dimensional structure of the Earth, thesedata have to be collected at many different locations on the surface of the Earth as wellas at many frequencies.For a two-dimensional conductivity model, the electric and magnetic fields inducedin the model by a plane-wave source divide into two decoupled modes. One mode involvesthe component of the electric field parallel to the strike direction (i.e., the direction inwhich the model is invariant) and the two components of the magnetic field perpendicularto the strike direction. This mode is known as the “E-polarisation” mode since the onlycomponent of the electric field within the model is parallel to the strike direction. It isalso referred to as the “transverse electric” (TE) mode since the only component of theelectric field in the model is perpendicular (or transverse) to the vertical direction (seesection 4.1.1, Weaver 1994). The other mode involves the component of the magnetic field30CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 31parallel to the strike direction and the components of the electric field perpendicular to thestrike direction. It is therefore called the “H-polarisation” or “transverse magnetic” (TM)mode since the only component of the magnetic field in the model is parallel to the strikedirection (and transverse to the vertical direction). This decoupling into two modes alsoapplies to any situation in which the source is invariant in the strike direction of themodel.In the exploration industry, controlled-source methods are used more often than themagnetotelluric method. A typical controlled-source method has either a current loopor a grounded line current as a source, and involves measurements of the frequency- ortime-dependence of the electric or magnetic fields. When considering a two- or three-dimensional target region, these measurements are either collected at many differentlocations relative to a stationary source, or at a fixed distance from the source for manydifferent source positions, or a combination of the two.In this thesis, the term “forward modeffing” is used to refer to the process of computing the values of the measurements that would have been obtained if the geophysicalexperiment had been carried out over the mathematical conductivity model rather thanover the Earth itself. One forward modeffing produces values for the measured quantitiesat every observation location for every source location. I shall further specify that oneforward modeffing, as considered in this thesis, corresponds to only one frequency or delay time: if measurements were collected at ten frequencies then ten forward mode]iingswould be required to compute values for all the quantities measured in the experiment.I shall also assume that, in the process of computing the predicted data for a particularfrequency or delay time, a forward modeffing produces the values of the electric fieldeverywhere in the conductivity model.As a final point, it is worth noting the correspondence between “sensitivities” and“Fréchet derivatives”. Sensitivities are obtained by first parameterising the conductivitymodel and then differentiating the equations that define the forward problem with respectto the model parameters (see sections 2.3 and 4.2). An alternative approach is to considerCHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 32a general perturbation to the model before any parameterisation is introduced. Theresulting effect on the observed quantity can be expressed as (Parker 1994):d [m(r) + Sm(r)] = d {m(r)] + (D(r),Sm(r)) + R, (3.1)where d represents the observed quantity and (.,.) represents the inner product. If theremainder R is such that jR/j6mj — 0 as — 0, then D is the Fréchet derivativeof d. If the model parameterisation is now introduced such that m(r) = m&(r)and Sm(r)=Smib(r), eq. (3.1) becomesd [m+ Sm] = d [m] + (D(r),(r)) 5m + R, (3.2)where m = (m1,... , m)’. For the problems considered in this thesis, the inner product (•,) is the volume integral over the whole of three-dimensional space. Hence,(D(r),(r)) = jDi(r)j(r)dv. (3.3)Comparing eq. (3.2) with eq. (1.1), it can be seen that the sensitivity is equivalent to thevolume integral of the product of the Fréchet derivative and the basis function:I= I D(r) (r) dv. (3.4)um j3.2 The need for approximate sensitivitiesThe Gauss-Newton algorithm described in Chapter 2 is the inversion methodologythat one would like to use to solve the two- and three-dimensional electromagnetic inverse problems. However, the computations required by some of the components ofthis algorithm become very time-consuming when applied to two- and especially threedimensional problems. The three most computationally intensive components are (1) theforward modeffing, (2) the generation of the sensitivity matrix and (3) the solution ofCHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 33the linear system of equations. Each of these operations has to be performed at leastonce for each iteration. Highly sophisticated and efficient forward modeffing schemes arenow being developed for two- and three-dimensional conductivity models (e.g., Everett &Edwards 1992; Unsworth, Travis & Chave 1993; Mackie, Madden & Wannamaker 1993;Livelybrooks 1993; MDona1d & Agarwal 1994; Newman & Alumbaugh 1995). The inversion of the linear system of equations has also received attention, and techniques suchas a generalised subspace method can be used to speed up the solution of the matrixequation without hindering the progress of the iterative inversion scheme (Oldenburg,WGiffivray & Ellis 1993; Oldenburg & Li 1994). The generation of the Jacobian matrixof sensitivities, however, has received little attention.The three main methods of calculating the sensitivities for the general electromagnetic inverse problem are the brute-force or perturbation method, the sensitivity-equationmethod and the adjoint-equation method (McGiffivray & Oldenburg 1990). The computation times for these methods are roughly equivalent to N x Mf, N x Mf and M0 x Mfforward modeffings respectively where N is the number of model parameters, Mf is thenumber of frequencies (or delay times) and M0 is the number of observation locations(Mf x M0 = M). The actual time taken by each of these three methods will depend ontheir implementation for a particular problem. It is sometimes possible to exploit featuresof a forward modeffing algorithm to reduce the time needed to generate the sensitivities.For example, Oristaglio & Worthington (1980) make use of the already-factored matrixfrom their finite-difference forward-modeffing program when calculating the sensitivitiesusing the sensitivity-equation method. However, the above estimates of the computationtimes for the three methods of calculating the sensitivities are reasonable when constructing an inversion scheme around an existing forward-modeffing program, and are a goodindication of how the computation times for the sensitivities are affected by the numberof model parameters and data.For surveys over two- or three-dimensional geological structures, the typical numberof measurements is of the order of one or ten thousand. The number of parameters required to adequately describe any model used in the inversion of such measurements willCHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 34also be of the order of one, ten or a hundred thousand. Given that forward-modellingprograms for such two- or three-dimensional conductivity structures currently have computation times of at least several minutes, it can be seen that calculating the sensitivitiesusing any of the accurate methods described above will be very time-consuming, if notprohibitive. It therefore becomes necessary to consider approximate methods for generating the Jacobian matrix of sensitivities if multi-dimensional electromagnetic inverseproblems are to be tackled.3.3 Previous forms of approximate sensitivities3.3.1 The one-dimensional problemThere has been one investigation of approximate sensitivities for the one-dimensionalelectromagnetic inverse problem. Boerner & Holladay (1990) demonstrated that thesensitivities with respect to the layer conductivities do not depend strongly on the modelwhen components of the electric or magnetic fields associated with only horizontally-flowing currents in a horizontally-layered model are considered. They illustrated this byconsidering a horizontal electric dipole source and measurements of the resulting verticalcomponent of the magnetic field, with both the source and measurement location onthe surface of the Earth. They computed the sensitivities of these measurements withrespect to the conductivity of a thin test layer for three models: a homogeneous halfspaceof conductivity 0.01 S m1, a halfspace of 0.01 S m1 in which a conductive layer of0.05 S m1 was embedded, and another halfspace of 0.01 S m1 in which a resistive layerof 0.002 S m1 was embedded. The plots of the sensitivities as functions of the depth ofthe test layer illustrate the weak dependence of the sensitivities on the model: all threecurves are smooth and continuous across the boundaries of the conductive or resistivelayer, they all exhibit the same general characteristics, and differ at most by 10%.Boerner & Holladay call the sensitivity with respect to the conductivity the “incremental” sensitivity. The sensitivity with respect to the logarithm of the conductivity,CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 35what Boerner & Holladay call the “fractional” sensitivity, is obtained by scaling the incremental sensitivity by the conductivity of the layer: c9/O(ln u) = o O/8o. This is thesensitivity that is applicable when the model is defined in terms of the logarithm of theconductivity, as is usually the case. Because of the weak dependence of the incrementalsensitivities on the structure in the model, the dominant influence on the fractional sensitivities is the scaling by the layer conductivities. Boerner & Holladay exploited this byinverting synthetic data using a simple, linearised procedure in which the fractional sensitivities were approximated by the incremental sensitivities calculated for a homogeneoushalfspace and scaled by the conductivity of the layers. They found that the convergenceof the inversion procedure was not hindered by the use of these approximate sensitivities,and that in some instances the rate of convergence was increased.Boerner & Holladay state that the weak dependence of the incremental sensitivitieson the conductivity model is not to be expected for measurements of the electric andmagnetic fields associated with current flow across the boundaries of the layers. However,the plots of Boerner & West (1989) show that, for a situation in which currents flow acrossthe layer boundaries, the fractional sensitivities are also strongly dependent on the scalingby the layer conductivities, and that the fractional sensitivities for a one-dimensionalmodel can be approximated by scaling the incremental sensitivities for a homogeneoushalfspace by the conductivities of the layers. (The situation Boerner & West consideredwas that of a grounded electric dipole source and measurements of the correspondingelectric field, although calculated in the wavenumber domain, at zero frequency and witha geometrical factor removed.)3.3.2 The two-dimensional problemApproximate sensitivities have been used to reduce the time required to invert magnetotelluric data for two-dimensional conductivity models. Oldenburg & Ellis (1993) andEffis, Farquharson & Oldenburg (1993) used the sensitivities from the one-dimensionalmagnetotelluric inverse problem as approximate sensitivities for the two-dimensionalCHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 36problem. For each observation location, a horizontally-layered conductivity structureis imagined in which the conductivities of the layers are equal to the conductivities ofthe cells in the two-dimensional conductivity model immediately beneath the observationlocation. One-dimensional sensitivities are calculated with respect to the layer conductivities. These sensitivities are then used in the two-dimensional inversion as the sensitivitiesfor the cells immediately below the observation location.It was found that iterative inversion procedures using these approximate sensitivitiesconverged to acceptable models when the data were either for the E- or il-polarisation,or determinant averages of the data for the two modes. These approximate sensitivitieswere not as successful, however, when used to jointly invert observations from both theE- and il-polarisations, especially when there were significant differences between thedata from the two modes. This is not surprising since the decoupling into two distinctmodes does not exist in the one-dimensional magnetotelluric problem, and so the one-dimensional sensitivities used as approximate sensitivities in the joint inversion were thesame for both modes.Smith & Booker (1991) developed a slightly better form of approximate sensitivitiesfor the two-dimensional magnetotelluric problem. They assumed that the horizontalgradient of the electric field in the two-dimensional conductivity model was small relativeto its vertical gradient. The resulting expressions for the approximate sensitivities arethe same as the expressions for the one-dimensional sensitivities used by Oldenburg &Effis and Ellis et al., but the electric field in the two-dimensional conductivity model isused rather than the electric field in the layered model imagined to exist beneath eachobservation location. Smith & Booker’s approximate sensitivities are therefore differentfor the E- and H-polarisations, and through the electric field from the two-dimensionalmodel, a small amount of information about the two-dimensional structure of the modelis incorporated in the sensitivities.The advantages of the above two forms of approximate sensitivity are that theyare quick to compute and are sufficiently accurate to enable, in many cases, iterative,linearised inversion schemes to converge to an acceptable model. The “Rapid RelaxationCHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 37Inversion” program of Smith & Booker (1991) is routinely used to invert data fromregional magnetotelluric surveys.The disadvantage of the above forms of approximate sensitivity is that they lead tosparse Jacobian matrices: the approximate sensitivity for a measurement is only definedwith respect to the conductivities of the cells immediately below the location of the measurement. There is no information, therefore, in these approximate sensitivities about thedependence of a measurement on the cells that are close to the location of measurementbut not immediately beneath it. Another consequence is that the discretisation of themodel is restricted to have each column of cells below an observation location.3.3.3 The three-dimensional problemIn principle, the approximate sensitivities of Smith & Booker (1991) described insection 3.3.2 can be extended to the three-dimensional magnetotelluric problem. Also,attempts have been made to calculate approximate sensitivities using the Born approximation (see Hohmann & Raiche 1988). However, the amount of work that has been doneon approximate sensitivities for three-dimensional problems is small.3.4 ConclusionsSections 3.3.1 and 3.3.2 suggest that approximate sensitivities need not be extremelyaccurate to enable a linearised, iterative inversion procedure to converge to an acceptablemodel. The work of Boerner & Holladay (1990) shows that the fractional sensitivities areheavily dependent on the scaling by the conductivity structure of the model, rather thanthe specific features of the incremental sensitivities. A moderate approximation of theincremental sensitivities, when scaled by the conductivities of the model, should thereforelead to a good approximation of the fractional sensitivities. The relative success of theapproximate sensitivities described in section 3.3.2, despite their somewhat poor levelof approximation, strengthens the argument that approximate sensitivities need not behighly accurate to be useful.CHAPTER 3: INTRODUCTION TO APPROXIMATE SENSITIVITIES 38It is also worth noting that in an iterative Gauss-Newton procedure for solving thenon-linear inverse problem (see eq. 2.20 and the surrounding description), the Jacobianmatrix of sensitivities is evaluated at the current model: J[m(’)]. Consider for the moment that the need for an iterative procedure has not yet been recognised. To find themodel that minimises the objective function, therefore, the objective function would simply be differentiated with respect to the model parameters and the resulting expressionsequated to zero. This would produce a matrix equation for the model parameters, m(s01),that minimise the objective function. However, this matrix equation involves the Jacobian matrix evaluated for this unknown model: J[m(s01)]. The matrix equation istherefore non-linear and intractable. (See the section “The Nonlinear Problem” of Constable, Parker & Constable 1987). Comparing this to the iterative approach, [[m(’)]can be thought of as an approximation to J[m(s01)], with only the dominant features ofJ[m(’)] influencing the direction in which the iterative procedure moves through modelspace, at least for the early iterations. This further suggests that approximate sensitivities might be effective in enabling an iterative, linearised inversion algorithm to convergeto an acceptable solution.It is the aim of this thesis to develop a general form of approximate sensitivity thatwill be of use in the solution of any electromagnetic inverse problem, especially thoseinvolving controlled sources.Chapter 4An Approximate Form of Sensitivity for theGeneral Electromagnetic Inverse Problem4.1 IntroductionThe three main methods of accurately calculating the sensitivities are the brute-force or perturbation method, the sensitivity-equation method and the adjoint-equationmethod (McGillivray & Oldenburg 1990). The computation times of these methods areroughly equal to N x Mf, N x Mf and M0 x Mf forward modellings respectively. (N isthe number of model parameters, Mf is the number of frequencies and M0 is the numberof observation locations). The inversion procedure envisioned for the multi-dimensionalproblem, just as for the procedure described in Chapter 2 for the one-dimensional problem, is one in which there are many more more model parameters than observations.The adjoint-equation method is therefore the most efficient method of calculating thesensitivities for such an under-determined problem. Because it is the natural methodfor an under-determined problem, and because it can be readily modified to give a useful approximation, the adjoint-equation method is the basis for the approximate formof sensitivity presented in this thesis. An abridged derivation and description of theadjoint-equation method are given below.4.2 The adjoint-equation method for accurately calculating the sensitivitiesThe adjoint-equation method has previously been applied to electromagnetic inverseproblems by, for example, Weidelt (1975) and Madden (1990). Also, the method used inChapter 2 to calculate the sensitivities for the one-dimensional electromagnetic inverseproblem considered in that chapter is a special case of the adjoint-equation method.39CHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 40The derivation presented here follows that of McGiffivray et al. (1994). Consider a domain, D, whose conductivity, o, varies only as a function of position, and whose electrical permittivity, e, and magnetic susceptibility, t, are constant. Consider also somegeneral source which can be defined in terms of electric and magnetic current densities Jand m = —iWILM. This description of the source terms is not that of McGifflvray etal., but that used by Ward & Hohmann (1988) where M is the magnetisation. Thisdescription is consistent with Appendices A and C in which a unit electric dipole sourcecorresponds to eI = 1 and a unit magnetic dipole source corresponds to M = 1. Theelectric and magnetic fields, E and H respectively, that are generated in the domain aregiven by the frequency-domain Maxwell’s equations:V xE = iWUH+Jm (4.1)VXH(O+W6)E+Je. (4.2)The electric and magnetic fields satisfy the general boundary conditionso(flxU) + ,13(ñxñxVxU) = Q (4.3)on the boundary, ÔD, of the domain D. U represents either E or H, c and /3 areconstants, ii is the unit normal to OD and Q is the appropriate electric or magneticsurface current density.For the purposes of the inversion, the conductivity is represented as a finite linearcombination of suitable basis functions:(r) = (4.4)where are the basis functions and are the coefficients. Substituting this representation of the conductivity into eqs. (4.1) and (4.2), and differentiating these equations withCHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 41respect to o, one obtainsOE . 8HV x = —iw-— (4.5)00k811 . OEV x = (tr + zwe)--— + ct’kE. (4.6)The partial derivatives OE/ôok and OH/ôok are the sensitivities of the electric and magnetic fields with respect to the coefficient ok. These sensitivities satisfy the homogeneousboundary conditionsa (ñx) + (ñxñxVx) = 0. (4.7)To arrive at a numerical method for calculating the sensitivities, a suitable solutionto eqs. (4.5) and (4.6) has to be found. With this in mind, consider an auxiliary problemin which auxiliary electric and magnetic fields, Et and Ht, are generated in the conductivity, o, by as yet undefined electric and magnetic sources J and J. These auxiliaryfields satisfy the following version of Maxwell’s equations:V x = —iw1uH + J (4.8)V x Ht= (o. + iw6)Et + J. (4.9)The auxiliary fields satisfy the homogeneous boundary conditionsa(ñxU) + ,Bt(nxnxVxUt) = 0. (4.10)If the boundary 9D extends to infinity, at and ,B need not be the same as a and /3. Ifthe domain D is finite, however, a = a and 13t = /3 for the analysis that follows.Using the vector identityV.(AxB)=B.(VxA)—A.(VxB) (4.11)CHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 42eqs. (4.5) and (4.6) can be combined with eqs. (4.8) and (4.9) to give.(t x — x Hf) = J. + J . — Et . E&k. (4.12)Integrating this expression over the domain D and using the divergence theorem givesF / OH UEI tEx——---——xHindsJ8D k 0k I= ID [Jt. 0(7k + 0o — E E k] dv. (4.13)It can be shown that, for the boundary conditions specified above, the left-hand side ofeq. (4.13) is equal to zero. Hence,ID (Jt. + jt. -) dv = ID E k dv. (4.14)Eq. (4.14) is the main result from which the sensitivities for the electric and magnetic fields can be obtained for a particular problem. It requires appropriate choice ofsources, J and J, for the auxiliary problem. For example, to obtain the sensitivity ofthe x-component of the electric field at an observation location r0, choose J =and J = 0. Substituting these expressions into eq. (4.14) gives= J E E bk dv. (4.15)0k r0 DFor this example, the auxiliary electric field, E, is now defined as the electric field in thedomain, D, due to an x-directed unit electric dipole at the location, r0, of the observationsof E. E is the electric field due to the particular source (described by e and m) usedin the geophysical survey. To obtain the sensitivity of the x-component of the magneticfield at r0, set J = 0 and J = S(r — r0) ê in eq. (4.14). The sensitivity is then givenby =J E .Ek dv, (4.16)DCHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 43where E is now the electric field in the domain due to an x-directed magnetic dipole ofstrength (—iwt) at r0.Eq. (4.14) can be obtained by a formal Green’s function solution of eqs. (4.5) and (4.6)(see the references given at the start of this section). The auxiliary electric field, E, inthis section is equivalent to the complex conjugate of the appropriate adjoint Green’sfunction for eqs. (4.5) and (4.6). I shall therefore use the term “adjoint field” rather than“auxiliary electric field” to refer to E.4.3 An approximate form of sensitivityTo calculate the sensitivities using the adjoint-equation method outlined in the previous section, the following are needed: (1) the electric field, E, in the domain due to thesource used in the geophysical experiment, (2) the adjoint field, E, in the domain due tothe appropriate dipole source at the observation location, and (3) the evaluation of thevolume integral (weighted by the basis function) of the scalar product of these two fields.The computation time required to numerically evaluate the volume integral is negligiblecompared to that for the other parts of the procedure. Also, in an iterative inversionprocedure the electric field, E, will already have been computed in order to calculate thedata misfit resulting from the previous iteration. No additional work therefore has tobe done to fulfil the first requirement above. However, for the second requirement, theadjoint field has to be explicitly computed in the conductivity model for a dipole sourceat each observation location in turn. It is this that accounts for the vast majority of thecomputation time for this method, and hence the statement that the computation timefor the adjoint-equation method is equivalent to that for M0 x Mf forward modellings.In this thesis, I present an approximate form of sensitivity that is much quicker tocompute than the adjoint-equation method: instead of computing the adjoint field, E,in the multi-dimensional conductivity model, an approximate adjoint field is computedthat is considerably quicker to obtain. The types of approximate adjoint fields given inthis thesis are the electric fields in a homogeneous halfspace, or in a horizontally-layeredCHAPTER 4: AN APPROXIMATE FORM OF SENSITIVITY 44halfspace, which are generated by the appropriate dipole source. The Born approximationof the adjoint field is also considered.Computing an approximate adjoint field will be quicker than computing the trueadjoint field in the multi-dimensional conductivity model. Repeating this calculation forthe dipole source at each observation location will lead to considerable time-savings inthe generation of the Jacobian matrix. The approximate adjoint field will obviously notcontain all the features of the true adjoint field in the multi-dimensional model. However, as I shall illustrate in the next chapter, the true adjoint field is dominated by itsdecay away from the dipole source. Suitable choice of conductivities for the homogeneousor layered halfspace result in an approximate adjoint field that has the same dominantbehaviour. Such an approximate adjoint field, when combined in the volume integrationwith the electric field, E, from the forward modelling (which is exact) produces approximate sensitivities that are sufficiently accurate to allow an iterative inversion procedureto converge to the desired result.Furthermore, the adjoint-equation method of exactly calculating the sensitivities iscompletely general. The approximate form of sensitivities presented here can thereforebe used in any electromagnetic inverse problem irrespective of source-receiver geometry.In the next Chapter, I shall investigate the forms of the approximate adjoint fieldfor the two-dimensional electromagnetic inverse problem and emphasise the comparisonsbetween the approximate adjoint fields and the true adjoint field for a two-dimensionalconductivity structure. In Chapter 6, approximate sensitivities for the two-dimensionalmagnetotelluric problem are calculated for a simple model and compared with the accurate sensitivities. The approximate sensitivities are then used in the inversion of asynthetic data-set and a field data-set. In Chapter 7, the approximate sensitivities forthe 2.5d problem and the three-dimensional controlled-source problem are calculated andcompared with the accurate sensitivities for simple conductivity models.Chapter 5Approximate Adjoint Fields for theTwo-Dimensional Problem5.1 IntroductionThe approximate form of sensitivity presented in this thesis is based on approximating the adjoint field required in the adjoint-equation method of accurately calculatingthe sensitivities. The rationale for this was given in Chapter 4. Here I investigate various possibilities for the approximate adjoint field in the context of the two-dimensionalproblem. The possibilities are the electric field computed in a homogeneous halfspacefor the appropriate dipole source, the electric field computed in a horizontally-layeredhalfspace, and the Born approximation of the adjoint field. The approximate adjointfields are computed and compared with the corresponding true adjoint field for a simpletwo-dimensional conductivity model.5.2 True adjoint fieldConsider the two-dimensional conductivity model shown in Fig. 5.1 comprising a conductive block of 0.1 S m1 in a more resistive background of 0.01 5 m1, and a conductivebasement of 0.1 S m1. (The mesh used for all calculations relating to this conductivitymodel is shown in Fig. 5.2.) Suppose measurements are made of the electric and magneticfields induced in this model by a source that is also invariant in the strike direction (forexample, the plane-wave source assumed in magnetotelluric experiments). For this purelytwo-dimensional problem, the adjoint field is the electric field due to a two-dimensionaldipole source at the measurement location (see Appendix B). This source is equivalent toan infinite line of point (in three-dimensional space) dipoles which is parallel to the strikedirection and which passes through the measurement location. Suppose a measurement45CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 46—080—0.95I,.-1.09—1.24—1 38—1 52—1.67—1.81—1.95—2.10log10 crFigure 5.1 The two-dimensional conductivity model used here and in subsequent chapters to illustrate the approximate adjoint fields and approximateadjoint sensitivities for the two-dimensional and 2.5d inverse problems.of the along-strike electric field, E, has been made at x0 = —3000 m on the surface ofthe model. From eq. (4.14) and Appendix B, the sensitivity of this measurement withrespect to the model parameter is given by= fEt.Ejds. (5.1),=—3OOO AHere, the adjoint field, Et, is the electric field due to a unit two-dimensional y-directedelectric dipole source (essentially an infinite line current in the along-strike direction) atthe observation location x1, = —3000 m.The true adjoint field for the two-dimensional model in Fig. 5.1 and for the twodimensional p-directed electric dipole source at = —3000 m was computed using theN 100001I15000(m)CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 47UFigure 5.2 The mesh used for allin Fig. 5.1.calculations associated with the model shownfinite-element program of Unsworth, Travis & Chave (1993). The adjoint field was computed at each node in the mesh shown in Fig. 5.2. The amplitude of this adjoint field isshown in Fig. 5.3a. The phase is shown in Fig. 5.4a. The frequency of the source was0.2 Hz. Note the distortion of the phase due to the conductive block, and the increaseddecay of the amplitude and the increase in the phase as the electric field penetrates theconductive basement.5.3 Approximate adjoint field: homogeneous halfspaceThe first possible form of approximate adjoint field that I present is the electric fieldcomputed in a homogeneous halfspace for the appropriate dipole source. The methodused for computing the electric field generated in a homogeneous halfspace by a two-dimensional dipole is briefly described in Appendix C. Fig. 5.3b shows the amplitudeof the electric field computed in a homogeneous halfspace of conductivity 0.028 S m1for a unit two-dimensional y-directed electric dipole source at x0 = —3000 m, and for afrequency of 0.2 Hz. The corresponding phase is shown in Fig. 5.4b. The above value of50001 flflflflN1500020000t’3 — (TI0 (TI 0 00 0 0 00 0 0 00 0 00 (11 I— — (‘30 0 Cl’ 00 0 0 00 0 0 00 0 0x(m)z(m) z(m)Figure 5.3 The amplitudes of the true and approximate adjoint fields for adirected electric line source over the two-dimensional conductivity model shownin Fig. 5.1. The location of the line source is shown by the triangle. Panel (a)shows the true adjoint field, and panels (b) to (d) show respectively the approximate adjoint fields computed in a homogeneous halfspace, in a layered halfspace,and using the Born approximation. The grey-scale represents log10 Ej.CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 48—5.80—6.17—6.53—6.90—7.27—7.63—8.00—8.36—8.73—9.10ct c)C C C C CC C C CC C C CU, C U C— — CUCHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS—10000—15000410.377.344.310.277.243.210.176.143.110.Figure 5.4 The phases of the true and approximate adjoint fields for thetwo-dimensional conductivity model shown in Fig. 5.1. The location of the linesource is shown by the triangle. Panel (a) shows the true adjoint field, andpanels (b) to (d) show respectively the approximate adjoint fields computed ina homogeneous halfspace, in a layered halfspace, and using the Born approximation. The grey-scale represents phase in degrees.49500001500010000—5000—10000—15000—2000020000150001000050000—5000—20000C-)z(m) z(rn)CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 50the conductivity was obtained by averaging the logarithm of the conductivities, weightedby across-strike area, of the model shown in Fig. 5.1.The approximate adjoint field shown in Figs. 5.3b and 5.4b obviously contains noinformation about the conductive block nor the conductive basement. However, the dominant behaviour of the true adjoint field is the decay in amplitude and increase in phaseaway from the source, which is exactly the behaviour of the approximate adjoint fieldcomputed in the homogeneous halfspace. The effect of the two-dimensional structure inthe conductivity model shown in Fig. 5.1 can almost be relegated to that of a perturbationto the field that would otherwise be present in the homogeneous halfspace. It is felt thatthis level of agreement between true and approximate adjoint fields should be sufficientfor the subsequent approximate sensitivities to be useful in an iterative inversion scheme.5.4 Approximate adjoint field: layered halfspaceThe second form of approximate adjoint field is that computed in a layered halfspace.An appropriate layered halfspace can be constructed by averaging the conductivities ofthe multi-dimensional model in each of a series of horizontal layers. The electric field generated in a layered model by a two-dimensional source can be computed using the methoddescribed in Appendix A. Doing this for a unit two-dimensional y-directed electric dipoleat x = —3000 m and for a frequency of 0.2 Hz gives the electric field shown in Figs. 5.3cand 5.4c. The layered halfspace used for this example was obtained by averaging thelogarithms of the conductivities, weighted by across-strike area, in each horizontal layerof the model in Fig. 5.1. Just as for the approximate adjoint field computed in a homogeneous halfspace (see Figs. 5.3b and 5.4b), there is, of course, no manifestation in thefield shown in Figs. 5.3c and 5.4c of the two-dimensional effects of the conductive block.However, the rapid decrease in amplitude and increase in phase of the true adjoint fieldin the conductive basement is well reproduced in the approximate adjoint field computedin the layered halfspace. This is not surprising since this feature of the two-dimensionalCHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 51model is replicated in the layered halfspace used for computing the approximate adjointfield.5.5 Approximate adjoint field: Born approximationThe third form of approximate adjoint field presented here results from the Bornapproximation. Consider the integral equation describing the true adjoint field, E, inthe multi-dimensional conductivity model (eq. 27, Hohmann 1988):Et(r) = E(r) + j G(r; r’) Et(r’) (r’) dv’, (5.2)where is the electric field calculated in some background conductivity model and Lio isthe difference between the background and true conductivity structures. G is the tensorGreen’s function for the background conductivity model. Eq. (5.2) is exact. The Bornapproximation for E is obtained by replacing E in the integrand by E:E(r) E(r) + j G(r; r’) . E(r’) (r’) dv’. (5.3)The elements of the tensor Green’s function,/GG(r;r’)= G?,,z ) , (5.4)\ G Jare such that is the x-component of the electric field generated in the background conductivity model by a unit y-directed electric dipole. For the E-polarisation case, eq. (5.3)reduces to a scalar equation and the only relevent component of the tensor Green’s function, for a homogeneous halfspace background model is given by Hohmann (1988).For the il-polarisation case, the relevent components of the tensor Green’s function (G,G & G) for a homogeneous halfspace are given by Lee & Morrison (1984).CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 52The Born approximation for the adjoint field was calculated for the conductivitymodel in Fig. 5.1. The background conductivity model in which and were caldulated was a homogeneous halfspace of 0.01 S m1. This value is the same as that for theresistive background in the model in Fig. 5.1, and was chosen to minimise the number ofcells over which the numerical integration in eq. (5.3) had to be carried out. The resultingapproximate adjoint field generated by a two-dimensional y-directed electric dipole sourceis shown in Figs. 5.3d and 5.4d. As before, a frequency of 0.2 Hz was used. From Fig. 5.3it can been seen that the Born approximation produces a poorer approximation of theamplitude of the adjoint field than that computed in the homogeneous halfspace. Thisis because of the difference between the background conductivity used to compute theBorn approximation and the conductivity of the homogeneous halfspace used to produceFig. 5.3b. However, the phase of the Born approximation (Fig. 5.4d) is somewhat betterthan that computed for the homogeneous halfspace above about 1000 m and does givean indication of the increase in phase that the true adjoint field exhibits in the conductive basement, although this is not nearly as well reproduced as in the approximationcalculated in the layered model.5.6 Computation timesThe true and approximate adjoint fields described in this chapter were computedon a Sun Sparclo workstation. The time required to compute each of these is listedin Table 5.1. For the given values, the adjoint field was computed at the 861 nodes inthe mesh shown in Fig. 5.2 (21 vertically and 41 horizontally) for the two-dimensionaly-directed electric dipole source. The approximate adjoint fields calculated in a homogeneous or layered halfspace have similar computation times which are significantly lessthan the computation time for the true adjoint field. The Born approximation takesconsiderably longer to compute than the true adjoint field. This is because the Green’sfunction needed in eq. (5.3), that is the electric field in the homogeneous halfspace dueCHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 53True adjoint field 75Approximate adjoint field: homogeneous halfspace 1layered halfspace 2Born approximation 25Table 5.1 Example computation times for the various types of two-dimensionaladjoint fields.to a two-dimensional electric dipole, has to be computed for dipoles at depths within themodel, not just for one at the surface.5.7 ConclusionsIn this chapter, I have presented three forms of approximate adjoint field. An example for each was computed for a simple two-dimensional conductivity model and comparedwith the true adjoint field. The approximate adjoint field computed in the homogeneoushalfspace replicates the dominant features of the true adjoint field: the decay in amplitude and increase in phase away from the dipole source. The approximate adjoint fieldcomputed using the layered halfspace matches the behaviour of the true adjoint field inthe layered features within the two-dimensional model. The Born approximation of theadjoint field exhibits some of the two-dimensional features present in the true adjoint field,although it does not match the amplitudes of these features. These points are furtherillustrated by Figs. 5.5 and 5.6: the differences between the logarithm of the amplitudeof the true adjoint field and the logarithm of the amplitude of the three approximateadjoint fields are shown in Fig. 5.5, and the differences between the phases are shown inFig. 5.6. It is further evident from these figures that all three approximate adjoint fieldsare good overall approximations of the true adjoint field.The approximate adjoint field computed in the homogeneous and layered halfspacesare significantly quicker to obtain than the true adjoint field. The Born approximation,however, is substantially slower to compute than the true adjoint field.CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 54In the next chapter, I shall use the three forms of the approximate adjoint fieldintroduced in this chapter to compute approximate sensitivities for the two-dimensionalmagnetotelluric problem.Figure 5.5 The differences between the logarithm of the amplitude of thetrue adjoint field shown in Fig. 5.3a and the logarithm of the amplitudes of theapproximate adjoint fields shown in Figs. 5.3b to 5.3d. The panels show thedifferences between the true adjoint field and (a) the approximate adjoint fieldcomputed in a homogeneous halfspace, (b) that computed in a layered halfspace,and (c) using the Born approximation.CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 550 1700.049—0.072—0.193—0.314—0.436—0.557—0.678—0.799—0.920——20000z(m)Sz(m)CHAPTER 5: 2d APPROXIMATE ADJOINT FIELDS 56Figure 5.6 The differences between the phase of the true adjoint field shownin Fig. 5.4a and the phase of the approximate adjoint fields shown in Figs. 5.4bto 5.4d. The panels show the differences between the true adjoint field and(a) the approximate adjoint field computed in a homogeneous halfspace, (b) thatcomputed in a layered halfspace, and (c) using the Born approximation.C.)z (in)90.073.356.640.023.36.7—10.0—26.7—43.3—60.0Sz(m)Chapter 6Approximate Sensitivities for the Two-DimensionalMagnetotelluric Inverse Problem6.1 IntroductionIn the previous chapter, I compared three forms of approximate adjoint field thatcould be used for a two-dimensional inverse problem. In this chapter, I use these adjointfields to calculate approximate sensitivities for the two-dimensional magnetotelluric problem. These sensitivities are compared with sensitivities calculated using the brute-force(or perturbation) method. In the final two sections of this chapter, the approximatesensitivities are used within a minimum-structure, Gauss-Newton algorithm to invertsynthetic and field data.6.2 Calculation of the sensitivities6.2.1 Sensitivities for apparent resistivity and phaseThe observed quantities in a magnetotelluric experiment are usually considered tobe the frequency-domain values of apparent resistivity, pa, and phase, ç. For a twodimensional Earth, for which the magnetotelluric problem can be split into two decoupledmodes, the apparent resistivity and phase are given bypaQ”) = and b(w) = phase () (6.1)for the E-polarisation mode, and bypa(W) = and (w) = phase () (6.2)57CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 58for the il-polarisation mode, where w is the angular frequency and t is the magneticpermeability. The strike direction is assumed to lie in the y-direction. The sensitivitiesfor the apparent resistivity and phase can be expressed in terms of the sensitivities forthe electric and magnetic fields (see Appendix D):= 2Pa {e () - e iL) } (6.3)and(1OE”\ (1OH’\= sm- (6.4)These relationships therefore enable approximate sensitivities for the apparent resistivityand phase to be calculated from the approximate sensitivities for the electric and magneticfields obtained using the method described in Chapter 4.6.2.2 Area integrationFor the two-dimensional magnetotelluric inverse problem, the sensitivities for theelectric and magnetic fields have the form (see section 5.2)=jEt.Ezi&jds (6.5)where F represents one component of the electric or magnetic field (either E, E, H orHr). For the typical situation in which the model comprises rectangular cells of uniformconductivity, that is the basis function is equal to 1 within the j” cell and equal tozero everywhere else, the area integration in eq. (6.5) becomesOF = J J E . E dz dx, (6.6)9o-j 0 ZOCHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 59Z (x0,z) (x1,z0)(a0,z1) (xj,z1)Figure 6.1 The notation used for the bilinear interpolation and integrationof the electric and adjoint fields over a cell in a two-dimensional conductivitymodel.where, for simplicity of notation, the th cell is assumed to extend from z0 to x andfrom z0 to z1 as shown in Fig. 6.1.The forward-modelled electric field E will usually only be known at the nodes in themodel (that is, at the vertices of the cell shown in Fig. 6.1). To calculate the integral ineq. (6.6), the real and imaginary parts of the components of E are bilinearly interpolatedover the rectangular cell:E(x,z) = c1xz + c2z + c3x + c4. (6.7)E represents the real or imaginary part of the component of the electric field. For ageneral cell, the 4 x 4 system of equations relating the values of E at the vertices of thecell to the coefficients c, i = 1,... , 4, in the above expansion was solved analytically. Theresulting expressions are then used to calulate the coefficients for a particular cell fromthe values of E at the vertices of that cell. The adjoint field, E is bilinearly interpolatedin an analogous manner. With the real and imaginary parts of the components of Eand E in the form of eq. (6.7), the integration with respect to x and z in eq. (6.6) couldalso be carried out analytically for a general cell. The resulting expression is a functionCHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 60of the coefficients in the bilinear expansions of E and E. For a particular cell, therefore,the integral in eq. (6.6) can be evaluated, using analytic formulae, from the values of theelectric and adjoint fields at the vertices of the cell.It was found that the horizontal discretisation of the mesh (see Fig. 5.2) used tocalculate the sensitivities in the next two sections was too coarse for the above techniqueto accurately evaluate the integral over the cells close to the dipole, especially for theH-polarisation mode. The adjoint field for this mode varies more rapidly near the sourcethan the adjoint field for the E-polarisation mode shown in Figs. 5.3 and 5.4. The cells inthe mesh were therefore sub-divided horizontally and the forward-modelled electric field(which does not vary rapidly horizontally) linearly interpolated over these sub-divisions.The adjoint field was explicitly computed at the vertices of the sub-divisions. This required an insignificant number of additional computations. The technique described inthe previous paragraph was then used to evaluate the area integration over each subdivision.6.3 Comparison of approximate and accurate sensitivities6.3.1 Homogeneous halfspacePrograms were written to compute approximate sensitivities for the two-dimensionalmagnetotelluric problem using each of the three forms of approximate adjoint field described in Chapter 5. The programs were first tested by computing the approximatesensitivities for a homogeneous halfspace. For this special case, the approximate adjointfields are not in fact approximate but exact, and so the resulting sensitivities are exact.Consider a homogeneous halfspace of conductivity 0.01 S m1 with the same horizontal and vertical dimensions as the model shown in Fig. 5.1. Suppose measurements ofthe E-polarisation apparent resistivity and phase had been made at x, = —3000m on thesurface of the model at a frequency of 0.2 Hz. Sensitivities were calculated for these observations using the brute-force (or perturbation) method. The model was divided into the800 cells (20 vertically and 40 horizontally) shown in Fig. 5.2. The conductivity of eachCHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 61cell was perturbed in turn by 10% (this perturbation is large, but was required to obtainaccurate values for the cells far from the observation location) and a forward modellingperformed to calculate the resulting changes in the apparent resistivity and phase. Thesensitivity with respect to the perturbed model parameter is then given by the changein the apparent resistivity or phase divided by the change in the model parameter. (Theforward modelling was carried out using the transmission-surface modelling program ofMadden 1972.) These sensitivities are shown in Figs. 6.2a and 6.2b. The sensitivity isplotted as a constant value over the area of the cell to which it refers.Approximate sensitivities were calculated for the homogeneous conductivity modeldescribed above. The forward-modelled electric field was computed using the program ofMadden (1972). All three forms of approximate adjoint field were used, and the resultingsensitivities compared with the brute-force sensitivities shown in Figs. 6.2a and 6.2b.The E-polarisation sensitivities calculated using the adjoint field computed in a layeredhalfspace (all the layers of which had the same conductivity) are shown in Figs. 6.2cand 6.2d. There is very good agreement between the brute-force sensititivities and thesensitivities calculated using the adjoint field in the layered halfspace. The approximatesensitivities calculated using the adjoint field computed in a homogeneous halfspace andthose calculated using the Born approximation of the adjoint were indistinguishable fromthose in Figs. 6.2c and 6.2d.A comparison was also made of the brute-force and approximate sensitivities forthe H-polarisation mode. The brute-force sensitivities are shown in Figs. 6.3a and 6.3b.The approximate sensitivities calculated using the adjoint field in a layered halfspace areFigure 6.2 (Following page) Sensitivities for the E-polarisation apparent resistivity and phase for a homogeneous halfspace of conductivity 0.01 S m1 andfor a frequency of 0.2 Hz. The observation location is indicated by the triangle. Panels (a) and (c), and the colour-bar on the left of the figure, showlog10ãlnp/Olno. Panels (b) and (d), and the colour-bar on the right, showlog8q5/6lncr. Panels (a) and (b) were produced by the brute-force method,and panels (c) and (d) were produced by the approximate-sensitivity programusing an adjoint field computed in a layered halfspace.CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 62o L) 0 0C C o— — —— CVI I ICV N — tf) CCD 0 -CV CV CV C) 0I I I I I20000150001000050000—5000—10000—15000—20000200001600010000—5000—10000—15000—20000I 1; 1o N C N C N CC’ I C’- 0 CV N 0-jCV CV CV C’) C’) C’) C’) C’) ‘ ‘0—5000—10000—15000—200000 L’0 0o 0o oz(m)CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 63shown in Figs. 6.3c and 6.3d. Again there is good agreement between the brute-force andapproximate sensitivities.6.3.2 Two-dimensional modelApproximate sensitivities were calculated for the two-dimensional model shown inFig. 5.1. As in the section above, measurements of the apparent resistivity and phasemade at Xc, = —3000 m on the surface of the model at a frequency of 0.2 Hz wereconsidered. The E-polarisation brute-force sensitivities are shown in Figs. 6.4a and 6.4b,and the approximate sensitivities in Figs. 6.4c and 6.4d. The approximate sensitivitieswere calculated using the approximate adjoint field computed in a layered halfspace. (Forthis example, the approximate adjoint field is that shown in Figs. 5.3c and 5.4c.)There is good agreement between the brute-force and approximate sensitivities shownin Fig. 6.4. A comparable, although slightly poorer, level of agreement was found for thetwo other forms of approximate sensitivity. It should be noted that, because the sensitivities shown in Fig. 6.4 are with respect to ln a, they depend directly on the conductivitymodel as well as the behaviour of the electric and adjoint fields within the model. (ThisFigure 6.3 (Following page) The H-polarisation sensitivities for the homogeneous halfspace. Panels (a) and (c), and the colour-bar on the left of the figure,show log10jOlnp/Oln4 Panels (b) and (d), and the colour-bar on the right,show logOq/Olno. Panels (a) and (b) were produced by the brute-forcemethod, and panels (c) and (d) were produced by the approximate-sensitivityprogram using an adjoint field computed in a layered halfspace.Figure 6.4 (Page 65) Brute-force and approximate sensitivities for the Epolarisation apparent resistivity and phase for the conductivity model shownin Fig. 5.1. The observation location is indicated by the triangle. The frequency was 0.2 Hz. Panels (a) and (c), and the colour-bar on the left of thefigure, show log10Olnp/O1ncrj, and panels (b) and (d), and the colour-baron the right, show log10Oq/Olnoj. Panels (a) and (b) show the sensitivitiesproduced by the brute-force method, and panels (c) and (d) show the approximate sensitivities calculated using the approximate adjoint field computed in alayered halfspace.CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 640 If) C) — C) f) 0C D C C- — If) C) C’)0 - — CIf CIf CIf C’) C’) C’)I I I I I II’0 00 0o 0If) C— C’)20000150001000050000—5000—10000—15000—20000-dC)z(m) z(m)ii0 C’- 0 C’- 0 C’- CC’) CC C) C’) CC C) C’) CD C) C’),_i -— C C’) C’) C’) C’) C’)I I I I I I ICHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MTI IiC LID OD ID 0cc Q C’) ID C C’) IL)C o———— C’2 CQ CIL C’I I I I I I I I I11o II) 0) ID C’) C’- — 11) 0C C’) C C’- — ‘— CQ CIL CL) C’) C’) C’) ‘ ‘ ‘I I I I I I I I6500IDz(m) z(m)CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 66can be seen from the relation O/O(ln o) = o 6/0o.) The dramatically increased sensitivity in the conductive block that results from this, and the enhanced sensitivity in theconductive basement, can be clearly seen in Fig 6.4. This is the two-dimensional manifestation of what Boerner & Holladay (1990) observed in their “fractional” sensitivitiesfor the one-dimensional problem (see section 3.3.1).The H-polarisation sensitivities were also calculated for the two-dimensional conductivity model. The brute-force sensitivities are shown in Figs. 6.5a and 6.5b, andthe approximate sensitivities calculated using the adjoint field in a layered halfspace areshown in Figs. 6.5c and 6.5d. The agreement between the brute-force and approximatesensitivities is similar to that for the E-polarisation sensitivities. Again, the two otherforms of approximate sensitivity gave a slightly poorer level of agreement.6.4 Computation timesTable 6.1 lists the computation times for the various methods of calculating theJacobian matrix of sensitivities for a small two-dimensional magnetotelluric example.The sensitivities were calculated for 80 data (observations of apparent resistivity andphase for both E- and H-polarisations at four locations and at five frequencies) withrespect to the conductivities of the 800 cells used to parameterise the model in Fig. 5.1.The approximate adjoint fields in both the homogeneous and layered halfspaces wereonly explicitly computed for one observation location and then translated horizontallyto calculate the sensitivities for the three remaining observation locations. This resultsin a considerable time-saving that becomes more significant with increasing number ofobservation locations. In a similar manner, for the Born approximation of the adjointFigure 6.5 (Following page) Brute-force and approximate sensitivities for theH-polarisation mode for the conductivity model shown in Fig. 5.1. Panels (a)and (c), and the colour-bar on the left of the figure, show log10Olnp/ôlncr,and panels (b) and (d), and the colour-bar on the right, show log10 q/Olnu.Panels (a) and (b) show the sensitivities produced by the brute-force method,and panels (c) and (d) show the approximate sensitivities.CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MTo a a a a a a a a ac— q c’ COa a a — — - C\ C’ 01 01I I I I I II ItJ *1 I 1o 0) 0 CO U) C) 01 — 00) 0) 0 C’) 0) C—— CU CU C’) C LI)I I I I I I I I67Iaa0U)z(m)aaa aa U)z(m)CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 68field, both and G, which were calculated in a homogeneous halfspace, were explicitlycomputed for only the first observation location and then translated horizontally forthe other observation locations. This reduces the computation time for this third form ofapproximate sensitivities, although it is still not less than that of the accurate sensitivities.It is not until there are tens of observation locations that the approximate sensitivitiesbased on the Born approximation of the adjoint field require less time to calculate thanthe accurate sensitivities.Accurate sensitivities: brute-force 500 mmadjoint-equation 2.5Approximate sensitivities: homogeneous halfspace 0.25layered halfspace 0.5Born approximation 5Table 6.1 Example computation times for the various types of sensitivities forthe two-dimensional magnetotelluric inverse problem.6.5 Two-dimensional inversion of magnetotelluric dataSo far in this chapter, I have been concerned with the comparison of the values ofthe approximate and accurate sensitivities. Although this is important, the usefulness ofthe approximate sensitivities ultimately depends on their successful performance withinan inversion routine. In this section, therefore, I test the behaviour of an iterative,linearised inversion procedure that uses the approximate sensitivities. The inversionstrategy discussed in Chapter 2 is used to invert synthetic and field magnetotelluric datato recover two-dimensional conductivity models.The sensitivities are just one of several components in an inversion procedure, and sothe effectiveness of the approximate sensitivities will depend on the details of the inversionstrategy used: convergence of one algorithm using approximate sensitivities does notmean that all algorithms will converge. However, testing approximate sensitivities in allCHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 69possible algorithms is clearly impractical. Here I consider an algorithm that is typicalof those that one would like to apply to the multi-dimensional electromagnetic inverseproblem.6.5.1 Synthetic dataApproximate sensitivities, calculated using the approximate adjoint field in a layeredhalfspace, were used to invert a set of synthetic data generated from the conductivitymodel in Fig. 5.1. Values of the E- and H-polarisation apparent resistivity and phasewere calculated at four observation locations (Xe —11000, —3000,5000 and 13000m)for five frequencies (1, 0.5, 0.2, 0.1 and 0.05 Hz). Gaussian random noise was added tothese data. The standard deviation of the noise added to the apparent resistivity was 5%,and the standard deviation of the noise added to the phase was 2 degrees. The data areshown in Fig. 6.6.The 80 synthetic data were inverted using an iterative, linearised, minimum-structureinversion procedure in which the system of equations at each iteration was solved usinga subspace technique. The details of the inversion algorithm are given in Oldenburg &Ellis (1993). The model norm that was minimised was the discrete equivalent of= fa3w(m — rn0)2 +JA 1/O(m_mo)2______cw Ox ) + wz ) ) ds, (6.7)where the model, m, was equal to ln o, and the reference model, m0, corresponded toa halfspace of conductivity 0.01 5 m1. The values of the coefficients were = 10—2,= 3 and c = 1. The weighting functions, w8, w and w were constant and equal tounity throughout the model, except for w, which was equal to 10 in the top four layersof the model. This was to suppress structure developing in the near-surface layers as aCHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MTFigure 6.6 The synthetic data (shown by the circles and error bars) generatedfrom the model shown in Fig. 5.1, and the predicted data (shown by the lines)from the final model produced by the inversion of the synthetic data. Thefilled circles and solid lines correspond to the E-polarisation data, and the opencircles and dashed lines to the H-polarisation data. The units of the apparentresistivity, Pa are fm and those of the phase, q, are degrees.70701021021 2102102Ct’qa: ,aE‘. -0.05 Hz—15000 0 15000j Tz—15000 0 15000rZ—15000 0 150000.5Hz—15000 0 150001Hz—I- 0.05 HzE—15000 0 15000700.1 Hz60 -J—15000 0 15000700.2 Hz60-5- f---;_;40 p—15000 0 1500070- 0.5 lIz—15000 0 1500070- 1 Hz:x(m)40—15000 0 15000 —15000 0 15000x(m)CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 71result of the large H-polarisation sensitivities for these layers (see Figs. 6.3 and 6.5). Thesubspace vectors were the steepest-descent vectorsvk = (Wm)’ Vm (6.8)where YLm is the weighting matrix associated with the discrete equivalent of eq. (6.7),and is the data misfit for the kt partition of the data. The data were orderedaccording to their individual misfit and then divided into, in this case, 20 partitions. Theapproximate sensitivities were calculated using the approximate adjoint field in a layeredhalfspace. The conductivities of the layers were chosen by averaging the logarithms ofthe conductivities within each layer of the model.The final model is shown in Fig. 6.7. The conductive block has approximatelythe correct amplitude and is situated at the correct depth and horizontal location. Theconductive basement has also been recovered. The predicted data produced by this modelare shown in Fig. 6.6.The values of the data misfit, model norm and Lagrange multiplier during the inversion are shown in Fig. 6.8. The data misfit decreases at the requested rate (/3 = 0.8) forthe first seventeen iterations. The inversion algorithm can then only decrease the misfitmore slowly than the requested rate. The target misfit of 80 is achieved after 31 iterationsand the misfit then remains constant for the remaining iterations. The model norm gradually increases before levelling off once the algorithm has achieved the target data misfit.The Lagrange multiplier steadily increases for the first twenty iterations, but then variesconsiderably for the next ten. This coincides with the iterations for which the algorithmcannot decrease the misfit at the requested rate. It appears as if the algorithm is havingto work hard to lower the misfit to the target value. Once the target misfit has finallybeen achieved, the Lagrange multiplier levels off and shows only small variations as thealgorithm attempts to further minimise the model norm.The results of the above inversion of synthetic data are what would be expected ifaccurate sensitivities had been used in place of the approximate sensitivities: the finalCHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT—0.80—095—1.09—1.24—1.38—1.52—1.67—1.81--1.95—2.10log10Figure 6.7 The final model produced by the inversion of the synthetic datashown in Fig. 6.6. The model used to generate these data is the one shown inFig. 5.1. The colour-bar is the same for both figures. The triangles indicate theobservation locations.72model reproduces the data to the desired level of misfit and has the typical smeared-out appearance of a model produced by a minimum-structure inversion algorithm. Thebehaviour of the data misfit, model norm and Lagrange multiplier (with the possibleexception of the oscillations between iterations 24 and 30) is also what would have beenexpected if accurate sensitivities had been used in the inversion procedure.As a final comment for this section, if the approximate sensitivities of Smith &Booker (1991) had been used to invert the synthetic magnetotelluric data considered inthis section, the model could only have contained four columns of cells, one beneath eachobservation location. The horizontal resolution of the final model would therefore havebeen considerably less than that in Fig. 6.7.N 10000IS— 15000x(m)CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 73a102 -0 10 20 30 40i o -io101 b1000 10 20 30 401 o310000 10 20 30 40L er a Li onFigure 6.8 The data misfit, q, model norm, 1m’ and Lagrange multiplier, ,as functions of iteration during the inversion of synthetic data discussed insection 6.5.1.6.5.2 COPROD2 dataTo conclude this chapter, the approximate sensitivities were used in the inversion ofa set of field data. The data-set chosen was a sub-set of the COPROD2 data-set collectedin southern Saskatchewan and Manitoba. A description of the COPROD2 data-set andtheir context is given by Jones (1993). The data are dominated by the effects of the NorthAmerican Central Plains (“NACP”) anomaly: a linear zone of high conductivity at midcrustal depths that runs north from South Dakota through Saskatchewan before turningeastwards under Hudson Bay. The anomaly appears in the COPROD2 data as low valuesof the E-polarisation apparent resistivity and high values of the E-polarisation phase. Incontrast, the H-polarisation measurements show little indication of the anomaly.CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 74The COPROD2 data have been used by the magnetotelluric community as a testdata-set for the comparison of two-dimensional inversion schemes (see Ingham, Jones &Honkura 1993). The considerable splitting of the E- and H-polarisation modes mentionedabove makes the COPROD2 data-set a challenging test for any inversion procedure, especially considering the size of the problem is close to the limit of what can presentlybe contemplated. All the inversion programs thus far applied to the COPROD2 data-sethave struggled to some extent, either because of the nature of the data or because oftheir considerable number.The sub-set of the COPROD2 data chosen for consideration here is the same asthat used by Ellis, Farquharson & Oldenburg (1993), except that stations to the east of118 km, and so relating to the Thompson belt (“TOBE”) anomaly, were omitted. Thetotal number of data was therefore 896. The error estimates were the same as those ofEllis et al. The observed values of apparent resistivity and phase are shown in Fig. 6.9.For the inversion carried out here, the model was divided into 30 layers and 62columns as shown in Fig. 6.10. The columns were arranged so there was one columnbeneath each observation location and an additional column between each pair of stations.There were therefore just over twice as many columns as there were stations. The resultingnumber of cells was 1860.The inversion algorithm was the same as that used in section 6.5.1. Approximatesensitivities were calculated using the approximate adjoint field computed in a layeredhalfspace. The conductivities of the layers of this halfspace were obtained by averagingthe logarithms of the conductivities in each layer of the two-dimensional model. A one-dimensional inversion of the complete data-set was carried out to produce the best-fittinghorizontally-layered model. This model, which has a conducting surface layer of 0.5 S m1that ramps down to a resistive basement of 0.01 S m1 below about 7 km, was then usedas the reference model in the two-dimensional inversion. This reference model was notincluded in the horizontal and vertical gradient terms (see eq. 6.7) in the weightingmatrix Lm The subspace vectors were again obtained by partitioning the data in termsof misfit, but 200 partitions were used here rather than the 20 in section 6.5.1.CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 751020.023 Hz 60 0.023 HzI I______—200 —100 0 100 —200 —100 0 1001020.047 Hz 60 0.047 HzI__I__I__I__I_ ____ _ _—200 —100 0 100 —200 —100 0 1001020.093 Hz 60 0.093 Hz10120100 0 I—200 —100 0 100 —200 -100 0 1001 20.19 Hz 60101 °100_____________20- 0.19 liz—200 —100 0 100 —200 —100 0 100x (km) x (km)Figure 6.9 The sub-set of the COPROD2 data considered in section 6.5.2.The observations are represented by the circles, crosses and error bars. Thelines indicate the predicted data for the final model produced by the inversion(see Fig. 6.11). The filled circles and solid lines correspond to the E-polarisationdata, and the crosses and dashed lines to the il-polarisation data. The units ofthe apparent resistivity, Pa are m and those of the phase, q, are degrees.CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 761 2100102604020—200 —100 0 100604020102CO 10100—2001 2COLU100—100—200 —1000 1000 100x (km)604020- 0.0059 Hz0—20060‘10 -200100—200 —1000 1000 100x (km)CO hi—200 —100 0 10000.0029 lIziOo I I—200 —100 0 1000.0029 1-Tz0—200 —1 00 0 1000.0059 HzI I I I0.012 Hz- 0.012 HzFigure 6.9 (contd.)CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 77010______2030-‘-‘ 405060I I I 0 ——— (31 0 0o a a 0o a ax (km)Figure 6.10 The mesh for the model used in the inversion of the COPROD2data. The final model is shown in Fig. 6.11. The vertical exaggeration is 2.9:1(this is different from Fig. 6.11).As mentioned at the start of this section, the COPROD2 data-set has proved to bechallenging to many inversion programs. This was also the case here. The main difficultywas that both the E- and H-polarisation sensitivities were large for the uppermost layers inthe model. The values of w3 and w for these layers were therefore increased to counteractthe development of excessive structure that these large sensitivities would otherwise cause.The process of choosing the most effective form of this additional weighting was empirical.I therefore ran the inversion program numerous times, each using a different form ofweighting in the near-surface layers. Once the misfit had stopped decreasing, I wouldalter the weighting to see if a further decrease in the misfit could be achieved. The finalresult of the inversion process was the model shown in Fig. 6.11. The predicted data forthis model are shown in Fig. 6.9, and the corresponding value of the x2 misfit is 1370.This is greater than the target value of 896. However, the final misfit is close to thistarget misfit, and the comparison between the observed and predicted data in Fig. 6.9 isgood.CHAPTER 6: APPROXIMATE SENSITIVITIES FOR 2d MT 78In summary, therefore, the inversion of the CQPROD2 data using the approximatesensitivities has produced a model that reproduces the observations to almost the desiredlevel of misfit. The final model also contains the multiple, distinct, mid-crustal conductorsthat have been produced by other inversions of the COPROD2 data.6.6 ConclusionsIn this chapter, I have calculated approximate sensitivities for the two-dimensionalmagnetotelluric inverse problem. These sensitivities were compared with brute-force sensitivities for a two-dimensional conductivity model. It was shown that the agreement between the approximate and brute-force sensitivities is good, especially when consideringthe sensitivity with respect to the logarithm of the conductivity. The level of agreement,as illustrated in Figs. 6.6 and 6.7, seems sufficient to enable iterative, linearised inversions of magnetotelluric data which use these approximate sensitivities to converge to anacceptable model. This was emphasised by successful inversions of a synthetic data-setand of the COPROD2 data.It was also shown that the approximate sensitivities are considerably quicker tocompute than even the most efficient method of accurately calculating the sensitivities.For the small example considered in section 6.4, a factor of five difference was obtained.This factor will increase with the size of the problem, since the approximate adjointfield is calculated for the dipole at only one observation location and then translated forall other observation locations. Calculating the accurate sensitivities using the adjointequation method requires the adjoint field to be computed for a dipole at each observationlocation. The difference-factor for the inversion of the COPROD2 data in section 6.5.2was estimated to be two orders of magnitude.Figure 6.11 (Following page) The final model produced by the inversion, usingthe approximate sensitivities, of the sub-set (shown in Fig. 6.9) of the COPROD2data. The vertical exaggeration is approximately 2.3:1.() I0.900.35—0.19—0.74—1.28—1.82—2.37—2.91—3.45—4.00c(km)log10aChapter 7Approximate Sensitivities for the 2.5d andThree-Dimensional Inverse Problems.1 IntroductionIn Chapter 6, I investigated approximate sensitivities in the context of the two-dimensional magnetotelluric problem. I presented a comparison of the approximate andaccurate sensitivities for a simple two-dimensional conductivity model. This comparisoncovered all the observed quantities relevent to the two-dimensional magnetotelluric problem: measurements of apparent resistivity and phase for both the E- and il-polarisations.It was shown that there was good agreement between the accurate and approximate sensitivities. It was also shown that the computation time for the approximate sensitivitieswas considerably less than for any of the accurate methods of calculating the sensitivities. Finally, in the last part of Chapter 6, I presented the results of two inversionsusing the approximate sensitivities, one of synthetic data and one of field data (a sub-setof the COPROD2 data). The approximate sensitivities performed successfully in bothinversions.In this chapter, I consider the approximate sensitivities for the 2.5d problem andthe three-dimensional controlled-source problem. The main goal is to construct computational algorithms to calculate the approximate sensitivities for these more realisticand widespread problems, and to determine how much faster the approximate sensitivities are to compute than the accurate sensitivities. A complete investigation similar tothat in Chapter 6, including test inversions using the approximate sensitivities, was notcontemplated. It was assumed that if a similar level of agreement to that in Chapter 6was obtained between the approximate and accurate sensitivities, then the approximate80CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 81sensitivities would perform as well in a 2.5d or three-dimensional inversion as they did inthe two-dimensional inversions of magnetotelluric data.In the first half of this chapter, I consider the approximate sensitivities for the 2.5dinverse problem in which the electric and magnetic fields vary in all three dimensionsbut the conductivity model varies in only two-dimensions (the vertical direction andone horizontal direction). The approximate sensitivities are compared with accuratesensitivities for a homogeneous halfspace and for a simple conductivity model, and thedifference in computation time between the approximate and accurate sensitivities isdemonstrated.In the second half of this chapter, I consider the approximate sensitivities for thethree-dimensional controlled-source problem. This is the most general geophysical electromagnetic inverse problem and is the one we ultimately want to solve on a routinebasis. However, the extent to which a comparison of approximate and accurate sensitivities could be made was severely restricted by the limitations placed on three-dimensionalforward modelling by current computer technology. The only feasible comparison wastherefore for a model in which a confined region had been parameterised into 5 x 5 x 5cells. Using this parameterisation, the approximate and accurate sensitivities were compared for a homogeneous model and for a simple three-dimensional conductivity structure.These comparisons also served to illustrate the considerable difference in computationtimes for the approximate and accurate sensitivities.1.2 Calculation of the sensitivities for the 2.5d problemTo calculate the sensitivities for the 2.5d problem, the general expression given byeq. (4.14) is first manipulated to take advantage of the particular nature of the problem(Unsworth & Oldenburg, 1995). As for the two-dimensional magnetotelluric problemdiscussed in Chapter 6, the model for the 2.5d problem is usually parameterised in termsof rectangular cells of uniform conductivity. These cells are infinite in extent in the strikeCHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 82direction (assumed as before to lie in the y-direction). The expression for the sensitivitycan therefore be written as= jfEt.Edsdy, (7.1)oo A3where F represents the appropriate component of the electric or magnetic field (Es, E,E, H, H or H) and A is the area of the cell (in the across-strike plane).The forward modeffing for the 2.5d problem was carried out using the program ofUnsworth, Travis & Chave (1993). This program essentially solves a two-dimensionalproblem for each in a sequence of along-strike wavenumbers k. The required electricor magnetic field is then obtained by performing the one-dimensional inverse Fouriertransform in k. The forward-modelled electric field, E, in eq. (7.1) can therefore bewritten asE(x,y,z) = _j E(x,k,z) dk, (7.2)where y8 is the y-coordinate of the source. E(x, k, z) emerges naturally from the programof Unsworth et al. The approximate adjoint field, E, can be expressed in a similarmanner:Et(x,y,z) =— j Et(x,1c,z) eio) dk, (7.3)where y0 is the y-coordinate of the observation location. Et(c, ks,, z) is easily obtainedwhen it comes to considering an approximate adjoint field. Substituting the above twoexpressions for the forward-modelled and adjoint fields into eq. (7.1) gives= j-- J J f A(k,k) ei eio) dk dk dy, (7.4)3 —00 —00 —00where the area integration, A, in the along-strike-wavenumber domain is defined asA(k,k) = f Et(x,k,z) .E(,k,z) ds. (7.5)CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 83Using the Fourier transform definition of the Dirac delta function:S(k) =—j:e’ dy, (7.6)the integration with respect to y in eq. (7.4) can be performed to give= (7.7)assuming, for simplicity, that y = y0 = 0. This equation reduces toOF 1[= (7.8)To calculate the sensitivities for the 2.5d problem, therefore, the area integration in thealong-strike-wavenumber domain (eq. 7.5) is evaluated at each of the wavenumbers required by the forward modeffing. The area integration at each wavenumber is carriedout using the same method as for the purely two-dimensional problem described in section 6.2.2. The sensitivity for a particular cell is then obtained by the integration overthe wavenumbers in eq. (7.8).7.3 Comparison of approximate and accurate 2.5d sensitivities7.3.1 Homogeneous halfspaceAs for the two-dimensional magnetotelluric problem, the program for computing theapproximate sensitivities was first tested on a homogeneous halfspace conductivity model.For this simple model, the sensitivities calculated using this program should be equal tothe accurate sensitivities.Consider a homogeneous halfspace of conductivity 0.01 Sm’ having the same dimensions as the model shown in Fig. 5.1 and parameterised using the mesh in Fig. 5.2.Consider a controlled-source electromagnetic experiment in which an along-strike electricdipole source (a point source in three-dimensional space) is located at x3 = —11000m and= = 0. Suppose measurements are made of E at x0 = l000m and y0 = = 0, andCHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d S 3d 84at a frequency of 0.2 Hz. The accurate sensitivities for this arrangement were calculatedby both the brute-force (or perturbation) method using the forward-modelling programof Unsworth, Travis & Chave (1993) and the adjoint-equation method (Unsworth & Oldenburg 1995). The sensitivities calculated using the latter method are shown in Figs. 7.laand 7.lb. The sensitivity is plotted as a constant value over the area of the cell to whichit refers.The program written to compute approximate sensitivites for the 2.5d problem wasused to calculate sensitivities for the homogeneous halfspace. The experimental configuration was the same as that described above. The adjoint field was that due to a y-directedunit electric dipole and was computed in a layered halfspace using the method describedin Appendix A. The resulting sensitivities are shown in Figs. 7.lc and 7.ld. There is veryclose agreement between the accurate sensitivities calculated using the adjoint-equationmethod and the sensitivities calculated using the approximate-sensitivity program (whichshould equal the accurate sensitivities for this special case of the homogeneous halfspace).7.3.2 Two-dimensional modelApproximate sensitivities were calculated for the two-dimensional model shown inFig. 5.1. The controlled-source experimental configuration used in the previous sectionwas also considered here. The accurate sensitivities calculated using the adjoint-equationFigure 7.1 (Following page) Sensitivities for the 2.5d problem described in section 7.3.1. The conductivity model was a homogeneous halfspace of 0.01 5m’.The source and observation locations are indicated by the open and solid triangles respectively. Panels (a) and (c) show log10Oln Ej/Oln °i, and panels (b)and (d) show log10 0 ç/0ln crj. IE and q are the amplitude and phase, respectively, of the electric field. The phase is measured in radians. The colour-barrefers to all four panels. Panels (a) and (b) were produced by the adjointequation method, and panels (c) and (d) were produced by the approximatesensitivity program using an adjoint field computed in a layered halfspace.CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3dC (‘2 CD C — C) CD CD Ca iD C U C CD —. CD — 1’-(‘1 (‘2 C) ) CD CD CD CDI I I I I I—5000—10000—1500085Q Co oo oCD 0— Cv—200000 0 0 0 0o o 0 0o o 0 0CD 0 CD 0— — Cv—20000—5000—1000015000c.)CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 86method are shown in Figs. 7.2a and 7.2b. The approximate sensitivities calculated usingan approximate adjoint field in a layered halfspace are shown in Figs. 7.2c and 7.2d.There is good agreement between the approximate and accurate sensitivities shown inFig. 7.2. Because these sensitivities are with respect to the logarithm of the conductivity,the direct effect of the conductivity structure is clearly evident, just as for the two-dimensional magnetotelluric sensitivities shown in Chapter 6. The agreement betweenthe approximate and accurate sensitivities for the 2.5d problem is better than that forthe two-dimensional magnetotelluric sensitivities. This is because the adjoint field is nowfully three-dimensional and so is even more dominated by its decay away from the dipolesource compared to its decay away from the line source in two dimensions.7.4 Computation times for the 2.5d problemThe computation times for the accurate and approximate sensitivities for a small2.5d example based on the model in Fig. 5.1 are listed in Table 7.1. These times are forobservations of the amplitude and phase of E at six observation locations at a frequencyof 0.2 Hz and for the 800 cells in the model. For this small example, the computation ofthe approximate sensitivities was quicker by a factor of five than the most efficient methodof calculating accurate sensitivities. This factor will increase further as the size of theproblem increases, just as for the magnetotelluric sensitivities discussed in Chapter 6 (seesection 6.4).Figure 7.2 (Following page) 2.5d sensitivities for the conductivity model shownin Fig. 5.1. The experimental configuration was the same as that considered insection 7.3.1. The source and observation locations are indicated by the openand solid triangles respectively. Panels (a) and (c) show log1öln E/Olnoi,and panels (b) and (d) show log108q5/Olnu. The colour-bar refers to all fourpanels. Panels (a) and (b) were produced by the adjoint-equation method, andpanels (c) and (d) were produced by the approximate-sensitivity program usingan adjoint field computed in a layered halfspace.CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3dC) D.? f) C’- o . CO CD U) C)C I! C IC) C CD CD ‘—uC’2 CQ CO Co U) IL) CD CDI I I I I I I Ip2000015000100005000—5000—10000—15000—2000087?jjWIz(m) z(m)CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 88Accurate sensitivities: brute-force 5000 mmadjoint-equation 25Approximate sensitivities: layered halfspace 5Table 7.1 Example computation times for the accurate and approximate sensitivities for the 2.5d inverse problem.7.5 Calculation of the sensitivities for the three-dimensional problemThe general expression for the sensitivity given by eq. (4.14) applies without simplification to the three-dimensional problem. It is usual for the model to be parameterisedinto cuboidal cells of uniform conductivity. For this situation, the sensitivity with respectto the conductivity of the th cell is given by= JET.Edv, (7.9)where F is any component of the electric or magnetic field, and V is the volume of theth cell. Most forward-modeffing programs produce the values of the electric field, E, atthe vertices of the cells. For E in this form, the integral in eq. (7.9) can be evaluatedby tn-linearly interpolating the real and imaginary parts of the components of B overthe volume of the cell (an extension to three dimensions of the method described insection 6.2.2):E (at, y, z) = c1 xyz + c2 xy + c3 yz + c4 xz + c5 x + c6 y + c7 z + c8, (7.10)where B represents the real or imaginary part of a component of the electric field. The 8 x8 system of equations relating the values of B at the vertices of the cell to the coefficientsc, i = 1,... , 8, was solved for a general cell using an analytic software package (WolframResearch, Inc. 1992). The adjoint field, E, is also tn-linearly interpolated over the cellin an analogous manner. With both B and Et in the form of eq. (7.10), the integration ineq. (7.9) can be performed analytically for a general cell to give an expression involvingCHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 89the coefficients of the tn-linear expansions. These analytic formulae are then used toevaluate the integral in eq. (7.9) for a particular cell from the appropriate nodal valuesof E and E.7.6 Comparison of approximate and accurate three-dimensional sensitivities7.6.1 Homogeneous halfspaceA program was written to compute the approximate sensitivities for the three-dimensional problem. The program was tested by comparing the sensitivities it produced for a homogeneous conductivity model with those calculated using the brute-forcemethod. The forward-modeffing program “SAMAYA” (Gupta, Raiche & Sugeng, 1989)was used for the brute-force calculation of the sensitivities. A grounded wire sourceextending from (50, —5, 0) m to (50, 5,0) m, and measurements of the vertical component of H at (0, 750,0) m were considered. The conductivity model was a homogeneoushalfspace of 0.01 Sm’. The region extending from x = 400 to 650 m, from y = —100to 150 m and from z 0 to 250 m was discretised into 5 x 5 x 5 cuboidal cells as shown inFig. 7.3. The conductivity of each of these 125 cells was perturbed in turn by 1%. Theresulting brute-force sensitivities for measurements of the amplitude and phase of H ata frequency of 100 Hz are illustrated in Figs. 7.4a and 7.4b. The sensitivities shown arethose for the shaded cells in Fig. 7.3. The sensitivity is plotted as a constant value overthe cell to which it refers.Sensitivities were calculated for the homogeneous halfspace using my approximatesensitivity program. The model parametenisation and experimental configuration werethe same as described above. The forward-modelled electric field, E, was computed usingSAMAYA (Gupta et al., 1989). The adjoint field was computed in a layered halfspace usingthe method described in Appendix A. Because observations of the vertical component ofthe magnetic field were being considered, the source for this adjoint field was a verticalmagnetic dipole. The resulting sensitivities are shown in Figs. 7.4c and 7.4d. There isFigure 7.3 Geometry for the three-dimensional example considered in sections 7.6 and 7.7. Panel (a) is the side view and panel (b) is the plan view.S indicates the y-directed grounded wire source and R indicates the observation location. The shaded cells are the ones for which the sensitivities are shownin Fig. 7.4. The numbers correspond to distances in metres.good agreement between these sensitivities and those calculated using the brute-forcemethod.Figure 7.4 (Following page) Sensitivities for the three-dimensional problemdescribed in section 7.6.1. The conductivity model was a homogeneous halfspaceof 0.01Sm1 Panels (a) and (c), and the colour-bar on the left of the figure, showlog10 IOln IHI/ôlnoi, and panels (b) and (d), and the colour-bar on the rightof the figure, show log10i8q/ôln. HI and 4’ are the amplitude and phase,respectively, of the vertical component of the H-field. The phase is measuredin radians. The sensitivities shown are those for the shaded cells in Fig. 7.3.Panels (a) and (b) were produced by the brute-force method, and panels (c)and (d) were produced by the approximate-sensitivity program using an adjointfield computed in a layered halfspace.CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d 3d 90a____z650OS50400b750250Os-100z xv4400R150650abN—6.20—6.35—6.49—6.64—6.78—6.92—7.07—7.21—7.35—7.50Ca0 LI) LI) 0—5.40—5.45—5.49—5.54—558—5.62—5.67—5.71—5.75—5.80y(m)y(m)CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 927.6.2 Three-dimensional modelBrute-force and approximate sensitivities were also calculated for a simple three-dimensional conductivity model. The survey configuration and model parameterisationwere the same as those in section 7.6.1. The conductivity model comprised a conductiveblock of 0.1 Sm’ in a background of 0.01 Sm’. The conductive block coincided withthe central square of nine cells in the third plane of cells down from the surface and thecentral square of nine cells in the fourth plane of cells below the surface. The brute-forcesensitivities are shown in Figs. 7.5a and 7.5b, and the approximate sensitivities calculated using the approximate adjoint field in a layered halfspace are shown in Figs. 7.5cand 7.5d. There is, in general, good agreement between the brute-force and approximatesensitivities shown in Fig. 7.5, although there are differences between the sensitivities ofthe phase for the cells outside the conductive block, especially in the second row (seeFigs. 7.5b and 7.5d). However, the resolution of the comparison is poor because of thelimited number of cells in the model, making it difficult to determine the significance orcause of these differences.7.7 Computation times for the three-dimensional problemThe computation times for the sensitivities calculated in sections 7.6.1 and 7.6.2 aregiven in Table 7.2. Six observation locations were considered at the single frequency of100 Hz. The computations were carried out on a Sun Sparcl0 workstation. An estimatedFigure 7.5 (Following page) Sensitivities for the three-dimensional problemdescribed in section 7.6.2. The conductivity model comprised a homogeneousbackground of 0.015m1with a conductive block of 0.15m coinciding with thecells of enhanced sensitivity shown here. Panels (a) and (c), and the colour-baron the left of the figure, show log10Oln jHI/ölno, and panels (b) and (a), andthe colour-bar on the right of the figure, show log10 ô qf/öln oi. The sensitivitiesshown are those for the shaded cells in Fig. 7.3. Panels (a) and (b) were produced by the brute-force method, and panels (c) and (d) were produced by theapproximate-sensitivity program using an adjoint field computed in a layeredhalfspace.ab—4.60—4.74—4.87—5.00—5.14—5.27—5.40—5.54—5.67—5.80CdI—5.10—5.44—5.77—6.10—6.44—6.77—7.10—7.44—7.77—8.10y(m)y(m)CHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 94time for the adjoint-equation method is also given. This estimate is equal to the timerequired to carry out six forward modeffings for this test-case using SAMAYA. Even forthis small example the approximate sensitivities are almost two orders of magnitudequicker to compute than the accurate sensitivities. This time difference will increase asthe number of observation locations increases for the same reason as that described insection 6.4.Accurate sensitivities: brute-force 260 mmadjoint-equation 12 mmApproximate sensitivities: layered halfspace 10 sTable 7.2 Example computation times for accurate and approximate sensitivities for the three-dimensional inverse problem.It is noteworthy that the approximate sensitivities for the three-dimensional problemare quicker to compute than those for the 2.5d problem for the same number of observationlocations and frequencies, and for the same number of layers in the model. To understandthis, consider the way in which the adjoint field is computed for the two problems (seeAppendix A). For the three-dimensional problem, the two-dimensional inverse Fouriertransform that gives the spatial dependence of the adjoint field is converted to a one-dimensional Hankel transform using eq. (2.10) of Ward & Hohmann (1988). This is thenevaluated using the digital filtering program of Anderson (1979b). For the 2.5d problem,the k and k dependencies are considered separately, as described in section 7.2. Theforward-modelled electric field and the adjoint field are combined in the (x, k) domain.The adjoint field is therefore required as a function of x and for a sequence (typicallyten) of values of k?,,. To obtain the adjoint field as a function of x for each value of k, aninverse cosine/sine transform has to be carried out. This is done using the digital filteringtechnique of Newman, Hohmann & Anderson (1986). The time required to evaluate oneHankel transform and one cosine/sine transform is similar. Because of this, and becausethe evaluation of the Hankel and cosine/sine transforms accounts for the majority of theCHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 95computation times for the two methods, the sensitivities for the 2.5d problem are roughlyten times slower to compute than those for the three-dimensional problem.The number of cells in the example considered in sections 7.6.1 and 7.6.2 is tinycompared to the number that would be required in a three-dimensional inversion of a realistic data-set. The number of cells required would be of the order of iOO. However, thenumber of cells used to produce the sensitivities shown in Figs. 7.4 and 7.5 was the largestfor which the brute-force sensitivities could feasibly be calculated. A parameterisationinvolving 8 x 8 x 8 cells was investigated, but a single forward modeffing for this numberof cells required 75 minutes. Carrying out 513 such forward modellings to calculate thesensitivities was considered impracticable. In addition to computation time, the amountof memory required by three-dimensional forward-modeffing programs severely constrainsthe number of cells that can be considered. For example, the model containing 8 cellsrequired 70 Mbytes of memory. This is approaching the maximum amount of memorytypically available on today’s computers.7.8 ConclusionsPrograms were written to compute approximate sensitivities for the 2.5d and three-dimensional inverse problems. For the 2.5d problem, the approximate sensitivities fora simple conductivity model were compared with the accurate sensitivities. The agreement was found to be better than that for the two-dimensional magnetotelluric sensitivities discussed in Chapter 6. The approximate sensitivities for the 2.5d problem shouldtherefore enable an iterative, linearised inversion procedure to converge to an acceptablemodel. Also, for the small example considered in this chapter, the computation of theapproximate sensitivities was quicker by a factor of five than the most efficient method ofcalculating the accurate sensitivities. This factor will increase further as the size of theproblem increases, just as for the magnetotelluric sensitivities discussed in Chapter 6.For the three-dimensional problem, the approximate sensitivities were computedfor a simple three-dimensional conductivity structure and compared to the brute-forceCHAPTER 7: APPROXIMATE SENSITIVITIES FOR 2.5d & 3d 96sensitivities. There was good agreement between the approximate and brute-force sensitivities. However, the number of model parameters for which the brute-force sensitivitiescould be computed was severely restricted by the computation time required by thethree-dimensional forward-modelling program. The comparison of the sensitivities wastherefore of a much lower resolution than for the two-dimensional magnetotelluric and2.5d problems.Although the example for the three-dimensional problem was severely restrictedin size, it did show that the approximate sensitivities for this problem are considerablyquicker to compute than the accurate sensitivities. The difference was roughly two ordersof magnitude for the small example presented in this chapter. This time difference willincrease significantly as the number of observation locations increases.Chapter 8SummaryThe multi-dimensional inversion of geophysical electromagnetic data is too computationally demanding to be carried out on a routine basis, even using present-day computertechnology. The aim of this thesis was to contribute a means of accelerating the multidimensional inversion process.In Chapter 2, I described an iterative, linearised inversion procedure for invertingtime-domain electromagnetic measurements for a one-dimensional conductivity model.This procedure is typical of the inversion strategy that one would ultimately like toapply to the multi-dimensional inversion of electromagnetic data. The most computationally intensive components of this inversion procedure are the forward modelling, thegeneration of the Jacobian matrix of sensitivities, and the solution of the linear systemof equations. In this thesis, I concentrated on developing a quick, approximate methodfor calculating the sensitivities.The approximate form of sensitivity that I presented is based on the adjoint-equationmethod. This method is the most efficient way of accurately calculating the sensitivitieswhen the inversion strategy is based on an under-determined problem. As described inChapter 4, the vast majority of the computation time required by the adjoint-equationmethod is due to calculating the adjoint field, E, in the multi-dimensional conductivitymodel for dipole sources at every observation location. The approximate sensitivitiesare obtained by replacing the true adjoint field by an approximate adjoint field that ismuch quicker to compute. Three possible forms of the approximate adjoint field wereinvestigated: the electric field computed in a homogeneous halfspace for the appropriatedipole source, the electric field computed in a horizontally-layered halfspace, and thatcomputed using the Born approximation.97CHAPTER 8: SUMMARY 98In Chapter 5, I investigated the three different forms of approximate adjoint field inthe context of the purely two-dimensional problem. It was shown that all three formsreplicated the dominant behaviour of the true adjoint field. The approximate adjointfield computed in a layered halfspace matched the features in the true adjoint field dueto one-dimensional structures in the conductivity model. The approximate adjoint fieldscomputed in the homogeneous and layered halfspaces were considerably quicker to compute than the true adjoint field. The Born approximation of the adjoint field, however,was slower. Although the computation time of the sensitivities using the Born approximation could be made smaller than that for the accurate sensitivities, the computationtimes for the sensitivities calculated using the approximate adjoint fields in a homogeneous or layered halfspace were even quicker. The approximate adjoint field in the layeredhalfspace is therefore the most useful of the three possible forms investigated here.In Chapter 6, I investigated the approximate sensitivities for the two-dimensionalmagnetotelluric problem. I compared approximate and accurate sensitivities for a simpleconductivity model. The level of agreement was good, and the approximate sensitivitiesusing the approximate adjoint field in either a homogeneous or layered halfspace wereconsiderably quicker to compute than the accurate sensitivities. The approximate sensitivities were used within an iterative, linearised procedure to invert both a syntheticdata-set and a field data-set. The inversion strategy was the same as that used for theone-dimensional inversion of time-domain electromagnetic data in Chapter 2, except thata subspace methodology was used to solve the linear system of equations. The approximate sensitivities performed successfully in the inversion of the synthetic data-set. Thefield data-set was a sub-set of the COPROD2 data-set which has proved challenging forthe many inversion procedures that have been applied to it. The inversion using theapproximate sensitivities did not achieve the final target misfit. However, the misfit thatwas achieved was close to the target value and the major features in the data were reproduced. The final model also contained the multiple mid-crustal conductors that otherinvestigations have produced. The inversion using the approximate sensitivities thereforeperformed as well as any of the other inversion procedures applied to the COPROD2 data.CHAPTER 8: SUMMARY 99To conclude the work for this thesis, approximate sensitivities were computed forthe 2.5d problem and three-dimensional controlled-source problem. These situations,especially the latter, are much more common in geophysics than the special case oftwo-dimensional magnetotellurics. These situations are also significantly more time-consuming. The approximate sensitivities for the 2.5d and three-dimensional problemswere computed for simple conductivity models and compared to the accurate sensitivities.The agreement for the 2.5d problem was found to be better than for the two-dimensionalmagnetotelluric problem. The agreement for the three-dimensional problem was alsofound to be good, although the number of model parameters used for the comparison wasseverely restricted by the computation time required by the three-dimensional forward-modelling program. Since the level of agreement between the approximate and accuratesensitivities was found to be similar, if not better, than that for the two-dimensionalmagnetotelluric problem, it was assumed that these approximate sensitivities would perform just as well in an iterative, linearised inversion procedure. The more importantresult from the comparisons of the approximate and accurate sensitivities was the difference in the computation times. This difference was considerable, especially for thethree-dimensional problem where a difference of two orders of magnitude was found,even for the small example that was considered. Also, the relative differences in computation times between the approximate and accurate sensitivities for the two-dimensional,2.5d and three-dimensional problems all increase as the number of observation locationsincreases.In this thesis, therefore, I have presented an approximate form of sensitivity that isappropriate for the multi-dimensional inversion of any geophysical electromagnetic dataset. The approximate sensitivities appear to be sufficiently good approximations of theaccurate sensitivities to enable an iterative, linearised inversion procedure to convergeto an acceptable model. And, as desired, the approximate sensitivities are considerablyquicker to compute than the accurate sensitivities, especially for the three-dimensionalcontrolled-source problem.ReferencesAnderson, W.L., 1979a. Programs TRANS_HCLOOP and TRANS_HZWIRE: calculationof the transient horizontal coplanar loop soundings and transient wire-loop soundings,USGS Open-File Report 79-590.Anderson, W.L., 1979b. Numerical integration of related Hankel transforms of orders 0and 1 by adaptive digital filtering, Geophysics, 44, 1287—1305.Asten, M.W., 1987. Full transmitter waveform transient electromagnetic modeling andinversion for soundings over coal measures, Geophysics, 52, 279—288.Bailey, R.C., 1970. Inversion of the geomagnetic induction problem, Proc. Roy. Soc.London, A315, 185—194.Boerner, D.E. & Holladay, J.S., 1990. Approximate Fréchet derivatives in inductiveelectromagnetic soundings, Geophysics, 55, 1589—1595.Boerner, D.E. & West, G.F., 1989. A spatial and spectral analysis of the electromagneticsensitivity in a layered earch, Geophys. J. mt., 98, 11—21.Brant, A.A., 1966. Introduction, in Mining Geophysics, Vol. 1, pp. 3—6, Soc. of Expl.Geophys., Tulsa.Brewitt-Taylor, C.R. & Weaver, J.T., 1976. On the finite difference solution of two-dimensional induction problems, Geophys. J. R. astr. Soc., 47, 375—396.Cagniard, L., 1953. Basic theory of the magneto-telluric method of geophysical prospecting, Geophysics, 18, 605—635.Constable, S.C., Parker, R.L. & Constable, C.G., 1987. Occam’s inversion: A practicalalgorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289—300.deGroot-Hedlin, C. & Constable, S., 1990. Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data, Geophysics, 55, 1613—1624.Dosso, S.E. & Oldenburg, D.W., 1991. Magnetotelluric appraisal using simulated annealing, Geophys. J. mt., 106, 379—385.Ellis, R.G., Farquharson, C.G. & Oldenburg, D.W., 1993. Approximate Inverse Mappinginversion of the COPROD2 data, J. Geomag. Geoelectr., 45, 1001—1012.Ellis, R.G. & Oldenburg, D.W., 1994. The pole-pole 3-D DC-resistivity inverse problem:a conjugate-gradient approach, Geophys. J. mt., 119, 187—194.Everett, M.E. & Edwards, R.N., 1992. Transient marine electromagnetics: the 2.5-Dforward problem, Geophys. J. mt., 113, 545—561.100REFERENCES 101Farquharson, C.G. & Oldenburg, D.W., 1993. Inversion of time-domain electromagneticdata for a horizontally layered Earth, Geophys. J. mt., 114, 433—442.Fullager, P.K. & Oldenburg, D.W., 1984. Inversion of horizontal loop electromagneticfrequency soundings, Geophysics, 49, 150—164.Gradshteyn, I.S. & Ryzhik, I.M., 1994. Table of Integrals, Series, and Products, 5th ed.,Academic Press, Boston.Gill, P.E., Murray, W. & Wright, M.H., 1981. Practical Optimization, Academic Press,London.Grant, F.S. & West, G.F., 1965. Interpretation Theory in Applied Geophysics, McGraw—Hill, New York.Golub, G.H. & Van Loan, C.F., 1989. Matrix Computations, The John Hopkins University Press, Baltimore.Gupta, P.K., Raiche, A.P. & Sugeng, F., 1989. Three-dimensional time-domain electromagnetic modeffing using a compact finite-element frequency-stepping method, Geophys. J. mt., 96, 457—468.Heiland, C.A., 1940. Geophysical Exploration, Prentice—Hall, New York.Hohmann, G.W., 1988. Numerical Modeling for Electromagnetic Methods of Geophysics,in Electromagnetic Methods in Applied Geophysics, Vol. 1, pp. 313—363, ed. Nabighian,M.N., Soc. Expl. Geophys., Tulsa.Hohmann, G.W. & Raiche, A.P., 1988. Inversion of Controlled-Source ElectromagneticData, in Electromagnetic Methods in Applied Geophysics, Vol. 1, pp. 469—503, ed.Nabighian, M.N., Soc. Expl. Geophys., Tulsa.Ingham, M.R., Jones, A.G. & Honkura, Y., (eds.) 1993. MT-DIW1 Special Section, J.Geomag. Geoelectr., 45, 931—1150.Jones, A.G., 1993. The COPROD2 Dataset: Tectonic Setting, Recorded MT Data, andComparison of Models, J. Geomag. Geoelectr., 45, 933—955.Jones, F.W., 1974. The Perturbation of Geomagnetic Fields by Two-dimensional andThree-dimensional Conductivity Inhomogeneities, Fag eoph, 112, 793—800.Jupp, D.L.B. & Vozoff, K., 1975. Stable Iterative Methods for the Inversion of Geophysical Data, Geophys. J. B. astr. Soc., 42, 957—976.Jupp, D.L.B. & Vozoff, K., 1977. Two-dimensional magnetotelluric inversion, Geophys.J. B. astro. Soc., 50, 333—352.Kaufman, A.A. & Keller, G.V., 1983. Frequency and Transient Soundings, Elsevier,Amsterdam.Lanczos, C., 1961. Linear Differential Operators, Van Nostrand, London.REFERENCES 102Lee, K.H. & Morrison, H.F., 1984. A solution for TM-mode plane waves incident on atwo-dimensional inhomogeneity, Lawrence Berkeley Lab., Rep. LBL-1 7857.Livelybrooks, D., 1993. Program 3Dfeem: a multidimensional electromagnetic finiteelement model, Geophys. J. mt., 114, 443—458.M8Donald, K.J. & Agarwal, A.K., 1994. Two and Three dimensional induction modelling on massively parallel supercomputers, presented at the J2th Workshop on EMInduction in the Earth, Brest, France.McGillivray, P.R. & Oldenburg, D.W., 1990. Methods for calculating Fréchet derivativesand sensitivities for the non-linear inverse problem: a comparative study, Geophys.Prosp., 38, 499—524.McGiffivray, P.R., Oldenburg, D.W., Ellis, R.G. & Habashy, T.M., 1994. Calculation ofsensitivities for the frequency domain electromagnetic problem, Geophys. J. mt., 116,1—4.Mackie, R.L. & Madden, T.R., 1993. Three-dimensional magnetotelluric inversion usingconjugate gradients, Geophys. J. mt., 115, 215—229.Mackie, R.L., Madden, T.R. & Wannamaker, P.E., 1993. Three-dimensional magnetotelluric modeling using difference equations - Theory and comparisons to integralequation solutions, Geophysics, 58, 215—226.Madden, T.R., 1972. Transmission systems and network analogies to geophysical forwardand inverse problems, Tech. Rep. 72—3, Dept. of Earth and Planetary Sciences, MIT.Madden, T.R., 1990. Inversion of low frequency electromagnetic data, in Oceanographicand Geophysical Tomography, pp. 379—408, eds. Desaubies, Y., Tarantola, A. & VinnJustin, J., Elsevier, Amsterdam.Madden, T.R. & Swift, C.M., 1969. Magnetotelluric studies of the electrical conductivitystructure of the crust and upper mantle, in The earth’s crust and upper mantle, Geophys. Monogr. 13, pp. 469—479, ed. Hart, P.3., Am. Geophys. Union, Washington,D.C.Menke, W., 1984. Geophysical Data Analysis: Discrete Inverse Theory, Academic Press,Orlando.Nabighian, M.N., 1988. (ed.) Electromagnetic Methods in Applied Geophysics, Vol. 2,Soc. Expl. Geophys., Tulsa.Nabighian, M.N. & Macnae, J.C., 1991. Time-domain electromagnetic prospecting methods, in Electromagnetic Methods in Applied Geophysics, Vol 2., pp. 427—479, ed.Nabighian, M.N., Soc. Expl. Geophys., Tulsa.Newman, G.A. & Alumbaugh, D.L., 1995. Frequency-domain modeling of airborne electromagnetic resonses using staggered finite differences, Geophysical Prospecting (accepted).REFERENCES 103Newman, G.A., Hohmann, G.W. & Anderson, W.L., 1986. Transient electromagneticresponse of a three-dimensional body in a layered Earth, Geophysics, 51, 1608—1627.Oldenburg, D.W., 1983. Funnel functions in linear and non-linear appraisal, J. geophys.Res., 88, 7387—7398.Oldenburg, D.W. & Ellis, R.G., 1993. Efficient inversion of magnetotelluric data in twodimensions, Phys. Earth Planet. mt., 81, 177—200.Oldenburg, D.W. & Li, Y., 1994. Subspace linear inverse method, Inverse Problems, 10,915—935.Oldenburg, D.W., McGifflvray, P.R. & Ellis, R.G., 1993. Generalised subspace methodsfor large scale inverse problems, Geophys. J. mt., 114, 12—20.Oristaglio, M.L. & Worthington, M.H., 1980. Inversion of surface and borehole electromagnetic data for two-dimensional electrical conductivity models, Geophys. Prosp.,28, 633—657.Parker, R.L., 1977. Understanding inverse theory, Ann. Rev. Earth planet. Sci., 5,36—64.Parker, R.L., 1980. The inverse problem of electromagnetic induction: existence andconstruction of solutions based on incomplete data, J. geophys. Res., 85, 4421—4428.Parker, R.L., 1994. Geophysical Inverse Theory, Princeton University Press, Princeton.Poddar, M., 1982. A rectangular loop source of current on a two-layered Earth, Geophys.Prosp., 30, 101—114.Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 1992. NumericalRecipes in FORTRAN: The Art of Scientific Computing, 2nd ed., Cambridge UniversityPress, Cambridge.Roach, G.F., 1982. Green’s Functions, Cambridge University Press, Cambridge.Rodi, W.L., 1976. A Technique for Improving the Accuracy of Finite Element Solutionsfor Magnetotelluric Data, Geophys. J. R. astr. Soc., 44, 483—506.Ryu, J., Morrison, H.P. & Ward, S.II., 1970. Electromagnetic Fields About a LoopSource of Current, Geophysics, 35, 862—896.Schultz, A., Kurtz, R.D., Chave, A.D. & Jones, A.G., 1993. Conductivity discontinuitiesin the upper mantle beneath a stable craton, G.R.L., 20, 2941—2944.Scriba, H., 1974. A Numerical Method to Calculate the Electromagnetic Field of aHorizontal Current Dipole, Pageoph, 112, 801—809.Smith, J.T. & Booker, J.R., 1988. Magnetotelluric inversion for minimum structure,Geophysics, 53, 1565—1576.REFERENCES 104Smith, J.T. & Booker, J.R., 1991. Rapid Inversion of Two- and Three-DimensionalMagnetotelluric Data, J. geophys. Res., 96, 3905—3922.Telford, W.M., Geldart, L.P. & Sheriff, R.E., 1990. Applied Geophysics, 2’ ed., Cambridge University Press, Cambridge.Torres-Verdin, C. & Habashy, T.M., 1993. Cross-Well Electromagnetic Tomography,Expanded Abstract, 3 International Congress of the Brazilian Geophysical Society,Rio de Janeiro.Unsworth, M.J. & Oldenburg, D.W., 1995. Subspace inversion of electromagnetic data -application to mid-ocean ridge exploration, Geophys. J. mt. (submitted).Unsworth, M.J., Travis, B.J. & Chave, A.D., 1993. Electromagnetic induction by a finiteelectric dipole source over a 2-D earth, Geophysics, 58, 198—214.Ward, S.H., 1990. (ed.) Geotechnical and Environmental Geophysics, Vols. 2 & 3, Soc.Expl. Geophys., Tulsa.Ward, S.H. & Hohmann, G.W., 1988. Electromagnetic theory for geophysical applications, in Electromagnetic Methods in Applied Geophysics, Vol. 1, pp. 131—311, ed.Nabighian, M.N., Soc. Expl. Geophys., Tulsa.Weaver, J.T., 1994. Mathematical Methods for Geo-Electromagnetic Induction, JohnWiley, New York.Weidelt, P., 1972. The inverse problem of geomagnetic induction, Z. Geophys., 38, 257—289.Weidelt, P., 1975. Inversion of two-dimensional conductivity structures, Phys. EarthPlanet. mt., 10, 282—291.Whittall, K.P. & Oldenburg, D.W., 1992. Inversion of Magnetotelluric Data for a One-Dimensional Conductivity, Geophysical Monograph Series, 5, Soc. Expi. Geophys.,Tulsa.Wiggins, R.A., 1972. The general linear inverse problem: Implication of surface wavesand free oscillations for Earth structure, Rev. Geophys. Space Phys., 10, 251—285.Wolfram Research, Inc., 1992. Mathemaiica, version 2.1.Wu, F.T., 1968. The inverse problem of magnetotelluric sounding, Geophysics, 33, 972—979.AppendicesAppendix A: Computation of the Electric and Magnetic Fields in aLayered Halfspace due to a Dipole SourceThe electric or magnetic field generated in a layered halfspace by an electric ormagnetic dipole source is required throughout this thesis: in Chapter 2 the verticalcomponent of the magnetic field is required on the surface of a layered halfspace for ahorizontal electric dipole source, and in Chapters 5 to 7 the approximate form of adjointfield that is found to be the most useful is the electric field generated in a layered halfspace.In this appendix, I outline the method that was used to calculate the electric and magneticfields generated in a horizontally-layered halfspace by an electric or magnetic dipolesource.The layered haifspace in which the electric or magnetic field is to be calculated ismade up of uniform layers of constant conductivity as shown in Fig. Ad. The TE andTM mode Schelkunoff potentials of Ward & Hohmann (1988) are used:A = Ae, (A.1)F = Fe, (A.2)where ê is the unit vector in the z-direction. In the j” layer, A and F satisfy thefollowing ordinary differential equations (section 4, Ward & Hohmann 1988):( — ) A = o,( - = o,where the tilde represents the two-dimensional Fourier transform, = k + k — i€W2 +iw,uoj, w is the angular frequency, p and e are the magnetic permeability and electrical105y xcp0•Figure A.1 Notation and coordinate system for the horizontally layered conductivity model of the Earth used in this thesis. is the depth to the bottomof the th layer, and and are the conductivity and thickness, respectively,of the th layer. S represents the dipole source.permittivity, and is the conductivity of the th layer. The definition of the two-dimensional Fourier transform used by Ward & Hohmann is also used here:andaDO aDO(k,k,z) = / / F(x,y,z) e_i k) dxdy, (A.5)i-CO i-DODO DOF(x,y,z)=U1APPENDIX A: EM FIELDS IN A LAYERED HALFSPAGE 106Sz=-ho= 0z=0z’z2zi-’tiziZN 2UN1 tN-1ZN1P(k,k,z) e”’ d1cdk. (A.6)APPENDIX A: EM FIELDS IN A LAYERED HALFSPACE 107The solutions for and areA(k,k,,z,w) = C(k,k,w) e z_1) + D(k,k,w) e_uj_zj_1), (A.7)= U(k,k,w) euj z_i) + Vj(k,k,w) e_Ui_Zj_). (A.8)In the region of free space above the layered halfspace, A and F satisfyA0(k,k,z,w) = Co(kz,ky,w)eu0z + Do(k,k,w)e_u0z, (A.9)Po(k,k,z,w) = U0(kz,ky,w)eu0z + V(kz,ky,w)e_u0z, (A.10)where ug = k + k.Firstly, consider the potential F. From Ward & Hohmann eqs. (1.152) and (1.153),the boundary conditions on F at z = arek, z=z_i,w) = k, z=z_i,w), (A.11)oP.____a2’k, z=z_1,w) = ‘ (ku,, k, z_—z_1,w). (A.12)Substituting eq. (A.8) into the above boundary conditions gives+ V = Ui_ie_1ti_1 + j_ieUj_huj_1, (A.13)uU — uVj = uj_iUi_ieui_lti1 — u_iVj_ieui_1ti_1 (A.14)= — is the thickness of the J’’ layer. These two boundary conditions can becombined in a matrix equation:(1 1 (U= ( e1ti1 e_Ui_lti_1 (A 15u—Ui) ) U_i&L_1t_1 _u_ei1i1) a_i)This can be rewritten as(v’) = e’ M (), (A.16)APPENDIX A: EM FIELDS IN A LAYERED HALFSPA GE 108where/ (1 + ___)e_2u_1ti_1 (1 —M — (A 17(‘—ti) (l+L)In eq. (A.16), the exponential term exp(u1t_has been factored out of the matrixMso that the propagation through the layers carried out below remains stable. Applyingthe boundary conditions at the surface of the halfspace and using eq. (A.lO) for thepotential F in the free space above the layered halfspace leads to() =where/(1i!a’ ‘i—-’_I\ u01 \ u0’1‘Al- (1-) (i+’))uo UoEqs. (A.18) and (A.16) can be used to propagate the boundary conditions through thelayers to produce an expression relating the coefficients of the potential P in the freespace above the halfspace to those in the basement halfspace:N+1 N N+1(E.°)= () ex(>uiti) ]JM, (E11). (A.20)There can be no upward-decaying solution for F in the basement halfspace, so UN+l = 0.It is assumed, at the moment, that the source is at a height z = —h above the layeredhalfspace. If this source is a unit x-directed electric dipole then, from eq. (4.137) ofWard & Hohmann,v—_____—uoh A2lo— 2k+k . ( . )For a unit x-directed magnetic dipole,v—-- k —uoh2(A.22)APPENDIX A: EM FIELDS IN A LAYERED HALFSPACE 109using eq. (4.106) of Ward & Hohmann. Eq. (A.20) is therefore a pair of simultaneousequations in the two unknowns U0 and VN+l. Once U0 is known, eq. (A.18) can be usedto obtain U1 and V1. Successive applications of eq. (A.16) then give the coefficients Uand Vj, and hence the potential F, in every layer.A modified version of the above approach was used for the potential A. The boundaryconditions on A at z = z_1 are, from Ward & Hohmann eqs. (1.182) and (1.183),A(k, k, z=z_1, w) = A1(kr, k, z=z_1, w), (A.23)1 oA. 1 oA.=(k,,k,z=z_i,w). (A.24)Because of the awkwardness, numerically, of the second boundary condition at the surfaceof the layered halfspace (above which o = 0), the source is assumed to lie within layer 1(at z = h). The final result for the source on the surface of the halfspace is obtainedby letting the source approach the surface from below. In layer 1, therefore, eq. (A.7) ismodified to contain a term representing the particular solution:A1 = G eU1Z + D1 e_U1Z + A e_U1_l. (A.25)Extending the analysis of Ward & Hohmann to determine explicit expressions for theirquantity A, the appropriate expression for A eq. (A.25) for a unit c-directed electricdipole is= 2 k+k’(A.26)andA — [LW — iW/2U1 ik A 27p ——2u k+k’for a unit x-directed magnetic dipole. The boundary conditions can now be propagated through the layers in a similar manner to the potential F to obtain two simultaneous equations in the two unknowns C0 and DN+l (D0 is zero since there is nowno downward-decaying part of the solution in the free space above the halfspace). ToAPPENDIX A: EM FIELDS IN A LAYERED HALFSPAGE 110obtain the potentials for y-directed electric and magnetic dipoles, the transformation(x,y) —* (y,—x) [ —* (k,_k)] can be used.Once the potentials A and F are known, the components of the electric and magneticfields can be calculated throughout the layered halfspace using eqs. (1.129) and (1.130)of Ward & Hohmann:E= a2A — (A.28)zo+iew 9xOz= 1 82A+ , (A.29)cr + iew Oyôz72=. (- + k ) A, (A.30)c7+iEW \UZ Jand8A 1 L92FH = + -—-- , (A.31)Oy iwu 9x8zOA 1 02F (A.32)OX iW[L oyaz72H = (-- + k ) F, (A.33)ZWI \OZ /where k2 = — iwuo. The above equations are in fact transformed to the (k,k)domain and used to obtain B and H from A and F, and then the transformation backto the (x,y) domain that is appropriate for the particular problem is carried out. Forthe purely two-dimensional case discussed in Chapters 5 and 6, the electric field for thetwo-dimensional dipole sources is obtained by setting the along-strike wavenumber, k,equal to zero before carrying out the inverse cosine/sine transform to recover the xdependence of the electric field. For the 2.5d case discussed in Chapter 7, the inversecosine/sine transform with respect to the across-strike wavenumber, k, is performed atdifferent values of the along-strike wavenumber, k. The resulting (z, k )-dependence ofthe electric field is then used in the calculation of the sensitivities. Finally, for the threedimensional case also described in Chapter 7, the two-dimensional Fourier transformwith respect to k and k is converted to a Hankel transform using eq. (2.10) of Ward &APPENDIX B: PURELY 2d SENSITIVITIES 111Hohmann, and this Hankel transform performed to obtain the full spatial dependence ofthe electric field. The numerical evaluation of the cosine/sine transforms was carried outusing the method of Newman, Hohmann & Anderson (1986), and the evaluation of theHankel transform was carried out using the method of Anderson (1979b).Appendix B: Adjoint-equation method for the purely two-dimensionalinverse problemFor a purely two-dimensional problem, that is, one in which both the conductivitymodel and the source are invariant in the strike direction, the adjoint field requiredto calculate the sensitivities becomes the electric field due to a two-dimensional dipolesource. This can be seen as follows. Consider the particular form of eq. (4.14) for a givenmeasurement and assume that the model and ID are invariant in the strike direction whichis chosen here as the y-direction. One obtains= f J J E(z,y,z) .E(z,z) (x,z) dxdydz (B.1)9o.j—00 —00=jj(x,z) E(x,z).{f00Ef(x,Y,z) dy} dxdz. (B.2)where F represents the component of interest of the electric or magnetic field. Without loss of generality, the adjoint field can be expressed as a two-dimensional Fouriertransform:E(x,y,z) = J J Et(k,k,z) dkdk. (B.3)-00 -00APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 112If this expression for the adjoint electric field is substituted into eq. (B.2), then theintegration with respect to y can be carried out, since= S(k). (B.4)This reduces the adjoint field to only its zero along-strike-wavenumber component:= jfbj(xz) E(a,z).{_j dk} dxdz. (B.5)The term within the braces is the electric field due to a two-dimensional source, that is,one that is invariant in the strike direction. Eq. (B.5) can therefore be rewritten as= j E(x, z) . E(x, z) (x, z) ds, (B.6)where the adjoint field, Et, is now due to a two-dimensional dipole source at the observation location (x0,z0). This is the appropriate expression for the sensitivities for a purelytwo-dimensional inverse problem.Appendix C: Electric fields due to two-dimensional dipole sourceson a homogeneous halfspaceThe first form of approximate adjoint field presented in Chapter 5 for the purely twodimensional inverse problem is the electric field generated in a homogeneous halfspace bya two-dimensional dipole source. To calculate the corresponding approximate sensitivitiesin Chapter 6 for the two-dimensional magnetotelluric inverse problem, this approximateadjoint field is required for both electric and magnetic dipole sources oriented in both thea- and y-directions on the surface of the halfspace.APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 113Electric dipole sourceTo develop an expression for the electric field generated in a homogeneous halfspaceby an x- or y- directed two-dimensional unit electric dipole, initially follow the derivation ofKaufman & Keller (1983) for the electric field generated by a point (in three-dimensionalspace) electric dipole. Consider an electric dipole oriented in the y-direction and supposethat both the dipole and the point at which the electric field is to be calculated are situatedwithin the conductive halfspace (see Fig. C.1). Assume that z is positive downwards.Kaufman & Keller make use of the vector potential A:A = (O,A,A) (C.1)such thatE = iwA + V(V . A). (C.2)w is the angular frequency, and [I and cr are the magnetic permeability and conductivityof the halfspace. The two non-zero components of the vector potential are:=_ j {e’l + J0(Ap) dk (C.3)= j u2+A e_UJ1(p) dA, (C.4)where u2 = )2 + iw,ucr and p2 = x2 + y2. I have followed the convention that Ward &Hohmann (1988) use for the time- to frequency-domain Fourier transform. This accountsfor the plus sign in the above expression for u. h is the z-coordinate of the dipole source.Using eq. (2.10) of Ward & Hohmann (1988) to convert the above Hankel transformsto two-dimensional Fourier transforms gives1 °° °°1 —= J J — { Iz-hI + + e_u} dk dk (C 5)=— L L A(u+ A) e_U ik e ky) dk dk, (C.6)APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 114o=Oz=OSz=h 00•Figure C.1 The coordinate system and geometry for computing the electricfield induced in a homogeneous halfspace by a two-dimensional dipole source, S.where now A2 = k + k. Substituting these two expressions into eq. (C.2) gives the threecomponents of the electric field resulting from the y-directed finite electric dipole:1 J kk (e’ u A82 _ u+A_u(z+h))+ e+k1cuuA(u + A)e} dk dk, (C.7)1 J J°° f (iwp — k (e_uIz_ + U —82 — u u) u+A2ku— A(u + A)e} dkdk, (C.8)and1 fofoIik /= — (u sgn(z — h) e_uIz_hI + u(u — A) e_u(z+h))_cru u+A2ik AV U+ ( + A) e(z+h)} dk dk. (C.9)APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 115Appendix B shows that the adjoint field required to calculate the sensitivities for a purelytwo-dimensional inverse problem can be obtained from the electric field for a finite dipolesource by considering only the zero along-strike-wavenumber component. If this is done inthe above equations the following expressions are obtained for the electric field generatedin a homogeneous halfspace by a two-dimensional y-directed electric dipole source:= 0, (0.10)B =— f e_UZ dlc , (0.11)2r J_u+kj= 0. (0.12)where it has been assumed that the source is at the surface of the halfspace (h = 0)and now u2 = k + iwo. From the above equations it is clear that for this type ofsource there is only an along-strike component of the electric field. This is exactly whatis required when calculating the sensitivities for the E-polarisation mode of the two-dimensional magnetotelluric inverse problem. Note that eq. (0.11) agrees, as it should,with eq. (4.208) of Ward & Hohmann for the electric field due to an infinite line current.To obtain expressions for the electric field produced by an x-directed two-dimensionalelectric dipole first consider the electric field produced by a finite x-directed electric dipole.This field is given by eqs. (0.7) to (0.9) after rotation of the coordinate axes correspondingto the transformation (x,y) —* (y, —x). This implies that (1c,k) —÷ (ks,, —kr) and(Br, B, E) —* (Er, —Er, Es). The electric field due to the two-dimensional dipole sourceis then obtained by considering only the components of the equations corresponding tozero along-strike wavenumber:___jU(e’ + e_u) dk, (0.13)= 0, (0.14)= f k (sgn(z — h) eZ + e_u(z+h)) dk. (0.15)APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 116The electric field for an x-directed two-dimensional electric dipole is therefore restrictedto lie in the plane perpendicular to the strike direction. The integrals in the aboveequations can be evaluated (Lee & Morrison 1984) to give analytic expressions for thefield components:= —- j- K(ik’r) +z)K1(ikr), (C.16)= K0(ikr) + --K1(ikr), (C.17)where k2 = —i wjio- and r2 = x2 + z2 for the electric dipole on the surface of the halfspace.K0 and K1 are the zeroth and first order Modified Bessel functions of the second kind.Magnetic dipole sourceTo develop expressions for the electric field induced in a homogeneous halfspaceby a two-dimensional magnetic dipole source initially follow the analysis of Ward &Hohmann (1988). Consider the Schelkunoff potentials A and F such thatE— 1 ÔAZ OFZ C18OxOz182A OFz z C19yuOyOz Ox’=+ k2) A, (C.20)where k2 = —iwiu. For a unit x-directed magnetic dipole= j°° / (e’ — e_u) ei(kx+kv) dk dk, (C.21)andF = i_: L (€_ulz_h + e_u) dk dk. (C.22)APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 117Here, as before, u2 = k + k + iw,uo and )2 = k + k. Substituting these expressionsinto eqs. (0.18) to (0.20) givesJ f (sgnz — h) — 1]e_uIz_ + 2A e u(z+h))8712 — j_ U +x dk dk, (0.23)/= j j { (sgn(z — + e_u)k2 U— e_u) } dk dk, (0.24)+ (eu1zh1 + +andj°° °° ik= 8 2 j ___(e_ul7_I — e_u(z+h)) dk dk. (0.25)71 j_ - UTo obtain expressions for the electric field generated by a two-dimensional z-directedmagnetic dipole on the surface of the halfspace consider the reduced forms of the aboveequations for k = 0:(0.26)iwI f Ue_UZ dk, (0.27)= 2irju+k= 0. (0.28)For a y-directed two-dimensional magnetic dipole use the transformation (x, y) —* (y, —z)in eqs. (0.23) to (0.25) before considering the reduced form of these equations. This gives= f e — u z e z ai, (0.29)271= 0, (0.30)= 0. (0.31)APPENDIX C: 2d DIPOLES ON A HOMOGENEOUS HALFSPACE 118Using eq. (3.914) of Gradshteyn & Ryzhik (1994), eq. (C.29) can be reduced towpkz=— —K1(zkr) (C.32)where k2 = —iw,uo and r2 2 +The electric fields described above were computed by evaluation of the Bessel functions or, if the expression for the electric field could not be reduced to one involvingBessel functions, by evaluation of the Fourier transform using the digital filtering code ofNewman, Hohmann & Anderson (1986).Appendix D: Sensitivities for the magnetotelluric apparent resistivityand phaseIn the magnetotelluric inverse problem, the data are usually considered to be valuesof the apparent resistivity, pa, and phase, q, wherep(w) = --- and (w) = phase (i). (D.1)E and H represent orthogonal horizontal components of the electric and magnetic fields.w is angular frequency and u is magnetic permeability. Since it is the apparent resistivityand phase that are the data in the inverse problem, the sensitivities are required for thesedata rather than for the electric and magnetic fields. However, the sensitivities for theapparent resistivity and phase can be obtained from those for the electric and magneticfields as follows. Differentiating the apparent resistivity in eq. (D.1) with respect to themodel parameter o gives(D28o wi H Ocr \jHj)APPENDIX D: SENSITIVITIES FOR MT 119— 2 E fi OEI IE OH D 3— w H UH H2 (Note that all the quantities in the above expression, the electric and magnetic fields aswell as the sensitivities, are to be evaluated at the observation location.Since the electric field in the frequency domain is a complex quantity it can bewritten asE = Ee4’E. (D.4)Treating both Ej and 4E as functions of the model parameters, differentiating eq. (D.4)with respect to o gives-= E —- (i) + e (D.5)Ocr ôcr Ocrj= iEea +EOjE{ (D.6)8cr E &r1 ÔE___1 ÔIEI= ä=ô (D.7)Equating the real and imaginary parts of eq. (D.7) gives___= () (D.8)and___= m (). (D.9)This argument can also be applied to the magnetic field resulting inôIH= Hj Re (D.1O)9uj \H 9c7jJand= m (). (D.11)APPENDIX D: SENSITIVITIES FOR MT 120Eqs. (D.8) and (D.10) can be substituted into eq. (D.3) to give a final expression for thesensitivity for the apparent resistivity:—- D12— w H UHI E Ou) HI H o)j .= 2Pa {e () - e £Z) }. (D.13)This expression can be used to calculate the sensitivity for the apparent resistivity sincePa is known at the observation location, and OE/Ou and OH/ôo can be calculated usingeq. (4.14).To determine the final expression for the sensitivity for the phase consider the ratioB — IEIei4’ — (D14HIHIe’H HBy definition,= —(D.15)and so—__— (D16)oo.j— oo= m () - m (E (D.17)using eqs. (D.9) and (D.11). The sensitivity for the phase can therefore be calculatedusing eq. (D.17) and eq. (4.14).