UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analysis of geomagnetic depth sounding data Stinson, Kerry James 1981-03-26

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

Item Metadata


831-UBC_1981_A6_7 S85.pdf [ 8.51MB ]
JSON: 831-1.0052943.json
JSON-LD: 831-1.0052943-ld.json
RDF/XML (Pretty): 831-1.0052943-rdf.xml
RDF/JSON: 831-1.0052943-rdf.json
Turtle: 831-1.0052943-turtle.txt
N-Triples: 831-1.0052943-rdf-ntriples.txt
Original Record: 831-1.0052943-source.json
Full Text

Full Text

ANALYSIS OF GEOMAGNETIC DEPTH SOUNDING DATA by KERRY JAMES STINSON B.Sc, Simon Fraser University, 1973 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Geophysics and Astronomy) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October, 1981 (c) Kerry James Stinson, 1981 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of (Si^-pA^s^c s QLW</ A sfv o *o i*i<j The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date Oct. X%/% > Abstract The electromagnetic induction problem is non-linear, and thus is very difficult to solve for all but the simplest symmetries. Because of this, quantitative modelling of the conductivity structure from geomagnetic depth sounding data is expensive and time consuming, and the possibility that the anomaly is produced by channelling of regionally induced currents may invalidate the results. For this reason traditional methods of analysis are generally qualitative in nature, with quantitative information estimated on the basis of simplified models of the anomaly. The theory and assumptions used in these traditional methods are studied in this thesis, and the range of their applicability is investigated. To avoid the current channelling complication, and to also get a linear relation between the model and the data, the problem is reformulated with the subsurface current density as the model parameter, rather than the conductivity. The disadvantage of this formulation is that models that fit the data are very non-unique. The character of this non-uniqueness has been explored using Backus-Gilbert appraisal, and by the construction of unconstrained models. The results indicate that reasonable resolution of the true model's horizontal features is possible, but that vertical resolution will be lacking. To i counter this, the infinite range of possible models is constrained by introducing expected physical features of the true model into the model construction algorithm.' This construction algorithm was tested using data generated from a variety of artificial models, and was successful in resolving both the horizontal and vertical positions of the major features in all of them. The algorithm was then used to determine the subsurface current structure for real data taken across the Cascade anomaly in Washington State. Table of Contents Page Abstract ii List of Tables viList of Figures viiAcknowledgements xiii Introduction 1 Chapter I: The Source Field 6 1.1 General Nature of Sources1.2 Mathematical Models of the Source Field 17 Chapter II: The Electromagnetic Induction Problem 26 2.1 General Source Field Over a Homogenous Half-Space 26 2.2 Uniform Source Field Over a Two-Dimensional Earth: Forward Modelling 48 Chapter III: Traditional Methods of G.D.S. Analysis 58 3.1 The Formulation and Separation of the Normal Field ....58 3.2 Visual Methods of Analysis 74 3.3 The Induction Tensor and Induction Arrows 78 3.4 Quantitative Methods Used in G.D.S 97 V Chapter IV: The Current Density Model 107 4.1 The Current Density Formulation 104.2 Uniqueness and Backus-Gilbert Appraisal 113 4.3 Construction: The Smallest Model 129 4.4 Constrained Model Construction Using Linear Programming 142 Chapter V: Analysis of G.D.S. Across the Cascade Anomaly ....160 Conclusions 188 Bibliography 192 Appendix A: Maxwell's Equations in a Conductor 212 Appendix B: Correlation of Price's Induction Solutions with Plane Wave Solutions 219 Appendix C: The Uniform Field Assumption 229 Appendix D: Separation of the External and Internal Fields ..234 Appendix E: Determination of the Induction Tensor Elements ..240 Appendix F: Properties of Current Distributions that Mimic a Line Current 246 vi Appendix G: Analytic and Numeric Integrations 251 Appendix H: Determining the Major Axis of an Ellipsoid 264 vii List of Tables Page Table 1.1 Skin Depths for a Homogenous Earth 8 Table 2.1 Parameter Ranges in G.D.S 39 Table 2.2 Values of £ for Varying T, A , with CT= .005 S/m 40 Table 2.3 Values of ft for Varying T, A , with (5 = .5 S/m 1 Table 2.4 Phase Difference Between the Anomalous Surface Magnetic Field and the Normal Surface Electric Field, for Both the Line Current and the Jones-Pascoe Results 57 Table 5.1 Station Positions ....182 vi i i List of Figures Page Fig. 1.1 The Steadystate Interaction Between the Earth's Dipole Field and the Solar Wind 10 Fig. 1.2 Relationship Between the Solar Wind-Magnetosphere Coupling Function and the Magnetospheric Substorm Index .. •• 13 Fig. 1.3 Polar Magnetic Substorm Models 5 Fig. 1.4 Geometry of the Source Field and the Earth 18 Fig. 1.5 Plane Wave Amplitude Spectrum of a Line Current Source 24 Fig. 2.1 The Ratios of the Complex Magnitudes of the Inducing and Induced Magnetic Field 37 Fig. 2.2 Three Representative Models Used for the Jones-Pascoe Forward Induction Program 53 Fig. 2.3 Results From the Jones-Pascoe Program for the Three Models of Fig. 2.2, at a Period of 5 Min 54 Fig. 2.4 Results From the Jones-Pascoe Program for the Three Models of Fig. 2.2, at a Period of 50 Min 55 Fig. 3.1 The Assumed Conductivity Model of G.D.S 59 ix Fig. 3.2 One Dimensional Surface Array Perpendicular to the Strike of a Two-Dimensional Earth 68 Fig. 3.3 Magnetometer Locations and Magnetograms for a Substorm in Aug., 1972 75 Fig. 3.4 Relation Between Induced Current and Inducing Field Presuming First Order Induction Only 77 Fig. 3.5 Contoured Amplitude and Phase of the Z Component of the Magnetograms of .Fig. 3.3, at a Period of 68.3 Min 77 Fig. 3.6 The Direction of the In-Phase Induction Arrow Near a Line Current Running North-South 88 Fig. 3.7 First Order Induction in Crossing Conductive Paths 89 Fig. 3.8 The Effects of Current Channelling through a Conductive Sphere 92 Fig. 3.9 A Variety of Two-Dimensional Current Density Models 9 Fig. 3.10 The Vertical Fields of the Current Models of Fig. 3.9 9Fig. 4.1 The Presumed Two-Dimensional Model 108 Fig. 4.2 Averaging Functions Calculated for the Two-Dimensional Current Density Problem.' 126 X Fig. 4.3 Determination of the Appropriate Weighting Factor to Offset Geometric Decay 133 Fig. 4.4 L^-Norm Smallest Model Construction 135 Fig. 4.5 L^-Norm Smallest Model Construction to Test Horizontal Resolution .138 Fig. 4.6 L2_Norm Smallest Model Construction to Test Vertical Resolution 140 Fig. 4.7 L,-Norm Linear Programming Model Construction 149 Fig. 4.8 L|-Norm Linear Programming Model Construction to Test Horizontal Resolution 152 Fig. 4.9 L| -Norm Linear Programming Model Construction to Test Vertical Resolution 153 Fig. 4.10 L|-Norm Linear Programming Model Construction Using Vertical Current Dike 154 Fig. 4.11 Lj-Norm Linear Programming Model Construction Using Dipping Current Dike 155 Fig. 4.12 Li-Norm Linear Programming Model Construction Using Jones-Pascoe Induction Program Output for a Single Conductive Path 157 Fig. 4.13 L(-Norm Linear Programming Model Construction Using Jones-Pascoe Induction Program Output for Two Vertically Separated Conductors 158 xi Fig. 5.1 Array of Magnetometers Used by Law et al to Study the Cascade Anomaly 161 Fig. 5.2 In Phase and Quadrature Arrows for Data from Array of Fig. 5.1 162 Fig. 5.3 The Array of Magnetometers which Measured the Data Used in this Thesis 163 Fig. 5.4 Magnetograms for the Magnetic Storm of Feb. 1980 164 Fig. 5.5 Magnetograms for the Magnetic Storm of Feb. 1980 165 Fig. 5.6 The Amplitude Spectra of the Anomalous Field at ORT 1 68 Fig. 5.7 The Amplitude Spectra of the Anomalous Field at MUD 169 Fig. 5.8 The Amplitude Spectra of the Anomalous Field at GRE 170 Fig. 5.9 The Degree of Polarization Between the Three Directional Components of the Anomalous Field 174 Fig. 5.10 The Ratio of the Smallest to the Second Smallest' Eigenvalues of the Polarization Ellipsoid 177 Fig. 5.11 The Induced Field from Two Conductors at Different Depths 179 xii Fig. 5.12 The Estimated Dip and Strike of a Line Current Model 180 Fig. 5.13 Lj-Norm Linear Programming Current Models Using the Real Data 183 Fig. 5.14 A Final Map of the Cascade Anomaly 187 Fig. B.1 Rotating to a New Horizontal Coordinate System 223 Fig. C.1 The Propogation Angle of the Transmitted Wave into a Half-Space of Uniform Conductivity 231 Fig. G.1 The Parameterized Current Density Model 261 xi i i Acknowledgements The stolid brick exterior of the Geophysics and Astronomy building at U.B.C. belies the warmth and enthusiasm that I have encountered within. I thank all the people of the department for the enjoyably conducive atmosphere that has sustained me during this thesis. My supervisor, Doug Oldenburg, deserves special mention; his insight, knowledge, and infinite eagerness have been drawn on throughout this thesis". Also worthy of note are the many long discussions I have . had with my officemate Tim Scheuer, on all aspects of life and geophysics. I would also like to thank Peter Fullagar for his suggestions on the inversion methods. I thank Gerald Hensel and Dr. John Booker of the University of Washington Geophysics Department for kindly supplying me with the Cascade anomaly data. I am immensely grateful to my wife Tracy for her support and patience during the course of this thesis. With luck she can take the candle out of the window now. 1 INTRODUCTION Geomagnetic depth sounding (G.D.S.) is one of the many electromagnetic methods used to determine the conductivity structure of the earth. In this particular method the electromagnetic measurements used are the surface values of the three directional components of the magnetic field. G.D.S. in practice is further subdivided into problems of global and'local extent, with the first evaluating the earth's gross conductivity structure using the assumption of radial symmetry, and the second determining more local structures using the plane earth assumption. This thesis will be concerned only with the local G.D.S. problem, so henceforth, unless specified otherwise, G.D.S. will refer only to this subdivision of the problem. In G.D.S. problems, the sources (or primary fields) for the electromagnetic induction into the earth have traditionally been natural large scale variations in the magnetosphere. Several types of natural source have been used in the past, as will be outlined in Chapter 2, and the ones most commonly used for G.D.S. studies will be discussed in somewhat more detail, also in that chapter. All electromagnetic induction problems are identical in initial formulation; from this point the intractability of the general problem has led to the branching into separate methods 2 along paths of different simplifying assumptions. In the magnetotelluric method, the assumption is usually that of a one-dimensional plane earth, with conductivity varying only with depth. In this one-dimensional case, the forward problem of electromagnetic induction can be solved quickly and cheaply, for both the plane layered case, (Cagnaird, 1953; Keller and Frischnecht, 1966) and for the continuous case (Oldenburg, 1979). The inverse problem is also quite well solved, and a variety of different methods are in the literature (Cagnaird, 1953; Becher and Sharpe, 1969; Bailey, 1970; Oldenburg, 1979; Fischer et al, 1980). The forward problem of the two-dimensional induction problem is more difficult and more expensive, although numerical solutions do exist (Jones and Price, 1970; Jones, 1970; Madden and Swift, 1969). There is also some ambiguity about the correctness of the solutions in each case, and the fastest way to do the problem numerically (Praus, 1975). Correspondingly, although methods for inversion of electromagnetic surface readings over two-dimensional structures do exist, (Jupp and Vozoff, 1977), they are very expensive and their iterative procedures can fail because of lack of convergence (Jupp and Vozoff ,1977). In general, most model constructions for the two-dimensional case are simply forward modelling attempts, with starting models either adjusted interactively (Porath et al, 1970), or adjusted randomly in Monte Carlo fashion (Anderssen ,1975; Woods, 1979) until a satisfactory fit with the 3 data is achieved. The three-dimensional problem easily involves an order of magnitude increase in difficulty over the two-dimensional case. The solution of the forward problem is an enormously expensive proposition (Lee et al, 1981), and to this author's knowledge there are at this time no working methods for true three-dimensional inversion. In this thesis, the three-dimensional case will not be considered. The assumption of vertical one-dimensionality is probably very good in some areas, such as in the deep ocean away from the ridges (Oldenburg, 1981) or in sedimentary basins (Vozoff, 1972). However, in many cases the reason an area will be of interest will be because of its anomalous electromagnetic responses, with this anomalous behavior implying an underlying structure more complex than the one- dimensional case. It is to handle these cases that the methods constituting G.D.S. have been devised. Unlike some electromagnetic induction methods, such as the global G.D.S. problem in which the assumption of radial symmetry is made, or controlled source studies in which the primary field is known, there are no universally applied assumptions in G.D.S. Rather, the avoidance of the intractability of the general induction problem is attained simply by concentrating on qualitative results, with little emphasis placed on rigorous quantitative evaluation. The intuition needed to arrive at these qualitative results 4 is generally obtained by reviewing the solutions for simple, but basic, forward problems. For this reason , Chapter II gives the results for induction of a general source field over a homogenous half space (in section 2.1) and also gives the results for three idealized two-dimensional structures which should be quite representative of the range of two-dimensional anomalies possible (in section 2.2). A review of the traditional methods used in geomagnetic depth sounding is presented in Chapter III, with many of the discussions drawing upon the results of Chapter II. It will be seen that both the qualitative and quantitative information obtained using these methods is crucially dependent on the horizontal variations in the source field, as well as upon induction outside the measurement array. Both of these factors are difficult to estimate well (Gough, 1973; Lilley, 1975). As well, many of the final qualitative results obtained from these methods are based on the intuition of the individual authors, where often it would be preferable to have an automated method to avoid individual biases. To avoid these problems, a different method of analysis of G.D.S. data for presumed two-dimensional anomalies is suggested in Chapter IV. In this approach the current density distribution in the earth is treated as the desired model rather than the conductivity structure. By formulating the problem to consider only the results of the induction process, (that is, the induced currents in the earth) rather than the total induction process, 5 the difficulties due to the non-uniformity of the source field and the non-local induction effects are avoided. As well, the problem in this formulation is linear, simplifying manifestly the inversion procedure. On the other hand, by ignoring the primary field altogether a great deal of model-limiting information is thrown away, and so we expect the non-uniqueness of our results to be increased. In fact, using only the two-dimensional assumption, the current density formulation will be shown to be hopelessly non-unique (section 4.2 and 4.3). However, by presuming certain physical characteristics expected of natural two-dimensional current structures, and incorporating these as constraints in a model construction routine (section 4.4) this non-uniqueness can be greatly reduced. In the final chapter (Chapter V) of this thesis, the new formulation described in Chapter IV is applied to real data taken along a linear array in southwestern Washington, over the Cascade Anomaly. Throughout this thesis, the attempt is made to present the material in as readable a fashion as possible. For this reason proofs or mathematical developments which are not considered essential for the continuity of the discussion are either referenced, if possible, or relegated to the Appendix. 6 Chapter I The Source Field I.1 General Nature of Sources The sources utilized for geomagnetic depth sounding are the naturally occurring variations in the earth's magnetic field. These can include the variations in the earth's internal field as well as the variations and disturbances in the external magnetosphere, but the most commonly employed are the latter. Because of the unpredictability of most magnetospheric disturbances, the 'choice' of source to be used is usually limited to post-procurement editing of data. For this reason most field experiments entail continuous recording of magnetic data for long periods, of up to months. The selection of source type will depend on the size and position of the array, as well as upon the depth one wishes to 'see' into the earth. In general, analysis of. the results will be simplified if the horizontal wavelengths of the source field are much larger than the dimensions of the array. Thus, for large arrays (200 km. by 500 km.) the sources used will be the large scale disturbances 7 such as the Dst, or polar magnetic substorms. For smaller arrays, it might be possible to use the more localized source fields, such as micropulsations. The desired depth of penetration would determine the range of initial amplitude of the source required, as well as the dominant period range in the source. The skin depth formula for a homogenous earth is: d = f T where T = period d = depth at which field original magnitude O = conductivity -2. If a conductivity of 10 S/m is assumed, which is perhaps a reasonable average for crustal and mantle conductivities (Brace, 1971), then the approximate skin depths can then be evaluated for the dominant period ranges (Garland, 1979, pg. 257)) of selected magnetospheric disturbances, (see Table 1.1). (1.1.1) has decayed to 1/e of 8 Table 1.1: Skin Depths for a Homogenous Earth Magnetospheric Di sturbance Dominant Period Skin Depth (d) Lightning -a 10 s. 0.159 km. Micropulsation 10 - 10s. 1.59-50.3 km Polar Substorm 1 01s. 159 km. Dst 4 10 s. 503 km. Diurnal Variations Solar Lunar 24 hrs. 25 hrs. 1480 km. 1540 km. Diurnal variations, long period (up to one year) external variations, and even the internally origined secular variations, have been used for global studies of the earth's gross radial conductivity structure (Chapman, 1919; Chapman and Price, 1930; 9 Lahiri and Price, 1939; Rikitake, 1950; Runcorn, 1955; McDonald, 1957; Eckhardt et al, 1963; Banks, 1969). However, this thesis deals only with the local G.D.S. problem, so attention will be focussed mainly on those disturbances with period ranges enabling investigation to depths of 500 km. Although micropulsations fit into this category they are not a common source in G.D.S. As well, they are diverse and complex in nature (Jacobs, 1970); for these reasons they will not be discussed here. All the external sources to be ' considered result from disturbances in the steady-state interaction between the earth's dipole field and the solar wind. In this steady state, or quiet time, the magnetic field frozen into the highly conducting plasma of the solar wind compresses the earth's dipole field (see Fig. 1.1). A shock front of highly compressed field lines marks the boundary between the regions of influence of the interplanetary magnetic field (IMF) and the earth's dipole field, and acts as an effective 'shield', excuding the ions of the solar wind from the earth (Nishida, 1978). However, it has long been accepted that a portion of these energetic charged particles must be finding entry into the earth's ionosphere (Rostoker, 1972). This was experimentally confirmed by the correlation between satellite measurements of the ion flux in the ionosphere and observed acitivity on the sun (Garland, 1979, pg. .253). The mechanism for the ion entry is still not well 10 Fig. 1.1 The steady state interaction between the earth's dipole field and the solar wind. The dotted line indicates the magnetopause, inside which the earth's field is confined. 11 understood or agreed upon, but it has become apparent that the direction and magnitude of the components of the IMF are important controlling parameters in the process (Akasofu, 1979). Satellite measurements of the IMF and coincident surface readings of the earth's field have indicated that a southward directed IMF is probably a neccessary (although not sufficient) condition for substorm activity (Rostoker, 1972). Perreault and Akasofu quantified this by showing that an empirical relation could be found between the IMF and the development of geomagnetic storms (Perreault and Akasofu, 1978; Akasofu, 1979). The rate of energy dissipation, u(t), was evaluated for 15 major geomagnetic storms, using measurements of ring current and aurora particle injection, and estimates of Joule dissipation in the ionosphere. It was found that this estimated energy dissipation u(t) could be closely duplicated by an empirically determined function €(t), dependent only on interplanetary parameters: €(t) sin (Joules/sec) (1.1.2) with: V(t) speed of the solar wind plasma B(t) total magnitude of the IMF 12 0(t) : Tan"'(IBy/Bzl) for Bz > 0 180 -Tan"1(IBy/Bzl) for Bz < 0 J?0 : estimate of linear dimension of the cross-section of the magnetopause (assumed = 7Re) The right hand side of eguation 1.1.2 is closely related to the Poynting vector flux of the interplanetary electromagnetic field, so that 6(t) can be regarded as the rate of energy coupling between the solar wind and magnetosphere (Akasofu, 1979). It is noted that if the IMF has only a northward component, then the coupling should be zero, and that for a completely southward directed IMF the coupling is a maximum VB*j£*. The correlation of fc(t) with substorm activity (as given by the substorm index AE) has subsequently been found to be quite good for values of €.(t) less than 10 J/s (Akasofu, 1979). (see Fig. 1.2). Once inside the magnetosphere, the bulk of the solar ion flux is carried to the magnetotail (Rostoker, 1972; Akasofu, 1979), with a certain amount of these charged particles spiralling along the earth's dipole field lines. At the poles the convergence of the field lines results in the particle's reflection so that they travel from pole to pole (Akasofu and Chapman, 1961). The centrifugal force due to the curvature of the field lines, as well as the inhomogeneity of the earth's field (Alfen, 1950,pg.14-23; Akasofu and Chapman, 1961) results in an eastward drift of the electrons in their pole to pole 13 ) Fig. 1.2 The relationship between the solar wind-magnetosphere coupling function and the magnetospheric substorm index AE for a storm in the month of July, 1974 (after Akasofu (1979)). 1 4 travel, and a westward drift of the protons. This drift gives the effect of a ring current at equatorial latitudes at about two or three earth radii (Akasofu and Chapman, 1961), with the result of this ring current being a dipole field which opposes the earth's internal dipole field. Low values of energy input into the magnetopause ( 6 < 10'" J/s ) correspond to extremely quiet conditions, with 'normal' energy dissipation processes in the magnetosphere maintaining a steady state. When the value of € increases to a critical value of io" J/s , it is suggested that the normal dissipative modes cannot handle the increased rate of energy being transferred to.the magnetotail (Akasofu,1979). The result is that a portion of the currents in the magnetotail are then diverted along magnetic field lines to the polar ionosphere. Equivalent current systems for this process have been suggested by a number of authors using ground based measurements, (Birkeland, 1908; Bostrom, 1964; Kisabeth and Rostoker, 1971; Kisabeth and Rostoker, 1977), with most recent models (for northern latitudes) featuring field-aligned currents (Birkeland currents) connected by a strong westward electrojet in the morning sector, and Birkeland currents connected by a strong eastward electrojet in the evening sector (see Fig. 1.3). This postulated diversion of magnetotail energy into currents in the polar ionosphere is generally accepted as the basic model for the polar substorm. It has been noted (Akasofu, 1979) that the strong correlation between the coupling rate €(t) and substorm 15 NOON (b) (a) «- J Fig. 1.3 Polar magnetic substorm models: (a) detailed model of auroral zone electrojets and connected field-aligned Birkeland currents from Rostoker (1978) (b) magnetospheric and field-aligned current model from Kisabeth (1975). 16 activity, both in growth and in decay (see Fig. 1.2) indicates that it is the energy coupling rate which controls all phases of the substorm. This is in contrast to the previous concept (for example, Rostoker, 1972) of the substorm as simply an energy unloading device activated when the energy density of the magnetotail reached some critical limit. In addition to the diversion of currents from the magnetotail to the polar ionosphere, there is also large scale injection of ions into the ring current system, with a resultant decrease in the main dipole field. This effect and its subsequent slow decay constitute the storm-time disturbance or Dst. The source most commmonly used for G.D.S. is the polar substorm. The main reasons for this are the strength (hundreds of nanotesla) and the uniformity (scale length of thousands of kilometers) of its disturbance field. In practice there are usually more than one overlapping substorms, with the final superposition constituting a polar storm. The entire storm is used as the source, with generally no attempt made to separate the individual substorms, or the Dst. I.2 Mathematical Models of the Source Field To facilitate a physical understanding of the induction process in the case of complex magnetospheric source fields, it is practical to mathematically formulate the true source field as a summation of elementary solutions of the wave equation, and then attempt the physical understanding in terms of a single solution. In the region 0 > z > -h, between the surface of the earth and the lowest current element of the source field (see Fig. 1.4), Maxwell's equations are: V-B = 0 (1.2.1) V-D = 0 (1.2.2) VxH = D (as 0*= 0) (1.2.3) 18 Fig. 1.4 Geometry of the source field and the earth 19 VxE = -B (1.2.4) Thus, presuming p = fi0 , and €=£<> everywhere in the region : > -h, using a time dependei of 1.2.3 and 1.2.4, we arrive at 0 > z  dependence of etW/t, and taking the curl V E = 9aE or VH = %H with: and V-E = 0 (1.2.5) (1.2.6) (1.2.7) (1.2.8) The solution of 1.2.5 (or 1.2.6) results in two independent solutions, corresponding to the TE mode (E4 = 0) and the TM mode (H4 = 0) (Budden, 1961, pg. 13-15,22-30). For the TE mode, the 20 elementary solution for the frequency w is: E(x,y,z,t) = Eo{ky,-kx,0}eLk«X eik»* e^V"* (1.2.9) and for the TM mode: H(x,y,z,t) = Ho{ky,-kx,0}e'^K e£K* V^V"' (1.2.10) where in both cases: ki +k! +k\c, = - <fc0 (1.2.11) Note that the parentheses {} denote the vector components. To completely describe the primary field of an arbitrary source would require summation of the elementary solutions of both the TE and TM modes over all possible k* and k^ values. (By virtue of 1.2.11, k2a is not independent). For simplicity in the discussion, consider a source that is producing TE mode waves only. (It will be shown in section 2.1 of Chapter II that this is the only mode that need be considered for induction in the earth). In this case the electric field vector of the primary field at any point (x,y) is: OO CO CU)f E(x,y,z,t) = etlA,t^ A(kx,ky) {ky,-kx,0} •oo -co 21 • e * e J J e dkxdky (1.2.12) The negative sign for the exponential involving the z component is chosen to ensure that the waves will always by propagating downward; k^ here is thus > 0 always. The x and y components of E are not independent because of 1.2.8, so the spectrum A(kx,ky) (denoted the 'angular spectrum' by Booker and Clemmow (1950)) can be found using the values of either Ex or Ey over an arbitrary plane surface (Booker and Clemmow, 1950). Using the z = 0 plane, dropping the time dependence for convenience, and considering Ex only, we have: 00 CO Ex(x,y,0) = ^ J A(kx,ky)ky etk** e^*1 dkxdky ~00"00 (1.2.13) This is a double inverse fourier transform, allowing us to write: A(kx,ky) = l/(4Trvky)^ Ex (x , y, 0) e^** e" 3^dxdy -co -co (1.2.14) It is apparent then, as shown by Booker and Clemmow (1950) and Wait (1954), that the angular spectra of complex sources can be calculated for known values of Ex on one plane, and then used for the calculation of the electric ( or magnetic ) field at any other position. It was indicated in section 1.1 of Chapter I that the equivalent model for a substorm would in some areas 22 resemble a horizontal line current; the angular spectrum for this model is thus of obvious interest and is calculated below. A line current of magnitude I in the y direction, oscillating at a frequency w, at a height h, has an associated electric field with only a component in the y direction (Landau and Lifschitz, 1960, pg.195; Wait, 1970, pg.23): Ey(x,z) = -i(jiewl)/(2fr) K0 {(<pa)'l[x*- + (z+hf ^} (1 .2.15) where K0 is the modified Bessel function of order zero. Expressing this as a summation of elementary solutions on an arbitrary plane z = constant: CO Ey(x,z) = j~A(kx)e"tkxXLdkx (1.2.16) (Note that as there can be no plane wave propagation in the y direction, then kx and k?0L are no longer independent in this case). The amplitude spectrum is given by: CO A(kx) = -(iHowI)/(41Vl-)|K0{ cPj'l[x2'+(z+h)1 ]'l]ei^K dx -oo = -(iju0wl )/[4Tr ( qVkx ) 3 e ( 1 .2.17) (from Gradshtyn and Rhysik , 1965, pg. 736) Thus, the final expression for Ey is: 23 oo Ey(x,z) = -(ho«I)/(4n)Je-^«+k^'UV^ -00 (1.2.18) For the line current example, a plot of the surface value of the elementary wave amplitude vs. the horizontal wavelength, 2 IT /kx, is given in Fig. 1.5 for a line current of period one hour, at a height of 100 km. It is important to note from this figure, that even for this very uncomplicated source an infinite superposition of waves of differing wavenumbers is produced. For each elementary solution (as denoted by equation 1.2.9 or 1.2.10) the total horizontal wavelength of the primary field will be given by: A = 2fT/(kxl + ky1 (1 .2.19) As seen from the line current example, most real fields will be composed of an infinite number of waves of different horizontal wavelengths. In G.D.S. the important wavelength will be the minimum value which still has a significant amplitude at the surface. This wavelength will be the value of the largest significant spatial non-uniformity, and is commonly called the scale length. It will be seen in Chapter 2 and Chapter 3 that the scale length of the primary field and its estimation are of fundamental importance in virtually all aspects of conventional G.D.S. analysis. The traditional range of values of X was given by Price (1962): 24 9.0 FIG. 1.5 Plane wave amplitude spectrum on the surface of the earth from a line current at height 100 km. The log of the amplitude is plotted against the log of the spatial wavelength X ( A in km.) at periods: 1 sec. (O) 20 hr. (A) 25 4x10**km. < }\ < 4x10s km. The maximum value corresponds to the circumference of the earth, whereas the minimum value was obtained from an estimate of the sharp drop-off point of the amplitude spectrum for a line current, as shown in Fig. 1.5. (Note that Price also considered a height of 100 km. for his line current, and used four times this height as his estimate of the cut-off point). Methods of experimentally estimating the scale length will be discussed in section two of Chapter III. 26 Chapter II The Electromagnet ic Induction Problem 2.1 General Source Field Over a Homogenous Half-Space To illustrate the major features of the induction process, including non-uniformity of the source field, the simplified case of an isotropic, homogenous earth is considered. As will be seen in Chapter III, the real conductivity structure of the earth will often be modelled as a normal, one dimensional structure containing a small anomalous region. To be able to intuit the variations in the 'normal' induced field of the one dimensional earth due to the inclusion of the anomalous portion, one must obviously first understand fully the normal field case. It is possible to solve for the normal field for both a layered earth (Cagnaird, 1953; Keller and Frischnecht, 1966) and a continuously varying earth (Oldenburg, 1979), but for simplicity only the homogenous earth case will be treated here. The extension of the theory to the full one-dimensional problem is relatively straightforward, and in any case, the major qualitative points of the discussion will be the same for both. 27 Much of the development that will be presented here is similar to that done by Price (1950,1962). However , slightly different assumptions will lead to a different physical description of the solutions and should provide more insight into the problem. It should be emphasized that the assumptions that Price makes are not invalid; it is simply that the form of his solutions lead to a conceptually different way of viewing the induction process. Presume two homogenous half-spaces as in Fig. 1.4 (disregarding the source field) , with the positive z direction downwards, and CT as denoted. As shown in Appendix A, for a time dependence of e^w' in a homogenous medium, Maxwell's equations can be put in the form: V E = <pE (2.1.1) where: 9 = iw^o-- vxfji0€.o (2.1.2) As well, E is non-divergent everywhere except at the boundary between the half-spaces (see Appendix A): V-E = 0 (2.1.3) 28 For the range of values of period range commonly used in second term in equation 2.1 comparison to the first term, space. In other words, the conductor will be negligible Using this and the zero conduct have for equations 2.1.1 and 2. 9 = 9a= -w^0£o z °" expected in the earth, and the G.D.S. (listed in Table 2.1), the .2 will be insignificant in within the conductive lower half displacement currents in the compared to the real currents, ivity in the upper half-space, we 1.2: < 0 (2.1.4) <f> = <f>e = iwju00- z > 0 (2.1.5) Using a separation of variables to solve equations 2.1.1 and 2.1.3, we presume a solution of the form: E(x,y,z) = Z(z)F(x,y)e''*ut where: (2.1.6) F(x,y) = {Fx,Fy,Fz} (2.1.7) (Note, the time term, eLV°^, will appear in all equations, so we 29 will drop it henceforth, resurrecting it only in the final solut ions). Using 2.1.6 in 2.1.1, with a separation constant - V* we get: ?- (^)/Z = -Vl= (0+ y^)/Fx (2.1.8) Similarly, using equation 2.1.6 to separate 2.1.3 with a separation constant - e<: (4JEL+ l£iL)/F2 = -(li)/E = -* (2.1.9) The solutions for the equations involving 2(z) in 2.1.8 and 2.1.9 are, respectively: E(z) = ei(^^^a (2.1.10) Z(z) = e** (2.1.11) 30 Obviously then, »<= MS)1*?)"1 (2.1.12) The left hand equation in 2.1.9 has two possible solutions: Solution 1 (2.1.13) and Fz = 0 (2.1.14) Solution 2 V1 3x (2.1.15) (2.1.16) where the second solution makes use of the last expression in equation 2.1.8. 31 It is noted that solution one requires the 'Z' component of the electric field to be zero; this corresponds to a TE mode solution (ie., E is in the boundary plane). On the other hand, upon taking the curl of the electric field for solution 2, it is found that the 'Z' component of the magnetic field is zero, corresponding to the TM mode solution (ie., H in the boundary plane). Thus, it is expected that the separation of the problem into solutions of Type 1 and 2 should be directly related to the separation of plane waves into TE and TM mode components, and this will be seen to be the case. In identical fashion to Price's development, the first type is solved by letting: Fx = >JfeQ (2.1.17) Fv = -1 ax (2.1.18) where P(x,y) is a scalar quantity. Thus, all terms on the right hand side of equation 2.1.8 lead to: + JlL + S)1P = o (2.1.19) 32 The function P(x,y) and the value of ^ can be shown to be the same in both half-spaces. Thus, the final form of the first type solution is: E,(x,y,z) . l^-$,0> U,e-<^\ B,e(*W'^ (2.1.20) for z < 0 B,(,,y.,) - {»-»?.. 0} C, ,-^*^"%* (2.1.21) for z > 0 with P(x,y) satisfying 2.1.19. (Note that in 2.1.21 the solution for 2(2) has only the negative exponential term to avoid unbounded solutions at 2 =*oo ). It should be pointed out here that each solution corresponding to a different "V (where S) can vary between 0 and oo) is as valid as any other; the total solution will be the summation (or integration) of all these 'elementary' solutions over the entire range of . Using Maxwell's equation relating H and the curl of E (A.13 in Appendix A), the magnetic field for the first type solution is found to be: 33 H,(x,y,z) - i(^+9J'V(vH -3' (2.1.22) for z < 0 .Cie-(^**.V"» (2.1.23) for z > 0 At the boundary of the half-spaces, z = 0, the tangential components of E, and all three components of H must be cont inuous. This leads to: B, = -A j• (1 - R)/(1 + R) (2.1.24) and C, = A,•(2R)/(1+R) (2.1.25) where: R = (2.1.26) 34 To get the same form as Price (1950), let: A, = A, 1 ( V + Va. ) /(Wjuo) (2.1.27) B, ' = -B, i ( Ol + qv )'7(wji0) (2.1.28) The magnetic field at the surface then becomes: H(x,y,o) = -{(A/ +Bl' )iL,(A; +B; )VL,(A; -b; )vs^y,} (2.1.29) and from 2.1.24, 2.1.27 and 2.1.28: B; = A,1 (1 - R)/(1 + R) (2.1.30) From the sign of the exponential in equation 2.1.22 (and using the value of <PA , (from equation 2.1.4) it is clear that A/ is the complex magnitude of the inducing or primary field, and B,1 is the complex magnitude of the induced or secondary field. It will be seen in Fig. 2.1 that the argument of (1-R)/(1+R) is always between 0 and This ensures that in equation 2.1.29 the magnitude of the horizontal components of the primary magnetic field will always be increased by the addition of the secondary horizontal field, whereas the magnitude of the primary 35 vertical field will always be decreased by the secondary vertical field. Substituting 2.1.4 and 2.1.5 into 2.1.26, the value of R is: R = [( S>v- w^o€e)/( iw/*6<r)],'^ (2.1.31) The separation constant, "v , will be shown later in this development to be equivalent to the,total horizontal spatial wavenumber, (kxl + ky1)'1, so that is directly related to the scale length, "X , discussed in section 1.2: > = 2TT/s> (2.1.32) Thus, the range of expected - values of "0 is readily obtained from the range of ?i arrived at in section 1.2; these are given in Table 2.1, along with the ranges of the other variables important in G.D.S. Even for the largest expected scale length, and for all periods greater than 5 seconds, it is found that: >)* » wl|uet0 (2.1.33) Thus, 2.1.31 can be approximated by: 36 R = I/O + ip )''*-(2.1.34) where: (2.1.35) If p, is very large, the magnitude of R (from equation 2.1.34) will be very small, so that from equation 2.1.30 |B|/A(| X 1. This corresponds to significant induction of the horizontal components into the earth, and almost total extinction of the vertical field. As well, the phase difference between A,' and B,' will be near zero. If j3 is very small, the magnitude of R will approach 1, so that | B,'/Aj is nearly zero, and B,1 is nearly ft/2 out of phase with A,'. This indicates near total reflection of the inducing field from the surface. The complete dependence of B('/Aj on j} is plotted in Fig. 2.1. A median value of |?> = 1 is suggested as the value of B above which there is significant induction into the earth. Using the ranges of S> and w in Table 2.1, and taking conductivities of .005 S/m (resistive crustal rock) and .5 S/m (conductive crustal rock) (Brace, 1971; Garland, 1979, pg. 277) the range of values of JS is calculated, and displayed in Tables 2.2 and 2.3. It is clear from the tables that the scale length has a sizeable effect on the value of 8. However, it must be noted from Fig. 2.1 that for all values of f> greater than 100, the magnitude and phase of |BJ/A' | will be approximately the 37 •4.0 0.0 4.0 LOGJp) 8.0 a.o 2.1 The ratio of the complex magnitudes of the inducing and induced magnetic field at the surface for the TE mode solution is B'/A'. Given is the dependence of this ratio on the dimensionless parameter p . a) MOD(B'/A') vs LOG,0 (fi) b) ARG (B'/A' ) vs LOGIO (p>) 38 same, so that for all ^ > 10*km. the surface values will be nearly the same. Also from Tables 2.2 and 2.3 it is clear that for the majority of the values of period, conductivity, and scaleiength there will be significant induction of the primary field into the earth for the TE mode. 39 Table 2.1: Parameter Ranges in G.D.S. 0" (crustal) .005 S/m - .5 S/m (separation parameter) 1.57x10 m - 1.57x10 m A (horizontal scale length) 400 km. - 40000km. T (period) 1s. - 72000s. (20hrs.) w = 2TT/T -s 8.7x10 Hz. - 6.28 Hz. / Table 2.2: Values of p for varying T, ?k . C = .005S/m T (km. ) 1 0s. 1m. 5m. 20m. 2hr. 1 Ohr. 4x1 04 4 .16x10 .27X10S 4-.53x10 .13x104 220 44 A-1x10 . 1 0x1 0S .17x104 .33xl03 83 1 4 2.8 5x 1 03 A-.25x10 .42x1 03 83 -21 3.5 .69 1x1 03 100 1 7 3.3 .83 . 14 .028 5x1 0Z 25 4.2 .83 .21 .035 .0069 T Table 2.3: Values of p for varying T,?v. 0" = .5S/m T (km. ) 1 Os. 1m. 5m. 20m. 2hr. 1 Ohr. 4x1 0* % .16x10 • 27X10-7 .53x10 .13x10 5 .22x1 0 4400 1x10 .10x101 .17x104 5 .33x10 8300 1 400 280 A 5x10* .25x10 .42x10* 8300 -21 00 350 69 1 x 1 0l . 1 0x1 0s" 1700 330 83 1 4 2.8 5x10* 2500 420 83 21 3.5 .69 T 42 Returning to the second type solution, the form of the electric field satisfying 2.1.8, 2.1.15, and 2.1.16 is found to be: Ez(x,y,z) - {^,^l^)Fz} A,e (2. 1.36) for z < 0 »3 (0X+<J>« for z > 0 (2.1.37) with Fz satisfying: (2.1.38) In similar fashion to the case for the first type solution, Fz(x,y) and are found to be the same for both the nonconductive and the conductive half spaces. As well, each value of S> between 0 and co merely represents one elementary solution with the complete solution again being the sum of all elementary solutions. Also in analogous fashion, the corresponding magnetic field for each elementary solution is found to be: 43 (2.1.39) for z < 0 Hx<xfy,z) - (i/w^)_^—{l&.F-^F0} Cxe (2.1 .40) for z > 0 Again using the boundary conditions at z = 0, we get: Bi - AX {[i -_|?_R]/[I +-|LR]} TO. Ya (2.1.41) CZ = AR {2/[l + AR]} (2. 1.42) where R is the same as defined previously, in equation 2.1.26, and: (2.1.43) To obtain the same format as for the type one solution, def ine: 44 (2.1.44) and (2.1 .45) Thus, the magnetic field at the surface becomes (using 2.1.39): H2(x,y,0) = {(k\ + B^)^.,-(A[ + Bl')^,0)} (2.1.46) with: - Al [(-^R - D/(-^-R + D] (2. 1 .47) When the amplitude of (A^ + Bj*) is evaluated over the range of values of c,w, and-0 given in Table 2.1, it is found that in all cases: lAi + Bj_' I = \C[\ V- 0 (2.1.48) The largest value of IC^I occurs when the period is 1 second, the conductivity is .005 S/m., and the horizontal wavelength is 40,000 km., and is still only 2.8 x 10 . Thus, for the range of parameters of G.D.S., both the transmitted field and the total 45 surface field are effectively zero, so that the TM mode will neither induce currents within the earth, nor be measureable at the surface. The TM mode can therefore always be ignored in G.D.S. The analogy between the first type solution and TE mode waves, and the second type solution and TM modes waves is completed in Appendix B. Also included is the derivation of the standard form of the Fresnel relations and Snell's law. It is found that the elementary solutions of the induction equation corresponding to different values of the separation constant , are in fact plane waves of differing total horizontal wavenumber. The main points of this analogy are given below. First type solution (TE Mode) E,(x,y,z,t) = {ky,-kx,0}A, e^e^^^e^1 + {ky,-kx,0}B, e^elVk'»lei-t (2. 1.49) H,(x,y,z,t) = (-1/wjUe) {( kza ) kx , ( kza ) ky ,-V1} I Aie^e J^e e 46 + d/w^e) {(kza)kx, (kza)ky, "v*1} B, eik*xetk^eik^* e iujt for z < 0 E,(x,y,z,t) = {ky,-kx,0} C, e**" e^V^6* eiu>± H,(x,y,z,t) = (1/wjuio) {(kze)kx,'(kze)ky,-^} C, e^VSV'^V"^ for z > 0 Here, A,, B,, and C, obey 2.1.24 and 2.1.25, Second type solution (TM Mode) E1(x,y,z,t) = (1/w€o) {(kza)kx,(kza)ky,- Vv} A, e * e J Je **• e (2.1.50) (2.1.51) (2.1.52) + (1/wO {(kza)kx, (kza)ky, ^v } B^e^V^V^e*"* (2.1.53) 47 H. 2(x,y,z,t) = {ky,-kx,0} A^^VjV^ e -{ky,-kx,0} B^V^e^* edl0t (2. 1.54) for z < 0 E\(x,y,z,t) = (l/w£6) { (kze)kx, ( kze) ky,->0V} (2.1.55) . H2U,y,z,t) = {ky,-kx,0} cie:k*Vk:i V:k* V"' (2.1.56) for z > 0 Here, At ,BX , and C1 obey 2.1.41 and 2.1.42. In these elementary solutions, the wavenumbers are shown (in Appendix B) to be related to the separation parameter by: Ckx* +ky'" = V (2. 1.57) kza = -i( WL = -i( 0r-w U0CT)"< (2.1.58) 48 kze = -i( % + ^)'K= -i( ^'L + iw/u6c)"r (2.1.59) The equivalence of the elementary waves of Chapter I (section 1.2) to the elementary solutions of the homogenous earth induction problem is thus complete. The separation constant is found to be related to the scale length discussed in the source field description: A = 2tr/v (2.1.60) Again, the complete solution will be a summation of all the elementary solutions, or in the continuous case, an integral over all horizontal wavenumbers. In this complete solution only the TE mode need be considered. 2.2 Uni form Source Field Over a Two-Pimensional Earth:  Forward Modelling A common simplification in G.D.S., and one used in most forward modelling algorithms, is that the source field is effectively horizontally uniform. Thus, the waves of the source field can be considered to be propogating vertically downward. The limits of this assumption are discussed in detail in Appendix C; in quick summary, the results indicate that for 49 periods less than 2 hrs. and for areas where the source field scale length is. > 5000 km. (such as at midlatitudes), this assumption will be valid. Jones and Pascoe (1971) have embodied this simplification in a finite difference approach to calculating the surface magnetic and electric fields over two dimensional structures, using the Gauss-Seidel iterative method. Their programs handle separately the situations where the electric vector of the primary field is parallel to the strike of the structure ( E-polarized) or the magnetic vector of the primary field is parallel to the strike ( H-polarized). However, in the case of H polarization, the symmetry along the strike of the structure ensures that the only non-zero component of the total magnetic field at the surface will be the horizontal component perpendicular to the strike. Further, for this polarization the continuity of any magnetic field vector across an interface guarantees that the values of this non-zero magnetic component will be the same at every surface position (Jones and Price, 1970). Thus, in G.D.S. the only interesting case is that of E polarization. Because of the expense involved in running the program of Jones and Pascoe (J-P), only a few selected conductivity structures have been studied. It is often considered that there are three major types of localized anomalies present in the earth, these being (Schmucker, 1970, pg. 78; Oldenburg, 1969, pg. 117): 50 1 ) near surface anomalies in the upper crust intermediate anomalies in the lower crust and uppermost 2) mantle 3) deep anomalies due to imbalances or undulations in the high conductivity layer of the upper mantle The three models shown in Fig. 2.2a,b,and c are representative of these three anomaly types. For each, the J-P forward programs have been used to generate the surface values due to an E-polarized inducing field. The "anomalous' field was then produced by subtracting from these the 'normal' field values due to the anomaly-free one dimensional structure. The results are given in Fig. 2.3a,b, and c for a period of 5 min., and in Fig. 2.5a,b, and c for a period of 50 min. If first order induction effects only are important, that is, mutual induction is insignificant, then the induced currents will be directly proportional to, and in phase with, the local electric field. Thus, for these localized conductive bodies used, the anomalous induced currents would resemble line currents in phase with the local electric field. From equation 2.1.52, the downgoing electric field in a half-space of conductivity can be written (considering only the TE mode): 51 E(x,y,z,t) = {ky,-kx,0} C(kx,ky) e'^e'^e'^V"'4 (2.2.1) where from equation 2.1.59: kze = -i (^,"+iw/Al>ff),'v (2.2.2) and from 2.1.57: S>V= kx1" + ky"*" (2.2.3) To allow for vertical .incidence, the vector components are normalized by division withV, giving: - i / \) ( StA\© • X + COS© • n \ E(xry,z,t) = {cos©,-sin©,0} C (©) e J -Ck4,4 Cult . e tc e (2.2.4) with: sin© = kx/S) ; cos© = ky/X> (2.2.5) Taking the strike of the two dimensional structure to be in the 'y' direction, then a vertically propogating (S) =0) E-polarized wave (sin© =0) of the TE mode (E% =' 0) will be represented by: 52 E(z,t) = {1,0,0} C' e"ik*c*elu,t (2.2.6) The portion of the spatial exponential which corresponds to propagation is obtained by taking the real part of kze from equation 2.2.2, so that the phase difference between the electric field at the surface and that at a depth Z0 will be given by: Ze (in radians) (2.2.7) If first order induction only is significant, the current along the anomaly and its corresponding magnetic field at the surface will be in phase with the electric field at the depth ZQ. The phase difference between the anomalous magnetic field and the electric field measured at the surface will thus be from equation 2.2.7. To test the approximation of first order induction the field of a line current at the same depth as the anomalous conductivity has been superimposed on the results from the forward induction, in Fig. 2.3a,b, and c, and Fig. 2.4a,b, and c. In all cases the magnitude of this current is normalized to the maximum value of the calculated horizontal field from the J-P induction program. For both the shallow and the intermediate anomalies, the match of the line current and J-P results is nearly exact for the horizontal field component. This should in turn mean that the vertical field match is as good, because the 53 SURFACE POSITION (KM-) 0 lOO 200 300 400 J—i i i i i SURFACE POSITION (KM) 0 100 200 300 400 o . in X-I Q_ . LU • o in . 1— 1 1 1 1 1 1 1 fb) SURFACE POSITION (KM1 0 100 200 300 400 J 1 1 I I L o . in 5T. X-h-CL UJ Oo in. 0 fl o b (<0 2.2 Three representative models used for the Jones-Pascoe (1971) forward induction program. The size of the anomalies is to scale. a) Shallow b) Midcrustal c) Mantle Undulation or 'bump' 54 100. 0 400,0 JOO.O Fig. 2.3 Results from the Jones-Pascoe program for the three models (a,b,c) of Fig. 2.2, at a period of 5 min. The solid line represents the J-P results, and —•©—- is the superimposed field of a line current at the depth of the anomalous conductivity. The magnitude of the line current is normalized to the maximum value of the Bx component. 55 Bx Bz 80,0 1 SO.0 200 0 320.0 100 0 3 0 SO 0 160.0 240.0 320.0 400.0 400.0 X X ig. 2.4 Results from the Jones-Pascoe program for the three models (a.,b,c) of Fig. 2.2, at a period of 50 min. The solid line represents the J-P results, and e is the superimposed field of a line current at the depth of the anomalous conductivity. The magnitude of the line current is normalized to the maximum value of the Bx component. 56 vertical field is always the negative Hilbert transform of the horizontal field for a two dimensional structure (see Appendix D). However, the vertical components are not nearly as similar as their horizontal counterparts. It is suggested that this indicates some minor failing in the J-P program, of either inherent nature, or due to the grid selection by the author. The matches for the third type anomaly are quite good for the horizontal component, but are very poor for the vertical component. In this case, however, the original assumption of a localized anomaly is not true, so that mutual induction is possible within the high conductivity zone beneath the undulation. Thus, the fit of the line current field is not expected to be as good in this case. To further check the validity of the line current approximation the phase difference <j> between the anomalous field and the surface electric field has been calculated from equation 2.2.7, and also directly from the J-P results (see Table 2.4). The values for the shallow and midcrustal anomalies are comparative in each example. For the mantle undulation there is a large difference between the two estimates of phase, and this is again indicative that in this case mutual induction is important. 57 Table 2.4: Phase Difference Between the Anomalous Surface Magnetic Field and the Normal Surface Electric Field, for Both the Line Current and the Jones-Pascoe Results. Anomaly Type Depth Period Phase eq.2.2.7 Phase J-P Shallow 9 km. 5 min. 1 .87° 1.74° Midcrustal 35.5 km. 5 min. 7.38° 5.3° Mantle 'Bump' 140 km. 5 min. 29.10° 21.0° Shallow 9 km. 50 min. .59° .46° Midcrustal 35.5 km. 50 min. 2.33° 2.44° Mantle 'Bump' 140 km. 50 min. 9.19° 35.93° In summary, for localized anomalies in a host material of comparative low conductivity, the J-P forward modelling results suggest that the line current approximation is in fact valid, indicating that the first order simplification for such anomalies is appropriate. 58 Chapter III Traditional Methods of G.D.S. Analysis 3.1 The Formulation and Separation of the Normal Field G.D.S. has been used for one dimensional as well as two and three dimensional sounding (Schmucker, 1970; Kuckes, 1973; Lilley, 1975; Woods, 1979), but the one dimensional problem will not be discussed here. The general presumption made in G.D.S. analysis, as originally suggested by Schmucker (1970), is that the earth model is basically plane layered, with only relatively small anomalous variations of conductivity (as indicated in Fig. 3.1a): Oj(x,y,z) = C^"+C*(x,y ,z) (3.1.1) The magnetic field within the j layer, Hj(x,y,z) can always be written: 59 Fig. 3.1 Possible Conductivity Models (a) The assumed conductivity model of G.D.S. is basically plane layered, with only small localized regions of anomalous conductivity. (b) A conductivity model that is plane layered except for a large discontinuous 'step'. (c) The abutment of two different layered structures (as at the land-sea boundary).. In both (b) and (c) the concept of a normal field breaks down if the array is close to, or spans the discontinuity. 60 Fig.3.1 AIR EARTH CO 00 AIR cr,' EARTH.-9}" 61 Hj(x,y,z) =Hj(x,y,z) + H^(x,y,z) (3.1.2) where Hj is the 'normal' field which would exist in the absence -i Ok. _i r of the anomalous conductivity, and Hj is the variation from Hj due to the anomaly. The general induction equation is equation A.19 from Appendix A: (3.1.3) Inserting equations 3.1.1 and 3.1.2 into 3.1.3, and using equation A.31 from Appendix A for the normal field in a region of constant conductivity, the anomalous field Hj can be expressed as a function of the normal field: 7V^ + -feT" x (Vx H>) - iw^OjH* = iv|-.<?HT " |^x ( § x Hp J (3.1.4) This equation, although never used directly, forms the basis for most of the methods that attempt to relate the designated normal field to the anomalous field in some statistical manner, and then from the ensuing relations estimate the anomalous structure. In cases where the extent of the anomalous variation is very large, the concept of a normal layered structure will not be applicable, and the defined 62 concept of a normal field will break down. Examples of this are when the measuring array spans a large subsurface conductivity step, or is in the vicinity of the abutment of two different layered structures (as in the case of the earth-ocean boundary), (see Fig. 3.1b,c) To avoid these problems of definition Lilley (1974) has suggested that the total field be formulated in terms of internal and external parts rather than anomalous and normal parts. Banks (1979) has also criticized Schmucker's approach on the grounds that the earth's crust is far too locally heterogenous for a normal field to ever truly exist. In view of the marked similarity of magnetograms in general over an array, this is probably too harsh an indictment; however, the existence of a true normal field should always be viewed with some reservat ions. If the concept of the measured field being composed of normal and anomalous components is accepted, its usefulness will then be dependent on the ability to separate these two components. Similarly, if Lilley's suggestion to use internal and external components is adopted, then it will be required that these two components be separable from the total measured field. As indicated in Tables 2.2 and 2.3 and Fig. 2.1, the induced vertical component over a conductive half-space almost completely cancels out the vertical component of the inducing field for nearly all values of period, conductivity, and scale length expected in G.D.S. Except for a combination of long period (> 2hr.), short scale length (< 5000 km.), and low 63 conductivity (< .005 S/m.), the vertical component of the normal field will be very small. Thus, in anomalous regions, the separation is often effectively already done for this component, as virtually all of the measured value will be anomalous, and internal. This forms the basis for the visual methods described in section 3.2 of this chapter, which use the unseparated data for preliminary analysis. To separate the normal and anomalous fields, two methods suggested by Schmucker (1970) are commonly employed. These methods differ only in their manner of separation of the horizontal components. In the first method a station presumed to be distant from the anomalous region (as would be indicated by the visual methods of section 3.2) is selected as the reference station. If the scale length of the observed field is much larger than the array size, then the horizontal components H and D of the field {H,D,Z} measured at this site will serve as the normal field. (Note that H is the component of the field in the direction of magnetic north, D is the component in the direction of magnetic east, and Z is the vertical component, where the positive direction of Z is downwards). In the case that this condition is not met, one can use the regional horizontal derivatives of the fields to calculate the first order corrections at any position (x,y) (where the position (0,0) is the reference site, 'x' is toward magnetic north, and 'y' is toward magnetic east): 64 HM(x,y) = H(0,0) + ("^).x + (^).y ^x j-(3.1.5) DN(x,y) = D(0,0) + (^L)-x + (I&»).y 3x *»j (3.1.6) From equations 2.1.20 and 2.1.36 in Chapter II, and the values of |B /A I given, the vertical electric field at the surface should always be nearly zero. Thus, the vertical component of the curl of the magnetic field should be near zero, allowing a check of the values of ^Ui/ax. and ^Hv/aij : liiii. - 22*. = iwecE, v o ^ ^ ^x ° * (3.1.7) The second method used to separate the normal and anomalous horizontal fields takes advantage of"the presumption that the scale length of the inducing field should be much larger than the spatial wavelengths of the induced anomalous field. Thus, spatial smoothing of the measured horizontal components over the array area should define the normal field. To define the normal vertical field, Schmucker uses the spatial derivatives of the normal horizontal field. From equations 2.1.19 and 2.2.23 of Chapter II it can be shown that at any frequency the value of the vertical magnetic field in the homogenous earth case is linearly related to the derivatives of 65 the horizontal components by a frequency dependent constant, C: 1 = c ( + ^M ) (3.1.8) where: C = l/(^+9e)'/z (3.1.9) In the plane layered case the relation is the same as in equation 3.1.8, but now C is a more complex, frequency dependent function of the conductivity structure (Schmucker, 1970, pg.15). (In both cases C is a measure of the skin depth of penetration). The value of C is calculated from some given model of the one dimensional conductivity structure. Using this, and estimates of the horizontal field gradients, it is then possible to evaluate the value of the vertical normal field, Z^. In practice, this method of determining Z^ is not very reliable for a number of reasons. The difficulties in accurately determining the one dimensional structure and the horizontal field gradients constitute the first obvious problem with this method. More fundamental, however, are the problems incurred because of the frequency dependence of C, and the time variations of the horizontal gradients. If the value of ZN is calculated in the frequency domain, then the time variations of the gradients cannot be taken into account, whereas in the time domain the frequency dependence of C cannot be introduced (Schmucker, 1970, 66 pg.16). To avoid these difficulties, it is this author's opinion that the normal vertical field ZN is just as reliably determined by assigning to it the value of ZN at the reference site. Whatever the method used to determine ZN, its spatial variation can safely be ignored. This follows from equation 3.1.8, as the gradients of Zw will now be second order corrections with respect to HN and Dw. Once a normal field {H^D^Z^} has been defined, the anomalous field {HA(x,y),DA(x,y),Zft(x,y)} at any station position can be calculated from the measured field, {H(x,y),D(x,y),Z(x,y)}: HA(x,y) = H(x,y) - Hw (3.1.10) Dfl(x,y)=D(x,y)-DNI 7 (3.1.11) ZA(x,y) = Z(x,y) - ZN (3.1.12) The problem of separating the internal and external fields was first solved by Gauss (1839) for a spherical earth. For the flat earth case,the separation method has been derived in a 67 variety of ways (Vestine, 1941; Siebert and Kertz, 1957; Weaver, 1963)). The derivation of the separation formulae for a two dimensional structure which is given in Appendix C is taken mainly from Weaver (1963). Consider a continuous one dimensional array running east-west perpendicular to the north-south strike of a two dimensional current density structure (see Fig. 3.2). By symmetry considerations the only non-zero components of the magnetic field will be D and Z. The relations between the internal and external components of D and Z (as derived in Appendix C) are: K(DX) = -zr (3.1 .13) K(DE) = ZE (3.1.14) K(ZX) = Dx (3.1.15) K(ZC) = -D£ (3.1.16) where the operator K is the Hilbert transform: 68 AIR Fig. 3.2 One dimensional surface array perpendicular to the strike of a two-dimensional earth. The current densities travel into (or out of ) the page. 69 oo r K(A(x)) = -1/1Y j A(u) /(x-u) du "* (3.1.17) with denoting the principal value of the integral. The measured data is: D(x) = Dr(x) + D£(x) Z(x) = Zx(x) + ZE(x) Thus, combining 3.1.13 - 3.1.19, we get Dx(x) = {D(x) + K[Z(x)]}/2 (3.1.18) (3.1.19) (3.1.20) Zj(x) = {Z(x) - K[D(x)]}/2 (3.1.21) The intrinsic shortcoming here lies in a fundamental property of the Hilbert transform. Because the denominator of the integrand in 3.1.17 is an odd function, the Hilbert transform of a constant is zero. Thus, for the normal field, this separation technique will never be able to separate the internal and external portions if the inducing field is uniform, 70 that is, if it consists of waves propagating vertically downward (corresponding to an infinite scale length). As seen in section 2.1 (arid as discussed by Weaver (1973)), this is an intrinsic property. In the case of a uniform field the ratio of the secondary to the primary field amplitude can have any value between 0 and 1, and there is no way to distinguish from measurements what the true value is. Thus, there will be problems separating the internal and external portions of the normal field if the scale length of the inducing field is greater than the dimensions of the array. It is clear from the above that estimation of the scale length from the surface readings will be important. In a layered structure, the crossing of a plane wave into the next layer will never result in a change of the original horizontal wavenumbers for the reflected or refracted wave (as seen for the earth-air interface in section 2.1, also, see Panofsky and Phillips (1962, pg. 196). Thus, the scale length should be calculable from either the separated normal field, or the separated external field. In the case of the homogenous earth, the inducing field ( Hp) consisting of only a single TE mode wave may be found by applying equation 1.2.4 to 1.2.9: Hp(x,y,0) = H0{-[kzkx],-[kzky],[kx1+ky' ] }et^,< e lk3^ (3.1.22) This would correspond to the separated external field. Using equation 2.1.29 the resultant total surface field (H) will be: 71 H(x,y,0) = H0{[-kzkx(l+T)],-[kzky(1+T)], [(kx1+kyl)d-T)]} elkx*e (3.1.23) where T is the ratio of induced to inducing complex amplitudes from equation 2.1.30. H corresponds to the separated normal field. The true' scale length is given by: A = 2TT/(kxl+ky V1 (3.1.24) Using the horizontal components of either the separated normal or external fields, an estimate of A is obtained from: > = ia(HZ+DVV[i(^L+^-) ] (3.1.25) where it is noted that H1, Da are H-H, D-D, and not |H|a,|D|*". The expression for the plane layered case is identical to equation 3.1.25, where the constant,C in this case is a complex function of the layered structure. In both cases C is a measure of the depth of penetration of the field. If the inducing field consists of many waves of differing spatial wavelengths, the non-linearity of equation 3.1.25 with respect to H and D will result in an incorrect value for A . However, it is expected that the value will still be a reasonable estimate of the smallest significant horizontal wavelength. It should be noted that as the complex form of the wave solution has been used in equation 3.1.22 and 3.1.23 that evaluation of X by this method would neccessarily be done in 72 the Fourier transform frequency domain. A similar estimator has been used which is applied to the time domain values of the horizontal normal or external field (Porath et al, 1971): V = 21YF/|VF| (3.1.26) F is the total horizontal field measured at the spatial position -a of the maximum gradient of F, and |v"F| is the magnitude of the maximum gradient. This estimator will be incorrect except under very fortuitous circumstances. Consider the expression for the complex field in equation 3.1.23. Because the estimator uses real time domain values, we use only the real part from this expression. Rotating into a new coordinate frame so that there is only one horizontal component F along ^, the value of F will be: F = FQ kHcos ( kH ) (3.1.27) where: Fc = -H^kz-Realf 1+T] (3.1.28) and kH is the total horizontal wavenumber. Using this to evaluate ~)\ in equation 3.1.26, we get: 73 = (2 TT/kH ) cot(k„ 1j ) (3.1.29) = A cot(kH^ ) (3.1.30) Thus, the estimate of the value will be incorrect by the factor cot(kH^), and as this can have any value between 0 and CO, the possible fluctuations of V from the true value X could be huge. It would be possible to use this estimator if the gradient and total field value were taken at positions a quarter scale length apart. However, as it is the scale length that is being determined this is not a very practical suggestion. For both of these methods the estimators are applied to either the separated normal or separated external fields, so that errors in the separations could result in errors in the scale length value. However, both methods of separation of the normal and anomalous fields will tend to error by not including enough small spatial wavelengths, whereas the separation of the internal and external fields is unable to separate the long spatial wavelengths. Thus, estimates of A made using both kinds of separated fields should provide bounds on the value of the scale length. 74 3.2 Visual Methods of Analysis As was mentioned in the previous section, because the induced vertical field over a normal layered earth always opposes the inducing vertical field, the anomalous amplitude will generally be much larger than the normal amplitude. The vertical component is therefore effectively already separated, by virtue of its being mainly anomalous and internal. This offers the rationale for doing initial preliminary analysis before separation, by using displays of the data itself. For example, consider the case of a buried linear feature with enhanced conductivity relative to the host rock. As seen in section 2.2 of Chapter II, the anomalous magnetic field at the surface will mimic the field of a line current. The vertical field will thus undergo a reversal in sign along a profile at right angles to the strike of the linear feature, with the zero field point directly over top of it. As seen in Fig. 3.3, this allows preliminary tracing of linear features; clearly a conductor runs roughly north-south between stations CHU and RAW, WIC and CUS, and REE and BAK, in the southern part of the array. Another possibility is that the time variations of the vertical component of the magnetic field will be closely correlated with one of the horizontal field components. This is seen in Fig. 3.3, in which , for example, the Z components of stations REE and WIC are closely correlated with their 75 Fig. 3.3 (a) Magnetometer locations for the magnetograms of (b). The array consists of 8 lines trending East - West, which are numbered from North to South. (b) Magnetograms for a substorm of August, 1972, from the southern part of the array (after Alabi et al, 1975). 76 associated Y components. The correlation indicates that it is the primary field of that horizontal component which is the dominant source of induction in the anomaly. Presume first order induction only, that is, the secondary effects of mutual induction by the induced field are ignored. The anomalous current will then be perpendicular to the correlated magnetic field component, as illustrated in Fig. 3.4. In accordance with the results of section 2.2 of Chapter II, the anomalous conductivity must then trend in this same direction. For the anomaly of Fig. 3.3, this trend will therefore be in the north-south direction, which is in agreement with the result of the Z reversals. An enhanced use of the unseparated data takes advantage of the feature noted in Chapter II in both sections 1 and 2, that the induced field due to a conductive region will be more closely in phase with the inducing fields the higher its conductivity is. Thus, by taking temporal Fourier transforms of the vertical component at every station in a two dimensional array, and contouring separately the amplitude and phase results, one will clearly see the paths of large highly conductive anomalies by the clustering of the contours in the phase diagram. In Fig.3.5, the path of the North American Central Plains anomaly through Saskatchewan and into the U.S.A. is clearly delineated in the contour plot of the phase of the Z component (Alabi et al, 1975). 77 Fig. 3.4 A plan view of an anomalous current in the earth is shown. Presuming only first order induction, the induced current in a localized conductor will be perpendicular to the magnetic vector of the inducing field. Fig. 3.5 The contoured amplitude and phase of the magnetograms of Fig. 3.3b, at a period of 68.3 min. 78 3.3 The Induction Tensor and Induction Arrows Presume that the normal field {HN,DN,ZN} has been defined, so that the anomalous field {H^,D^,Zft} can be separated from the total field {H,D,Z}. The methods of analysis using induction arrows are then based on the assumption that for any frequency w, the anomalous field components are related to the normal field components by the induction tensor, I: The elements of I may be complex, to allow for possible phase differences between the anomalous and normal fields. Obviously this relation can always be defined at each time t, since there are only three equations for nine unknowns. Unfortunately, this admits the possibility that any scheme used to compute the tensor elements will result in the values being (3.3.1) where: I = (3.3.2) 79 time dependent. If the induction tensor elements can fluctuate with time this indicates the tensor is a function of the inducing field as well as the underlying conductivity. The ensuing problem of trying to separate the two influences in the tensor will severely limit its usefulness in resolving the earth's structure. Thus, the tensor must be reasonably independent of the inducing field, and thus of time, to be of value. If the source field at each frequency w was comprised of only one elementary electromagnetic • wave, characterized by wavenumbers (kx,ky,kz), then obviously the linearity of equation 3.1.4 ensures that the values of the induction tensor will be constant,even with time varying amplitudes. For a single wave, and a given conductivity structure, the relation of 3.3.1 would be: values of the normal field components for an inducing wave of unit amplitude at wavenumbers (kx,ky). A(kx,ky) is the actual time varying amplitude of the inducing wave, and I'(kx,ky) is the induction tensor for the wavenumbers (kx,ky). The values of the elements in I'(kx,ky) will depend on the interaction of the / i' (kx,ky)A(kx,ky,t) / Hw(kx,ky) D^(kx,ky) V Z*N(kx,ky) (3.3.3) The vector components {HN(kx,ky),DN(kx,ky),ZN(kx,ky)} are the 80 transmitted wave with the anomaly, and so will be dependent on the spatial characteristics of the wave as given by (kx,ky). Thus, I'(kx,ky) will be different for each (kx,ky). Because the normal field of a wave at a plane boundary (incident plus reflected) will always have the same magnitude as the transmitted wave, the actual magnitude of the inducing wave does not enter into I'(kx,ky). In the case of more than one wave in the source field the total anomalous field will be the integral over all possible wavenumbers: CO CD where now A(kx,ky) is an amplitude density in wavenumber space. The total normal field is given by: Thus, to be able to represent the total anomalous field as given in equation 3.3.4, by an expression of the form of equation 3.3.1 in which I is source ( and thus time), invariant requires that: (3.3.4) ( (3.3.5) 81 00 co -oc -co = j^[l-I *(kx,ky) ] A(kx,ky,t) ( Hw(kx,kyAdkxdky )w(kx,ky) J !N(kx,ky)/ (3.3.6) Because the values of A(kx,ky,t) are arbitrary, this can only be satisfied if for all (kx,ky): I = I1(kx,ky) (3.3.7) As was pointed out earlier, this is not in general true. Thus there is nothing that requires a priori that the induction tensor I relating the normal and anomalous fields in 3.3.1 will be source independent. A variation of equation 3.3.1 was suggested by Dragert (1973, pg. 41 - 42). Taking the spatial Fourier transforms of the array measurements of both the normal and anomalous fields at a particular frequency, w, one would obtain their respective complex amplitudes at each possible pair of wavenumbers kx,ky. The assertion was then that the relation of 3.3.1 would be true when considered independently at each wavenumber pair (meaning that I, the induction tensor, would be a function of kx,ky). The suggestion is not correct, because a single elementary wave of arbitrary wavenumbers kxo,kyo, can generate an entire horizontal spectrum of waves in the anomalous field. Consider the example, treated in section. 2.2 of Chapter II, of a single downgoing wave impinging on an earth that is homogenous except for a 82 buried cylinder of high relative conductivity. As was seen, the anomalous field of this example to first order resembled that of a buried line current. It was shown in section 1.2 of Chapter I that a line current has energy at an infinite number of wavenumber values, so that the spatial transform of the anomalous field would be non-zero at an infinite number of wavenumber values. On the other hand, as was shown in section 2.1 of Chapter II, the normal field due to a single impinging wave will only have one non-zero wavenumber component, at the wavenumbers, kxo,kyo, of the original inducing wave. Thus, to relate the anomalous field components at each wavenumber pair kx,ky would be impossible at all values for which the normal field value was zero and the anomalous field was not. It has been shown then, that the relation of equation 3.3.1 is not universally true. However, it has been proposed that for arrays at latitudes mid-way between the polar areas and the equator, the source field will be effectively horizontally uniform for the ranges of parameters used in G.D.S. (Lilley and Bennett, 1973; Banks, 1973), with estimates of the scale length being on the order of 5000 to 10,000 km. (Gough, 1973; Banks, 1973; Lilley, 1973; Mareschal, 1981). Under these conditions it is claimed that the induction tensor will be independent of the source variations. Consider the effect on the induction tensor formulation in the uniform field situation. In the limit as the kx and ky values of a wave go to zero, the direction and mode (TM or TE) of the wave no longer serve to distinguish the 83 direction of the electric and magnetic fields in the horizontal plane. For a vertically downgoing wave there is no distinction between the TE and TM modes, and the only wavenumber is kz, yet there are an infinite number of possible rotations of the wave components about the direction normal. It is found, however, that any vertical wave can be decomposed into a sum of two waves: one with its electric vector along some given horizontal direction (E polarized), the other with its magnetic vector along this same direction (H polarized) (Jones, 1971). Thus, the uniform field case still involves two wave types with the consequent relation between the anomalous and normal field now being: (3.3.8) where the subscripts E and M refer to E polarized and H polarized respectively. Equation 3.3.8 has been simplified because of the zero values of Dwe and HNh due to our choice of reference axes, and because of the zero value of the normal vertical field in the uniform field case. Also from this , we have the total normal field being given by {HNE'DNM'0}- Tnus, we can rewrite 3,3.8 as: 84 / Hft Here, the superscripts E and M indicate the original tensor that the elements are from. Thus, in the case of a perfectly horizontally uniform field the induction tensor will be independent of the source field, and should thus contain retrievable information about the underlying conductivity structure. What needs to be determined then, is whether a sufficient degree of uniformity for induction tensor invariance is attained at scale lengths in the range 5000 - 10,000 km. This is considered in detail in Appendix C. The discussion indicates that for the scale lengths suggested for midlatitudes •, and for periods greater than 2 hrs., the.values of I'(kx,ky) will in fact be nearly identical. Thus, the relation of equation 3.3.1 will be valid for these parameter ranges. It should be noted that it is because the tensor relates the normal field to the anomalous field that the magnitude of the transmitted wave need not be considered. This is not true of the similar relation between the external and internal fields. If these portions of the surface field are related by a tensor in identical fashion to that in equation 3.3.1 (as done by Lilley (1974)), the significant changes in the ratio of the transmitted wave to the external field at all values of (kx,ky) will always result in a 85 tensor that is time dependent. Consider that the induction tensor is in fact independent of time. One can then define the relation between the normal and anomalous field at each frequency w as in equation 3.3.1, except for an uncorrelated 'noise' term: (3.3.10) The estimated values of the elements of I will be those which minimize the power of the uncorrelated terms (Schmucker, 1970, pg. 20 - 21). By minimizing the power of each component of the noise vector with respect to the real and imaginary parts of each of the tensor elements, three independent sets of linear equations are obtained, one for each'column of I (see Appendix D for the complete method). As an example, the set of equations for the third column is: S%«-*KJ J (3.3.11) where: 86 S = / SHwHw S^jHN S?MHM f (3.3.12) is the matrix of auto and cross powers. The power of two signals A,B of length To is defined as: SflB(w) = A(w)B*(w)/ To (3.3.13) The most important terms of the induction tensor are those of the bottom row, corresponding to the relation of to the normal field: (3.3.14) Presuming that the contribution due to induction by the normal vertical component will in general be very small, we have: (3.3.15) For each frequency w, two induction arrows can be defined, corresponding to the portion of Z^ which is in-phase with and 87 U(w) = -{x-Real[CjH(w) ] + y-Real[C^(w) ]} (3.3.16) and the portion which is out of phase with HN and DN: V(w) = {x•Imag[CfrH(w)] + y•Imag[C*&(w)]} (3.3.17) where x is the unit vector in the direction of magnetic north and y is the unit vector in the direction of magnetic east. The negative sign on the left hand side of equation 3.3.16 is commonly introduced to maintain • the convention set by Parkinson's original development of induction arrows (Parkinson, 1959, 1962). Also, in the case of a localized conductive path, if first order induction only is significant, the sign convention will ensure that the in-phase induction arrow U will point towards this conductive anomaly. This is illustrated in Fig. 3.6 for a horizontal line current running north-south at depth zc. It is generally presumed that the in-phase induction arrow at a frequency w will point towards nearby current concentrations induced at that frequency , as seemingly indicated in Fig. 3.6. (Banks, 1973; Beamish, 1977; Garland, 1979, pg. 270). However, it is not clear that this can be accepted in general. Consider the horizontal conductive anomaly shown in Fig. 3.7. This anomaly is not offered as a possible real earth model, but simply as a test of the generality of the induction arrow's properties. Using the first order induction 88 H Polarized Hn za< 0 C£1> < 0 E ^Polarized Ed ~*"6 n C2B» ClH Z^>0 c« >o -« • A X Fig. 3.6 The direction of the in-phase induction arrow near a line current running North - South. 89 'n c, >0 >0 J&= rC|0| Drj JH=rC2aaHn C,<0 C2<0 a 2 — •.V X Fig. 3.7 First order induction in crossing conductive paths. The ratios of currents in the two conductors will depend on both the conductivity and the magnitude of each component of the inducing field. However, the coefficients of the induction tensor will depend only on the conductivity. 90 approximation throughout (that is, ignoring self-induction), we will have the H component of the normal field inducing current in the east-west direction and the D component of the normal field inducing current in the north-south direction. The final relation for ZA in terms of HW and D^, at any surface position (x,y) is: Zft(x,y,w) = r (w) • [C2(x,y) CTZHn + C((x,y)G, DW] (3.3.18) C, and Cj. are spatially dependent constants, with their sign as indicated in Fig. 3.7, and r(w) is a frequency dependent complex constant relating the normal components H^D^ to their associated electric fields at the depth of the anomalous conductivities. Thus, the in-phase induction arrow in this example is: U = -{[r(w)<rlc1(xfy)]x + [r (w)tr, C, (x,y) ]y} (3.3.19) The direction of U will depend on the position of the station through the values of C, and . More significantly, the direction of U will depend on the comparative values of cr, and C"2, rather than on the relative strengths of the currents in each conductor. If the station was situated such that |C,| and \CZ\ were equal (which would be true along the bisectors of the axes), then for <"i>>C"i, U would point towards the conductor lying along the y axis, and for <^i»CT2f U would point towards 91 the conductor lying along the x axis. The relative amplitudes of the current densities in these cases could take on any value, depending on the amplitudes of Hw and . This clearly illustrates that the induction arrows do not necessarily point towards current concentrations, but rather, from the examples considered, they point generally towards the area of high comparative conductivity. Up to this point, only the case where local induction effects were dominant in the production of the anomalous field has been considered, so that there was a neccessary correlation between the locally measured normal and anomalous fields. Another possibility for anomalous behavior is current concentration, or channelling (Whitham and Andersen, 1965; Dyck and Garland, 1969; Gough, 1973), in which a uniform current flow due to large scale regional induction is channelled through a localized high conductivity zone (see Fig. 3.8). The anomalous magnetic field in the channelled current case does not depend on local induction, thus, it is not immediately clear how the calculated values of the induction tensor, and therefore, the induction arrow, will relate to the position and orientation of the perturbing conductive body. In ignoring induction effects, we reduce the problem to one of quasi-static direct current (DC) flow. The pertinent Maxwell's equations, after elimination of the time derivative terms, are (from equations A.12, A.13, A.14 and A.27 in Appendix A): 92 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i i i i i i i i i i /ffii > « » \ u~ K 1 ' 1 1 ::: :f Y::: ,, ,^ ;:; II/// i i i / i i i i i i i i i i i 1111 1 1 1 1 1 1 1 1 1 1 Cb) Fig. 3.8 The effects of current channelling through a conductive sphere. The arrows give the magnitude and direction of the electric field outside the sphere, for the two orthogonal directions of the regional current flow: (a) regional current flows East - West (b) regional current flows North - South 93 V-E = 0 (3.3.20) V-H = 0 (3.3.21 ) VxE = 0 (3.3.22) VxH = (TE = J (3.3.23) The curl-free nature of the electric field enables us to write the electric field as the gradient of a scalar potential: E = -VV (3.3.24) which when coupled with equation 3.3.20 gives\ us Laplace's equat ion: (3.3.25) In this quasi-static case, we also have from equation A.20 in Appendix A the boundary condition ensuring the continuity of the 94 normal component of current flow across a boundary: V-(flTE) = V- J = 0 (3.3.26) Equations 3.3.24, 3.3.25 and 3.3.26 govern the DC problem. In the following discussion it will be presumed for simplicity that the source field is perfectly horizontally uniform. All the arguments pertaining to the validity of this assumption in practice and its effect on the invariance of the induction tensor will apply here as well, but will not be discussed. We will also presume first order induction effects only, so that the north-south currents will be proportional to D^ only, and the east-west currents will be proportional to only. The linearity of the equations involved allows us to consider separately the cases in which the regional current is only north-south, or only east-west. Any other direction of regional current flow can always be considered as a superposition of these, with the resultant current vector at any point being the sum of the vectors corresponding to each of the two solutions. In Fig. 3.8, a sphere of radius 'a' and conductivity imbedded in a host rock of lower conductivity <3", serves as the model. Equations 3.3.25 and 3.3.26 can be solved in spherical coordinates, for the boundary condition of a uniform electric -a field E0 at large distances away from the sphere. The symmetry of the problem ensures that the solution for one direction of 95 this regional electric field can be simply rotated by 1"T/2 to give the solution for the orthogonal direction. As an example, the solution for the potential outside the sphere for the uniform field in the 'x' direction is (Telford et al, 1976, pg. 647-649): V(x,y) = -E0 x{l-[a3(b-D]/[(b+2)(xX+yl)Vl]} (3.3.27) where 'b' is the conductivity ratio 0"2/cT, . From equation 3.3.35, the directional derivatives of the • potential give us the components of the electric vector at any point outside the sphere. For the uniform field in the 'x' direction: Ex(x,y) = EQ{1 + [ai(b-l)(2xl-y2)]/[(b+2)(x1+yIf/l]} (3.3.28) E,j(x,y) = E0{[aS(b-1)(3xy)]/[(b+2)(xl+yl)S/i]} (3.3.29) The electric field vectors for each of the two uniform field directions are shown in Fig. 3.8. It is clear that the current concentration effects on the outside of the sphere are most pronounced in the regions along the axis of symmetry parallel to the direction of the regional current flow (regions 96 1 and 3), and are at a minimum in the regions along the radial line orthogonal to this (regions 2 and 4). Thus, in regions 1 and 3 the anomalous vertical .magnetic field will be most dependent on D^fin regions 2 and 4 it will be most dependent on HN, and in the intermediate regions there will be a continuous grading, with the TT/4 bisectors marking the lines of equal dependences. Thus, from this qualititave discussion, the in-phase induction arrows will always point generally towards the sphere. In fact, symmetry will require that the arrows will always point towards the very center. In quick summary of the main points regarding the induction tensor and the induction arrow we have: (1) Relating the anomalous magnetic field vector to the normal field vector via the induction tensor is not universally correct. (2) In the case of a horizontally uniform field the induction tensor will correctly relate the anomalous and normal fields, to first order. The induction arrows will point in general towards regions of anomalously high conductivity, which will not neccessarily be the regions of highest current concentration. The degree of uniformity required for the invariance of the tensor will be attained at values of the scale length encountered in practice at midlatitudes, for periods greater than ~ 2 hrs. (3) The uncertainty of the invariance of the induction tensor 97 will apply to the case of current channelling as well. However, if the field is uniform,such as at midlatitudes, the induction arrows will point also toward the anomalous high conductivity zone. 3.4- Quantitative Methods Used in G.D.S. Quantitative modelling of the earth's conductivity for areas containing presumed anomalous structure has received much less emphasis in the 'past than the corresponding qualitative modelling. The true inverse problem of finding a model that fits the data in either a 'one-shot', or an iterative prodedure, has not been solved to this author's knowledge. The other approach to this is to select a possible model, either based on qualitative information and intuition, or simply picked randomly, and then use forward modelling routines to see if it fits the data. If model selection is done interactively with the forward modelling then the model can be adjusted until a good fit is obtained (Porath et al, 1970; Gough, 1973; ; Dragert, 1973, pg.94-96). However, the final model is always non-unique, so that it is very possible that the biases of different individuals doing the model perturbations could result in very different final solutions. As well, it has been found in some cases that a model to fit the data cannot be found (Whitham and 98 Andersen, 1965; Porath et al, 1971). (The suggested reason for this was that the anomalous field was due to current channelling, so that the areas considered in the model studies were simply not large enough to encompass the true model.) The alternative to the interactive approach would be the 'Monte Carlo' method, in which models are selected randomly for subsequent trial in the forward modelling routine (Cochrane and Hyndman, 1970). For the two dimensional and three dimensional cases this would clearly be prohibitively expensive. As well, it is not clear that for arbitrarily • complex models that the existing forward modelling routines will give accurate results. Thus, either the models must be made simpler (meaning they no longer are completely arbitrary) or the forward model results must be considered questionable. To avoid these difficulties, and yet still obtain estimates of some of the basic quantitative values, such as the depth and lateral extent of the anomalous conductivities, certain approximate measures have been devised. If the visual data displays and induction arrows consistently indicate a long narrow anomaly, then it is common to model the internal source as a line current. Presume that the current, of magnitude I, is parallel to the earth's surface and in the negative 'y' direction, and that the measuring array is perpendicular to the strike of the current (see Fig. 3.9). Then, by Ampere's law, the magnitudes of the directional components Bx,Bz at any position 'x' along the array will be given by: 99 3.9 A variety of different two-dimensional current density models, with the currents travelling into or out of the page. In all cases the station array is at the surface, perpendicular to the strike of the currents. The different models are indicated by: (a) <> (b) LU (c) X (d) X 3.10 The dependence of the 'peak to -peak' width of the Bz component, for the current density models of Fig. 3.9. The width's are seen to vary with the depth of the anomaly as well as with its lateral extent. 100 Fig. 3.9 Fig.3.10 87.5 100.0 Station Position 101 Bx(x) = ji0Iz0/{2fT [ (x-x0 )* +z*]} (3.4.1) Bz(x) = ju0I(x-xo)/{2fT [(x-x0)+z*]} (3.4.2) where (x0,z0) is the position of the line current (see Fig. 3.10 for sample plots of Bx(x) and Bz(x)). The depth of the line current can be found from either the distance between the point at which the Bx field is a maximum, '(xo,0), and the points at which it is half the maximum value (xi/,0): (3.4.3) or from the separation between the positions (x.^,0), (x_,0) of the positive and negative peaks of the Bz component: zo = 1/2|X_-XJ (3.4.4) It is sometimes claimed that this depth represents the maximum depth possible to the anomalous current, and that any other current distribution giving the same data values on the surface must be shallower. This statement is not strictly correct, as is shown in the following argument. Consider a distribution of currents along a sheet parallel 102 to the earth's surface, at a depth Z|. If the currents in this sheet are constrained to flow only in the positive or negative 'y' direction, with their magnitudes not varying with 'y', then a lineal current density j(x) is sufficient to describe them. It will be shown that for z( < z0 it is always possible to find a j(x) such that at the surface the field of the line current at (x0,z0) is completely duplicated. Without loss of generality we can simplify the algebra of the problem by setting the • coordinate origin on the surface directly above the line current. Thus, the field components of the original line current reduce to: Bx(x) = u0Iz0/[21T (z0+xj) ] (3.4.5) Bz(x) = p0Ix/[2TT (Zo+xJ-) ] (3.4.6) i i The surface fields, Bx , Bz , due to the sheet distribution of currents are obtained using the Biot-Savart law, and are: CO Bx' (x) = (/L»o/2*nr) ^{[j(u)z, ]/[(x-u)l+zl1]}du ~°° (3.4.7) 103 co Bz\x) = (^o/2TT )f {[j(u)(x-u)]/[<x-u)l+zj"]}du (3.4.8) Thus, for the fields of the sheet distribution of current to mimic that of the original line current requires that: oo y [ j(u)z]/[ (x-u)2+z^]du = Iz0/(z*+xl) ~°° (3.4.9) and CO ^ [ j(u)(x-u)]/[(x-u)l+z|l]du = lx/(z*+xr) ~°° (3.4.10) Taking the spatial Fourier transform of both sides of equation 3.4.9 and 3.4.10 (where the Fourier transform, 3~ , is defined as in equation C.5 in Appendix C), noting that the left hand side in both cases is a convolution, and then applying the Faltung (or convolution) theorem, we arrive at: J[j(x)]. J [z,/(x*+z.1)] = I [z0/ (z0l+xv) ] (3.4.11) 3 tjU)]-} [x/(x1+z,1')] = lJ[x/(zflV)] (3.4.12) Using the general Fourier transform expressions (from Appendix C, equations C.3 and C.4): 104 }[a/(x2+a1)] = Tre"1!1" (3.4.13) J [x/U^+a1 )] = -Tfie'^sgn^) (3.4.14) both equations 3.4.11 and 3.4.12 result in: (3.4,15) where: (3.4.16) Taking the inverse Fourier transform of equation 3.4.15 we arrive at the final result: j(x) = I<$V[TT( J^+x1)] (3.4.17) Thus, we can duplicate the field due to any line current with a distribution of current on a sheet situated between the line current and the surface. Correspondingly, this suggests that if such a sheet current was present at the same time as the line current, but with its direction of flow reversed, then the field at at the surface would be null. Thus, it is conceivable 105 that current distributions could exist of any magnitude and depth, but which in total give rise to no surface readings. If such an 'annihilator* distribution was present, then currents could conceivably be flowing undetected at greater depths than that of the original line current. It is also shown in Appendix F that even when all the currents in any model are restricted to being in phase, (so that 'annihilator' distributions are impossible) there can still be currents in the model that are deeper than the line current, with the model duplicating the surface data of the line current. The correct statement about the possible depth of the anomalous currents is that if the surface field readings can be duplicated by a line current at a certain depth, then there must be some current at that depth or shallower. This statement is proved in Appendix F. Another method of estimating the depth of an anomaly is by considering the smallest period at which the anomalous field is still significant, and then using this frequency in a skin depth calculation. One problem involved with this method is that of estimating the one dimensional conductivity structure for use in the skin depth calculation. The other problem arises because of the possibility of channelled currents. It is very possible that the channelling of the regional current system will not be strictly horizontal, so that the channelled current under the array will be at a different depth than that at which the regional system was induced. The dominant periods of the anomalous ,fields in such a case would thus be indicating the 106 depth of the regional current system, and not that of the channelled currents. Estimates of the lateral width of the anomaly might be made from the spatial extent of the anomalous field, as indicated on the contour maps of the vertical field's phase, or on contour maps of the ampitude of the horizontal anomalous field. However, the calculated fields (see Fig. 3.10) for the different current density models of Fig. 3.9 illustrate that the indicated width is a function of depth, as well as the actual lateral extent. Thus, the true width may be difficult to estimate in this manner, particularly for deep anomalies. A method used to calculate the extent of anomalous bodies has been to select a starting model that qualitatively fits the data, and then perturb the size of the anomaly until a best fit with the data is attained (Porath et al, 1970). As discussed previously, this type of interactive modelling is open to problems of human bias, as well as being expensive to perform. In summary of this section, it is clear that the intractability of the induction problem for two and three dimensional models severly limits the extent of quantitative analysis possible. The methods and estimators used have been shown to be of questionable validity. 1 07 Chapter IV The Current Density Model 4. 1 The Current Density Formulat ion Determining the earth's conductivity as a function of spatial position is the- ultimate goal in G.D.S. analysis. However, for reasons to be outlined shortly, a more useful approach in many cases is to consider the current density as the model parameter to be determined. In an initial simplification, only two dimensional models will be allowed, with all current flowing only in the 'y' direction. Although this two dimensional constraint will obviously not be applicable in all real examples, a sufficiently large number of elongated anomalies in G.D.S. do appear to have this degree of symmetry, making it a reasonable initial simplification (Gough, 1973). The symmetry assumptions ensure that a linear array perpendicular to the strike of the structure is sufficient to measure all non-redundant information available at the earth's surface (see Fig. 4.1). It will be presumed that the separation of the anomalous and normal fields measured along this array has been performed, 108 AIR Stations-v^ —•—x X *"* s\ /\ 7* X X X EARTH Densities Fig. 4.1 The presumed two-dimensional model. The conductivity is invariant in the Y direction (out of the page), and all currents are assumed to flow in this direction. A linear array of magnetometers perpendicular to the strike of th« model is sufficient to record all non-redundant information available. 109 so that henceforth all data will be considered to be purely anomalous and internal. The general expression used to calculate the magnetic field due to the current in a given volume V is the Biot-Savart law (Panofsky and Phillips, 1962, pg. 125): where r is the difference vector between the position of the (4.1.1) line^ current, (x*,y',z') and the position at which B is measured, (x,y,z ) : r = x, (x-x1 ) + X;, (y-y( ) + x^z-z1 ) (4.1.2) The components of B(x,y,z) are thus: Bx(x,y,z) = uD/4Try[[{[ (z-z1) J3-(y-y')j?]/[ (x-x1 )^ V + (y-y' )"" + (z-z' )x 3 dx'dy'dz (4.1.3) By(x,y,z) = ^o/4rr[[f {[(x-x1 )j4-(z-z' ) jx ]/[ (x-x ' V 110 +(y-y ) +(z-z ) J dx dy dz (4.1.4) Bz(x,y,z) = ^o/4tv]'([{[(y-y,)jx-(x-x, ) j,, ]/[ (x-x' )Z V + (y-y' )z+(z-z' J7" ] dx'dy' dz' (4.1.5) For the simplified current density system considered, j will be in the 'y' direction only, and its amplitude will not vary in that direction. The values of B(x,y,z) will be taken only at z = 0. Using these simplifications, and integrating from y =-cr> to y =+co , the components of B(x,y,0) become: Bx(x,y,0) = /io/2TT^{[-z'j^]/[ (x-x1 f+z,Z-]}dx'dz' S (4.1.6) By(x,y,0)=0 (4.1.7) Bz(x,y,0) = /V2tt^{[(x'-x)j l/IU-x^+z^Hdx'dz' S (4.1.8) where S is the area of interest outside of which will be presumed to be zero. In future, as Bx(x,y,0) and Bz(x,y,0) do 111 not depend on y or z, they will simply be referred to as Bx(x) and Bz(x). The major advantage of inverting the data to find a current density model rather than a conductivity model is the avoidance of the intractability of the non-linear induction equation. In the new formulation the model parameter j(x,z) is linearly related to the data as shown in equations 4.1.6 and 4.1.8. The forward problem expressed in these equations is easily and cheaply solved,as is shown in section 4.3 of this chapter. As well, the solution of the general linear inverse problem is very well understood, with a variety of methods and techniques available in the literature (Backus and Gilbert, 1967,1968; Wiggins, 1972; Oldenburg, 1976; Parker, 1977). Another advantage to using a current density model is that its very formulation avoids the problem of current channelling. The region of interest of the model can be defined precisely by using the geometric decay of the magnetic field away from a line current to determine the boundaries beyond which currents will have negligible effect at the measuring stations. The problem of having to solve over a large and poorly determined region of interest to ensure that the effects of channelling of the regional current flow are included (Porath et al, 1971; Gough, 1973) will simply not be encountered. The inherent disadvantage in using the current density model is its non-uniqueness. It has been pointed out earlier in section 3.4 of Chapter III that the field due to a line current 1 12 can be duplicated by an infinite number of different current distributions. As well, it is possible for 'annihilator distributions' to exist which will have non-zero current density values, but which will not give rise to any measurable surface magnetic field. As formulated, the current density inverse problem can only utilize the internally origined portion of the measured magnetic field, effectively ignoring the external field required to induce the internal current system. The disregarding of this information explains the difference in uniqueness difficulties between the conductivity and current density model approaches. The current density model has been used before (Banks, 1979; Woods, 1979), but the uniqueness problem has led the authors to constrain the current models to an infinitesmally thin sheet at some constant depth. The thin sheet formulation will probably give a good indication of the horizontal position of the actual earth currents, but it cannot give any information as to their depths. As well, the depth of the thin sheet must be chosen with some caution, for as shown in Appendix F, if the thin sheet is deeper than the actual current, then no current distribution along it will be possible that fits the data. The current density inverse problem left simply as originally stated is probably hopelessly non-unique. This will be illustrated in section 4.2 by Backus-Gilbert type appraisal, and in section 4.3 by the construction of a variety of dissimilar models, all unlike the original model, but all of 113 which match the data. However, by imposing certain expected physical features of the earth model as constraints in the construction process, it is possible to greatly restrict the range of allowable models. It will be shown in section 4.4 that 'true' models which match these physical requirements are closely recovered in this constrained construction. 4.2 Unigueness and Backus-Gilbert Appraisal The statement that a particular inverse problem is non-unique is rather nebulous, as there are varying 'types' and 'degrees' of non-uniqueness. In certain problems the uniqueness of the data set produced by a particular model can be proved, such as in the global G.D.S. problem (Bailey, 1970) or in horizontal loop electromagnetic sounding over a stratified earth (Fullagar, 1981). However, in both cases, to obtain the unique model requires an infinite amount of perfectly accurate data. As this condition can never be reached in practice, there will always be a certain range of acceptable values for the model, even though theoretically there is a one to one correspondence between models and data sets. The non-uniqueness incurred in the above-mentioned cases would be expected to be 'less' than that in cases such as the current density problem, in which an infinite amount of perfect data would still not guarantee a 1 1 4 unique model upon inversion. As well, even though the total model may not be unique, there may be portions of it, or features in it, which are required in all models that fit the data. Backus and Gilbert (1967,1968,1970) have outlined the appraisal procedure to quantify these different aspects of non-uniqueness. In a linear problem, the model and data can always be related by a Fredholm integral equation of the first kind (Parker,(1977): where: e- : is the i datum g- : is the kernel function associated with the itk datum m : is the model R : is the region of interest outside which the model is considered to be zero (in this case R is a surface) Consider that there are N data, with N associated linearly independent kernel functions. A function, designated an averaging function, is now constructed from a linear combination of the kernel functions: (4.2.1) 115 N A(r,T0) = £4i(r0 )gL(r) (4.2.2) such that the function A is unimodular (ie. it has unit area) and is as 'close' as possible to a Dirac delta function at some designated position r0. The measure of closeness to a delta function can be defined in a number of different ways (Backus and Gilbert, 1970; Oldenburg, 1976); the measure used here is the first Dirichlet criterion (Oldenburg, 1976): S(r0) = £$[A(r,r0) - <f( r-tQ ) ]*ds R (4.2.3) The average of a model at the position rQ is defined as: <m(r0)> =Jrm(r)A(r,r^)ds R (4.2.4) If the averaging function A(r,rc) was exactly a Dirac delta function, then: <m(r0)> = $}m(r)£(r,r0 )ds = m(r0 ) (4.2.5) —k so that the value of the model at the position rc would have been exactly recovered. In any event, independent of the form of A(r,r^,), the values of the averages <m(r0)> are unique to the problem, as they are dependent only on the values of the coefficients and on the data: 1 1 6 <m(r^)> = f^m(r )A(r (r6 )ds = $£m(r)£o<. (r0)g- (r)ds M = ^ ^ (re )[f[m(r")g-(r)ds] = (ro)e-(4.2.6) Thus, the averages must be the same for all models, including the true model, and therefore completely codify our unique information about the problem. If the averaging function for a position rc is very close to a delta function at that position, then its unit area ensures that the model average <m(rc)> is probably very close to the true.model value at that point. As well, as all models have the same unique value of the average then the range of model values at this position will be neccessarily limited. By constructing averaging functions at a variety of positions in the region of interest R, one can quickly determine those positions at which all models will be similar, and those at which the range of possible model values will be less constrained. —X In practice, because the averages <m(r&)> will always be calculated as a linear combination of the data (as in equation 4.2.6), the errors in the data will introduce error into the averages. Let the observed value of the data be e^, the true value, e;* , and the error in the observed value, <fet, such that: 117 (4.2.7) Similarly, we will have for the averages: <m(rD)> = <m(r0 )>* + £<m(r^)> (4.2.8) The assumption will be made that the errors are Gaussian, so that the expected value and covariance of the error terms are: E[<fe-] = 0 (4.2.9) COVfcfe- Jej ] = E[cf e Lfe - ] (4.2.10) We will also presume the errors are uncorrelated so that the covariance matrix, C, reduces to a diagonal matrix, with the diagonal elements being given by the square of the standard error 0"c of each datum e^: Using 4.2.7 in equation 4.2.6, we have: (4.2.11) 118 <m(rD)> =£^-eL KJ hi i-l t= i = <m<r0 )> + £*:«re<. (4.2.12) Thus N <f<m(r0 )> = ^oiLfe^ (4.2.13) The expected value of the error in the average is: E [ <m ( r0 ) > ] = E [ <f <m ( ) > ] = E[fe: ) = 0 The variance of the average then becomes VAR[<m(r0 )>] = VAR[<f<m( r0 ) > ] (4.2.14) 1 19 N N = E[{£*JeL)(£o<j£eS)] M N/ M N «-x,J5' i=, (4.2.15) The future algebra of the averaging function determination will be greatly simplified if we now reformulate the problem by normalizing with respect to the standard errors, c?c : ec(r) = e'(r)/^ = f${ [m(r )g{.(r) ]/oi }ds = ^ m(r )GC (r )ds R (4.2.16) VAR[<m(rc)>] = ^<^-VAR[ e^ ( r ) ] (4.2.17.) Backus and Gilbert (1970) show that resolution (as quantified by any suitable 'delta-ness' criterion) and accuracy are mutually antagonistic properties, so that any improvement in one neccessarily degrades the other. Thus, in each application one must choose the coefficents °<i so as to obtain the optimum 'trade-off between resolution and accuracy. To allow a 1 20 continuously variable trade-off the parameter © is introduced. The choice of G ,between 0 and TT/2, will determine the relative emphasis on resolution or accuracy, when the coefficients are calculated by minimizing the following objective function with respect to each o<L: Y(r0) = cos©S(r^) + sin9 VAR[<m(r0 )>] (4.2.18) Backus and Gilbert (1970) prove that the averaging function obtained in this way (where S would be any suitable 'closeness' criterion) will have the lowest possible value of the variance for a given value of the resolution. The trade-off curve of resolution and variance as a function of & thus obtained will be the optimal curve for the problem. A value of © = 0 will strictly minimize the 'closeYiess' of A(r,r0) to a delta-function, and the values of °<L determined will result in maximum resolution, but with a consequent maximum in the variance. On the other hand, 9 = Tr/2 will ensure that the vari-ance is a minimum, but now with resolution also at its minimum. We add the unimodular condition to the problem by the method of Lagrange multipliers, so that our objective function becomes: 121 V(r0) = 00565(4) + sin£vAR[<m(r0 )>] + 2B [ 1 ( r )ds ] R (4.2.19) with the values of °(i satisfying: 1 ^Gi(r)]ds = £*<L ^G-(r)ds (4.2.20) Minimizing with respect to each gives: Vfy3c<;= 0 = cos©(}S(re )/XL ) + sin6(V*<C){VAR[S<in(r0' )>]} - 2|ajy[}A(rfr0 )/M,)ds R (4.2.21) Substituting the expressions for S(rD), VAR[£<m(re )>], and A(r,r0) from equations 4.2.3, 4.2.15,and 4.2.2 into equation 4.2.19 we obtain for each o(; : ' -M cos0{ [£<*j f^Gj (r)G^(r)ds] - GL(r0)} + sine-*: -p^G^rJds = 0 R (4.2.22) To use simpler matrix notation, let: 122 ^G0(r)GL(r)ds = Pj = R (4.2.23) be an element of the NxN inner product matrix P , and: ^GL(r)ds = rjL (4.2.24) —* ~i be an element of the 1xN vector U, and G^(r0) and °<i be elements of the 1xN vectors G0 , and <\ respectively. Thus, considering the minimization with respect to each of the <<i's as in 4.2.21 the final matrix relation is: cos0[^-P-Go] + sin9-<* - pu = 0 (4.2.25) The unimodular constraint is now: _» -» o( • U = 1 (4.2.26) The inner product matrix, P , is obviously symmetric, and is also positive definite: 1 i P- xT =^^X- P; X = ii Jf(GL(r)xl)(GJ (r)Xj )ds L j R 123 = ff [<CGi (r)x- X^Gj (r)Xj )]ds R. L j = (^GL(r)xL )2 ds > 0 (4.2.27) Thus, r can be diagonalized (Parker, 1977): r = RART (4.2.28) where R and R are orthogonal NxN matrices such that: R- = RT (4.2.29) and A is an NxN diagonal matrix containing the eigenvalues of P , where all eigenvalues are greater than zero. Using this expansion, multiplying both sides by R and utilizing equation 4.2.29, 4.2.25 becomes: cos9[£-R-A - G • R] + sin6*R - B U-R = 0 (4.2.30) The matrix P can be considered to both 'rotate' and 'stretch' the components of an arbitrary vector; the decomposition of P in equation 4.2.27 separates these two operations, with the diagonal matrix A being pure 'stretching', and R and RT being pure rotations. We denote the resultant vectors of the rotation of <*,G0, and U by R as: 124 e< = oC- R (4.2.31) (4.2.32) U = U-R (4.2.33) Using these relations and collecting terms in equation 4.2.30, the rotated coefficients for the averaging function are found to be: c< = [BU + cos©G0]-D (4.2.34) where D is the diagonal matrix: D = cos 6 A + sinGl ( (4.2.35) Rotating two vectors will not change their inner product, so the unimodular condition of equation 4.2.26 is not changed: 125 oC . U = << • TJ = 1 (4.2.36) This equation allows the value of the Lagrange multiplier |3> to be calculated: p = [ 1 -cos6c?0- D '• U]/[U- D-1- U] (4.2.37) The values of the averaging function coefficients are then found by rotating equation 4.2.34 with the matrix RT: c< = c<- R = tpU + cosSS0 ]D • Rr (4.2.38) The utility of the diagonalization approach is now apparent, as the decomposition of the inner product matrix need be done only once for all values of © and r^, whereas otherwise a matrix inversion would be required for each different value. In the two dimensional curre'nt density problem the two forms of the kernels relating the data and the model for different station positions are given in equations 4.1.6 and 4.1.8. For a given array of stations over a region of interest R, the calculation of the averaging function at a position (x0,zo) and any © proceeds in the same manner as just outlined. The exact methods used to calculate the integrals neccessary for determination of the elements of the inner product matrix P and the elements of the vector U are given in Appendix G. Displayed in Fig. 4.2 are the averaging functions calculated in the above 1 26 Fig. 4.2 Averaging functions calculated for the two-dimensional current density problem. The averaging functions have been calculated at maximum resolution ( Q = 0) on the trade-off curve, and are thus as 'close' as possible to a Dirac delta function at the points: (a) (1,25) (b) (10,25) Fin. 4.2 1 27 (a) AVERAGING FUNCTION XLOC ~ 25 ZLOC ~ ] 6 = 0.000 Cb) AVERAGING FUNCTION XLOC = 25 ZLOC = 10 9 = o .ooo 128 manner for two positions in the region of interest R. Both averaging functions have been calculated at the highest resolution, 0=0. It is noted immediately that these functions are very spread out, particularly with respect to depth, z,and thus show very little resemblance to the delta function they are to emulate. The resolution at any other value of © would be even worse. It is also clear that the peak of the averaging function is much more localized in the 'x' direction than in the 'z' direction, indicating that the horizontal resolution of the true model's features will be superior to the resolution in depth. A more critical feature of the averaging functions is the shift of the peak of the function toward the surface away from its called for point,this being particularly noticeable in Fig. 4.2b. This large bias in the position of the averaging functions clearly invalidates the concept of the average as being a moment of the model around the desired point. Thus, in these cases, the actual averages would hold little meaning. A"s well, it further points out the extreme lack of depth resolution. The bias is due to the decrease in the required magnitude of currents in models fitting the data which are concentrated closer to the surface. This is a fundamental property of this current density formulation, and as such, is unavoidable through changes such as array design, etc. The computed averaging functions for the current density problem have thus shown that the non-uniqueness of the problem is indeed severe. However, although the resolution of the true 1.29 model's features with respect to depth is indicated to be extremely poor, there should be fair horizontal resolution. 4.3 Construction: The Smallest Model A customary first model to construct when inverting data is one with minimum structure. Here, that first model will be the smallest model, in a least squares sense, which fits the data. The objective function to be minimized is: (£(J) = ^ | J(x' ,z' ) |Zdx'dz' (4.3.1) To ensure that the current density model J fits the data, the data equations are added to the objective function via the Lagrange multipliers, pi: <f) (J) = JJ | J (x ' , z ' ) | dx'dz' R N + 2^^ [e. - JJG- (X' ,Z' )J(X* ,Z* )dx'dz' <--' R (4.3.2) However, the smallest model found by minimizing Cj^(J) in equation 4.3.2 would not be physically reasonable, as the 'smallest' criterion would ensure that all the currents would congregate in the uppermost portion of the region R. To offset 130 this effect, a weighting function, W(x,z) is introduced into the objective function: <^(J) = ^ W(x* ,z' ) | J(x' ,z' ) | dx'dz' + 2^ |3- [eL- G^(x',z')j(x' ,z')dx'dz'] (4.3.3) The exact form of the weighting will be dealt with later. Minimizing the objective function with respect to an arbitrary infinitesmal perturbation of the model gives: (f>(J+£j) - (J) = 0 = 2 Jfw(x' , z' )J(x' ,z' )«fj(x' ,z' )dx'dz' -2 £ B'^jd'J(x' , z' )GC (x' , z' )dx'dz' * (4.3.4) As cfj is arbitrary, this requires that: J(x,z) = £ )?• G: (X,Z)/W(X,Z) (4.3.5) However, the model must also satisfy the data, so that for all values of 'j': 131 ej = J(x' ,z' )Gj (x' ,z' )dx'dz' N rr = £ B- )\ [Gj (x' ,2')Gc(x',2')/W(x',2')]dx'd2' (4.3.6) or in matrix form: (4.3.7) where: Rl' = = ^ [ G^ (x ' , 2 ' )G^ (x ' , z ' ) /W(x ' , 2 1 ) ]dx'dz' (4.3.8) and p is the vector containing the Lagrange multipliers, fi . 1^ is symmetric and positive definite, so that P ' can be diagonali2ed: e = B • R-ART I A. "V <"V-(4.3.9) Thus, the values of fi are easily found: B = e • R A 'R (4.3.10) allowing the model J(x,2) to be calculated from equation 4.3.5. The required form of the weighting function W(xrz) must now be determined. Consider two line currents, I, ,Ia, of equal magnitude, but at different positions,(x,,z,) and (x2,z2) 1 32 respectively (see Fig. 4.3). The ratio of the Bx component of the first current to that of the second at a station (xo,0) is: Bx,/Bx2 = {[(x1-x0)2+z11 n/z^.U./Kx.-x^+z,1 ]} (4.3.11) This ratio could be made unity by multiplying the currents I, ,I2 v X by the weight factors W, , and Wj : W* (x,z) = zi/[ (x- -x0 )l +zL* ] (4.3.12) where i = 1 or 2. Following the same reasoning, the weighting function for the Bz component is found to be: WL*(x*,z) = (x--x0)/[ (xL-xe )l+Zi*] (4.3.13) Thus, the weighting function is different for the two components. As well, for each station position xD the weighting functions would change. Clearly there is no universal weighting function which will apply in all cases. For this reason an approximate weighting function to offset the depth effects only is suggested: W(x , z) = z * (4.3.14) where 'z' is the depth of the line current, and is some weighting factor. At large values of the station position x0, 133 Fig. 4.3 Determination of the appropriate weighting factor to offset the geometric decay in the magnetic field away from a line current. 134 the asymptotic values of the weghting functions of equations 4.3.5 and 4.3.6 are proportional to 1/z and 1/z1 respectively, so that a reasonable value of o< in equation 4.3.14 should be between 1 and 2. The numerical integration technique used to obtain the elements of P for the averaging functions of the previous section is easily adapted to include the variable weighting function in equation 4.3.14 (see Appendix G), so that the calculation of the smallest models for a given data set requires little extra programming. The artificial data to be used was generated from parameterized current density models (see Fig. G.1). The field components, Bx,Bz, at any station position due to a rectangular grid element of constant current can be calculated in closed form (see Appendix G). Because of the linearity of the equations involved, the contributions to the surface field from each grid element can then be summed, giving the resultant field due to the current model. Using these data sets, 'smallest' models were then constructed using the outlined method, for different values of the weight factor, <X . The first 'true' model used to generate artificial data was that shown in Fig. 4.4a. The models constructed using this data, for weightings of o( = 0, c< = 1, and <K = 2, are in Fig. 4.4b,c,and d respectively. In the unweighted case (°< = 0), the currents, as expected, are concentrated at the surface. However, the calculated model does indicate the proper horizontal position of the original model. In the second and third examples 135 Fig. 4.4 L^-Norm smallest model construction. The various models are: (a) the true model (b) model constructed at *< = 0. (c) model constructed at °< = 1 . (d) model constructed at <<= 2. co 1 37 the weightings are increased from oi = 1 to «< = 2, and the currents in the models in each case are 'pushed' deeper. Unfortunately, the models still show little resemblance to the true model, although they again have accurately indicated the 'true' horizontal position. As well, there is clearly no indication in the calculated models of the depth at which the true model might be. To further check the apparent ability to horizontally resolve features, the model in Fig. 4.5a was used as the true model, with the calculated models for the different values of °i being in Fig. 4.5b,c,and d. Once again the unweighted model is concentrated at the surface. It is noteworthy in this case that the horizontal features of the model are quite poorly resolved, which suggest that the sheet current models used by Banks (1979) arid Woods (1979) might run into similar difficulties. The weighted model constructions for this example quite clearly show the horizontal positions of the true model currents, but again give no indication of the depth. The final model construction was a further test of the presumed lack of vertical resolution. The true model is shown in Fig. 4.6a, and the weighted models are in Fig. 4.6b,c,and d. The results confirm the previous conclusions, as the two line currents are not resolved at any weighting value. The weighted smallest model constructions done here have substantiated the predictions made from the Backus-Gilbert appraisal. The general horizontal features of the true model are all reproduced to some extent in the calculated model, whereas 138 Fig. 4.5 L^-Norm smallest model construction. The various models are: (a) the true model (b) model constructed at <* = 0. (c) model constructed at <* = 1 . o (d) model constructed at <* = 2. 139 140 Fig. 4.6 L£-Norm smallest model construction. The various models are: (a) the true model (b) model constructed at «< = 0. (c) model constructed at °C = 1 . (d) model constructed at = 2. GO 142 the vertical features are not resolved at all. 4.4 Constrained Model Constructions Using Linear Programming It has been shown in the previous sections that the surface data do not limit the range of possible models enough to give vertical resolution of the true model's features. To further restrict the range of permissible models, it is suggested that expected physical features of the model be incorporated, or favoured, in the model pohstruction. In a great number of cases the anomalies which attract attention in G.D.S. are those which appear to be very localized, as these suggest interesting geophysical structures such as geothermal hotspots, ancient craton boundaries, fault zones, etc. This localized type of model consisting only of a few large model elements will not be favoured by the least-squares L2-norm used in the smallest model calculation, which is clear from the flattened, spread-out models calculated in the previous section. The norm which most favours construction of the desired sparse, localized models is the L,-norm (Levy and Fullagar, 1981). Another physical feature is that no currents in the model will be expected to be more than IT/2 different in phase. From equation B.9 in Appendix B, the complex vertical wavenumber in a homogenous earth is: 143 = i( ^+ iw^cr)"1-(4.4.1) where is the horizontal wavenumber, 2rt"/;\ (as defined in section 2.1 of Chapter II). The real part of kJe will define the oscillatory or 'travelling' nature of the wave, so that the depth at which the wave has travelled a quarter wavelength from the surface is given by: z'/4 = ™/2 {REAL[i(^>v+ ivft^f1]} (4.4.2) The value of z y+ will increase with increasing "v , so it will always be true that: z,4 > IT/2 (2/w/U0CT )"1 = 1 . 57 <f (4.4.3) where tT is the skin depth. Thus, the induced currents that are more than out of phase with those at the surface will be insignificant, as they are deeper than the skin depth of the inducing wave. (It is noted that in the case of currents that are channelled large distances vertically upward, this argument will break down.) To reiterate, the two expected physical features of the model are: 1. The current model will be sparse, that is, it will 144 consist of very localized current elements. 2. No currents in the model will be more than fr/2 out of phase. Both of these constraining features are easily applied using linear programming, with an Lt-norm objective function to be minimized. The current density model is parameterized into an N x M grid of rectangular elements, with each having a constant current within them (see Fig. G.1). The contribution to the surface field at (xk,0), due to a current J:J in the grid element (i,j), can be calculated (see Appendix G), and is linearly proportional to J;; : Presume now that the currents in all grid elements are in phase. ( The possibility of phase differences will be considered later). The linearity of the current density equations then allow the total field at any surface position (xk,0) to be calculated from the sum of each grid contribution: BM '<*k> (4.4.4) (4.4.5) 145 N tt B*(xk) = i $ J. • ki{ (x. ) (4.4.6) (4.4.7) Considering all the surface stations (x^: k = 1fD and both directional components results in an L x (2NM) set of linear equations, which in most cases will be underdetermined (that is, 2NM > L). A model to fit these equations can be found using linear programming. This model will be the smallest model in the L,-norm sense, if the model minimizes the objective function: <£(j) £ IJC, i (4.4.8) As well, because the L, -norm results in sparser, more localized models than the least-squares norm (L^-norm), the calculated model will be in accordance with the first expected physical feature. Linear programming solves the underdetermined problem with only positive values of the variables allowed, which is why the formulation to this point has presumed the grid element currents are all in phase. The desired condition was a much weaker one; that there be no currents in the model more than tr/2 out of phase. Presume now that the currents in each grid element (i,j) at each frequency 'w' are not in phase, with their complex time 146 domain representation being: (t) - |Jy | .«««**^ (4.4.9) and the corresponding frequency domain representation being: Jij (w) = |Jij | cos + i|Jij | sin©^ (4.4.10) for w > 0. If no currents are more than Tr/2 out of phase, then it is ensured that a minimum value of <9cj can be found such that [ ©:J - G^-.rv] > 0 always. Phase shift all the currents by this minimum value of B : (t) = Jcj (t) e"18^ = |Jij I e J (4.4.11) In the frequency domain the imaginary and real components of each current element (i,j) are: Jij (w) = | cos( ©Lj -<9~In) + i|J^| sin ( ©cj -(4.4.12) Thus, as COS(S;J and sin( G^-tSU,'*.) are guaranteed to be greater than zero, then the values of the real and imaginary parts of all currents (w) will be greater than zero. The linear programming approach as outlined can then be applied 147 separately to the real and imaginary portions of the phase-shifted data,, and the constructed models will be in accordance with the second expected physical feature. The final real and imaginary portions are then recombined to give the amplitude and phase of each model element. The original phase-shift of the data is subtracted from the phase values to give the final phase representation of the model. Using the unweighted L(-norm 'smallest' criteria as in equation 4.5.8 will result in the same unphysical shallow models as in the previous least-squares construction. To offset this the weight function W(z) is again included in the objective funct ion: As well, in practice there will be some error in the data, so to fit the data equations of 4.5.6 and 4.5.7 exactly would be overfitting the data. This is handled by expanding each equation into two inequality constraints. For example, equation 4.5.6 becomes: (4.4.13) BK(xk) +*\(xk) >t£ J^A*j(xk) (4.4.14) 1 48 BK(xk) -^Bx(xk) <£ £ JL-ACi(xk) (4.4.15) where cTB)C(xk) is the estimated standard deviation in Bx (x^). The model would then be required to fit the data only to within the standard deviation of each data point. This completes the formulation of the -norm linear programming construction. The linear programming construction was applied to the same suite of true models as used previously. Because the method of calculating the data from the parameterized true model is identical to the method used to calculate the contribution of each grid element in the construction, the grid meshes used in each were made dissimilar, with the grid size in the data calculation being 20 x 20, and that in the construction, 16 x 16. The first set of true and calculated models are shown in Fig. 4.7. In all cases, only the amplitude current density model is given. In Fig. 4.7, the calculated models are constructed at various values of the weight factor c< , for surface fields calculated using the true model. Gaussian white noise with a standard deviation of 10% of the maximum field value was added to the data. Because of the rapid geometric decay in the field of the anomalous current, this noise level will be very high at stations to either side of the anomaly, and is thus unreasonably large. The unweighted model is concentrated near the surface as expected. However, for both of the weighted constructions (o< = 1, = 2) the calculated models are very similar to the true Fig. 4.7 Lj-Norm linear programming model construction. 10% noise was added to the true model data.' The various models are: (a) the true model (b) model constructed at = 0. (c) model constructed at °( - 1 . (d) model constructed at o< = 2. o 151 model, with the horizontal and the vertical features recovered. When the grid size is taken into account, the total current in the large peak is also found to be close to that of the true model. It is an important point that both values of the weighting produce similar results, as this removes the problem of determining the 'best' value of o( for each data set to be inverted. In fact, values as large as 10 were tried for ^< , and caused little difference in the constructed model, indicating great stability with respect to the weighting. In all subsequent models this stability is taken advantage of, with the value of °< f ixed at 1. The horizontal and vertical resolution of this construction were then checked using data from the true models of 4.8a and 4.9a, to which 1 % noise was added. The construction was done using a weighting of oC = 1, and the results are in Fig. 4.8b and 4.9b. The resolving ability in both the horizontal and the vertical direction is evident. The true models used to this point have all matched the first presumed physical feature, in which the true model was assumed to be sparse, with only a few localized current elements. To check the limit of the L(-norm construction with respect to less localized true models, the models of Fig. 4.10a and 4.11a were used to generate the surface data. The constructed models in Fig. 4.10b and Fig. 4.11b do show the tendency of.the I^-norm to over-localize the current elements. However, the data still supply enough information so that the 1 52 fil 4.8 Linear programming model construction to test horizontal resolution. The constructed model was calculated with a weighting of <X = 1. 1 % noise was added to the true model data. (a) the true model (b) the constructed model 1 53 Fig. 4.9 Linear programming model construction to test vertical resolution. The constructed model was calculated with a weighting of °< = 1 . 1% noise was added to the true model data. (a) the true model (b) the constructed model 154 Fig. 4.10 Linear programming model construction to test the algorithm with less localized true models, in this case a vertical current dike. The weighting was o< = 1 , and 1% noise was added to the true model data. (a) the true model (b) the constructed model 155 Fig. 4.11 Linear programming construction to test the algorithm with less localized true models, in this case a dipping dike. The weighting was c< = 1, and 1% noise was added to the true model data. (a) the true model (b) the constructed model 156 main features of both dike-like anomalies are recovered, including their dip. In summary, this linear programming modelling approach is found to be stable with respect to both the weighting function and to reasonable amounts of noise in the data. It has good vertical as well as horizontal resolution when the true models consist of a small number of localized model elements. Even for dike-like anomalies which are localized only in one direction, the construction still recovers the major features of the true model. The data sets used to test the linear programming construction algorithm have to this point been generated from current density models, with all current elements in these models being in phase. A more realistic test is to use as the starting data the output from the Jones' forward induction program, where now the true' models would be conductivity models. This will in particular check the construction algorithm's ability to handle data containing various phase shifts, and will also test the correlation between the constructed current structure and the original conductivity model. The first conductivity model used to generate data was that in Fig. 4.12a, consisting of a localized conductor in a layered structure, with a horizontal position of 189 km. and a depth of 35 km. Because of the inconsistency between the two directional components of the output in Jones' program (as noted in section 2.2 of Chapter II), only the Bx component has been used for 157 SURFACE POSITION (KMi O (OO 200 300 400 Fig. 4.12 Linear programming model construction using as starting data the output from the Jones-Pascoe induction program. The models are: (a) the true conductivity model (b) the constructed current density model 158 SURFACE POSITION (KM) O 100 200 300 400 >99, 'yp-Fig. 4.13 Linear programming model construction using as starting data the output from the Jones-Pascoe induction program. The true conductivity model in this case consists of two vertically displaced conductors. (a) the true conductivity model (b) the constructed current density model 159 these inversions. The constructed amplitude current density model is shown in Fig. 4.12b. The basic structure of this current density model is a close reproduction of the original conductivity model, and the depth (34.6 km.) and horizontal position (188.1 km) also match that of the original model. A second conductivity model was used to check the vertical resolution, with two localized conductors at 9 and 35.5 km. depth, with both at a horizontal position of 189 km. (see Fig. 4.13a). Again, the constructed current model (in Fig. 4.13b) has resolved the major features of the conductivity model, and has accurately delineated the horizontal position (188 km.) of the original conductors. As well, a reasonable estimate of the depths of the conductors (8.3 km., 30.2 km.) has been made. Chapter V Analysis of G.D.S. Across the Cascade Anomaly G.D.S. data from an array of magnetometers spanning the Cascade anomaly of Washington State (Law et al, 1980) has been obtained from the University of Washington's Geophysics Section. The Cascade anomaly was originally detected by a linear array of stations in southwestern Washington (see Fig. 5.1) with the induction arrow responses indicating a localized north-south conductor between the stations KOS and WHI (see Fig. 5.2). The coincidence of this conductive path with the Cascade volcanic belt generated significant interest, prompting further investigations to trace the course of the conductive path to the north and to the south. One such investigation recorded the surface magnetic field of a polar magnetic storm in February of 1980, at the array sites indicated in Fig. 5.3. The data from this event (see Fig. 5.4 and 5.5) is used here to demonstrate the L,-norm linear construction programming routine described in section 4 of Chapter IV. The discussions of appraisal and construction in Chapter IV were all' based on the presumptions that the normal and anomalous 161 Fig. 5.1 The array of magnetometers used by Law et al (1980) to study the Cascade anomaly of Washington State. 162 Fig. 5.2 In-phase and quadrature induction arrows for 30,300, and 3000s. (After Law et al, 1980). 163 Fig. 5.3 The array of magnetometers which measured the data used in this thesis (after Hensel, 1980). The approximate continuation of the Cascade anomaly (Hensel, 1980) is marked as a dashed line. The study area is indicated by a rectangle in Fig. 5.1. 1 64 (a) . 5.4 Magnetograms for the magnetic storm of February, 1980. The measuring stations are: (a) NIS (b) ORT 1 65 (A) Fig. 5.5 Magnetograms for the magnetic storm of February, 1980. The measuring stations are: (a) MUD (b) GRE fields had already been separated, and that the subsurface structure was two dimensional. The linear programming construction further assumed that the internal currents producing the anomalous field were very localized and that none of these currents were more than 'TT/2 different in phase. Thus, before any construction of models is done, the separation of the total field into its normal and anomalous portions must be performed, and then all the assumptions must be checked using the anomalous field values. The normal field defined in Chapter III was the total surface field that would be induced over the regional one dimensional structure. However, the stations in this array are all very close to the land-sea boundary, so that the geomagnetic coast effect (Parkinson, 1959; Everett and Hyndman, 1967) will be pronounced at all of them. Although it is of interest and is often studied, the coast effect in this experiment simply obscures the desired response from the inland Cascade conductor. For this reason it is advantageous to include the anomalous field of the coast effect in the defined normal field, so that upon subtracting this normal field from the total field only the anomalous portion due to the Cascade conductor remains. Law et al (1980) have done this using the calculated response of the coast effect at each inland station.(Everett and Hyndman, 1967). However, for the array supplying the data for this thesis the coast effect will be complicated by the proximity of the Puget Sound. As well, the two outermost stations in this array, NIS 167 and GRE, are only separated longitudinally by about 80 km. Because of these points, the simplification will be used here that the coast effect is constant across the array. The strong similarity between the magnetograms at the stations which are most distant from the presumed location of the Cascade anomaly (NIS and ORT), indicate that this simplification is reasonable. The field at the outermost station, NIS, is designated the normal field. The anomalous field at the remaining three stations is then calculated as the difference between the individual measured fields and this normal field, and both the one-dimensional structure response as well as the coast effect should be removed. The first step in the analysis of the anomalous field (henceforth called simply the field) at the stations ORT, MUD, and GRE will be the calculation of the spectra for each directional component. The first 23.75 hours of the digitized magnetograms were used, which was 4096 points at a 20 second digitizing interval. The mean and linear trend were removed from each signal and the Fast Fourier Transform was applied. The amplitude spectra are given in Fig. 5.6, 5.7, and 5.8. The original measurement of the magnetic field was digital also, with each count representing .25 nanotesla (nt.). Thus, an intuitive 'threshhold level', below which the amplitude spectra might be considered to be merely noise, would be .125 nt. (The half count value is taken because there is an equal amplitude of signal at the corresponding negative frequency). If this level 168 H SPECTRUM . D SPECTRUM Fig. 5.6 The amplitude spectra of the anomalous field at ORT-(a) H (b) D (c) Z 169 H SPECTRUM D SPECTRUM CO Fig. 5.7 The amplitude spectra of the anomalous field at MUD-(a) H (b) D (c) Z H SPECTRUM D SPECTRUM Z SPECTRUM The amplitude spectra of the anomalous field at GRE (a) H (b) D (c) Z 171 is used, the meaningful range for most of the spectra would only include periods greater than ~1 hour. However, the strict use of this 'threshhold' level is not appropriate, because the inducing fields are more correctly modelled as a sum of frequency components with time varying amplitudes; b(t) = £ a•(t) cos(w-t) (5.1) rather than as a sum of frequency components with constant amplitude, as modelled by the Fourier transform. The time variations of the amplitude, aj^t), will spread the energy of each frequency component of the signal over a range of frequencies in the Fourier transform representation. As an example, for a signal b(t) as in equation 5.1, but with a single frequency component, w0 , its frequency spectrum will be: B(w) = [Ac(w)/2] «> [<f(w-w0)] + [AQ(w)/2] ® [cf(w+wc)] (5.2) with: A6Cw) = J[ac(t)] (5.3) Thus, the threshold level can probably be taken to be much lower than half the count rate. The possibility that the model of equation 5.1 is more physically representative of the source suggests that complex demodulation (Banks, 1975) may be a more 172 reasonable approach to determining spectra than the Fourier transform approach, but it has not been used here. The assumptions embodied in the construction routine of section 4 in Chapter IV can now be checked. A measure of the degree of polarization of the magnetic field vector at any frequency 'w' can be obtained from (Samson, 1977): p(w) = [nTr(SZ) - (TrS? ]/[(n-1)(TrS)^ ] (5.4) where p(w) is the degree of polarization at 'w', 'n' is the number of vector components, and S(w) is the smoothed spectral matrix of the magnetic vector. The unsmoothed spectral matrix is given by: S' (w) with S'^-g, (w) defined as: S^fc (w) = A(w)B*"(w) (5.6) (A(w) and B(w) are the complex magnitudes of the spectra for each vector component). The value of the degree of polarization can vary between 0 and 1, corresponding to completely random S^ (w) SiD (w) S^ (w) \S4H(w) S', (w) S^(w), (5.5) for two arbitrary vector components A and B being 1 73 polarization and to a completely polarized signal. Using the spectra previously determined, the spectral matrix S'(w) was calculated, and subsequently smoothed with a 3 point filter (.25,.5,.25) to give S(w): Real[SAB(w)] = Real [ S 'AE (w) ] ® (.25,.5,.25) (5.7) Imag[SftB(w)] = Imag[S'(qgw(w) ] ®. (.25,.5,.25) (5.8) The smoothing of the spectral matrix before the degree of polarization is calculated is based on the implicit assumption that the field is more correctly modelled by equation 5.1, rather than by the Fourier transform expression. If the unsmoothed spectral matrix was used, the degree of polarization from equation 5.4 would be 1 at every frequency, in accordance with the constant amplitude and constant phase assumption of the Fourier transform model. The values of p(w) are plotted in Fig. 5.9 for the stations ORT, MUD, and GRE. The magnetic vector at all stations is strongly polarized for nearly all frequencies, with the vector at MUD and ORT being almost completely polarized. This allows us to approximate the frequency components of the magnetic vector at each station position 'x' by a completely polarized signal (Born and Wolf, 1975, pg. 32): 174 Ca) PERIOD (MIN.) 50 25 16 7 4- r—-I I I M O.OO CYC.I F.C/MIN 10 _L_ ft) «•• i.n • •M O.H | .Qg CYCLES/ MIN. • •I t.ll CO 50 I PERIOD (MIN.) z? 'V v> CYCLES/MIN. Fig. 5.9 The degree of polarization between the three directional components of the anomalous field as a function of period. (a) ORT (b) MUD (c) GRE 175 H(x,w,t) = HQ(x,w.)- [cos(wt+ <}>H(x,w) ) ] (5.9) D(x,w,t) = De(x,w) • [cos(wt + <pt(x,w) ) ] (5.10) Z(x,w,t) = ZQ(x,w) • [cos(wt+<r\(x,w) ) ] (5.11) As indicated, H0, De, Z0,9H rand 9* are all independent of time. The varying position of the vector described by these three components will trace out the surface of an ellipsoid, with principle axes having lengths A,B, and C, and directions the polarization of the vector. The procedure to obtain the principle axes lengths and directions is detailed in Appendix H. Basically, the three component equations of 5.8 - 5.10 are converted to a quadric equation representing an ellipsoid, which is then simplified by diagonalizing the matrix of coefficients. The relative lengths of the principle components A,B, and C are then given by the eigenvalues (Xi.) of the coefficient matrix, which are the non-zero elements of the diagonal matrix: . These parameters completely determine the nature of 176 A i/>;'-(5.12) B «< 1/-X,,"1 (5.13) C « 1/W1 (5.14) The direction of each principle axis is given by the eigenvector associated with its eigenvalue. The spectra for each station were smoothed using the three point filter (.25,.5,.25),the eigenvalues were determined, and the ratios of the smallest to second largest eigenvalues were calculated at each frequency. The results given for the periods from 20 min. - 4 hrs. (see Fig. 5.10 ) indicate that at both ORT and MUD the anomalous field vector is effectively linearly polarized. The high degree of polarization found in the anomalous field at all three stations, and in particular ORT and MUD, could be due only to either a highly polarized inducing field, or to a direction of constant symmetry in the conductive structure. From the original magnetograms in Fig. 5.4 and 5.5 it appears that the inducing field is in fact polarized, at least for the long periods. However, even for a linearly polarized 177 PERIOD (MIN) Ca) e.i in Cb) CYCLES.'MIN. CO o.D li oi CYCI.ES/MIN. O.I a ii Fig. 5.10 The ratio of the smallest to the second smallest eigenvalues from the matrix of coefficients of the polarization ellipsoid, as a function of period. (a) ORT (b) MUD (c) GRE 178 inducing field, if the conductivity structure is three-dimensional the anomalous field will not be expected to be linearly polarized, because of the probable phase shifts between the directional components. Thus, it is felt that the linear polarization at ORT and MUD indicates that the structure has a direction of constant symmetry. Further to this, as shown in Fig. 5.11, the fact that two stations have a linearly polarized signal suggests that the current must be localized at a single depth, as currents at different depths would introduce a phase shift between the directional components at at least one station. If the field is presumed to be due to a line current, the vector directions of the longest principle axis at any two stations must both be perpendicular to the line current, and also must be perpendicular to the normals to the line current which intersect the stations. This would completely determine the location of the line current, including the dip and strike, if the vector directions were perfectly accurate. Using the directions of the largest principle axis from the associated eigenvectors, the dip and strike of the conjectured line current was calculated using the ORT and MUD eigenvectors. The mean values and standard deviation of the dip and strike of the proposed line current path were found from the 1 hr. - 20 hr. period range to be (see Fig. 5.12): DIP = -1.001° + 7.117° 179 AIR Stations A,/ H, ^FkX.B < I = W EARTH X \ \ lei<P> Is = ir » > Fig. 5.11 Two lineal conductors are shown at different depths. As shown in section 2.2 of Chapter II, the resultant field from each of these conductors can be modelled as a line current, as shown at stations A and B. As the conductors are at different depths, there will be a phase difference between their respective fields, so that both stations cannot have a linearly polarized signal. 180 DIP 50 Period (min.) 25 16.7 12.5 10 0.12 STRIKE o Cycles/mm. 0.]<! Fig. 5.12 The estimated dip and strike that a line current would need to match the field direction of the largest principal axis of the polarization ellipsoid at stations ORT and MUD. 181 STRIKE = -15.029° + 5.825° where the dip is downward in the direction of the strike, and the strike is measured clockwise from north. The dip is small enough to be ignored, thereby allowing the use of the two-dimensional approximation. The assumptions of Chapter IV have been validated for this data set, except for the assumption of the maximum phase difference. The check of the phase differences was incorporated in the modelling routine itself, and it was found that the assumption was not violated. The remaining step is to rotate the vector data and station positions to a new coordinate system aligned along the estimated strike of the structure, using the station NIS as the origin of the new coordinate system. The error in the angle of the strike will introduce error into both the station position, and the rotated magnetic data, with the resultant error in the magnetic field values easily incorporated into the linear programming method. There is no easy way to incorporate the error in the station position into the linear programming, but its relative value will be small so it has simply been ignored. The final station positions (x) are given in Table 5.1. 182 Table 5.1: Station Positions Station Position NIS 0.0 km. ORT 33.2 km. MUD 57.7 km. GRE 80.9 km. Using the one-dimensional station positions (x), the component of the rotated field perpendicular to the proposed strike, and the vertical component, the construction of models was performed using the linear programming method of Ch. IV. For the longer periods (1 - 4 hrs.) the constructed models consistently require a localized current in the depth range 11 -21 km., at a station position, x = 50.6 + 2.8 km. (see Fig. 5.13a,b,and c). Varying the size of the region of interest, or the 'tightness' of fit to the data can result in the removal or appearance of the currents at the edges of the region, but never 183 Fig. 5.13 Linear programming current models calculated using the real data at a variety of periods. The periods of each model, with the indicated position of the central line current, are: (a) 4 hr. Depth 16.9 km. Station Position ... 46.5 km. (b) 2 hr. Depth 14.1 km. Station Position...50.7 km. "(c) 1 hr. Depth 21.0 km. Station Position ... 52.5 km. (d) 20 min. Depth 11.0 km. Station Position ... 52.5 km. Fig.5..IT 4 hr 2 hr 1 hr 20 min. 185 results in the disappearance of this central current. This indicates that the other current elements are strictly artifacts of the noise and the modelling parameters, whereas the central current is required to satisfy the data. It is noted that the depth of this current (15.7 ± 4.2 km.) is in good agreement with the value of 17+8 km. estimated for the southern section of the Cascade anomaly, using a best fitting line current model for a period of 1 hr. (Law et al, 1980). For the linear programming models constructed at higher frequencies, the central current, although still evident, becomes less and less significant (Fig. 5.13d). For periods less than 15 min. it was not always possible to find a model that matched the data, indicating that either the two-dimensional symmetry assumption was no longer valid at these periods, or that the currents were perhaps travelling in a different direction. The induction arrow results of Law et al (see Fig. 5.2) are in accordance with this, as they indicate the disappearance of the north-south currents between KOS and WHI for periods less than ~1 hr. The appearance of the Cascade anomaly currents only for the longer periods ( 30 min. - 4 hrs.) is at odds with its shallow depth. As indicated in Table 1.1, the skin depths for these periods at a conductivity of .01 S/m (which is very high for shallow crustal rocks) are on the order of hundreds of kilometers. This suggests that the currents are probably due to channelling of a regional current system, and are not due to 186 local induction. The proposed direction of the segment of the Cascade anomaly studied in this thesis is shown in Fig. 5.14, along with the southern section from Law et al, and an approximate northern section from Hensel (1981). The implication of this 'complete' path is that the Cascade anomaly currents are simply being channelled along some relatively conductive path into the Puget Sound. This appears to rule out the possibility of a conductive conduit connecting the chain of Cascade volcanoes. 187 Fig. 5.14 A -final map of the Cascade anomaly showing the approximate positions of portions as estimated by: (a) Law et al (1980) » » • »— (b) This thesis —I—i—i—i— (c) Hensel (1981) 188 CONCLUSIONS A detailed investigation has been made of the assumptions and limitations inherent in the traditional methods used in G.D.S. The approximation of an induction tensor linearly relating the normal and anomalous fields is found to be probably valid at midlatitudes, for periods greater than ~2 hrs. However, as shorter periods are used, or in regions where the inducing field cannot be considered horizontally uniform, its validity will degrade, and ultimately will fail. The induction arrows derived from the induction tensor will share the same limitations. As well, it has been shown using very simple examples that the arrows will not neccessarily point towards current concentrations, but rather will point towards high relative conductivities, for both induced and channelled types of anomalies. The various quantitative measures commonly employed in G.D.S. to determine depths of currents, lateral extents of anomalies, and scale lengths, have been shown to be of limited usefulness, and in fact, one estimator for the scale length (Porath et al, 1971) is revealed to be erroneous. It is also concluded that quantitative modelling of the conductivity structure at this stage in its development is not always the most practical way to proceed, as the modelling is expensive and 189 time consuming because of the non-linearity of the induction formulation. As well, because of the possibility that the anomaly is due to channelling of regionally induced current systems, the extent of the region over which the modelling is to be done is always in doubt. The failure of modelling efforts in certain studies has been blamed on this inability to find the proper region of interest (Whitham and Andersen, 1965; Porath et al, 1971). To avoid these difficulties and to also put the problem in linear form, it is suggested in Chapter IV that the anomaly be modelled in terms of current density rather than conductivity, with an initial simplification of two-dimensionality. The region of interest of the model is easily determined in this formulation, and the linearity allows for fast and inexpensive computations for both the forward and inverse problems. The major disadvantage of the current density approach is the non-uniqueness that is inherent in this approach. The extent and type of this non-uniqueness is explored using the averaging functions of Backus - Gilbert appraisal (1967,1968,1970). It is found that the surface data will constrain the range of possible models such that resolution of the major horizontal features of the true model will be apparent in all constructed models. However, the averaging functions indicate that constraining the model construction with the data alone will result in no resolution of the true model's vertical features. These conclusions were confirmed by the construction of the 190 L^-norm weighted smallest models, with the different models that fit the data all correctly indicating the horizontal positions of the current elements in the true model, but with none of the models giving an accurate vertical placement. To overcome this uniqueness difficulty, certain expected physical features of the true model were incorporated or favoured in the model construction. It was suggested that many anomalies would be localized and would be sparsely distributed, and so a weighted L,-norm objective function was introduced for the construction of a parameterized current density model using linear programming. The currents were also presumed to be due mainly to first order induction, so that no significant currents in the model would be more than Tr/2 different in phase. With these constraints and model features incorporated in the linear programming construction, it was found that both vertical and horizontal resolution of the features of localized, sparsely distributed true models was now possible. This ability of the construction algorithm to recover the major features of the true model was found to be stable with respect to the weighting factor (for << £ 1), which allowed this factor to be fixed at << = 1 for all subsequent modelling. As well, the success of the algorithm persisted in the presence of reasonable amounts of noise (up to 10% of the maximum data value for simple models; and 1-2% of the maximum data value for more complex models). Even when data was inverted using as true models current density configurations that were localized only in one direction, the 191 resultant constructed models were still a good replication of the true model. In the final chapter the linear programming construction was applied to real data measured at stations crossing the Cascade anomaly in Washington State. The data was checked, and found to be in good accord with the assumption of two-dimensionality. As well, the nearly linear polarization of the signal at two stations indicated that the anomaly probably consisted of only a single localized current. After rotating the data and station positions into the proper reference frame, current density amplitude models were constructed over a range of periods from 20 min. to 4 hrs. and it was found that all models required a localized current near the station MUD, at a depth of from 11-21 km. The mean value of depth of 15.7 + 4.2 km. from all the models is in good agreement with the estimated depth of Law et al (1980) for the southern portion of the anomaly. The three suggested segments from Law et al (1980), this thesis, and Hensel (1981) show good continuity, and indicate that the current path does not extend to the northern Cascade volcanoes. The shallowness of the significant currents in the models in comparison with the skin depths at their respective periods indicates that the anomaly is probably due to the channelling of regional currents through a local high conductivity feature. 1 92 BIBLIOGRAPHY Akasofu, S-I., (1979), 'Dynamics of the Magnetosphere', Akasofu (ed.), D. Reidel Publ. Co., pg. 447-460. 'What is a Magnetospheric Substorm' Akasofu, S-I. and Chapman, S., (1961), J. Geophys. Res., 66 , pg.1321-1350. 'The Ring Current, Geomagnetic Disturbances, and the Van Allen Radiation Belts' Alabi, A.O., Camfield, P.A., and Gough, D.I., (1975), Geophys. J. R. astr. Soc, 43 , pg. 81 5-833. 'The North American Central Plains Conductivity Anomaly' Alfen, H., (1950), 'Cosmical Electrodynamics', Oxford University Press, Oxford. Anderssen, R.S., ( 1975), Phys. Earth Planet. Inter., J_0 , pg. 292-298. 'On the Inversion of Global Electromagnetic Induction Data' Backus, G. and Gilbert, F., (1967), Geophys. J. R. astr. Soc, J_3 , pg. 247-276. 'Numerical Applications of a Formalism for Geophysical Inverse Problems' Backus, G. and Gilbert, F., (1968), Geophys. J. R. astr. Soc, 16 , pg. 169-205. 'The Resolving Power of Gross Earth Data' Backus, G. and Gilbert, F., (1970), Phil. Trans. Roy. Soc, A266 , pg. 123-192. 'Uniqueness in the Inversion of Inaccurate Gross Earth Data' Bailey, R.C., (1970), Proc. Roy. Soc. Lond., A315 , pg. 185-194. ' Inversion of the Geomagnetic Induction Problem' 194 Banks, R.J., (1969), Geophys. J. R. astr. Soc, _T7 , pg. 457-487. 'Geomagnetic Variations and the Electrical Conductivity of the Upper Mantle' Banks, R.J., (1973), Phys. Earth Planet. Inter., 7 , pg. 339-348. 'Data Processing and Interpretation in Geomagnetic Deep Sounding' Banks, R.J., (1975), Geophys. J. R. astr. Soc, 43^ , pg. 83-101. 'Complex Demodulation of Geomagnetic Data and the Estimation of Transfer Functions' Banks, R.J., (1979), Geophys. J. R. astr. Soc, 56 , pg. 139-157. 'The Use of Equivalent Current Systems in the Interpretation of Geomagnetic Deep Sounding Data' 195 Beamish, D., (1977), Geophys. J. R. astr. Soc, 50 , pg. 311-332. 'The Mapping of Induced Currents Around the Kenya Rift: A Comparison of Techniques' Becher, W.D. and Sharpe, C.B., (1969), Radio Science, 4 , pg. 1089-1094. 'A Synthesis Approach to Magnetotelluric Exploration' Birkeland, K., (1908), 'The Norwegian Aurora Polaris Expedition 1902-1903, Vol. 1 , Section 1 , Aschhoug, Christiania. Booker, H.G. and Clemmow, P.C., (1950), Proc I.E.E., 97 , Part III, pg. 11-17. 'The Concept of an Angular Spectrum of Plane Waves, and its Relation to that of Polar Diagram and Aperature Distribution' Born, M. and Wolf, E., (1975), 'Principles of Optics', Pergamon Press, Toronto. 196 Bostrom, R., (1964), J. Geophys. Res., 69 , pg. 4983. 'A Model of the Auroral Electrojets' Brace, W.F., (1971), 'The Structure and Physical Properties of the Earth's Crust', ed. J.G. Heacock, Geophysical Monograph 4, American Geophysical Union. Bracewell, R.W., (1965), 'The Fourier Transform and its Applications', McGraw-Hill, New York. Budden, K., (1961), 'The Wave-Guide Mode Theory of Wave Propogation', Prentice-Hall, N.J. Cagnaird, L., ( 1 953), Geophysics, J_8 , pg. 605-635. 'Basic Theory of the Magnetotelluric Method of Geophysical Prospecting' Chapman, S., (1919), Phil. Trans. Roy. Soc. London, A218 , pg. 1-118. 'Solar and Lunar Diurnal Variations of Terrestrial Magnetism' 197 Chapman, S. and Price, A.T., (1930), Phil. Trans. Roy. Soc. London, A229 , pg. 427-460. 'The Electrical Magnetic State of the Interior of the Earth as Inferred from Terrestrial Magnetic Variations' Cochrane, N.A. and Hyndman, R.D., (1970), Can. J. Earth Sci., 6 , pg. 1208-1218. Dragert, H., (1973), Ph.D. Thesis. 'Broad-band Geomagnetic Depth-Sounding Along an Anomalous Profile in the Canadian Cordillera' Dyck, A.V. and Garland, G.D., (1969), Can. J. Earth Sci., 6 , pg. 513-516. 'A Conductivity Model for Certain Features of the Alert Anomaly in Geomagnetic Variations' Eckhardt, D., Larner, K., Madden, T., (1963), J. Geophys. Res., 68 , pg. 6279-6286. 'Long-Period Magnetic Fluctuations and Mantle Electrical Conductivity Estimates' 198 Everett, J.E. and Hyndman, R.D., (1967), Phys. Earth Planet. Inter., J_ , pg. 24-34. 'Geomagnetic Variations and Electrical Conductivity Structure in South-Western Australia' Fischer, G., Schnegg, P.-A., and Peguiron, M., (1980), to be published. • 'An Analytic One-Dimensional Inversion Scheme' Fulks, W., 'Advanced Calculus', John Wiley and Sons Inc., U.S.A. Fullagar, P.K., (1980), Ph.D. Thesis. 'Inversion of Horizontal Loop Electromagnetic Soundings Over a Stratified Earth' Garland, G.D., (1979), 'Introduction to Geophysics, (Mantle, Core and Crust)', W.B. Saunders Co., Toronto. 199 Gauss, C.F., (1838), Allgemeine Theorie des Erdmagnetismus, reprinted in Werke, Band 5, pg. 121-193. Gough, D.I., (1973), Geophys. J. R. astr. Soc, 35 , pg. 83-98. 'The Interpretation of Magnetometer Array Studies' Gradshtyn, I.S. and Rhysik, I.M., (1965), 'Tables of Integrals, Series and Products, Academic Press, New York. Hensel, G., (1981), Private correspondence. Jacobs, J.A., (1970), 'Geomagnetic Micropulsations', Springer-Verlag, New York. Jones, F.W., (1970), Geophys. J. R. astr. Soc, 22 , pg. 17-28. 'Electromagnetic Induction in a Non-Horizontally Stratified Two-Layered Conductor' 200 Jones, F.W., (1971), Phys. Earth Planet. Inter., 4 , pg. 417-424. 'Electromagnetic Induction in a Two-Dimensional Model of an Asymmetric Two-Layered Conductor' Jones, F.W. and Pascoe, L.J., (1971), Geophys. J. R. astr. Soc, 24 , pg. 3-30. 'A General Computer Program to Determine the Perturbation of Alternating Electric Currents in a Two-Dimensional Model of a Region of Uniform Conductivity with an Embedded Inhomogeneity' Jones, F.W. and Price, A.T., (1970), Geophys. J. R. astr. Soc, 2_0 , pg. 317-334 . 'The Perturbations of Alternating Geomagnetic Fields by Conductivity Anomalies' Jupp, D.L.B. and Vozoff, K., (1977), Geophys. J. R. astr. Soc, 50 , pg. 333-352 'Two-Dimensional Magnetotelluric Inversion'. 201 Keller, G.V. and Frischnecht, F.C., (1966), 'Electrical Methods in Geophysical Prospecting', Pergamon Press, New York. Kisabeth, J.L., (1975), Phys. Earth Planet. Inter., J_0 , pg. 241-249. 'Substorm Fields In and Near the Auroral Zone' Kisabeth, J.L. and Rostoker, G., (1971), J. Geophys. Res., ]_6 , pg. 6815. 'Development of the Polar Electrojet During Polar Magnetic Substorms' Kisabeth, J.L. and Rostoker, G., (1977), Geophys. J. R. astr. Soc, 49 , pg. 655-674. 'Modelling of Three-Dimensional Current Systems Associated with Magnetospheric'Substorms' Kuckes, A.F., (1973), Geophys. J. R. astr. Soc, 32 , pg. 1 19-131. 'Relations Between Electrical Conductivity of a Mantle and Fluctuating Magnetic Fields' 202 Lahiri, B.N. and Price, A.T., (1939), Phil. Trans. Roy. Soc. London, A237 , pg. 509-540. ' 'Electromagnetic Induction in Non-Uniform Conductors and the Determination of the Conductivity of the Earth from Terrestrial Magnetic Variations' Landau, L.D. and Lifschitz, E.M., (i960), 'Electrodynamics of Continuous Media', Pergamon Press, New York. Law, L.K., Auld, D.R., and Booker, J.R., (1980), J. Geophys. Res., 85 , pg. 5297-5302. 'A Geomagnetic Variation Anomaly Coincident with the Cascade Volcanic Belt' Lee, K.H., Pridmore, D.F., and Morrison, H.F., (1981), Geophysics, 46 , pg. 796-805. 'A Hybrid Three-Dimensional Electromagnetic Modelling Scheme' 203-Levy, S. and Fullagar, P.K., (1981), Geophysics, £6 , pg. 1235-1 243. 'Reconstruction of a Sparse Spike Train from a Portion of its Spectrum and Application to High Resolution Deconvolut ion' Lighthill, M.J., (1958), 'Introduction to Fourier Analysis and Generalized Functions', Cambridge at the University Press. Lilley, F.E.M., ( 1 974 ), Phys. Earth Planet. Inter., 8 , pg.. 301-316. 'Analysis of the Geomagnetic Induction Tensor' Lilley, F.E.M., (1975), Phys. Earth 231-240. 'Magnetometer Array Studies: A of Observed Fields' Planet. Inter., H) , pg. Review of the Interpretation Lilley, F.E.M. and Bennett, D.J., (1973), Phys. Earth Planet. Inter., 7 pg. 9-14 'Linear Relationships in Geomagnetic Variation Studies' 204 McDonald, K.L., (1957), J. Geophys. Res., 62 , pg. 117-141. 'Penetration of the Geomagnetic Secular Field Through a Mantle with a Variable Conductivity' Madden, T.R. and Swift Jr., CM., (1969), In: 'The Earth's Crust and Upper Mantle', P.J. Hart (ed.), Am. Geophys. Monograph 13, pg.469. 'Magnetotelluric Studies of the Electrical Conductivity Structure of the Crust and Upper Mantle' Mareschal, M., (1981), to be published. ' Source Effects and the Interpretation of Geomagnetic Sounding Data at Sub-Auroral Latitudes' Nishida, A., (1978), 'Geomagnetic Diagnosis of the Magnetosphere', Springer-Verlag, New York. Oldenburg, D.W., (1969), M Sc. Thesis 'Separation of Magnetic Substorm Fields' 205 Oldenburg, D.W., (1976), Geophys. J. R. astr. Soc, £4 , pg. 413-431. 'Calculation of Fourier Transforms by the Backus-Gilbert Method' Oldenburg, D.W., (1979), Geophysics, 44 , pg. 1218-1244. 'One-Dimensional Inversion of Natural Source Magnetotelluric Observations' Oldenburg, D.W., (1981), Geophys. J. R. astr. Soc, (65 , pg. 359-394. 'Conductivity Structure of Oceanic Upper Mantle Beneath the Pacific Plate' Panofsky, W.K.H. and Phillips, M., (1962), 'Classical Electricity and Magnetism', Addison-Wesley Publ. Co., Reading, U.S.A. Parker, R.L., (1977), Ann. Rev. Earth Planet. Sci., 5 , pg. 35-64. 'Understanding Inverse Theory' 206 Parkinson, W.D., (1959), Geophys. J. R. astr. Soc, 2 , pg. 1-'Directions of Rapid Geomagnetic Fluctuations' Parkinson, W.D., (1962), Geophys. J. R. astr. Soc, 6 , pg. 441-449. 'The Influence of Continents and Oceans on Geomagnetic Variations' Perreault, P. and Akasofu, S-I., (1978), Geophys. J. R. astr. Soc . , 5_4 , pg . 547 . Porath, H., Oldenburg, D.W., and Gough, D.I., (1970) Geophys. J. R. astr. Soc, J_9 , pg. 237-260. 'Separation of Magnetic Variation Fields and Conductive Structures in the Western United States' Porath, H., Gough, D.I., and Camfield, P.A., (1971), Geophys. J. R. astr. Soc, 20 , pg. 387-398. 'Conductive Structures in the North-Western United States and South-West Canada' 207 Praus, 0., (1975), Phys. Earth Planet. Inter., _1_0 , pg. 262-270. 'Numerical and Analogue Modelling of Induction Effects in Laterally Non-uniform Conductors' Price, A.T., (1950), Quart. J. Mech. Appl. Math., 3 , pg. 385-410. 'Electromagnetic Induction in a Semi-Infinite Conductor with a Plane Boundary* Price, A.T., (1962), J. Geophys. Res.,, 67 , pg. 1907-1918. 'The Theory of Magneto-Telluric Methods when the Source Field is Considered' Rikitake, T. , (1950), Bull. Earthquake Res. Inst., Tokyo University, 28 , pg. 45-100, 219-262, 263-283. 'Electromagnetic Induction within the Earth and its Relation to the Electrical State of the Earth's Interior' Rostoker, G., ( 1972), Rev. Geophys. Space Phys., J_0 , pg. 157-211. 'Polar Magnetic Substorms' 208 Rostoker, G., (1978), J. Geomagn. Geoelectr., 30 , pg. 67-107. 'Electric and Magnetic Fields in the Ionosphere and Magnetosphere' Runcorn, S.K., (1955), Trans. Am. Geophys. Union, 3_6 , pg. 191-1 98. 'The Electrical Conductivity of the Earth's Mantle' Ryshik, I.M. and Gradstein, I.S., (1963), 'Tables of Series, Products, and Integrals', Veb Deutsher Verlag Der Wissenschaften, Berlin, 1. Samson J.C, (1977), Geophys. J. R. astr. Soc., 5J_ , pg. 583-603. 'Matrix and Stokes Vector Reperesentations of Detectors for Polarized Waveforms: Theory, with Some Applications to Teleseismic Waves' Schmucker, U., (1970), 'Anomalies of Geomagnetic Variations in the Southwestern United States', Monograph, Scripps Inst. Oceanography 209 Siebert, M. und Kertz, W., (1957), Zur Zerlegungeines Lokalen Erdmagnetishen Feldes in Aussern and Innern Anteil', Nachr. Akad. Wiss. Gottingen, II. Math-Physik. K1., pg. 87-112. Sokolknikoff, I.S., (1951), 'Tensor Analysis, Theory and Applications', John Wiley and Sons, Inc., New York. Strang, G., (1976), 'Linear Algebra and its Applications', Academic Press, New York. Telford, W.M., Geldart, L.P., Sheriff, R.E., Keys, D.A., (1976), 'Applied Geophysics', Cambridge University Press. Vestine, E.H., (1941), I.Terrest. Magnetism Atmospheric Elec, 46 , pg. 27-41. Vozoff, K., (1972), Geophysics, 37 , pg. 98-141. 'The Magnetotelluric Method in the Exploration of Sedimentary Basins' 210 Wait, J.R., ( 1 954), Geophysics, J_9 , pg. 281-289. 'The Relation Between Telluric Currents and the Earth's Magnetic Field' Wait, J.R., (1970), 'Electromagnetic Waves in Stratified Media', Pergamon Press, Toronto. Weaver, J.T., ( 1963), 3. Geophys., 29 ,'pg. 29-36. 'On the Separation of Local Geomagnetic Fields into External and Internal Parts' Weaver, J.T., (1973), Phys. Earth Planet. Inter., 7 , pg. 266-281 . 'Induction in a Layered Plane Earth by Uniform and Non-Uniform Source Fields' Wiggins, R.A., (1972), Rev. Geophys. Space Phys., H) , pg. 251-285. 'The General Linear Inverse Problem: Implication of Surface Waves and Free Oscillations for Earth Structure' 21 1 Woods, D.V., (1979), Ph.D. Thesis. 'Geomagnetic Depth Sounding Studies in Central Australia' 212 Appendix A Maxwell's Equations in a Conductor Maxwell's equations in their general form are: V-D = V- B = 0 VXE = - ££. (A.1 ) (A.2) (A.3) VxH = j + IB. (A.4) If we presume that all time dependences are of the form eLUjt and that 6,^. are equal to €0, f*o everywhere, then the equations simplify to: 213 V-D = (A.5) V-B = 0 (A.6) VxE = -iwB (A.7) V xH = J. + iwD (A.8) Employing the constitutive relations for an isotropic medium: 5* = crE (A.9) D = €QE (A.10) B = fi0H (A.11) 214 we arrive at the starting form of Maxwell's equations for all problems in this thesis: V-E = _/>/€o (A.12) V-H = 0 (A.13) VxE = 1WI (A.14) V xH = ( <y + iwto ) E (A.15) All induction problems that will be dealt with in this thesis will consist of a non-conductive half-space (air) and a conductive half-space (earth) with the source in the non-conductive half-space at a distance from the boundary plane (as shown in Fig. 1.5). The Z axis will always be perpendicular to the boundary, with the positive direction downward. Within the earth, the lowest conductivity expected ( ~ 10 5 S/m) is still much greater than the value of w€0 corresponding v 215 to the highest frequency used in G.D.S. (the frequency range is shown in Table 2.1): = 10~5 S/m » w^K€0 = 5.6x10~lZ S/m (A.16) Thus, within the earth we can always neglect the displacement current term, so that A.15 becomes: VxH = GE for z > 0 (A.17) Taking the curl of both sides of equation A.17 and utilizing equation A.13 on the left hand side of the subsequent equation, we arrive at: VlH = -7t7 xE + cr ( $ xE) (A.18) Upon replacing E from A.17, and VxE from A.14, we have: v"*H = -^pX(vxH) + iwju^H (A.19) This is the general equation for the magnetic field in an inhomogenous earth. The corresponding result for the electric field is obtained in analogous fashion by taking the curl of both sides of A.14 and then simplifying with A.12 and A.17: 216 v^E = iw^ffE + (V/> )/€0 (A.20) Going back to equation A.15, taking the divergence of both sides (and noting that the divergence of the curl of a vector is always zero), it will always be true that: "V • [ ( Cr + iw€0 )E] = 0 (A.21) Thus, by application of the divergence theorem, the vertical component of the vector ( O" + iwde)E- is continuous across the boundary plane separating the halfspaces, so that at the boundary: (OkJr + iw€o)E4(0;^= (CWn, + iw€0 )EMca,^ (A.22) Using A.16 and the zero conductivity of the air, we have: (A.23) Again from A.16, it is apparent that at the boundary plane of the two half-spaces, the vertical electric field in the earth and consequently the vertical currents, will be insignificant. If the conductivity varies only with depth after the boundary, then by symmetry considerations, E-j and Jt in the earth will always be insignificant, regardless of the source type. Considering equation A.21 again, expanding the terms, 217 utilising A.12, and finally converting ' iw' back to ^ , leads to: 1|. -,5a-E, (A.24) If O" varies only in the Z direction, and Ej from the previous discussion is insignificant in the earth, then the source term for charge creation on the right hand side of A.24 is negligible, leading to the equation governing the rate of decay of existing free charge: = 0 (A.25) The solution of this is: (A.26) Even for the most resistive rock ( O" = 10~? S/m) the half-life for free charge is only ~10 s. Effectively then, there will be no free charge in the earth if cr = cr(z) only (except at the boundary ). As well there can be no free charge in the non conducting half-space. Thus, everywhere except at the boundary: V' E = 0 (A.27) 218 Simplification of A.20 is now possible because of the lack of free charge, leaving: V1 E = <p E (A.28) where: 9 = -y>^aafic (A.29) z < 0 (air) 9 = iw /u0cr (A.30) z > 0 (earth) Again, these results are for the case where G" is a function of depth only. If we now further stipulate that the variation of cr with depth occurs in a layered fashion, with <3~(z) constant within each layer (giving VO~ = 0 within each layer), then A. 19 simplifies to: (A.31) within the n^ layer, where 9*. is the value of ^ in that layer. 219 Appendix B Correlation of Price's Induction Solutions with Plane Wave As seen in section 2.1 of Chapter II, the required non-divergence of the electric field in the homogenous earth induction problem can be satisfied in two ways, leading to two independent types of solutions. Solutions of the first type will be correlated here with waves with a transverse electric field (TE) and the second type will be correlated with waves with a transverse magnetic field (TM). Type one solutions The first type solutions for the electric field are given by equations 2.1.20 and 2.1.21 in Chapter II: Solutions E . (x , y , z ) = -CS>* + qO"x* + B (B.1 ) z < 0 220 E,<x,y,z) = [C, e-^l+^V'lai (B.2) z > 0 with P(x,y) satisfying (from 2.1.19): + ILZ. + S^P = 0 "ax1 ^VJ'-(B.3) and with A,,B,, and C4 satisfying (from 2.1.24, 2.1.25, and 2. 1 .26): B, = -A|•[(1-R)/(1+R)] C, = A ,[2R/( 1+R) ] where: R = [( ^ +?a>/( ^ +9e)] The elementary solution to B.3 is: p(x,y) = eik** eU3* with: (B.4) (B.5) (B.6) (B.7) 221 kx 2 + kyZ = (B.8) Lett ing: kze = -i ( + <j>e j''2" (B.9) and: kza = -i ( 0* + <Va. (B.10) and using B.7, then the elementary solutions B.1 and B.2 take on standard plane wave format (where it is noted that the time dependence has been resurrected): E,(x,y,z) = i{ky,-kx,0} [A, e~a**e ^V^U eiu,t + B, eiki^ eCk*K elk3U iult (B.11) z < 0 ^/ \ • (, , «•> „ kj_^ <-kKX ik¥w coot E,(x,y,z) = i{ky,-kx,0} C, e K e e He (B.12) z > 0 To simplify these results we rotate by an angle -o\ about the Z axis to a new horizontal coordinate system ( f , Y) ) (see 222 Fig. B.1) such that the electric vector is only along |? , and horizontal propogation is only in the direction. The equation for the electric field becomes: E,(^ , r> f2) «i{S> ,0,0) [A, e^2 e'^eCwt + B, e *°- e L e z < 0 (B.13) E,(^,Y£,z) = i{"9,0,0} C| e~ck*V^ eiujt (B.14) z > 0 In terms of the incident angle 9" and the transmitted angle ©' (where the angles may be complex): ka sin© = \) (B.15) ka cos & = kza (B.16) for z < 0 where: 223 Fig. B.1 Rotating by an angle - o< about the Z axis to a new horizontal coordinate system. The horizontal direction of propagation is now purely in the r? direction. 224 ka = ( S)1- + kza* ) is the total wavenumber in the air. Also: ke sin = V) (B.17) (B.18) i ke cos c? = kze for z > 0 where: (B.19) ke = ( V*1" + kze* ) (B.20) is the total wavenumber in the earth. The equivalence of in the earth and the air, combined with B.15 and B.18, results in Snell's law (Panofsky and Phillips, 1962, pg.196-197): ka/ke = sin ©/sin © (B.21) Also, equating the first term in B.13 with the incident wave amplitude Ei, and the second term with the reflected wave amplitude Er, we have: 225 Er/Ei - B,/A, = [(-^i)"1- i] (B.22) Upon substitution using B.9, B.10, B.16 and B.19, we arrive at the usual Fresnel relation for the electric field for the TE mode (Panofsky and Phillips, 1962, pg.198): Er/Ei = [cos© - (ke/ka)cos©']/[cos© + (ke/ka)cos ©'] (B.23) In similar fashion we have for the transmitted wave amplitude, Et: Et/Ei = C,/At - 2'( )''V[1 + (*l±*S^)'i] (B.24) which after substitution becomes the second Fresnel relation for the TE mode (Panofsky and Phillips, 1962, pg.198): Et/Ei = 2 cos B/[ (ke/ka)cos©' + cos© ] (B.25) Type two solutions The second type solutions for the magnetic field are given by equations 2.1.39 and 2.1.40 in Chapter II: Hi(Xry.z) = (i/wp0) ^ {1E*,-Mi,,0} [A 226 (B.26) for z < o st<*.y.*> • <i/v.)s^3Ilitia..-2a.o) c,.-"1^"* (B.27) for z > 0 with FE(x,y,z) satisfying (from 2.1.38): ^1 * and with A2,BZ,C2 satisfying (from 2.1.41 and 2.1.42) To. Ya. q = A2-2 /[i + (-J5-)!*] Again, the elementary solution to B.28 is: Fa(x,y) = eik*K eikSl with: (B.28) (B.29) (B.30) (B.31 ) 227 kx2, + ky* = ^ (B.32) Using the definitions for kze and kza as in B.9 and B.10, and substituting the result of B.31 into B.26 and B.27 we obtain the standard plane wave format for the magnetic field in the type two case: Hz(x,y,z) = -(i/w«0) -7^— {ky,-kx,0} >[hze e l-'e e ~Bze e 33e e ] (B.33) for z < 0 Ht(x,y,z) = -(i/W/i0) -j^- {ky,-kx,0} Cze:t"x e'^e ik^1 e 1 wt (B.34) for z > 0 Rotating into the new coordinate system ), in identical fashion to that done for the type one solution, the magnetic field becomes: 228 •B2e fe e ] for z < 0 (B.35) Hz(f ,7 ,z) = -(i/w^--^-{S>,0,0} qe^V^V^ (B.36) for z > 0 • Following the same procedure as for the type one solution, Snell's law is again recovered, and also the Fresnel relations for the TM mode (Panofsky and Phillips, 1962, pg.198): Hr/Hi = -Ba/A^ = [(ke/ka)cos© - cos ©' ]/[(ke/ka) cos © + cos6 ] (B.37) Ht/Hi = Cz/hi = [2(ke/ka)cos6 ]/[ (ke/ka)cos© + cosS'] (B.38) 229 Appendix C The Uniform Field Assumption A common assumption in both G.D.S. and magnetotellurics is that the electromagnetic plane waves comprising the primary field can all be considered to be propagating vertically downward, so that the field is horizontally uniform. This assumption inherently implies that all waves of significant amplitude in the source field must adhere to three conditions: (1) The angle from the vertical of the direction of propagation of the wave transmitted into the earth is near zero. (2) The horizontal wavelength of the transmitted wave is much greater than the lateral extent of any anomaly, so that there will be little horizontal variation in the inducing field across the anomaly. (3) The complex magnitude of the transmitted wave (for a normalized incident wave) is the same for all waves. 230 From equation 2.1.51 or 2.1.52, and using equation 2.1.57, the downgoing electromagnetic wave in the earth can be expressed as: -» - i\>e -ifeui cut F = F0 e i e e (C. 1 ) where F is the total horizontal distance and F can be either the electric or magnetic vector. Equation C.1 can be rewritten as a product of decaying and propogating parts of the wave: (C.2) Note that from equation 2.1.59 the imaginary part of k %e is always negative, so that e I dogs -n facfc represent A decay with increasing depth. The propogation in the P direction is represented by eL 5 , and that in the z direction by -tRccJCk^H , _ e , so that the angle of propogation, measured from the vertical is given by: 9 = Tan"* [v»/-Real(kze) ] = Tan"1 [V/imag( 1 + ip) ] (C.3) with the definition of "p the same as in equation 2.1.35. The dependence of B on p is plotted in Fig.C.1. It is seen that for p > 100, all waves propogate virtually vertically downward. Using the expected physical range of given in Tables 2.2 and 2.3, and interpreting for a midrange conductivity of .05 S/m, it 231 a rsi. (N CD UJ Q a .to CLZ~ I— LJ X a oo a a 1 .0 ^ r 3.0 5.0 LOG rBETfl) 7.0 9.0 The propagation angle from the vertical of the transmitted wave into a half-space of uniform conductivity. The angle is plotted as a function of the dimensionless parameter, jg . 232 is seen that for large values of the horizontal wavelength ( ?\ > 3 5x10 km.) this magnitude of B is attained for all periods less than ~ 2 hours. Thus, within these parameter boundaries the first condition will be satisfied. As well, the large horizontal wavelengths required will certainly guarantee that the second condition is also satisfied. The remaining condition, if met, will ensure that the magnitude of the anomalous response will be the same for all significant waves in the source field. The relations between the complex magnitudes of the incident ' (A), reflected (B), and transmitted (C) waves are defined in equations 2.1.24 - 2.1.30, with the results being: B = f ( B )• A (C.4) C = A + B = [ 1 + f(p ) ] A (C.5) The amplitude of f(B ) is plotted in Fig.2.la. Between B = 100 and B-»oo , |f(p> )| varies by a factor of 2, so that condition (3) is not satisfied for this range of p . However, in some methods of analysis used in G.D.S., which will be discussed in Chapter III, one is concerned only with the relation between the 'normal surface field, and the addition to the 'normal' surface 233 field due to an interior anomalous region,(that is, the anomalous field). The 'normal' field is the total surface field that would exist in the absence of any anomalous variations from a one dimensional conductivity structure. In this case only the ratio of C to (A + B) must remain constant for all waves to act alike, and this is always guaranteed by the continuity of the magnetic field across a boundary, as expressed in equation C.4. Thus, in this case the third condition is not required, so that the uniform field assumption is valid for the parameter ranges satisfying the first two conditions. It should be emphasized that this works only for methods using the relation between the total 'normal' surface field and the anomalous surface field. It will not be true for methods which relate the magnitude of the incident field to the anomalous field, as that will require that condition (3) be fulfilled. 234 Appendix D Separation of the External and Internal Fields In illustration of the more general separation formulas derived by Siebert and Kertz (1957) and Weaver (1963), the simplified case of a two dimensional earth will be considered. Presume that neither the inducing field nor the earth structure varies in the 'y' direction (see Fig. 4.1), and also that the station array is along" the surface at right angles to the y axis. Thus, the magnetic field will only have directional components X£,Xj and ZC,ZT in the 'x' and 'z' directions respectively. The subscripts 'E' and 'I' will denote whether the field component is of external or internal origin. Let a single external line current, I, at the position x=0,z=-H, be the source of the fields at the surface. The measured field at any array position x, from Ampere's law will be: Xc(x) = (/i0I/2Tr)- [H/(Hl+x*)] (D.1 ) ZE(x) = -(/i0I/2Tr) [x/(Hl+x1)] 235 We will need to know the general (Ryshik and Gradstein, 1963, pg. } [a/(x*+al>] =TTe"^la (D.2) Fourier transform expressions 250) : (D.3) }[x/(xl+a1)] = -frie~^lasgn(^) (D.4) where the Fourier transform, J , of a function of x, f(x), is defined here to be: oo } [f (x)] = |*f (x)e'^X (D.5) and the 'sgn' function is 1 or -1 depending on the sign of its argument. Using these, the spatial Fourier transform of the measured fields along the surface will be: J[XE(x)] = (jUeI/2fr) [Tre~'flH] (D.6) }[ZE(x)] = (ji0I/2n) [TTeWsgnflj)] Thus: (D.7) 236 }[ZE(x)]/ J[XE(x)] = isgn(^) = jf [-1/(tTx)] (D.8) (from Lighthill, 1958, pg. 43). We can rewrite equation D.8: 3 [ZE(x)] = J [-l/(tfx)] J [XE(x)] (D.9) Then, using the Faltung, or convolution theorem,(Bracewell, 1965, pg.25) we can rewrite the relation between XE and ZE in the frequency domain as a convolution in the time domain: CO ZE(x) = f XE(x*)•1/[TT(x-x')]dx' -co (D.10) The convolution with the function [-1/(TTx)] is the Hilbert transform (Bracewell, 1965, pg. 267). Denoting this operation by the symbol 'K', D.10 becomes: ZE(x) = K[XE(x) ] (D.11) Noting that the reciprocal of sgn(^) is still sgn(^), we could also have rewritten equation D.8 in the form: J rx£(x)] = which results in: -isgn(^) J [ZE(x)] (D.12) 237 XE(x) = -K[Z£(x)] (D.13) The linearity of Ampere's law allows us to extend this result for one external line current, to the general case of an arbitrary number and distribution of external line currents. Consider now that the source is a line current at a depth H below the surface. From Ampere's law the field X j, Z-^ at the surface will be: XjU) = -(/i0I/2tr) [H/(x*+Hx)] (D.14) Zj(x) = -(^0I/21T) [x/(x1+H1)] (D.15) Following the same arguments as before, these equations lead to: XI(x) = K[Zx(x)] (D.16) Z1(x) = -K[Xx(x)] (D.17) Again, this result can be extended to an arbitrary distribution of internal currents. 238 The total fields measured at the surface will be: X(x) = X-j-(x) + XE(x) (D.18) Z(x) = Zx(x) + ZE(x) (D.19) Taking the Hilbert transform of the measured total data gives: K[x(x)] = K[XJ(X)] + K[XE(x)] = -Zx(x) + ZE(x) (D.20) and: K[Z(x)] = K[Zx(x)] +K[ZE(x)] = XX(x) - X£(x) (D.21) Thus, the portion of the measured field that is due to internal sources can be separated: XT(x) = (K[Z(x)] +X(x)}/2 (D.22) 239 Zj-U) = (Z(x) - K[X(x)]}/2 (D.23) It should be noted that the Hilbert transform of a constant gives a zero result/ so that in the case of a constant field of either internal or external origin, equations D.22 and D.23 will not be able to separate the components, and will merely divide the total field into equal parts of each component. 240 Appendix E Determination of the Induction Tensor Elements The following method of obtaining values for the induction tensor elements is due to Schmucker (1970). Other methods have been suggested by Everett and Hyndman (1967), and Woods (1979, pg.53 ) . The anomalous and normal fields at any frequency 'w' are considered related by the induction tensor, but with the possibility that a portion of the anomalous field cannot be correlated with the field: / CHH r DA = Cl>D ZA / Vc*u H, (E. 1 ) The values of the tensor elements are complex, with the real and imaginary parts represented here by superscript I or R as shown: C = CR + irx HH HH H H (E.2) The auto and cross powers of two signals a(t) and b(t) of equal length Tc will be defined as: 241 SflB = A(w)B*(w)/T0 (E.3) where A(w) and B(w) are the Fourier transforms of a(t) and b(t) respectively, and the star '*' denotes the complex conjugate. The desired values of the tensor are those which minimize the uncorrelated terms , Dr and . Consider the last anomalous vector component: zfl " C*HHN + CfcD dN + ZN +Zf (E.4) The values of Zc and Zj are then: and: Z£ = "(CJH + ic\H).HN - (C.J0 + iC^)-D -(C^ + iC^ )• Zw + zA Z* = - icJH).Hj - (C^ -icJD).DJ (E.5) (C4% - iC^)-ZN + ZA (E.6)' where we have expanded the tensor terms into both real and imaginary terms. Z^, and thus the autopower of Zr, is dependent on six independent variables, as seen in E.6, so that the minimum must occur when the derivative of is zero with respect to each of them: 242 ^ii£lL = 0 (E.7) (E.8) S^L?£ = 0 ; *Sit}f = 0 (E.9) Using the chain rule (where 'u' is any of the independent variables) we have: tlmL= (Zr/TC)^L+ (Zr/Tc)l^-} U, (E.10) Evaluating this using equations E.5 and E.6 for each of the cases, six separate equations are obtained: (E.11 ) "Z*HM + 1S H* = 0 (E.12) 243 4DN + VD* = 0 (E.13) -Z*DK + ZfD* = 0 (E.14) zJzN + Zf Z* = 0 (E.15) -Z*ZN + ZfZl = 0 (E.16) For each of the pairs (E.11 and E.12), (E.13 and E.14), and (E.15 and E.16) the only possible solutions are either that HN,DN, and ZN are zero, which maximizes Z^, or that: 4H*= zrH*= 0 'H** ° (E.17) Z*Dw = zrOj = 0 , Dw * 0 (E.18) 244 Z£ZM - % Z * « 0 , ZN * 0 (E.19) This is in fact a restatement of the original problem, as each of E.17, E.18 and E.19 is merely stating that Z is uncorrelated with the normal field, with: SH„^ = o ; s^Hw = o (E.20) (E.21) Inserting the expressions for lg each of the above equations, equations: C*tfSHioHw + C*DSDN HN + C2* (E.22) if. and Zc from E.5 and E.6 into we arrive at the set of linear 2N H* = Stk Hw (E.23) 245 (E.24) (E.25) which allow us to calculate the values of the three induction tensor elements: f 2H (E.26) where S is given by: S = SHMHN ST>MH- S2«/HW u 2w DM Sr\*\„ SD«iN Siw?k), (E.27) The values of the induction tensor elements for the other two rows would proceed in identical fashion. 246 Appendix F Properties of Current Distributions that Mimic a Line Current Presume that a particular two-dimensional current distribution that fits the observed surface readings Bx,By along a linear array at a frequency w, is a line current at a depth zQ, as in Fig. 3.7. The ability to fit the data with this line current model then imposes certain constraints on all other possible models that fit the data. As many actual anomalies have surface magnetic fields that closely resemble those of a single line current, these model constraints are of obvious interest. Consider a current distribution, j(x,z), which is non-zero between z=0, and z=co. For this distribution to duplicate the surface field due to the line current, the values of Bx and Bz must match at all surface positions 'x' along the array. Because of the linear relationship between the two components (as shown in Appendix D), it is sufficient to consider only one of the components. Using the Bx component, a suitable model j(x,z) will thus satisfy: 247 (Ue/2fr) lz6/(xl+ze1-) OO OD = (f./2tt)^ {[ j(u,v) -v]/[ (x-u)z+v*]}dudv (F.1) Noting the convolutional form of the right hand side we can rewrite equation F.1 as: CX) Iz0/(xl+zl0) = ^ [j(x,v) ® (v/(xi+v"1))] dv (F.2) Taking the Fourier transform with respect to 'x' of both sides of F.2 gives: ITT e"'!Uo = ^ [J(^,v) TT e 3 ] dv (F.3) where: CO -tea du J(<£,v) = ^ j(u,v) e~L^ -00 Equation F.3 simplifies to (F.4) oo Ie lV*' = ^ j(^v) e dv (F.5) Consider now a model that has no currents at, or shallower than the line current depth, that is, j(x,z) is non-zero only between z0+€ and oo , where € is an arbitrarily small number. Using this, moving e J ° to the right hand side in F.5, and 248 introducing a change of variables, v = (v-z), equation F.5 becomes: co I e (F.6) The upper limit of integration in F.7 must actually be finite because of the finite 'depth' of the earth; denote this finite limit by 'R'. Now, by the mean value theorem (Fulks, 1969, pg.125) a value of the integrand will be able to be found at some position v within the range of integration such that: This must hold for all values of ^, which means that in the limit as ^ goes to co , J(^,v) must increase monotonically to infinity. Thus, the requirement that the entire current distribution be at a greater depth than that of the line current results in physically untenable models. In other words, all realizable models that fit the data must have some current at, or shallower than the line current depth. Return to the general problem again, where the current distribution that will mimic the line current can be non-zero anywhere between z = 0 and z= co. The required relation between this current distribution and the line current is then given by equation F.5. If the current distribution J(^,v) has elements which are not in phase with the original line current, then that -iriu-v) e 5 (F.7) 249 model will have both real and imaginary portions, JR(^,v), Jj(^,v) respectively. Designating the line current as the zero phase position, then I will be pure real. This leads to equation F.5 becoming: ^l0= I JR(^V) e^dv ° (F.8) and , CD 0 o = $ Jx(^,v) e~l5,trdv (F.9) In the case of F.9, as e ^'^is always greater than zero, the neccessary zero value of the total integral requires that either Jr(^,v) is everywhere zero, or that Jj.(^,v) has elements of both negative and positive sign. If we now restrict our range of possible models to include only those where the current elements in the model are different in phase by less than TV/2, then the only way to satisfy equation F.9 will be to have Jx(^,v) zero everywhere. Thus, applying this constraint forces all current elements to be in phase with the original line current. Although this precludes the possibility of 'annihilator' distributions in the model, it does not preclude currents at depths greater than that of the line current. An example of such a distribution is one consisting of a line current of magnitude I0/4 at a depth 2z0, plus a lineal current density: 250 j(x) = I0{(z0-2l )/[xl + (ze-z, r" ] - (2ze-z, )/4[x2+(2z0-z, f ]} (F.10) along a depth z,, where z, < zQ. In a case such as that given, where all the currents are in phase and the current distribution contains only one line current I at a depth z, greater than zc, then it must be true that: zoxo ^ Mi (F.11) This is required to ensure that the X component of the original line current field is always greater than or equal to the X component of any portion of the mimicking current distribution at asymptotically large values of x. 251 Appendix G Analytic and Numeric Inteqrations In section 2 of Chapter IV the inner products of the current density kernels are to be computed. The kernels correspond to the contribution to the surface magnetic field at station position (x,0), due to a current density at the position (x',z'). There will be two types of kernels, corresponding to the vertical (z) and horizontal (x) components of the field (as given in equations 4.1.6 and 4.1.8): Gx(x,x',z*) = /V2TT {-z'/[ (x-x')2+z*2-]} (G.1) Ga(x,x',z') = /V2ir {(x'-x)/[(x-x' f+z'1]} ) (G.2) This allows for six types of inner product integrals: ' [jk = ^Gx (xj ,x' ,z' ) Gx(xk,x',z') dx'dz' S 252 /*o%rra ^ z'2- dx'dz' [ (X. -X' )l+z' X] [ (X. -X' J^+Z'1"] s J k (G.3) with either x^ H x^, or XJ = x^. = KG*(xj ,x' ,z' } G4(xk,x',z') dx'dz S ^oV4TTl ff (x'-xk) (x'-x; )y I (XJ -x' )l+z- v -x.i ) dx'dz'  ] [ (xt-x« )* + z«*] s (G.4) with either x; * x, , or x: = x, . Ijk = JjG4(xk ,x' ,z* ) GK(XJ ,x' ,z') dx'dz' 5 = ^v4Tri rr -z'-(x'-xk) dx'dz'  JJ [ (XJ -x' )«-+z'»- ] [ (xk-x ' )*-+z ' *• ] s (G.5) with either xj \ x^, or XJ = x^. The range S will be designated A < x' < B; €: < z' < D, where 6 must always be greater than zero. Using standard integrating techniques, the first integration of these inner products can be carried out with respect to either of the variables, resulting in the following analytic forms. 253 First integration with respect to z : r]k = ^o/4ir* $ dx' {(x'-xk)/[ (x,-xk)i-(x'-xj )z ] A -{Tan"1 [D/(x'-xk) ]-Tan~l [€./( x '-xk) ]} (x'-Xj )/[ (x'-xj )1-(x'-xk)1 ] • {Tan-1 [D/(x'-xj ) ]-Tan"' [t/(x'-xj )]}} (G.6) x: * xu 3 • Pjk = Po/srr1^ dx' { i/(x'-xk) * ~. -I .[Tan (D/(x'-xk)) - Tan ' (€/(x'-xk))] + (€/[ (x'-xk)Z + - D/[ (x'-xk)i+D1])} (G.7) Xj = xk B r$k = poV4ftr _f dx' {(x'-xj )/[ (X'-XJ ^-(X'-XL.)2-] * {Tan"1 [D/(x'-xk) ]-Tan"1 [€/(x'-xk)]} (x'-xk)/[ (X'-XJ )i-(x'-xk)2 ] 254 -i •{Tan [€/(x'-xj ) ]-Tan 1 D^/( x '-XJ ) ]}} (G.8) X; * x. Pjk = ^ol/8Trl £ dx' { 1/(x'-xk) [Tan'1 (D/(x'-xk)) - Tan"' (V(x'-xK))] + (D/[ (x'-xk)2+D ] - €/[ (x'-xk)* + £*])} (G.9) xj = xk 3P rf/8T['L § (x'-xk)/[(x,-x^.)i-(x'-xk)z] * ln{[ (x'-xj )* +6xir (x'-xO* +D*-])dx' (J (x'-x; )l+DM[ (x'-xk)1 + 6*]} (G.10) Xj * xk f]k = MoV8irl 5" (x'-xk)1 dx • d/[ (x'-xk)1+€l] - l/[ (x'-xk)l+D1]} (G. 1 1 ) xj = xk In all cases, the results for XJ \ xk can be shown to have as their limiting values the XJ = xk expressions, when x: xk« 255 First integration with respect to x 1 £al f z'ldz' . lni 4tttJ 8(z*1- +d'"). d I [ (XK -A)1 +z 1 ][(xi -B )*" +z [<xK -B)t+z M[(x- -A)1 +z 2 Tan' / [ (x«-A) (XJ -A)-Z'"M z' ^ [z'(xK-A)+z'(XJ-A)] -Tan (XK-B) (xi -B)-z'* 1 z* (xk-B)+z' (XJ -B) ] XJ * xk (G.12) P; Jk £l_ C dz'f (XK 8rrv J l[z't + € -A) (XK-B) (xk-A)1 ] [ z '*• + (x, -B)z ] + J_ Tan z' -i |(x^-A)j - i_Tan"' |(xR-B)jj (G.13) n hi 4fT J_ In 8d [(x;-A)x+Z'*][(x*-B)>2'M [(xk-A)l+z'1][(xj-B)z+z,x] 256 + J_ 4z' Tan"' /[ (xk-A) (xj -A)-z'1 ] ([z'(xk-A)+z'(x--A)] - Tan"1 [ (xn-B) (xj-B)-z'1 1 [z' (xk-B)+z' (X>i -B) ] 1 /d In 8(z'1 +d* [ (xk-A)z+z'M [JxizJillz^.]1 [(xk-B)l+z*s-][(xj-A)-'-+z,;L] 2d_* z' Tan"1 / [ (xk-A) (xj -Al-z^ ] 1 [z' (xK-A)+z' (XJ -A) ] Tan"1 / [ (XK-[z XK-B) (xi -B)-z',• 1 )]\\ ' (xK-B)+z' (XJ -B) )J\ j (G.14) dz' 1 Tan (xk-A1 - (xk-A) |z' z' [z"- + (xk-A)'- ] JI_Tan"' IXK-BI I Z7 (xk-B) , [z,J- + (xk-B)l]J (G.15) 257 r, Ho 4TP ' fz'dz'f^ [Tan"' ([ (x*-A-d)Z -d*+z'z]^ J [4z'd L \ 2zTd J - Tan"1 ([ (xK-B-df -d 2z'd ^)] d C_l In f[ (XK-A^+Z'1 ][ (xo -B)*+z,a" ^ 8(z"-+dx) \ d \[(xk-B)l+z,i][(xj-A)1'+z,J-]y + 2_ z' Tan -I /[ (X'K-A) (xi -A)-z'r 3 Uz- (x. -k A)+z' (Xj-A) ] - Tan /[ (xk-B) (xj -B)-z't- ] \1) Uz' (xk-B)+z' (x- -B) (G.16) xj * xk 8vr [ (x^-B^+z'1 ] [ (xk-A)l+z"-] (G.17) 258 d = <xk-xj)/2 (G.18) The integration with respect to the final variable for either of the two possibilities must be done numerically, due to the 'Tan-1 ' and 'In' terms in the integrand. The method of numerical integration that was used involved sampling the integrand function between the limits of integration, and then fitting a cubic spline to these values. The cubic spline estimates the first, second and third derivatives of the integrand curve at every sampled point, so that the integral can be readily estimated: The sections of the integrand curve which changed rapidly, and thus were difficult to spline, were located and sampled on a much finer basis. The integrated results of the splining technique were compared with results using various U.B.C. computer library integrating routines, and were always within .0001. Also, the numerical integration was done using the integrands from both the first integration with respect to 'x' and with respect to 'z', with the results agreeing within .0001. The particular utility of integrating with respect to 'x' first, is that it allows the simple introduction of the weighting f(x')dx' = Df(x) + Daf ' (x) + D3- f ' ' (x ) + D*- f ' ' ' ( x ) 2! 3! 4! (G.19) 259 factor, z , into the integrand for the weighted smallest model construction of Chapter IV. The second result needed for the calculation of the averaging functions in Section 2 of Chapter IV is the integral of each of the kernels over the region of interest. There are only two possibilities, and both have analytic solutions: B T> 1 Uk = ^ GK Uk,x' ,z' ) dx'dz' A £ = f(xk-B). ln/[D*+(B-XKf ]-\ 2TT{ 2 Ufel+(B-xk)1 ]/ (XK-A), Inf [DHtA-x.,)1 ] 2 \[€z-+(A-xk)i ] D'Tan"1 /XR-B\ - €.Tan~' D« Tan"1 (G.20) B T> Gs(xk,x',z') dx'dz* A € 260 t 2TT - (l [ D IT (.2 . In \PX +a V - £ In /ilV €»• +al/J where; + a [Tan"' (6/a) - Tan"1 (D/a)] - b [Tan"' (€/b) - Tan'1 (D/b) ] ^ (G.21 ) xk- A (G.22) b = xk- B (G.23) For both the forward modelling routine and the linear programming construction, the contribution to the surface field from a rectangle of constant current density is required (see Fig. G.1). The contributions for each component at a station posit ion 1 Xj_ ' are; -z'dx'dz' [ (x, -x' )*+z"- ] = (/X./2TT) Jjkl^ (G.24) 261 ^-Surface J., J2i Parameterized current density model. The current density is constant with respect to position within each grid element. 262 BzU,; ,0) = ^k 2rr (x'-xi) dx'dz' J [ (X; -X* )*• + Z"-] = (/lo/21t ) JA Ij'k The values of I:. and I;, are: Jk Jk (G.25) \lk Tan -xj - Tan Xt_" z kti Tan xc -x '• k + t - Tan -i xc k-H Xi-X; ln [ZK1+ (XC -x-,'-)2- ] + (x,-xj) ] k+( xc. -x Jtl In [zkli + (xi -x;Vi )* 3 [z- + (x£-xjt,)M (G.26) <xt~xj+iHTanY zUi ]- .Tan"'/ zV \] \x. -x' ' Iv -V • I XC X jVly 263 (x- -X-' )[Tan" \X; -X: + zW, . In /[ (xt-x'jV,)* +z&3\ 2 ^[(xt-x] V+z^]/ + z_k . In A (xi-xl )» +zlx] \ 2 \Kxt-XjV, )l+z^]J (G.27) 264 Appendix H Determining the Major Axis of an Ellipsoid Let the linear transformation of a three dimensional vector x be represented by a real valued matrix A: y = Ax (H.1) The square of the length of the new vector, y, is: yTy = (Ax)TAx = xT(ATA)x = x Bx (H.2) The matrix B will be real and symmetric. The endpoints of all the position vectors x for which the squared length of y is a constant will map out a quadric surface (Sokolnikoff, 1951; Strang, 1976), with the governing equation being: yTy = C 265 = (x, ,xz,xz) / b„ ba bn\ / x, bi, b27. b2* j I xj, a 3 xb" b" bV \XV (H.3) For C > 0 the quadric surface is an ellipsoid. From equation H.2, the matrix JB will be positive definite as well as symmetric, which ensures that we can always decompose B: B = SAST (H.4) where S is an orthogonal matrix: s- = sT (H.5) and A is a diagonal matrix, wiilv. all values of the diagonal greater than zero: A (H.6) Thus, equation H.3 becomes: C = xT(SAST) x 266 = (xTS)-A'(STx) (H.7) Defining a new coordinate system by: I - sT i (H.8) (where the length of £ is guaranteed to be the same as that of x by the orthogonality of S) then equation H.7 reduces to: (H.9) Thus, in the rotated system, the equation of the ellipsoid becomes: c = + + (H.10) so that the unit vectors , ^ lie along the principal axis of the ellipsoid. The maximum lengths of the ellipsoid in each of the axis directions (using the fact that }V , Alfand are all greater than zero) are: ^( max = (C/ >vj) 1 (H.11) 267 ^m** = (C/\)'/2-(H.12) ^mA.K= (C/^)''* (H.13) Thus, a measure of the relative lengths of the principal axis of the ellipsoid is given by the ratios of the square roots of the eigenvalues: (H.14) ^•j max ^ i (H.15) As well, the directions of the principal axis of the ellipsoid in the original space are given by the rows (or eigenvectors) of S: (H.16) 268 (H.17) ^•i _ ^SlS 'S*3 'S33 ) (H.18) where the orthogonality of S ensures the unit length of the eigenvectors. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items