Finite difference time domain modeling of cross-hole electromagnetic survey data

Cross-borehole electromagnetic (EM) surveying is an emerging technology in the mining industry that can offer resolution of geoelectric structures superior to that from other EM methods. However, existing interpretation techniques often employ physical models of electromagnetic wave propagation that are inappropriate for the instrument parameters and electrical properties commonly encountered in the metalliferous mining environment. Lacking appropriate numerical models, we are unable to properly assess imaging methods or to design and interpret surveys. In this thesis, I address such issues relating to the delineation of massive sulfide mineralization using cross-borehole EM techniques. Specifically, I discuss the use of the finite-difference time-domain (FDTD) method in modeling the cross-borehole electromagnetic response of a perfect electric conductor (PEC), which is a good approximation to a target of high conductivity. This model type has particular application to massive nickel sulphide ore bodies, such as those found in the Kambalda, Sudbury Basin, Norilsk, and Victoria Island regions. A key development of the thesis is the implementation of perfectly matched layer (PML) absorbing boundary conditions in a Yee FDTD algorithm. This revised algorithm permits the realistic modeling of the electromagnetic response of conductors in the frequency range transitional between the diffusion and wave domains. Such modeling reveals diffraction phenomena that are not included in attenuation-based interpretation schemes commonly used in Radio Imaging Method (RIM) interpretation. These phenomena can be expected to lead to artifacts in ray theory interpretations. Proper placement and frequency of transmitters and receivers are issues that need to be rigorously examined prior and subsequent to any cross-borehole EM experiment. Using plate type models, we illustrate how the sensitivity of receivers to perturbations in the vertical and lateral extent of the plates is a function of the receiver placement and frequency of excitation. Of particular interest here is how rapidly sensitivities can vary with spatial position and frequency. For this examination, we used the FDTD code to calculate data sensitivity maps for perturbations to three-dimensional (3D) models, and an analytic solution for the fields due to a perfectly conducting half-plane developed by Weidelt to estimate parameter uncertainties for noisy 2.5D model data generated using an array of transmitters. The thesis concludes with a discussion of weighting of transmitters to optimize resolution of a perturbation about an a priori model. For the test case, the optimal weights improve the resolution of the top edge of the half-plane by 20%. The theory could also be used to design a phased array device optimized for high resolution studies. The Weidelt code was also used in this study. Despite advances in the development of absorbing boundary conditions, the FDTD method is still an expensive solution to the 3D EM modeling problem. Such was the rationale for using a fast analytic solution to the 2.5D half-plane induction problem to orient survey design methods.

ELECTROMAGNETIC by David Michael Johnson Geophysics Department of Geology and Geophysics The University of Utah June 1997 FINITE DIFFERENCE TIME DOMAIN MODELING OF CROSS-HOLE ELECTROMAGNETIC SURVEY DATA by A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science in Copyright © David Michael Johnson 1997 Rights Reserved © All THE UNIVERSITY OF UTAH GRADUATE SUPERVISORY COMMITTEE APPROVAL of a thesis submitted by This thesis has been read by each member of the following supervisory committee and by majority vote has been found to be satisfactory. >w H . iwi AIL. C - i,w_. /i TChhaaiirr-: AAllaann CC . TTrriirpmp 1 » I./997 6uJdL LA Erich U. Petersen PC. 09. % HCOLJ 2Ucnj*/ Michael Zhdanov THE UNIVERSITY OF UTAH GRADUATE SCHOOL David Michael Johnson '- . \ ~ y) 19~7 Chair: Alan C. Tripp n DC () 9. /.J. Michael Zhdanov SCHOOL FINAL READING APPROVAL I have read the thesis of David Michael Johnson in its final form and have (2) its illustrative materials including figures, tables and charts are in place; and (3) the final manuscript is satisfactory to the supervisory committee and is ready for submission to The Graduate School. Date ' AAllaann CC.. TTrriipppp \ Y Chair, Supervisory Committee Approved for the Major Department Bartle^ Chair Approved for the Graduate Council ' ' " "AnnW.Hart THE UNIVERSITY OF UTAH GRADUATE SCHOOL To the Graduate Council of the University of Utah: found that (1) its format, citations and bibliographic style are consistent and acceptable; chans Alan C. Tripp l""( Chair, Supervisory Committee Depanment John M. Banley AnnW. Han Dean of The Graduate School Cross-borehole electromagnetic (EM) surveying is an emerging technology in the mining industry that can offer resolution of geoelectric structures superior to that from other EM methods. However, existing interpretation techniques often employ physical models of electromagnetic wave propagation that are inappropriate for the instrument parameters and electrical properties commonly encountered in the metal­liferous mining environment. Lacking appropriate numerical models, we are unable to properly assess imaging methods or to design and interpret surveys. In this thesis, I address such issues relating to the delineation of massive sulfide mineralization using cross-borehole EM techniques. Specifically, I discuss the use of the finite-difference time-domain (FDTD) method in modeling the cross-borehole electromagnetic response of a perfect electric conductor (PEC), which is a good ap­proximation to a target of high conductivity. This model type has particular appli­cation to massive nickel sulphide ore bodies, such as those found in the Kambalda, Sudbury Basin, Norilsk, and Victoria Island regions. A key development of the thesis is the implementation of perfectly matched layer (PML) absorbing boundary conditions in a Yee FDTD algorithm. This revised al­gorithm permits the realistic modeling of the electromagnetic response of conductors in the frequency range transitional between the diffusion and wave domains. Such modeling reveals diffraction phenomena that are not included in attenuation-based in­terpretation schemes commonly used in Radio Imaging Method (RIM) interpretation. These phenomena can be expected to lead to artifacts in ray theory interpretations. Proper placement and frequency of transmitters and receivers are issues that need to be rigorously examined prior and subsequent to any cross-borehole EM experiment. Using plate type models, we illustrate how the sensitivity of receivers to perturbations in the vertical and lateral extent of the plates is a function of the receiver placement ABSTRACT borehole is technology in resolution to methods. However, difference Vee conductors are schemes commonly used in Radio Imaging Method (RIM) interpretation. These phenomena can be expected to lead to artifacts in ray theory interpretations. Proper placement and frequency of t ransmitters and receivers are issues need to be rigorously examined and subsequent to any cross-borehole EM experiment. Using plate type models, we illustrate how the sensitivity of receivers to perturbations in the vertical and lateral extent of the plates is a function of the receiver placement vary with spatial position and frequency. For this examination, we used the FDTD code to calculate data sensitivity maps for perturbations to three-dimensional (3D) models, and an analytic solution for the fields due to a perfectly conducting half-plane developed by Weidelt to estimate parameter uncertainties for noisy 2.5D model data generated using an array of transmitters. resolution of a perturbation about an a priori model. For the test case, the optimal weights improve the resolution of the top edge of the half-plane by 20%. The theory could also be used to design a phased array device optimized for high resolution studies. The Weidelt code was also used in this study. FDTD rationale for using a fast analytic solution to the 2.5D half-plane induction problem to orient survey design methods. and frequency of excitation. Of particular interest here is how rapidly sensitivities can FOTO 3~) 2.50 The thesis concludes with a discussion of weighting of transmitters to optimize Despite advances in the development of absorbing boundary conditions, the FOTO method is still an expensive solution to the 3D EM modeling problem. Such was the 2.50 v ABSTRACT iv LIST OF FIGURES vii ACKNOWLEDGEMENTS ix 1. INTRODUCTION 1 1.1 Cross-borehole EM techniques 1 1.3 resolution 2. FDTD MODELING 8 transmitters 2.5 Modeling perfect electric conductors in resistive hosts 20 3. RADIO IMAGING METHOD MODEL STUDIES 29 3.1 3.2 Sensitivity of RIM measurements to geophysical model parameters . . 32 3.3 3.4 45 APPENDICES CONDITIONS . . CONTENTS ABSTRACT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . ..... . . .. .... . . .... . ... . ... .. .. . . . ... . . . .... . ........ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . IX INTRODUCTION... ........ .. ........ ... ... . ..... . ... . . . .. . . .. . . .. .... . . . . . . . . . . . . . . . 1.2 Metalliferous mining applications of cross-borehole EM . 3 1. 3 Numerical electromagnetic modeling .. . . . . . . . . . 4 1.4 Cross-borehole EM survey design for optimal resolut ion . 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . 2.1 The need for wavefield calculations . . 8 2.2 Outline of FDTD modeling technique . . . . . . . . . . . 10 2.3 Discretization of dipole t ransmitters. . . . . . . . . . . . 11 2.4 Optimization of the PML absorbing boundary conditions 16 . . . . . . . . . . . . . . . . . . . . The Radio Imaging Method (RIM) .... . .. . . . . . ... . . . Confidence limits of parameter resolution . . .. . .. ... . . .. . Optimization of transmitter arrays .... .. ... . . . .. ... . 29 37 4. CONCLUSIONS . . . ... . .. . . .. . .. . . . . .... . . . . .... . . . .. . . ........ . . . . .. ... 48 53 A. PML ABSORBING BOUNDARY CONDITIONS. . . . . . . . . . . . . . . . .. 53 B. DISCRETE PML ABSORBING BOUNDARY CONDITIONS . . .. 60 C. PERFECTLY CONDUCTING HALF-PLANE MODEL .. ... ....... 72 REFERENCES . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 75 configuration conductivity - 0.0001 S/m, relative electric permittivity er = 10, relative permeability fir 1, and propagation distances of 50 m and 100 m as functions of frequency. 9 2.2 Error introduced by the diffusive approximation for a medium with conductivity a - 0.001 S/m, relative electric permittivity er = 10, relative permeability /j,r = 1, and propagation distances of 50 m and of frequency. Yee Yee, 1966) 2.4 Schematic diagram showing contour of integration C enclosing surface S. The contour integral of the impressed magnetic fields on C is equal to the flux density of equivalent impressed electric current flowing through S 2.5 Comparison of FDTD and analytic solutions for a pulsed magnetic dipole in a 0.001 S/m whole-space at varying distances from the trans­mitter. The cell size is 5 m, and the receiver is placed at (a) 5 cells, transmitter results, data, (b) Optimal profiles for alpha and beta 20 2.7 Cross section of the PEC half-plane model. The conductor extends infinitely to the right and perpendicular to the page 24 2.8 Comparison of the analytic result for a perfectly conducting half-plane due to Weidelt (pers. comm.) with the numerical solution obtained using the FDTD algorithm of Wang and Hohmann (1993) with PEC boundary conditions enforced explicitly at each time step. The calcu­lated impulse responses are compared at delay times of 10-6, 10~5 and 10~4 seconds 25 2.9 Plate model used to represent a perfectly conducting half-plane. . . . 26 2.10 Comparison of the vertical magnetic field response of the model de­picted in Figure 2.9 with the analytic solution for the analogous half-plane model 27 3.1 Perfectly conducting plate models 33 LIST OF FIGURES 1.1 General cross-borehole electromagnetic survey configuration. 2 2.1 Error introduced by the diffusive approximation for a medium with a = Slm, Er = J.tr = . . . . . . . . . . . . . . . . . . . .. = Slm, Er = J.tr = 100 m as functions offrequency. . . . . . . . . . . 10 2.3 The Vee staggered grid (adapted from Vee, 1966). . . . . . . . . . . . 11 S. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 13 Sim (b) 10 cells, (c) 15 cells, and (d) 20 cells from the transmitter. . . . . 15 2.6 Optimized PML profile results. (a) Comparison of observed and fitted beta. ... . . . . . . . . .. page. . . ..... . . 10- 6 , 10- 10-4 seconds. . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . model. . . . . . . . . . . . . . models. B at clipped at -500 dB. The diffractive nature of the field behavior is evident. 35 been clipped at -30 dB. Contours at -5, -10, -15 and -20 dB are shown. The region of maximum sensitivity forms a wedge terminating at the plate field to the unperturbed model response at 500 kHz, where the minimum has been clipped at -30 dB. Contours at -5, -10, -15, and -20 dB are shown. The region of maximum sensitivity is more compact than that shown in Figure 3.3 37 3.5 Vertical magnetic field amplitude perturbations in decibels relative to been clipped at -30 dB. Contours at -5, -10, -15, and -20 dB are shown. The zone of maximum sensitivity has a greater vertical width than in example 3.6 Faulted and intact perfectly conducting plate models. The dashed line shows the position of the upper 30 m of the orebody in each case. The orebody been clipped at -30 dB. Contours at -1, -5, -10, -15, and -20 dB are taken above the horizontal fault plane and behind the upper faulted plate segment contain little information regarding the model perturbation 40 3.9 Perfectly conducting vertical half-plane lying between two vertical bore­holes. An array of VMD transmitters is placed in bore-hole 1, and an array of magnetic field receivers is placed in bore-hole 2. The half-plane extends to infinity downward and perpendicular to the plane containing the bore-holes 42 3.10 Schematic diagram showing two half-plane models in which the vertical position of the edge is perturbed by 20 m 43 Vlll 3.2 Horizontal magnetic field amplitudes for model Bat 2 MHz expressed in decibels relative to the maximum value. The minimum has been 3.3 Horizontal magnetic field amplitude perturbations in decibels relative to the unperturbed model response at 2 MHz, where the minimum has top edge of the perturbed plate. . . . . . . . . . . . . . . . . . . . .. 36 3.4 Horizontal magnetic fie ld amplitude perturbations in decibels relative 3.3. .. . . . ........ . .. . . . . .. . .. . . the unperturbed model response at 2 MHz, where the minimum has the previous example. . . . . . . . . . . . . . . . . . . . . . . . . . .. 38 dipole transmitter lies opposite the top edge of the orebody. ... . . 39 3.7 Vertical magnetic field amplitude perturbations in decibels relative to the unperturbed model response at 2 MHz, where the minimum has shown. The region of maximum field perturbation is broad. Measure­ments perturbation. . .. . ........................... 3.8 Vertical magnetic field amplitude perturbations in decibels relative to the unperturbed model response at 500 kHz, where the minimum has been clipped at -30 dB. Contours at -1, -5, -10, -15, and -20 dB are shown. The range of field sensitivity is decreased relative to that shown in Figure 3.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 41 half­plane holes. ..... .. .. . .......... . . . . m. . . . . . . . . . . . . . .. 15, and 20%. The parameter is estimated using data recorded at each . . . the 99% confidence region, given uniform data standard errors of 10, 15, and 20%. (a) Uncertainties for the vertical position, (b) Uncer­tainties for the lateral position. The uncertainty ranges for the lateral position are approximately twice as large as the uncertainty ranges for position array, arrays, 3.10 B.l Field components used to calculate at points lying on the interface mesh B.2 Field components used to calculate Ex at points lying on the interface mesh mesh C.l Projection of the half-plane model into the y-z plane illustrating the systems 3.11 Estimated uncertainties for the height of the half-plane edge within the 99% confidence region, given uniform data standard errors of 10, receiver for the set of 27 transmitters, and the uncertainty is plotted for each receiver number. Data recorded above the plate edge yield greater parameter uncertainty ranges. The parameter uncertainties associated with these data sets also increase more rapidly with noise level. 44 3.12 Estimated uncertainties for the position of the half-plane edge within position. the vertical position. . . . . . . . . . . . . . . . . . . . . . . . . . .. 44 3.13 Optimization of sensitivities for a weighted transmitter array. (a) Sen­sitivities of data for the uniformly and optimally weighted transmitter arrays. (b) The set of optimal weights applied to each transmitter. The model geometry is shown in Figure 3.10. . . . . . . . . . . . . .. 47 1 Hy between the PML and the interior region of the FD mesh. ..... . 68 between the PML and the interior region of the FD mesh. . . . ... 69 B.3 Field components used to calculate Ez at points lying on the interface between the PML and the interior region of the FD mesh. .. . . .. 69 1 Cartesian and cylindrical coordinate systems. ... . ......... 72 IX Financial support for this work was provided by the Consortium on Electromag­netic Modeling and Inversion at the Department of Geology and Geophysics, Univer­sity of Utah, and by Western Mining Corporation Ltd. I also thank my advisor, Dr. Alan C. Tripp, for his valuable collaboration on technical matters, and for his editorial direction. The time Dr. Tripp spent super­vising this work was funded by the DOE. I am grateful to Dr. Michael S. Zhdanov, Dr. Gerard T. Schuster and Dr. Erich U. Petersen for serving on my supervisory committee. discussions on FDTD modeling. wholespace. Oleg estimation code, and I thank him for several useful discussions during the course of my research. Su-percomputing Institute at the University of Utah. Dr. Kurt Strack and Mr. Berthold Kriegshauser are also to be thanked for making computing facilities at Western Atlas ACKNOWLEDGEMENTS I am grateful to Dr. Cynthia Furse and Dr. Om Gandhi at the Department of Electrical Engineering, University of Utah for providing me with their implementation of the Yee FDTD modeling algorithm, from which I developed the code used in this project. Dr. Cynthia Furse and Dr. Gianluca Lazzi are to be thanked for many useful Dr. Peter Weidelt at the Institut fur Geophysik und Meteorologie in Braun­schweig, Germany generously provided me with an unpublished solution for the tran­sient electromagnetic response of a perfectly conducting halfplane in a conducting I wish to thank Dr. Elena Cherkaeva for consultation regarding optimal transmit­ter arrays. Mr. OIeg Portniaguine provided valuable assistance with his parameter research . I gratefully acknowledge the use of the computer facilities provided by the Utah Su­percomputing Institute at the University of Utah. Dr. Kurt Strack and Mr. Berthold Kriegshiiuser are also to be thanked for making computing facilities at Western Atlas Logging Services available. 1. INTRODUCTION placed above, on, or below the ground surface. By placing arrays of transmitters and case of underground mine exploration, a cross-borehole EM survey can yield better in separate boreholes, in order to image the intervening geoelectric structures (Fig­ure 1.1). Existing methods of cross-borehole EM surveying fall into two main cate­gories. One group of systems, termed cross-borehole radar, employs electric dipoles to transmit an electromagnetic pulse with a center frequency in the 10-200 MHz range. Receivers are deployed to measure either the reflected (Fullagar and Livelybrooks, 1994) or hole-to-hole (Olhoeft, 1988) time-domain signal. On the other hand, the ra­dio imaging method (RIM) group of cross-borehole systems employ either a magnetic or electric dipole to transmit a continuous-wave signal. The associated frequencies range from several kilohertz (Wilt et al., 1995) to tens of megahertz (Nickel and Cerny, INTRODUCTION 1.1 Cross-borehole EM techniques Electromagnetic transmitters and receivers used in geophysical surveys can be receivers closer to the subsurface region of interest, we can attain higher resolution of the geoelectric structure, subject to an appropriate survey design. Therefore, in the resolution than, say, a surface-to-borehole survey. In many areas, the basement rocks that contain the target of interest are overlain by conductive, inhomogeneous lay­ers of porous sediments and/or weathered rock, which scatter and absorb EM waves transmitted from an above-ground source. Placing the transmitter a sufficient dis­tance beneath the ground surface reduces contamination of the signal by this geologic noise. Underground deployment of receivers also ameliorates these conditions, and reduces contamination of the data by atmospheric and cultural sources of electromag­netic noise. Subsurface transmitters and receivers are restricted to small moments and coarse lateral sampling and so constitute major disadvantages in this method. In this thesis, I consider the case where transmitters and receivers are placed wave aI., Figure 1.1: General cross-borehole electromagnetic survey configuration. 1989), depending upon the application, so that there is an overlap in the frequency ranges of the radar and RIM type systems. Cross-borehole electromagnetic techniques have been applied to a variety of geo­logical and engineering problems; including tunnel detection (Lytle et al., 1979; Ol­hoeft, 1988), reservoir characterization and monitoring of enhanced oil recovery (Wilt et al., 1995), and coal seam delineation (Stolarczyk et al., 1988). However, application of the method to the metalliferous mining environment is still in its infancy. 2 -/-,-R­X I Tx tunnel aI., 01- hoeft, 1995) , aI. , of the to the metalliferous is still its detect any mineralization which has hitherto eluded discovery by drilling. In this case, the objective is to estimate the position of the target with sufficient accuracy to enable further investigation by drilling. On the other hand, we may wish to resolve some specific geometric feature of a known body of mineralization. For instance, we may wish to define the down-plunge limit of a plate-like "shoot" of massive sul­fides, or we may wish to define the limits of a zone where mineralization has been stoped out by igneous intrusions. In these cases, we usually possess detailed a priori knowledge of the target, and are seeking specific information regarding its geometry. This information could be obtained by drilling. We wish, however, to design surveys and interpretation techniques that provide information of comparable accuracy at a significantly lower cost than drilling. and by Nickel and Cerny (1989) to define the edges of sulfide-rich zones. Data are model for each transmitter-receiver pair. This estimated attenuation is assumed to represent the integrated attenuation constant along a straight ray connecting the transmitter and receiver. Conductive sulfide mineralization lying on the ray path is expected to give rise to an increased attenuation. This model of electromagnetic wave propagation also forms the basis of tomographic reconstruction techniques commonly applied to RIM data, and will be discussed in detail in a later section. Stolarczyk (1992), Thomson et al. (1992), Wedepohl (1993), and Pitts and Kramers (1996) present case histories illustrating the use of tomographic reconstruc­tion to image geoelectric structure in a cross-borehole plane. In each case, the ampli­tude or phase of the axial electric or magnetic field component is measured at a set 1.2 Metalliferous mining applications of cross-borehole EM 3 This thesis examines the application of cross-borehole EM techniques to delin­eation of massive sulfide mineralization in an underground mining environment. In this case, the objectives of geophysical surveying are usually very specific. On the one hand, we may wish to conduct a reconnaissance of a given region, in order to Simple "radio shadow" techniques have been employed by Rao and Rao (1983) reduced by calculating the signal attenuation relative to an estimated background 1992) , 4 of points in one borehole for a set of transmitter deployments in another borehole. These data are used to recover a distribution of attenuation or phase constants. All of these authors seek to validate their recovered models by comparison with geologic data obtained from boreholes intersecting the image plane defined by the transmit­ter and receiver boreholes. In all cases, inconsistencies exist between the predicted and actual mineralization geometry, and the comparisons do not allow us to evaluate the accuracy of reconstructions in areas where no drilling information is available. These case histories therefore provide an unsatisfactory basis upon which to assess the usefulness of the survey and interpretation techniques. Numerical modeling can be used to more effectively evaluate the resolving power orebody 105 - 106 S/m. 1988). We will see that it is necessary to include both conduction and displacement currents when modeling high-frequency cross-borehole EM surveys. 1.3 Numerical electromagnetic modeling As mentioned above, a common assumption made in the interpretation of cross-borehole EM data is that interactions between electromagnetic waves and the geo­electric structure can be described by a simple ray-theory model (e.g., Stolarczyk, 1992). The amplitude or phase of the recorded field component is assumed to reflect the integrated attenuation or phase constant along a ray connecting the transmitter of points in one borehole for a set of transmitter deployments in another borehole. These data are used to recover a distribution of attenuation or phase constants. All of these authors seek to validate their recovered models by comparison with geologic data obtained from boreholes intersecting the image plane defined by the transmit­ter and receiver boreholes. In all cases, inconsistencies exist between the predicted and actual mineralization geometry, and the comparisons do not allow us to evaluate the accuracy of reconstructions in areas where no drilling information is available. These case histories therefore provide an unsatisfactory basis upon which to assess the usefulness of the survey and interpretation techniques. of cross-borehole EM surveys and the accuracy of imaging techniques. But numerical models are useful only insofar as they incorporate appropriate petrophysical informa­tion. The sulfide mineralogy in a number of important classes of ore body is dominated by highly conductive pyrrhotite, giving ore conductivities of about - 106 Sjm. Examples include the Kambalda (Gresham and Loftus-Hills, 1981), Sudbury Basin (Stanton, 1972), and Voisey Bay nickel deposits. All of these orebodies occur within highly resistive metamorphic rocks, so that the conductivity contrast between ore and host rock can be as high as 10 orders of magnitude. In the case of highly resistive rocks, dielectric polarization effects become significant at high frequencies. Common rock-forming minerals such as quartz and feldspar have relative dielectric constants in the range 4-8, whilst water has a relative dielectric constant of about 80 (Keller, 1988). We will see that it is necessary to include both conduction and displacement currents when modeling high-frequency cross-borehole EM surveys. cross­borehole g. , 5 and receiver. Nekut (1994) has described the conditions under which this model can be used as a basis for interpretation. Two key assumptions are that measurements are taken in the far-field (i.e., \kr\ > 1), and that conductivity contrasts are small. Both of these assumptions are violated in most metalliferous mining applications of cross-borehole EM, where wavelengths in the resistive host rocks are long relative to the scale of the survey, and conductivity contrasts are up to 10 orders of magnitude. Fullagar and Pears (1995) applied a ray-trace tomography code to synthetic data gen­erated using a conductive sphere model. They found that an image rich in artifacts, but approximately resolving the sphere location, could be generated for cases where the sphere lay between the survey boreholes. However, placing the sphere outside the survey plane resulted in spurious resistive features in the tomogram. These results demonstrate that data that record scattering from objects lying outside the survey plane may not be accurately interpreted using the ray-theoretic model. Yee and receiver. Nekut (1994) has described the conditions under which this model can be used as a basis for interpretation. Two key assumptions are that measurements are taken in the far-field (i. e., Ikrl ~ 1), and that conductivity contrasts are small. Both of these assumptions are violated in most metalliferous mining applications of cross-borehole EM, where wavelengths in the resistive host rocks are long relative to the scale of the survey, and conductivity contrasts are up to 10 orders of magnitude. Fullagar and Pears (1995) applied a ray- trace tomography code to synthetic data gen­erated using a conductive sphere model. They found that an image rich in artifacts, but approximately resolving the sphere location, could be generated for cases where the sphere lay between the survey boreholes. However, placing the sphere outside the survey plane resulted in spurious resistive features in the tomogram. These results demonstrate that data that record scattering from objects lying outside the survey plane may not be accurately interpreted using the ray-theoretic model. Other traditional numerical interpretation methods also may be inadequate for modeling the full wavefield EM response of a highly conductive body in a lossy dielec­tric medium. Early workers (Lajoie and West, 1976; Hohmann, 1983) recognized that the pulse subsectional basis functions commonly employed in method-of-moments so­lutions fail to accurately represent the scattering currents in an inhomogeneity as the conductivity contrast increases. For instance, Newman et al. (1986) claim that their 3D integral equation solution is accurate for conductivity contrasts up to 300 : 1. In an effort to overcome these numerical instabilities, San Filipo et al. (1985) introduced specialized basis functions to represent the curl-free and divergence-free components of the scattering currents. However, the improved accuracy of this solution is achieved at the expense of reduced geometric flexibility. Recently, Liu and Lamontagne (1995) presented a surface integral equation solution for discrete, homogeneous ellipsoidal targets, and claimed accurate results for high conductivity contrasts over a broad range of frequencies. These results are encouraging. However the application of this method to more general geometries is expected to involve labourious programming. In contrast, finite-difference techniques employing the Vee (1966) method of dis­cretization for electric and magnetic field components possibly offer a means for accu- 6 rately modeling high-contrast targets. Boundary conditions across physical interfaces are automatically satisfied, removing a major source of numerical error. Alumbaugh and Newman (1996) implemented a finite-difference frequency-domain solution that permits a high degree of geometric flexibility by making the simplifying assumption that the highly conductive objects are perfect electric conductors (PECs). This ap­proach may be useful in the case of a massive sulphide body composed mainly of a highly conductive mineral, such as pyrrhotite, where the contact between massive ore and host rock is sharp. The footwall contact of a Kambalda-style nickel sulphide orebody is an instance when the PEC model may be applicable. This thesis will concentrate on the FDTD method of calculating the full wavefield EM response of a perfect electric conductor in a conducting host. 1.4 Cross-borehole EM survey design for existing boreholes and underground mine workings in order to gain access to the region of interest. However, many surveys carried out in this manner will fail to resolve the features of interest. For instance, an electric dipole transmitter deployed in a borehole that is perpendicular to a conductive plate-like body will be null-coupled to the target, giving no scattered field response. Forward modeling which incorporates a priori information can be used to design a survey in which the data have maximum sensitivity to the features of interest. As an example, we will use the modeling code described in this thesis to map the sensitivity of frequency domain EM data to perturbations in the position of the edge of a conductor. In planning a survey to estimate this parameter, we would take measurements in the regions of maximum sensitivity. These maps of the Frechet derivatives for a given model parameter are useful, insofar as they predict the resolving power of noise-free data. We also used the parameter covariance matrix to estimate parameter uncertainties for various noisy measurement sets. Arrays of transmitters are employed in cross-borehole EM, in order to achieve optimal resolution Cross-borehole geophysical surveys are often opportunistic in their planning, ex­ploiting data to perturbations in the position of the edge of a conductor. In planning a survey to estimate this parameter, we would take measurements in the regions of maximum sensitivity. These maps of the Frechet derivatives for a given model parameter are useful, insofar as they predict the resolving power of noise-free data. We also used the parameter covariance matrix to estimate parameter uncertainties for various noisy measurement sets. shown, in the context of DC resistivity surveying, that it is possible to design a weighted array of transmitters which optimizes the sensitivity of a given set of mea­surements to perturbations from a given a priori geoelectric model. We apply the same method to design an array of borehole transmitters that produces optimal defi­nition of the position of a conductor edge. This approach has particular application to cross-borehole EM surveying in the mine environment, where we are often required to define particular parameters of the orebody geometry, and useful a priori information is available to aid survey planning. 7 varying illuminations of the region of interest. Cherkaeva and Tripp (1996) have ore body for the frequency range and physical properties commonly encountered in high resolu­tion hardrock discretization of the exact Maxwell's equations, with specialized time-stepping equa­tions for the impressed magnetic and electric currents. The Perfectly Match Layer (PML) absorbing boundary condition is used to deal with the serious problem posed by truncation of a finite-difference mesh in the low-loss wave propagation regime. The parameters of the PML must be optimized for each model. These theoretical a perfect electric conducting (PEC) target. 2.1 The need for wavefield calculations modeler. This is so within an FDTD context because a purely diffusive mode demanding Yee algorithm. Unfortunately the present exploration application requires a full wavefield treatment, as the following discussion shows. Consider two electromagnetic plane waves propagating with exact ki = -\j/J.eu>2 - ifiaco, and diffusive ft2 = - y - ifiacj, wavenumbers. these plane waves propagate a distance then the relative error in 2. FDTD MODELING Neglecting displacement currents (i.e., the diffusive assumption) is inappropriate imaging in hard rock mining environments. Therefore, FDTD modeling requires and practical considerations are addressed below, with application to scattering from The diffusive assumption, where justified, greatly eases the burden of the electro­magnetic of propagation permits the use of the unconditionally stable modified Du Fort-Frankel FDTD algorithm (Wang and Hohmann, 1993), rather than the more computationally Vee k, = -J fJ,EW2 - ifJ,CIW, diffusive If d Error = 3-ik\d -ik2<l\ -ikid\ (2.1) = S/m, electric permittivity er = 10, and relative magnetic permeability fir = 1. If the propagation distances are 50 m and 100 m, then the relative errors as functions of frequency are described by the curves in Figure 2.1. The errors at 100 kHz are. about 10% for d 50m and 20% for d 100m. When a 0.001 S/m, then the error at 100 kHz is approximately 4% for d = 50m and 8% for d = 100m (Figure 2.2). Clearly, the diffusive approximation is not appropriate for high resolution modeling above 100 kHz. diffusion and propagation regimes, so that both conduction and displacement currents must 10 Frequency (Hz) Figure 2.1: Error introduced by the diffusive approximation for a medium with con­ductivity a = 0.0001 S/m, relative electric permittivity er = 10, relative permeability Hr = 1, and propagation distances of 50 m and 100 m as functions of frequency. amplitude introduced by the diffusive assumption is given by = le-iktd _ C ik2d l le-iktdl 9 For example, consider a host medium with conductivity a 0.0001 81m, relative ir = J.lr = are = = = 81m, kHz is approximately 4% for d = 50m and 8% for d = 100m (Figure 2.2). Clearly, the diffusive approximation is not appropriate for high resolution modeling above 100 The models discussed in this thesis lie within the transition between the diffusion be included in order to accurately simulate the EM response. 10' r-:""'TCf.T..r. .C . :~==0=·-'.-.' T ··="':'·[·-·:·-T.. " .- 'J'!'.:..:.= ..:.. : .= .. :..: .: .:f ,,·T.. ·j 102 ..... '.:::: ......... ":::'.::::: c .. :::, .,.·.,::'·"·""": ..... ::: ~: d;S()m ::: : ::::: .. .. ", ... ', .. , .. : .. :.: :.:-: .... • •• " / , 0 • . . ., .( . . .. ;, ,: .J': . y( .. "' .... ," ." : .... : .. . ; . ; ,',,' ~;.~ ~ : d~16bitL . .......... ."..'.,. .......... . . . _. . . .... ; . ,.; ; con­ductivitya = 81m, ir = relatIve permeablhty J.lr fun ctions 10 10° Frequency (Hz) Figure 2.2: Error introduced by the diffusive approximation for a medium with con­ductivity o = 0.001 S/m, relative electric permittivity er = 10, relative permeability Hr 1, and propagation distances of 50 m and 100 m as functions of frequency. Inclusion of both electric displacement and conduction currents demands that the -//- OH dt' (2.2) „ dE V x H = aE + e-, (2.3) where // is the time-invariant magnetic permeability, e is the time-invariant electric permittivity, and a is the time-invariant electric conductivity. The modeling technique used in this thesis is a straightforward implementation the Yee (1966) method, as described by Taflove and Umashankar (1989). The model is discretized into a uniform cubic grid with a nodal spacing of Ax. The electric and 10° . . . ... : .. '-; :, ':::::( : ... ~ ... . .. J ....... , .;.~.,,; . .... , - . .. '.' ;;....;.;.;..: d~rrl :: ;c d~16oiri .. . ".,.~ ,., ,,' .. ; .. ', .. " .. , 10 a Sl m, fr J1r = frequen cy. 2.2 Outline of FDTD modeling technique exact form of Maxwell's equations be used. These are and aH 'V x E = -H­... at ' (2.2) (2.3) where J1 f electric permittivity, and a is time-electric conductivity. The modeling technique used in this thesis is a straightforward implementation of the Yee (1966) method, as described by Taflove and Umashankar (1989) . The model is discretized into a uniform cubic grid with a nodal spacing of 6.x. The electric and magnetic fields are sampled in a staggered fashion at half-cell intervals on this grid (Figure 2.3), and are updated on staggered time-steps of duration At. The space and time derivatives in Equations 2.2 and 2.3 are approximated by standard second-order central differences. matched layer (PML) material absorbing boundary condition (ABC) formu­lated Appendix A, and the implementation is discussed in Appendix B. Frequency-domain model data were derived from the time-domain data using a discrete Fourier transform (DFT). The DFT has been shown by Furse and Gandhi (1995) to be the most efficient method for computing frequency domain data from the output of an FDTD code for multiple frequencies. transmitters Infinitesimal magnetic and electric dipoles can be used to model the compact transmitters employed in crosshole EM experiments. order to introduce trans­mitters impressed magnetic and electric currents in a local region. this region, Maxwell's Ax- Figure 2.3: The Yee staggered grid (adapted from Yee, 1966). 11 !:J.t. differences. Field values on the boundaries of the model region are calculated using the per­fectly by Chen et al. (1996). The principles underlying the PML are discussed in 2.3 Discretization of dipole transmitters In to the FDTD mesh, we must modify the time-stepping equations to include In E,f----+--~ _---Ill---- Vee Vee, (9E V e- + aE + J\ (2.4) (9H V x E = - / i - - M\ (2.5) J* IvT These expressions can be discretized to yield time-stepping equations which model the transmitter. x, y be discretized by controlling the corresponding field component at that node. This write the differential form of Faraday's law as V X E = - ^ - ^ ' (2' 6) H* fj,^- Ml space and time: Hr**\ak = flr-°-8iy*+^wi«*-^i«-y*- 12 equations become aE . \7 x H = fa-t + erE + J ' ' and aH . \7 x E = -,H., a-t M' ' (2.5) where Ji is an impressed electric current and Mi is an impressed magnetic current. Electromagnetic field components within the FDTD mesh are defined at discrete points in space, and are assumed to vary linearly between adjacent nodes. A dipole oriented in the x , or z-direction and placed at a node in the FDTD mesh can effectively replaces an infinitesimal current element with a current density distributed over the finite volume of one cell. The field component must be scaled to simulate the effect of an infinitesimal dipole. Buechler et al. (1995) derived a time-stepping formula for the electric field in order to approximate a Hertzian electric dipole transmitter. Using a similar argument, we can derive an approximation for a magnetic dipole. We (2.6) where Hi is the impressed magnetic field intensity. The term jJ, a~i is an impressed magnetic displacement current which corresponds to Mi in Equation 2.5. Then, for example, the discrete time-stepping equation for Hz at the source location can be written by approximating the derivatives in Equation 2.6 by central differences in Hn+O . 5 I" z t)k Hn- O.5 1 + 6.t {En I··k _En I· j ·k- z ijk A x 1) X t - J w.XjJ, Ey \ijk +Ey \ij-ikj - At 13 where the impressed magnetic displacement current is evaluated at time-step n. by a FDTD cell by adding an impressed magnetic field -At^- at every time H* Consider an impressed magnetic field intensity FT defined on a closed contour enclosing an area The equivalent electric current density J* is related to the IF f IT Jc dl f r • ds. Js This expression is completely general, so that the magnetic field may be nonzero at only one point on the contour C, as in the case of an infinitesimal linear magnetic current element. Now we evaluate these integrals for the discrete case of a contour surrounding a Yee cell (Figure 2.4). Assuming that the field at a node represents the average field over the corresponding side of the square contour, the integrals become Ax . 0 Hz=H 'i Hz=0 Ax Enclosed surface S Hy=0 Figure 2.4: Schematic diagram showing contour of integration C enclosing surface S. The contour integral of the impressed magnetic fields on C is equal to the flux density of equivalent impressed electric current flowing through S. (2.7) n . We have introduced an impressed magnetic current within the finite volume de­fined -.6.ta.rr' step. We now determine the relationship between the quantity Hi and the moment m of the dipole we are attempting to simulate. Hi C S. Ji magnetic field intensity Hi by Ampere's law (2.8) Vee _ ___ ll, _ __ _ Hy=O .... -..... .. .. . . i / :........- ... . Contour C Hy=O HlAx = J1 Ax2 = J, (2.9) flowing dW 1 dl , ~df = Aim2' <2-10) rrn+0.5 I _ nrn-0.5 l , _~_ r r n i cm i nz Ujk - nz \ijk + ^ T i&x Hjk -&x U-ljk ~ Ey \ijk +E^ UJ-IA} - -^-QI- (2 - n) If we assume that the current I flows in a square contour of dimension Ax containing the electric field nodes that surround the transmitter location, then the magnitude of the corresponding magnetic dipole moment is simply m = I Ax2, so that | | = -^p^fjr- Therefore, we can rewrite the time-stepping equation with the impressed field given in terms of a time-varying magnetic dipole moment rj-n+0.5 I rjra-0.5 I i / p n I cm | tiz \ijk = tiz \ijk + ^ ~ i&x Hjk -&x U-ijk - Ey Ujk +Ey \ij-ik) ~ ~^-fo- (2"12) S/m, in time as m(t) = yje-'»-*>ax, (2.13) 14 where I is the total current flowing through the surface S. Differentiating with respect to time, we have 8W 8I at t:.x 8t z, (2 .10) where z is the unit vector describing the dipole orientation. Using this result, we can rewrite Equation 2.7 as H n+O . 5 1· · z t]k H n- O.5 I t:.t {En I En I - z ijk + 6.XJ..L x ijk - x i-ljk- En I En I } t:.t 81 y + y ij-lk t:.X 8t · (2.11) t:.x = 1 t:.x2, ~; = ",~, ~';' . n-O.5 I t:.t {En I En I Hz ijk +~ x ijk - x i-ljk- '-"XI' n n } t:.t 8m Ey lijk +Ey lij-lk - t:.x3 at . (2.12) This method of incorporating finite sources was tested by calculating the magnetic field for a pulsed horizontal magnetic dipole, contained in a conducting whole-space. The conductivity of the whole-space was 0.001 81m, the relative electric permittivity and relative magnetic permeability were equal to one, and the dipole moment varied = /!ie-O(t-tO)' x, (2.13) 0 = ^, 2 ' 3.4i Here uc 2irfmax with /max chosen to be 0.3 MHz. These model parameters were chosen to simulate a RIM-type transmitter operating in a moderately resistive host rock. The dipole transmitter was oriented in the x-direction, and placed at the center of a 50 x 50 x 50 cell mesh, with collinear magnetic field receivers placed at points along a line passing through the transmitter. The mesh resolution Ax was 5 m. diffusive impulse response (Ward and Hohmann, 1988) with the transmitter moment x10 E 3-0.5 S0.5 FDTD q(10-'2 0.4 0.6 0.8 O^Q-13 0.2 0.4 0.6 0.8 x10 10 from source Figure 2.5: Comparison of FDTD and analytic solutions for a pulsed magnetic dipole S/m transmitter. size m, cells from the transmitter. in which and 2 f} = We 2 ' to = 3.4~. 15 We = 27rlmax Imax ~x For comparison (Figure 2.5), an analytic solution was obtained by convolving the waveform. The model parameters were chosen such that the diffusive approximation is appropriate. The numerical results calculated using the FDTD algorithm were found to be in error at small distances from the transmitter, but converged to the correct x 1 !\ ' , - - FOlD result 5 , 5 cells from source / - Analytic solution 0< t 10-12 0.2 OA 0.6 0.8 1 1.2 1.4 !o:[ A : :,0 cells "o~ sou,:'0 1 0 Ck 10-13 0.2 0.4 0.6 0.8 1.2 1.4 ([ ~ : ' , x 10 1 15 cell, fm~ sou,ce : q 10-14 0.2 OA 0.6 0.8 1 1.2 1.4 fl ~ :20 cell, fm~ ,ou,:'0 1 00 0.2 OA 0.6 0.8 1.2 lA Time (sec) x 10-5 in a 0.001 81m whole-space at varying distances from the transmItter. The cell sIze is 5 m , and the receiver is placed at (a) 5 cells, (b) 10 cells, (c) 15 cells, and (d) 20 solution at greater distances. This result is not surprising, since we are approximating an infinitesimal dipole with discretized impressed fields. However, the numerical model is sufficiently accurate at distances greater than 10 cells from the transmitter. Note that the field computed 5 cells from the edge of the mesh appear to be unaffected by the presence of the mesh boundary. methods to geophysical modeling is the simulation of an unbounded region within the earth with a computational domain of finite extent. To address this issue, Berenger (1994) developed the Perfectly Matched Layer (PML) technique, which a modified set of Maxwell's equations and boundary conditions are used to create an unphys-ical zone which electromagnetic waves are attenuated rapidly without reflection. Berenger's formulation of the PML assumes that the model host medium has infinite resistivity. Chen et al. (1996) have developed a PML suitable for truncating a volume which materials of finite conductivity impinge upon the boundary. These PML ab­sorbing boundary conditions can be used to effectively eliminate spurious reflections from the outer surface of the computational mesh. However, the effectiveness of the PML technique relies upon an appropriate choice of stretching factor profiles within the absorbing boundary layer. of Cartesian coordinates {x, y, z) -v {xsx, ysy, zsz}, 2.14) where the "coordinate stretching factors" sx, sy and sz are complex numbers. The authors define a stretching factor as s(p) = hi , 2.15J ujre cue 16 solution at greater distances. This result is not surprising, since we are approximating an infinitesimal dipole with discretized impressed fields. However, the numerical model is sufficiently accurate at distances greater than 10 cells from the transmitter. Note that the field computed 5 cells from the edge of the mesh appear to be unaffected 2.4 Optimization of the PML absorbing boundary conditions A fundamental problem in the application of finite-difference time-domain (FDTD) in unphys­ical in in The PML technique developed by Chen et al. (1996) is based upon a transforma­tion (2 .14) s., Sy Sz s s(p) = PCp) + ia(p) , Wc€ WE (2.15) 17 uc elements of the Cartesian basis. The functions a(p) and 0(p) control the rate of attenuation within the PML in the direction p. Choosing a(p) = 0 and /?(p) = uce reduces the PML equations to Maxwell's equations, so that the rate of attenuation within the medium is controlled by its physical properties. In resistive media, this rate of attenuation may be too small to guarantee that all EM energy is dissipated before reaching the outer edges of the computational domain. Greater rates of attenuation can be attained by choosing larger values of a and /?. However, an attenuation rate that is too large gives rise to spurious numerical reflections. A compromise must be reached between insufficient attenuation, which causes EM waves to reach the outer boundary of the mesh and reflect back into the region of interest, and excessive attenuation, which causes spurious reflections to occur at interfaces within the PML. The performance of the PML absorber can be optimized by the application of parameter estimation. Let the exact solution to a given EM propagation problem be given by d°, which we will consider as observed data to be fitted by a set of predicted data dp produced by the FDTD numerical model. Now define A as the forward modeling operator that maps a set of PML coefficients m onto a set of predicted data d": dp = A(m). d° dp mapr. this case, ma T is simply our initial model. Thus we seek to minimize the regularized misfit functional 0(m?) = ||d' - d°||2 + «||m - mopr||2 = ||A(m') - d°||2 + K||m - mapr||2. (2.17) where We is the source centre frequency, and the vector p corresponds in turn to the ;3(= ;3(= WeE Maxwell 's ;3. dO, dP dP : dP m) . (2.16) Having fixed the physical parameters of the model, the PML coordinate stretching factors define the set of free model parameters. The error due to numerical reflections from the PML is equal to the misfit between the observed data dO and the predicted data dP which correspond to the model parameters m. We wish to minimize this misfit, whilst maintaining a model m which is close to some a priori model m apr · In this case, mapr is simply our initial model. Thus we seek to minimize the regularized functional 18 Here K is a regularization coefficient, which controls the degree to which minimization of the misfit functional is traded for minimization of variation in the model parameters from the a priori model. It is more important to obtain a low value of the misfit functional for a given background medium than it is to maintain a small degree of variation from the a priori coordinate stretching coefficients. Therefore, we can permit K to become small as we approach the desired minimum misfit. This problem can be solved using parameter estimation techniques. We applied a parameter estimation technique implemented by Portniaguine and Zhdanov (1995), which uses a regularized Newton method, with Frechet derivatives calculated by fi­nite differencing. In practice, this particular inverse problem is highly nonlinear, and therefore difficult to solve using this technique. Attempts to solve the inverse problem for independently varying coordinate stretching factors failed because the forward modeling algorithm is unstable with respect to perturbations in these pa­rameters, resulting in the Frechet derivatives being undefined. However, we can deal with this problem by parameterizing the real and imaginary parts of the PML co­ordinate stretching factor profiles in terms of simple functions, such as sines, whose shapes incorporate desirable properties, such as smoothness. In particular, the PML equations should match Maxwell's equations at the interior boundary of the layer, whilst the stretching factors may increase towards the outer boundary. In practice, it may be necessary to introduce a DC offset in the stretching factor profiles, to account for discretization of the PML. Thus we write ( \ ana a{n) = - ( 717T 1 + sin I- 7T ,NPML + b0 (2.18) and ag r . ( nir 7T 0(n)=u,ce + f [1 + ^ H N P M L _ 2 + b0, (2.19) where NPML is the PML thickness expressed as a number of cells, and n is the number of cells between the given node and the interior boundary of the PML. The variables aa, ba, ap and b3 are the new model parameters. If, t raded misfi t If, to become small as we approach the desired minimum misfit . t echnique t his part icular part icular, for discretization of the PML. Thus we write a", [ ( mr IT)] a(2 l +sin NPML -2" + b", 2. 18) a"" b"" ap bp 19 A ID EM wave propagation problem provides a forward modeling operator A which can be simply and rapidly computed. The PML coordinate stretching factors estimated using such a model may be applied to a 3D modeling problem with similar physical and numerical parameters. We constructed a ID FDTD model, in which a uniform plane-wave is discretized by staggered sampling of the field components Hz and Ey along a line in the x - direction. The fields are assumed to be invariant in the directions perpendicular to this line. A plane-wave pulse was introduced by forcing the vertical magnetic field in a mesh, representing plane perpendicular to the line. field by: t) = -26(t-t0)^e-^\ (2.20) V lu where 9 w2/2, and the center frequency is 1 MHz. The relative electric permittivity er 10, the magnetic permeability \i - fj,0 and the conductivity a = 0.0001 S/m. The effect of a boundless host medium was simulated by extending the finite-difference mesh a sufficient distance from the transmitter-receiver pair, such that no numerical reflections from the mesh boundaries were observed at the receiver during the time interval of interest. This model provided the ideal, reflection-free set of observed data d°. The transmitter-receiver pair were then moved to a point near the interface between the PML and the interior region of the mesh. The model data generated in this configuration constituted the predicted data 6P which we seek to fit to the observed data in a least-squares sense by adjusting the PML coordinate stretching factor profile. The parameter estimation program PAREST (Portniaguine and Zhdanov, 1995) was used to minimize the misfit functional given in Equation 2.17 by adjusting the parameters aa, 6Q, aB and bfi in Equations 2.18 and 2.19. The magnitude of the reflected pulse was reduced to negligible levels within a few iterations (Figure 2.6). The 3D model fields shown in Figure 2.5 were generated using an optimized set of PML coordinate stretching factors calculated by the ID modeling approach described 19 A 1D EM wave propagation problem provides a forward modeling operator A which can be simply and rapidly computed. The PML coordinate stretching factors estimated using such a model may be applied to a 3D modeling problem with similar physical and numerical parameters. We constructed a 1D FDTD model, in which a uniform plane-wave is discretized by staggered sampling of the field components Hz and Ey along a line in the x - direction. The fields are assumed to be invariant in the directions perpendicular to this line. A plane-wave pulse was introduced by forcing the vertical magnetic field at one point in a FDTD mesh, representing a plane perpendicular to the line. The magnetic field in this plane is given by: Hz(t) = -2(}(t - to)l!ie-O(t-tol', (2.20) () = w~/2 , €r = lO, J.l = J.lo and the conductivity 17 = 0.0001 S/ m. finite­difference receiver dO. receiver data generated configuration dP seek to fit to the observed data in a least-squares sense by adjusting the PML coordinate stretching factor profile. The parameter estimation program PAREST (Portniaguine and Zhdanov, 1995) was used to minimize the misfit functional given in Equation 2.17 by adjusting the parameters a", b", a{3 and b/3 in Equations 2.18 and 2.19. The magnitude of the reflected pulse was reduced to negligible levels within a few iterations (Figure 2.6). The 3D model fields shown in Figure 2.5 were generated using an optimized set of PML coordinate stretching factors calculated by the 1D modeling approach described s 1 .£" 0.5- c S or .93 i-0.5 c CO -1 | 0.002 Wealsolution Solution with 1500 Time step Cell index results, outer mesh boundary is good. resistive hosts The geophysical models discussed in this thesis consist of highly conductive mas­sive sulphide bodies situated within a weakly conducting host rock. In this case, the electromagnetic fields attenuate very rapidly within the conductor. For instance, if we consider electromagnetic fields varying with a frequency of 1 MHz incident upon sulphide mineralization (e.g., massive pyrrhotite) of conductivity 104 S/m, then the skin depth within the conductor will be 5 mm. Rather than attempting to explicitly model this rapid attenuation of the fields, we can approximate the object as a perfect electric conductor (PEC), where we assume that the fields vanish within the object. Several existing algorithms can be considered as potential solutions for this model- 20 E ~ f >- 1- ~~a\.soluti\W l ~ 0.5 ... .. a U Ion WI PML ABC :"<s= 0 '0 V Q; ;;: i"<=L 0.5 :'':";": -1 0 500 1000 2000 2500 3000 0.01 10 " 0.008 ~ :g 0.006 . <> ..... . .. '" ' " ~ 0.004 B i!i"i 0.002 0 1 2 3 4 5 6 7 8 9 10 11 Cell index Figure 2.6: Optimized PML profile results. (a) Comparison of observed and fitted data, (b) Optimal profiles for alpha and beta. above. Agreement between the analytic and numerical results within five cells of the 2.5 Modeling perfect electric conductors in 104 Slm, fields , of a perfectly conducting plate in free-space. This approximation, although useful, has a limited accuracy for conducting host media. Moreover, the thin plate model is lacking in geometric flexibility. Therefore, generalization of the algorithm to con­ducting host materials does not seem worthwhile for this application. Moreover, integral equation codes, which permit simulation of a conducting host, can become numerically unstable for high target-host conductivities (Hohmann, 1983; San Fil-ipo et al., 1985). In addition, our preliminary calculations suggested that neglecting displacement currents was not justifiable for the electrical properties and transmit­ter frequencies of interest. These considerations narrowed the field of candidates for codes with which to perform our modeling. target geometries afforded by the FDTD approach. outward propagating electromagnetic waves highly attenuate before reaching these perfectly reflecting surfaces. Expansion of the mesh beyond the region of interest is performed using gradually increasing cell dimensions, which is feasible in the diffu­sive regime because field variations become smoother with increasing distances from transmitters and scatterers. This code is appropriate for calculating the responses for lower frequencies used in this study. 21 ing problem. Annan (1974) has presented a numerical approximation for the response Fil­ipo aI. , I decided to further restrict the field of candidate algorithms to the class of FDTD methods which use the Yee staggered grid scheme to discretize the fields. One of my reasons was the ease with which FDTD algorithms can be programmed. Another is the apparent accuracy of the technique for full wavefield calculations (Wang and Tripp, 1996; Taflove and Umashankar, 1989). A third reason is the versatility in Three modeling codes were initially identified for experimentation and possible customization to the geophysical case at hand. The first code (TEM3DL), described by Wang and Hohmann (1993), was developed for the diffusion regime, and uses a modified version of the Du Fort-Frankel method to step Maxwell's equations in time on a Yee staggered grid. A homogeneous Dirichlet boundary condition is imposed on the fields at the subsurface boundaries of the computational mesh. The mesh boundaries are placed a sufficient distance from the region of interest to ensure that diffu­si ve lower frequencies used in this study. The second code, written by Dr. Cynthia Furse and others in the Department of Electrical Engineering, University of Utah, is an implementation of the Yee (1966) method using second-order accurate central difference approximations to the tempo­ral and spatial derivatives in Maxwell's equations. Derivation of the time-stepping equations is discussed in several texts (e.g., Taflove and Umashankar, 1989). Both electric conduction and displacement currents are included in the model. This code was originally developed for models contained within a region of free-space, and used Mur (1981) second-order absorbing boundary conditions (ABC) on the mesh bound­aries. This ABC is inappropriate for simulating broad-band EM responses of targets in dispersive media, since it is based upon a wave-equation in which phase velocity is assumed constant. An accurate ABC for lossy host materials was required. Attention has recently been focussed on Perfectly Matched Layer (PML) absorbing boundary conditions, which reportedly lead to superior accuracy (Berenger, 1994; Chew and Weedon, 1994; Chen et al., 1996). I decided to add a PML ABC for lossy media to the Yee algorithm implemented by Furse et al. (pers. comm.). Differing formulations of the PML for lossy media have been published by Fang and Wu (1995) and by Chen et al. (1996). I chose the latter for its ease of implementation. The third code (EMWAVE), described by Wang and Tripp (1996), is another Yee algorithm in two ways. Firstly, the spatial derivatives are approximated by special optimized second-order divided differences, which lead to better numerical dispersion characteristics. Secondly, the computational mesh is terminated using Liao absorb­ing boundary conditions, which are described by Chew (1990). A major disadvantage of this code is the requirement that transmitters, scatterers and receivers be placed a large distance from the mesh boundaries, so that the region of interest must be padded with extra cells (Wang and Tripp, 1996). In contrast, the PML absorbing boundary conditions have been shown to perform well in cases where the transmit­ter is placed close to the mesh boundary (Chen et al., 1996). An optimal solution might incorporate PML absorbing boundary conditions and the optimized second- 22 The second code, written by Dr. Cynthia Furse and others in the Department of Electrical Engineering, University of Utah, is an implementation of the Vee (1966) method using second-order accurate central difference approximations to the tempo­ral and spatial derivatives in Maxwell's equations. Derivation of the time-stepping equations is discussed in several texts (e.g., Taflove and Umashankar, 1989). Both electric conduction and displacement currents are included in the model. This code was originally developed for models contained within a region of free-space , and used Mur (1981) second-order absorbing boundary conditions (ABC) on the mesh bound­aries. This ABC is inappropriate for simulating broad-band EM responses of targets in dispersive media, since it is based upon a wave-equation in which phase velocity is assumed constant. An accurate ABC for lossy host materials was required. Attention has recently been focussed on Perfectly Matched Layer (PML) absorbing boundary conditions, which reportedly lead to superior accuracy (Berenger, 1994; Chew and Weedon, 1994; Chen et aI., 1996). I decided to add a PML ABC for lossy media to the Vee algorithm implemented by Furse et al. (pers. comm.). Differing formulations of the PML for lossy media have been published by Fang and Wu (1995) and by Chen et al. (1996). I chose the latter for its ease of implementation. full wavefield numerical solution to Maxwell's equations, incorporating both electric conduction and displacement currents. The code differs from the conventional Vee aI., order finite difference scheme described by Wang and Tripp (1996). However, the implementation of this algorithm would be a major undertaking. Therefore, I con­centrated upon development of the standard Yee FDTD modeling code using PML absorbing boundary conditions, rather than adapting the EMWAVE code. In all of these modeling codes, we impose the boundary condition Eten on the tangential electric fields at the surface of the PEC. When the fields are dis­cretized using the Yee method, continuity of the normal component of magnetic induction across this interface is automatically satisfied, as these components are calculated using values of Eton on the surface. Wang- Hohmann), FJ (Furse-Johnson) and WT (Wang-Tripp) codes. addition, we coded an analytic solution for the diffusive impulse response of a perfectly conducting half-plane in a conducting host, as derived by Weidelt (pers. comm.). This solution is discussed in Appendix C. The half-plane model can be used to simulate the EM response obtained near the edge of a large, plate-like body of massive sulphide min­eralization, as found in the Sudbury Basin and Kambalda nickel provinces. It is also a useful 2.5D model for use in validating the model responses calculated using the FDTD codes. order to verify that the electromagnetic responses of perfectly conducting bodies can be accurately estimated using the FDTD method, I first made a comparison between results calculated using the Weidelt half-plane solution and the WH FDTD algorithm. The model used in this comparison (Figure 2.7) consisted of a perfectly conducting half-plane lying within a 0.001 S/m medium. The magnetic transmitter dipole was oriented in the horizontal direction perpendicular to the edge, and placed 50 m above and to the left of the half-plane, in order to achieve good coupling with the target. The transmitter waveform consisted of a unit magnetic moment terminated at time t 0. A series of receivers measuring the vertical magnetic field were placed 23 order finite difference scheme described by Wang and Tripp (1996). However, the Etan = 0, (2.21) B Etan For simplicity of exposition, we denote these algorithms as the WH (Wang­Hohmann), In half­plane FDTD codes. In Sim target. The transmitter waveform consisted of a unit magnetic moment terminated at time t = O. A series of receivers measuring the vertical magnetic field were placed z=-50m -» z=-40m "^ receivers 0m 50m PEChalfplane-rj = 0.001 S/m z=50m to the right and perpendicular to the page. half-plane, before updating the magnetic field at each time step. The numerical results obtained using a cell size of 1 m within the region of interest agree well with the closed-form result, particularly at late times (Figure 2.8). The z component of the magnetic field response was used in this comparison because the whole-space contribution is equal to zero for this set of receivers, giving a response which is entirely due to the half-plane. The discrepancy between the analytic and numerical results decreased with decreasing cell size. These results give us confidence that the analytic result is correct, and that the diffusive impulse response of a perfectly conducting half-plane can be accurately computed using the FDTD approach. Since the comparison suggests that the analytic solution is correct in the diffusion regime, and an accurate response for a PEC can be obtained using an FDTD solution to the Maxwell's equations in which displacement currents are neglected, we can now test the wavefield FDTD modeling algorithms against the diffusive half-plane model. The Yee algorithm has been exhaustively tested for models involving 24 transmitter x .• -_ recei.v ers Y~ •••••• z=Om • SOm • ••••• (PJ E=e 0 h.0a0l1fp lSan/me ---./' ••• • Figure 2.7: Cross section of the PEC half-plane model. The conductor extends in­finitely below the transmitter, recording the time derivative of the magnetic field intensity (i.e., the impulse response). The perfectly conducting half-plane was incorporated in the FDTD model by explicitly setting the electric field components to zero at nodes lying within the half­plane, m within the of interest agree well the closed­form can be accurately computed using the FDTD approach. diffusion full half­plane Vee 0.01 Time= 10"6 seconds Figure 2.8: Comparison of the analytic result for a perfectly conducting half-plane explicitly at each time step. The calculated impulse responses are compared at delay 10- 6, 10- 5 10- 4 scattering from perfectly conducting objects in the wave-propagation regime (a survey of the relevant literature is provided by Shlager and Schneider, 1995). We must now examine its performance in the FJ implementation for scattering in the diffusion regime. The strategy is thus to rigorously test the code for both end regimes, as a necessary condition for correctness for all frequency regimes. The full wavefield FJ code was used to simulate the response of a perfectly conduct­ing plate measuring 200 m in length and 150 m in width (Figure 2.9) within a medium of conductivity 0.001 S/m. The finite difference mesh consisted of 100 x 150 x 150 cubic cells of dimension Ax 2 m. A y-directed magnetic dipole transmitter was placed 100 m above the surface of the plate and 50 m from its edge, with a vertical line of receivers placed below the transmitter. The transmitter waveform consisted of a Gaussian pulse whose amplitude spectrum decreased to half of its peak value at 25 Time = 10-6 seconds 0.Q11~------'------=::::=='===~-=---~_-'-_'--1 -- -- - - 0.005 - - Weidelt - Wang-Johnson -~.~0--4--~3~0~---2~0~~-?tO~=.==0=====tCO~--2LO----3~0----4~0---==50 OX 10 TIme = 10-5 seconds - - - "-"' -6~~~~~~~~-~-~-~-~-~ -40 -5 30 -20 - to Time ~ t 04 se~<i,nds 20 30 -1 x 10 40 50 -1.2 -1.4.~----:';:----::';:----=:;;=:==::::::=---L---'------'---- -40 -30 -20 -10 0 10 20 30 40 50 Z Coordinate (m) due to Weidelt (pers. comm.) with the numerical solution obtained using the FDTD algorithm of Wang and Hohmann (1993) with PEC boundary conditions enforced times of 6 , 5 and seconds. F J diffusion F J ~x = of a Gaussian pulse whose amplitude spectrum decreased to half of its peak value at 100 150 cells Mesh boundary 150m the closed-form impulse response solution with the transmitter waveform used in the FDTD simulation. The host resistivity and transmitter frequency range of the model were chosen so that the diffusive assumption did not contribute significant error to the analytic result used in the comparison. Comparison of the vertical magnetic field responses (Figure 2.10) for the receiver opposite the edge of the plate suggests that the finite-dimensional plate is an imperfect representation of an infinite half-plane. Hanneson (1981) found that, for a horizontal loop EM prospecting system with a loop separation L, a thin vertical plate can be considered effectively infinite when its length is > 3L and its width > 1.5L. Using this criterion, a transmitter-receiver separation of 100 m would require a plate of length 300 m and width 150 m for the simulation of a half-plane response. Unlike the WH code, which models the fields within a large volume with cells of increasing size, the FJ code uses a uniform mesh to model a restricted region, so that the dimensions of the plate which can be modeled using the FJ code are considerably smaller than those of the plate simulated by the WH code. The response of the plate modeled using the FJ algorithm decays more rapidly at late times than the response of the true half-plane. However, since early time and peak scattered field values for the two models are in good agreement, we suppose that the FDTD fields for a physically larger 26 1 ()() cells / Transmitter 150 cells Figure 2.9: Plate model used to represent a perfectly conducting half-plane. 300 kHz. The analogous half-plane model response was also calculated by convolving ;:: width;:: F J F J those of the plate simulated by the WH code. The response of the plate modeled using the F J algorithm decays more rapidly at late times than the response of the true half-plane. However, since early time and peak scattered field values for the two models are in good agreement, we suppose that the FDTD fields for a physically larger 27 0.5 0 | - 0 . 5 $ i<o -1 c a> •£ £-1.5 iel o •S -2 c O) CO E g-2.5 ••s CD > o •^ -3 -3.5 -4 X . - 10" t5 i \ l\ l\ t\ \ J \ / i i i i / i / i / i / if i i 11 ii ii i il i i ., - " - > • - / ^--~ / • • ^ ^ ^ "^ / ^^ / / / / Analytic solution FDTD i i i 0.2 0.4 0.6 0.8 1.2 x10 1.4 •5 Figure 2.10: Comparison of the vertical magnetic field response of the model depicted Yee a perfectly conducting body in the diffusion regime. the diffusive impulse response of a perfectly conducting body. The WT and FJ codes can be used to calculate full wavefield EM responses, and our numerical experiments suggest that the FJ code can accurately model PEC responses in the diffusion regime. Since the FJ code has been validated for free-space EM scattering from PEC targets, the necessary conditions for accuracy in the transitional frequency range have been established. The WT algorithm has a similar formulation to the FJ code, and is therefore expected to produce simulations of equal or better accuracy for PEC models. Moreover, the WT code has superior numerical dispersion characteristics (Wang and Tripp, 1996). However, the Liao absorbing boundary conditions depend upon the x 10-15 0.5 Olr------..... E-0.5 3- ~ -~ J!! <= :;; -1 .5 ~ .9 (j) -2 t .. -2.5 .\1 t: >" -3 I ' - / I I I I f I ,, ,, ,, ,­/ / --- - Analytic solution - - FOTD result ~o~ --~0.2~ --~0.4- -~~0.6- -~~0.8- -~----1~.2 --~1. 4 Time (sec) x 10-5 in Figure 2.9 with the analytic solution for the analogous half-plane model. plate would provide a better match to the half-plane response. Hence, we tentatively conclude that the Vee algorithm can be used to accurately estimate the response of In summary, the WH code provides an efficient and accurate means of calculating F J F J F J established. The WT algorithm has a similar formulation to the F J code, and is therefore expected to produce simulations of equal or better accuracy for PEC models. Moreover, the WT code has superior numerical dispersion characteristics (Wang and Tripp, 1996). However, the Liao absorbing boundary conditions depend upon the 28 establishment of a sufficiently large separation between the mesh boundaries and any transmitters, scatterers or receivers placed therein. 28 t ransmitt ers, The Radio Imaging Method or RIM (Stolarczyk, 1990) belongs to a class of cross-hole for use in coal mines as a means of detecting various mining hazards. However, attempts have recently been made to apply the technique to delineation of base metal mineralization (Thomson et al., 1992; Stolarczyk, 1992). A similar EM system has been developed in South Africa for this purpose (Wedepohl, 1993). The data reduc­tion and interpretation methods applied in each case are similar. order to examine the assumptions underlying RIM interpretation, let us con­sider er, e ix containing a harmonic vertical magnetic dipole of angular frequency u> and moment M = Me-^'z, where is the peak magnitude of the dipole moment and z is the vertical unit vector. The magnetic field components in polar coordinates are given by HT = -^eikr(l ikr) cos0, and He = -^eikr(l - ikr - k2r2) sin 0, (3.3) 47rr3 3. RADIO IMAGING METHOD MODEL STUDIES 3.1 The Radio Imaging Method (RIM) cross­hole imaging methods wherein variations in signal amplitude are used to reconstruct a distribution of attenuation rates within an image plane. RIM was originally devel­oped aI. , In a uniform medium of conductivity a, electric permittivity f and magnetic per­meability /1 w (3.1) M M ·k Hr --e' T(l - cosB, 27fr3 (3.2) (3.3) where k2 io\xu u2efx. In the far zone of the dipole (\kr\ > 1), we have Hr « 0, (3.4) and Following Stolarczyk (1992), let us examine RIM data reduction and interpretation theory. Commonly used interpretation techniques for the Radio Imaging Method assume a simplified ray-theoretic model of EM propagation based upon Equation 3.5, in which the amplitude of the signal in a homogeneous medium is given by , sin 9r = A0 -e'ar, 3.6 A0 9T ^sm(k). a/toe 3> Q LJ/J,(7 = 1/5, 5 r _ 1 r~1/2 in coal mining applications of RIM. The measured signal at a dipole receiver has a similar angular dependence and is given by ^ sinflrsinfl^-^ 3.8) 30 k2 = iaJ.Lw + W 2fJ.L. (Ikrl » (3 .5) 1992) , A -A ()T - or - o--e , r (3.6) where Ao denotes a "source strength", commonly determined by empirical methods (e.g., McGaughey and Stolarczyk, 1991), ()T is the polar angle with respect to the transmitter axis and a is the attenuation constant of the medium equal to S'm(k) . If the field behavior in the medium is entirely diffusive (a /Wf » 1), then the attenuation constant is equal to a = Jw~a = 6, (3.7) where 6 is the skin depth. The r- 1 dependence (or "spherical spreading term") in Equation 3.6 is replaced by a r- 1/ 2 cylindrical spreading term in the case of in-seam guided waves encountered The measured signal at a dipole receiver has a similar angular dependence and is given by A sin ()T sin () R - or A = 0 e , r (3.8) 6R subtended receiver axis. A = s i n ^ s i n ^ jL a dl 3-9) where JL is the line integral of the attenuation constant along the raypath linking the transmitter and receiver. r T= f adl - \n(A/A0) + ln(sin 9R) + ln(sin 9T) - ln(r). (3.10) signal amplitude, and the remaining terms are geometric corrections for dipole cou­pling the emf induced in the receiver coil. The emf is related to the magnetic field H by emf = IMJNOH • n, (3-11) a field. r and converting Equation 3.10 to decibels we have T = -0.115 [20\og(A/AR) + 20\og(AR/AQ) 3 12N - 20 log(sin 9R) 20 log(sin 9T) 20 log(r)]. Data reduction for tomographic imaging of the signal amplitude data is carried out using Equation 3.12. The instrument output has been designed around the RIM 31 where OR is the angle sub tended by the ray joining receiver and transmitter and the receiver axis. In an inhomogeneous medium the measured data are approximated by A -A sin Or sin OR - J adl - 0 e L r ' (3 .9) h adl a L Taking natural logarithms of Equation 3.9 and rearranging gives an expression for the net attenuation 1" given in Nepers: 1" = i = In(Ao) + In (OR) + In(sinOT) In(r). On the righthand side of this expression, the first term contains the normalized and spherical spreading. The raw data consist of the amplitude and phase of J1wNaH· n, (3.11) where N is the number of turns in the receiver coil, is the area of each turn and n is the unit normal vector for the receiver turns. The emf is thus proportional to the magnetic field . Therefore, the net attenuation 1" is calculated using the emf in place of the magnetic field in Equation 3.10. The signal level is recorded in decibels relative to a reference amplitude AR (usually one nanovolt). Using this form of normalization 1" - 20Iog(A/AR) + 20Iog(Ao) log (OR) - 20Iog(sinOr) + 20Iog(r)]. (3.12) data reduction and interpretation procedures, and extracting the type of normalized magnetic field data yielded by other EM equipment is impossible. For example, no information regarding current magnitudes in the transmitter coil is contained in the recorded data, and the "source strength" must be estimated from data acquired in a relatively homogeneous subdomain within the survey area. Although our numerical techniques can be used for RIM interpretation, the amount of information yielded by our calculations exceeds that in the RIM data. in a manner that maximizes the sensitivity of measurements to the geological pa­rameters which maps data uncertainties into parameter uncertainties (Kriesgshauser et al., 1996). However, we can also perform simple forward modeling experiments to estimate ideal, noise-free sensitivities for measurements taken within a given region, which indicate where data must be measured in order to obtain optimal resolution under ideal conditions. In this case, we compare the data generated by one model with those generated by a similar model, in which the parameters of interest have been perturbed. In performing the actual survey, we should collect the measurements which are most sensitive to perturbations in the parameters of interest. Whether or not these measurements are sufficiently sensitive for the parameters to be resolved to the desired level of accuracy is a question that can be answered only with reference to the predicted noise model and the instrument specifications. Consider the case where a plate-like body of massive sulfide mineralization has been intersected by several drill holes. We know the position of the body, and would 32 relatively homogeneous subdomain within the survey area. Although our numerical techniques can be used for RIM interpretation, the amount of information yielded by 3.2 Sensitivity of RIM measurements to geophysical model parameters The ability of geophysical measurements to resolve geological features is limited by such factors as noise and the physical parameters of the survey and geology. We can usually control only the survey design, subject to physical access restrictions and the limitations of our instruments. Therefore, it is important to design the survey of interest. In order to determine the parameter sensitivities accurately, we might perform a rigorous sensitivity analysis based on singular value decomposi­tion, Kriesgshiiuser aI., sufficiently to the predicted noise model and the instrument specifications. like to determine the location of its edges, for the purposes of ore reserve calculation and mine planning. It is possible to constrain the location of the edges by drilling more holes, although it may be more cost-effective to use geophysical methods. A cross-borehole electromagnetic survey can be carried out using a magnetic dipole transmitter operating at several different frequencies. For a given transmitter lo­cation, we should determine those locations where our measurements will be most sensitive to the position of the plate edge, in order to optimize the survey design. FJ and z) = 90m), S/m, 6.4e0. 150 120 ? 90 N c o W 30 . n Model B plate edge Model A plate edge 30 60 Horizontal distance x (m) 120 150 Figure 3.1: Perfectly conducting plate models. 33 like to determine the location of its edges, for the purposes of ore reserve calculation and mine planning. is possible to constrain the location of the edges by drilling more holes, although it may be more cost-effective to use geophysical methods. A cross-borehole electromagnetic survey can be carried out using a magnetic dipole We used the F J code to calculate the electromagnetic fields at 300 kHz, 500 kHz, 1 MHz and 2 MHz for two perfectly conducting plate models, in which the position of the upper plate edge differed by 30 m (Figure 3.1). In both cases, a horizontal magnetic dipole transmitter was placed at the position (x, z) (30m, gOm), the conductivity of the host rock was 0.001 81m, and the electric permittivity was 6.4£0· g 90 C .g ~ ~ Ui 60 o o Horizontal magnetic dipole transmitter - ,• , 90 34 The fields were recorded in a vertical plane perpendicular to the plates and containing the transmitter. The plates extended 50 m to either side of the survey plane. We will denote the model in which the upper edge of the plate lies at z = 60 m as A, and the model in which the edge lies at z = 90 m as B. A conventional RIM survey collects only amplitude data, which are usually ex­pressed in decibels relative to some arbitrarily chosen signal level. Therefore, it is useful to examine the model amplitude responses and sensitivities. For the purposes of visualization, the amplitude data were transformed as follows: F1 = 2 0 1 o g - ^ , (3.13) -fmax where is the field amplitude and Fmax is the maximum value. The values of the transformed field F 1 were clipped at a lower bound of Fl = -500 dB. The maximum value is actually meaningless, since field components computed within 10 cells of the source are inaccurate. This transformation, although somewhat arbitrary, facilitates image display. diffracted a shadow portion of the plate. Two sets of model data were used to calculate the sensitivity of the measurements to the position of the plate's upper edge. The observed field can be expressed as a sum of a primary field component, due to the dipole in a homogeneous whole-space, and a secondary field component due to scattering from the target. Since we want to examine the difference in response between the two models, it is useful to express the field perturbation due to the change in the plate edge position in decibels relative to the unperturbed model fields, so that FA - FB (* u) F = 2 0 1 o g - 7 T T 7 , ^3.14) 34 The fields were recorded in a vertical plane perpendicular to the plates and containing the transmitter. The plates extended 50 m to either side of the survey plane. We will denote the model in which the upper edge of the plate lies at z = 60 m as A, and the model in which the edge lies at z = 90 m as B. A conventional RIM survey collects only amplitude data, which are usually ex­pressed in decibels relative to some arbitrarily chosen signal level. Therefore, it is useful to examine the model amplitude responses and sensitivities. For the purposes of visualization, the amplitude data were transformed as follows: F F' 20 log z:;--, rmax (3.13) F field field F' F' 500 image display. The transformed amplitudes of the horizontal magnetic field components Hx at 2 MHz for model B (Figure 3.2) display a distinct shadow zone behind the plate, particularly in the lower portion of the target. Electromagnetic waves are diffracted by the edges of the plate, creating a zone of partial shadow behind the upper portion of the plate. Two sets of model data were used to calculate the sensitivity of the measurements to the position of the plate's upper edge. The observed field can be expressed as a sum of a primary field component, due to the dipole in a homogeneous whole-space, and a secondary field component due to scattering from the target. Since we want to examine the difference in response between the two models, it is useful to express the field perturbation due to the change in the plate edge position in decibels relative to the unperturbed model fields, so that I FA - FBI F = 20 log FA + f ' (3. 14) 35 140 120 100 ¥ Elevation z ( 05 00 o o 40 20 - M _ ! i i 1_ 1 1 » -» * 0 -20 -40 -60 -80 -100 -120 -140 -160 -180 -200 20 100 x (m) 120 140 Figure 3.2: Horizontal magnetic field amplitudes for model B at 2 MHz expressed in relative where FA and FB are the amplitude responses of unperturbed model A and the e FA 0. observational error, F useful edge. F at 2 MHz (Figure 3.3) indicates the data resolve the position of the plate edge lie within a wedge-shaped zone radiating from the top edge of the plate in model B, and lying on the far side of the plate relative to the transmitter. If the transmitter frequency is lowered to 500 kHz, then the zone of maximum sensitivity becomes more focussed. For instance, the region defined by the -5 dB sensitivity contour at 2 MHz (Figure 3.3) remains open at the lower edge of the figure, whilst the contour at 500 kHz (Figure 3.4) only extends about 30 m horizontally from the plate. The required sensitivity threshold of the instrument decreases more rapidly with distance from the target, and optimal positioning of receivers becomes more important as the frequency 3.) 0 -20 -40 -60 -80 -100 -120 -140 -100 -180 20 40 60 80 100 -200 Horizontal distance )( (m) decibels relat ive to the maximum value. The minimum has been clipped at -500 dB. The diffractive nature of the field behavior is evident. perturbed model B respectively, and € is a small number chosen to prevent F from becoming undefined when = O. For an instrument of given dynamic range and an estimated level of observational error, we can use maps of F to decide whether a given measurement is likely to yield useful information regarding the vertical position of the plate edge. The map of F at 2 MHz (Figure 3.3) indicates that the data that best resolve the position of the plate edge lie within a wedge-shaped zone radiating from the top edge of the plate in model B , and lying on the far side of the plate relative to the transmitter. If the transmitter frequency is lowered to 500 kHz, then the zone of maximum sensitivity becomes more fo cussed. For instance, the region defined by the -5 dB sensitivity contour at 2 MHz (Figure 3.3) remains open at the lower edge of the fi gure, whilst the contour at 500 kHz (Figure 3.4) only extends about 30 m horizontally from the plate. The required sensitivity threshold of the instrument decreases more rapidly with distance from the target, and optimal positioning of receivers becomes more important as the frequency --30 the unperturbed model response at 2 MHz, where the minimum has been clipped at -30 dB. Contours at -5, -10, -15 and -20 dB are shown. The region of maximum sensitivity forms a wedge terminating at the top edge of the perturbed plate. receivers (compare with Figure 3.3). faulting, this feature. Using the FJ code, the responses of two models were computed; one in which a plate of 75 m vertical extent is excited by a vertical magnetic dipole, and another in which the top 30 m of the plate is displaced horizontally 20 m away from the transmitter (Figure 3.6). As in the previous example, the width of the plate was 100 m. Amplitude data were modeled at 2 MHz and 500 kHz, and the sensitivities in a 36 o -5 iD -10 g :E-g. N. ~: I .>. -15 W -20 -25 20 40 60 80 100 120 140 Horizontal distance x (m) Figure 3.3: Horizontal magnetic field amplitude perturbations in decibels relative to is lowered. Replacing the 2 MHz horizontal magnetic field transmitter with a vertical mag­netic dipole and measuring the vertical magnetic field component results in a marked change in the sensitivity map (Figure 3.5). The area enclosed by the -10 dB sensitiv­ity contour is now greater than the area enclosed for the horizontal transmitters and Now consider an instance where a part of the original orebody has been displaced by fau lting, and we wish to determine a survey design that will allow us to resolve F J --30 140 Figure 3.4: Horizontal magnetic field amplitude perturbations in decibels relative to the unperturbed model response at 500 kHz, where the minimum has been clipped at -30 dB. Contours at -5, -10, -15, and -20 dB are shown. The region of maximum sensitivity is more compact than that shown in Figure 3.3. vertical slice containing the transmitter were calculated as in the previous example. For a transmitter frequency of 2 MHz, the data behind the lower portion of the orebody can be used to resolve the throw of the fault (Figure 3.7), whereas the data acquired at 500 kHz are comparatively insensitive to the offset (Figure 3.8). The lateral resolution of the survey deteriorates rapidly with decreasing frequency. 3.3 Confidence limits of parameter resolution Sensitivities of data to perturbations in a model parameter (i.e., Frechet deriva­tives) can be plotted as a function of measurement location, in order to determine the locations where the data are most sensitive to the model parameter of interest. However, these maps are of limited usefulness insofar as they do not consider mea­surement error. If a certain amount of noise is introduced, then the degree to which each datum can be used to resolve a parameter depends upon the magnitude of that datum relative to the error, as well as the magnitude of the Frechet derivative. We 37 o -5 iD -10 :[ :!l- '8N" ~¥' I -15 ~ > ~ w .w 0 '" -20 -25 20 40 60 80 100 120 140 Horizontal distance x (m) perturbat ions 10 , 1 5, sensit ivity e. , sensit ive 40 Horizontal distance x (m) -Figure 3.5: Vertical magnetic field amplitude perturbations in decibels relative to the unperturbed model response at 2 MHz, where the minimum has been clipped at -30 dB. Contours at -5, -10, -15, and -20 dB are shown. The zone of maximum sensitivity has a greater vertical width than in the previous example. then require a means of mapping data errors into parameter uncertainties. Assuming that perturbations in data are linearly related to perturbations in model parameters through the Jacobian matrix, we can use the Singular Value Decompo­sition (SVD) of the Jacobian to estimate the parameter covariance matrix from standard Martin, v confidence intervals for each parameter p3 are then given by 5pi = ±i /M (3.15) (Press et al., 1992, 690-691), where A\2 V is the x2 variable of v degrees of freedom for the desired confidence level. Using these estimates of the parameter uncertainties, the is discussed below. 38 o -5 CD -10 I ~ Nc ~" I ~0 -15 .. m ;!; ~ ~ 0 '" -20 -25 -30 60 80 100 120 140 Horizontal distance )( (m) magnet ic fi eld amplit ude unpert urbed 10 , vert ical perturbat ions t hrough mat rix, t he t he J acobian C the data st andard errors (Mar tin, 1971). For a model containing free parameters, Pj aI. , 691 ), 6 X~ X2 fr eedom Csing we can evaluate t he usefulness of data acquired at different locations in resolving a given parameter. An application of this method to cross-borehole EM survey design 39 150 E o •s m 60 30 Vertical 20m1 Horizontal shows the position of the upper 30 m of the orebody in each case. The dipole trans­mitter lies opposite the top edge of the orebody. Cross-borehole EM surveys usually consist of measurements taken at a number of receiver locations, using several different transmitter deployments. This is done to record data for several different illuminations of the target. We need to determine the locations at which we should place receivers in order to achieve optimal resolution of a parameter of interest. Consider the case where two boreholes have been drilled on ei­ther side of a large sheet of highly conductive sulfide mineralization (Figure 3.9). The position of the sheet edge is to be determined using cross-borehole electromagnetic data. Model data were generated using the analytic solution for the PEC half-plane due to Weidelt (pers. comm.) for two cases, in which the position of the half-plane edge differed by 20 m vertically (Figure 3.10) with the conductivity of the host ma­terial fixed at 0.001 S/m. In each case, an array of receivers with a spacing of 5 m 150 1------------------ 120 g 90 t Model A , , , , N , , c , , , , .2 VenicaJ magnetic ,, ,, ,, ,, :;; dipole transmitter , , 20m I , > " , , , , US ,, ,, ,, ,, , , , , Model B o o 30 60 90 120 150 HorizontaJ distance x (m) Figure 3.6: Faulted and intact perfectly conducting plate models. The dashed line ore body 81m. 40 140 120 100 LU 60 20 40 60 80 100 Horizontal distance x (m) 120 140 Figure 3.7: Vertical magnetic field amplitude perturbations in decibels relative to the unperturbed 1, field broad. model perturbation. was placed in one bore-hole. Magnetic field responses were calculated for a vertical magnetic dipole (VMD) transmitter operating at 300 kHz and placed in turn at 5 m spaced positions in the other bore-hole. These models emulate the tomographic mode of cross-borehole EM surveying in common use. Frechet derivatives were estimated by differencing the model results. For each receiver location, a Jacobian matrix was formed from the sensitivities corresponding to each transmitter location. Assuming standard errors of 10%, 15% and 20% for each measurement, the 99% confidence intervals for the position of the half-plane edge were estimated at each receiver location using the method described above (Figure 3.11). This mapping of data errors into parameter uncertainties shows that receivers placed well above the sheet edge do not record data that can be used to closely constrain our estimate of its position. It is obvious that data must be recorded at positions below the edge of the half-plane. Parameter uncertainties corresponding to receivers within the "shadow 40 120 100 80 I 60 N i 40 w 20 o -20 20 40 60 80 100 140 Horizontal distance )( (m) Vert ical pert urbations unpert urbed model response at 2 MHz, where the minimum has been clipped at -30 dB. Contours at -1 , -5, -10, -15, and -20 dB are shown. The region of maximum fi eld perturbation is broad . Measurements taken above the horizontal fault plane and behind the upper faulted plate segment contain little information regarding the model perturbation. operat ing spaced posit ions of cross-common estimated by differencing the model results. each receiver location, Jacobian formed from the sensitivities corresponding each transmitter location. Assuming standard errors of 10%, 15% and 20% for each measurement, the confidence intervals for the position of the half-plane edge were estimated at each receiver location using the method descri bed above (Figure 3.11 ). This mapping of data errors into parameter uncertainties shows that receivers placed well above the sheet edge do not record data that can be used to closely constrain our estimate of its position. It is obvious that data must be recorded at positions below the edge of the half-plane. Parameter uncertaint ies corresponding to receivers within the "shadow 41 -20 perturbations unperturbed model response at 500 kHz, where the minimum has been clipped at -30 1, for all receiver positions, and optimal positioning of the measurements becomes more critical (compare Figure 3.11 and Figure 3.12a). Estimates of the lateral position of the conductor (Figure 3.12b) have higher uncertainties by a factor of about 2 overall. Both vertical and lateral resolution are optimal for the receiver lying at the same height as the edge of the half-plane in the unperturbed model. In designing a survey to estimate the vertical and lateral position of the conductor to within a few meters, §: N C i w 20 40 60 80 100 120 140 Horizontal distance x (m) -11 120 100 80 :igii. 60 40 20 0 Figure 3.8: Vertical magnetic field amplitude pert urbations in decibels relative to the kHz , dB. Contours at -1 , -5, -10, -15, and -20 dB are shown. The range of field sensitivity is decreased relative to that shown in Figure 3.7 zone" of the plate are fairly uniform, increasing slightly as the signal magnitude decreases with depth. If both the vertical and lateral positions of the half-plane edge must be estimated from the data, then the parameter uncertainties are expected to be somewhat greater. Allowing both parameters to vary, and assuming uniform standard errors of 10%, 15% and 20% for each measurement, the 99% confidence intervals for the vertical and lateral positions of the edge were estimated (Figure 3.12) using the method described above. Our ability to resolve the vertical position of the conductor edge decreases 42 Borehole 1 Borehole 2 -•, Transmitter Figure 3.9: Perfectly conducting vertical half-plane lying between two vertical bore­holes. An array of VMD transmitters is placed in bore-hole 1, and an array of magnetic field receivers is placed in bore-hole 2. The half-plane extends to infinity downward and perpendicular to the plane containing the bore-holes. 42 2 I Receivers Figure 3.9: Perfectly conducting vertical half-plane lying between two vertical bore­holes. An array of VMD transmitters is placed in bore-hole 1, and an array of magnetic field receivers is placed in bore-hole 2. The half-plane extends to downward and perpendicular to the plane containing the bore-holes. Tx#l Tx#10 Tx#20 Bore-hole 2 Rx#l - z=70m VMD transmitters ^ 20 m I 50 m Model B Model A 50 m - PEC half-plane Rx#10 z=20m z=0m Rx#20 Rx#27 J Figure 3.10: Schematic diagram showing two half-plane models in which the vertical position of the edge is perturbed by 20 m. we should plan to measure data within the range of depths where the edge of the conductor is expected to lie. Data measured above the conductor have very poor sensitivity to the parameters of interest, especially in the presence of noise. The slight increase parameter uncertainty estimates at depths below the plate edge is more marked in this case than it was when only the vertical position was allowed to vary (compare Figures 3.11 and 3.12). The deleterious effect of noise upon resolution becomes more acute as the number of parameters to be estimated increases. There are two main problems with this method of estimating uncertainties, which both lead to estimates of the parameter uncertainties which are unrealistically small. Firstly, it relies on the assumption that the mappings between perturbations in data and parameters are linear. In electromagnetic models, this assumption usually holds only over a small range of parameter variations. Secondly, the use of uniform data relative errors implies that we can measure signals over an arbitrary dynamic range 43 Bore-hole 1 1 ~'MO","~'~ ~ 1 z.=70 m Magnetic field receivers Tx #I Rx #- 20m - z:;;:{) m 50m 50m Tx #20 - PEe Tx#27 I Rx #3_10: in t his c 1 83- §2 ' \ . \ V \ \ \ \ \ \ \ \ v \ * \ x- ^ \ x- v \ ^ N \ \ \ \ '\ x \ 'NA \ x. >• ' 15% - 20% • 15 25 30 3.11: level. 15 30 Figure 3.12: Estimated uncertainties for the position of the half-plane edge within the 99% confidence region, given uniform data standard errors of 10, 15, and 20%. (a) Uncertainties for the vertical position, (b) Uncertainties for the lateral position. the uncertainty ranges for the vertical position. 6;,-----~----~----~----~----~---- 5 , , , , , :[4 \ \ :r " ~ \ \ . , , " ~3 § I,2 If , , .. . ' , '.' ., ,, , , , \ ., ' , , '.' .,, , , , ' - 10% standard error ._ .. 15% - - 20% " ........ - - - - - - - - - - °0L-----~----~------~-----L----~L---~ 5 10 15 ~ ~ ~ Receiver number 44 Figure 3.11 : Estimated uncertainties for the height of the half-plane edge within the 99% confidence region, given uniform data standard errors of 10, 15, and 20%. The parameter is estimated using data recorded at each receiver for the set of 27 transmitters, and the uncertainty is plotted for each receiver number. Data recorded above the plate edge yield greater parameter uncertainty ranges. The parameter uncertainties associated with these data sets also increase more rapidly with noise 201r-----~----_.------~----~----_.------ , I , ' _ _ _ .-.e 10 ,. , , g 5 ... "_, - -- --'":"_-- 001-- ----sL----=='o=':'::::.':::==,;:=·= -=-=O~:~O::=~2S~=--J30 Receiver number 401 ,-----~----~----r===========~---, 1 -- 10% slandard error l -- -· 15% ~30 ._. .... - - 20% . 6 " " ~20 -" , "" ~ § 10 <'::--""" :> IL-~~' ' =~ =-=--~- - =~.~-- ~:-:=- -:=:;-:~. ---.~ °0 5 10 15 20 25 Receiver number position. The uncertainty ranges for the lateral position are approximately twice as large as 45 3.4 Optimization of transmitter arrays survey, Cherkaeva and Tripp (1996) have shown that it is possible to find an applied current distribution which maximizes the sensitivity of measurements to a geologic feature of interest. The method derived by these authors is briefly outlined below, with particular application to the cross-borehole EM problem. A simple but practical example is used to illustrate the method. Consider a set of electromagnetic measurements taken in the vicinity of an array Let each transmitter have a moment of magnitude u>i, so that Faw = d, where Fa is an operator which maps the vector of transmitter moments w = [toi, w2, w3 ... wN] with as the set of weights. + 5 F*+*w = d\ where F°+5 is the perturbed operator which maps the vector of moments w onto the perturbed data ds. - d 5 . I2, 45 to that level of accuracy. In practice, geophysical instruments have a limited dynamic range. Given sufficient a priori information concerning the geological setting of a geo­electric possihle of dipole transmitters and a given distribution of physical properties denoted as a. Wi, so that (3.16) FU = WI, W2, W3 . .. W N J onto the vector of data d. We can also consider the measured data as a weighted sum of responses to a set of unit dipoles placed at each transmitter loca­tion, w If the distribution of physical properties is perturbed, then we can write a second operator equation for the perturbed geophysical model denoted a + () as (3.17) FU+o dO. One measure of the sensitivity of the measured data to the perturbation is the L2 norm of the data residual d - dO. We denote this quantity as [2 , which can be considered a measure of the information content of the data with respect to the T2 = (F°+'-F°)w . (3.i8) w / . The measured response scales with w, and we wish to find the geometry of the / ||w||2 = Mi = max = max (Fa+5 - Fa)w . (3.19) w w : ||w||2 = 1 w is the eigenvector corresponding to the maximum eigenvalue Mt. An example of how this technique might be applied to optimization of cross-borehole EM survey resolution is discussed below. 3.10). is possible to form a weighted sum of the data recorded at each receiver the SVD of the impedance matrix, and the sensitivities of the resulting data are compared with the sensitivities corresponding to a uniform weighting (Figure 3.13). The maximum sensitivity is improved by about 20% by using optimal weightings, as opposed to a uniform weighting. The optimized weightings have the effect of emphasizing those transmitters for which the measured data have the greatest sensitivity to the model parameter of interest. Conversely, data that are unimportant are deemphasized, so that there is less opportunity for errors in these measurements to contribute uncertainty to the estimation of the model parameter. Thus, the distribution of weightings shown in Figure 3.13b has a shape similar to that of the optimal sensitivity distribution shown 46 perturbed model parameter. Thus we write (3.18) We would like to find the set of moments which maximizes the magnitude of I. optimal transmitter moments, so we optimize I subject to IIwl12 = 1, where MI = I(w) = 11(F'7H - puHI . w w: IIwll2 (3.19) This is an eigenvalue problem, wherein the optimal set of transmitter moments MI' Consider the problem of detecting the edge of a conductive sheet discussed pre­viously. We will again use the Weidelt code to generated the model data. Data are recorded at each receiver location using a number of different transmitters (Fig­ure It location which has optimal sensitivity to the parameter of interest. Considering only amplitude data, the optimal weights for this transmitter array were estimated us­ing 3.13) . . x10" .£•0.8 g0.6 CO CD 3 o en §I 1^ - ^ ' - 1 - / N / / \ \ / / • [ • 1 - Optimized weighting Uniform weighting - • 10 15 20 Receiver number 20 Figure 3.13: Optimization of sensitivities for a weighted transmitter array, (a) Sensi­tivities arrays, The Figure 3.10. 47 1 X 10-' .~O. 6 , - - ~ , , , , - i 0.6 , , • "5 0.4 , -0 • ~0.2 0 0 5 10 15 20 25 30 Receiver number 0.4 '" + + + + ~ 0.3 + + ~ + + .~ + + 3: 0.2 + + ".g + + + 8 0.1 + + + + + + + + + + + + 0 0 5 10 15 20 25 30 Transmitter number t ransmitter array. Sensi­t ivit ies of data for the uniformly and optimally weighted transmitter arrays. (b) The set of optimal weights applied to each transmitter. The model geometry is shown in in Figure 3.13a. CONCLUSIONS might include topics ranging from instrumentation and survey design to dielectric host. Given the range of frequencies and host-rock physical properties should include both diffusion and wave propagation effects. Therefore, a numerical modeling algorithm was adapted to treat the general case of full wavefield propagation. I chose a commonly used finite-difference time-domain solution to Maxwell's equations for its simplicity and demonstrated ability to model geometrically complex structures. in a conducting host material is a difficult numerical problem. Numerical instability results when attempts are made to explicitly model the behavior of fields within such targets. However, when the simplifying assumption of infinite target conductivity is made, accurate model responses can be calculated using FDTD methods. Accuracy of the code in the limiting cases of diffusion and wave propagation is a necessary condition for correctness in the transitional regime. The accuracy of the FDTD method for metallic targets in free space has been thoroughly tested by other authors (e.g., Taflove and Umashankar, 1989). There is good agreement between the analytic solution for the diffusive impulse response of a perfectly conducting half-plane and the numerical solution obtained using the modified Du Fort-Frankel method of Wang and Hohmann (1993) for a very large plate. The full wavefield FDTD algorithms 4. CONCLUSIONS A study of cross-borehole EM methods for the delineation of massive sulfide min­eralization imaging and interpretation. I have chosen to address the application of numerical modeling to the design of cross-borehole EM surveys. Selection and design of accurate numerical modeling techniques is a prerequisite for survey design and for other steps in the interpretation process. The class of models studied in this thesis is characterized by highly conductive bodies lying within a lossy commonly encountered in the hardrock mining environment, an accurate simulation Simulation of broad-band electromagnetic responses of highly conductive objects (e.g., Taflove and Umashankar, 1989). There is good agreement between the analytic the numerical solution obtained using the modified Du Fort-Frankel method of Wang and Hohmann (1993) for a very large plate. The full wavefield FDTD algorithms 49 model a region of comparitively restricted extent, so that it is prohibitively expensive to model an effectively infinite plate. However, comparison of the analytic half-plane model response with the full-wavefield FDTD result for a plate suggests that good agreement could be achieved for a sufficiently large plate. Our numerical experiments in the diffusion regime suggest that the full wavefield code is adequate for modeling transitional regime responses. of the computational mesh. This problem is especially acute in the case of a highly resistive host medium, where attenuation of electromagnetic waves is minimal so that spurious reflections from the mesh boundary can be large. Wang and Tripp (1996) successfully applied a modified Liao absorbing boundary condition (ABC) on the mesh boundary. However, these authors caution that it is necessary to place the source, receivers and scatterers a large distance from the boundary when the Liao ABC is used, creating a large computational overhead. I chose to implement perfectly matched layer (PML) absorbing boundary conditions, because Chen et al. (1996) have demonstrated that accurate results can be obtained in cases where the source is placed close to the PML mesh boundary, thus allowing smaller FDTD meshes to be used. The performance of the PML ABC depends upon an appropriate design of the perfectly matched media properties, and I applied a parameter estimation technique to choose an optimal PML profile. Using the optimized PML, I obtained accurate results for host media with conductivities as low as 0.001 S/m. 1988), which results from the difference between the numerical and physical wave velocities. The errors contributed by numerical dispersion grow with increasing prop­agation distance and frequency for a given mesh resolution. For the standard second-order accurate Yee algorithm used in this study, comparison with analytic solutions suggest that numerical dispersion errors are not significant when a 0.001 S/m medium is discretized at a 1 m mesh resolution, for reasonably large propagation distances (e.g., 100 m) and frequencies up to about 2 MHz. However, the mesh resolution demanded by more resistive host materials (e.g., 10~4 S/m) becomes prohibitively model a region of comparitively restricted extent, so that it is prohibitively expensive to model an effectively infinite plate. However, comparison of the analytic half-plane model response with the full-wavefield FDTD result for a plate suggests that good agreement could be achieved for a sufficiently large plate. Our numerical experiments A subsidiary problem encountered in FDTD modeling is the correct truncation 81m. An important source of error in FDTD modeling is numerical dispersion (Taflove, second­order 81m g. , 10- 4 81m) 50 fine as frequencies approach the range of interest (>100 kHz) for this study. If the mesh resolution is made sufficiently fine to practically eliminate numerical disper­sion, then the PML becomes physically thin and ceases to provide effective absorp­tion of outgoing energy. Wang and Tripp (1996) devised an optimized second-order differencing scheme which is claimed to have superior dispersion characteristics com­pared to the standard Yee algorithm. I conducted comparisons that demonstrated that the Wang-Tripp discretization produces significantly more accurate results in cases where the host medium conductivity reaches 10"4 S/m. However, the transmit­ters, scatterers and receivers must be placed a sufficient distance from the Liao ABC mesh boundaries. An optimal FDTD modeling code might incorporate the optimized second-order finite differences and PML absorbing boundary conditions, which allow the mesh boundaries to be placed closer to the region of interest. Yee a chosen parameter of a geoelectric model. By expressing the data perturbation as a proportion of the unperturbed model response, a map of the sensitivities can be created for a given transmitter deployment. I considered resolution of the vertical position of the upper edge of a vertical plate. For a vertical plate illuminated by a horizontal magnetic dipole aligned with its upper edge, the zone in which horizontal magnetic field measurements provide optimal resolution is wedge shaped, lying on the opposite side of the plate to the transmitter and tapering towards the upper edge. This zone of optimal receiver locations becomes smaller as the frequency is lowered - in order to attain a given level of sensitivity to the parameter of interest, receivers must be placed closer to the region in which the perturbation is taking place. Per­forming a similar experiment, in which I defined sensitivities to the lateral position of a faulted block of ore, I found a more pronounced degradation in lateral resolu­tion with decreasing frequency. The modeling code I have adapted can be used to determine receiver positions for optimal parameter resolution in arbitrarily complex models. Maps of the parameter sensitivities are useful insofar as they refer to noise-free line as frequencies approach the range of interest (>100 kHz) for this study. If the mesh resolution is made sufficiently line to practically eliminate numerical disper­sion, then the PML becomes physically thin and ceases to provide effective absorp­tion of outgoing energy. Wang and Tripp (1996) devised an optimized second-order differencing scheme which is claimed to have superior dispersion characteristics com­pared to the standard Vee algorithm. I conducted comparisons that demonstrated that the Wang-Tripp discretization produces signilicantly more accurate results in cases where the host medium conductivity reaches 10-4 S/m. However, the transmit­ters, scatterers and receivers must be placed a sufficient distance from the Liao ABC mesh boundaries. An optimal FDTD modeling code might incorporate the optimized second-order linite differences and PML absorbing boundary conditions, which allow the mesh boundaries to be placed closer to the region of interest. I have used my adaptation of the second-order Vee FDTD algorithm to demon­strate a method for determining the optimal placement of receivers in order to resolve lield _ deli ned with decreasing frequency. The modeling code I have adapted can be used to determine receiver positions for optimal parameter resolution in arbitrarily complex data. In order to determine the influence of noise on optimal survey design, I used the parameter covariance matrix estimated from the singular value decomposition (SVD) of the Jacobian to project data errors into parameter uncertainties, and thus find the transmitter-receiver combinations that resulted in minimal parameter uncertainty in the presence of noise. I considered the frequency domain response of a perfectly conducting vertical half-plane in a moderately resistive host. This model has practical application to various styles of mineralization, and can be computed using a closed-form expression. If the vertical position of the half-plane edge is to be resolved, then receivers placed below the edge provide almost uniformly good resolution, with uncertainties increasing slightly as the noise to signal ratio increases with depth. If estimates of both the vertical and lateral position are required, then the zone of optimal receiver placement becomes restricted to an interval below the edge of the plate. Uncertainties in lateral position are a factor of 2 greater than the uncertainties in vertical position for the 300 kHz signals considered in this model. 51 data. In order to determine the influence of noise on optimal survey design, I used the parameter covariance matrix estimated from the singular value decomposition (SVD) of the Jacobian to project data errors into parameter uncertainties, and thus find the transmitter-receiver combinations that resulted in minimal parameter uncertainty in the presence of noise. I considered the frequency domain response of a perfectly conducting vertical half-plane in a moderately resistive host. This model has practical application to various styles of mineralization, and can be computed using a closed­form expression. If the vertical position of the half-plane edge is to be resolved, then receivers placed below the edge provide almost uniformly good resolution, with uncertainties increasing slightly as the noise to signal ratio increases with depth. If estimates of both the vertical and lateral position are required, then the zone of optimal receiver placement becomes restricted to an interval below the edge of the plate. Uncertainties in lateral position are a factor of 2 greater than the uncertainties in vertical position for the 300 kHz signals considered in this model. When data are collected using a set of transmitters, the measurements at each receiver can be combined as a weighted sum such that the sensitivity of the resulting datum is optimized for a given parameter of the geoelectric model. I illustrated this technique using the PEC half-plane model, considering the problem of estimating the vertical position of the conductor edge. The set of optimal weights is equal to the elements of the eigenvector corresponding to the maximum eigenvalue of the impedance matrix that is formed by differencing the responses of perturbed and unperturbed models. The distribution of optimal weights for this problem is similar to the distribution of sensitivities for the resultant weighted summations. Thus the solution of the eigenvalue problem naturally leads to a set of weights which emphasizes the transmitters producing the greatest sensitivities to the parameter of interest. I have used several modeling techniques in this thesis. The analytic solution for the PEC half-plane (Weidelt, pers. comm.) proved useful for validating numerical solutions and for orienting studies of resolution and survey design. This code could form the basis of an inversion code for specialized application to problems involving large, plate-like conductors. The full wavefield FDTD modeling codes, which include 52 the WT and FJ programmes, appear to provide accurate solutions for PEC targets and are most suitable for 3D problems in which the fields exhibit wave-like behavior. Further work is required to compare the accuracy and efficiency of the WT and FJ codes. Based upon the examples presented by Chen et al. (1996), it is expected that the PML absorbing boundary conditions incorporated in the FJ code will allow the use of smaller finite-difference meshes than those allowed by the Liao ABC of the WT code. On the other hand, results of limited numerical experimentation have suggested that the optimized second-order differencing scheme of the WT code per­mits the use of larger cells. Therefore, it may be profitable to add PML absorbing boundary conditions to the WT code. In cases where purely diffusive behavior are to be modeled, the WH code provides a more efficient method of solution, owing to the conservation of memory afforded by expanding cell sizes and the efficiency with 52 and targets suitable behavior. F J condit ions F J which late-time responses can be calculated using expanding time-steps. However, simple comparisons of plane-wave propagation solutions using exact and diffusive ex­pressions for the wavenumbers imply that the diffusive approximation is unsuitable for most high resolution cross-borehole EM applications. expense. The PML technique of Berenger (1994) attempts to maximize accuracy whilst minimizing cost by embedding Maxwell's equations in an expanded set of equations. The additional degrees of freedom of this expanded set of equations are chosen by the modeler to accomplish two goals: 1. normal to the interface without giving rise to reflections. In fact, the without giving rise to reflections. The subsequent discussion in this section demonstrates one method in which such a PML may be constructed. The presentation follows Chen et al. (1996) and Chew and Weedon (1994). APPENDIX A PML ABSORBING BOUNDARY CONDITIONS Terminating a finite difference mesh by simulating a medium of infinite extent has always involved a compromise between computational accuracy and computational The interface between the conventional media, in which Maxwell's equations govern electromagnetic phenomena, and the PML is reflectionless. 2. The attenuation of a wave travelling into a PML can be controlled in the di­rection attenuation can be adjusted so that the wave amplitude at the outer boundary tends to zero, greatly facilitating truncation of the computational domain. Thus Maxwell's equations govern electromagnetic phenomena within the measure­ment region, which is surrounded by a PML "sponge zone" in which loss can be con­trolled 54 l {x, y, z} -• {xsx, ysy, zsz], (A.l) sx, sy sz this transformation, we can write an extension of Maxwell's equations in the frequency {e~'lwt dependence): V s x = iu/jH (A.2) V, x = icoeE + (A.3) Vs • eE = p (A.4) V, • //H = 0 (A.5) - - - - -L - - sx dx sy dy szdz a, e, p, p of conductivity, electric permittivity, magnetic permeability and charge density re­spectively. These equations reduce to the conventional Maxwell's equations when SX 1 Sy SZ. Let us now choose appropriate values of the coordinate stretching factors. Our first aim is to achieve a reflectionless interface between the PML and another medium, in which electric and magnetic fields are governed by the conventional Maxwell's equa­tions. This goal is facilitated by deriving a dispersion relation for the generalized 54 A.1 Generalized Maxwell's equations Let us introduce a transformation of Cartesian coordinates (A.l) where the "coordinate stretching factors" Sx> Sy and Sz are complex numbers. Under domain (e- iwt time dependence) : \7, E iW/lH \7, H -iw€E + aE \7, . €E = P (A.4) \7,'/lH =O (A.5) where 10 10 10 \7, =x-- +y-- +Z--. (A.6) SX ox Sy oy Sz oz The variables €, /l and correspond to the conventional physical parameters Sx = 1 = = Sz. A.2 PML boundary conditions 55 Maxwell's equations in a homogeneous medium. This dispersion relation can then be written in parametric form to give independent expressions for the scalar components of the wavenumber vector k = kxx + kyy + kzz. Values of the medium parameters in these expressions are chosen such that the phase matching conditions for the trans­verse components of the wavenumber vectors are met at the interface. This choice of parameters gives rise to a reflectionless boundary. Neglecting harmonic time dependence, a general Cartesian coordinate plane wave solution of Equations A.2 to A.5 has the form E = 'EIQe ikxx+ikyy+ikzz /A 71 and H = a0eikxX+ikyy+ik'z. (A.8) We derive a modified homogeneous Helmholtz equation in H by applying the stretched curl operator Vsx to Equation A.3 and substituting Equation A.2 on the righthand side: V s x V , x H = (a - iuje)Vs x E = (er - itoe)iu>p,H, so that V s x V s x H - K2H = 0, (A.9) where K2 = u2p,e + iup,a. Assuming solutions of the forms given in Equations A.7 and A.8, then Vs -> iks where ks = x^* + y^- + z^-, so that Equation A.9 becomes 3 y = (ks-ks)H-ks(ks.H). lA>iUj However Equation A.5 gives ks • H = 0, so that K2 = ks • ks, (A.ll) 55 Maxwell's equations in a h omogeneous medium . T hi s di spersion relation can then be written in parametric form to give independent express ions for t he scalar components of the wavenumber vector k = kx x + kyY + k,z. Valu es of the medium parameters in these expressions are chosen such t hat t he phase matching co nditions for the trans­verse components of the wave number vectors are met at t he interface. T his choi ce of parameters gives rise to a reflectionless boundary. Neg lec ting harmonic t im e depend ence, a general Cartesian coo rdina te plane wave solution of Equations A.2 to A.5 has the form (A .7 ) a nd (A .8) We derive a modified homogen eo us Helmh oltz equation in H by applying the stretched curl operator \7 , x to Equation A.3 and substituting Equation A.2 on t he righthand sid e: so that \7 , x \7 , x H = (0" - iWf) \7 , x E = (0" - iWf)iw J.LH , whe r e K,2 = w2 J.L f + iw J.LO". (A.9) Assuming solutions of t h e forms give n in Equations A.7 and A.8 , then \7 , ---t i k s where k = x & + Y""- + z &, so that Equatio n A.9 becomes s Sz: s~ Sz K,2 H = -k, x k s x H = (k, . k,) H - k,(k , . H ). (A.IO) However Equation A.5 gives k , . H = 0, so that K,2 = k s . k " ( A.ll) 1.2 i.2 12 w juje + i w / i a = - f -+ - £ - r - - | . S x Sy S2 kx = KSX 9 <f>, (A. 13) ky = KSy sin 9 sin <j> (A. 14) and kz = KSZ 9. (A. 15) can derive relations between the coordinate stretching factors in two media which matching the stretching factors for coordinates tangential to the interface gives us RTE = = R™. chosen arbitrarily. = Region 1 has the set of physical properties {o"i,ei,//i} and stretching factors {six,Siy,Siz}. Region 2 has the set of physical properties {a2,e2,p2} and stretch­ing factors {s2x, s2y, s2z}. The interface is illuminated by the plane wave given by Equations A.7 and A.8. We choose ax = a2,ex = e2 and p,i = p,2, so that KX = K2. In order to satisfy klx = k2x kly = k2y. Equations A.13 and A.14 give six sin 9i cos </>i = s2x sin 92 cos fa (A. 16) and 5iv sin 0i sin fa = s2y sin 92 sin fa. (A.17) 56 or (A.12) This is the equation of an ellipsoid, which can be parameterized by writing kx = KSx sin £I cos ¢, A.13) ky = KSy sin £I sin ¢ (A.14) k z = KSz cos £I. Using this set of expressions for the components of the wavenumber vector, we allow us to satisfy phase matching boundary conditions. Then we will show that RT E = 0 = RT M. The coordinate stretching factors for the normal direction can be Hence, consider an interface lying in the plane z constant, separating two me­dia. {O"I, EI, Mt} {Slx ,Sly, S l z }. {0"2 ,E2 , M2} {S2"> S2y, S2.}. S. 0"1 = 0"2, EI = E2