@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Earth, Ocean and Atmospheric Sciences, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Butler, Karl Edward"@en ; dcterms:issued "2009-03-20T18:39:36Z"@en, "1996"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """Seismoelectric effects are electromagnetic signals that arise when seismic waves stress earth materials. Their measurement is challenging because they are typically much weaker than the ambient electromagnetic noise. For this study, specialized instrumentation, field methods, and data processing techniques were employed to eliminate acquisition artifacts, and optimize signal-to-noise ratios. Experimental data were acquired which demonstrate clearly that seismoelectric effects can be measured in the field, and used to map shallow boundaries in porous sediments. Two types of seismoelectric signals were observed during field experiments at Ffaney, BC. The primary response was generated as compressional waves impinged upon a boundary between road fill and impermeable glacial till. Sledgehammer and blasting cap seismic sources positioned up to 7 m away from the boundary induced seismoelectric conversions with amplitudes of up to 1 mV/m, which were measured at the surface with grounded dipole receivers. This response arrived simultaneously at widely separated receivers, and was the dominant signal observed at near offsets. Recordings taken by receivers farther from the shotpoint were dominated by a second type of seismoelectric arrival which originated in the immediate vicinity of each dipole. Attempts to use seismoelectric conversions to map boundaries at other sites were not successful, but signals like the secondary arrivals at Ffaney were observed in all cases. Electrokinetic effects, arising from motion between the pore liquid and solid frame, are considered the most likely mechanism for the seismoelectric responses presented here. A simple conceptual model for seismically-induced electrokinetic effects can account for the two types of arrivals observed in field data. The model predicts that charge separations induced by a compressional wave in porous media produce electric fields that are observed (i) as the seismic wave passes by a receiver, and (ii) when the distribution of charge associated with the seismic wave is altered by a boundary or other inhomogeneity. In principle the boundary may separate regions with differing elastic or electrical/electrokinetic properties, permeabilities, or pore fluids since all of these factors influence charge transport. The experimental results from Haney, BC, support this model and other more elaborate theories for seismoelectric effects of electrokinetic origin."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/6270?expand=metadata"@en ; dcterms:extent "11833657 bytes"@en ; dc:format "application/pdf"@en ; skos:note "S E I S M O E L E C T R I C E F F E C T S OF E L E C T R O K F N E T I C ORIGIN by K A R L E D W A R D BUTLER B.Sc , Queen's University, 1987 M.Sc , The University of British Columbia, 1991 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF DOCTOR OF PHILOSOPHY in THE F A C U L T Y OF G R A D U A T E STUDIES (Department of Earth and Ocean Sciences) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A August 1996 © Karl Edward Butler, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of £ » v £ ^ ^Ut^CC4> The University of British Columbia Vancouver, Canada Date DE-6 (2/88) ii Abstract Seismoelectric effects are electromagnetic signals that arise when seismic waves stress earth materials. Their measurement is challenging because they are typically much weaker than the ambient electromagnetic noise. For this study, specialized instrumentation, field methods, and data processing techniques were employed to eliminate acquisition artifacts, and optimize signal-to-noise ratios. Experimental data were acquired which demonstrate clearly that seismoelectric effects can be measured in the field, and used to map shallow boundaries in porous sediments. Two types of seismoelectric signals were observed during field experiments at Ffaney, BC. The primary response was generated as compressional waves impinged upon a boundary between road fill and impermeable glacial till. Sledgehammer and blasting cap seismic sources positioned up to 7 m away from the boundary induced seismoelectric conversions with amplitudes of up to 1 mV/m, which were measured at the surface with grounded dipole receivers. This response arrived simultaneously at widely separated receivers, and was the dominant signal observed at near offsets. Recordings taken by receivers farther from the shotpoint were dominated by a second type of seismoelectric arrival which originated in the immediate vicinity of each dipole. Attempts to use seismoelectric conversions to map boundaries at other sites were not successful, but signals like the secondary arrivals at Ffaney were observed in all cases. Electrokinetic effects, arising from motion between the pore liquid and solid frame, are considered the most likely mechanism for the seismoelectric responses presented here. A simple conceptual model for seismically-induced electrokinetic effects can account for the two types of arrivals observed in field data. The model predicts that charge separations induced by a compressional wave in porous media produce electric fields that are observed (i) as the seismic wave passes by a receiver, and (ii) when the distribution of charge associated with the seismic wave is altered by a boundary or other inhomogeneity. In principle the boundary may separate regions with differing elastic or electrical/electrokinetic properties, permeabilities, or pore fluids since all of these factors influence charge transport. The experimental results from Haney, B C , support this model and other more elaborate theories for seismoelectric effects of electrokinetic origin. iii Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgements viii Chapter 1: Introduction 1 1.1 Background 3 1.2 Electrokineric Phenomena 9 1.2.1 The Electrical Double Layer 9 1.2.2 The Streaming Potential 12 1.2.3 Charge Accumulation and the Physical Nature of the Streaming Potential 17 1.3 Electrokinetic Effects as Coupled Flows 20 1.3.1 Theory of Coupled Flows 20 1.3.2 Sources of Streaming Potentials in the Earth 22 Chapter 2: Field Methods 26 2.1 Introduction 26 2.2 Instrumentation 26 2.2.1 Overview 26 2.2.2 Seismic Sources and Triggering 28 2.2.3 The Grounded Dipole Sensor 31 2.3 Sources of Noise 38 2.4 Field Methods 42 Chapter 3: Data Processing Methods 46 3.1 Introduction 46 3.2 Lowpass and Flighpass Filtering 47 3.3 Remote Reference Subtraction 52 3.4 Removal of Powerline Harmonic Noise 54 3.3.1 Block Subtraction 56 3.3.2 Sinusoid Subtraction 59 3.3.3 Comparison with Notch Filtering 67 3.3.4 Conclusions 70 Chapter 4: Field Experiments 71 4.1 Introduction 71 4.2 Experiments at Haney, B C 74 4.2.1 Site Description 74 4.2.2 Results and Discussion 75 4.2.3 The Seismoelectric Conversion Mechanism 87 4.2.4 Summary of the Haney Experiments 90 4.3 Experiments at Base Borden 91 4.4 Experiments on the Fraser Delta 94 4.4.1 The M O T Tower Site, Richmond 95 Chapter 5 : Modelling of Electrokinetic Effects Induced by Seismic Waves 103 5.1 Introduction 103 5.2 Seismic Wave Propagation in Porous Media 104 5.3 Conceptual Modelling of Seismically Induced Electrokinetic Effects 110 5.3.1 Charge Accumulations Induced by Compressional Waves I l l 5.3.2 Overview of the Conceptual Model 114 5.3.3 Detailed Analysis of the Conceptual Model 117 5.4 Quantitative Models 132 5.4.1 Frenkel's Suggestion 132 5.4.2 A Quasistatic Theory for Electrokinetic Coupling in the Elastic Wave Equations 134 5.4.3 Fully Coupled Seismoelectric Waves 138 5.4.4 Order of Magnitude Estimates for Electrokinetic Signal Strength 140 Chapter 6: Conclusions 144 References 148 Appendix A; Notes on Fortran Programs for Data Processing 154 Appendix B: Mathematical Detail for Sinusoid Subtraction 156 Appendix C: Corrections to Typographical Errors in the Literature 157 Appendix D: Physical Properties of a Fraser Delta Sand 159 V List of Tables 1.1 Examples of linear relations between fluxes and their direct driving forces 21 3.1 Main programs used to process seismoelectric records 46 4.1 Chronological list of field studies carried out (in whole or in part) for this research 73 4.2 Specifications of geophones used at the M O T Tower site 98 CI Differences between symbols used in this thesis and those used by Neev and Yeatts. ..157 D I Physical properties of a Fraser Delta sand 159 D2 Elastic and inertial constants for a Fraser Delta sand 160 List of Figures 1-1 The seismoelectric exploration method 1.2 The electrical double layer 1.3 Schematic illustration of streaming potentials generated by steady flow through a capillary tube 1.4 Example of a fluid flow that is not accompanied by a streaming potential 1.5 Illustration of a simple case for which the charge accumulation giving rise to a streaming potential can be calculated analytically 2.1 Circuit diagram for the T-box preamplifier 2.2 Equivalent circuit for a grounded dipole connected to the input (passive filter) stage of the T-box preamplifier 2.3 Amplitude and phase response of the T-box preamplifier 2.4 Intrinsic noise of the T-box and V-box amplifiers as measured at their outputs 2.5 Two useful source-receiver geometries for the acquisition of seismoelectric data with an eight channel recording system 3.1 Step responses for two types of highpass and lowpass digital filters 3.2 Noise removal by subtraction of a remote reference 3.3 Removal of powerline harmonics by block subtraction 3.4 Errors in amplitude and phase estimates as a function of the length of the estimation window in a record containing two harmonic signals 3.5 Removal of powerline harmonics by sinusoid subtraction 3.6 Comparison of three processing techniques for the suppression of harmonic noise 4.1 Eight-channel record showing seismic and electrical responses to a single sledgehammer blow on the surface at the Haney site 4.2 Cross-road traverse at Haney, B C 4.3 Seismoelectric response to detonation of a blasting cap in a borehole at 8 m depth 4.4 Seismoelectric response vs. shot depth in borehole 94A at Haney, B C 81 4.5 Seismoelectric response vs. shot depth in borehole 92B 82 4.6 Comparison of seismic and seismoelectric arrivals vs. dipole offset at Haney, B C 84 4.7 Linear and log-log plots of seismoelectric amplitude vs. dipole offset 85 4.8 Confirmation of the linearity of the seismoelectric response at Haney, B C 87 4.9 Radial symmetry of the seismoelectric response at Haney, B C 89 4.10 Plan view of the experimental layout at Base Borden 92 4.11 Comparison of seismic and seismoelectric arrivals vs. dipole offset at Base Borden.... 93 4.12 Location map for the MOT tower site on the Fraser Delta, Richmond, B C 95 4.13 Near-surface geology and experimental layout at the M O T tower site 97 4.14 Comparison of seismic and seismoelectric arrivals vs. dipole offset at the M O T site. 100 5.1 Conceptual modelling of seismically-induced electrokinetic effects 116 5.2 Spherical charged wavefronts emanating from a seismic point source at the surface of a homogeneous halfspace 119 5.3 Diagram showing the variables involved in the equation for the electrostatic potential due to a uniformly charged spherical cap 119 5.4 Model profiles for the electric field generated at the surface by reflection of a compressional seismic wave at five different depths 123 5.5 Diagram of the vertical dipole source used to approximate a pair of charged caps 125 5.6 Comparison of electric fields due to a vertical dipole and a pair of charged caps ........ 125 5.7 Comparison of measured and modelled profiles for seismoelectric amplitude vs. offset 128 5.8 Charge distribution induced by a critically refracted compressional wave 129 5.9 Electric field profiles generated by critically refracted charged wavefronts 131 Acknowledgements The field experiments at the core of this thesis would not have been possible without the assistance of past and present members of the Geophysical Instrumentation Group -particularly R. Don Russell, Anton Kepic, Michael Maxwell, and Barry Narod. We have shared many exciting and memorable experiences in pursuit of seismoelectric phenomena. I appreciate the expertise they provided on numerous trips to the field, and their many helpful suggestions. Special thanks to Anton for his advice on matters of instrumentation, for accompanying me on more trips to the field than anyone else, and for countless conversations - on matters of varying importance - in the office we shared during our student careers. My supervisor, Dr. R. Don Russell, was a steady source of support and good advice -always generous with his time, and mindful of my ideas. In particular, I would like to acknowledge his helpful suggestions on the topic of modelling electrokinetic effects. My respect for him as a scientist and advisor continues to grow after two postgraduate degrees under his supervision. I also acknowledge the guidance provided by my other committee members - Dr.'s G.K.C. Clarke, E .V. Jull, R J . Knight, and D.W. Oldenburg. Their dedication to reviewing research that was not directly related to their fields of study is greatly appreciated. I am indebted to many other people for assistance with field work. These include fellow graduate students Rob Luzitano, and Kevin Jarvis, and Dr. Jim Wright of Memorial University of Newfoundland, who joined our group during his 1993 sabbatical leave atUBC. The cooperation of the staff at the U B C Malcolm Knapp Research Forest in Haney , BC, and of their director Dr. Don Monro, is greatly appreciated. For access to sites at Base Borden, Ontario, I thank David Redman of the University of Waterloo. Dr. Jim Hunter and Dr. John Luternauer of the Geological Survey of Canada (GSC) were instrumental in helping me find sites suitable for experiments on the Fraser Delta. I appreciate the keen interest Dr. Hunter took in this research, and thank him and Dr. Susan Pullan for the loan of seismic equipment including an array of borehole hydrophones. Financial support for this research was provided by NSERC and five industrial partners in the form of a Cooperative Research and Development Grant (CRD-0094237) to R.D. Russell. The partners were BHP-Utah Mines Ltd., Cominco Ltd., C R A Exploration Pty. Ltd., Lamontagne Geophysics Ltd., and Placer-Dome Inc. Additional funding was provided by a NSERC Operating grant to R.D. Russell (OGP0000720), and NSERC and U B C scholarships to the author. Family and friends have been a great source of support. I thank my parents, Gordon and Nora Butler, for always encouraging me to take on whatever challenges interested me. Most of all, I thank my wife Becky Mills. She has shared in my excitement, has helped me keep difficulties in perspective, and has been a constant source of encouragement and inspiration throughout my years at UBC. Chapter 1 Introduction There are still old knots that are unrecorded, and so long as there are new purposes for rope, there will always be new knots to discover. C.W. Ashley, 1944, The Ashley Book of Knots The use of seismic sources to excite electromagnetic responses in the earth has been studied sporadically since the 1930's. However, the potential of such 'seismoelectric' methods for application in geological exploration has remained largely unknown. This is due in part to the fact that the various conversion mechanisms are only partially understood. More importantly, the body of empirical evidence showing that seismoelectric signals can be used to map geologic structures or provide other useful information is still relatively small. The major experimental difficulty is the identification of converted signals in the presence of much larger ambient electromagnetic noise. In recent years it has become possible to combat this problem through the use of recording systems with high dynamic range and digital processing techniques. The principal contribution of this thesis is a demonstration that conversions of seismic to electromagnetic energy can be measured in the field and used as a geophysical tool to map shallow boundaries. At our Haney test site, near Vancouver, Canada we have measured seismoelectric conversions that clearly originate at a shallow interface 1 to 3 m below the surface. The response exhibits the characteristics expected of a seismoelectric conversion, and can be used to map a subsurface boundary between road fill and glacial till. This study constitutes part of an ongoing effort by researchers at the University of British Columbia (UBC) to investigate the potential of seismoelectric methods in the search for minerals and other geologic targets. Early U B C studies, beginning in the late 1970's focused on the measurement of piezoelectric signals generated by the passage of a seismic waves 1 Chapter 1: Introduction 2 through quartz veins - a common environment for gold deposits. By 1990, the interests of the group had expanded to include the measurement of responses generated by an intriguing non-linear mechanism in sulphide ore bodies. In 1991, the remarkable seismoelectric conversions at Haney, B C , were discovered, inadvertently, by the author and his colleagues during a test of new instrumentation. Further investigations at the site revealed that piezoelectricity was an unlikely candidate for the conversion mechanism; an electrokinetic mechanism - involving relative motion between the solid and fluid components of a porous solid - was considered more likely. This thesis was motivated by the desire to determine the origin of the seismoelectric effect at Haney, and investigate its potential utility as a tool for mapping subsurface boundaries in porous sediments. My objectives remained much the same throughout this study. They were (i) to confirm that the seismoelectric responses measured in the field were real and indeed generated by the stressing of geologic materials; (ii) to develop data acquisition and processing techniques that would enable small seismoelectric signals to be resolved in the presence of much stronger ambient noise; (iii) to identify the geologic targets/materials responsible for observed signals; (iv) to develop a conceptual model to explain the origin of seismically-induced electrokinetic effects, and (v) to determine if the electrokinetic mechanism is a reasonable explanation for the signals I observed at Haney and elsewhere. The steps taken to satisfy each of these objectives are detailed in six chapters. In Chapter 1,1 introduce the concept of seismoelectric exploration, and discuss the more widely known mechanisms for seismoelectric conversion. The remainder of the chapter is a review of electrokinetic phenomena associated with steady flows of fluid through porous media. That topic has been studied for more than 100 years, and the concepts that have evolved are Chapter 1: Introduction 3 valuable for understanding the origin of electrokinetic effects induced by seismic waves. Chapter 2 is a discussion of the instrumentation and methods used to acquire seismoelectric data in the field. The topic deserves attention because there are no commercially available instruments or standardized techniques for that purpose. Most instrument development was carried out by other members of the U B C seismoelectric group, but careful tests of equipment and field methods were required to ensure that the measurements I made in the field were reliable, and free from acquisition artifacts. Data processing algorithms, developed to remove noise from field data are presented in Chapter 3. The most useful of these is a technique by which harmonic noise can be subtracted from a record without affecting the underlying signal; signal-to-noise ratios were improved by factors of up to 45 dB through this process. Chapter 4 contains the important results of field experiments carried at Haney, B C , and two other sites. These field trials confirm that seismoelectric effects can be measured in porous sediments where electrokinetic phenomena are likely to be induced. In Chapter 5,1 present a conceptual model for seismically-induced electrokinetic effects that accounts for the two types of arrivals observed in the field records. Quantitative models derived by other researchers are reviewed, and one of them is used to calculate order-of-magnitude estimates for seismoelectric conversions that are comparable to the signal strengths measured in the field. Finally, conclusions and recommendations for future research are given in Chapter 6. 1.1 BACKGROUND Seismoelectric exploration, as illustrated in Figure 1.1, involves the use of a seismic source, and electric or magnetic field receivers. In earth materials, electromagnetic signals travel much faster than seismic waves. As a result, the delay between the shot moment and reception of the electromagnetic response is due primarily to the time taken by the seismic wave to travel from the shotpoint to the target where the seismoelectric effect is produced. The product of this delay with the seismic velocity gives the distance from shotpoint to target, and data from different shotpoints can be used to locate the target. In addition, the magnitude Chapter 1: Introduction 4 of a seismoelectric response, coupled with an understanding of the mechanism that produced it, may provide information on the physical properties of the target material. Electromagnetic Sensors • • Target I -•• I Figure 1.1 Schematic illustration of the seismoelectric exploration method. A seismic wave produced by an explosion or other source propagates through the earth and stresses a seismoelectric target. The target converts some of the seismic energy to an electromagnetic response which is detected by electric or magnetic field receivers. There are at least four seismoelectric effects of interest in exploration geophysics: (i) the modulation, by seismic stress, of the resistivity of a volume of earth through which steady currents flow; (ii) seismically induced electrokinetic effects or streaming potentials; (iii) the piezoelectric effect; and (iv) highly non-linear processes that generate radio-frequency impulsive responses in sulphide-rich rocks. A summary of the experimental evidence for these four phenomena - with emphasis on the electrokinetic mechanism - is given blow. Measurements of resistivity modulation were originally reported by Blau and Statham (1936), and Thompson (1936, 1939). Further field studies were carried out by Long and Rivers (1975) and by Soviet researchers who called it the /-effect (e.g. Ivanov, 1949). Field experiments generally involved applying a constant voltage across a pair of electrodes inserted into the ground, and detonating an explosive charge nearby. The arrival of seismic waves Chapter 1: Introduction 5 beneath the electrodes caused current fluctuations in the electrode circuit which were attributed to changes in the resistivity of the volume of earth between the electrodes as it was alternately compressed and dilated by the incoming seismic pulse. Blau and Statham (1936) proposed that suitable arrays of electrodes could be used to discriminate against surface waves and enhance the detection of the vertically travelling seismic reflections most useful in seismic reflection prospecting. However, conventional seismometers proved better suited to the task and the idea was never applied in practice. No currents were injected into the ground during the experiments carried out for this thesis. Nonetheless, the possibility that resistivity changes might produce voltage fluctuations in areas where strong telluric currents flow through the ground must be kept in mind when making seismoelectric measurements in the field. Piezoelectric effects in rocks were first studied in the laboratory by Parkhomenko and Volarovich beginning in 1953. Parkhomenko (1971) reported that samples of vein quartz and pegmatitic quartz generally exhibited much larger piezoelectric effects than samples of other quartz-bearing rocks including quartzites, granites, gneisses, and sandstones. She concluded that vein and pegmatitic quartz were much more likely to have the preferential alignment of crystals required to generate a coherent piezoelectric polarization. Encouraged by the laboratory results, Soviet scientists began making field measurements of the piezoelectric effect in rocks during the late 1950's (e.g. Volarovich et a l , 1962; Neishtadt et al., 1972; Sobolev and Demin, 1980; Kondrashev, 1980). They carried out numerous case studies and claimed that ore bodies rich in the piezoelectric minerals quartz, and sphalerite (zinc sulphide) could be detected by the method of seismoelectric exploration described above. Quartz veins - a common environment for gold deposits - seem to have been the main target of interest. Interest in seismoelectric phenomena at U B C was boosted by the 1983 visit of G.A. Sobolev and V . M . Demin during which measurements were made at two Canadian mines using Soviet equipment (Sobolev at al., 1984). Subsequent field experiments by U B C researchers at sites in both Canada and Australia supported Soviet claims that quartz veins could be detected by this Chapter 1: Introduction 6 method (Butler, 1991; Russell et al, 1992). Furthermore, theoretical estimates of expected piezoelectric signal strengths (Russell and Barker, 1991) were found to be in reasonable agreement with the p.V/m to mV/m field strengths measured in the field. The discovery of unusual radio-frequency seismoelectric responses from sulphide ore bodies was made by Sobolev et al. (1980, 1982) in the late 1970's while performing piezoelectric surveys. The phenomenon, known as PRRER (pulsed radio wave range electromagnetic radiation) or RPE (the radio pulsed effect), differs from other seismoelectric effects in that it is highly non-linear. For example, the signals are generated only if the seismic energy exceeds a certain threshold; they also contain frequencies up to a few M H z while the seismic waves which generate them contain frequencies no higher than a few kHz. The existence of RPE, its non-linear nature, and potential for use in mineral exploration have been confirmed by field trials carried out in three sulphide mines by the U B C seismoelectric group (Maxwell et al., 1992; Kepic et al., 1995; Kepic, 1995). Theoretical models for the effect have been proposed by Sobolev et al. (1982), and by Kepic (1995) who reviewed several candidates. Kepic's preferred mechanism involves the opening of small cracks in the ore body, charge separation across the crack (triboelectricity), and subsequent electric discharge across the crack gap. Nonetheless, at this time, the potential mechanisms remain rather speculative. In 1939 and 1940, Ivanov presented convincing field measurements of a seismoelectric effect in sedimentary materials and proposed an electrokinetic mechanism for the conversion. He used a simple, elegant test based on signal polarity to demonstrate that the effect differed from the modulation of telluric currents by seismic stress (the /-effect) which had been discovered earlier. He suggested the names 'seismoelectric effect of the second kind' and 'E-effect' to emphasize that point. Electrokinetic effects arise because of the electric double layer that exists at a solid-liquid interface. The double layer consists of a layer of ions adsorbed on the solid matrix, and a parallel, diffuse layer of counter-ions in the pore fluid. At least part of the diffuse layer is free to move with the pore fluid. Thus, the flow of fluid relative to solid Chapter 1: Introduction 1 allows for the possibility of charge separation and the development of an (electrokinetic) electric field known as the streaming potential. Inspired by Ivanov's (1940) observations, Frenkel (1944) made the first attempt to explain seismically-induced electrokinetic effects in a quantitative fashion. He showed that the propagation of a F-wave through porous, elastic media is accompanied by relative motion between the solid and pore fluid, and that this relative motion can give rise to dynamic streaming potentials. Neev and Yeatts (1989) reached similar conclusions by modifying Biot's (1956) equations for seismic wave propagation in porous media to account for electrokinetic coupling at low frequencies. Thompson (1990) claimed that the effect could also be used as an exploration technique to detect permeable formations, and oil/water or gas/water contacts in sedimentary rocks. More recent and comprehensive theoretical work (Pride, 1994; Pride and Haartsen, 1994; Haartsen, 1995) predicts that seismic body waves in porous media produce two types of seismoelectric effects through electrokinetic coupling: (i) non-radiating electromagnetic fields that are trapped within the seismic pulse and travel at the seismic velocity, and (ii) electromagnetic fields that are generated by seismic waves traversing boundaries where there is a change in elastic properties, permeability, or pore fluid chemistry (because of its effect on electrical properties such as the zeta potential). The seismoelectric responses recorded by Ivanov (1939, 1940) were not clearly linked to a subsurface interface. Rather, they were observed whenever a seismic wavefront arrived near the grounded dipole sensor. The delay between the moment of shot detonation and reception of an electromagnetic response depended on the shot-dipole separation rather than the distance from the shot to any particular interface or target. Since the time of Ivanov, there have been several other reports of seismoelectric effects observed in the field and attributed to electrokinetic phenomena. Martner and Sparks (1959) documented clear seismoelectric responses from the base of the seismic weathered layer. Broding et al. (1963) measured the response along a 35 m profile in a borehole and found that it peaked opposite a sandy loam/shale interface. Borehole measurements by Parkhomenko and Gaskarov (1971) showed Chapter 1: Introduction 8 that seismoelectric responses in limestone were consistently stronger than those observed in clays. Migunov (1987) reported that electrical signals were generated when seismic waves reached the boundaries of kimberlite pipes. Maxwell et al. (1992) found that piezoelectric responses from a quartz vein were accompanied by other responses originating in the surrounding sediments and host rock. Large scale field experiments have been carried out by Exxon researchers Thompson and Gist (1991, 1993). They reported that they were able to detect electrokinetic conversions from boundaries between impermeable rocks and permeable water-saturated sands at depths of up to at least 300 m. The seismoelectric measurements presented in this thesis, and in two related publications (Butler et al., 1994, 1996) are also attributed to electrokinetic phenomena. Parkhomenko (1971) has summarized laboratory studies of seismically-induced electrokinetic effects made on samples of porous rocks and soils by Soviet scientists in the 1950's and 60's. Most measurements were made using ultrasonic piezoelectric transducers to generate compressional waves, and electrodes attached to the sample faces to measure electric potential differences. Unlike piezoelectric effects - measured in the same way on samples of rocks such as vein quartz - the polarity of the effect observed in these porous samples exhibited no dependence on sample orientation. Parkhomenko reported that no responses were observed in rocks that were completely dry. As the water content was increased from zero the seismoelectric response increased rapidly. However, beyond a certain, relatively low, water content (which varied with rock type) the magnitude of the response remained essentially constant or even decreased as more water was added. The strong dependence on water content at low saturations is intuitively satisfying in that it is the motion of water in the electrical double layer, which forms in a thin zone close to the pore walls, that is responsible for electrokinetic effects. Once this layer has formed, the addition of more water to fill the pores evidently has little effect on the electrokinetic response. Gaskarov and Parkhomenko (1974) also investigated seismoelectric effects in rocks containing frozen pore water, and in Chapter 1: Introduction 9 rocks that had been saturated with gasoline and petroleum after having being dried to a constant weight. The measured effects in both cases were negligible compared to those observed in the same rocks with water in their pores. The authors suggested that this difference could be exploited to search for permafrost and oil/water contacts by measuring seismoelectric effects as a function of depth in a borehole. Electrokinetic phenomena have been studied extensively because of their importance in physical chemistry, and for the information they provide about flows of water in the earth. Before addressing the more complicated phenomenon of seismically-induced electrokinetic effects in Chapter 5, it is instructive to review the important concepts already known from study of the electrokinetics of steady flows. 1.2 ELECTROKINETIC PHENOMENA 1.2.1 The Electrical Double Layer Electrokinetic effects are the result of a coupling between fluid flow and electric current flow in a medium containing both solid and liquid components. They arise because of a charge distribution known as the electrical double layer that exists at a solid-liquid interface. The solid surface comprises one half of the double layer. It carries a surface charge which arises from the dissociation (ionization) of surface groups, the adsorption of ions from solution, or, in the case of clays, from crystal lattice defects (Hunter, 1987, p. 374; Israelachvili, 1991, p. 213). The other half of the electrical double layer consists of ions in the liquid, and several models have been proposed to account for the ion distribution. The currently accepted model is due largely to Stem (1924) (Ishido and Mizutani, 1981; Morgan et a l , 1989). Figure 1.2 schematically depicts the ion distribution as well as the electric potential profile in the double layer according to the Stem model. The solid surface is here assumed to be negatively charged. To preserve electroneutrality, the net charge in the liquid must be of 1 0 Figure 1.2. Stern model for the charge and electric potential distribution in the electrical double layer at a solid-liquid interface. In this example (as in most rock-water systems), the solid surface is negatively charged and the mobile counter-ions in the diffuse zone are positive. Chapter 1: Introduction 11 equal magnitude but opposite sign to the surface charge. Some of the cations in the liquid are adsorbed on the surface of the rock in a region known as the Helmholtz layer. At greater distances the cations are attracted to the surface by electrostatic forces only and spread out into a diffuse layer known as the Gouy-Chapman diffuse zone. The distribution of ions within this zone is governed by the Poisson-Boltzman equation which accounts for the balance between electrostatic, and thermal diffusional forces. For a symmetric electrolyte, the solution to this equation in one dimension (i.e. for a broad planar interface) is well-known and yields an electric potential profile that decays approximately exponentially with distance from the Helmholtz layer (e.g. see Mitchell, 1993 or Pride and Morgan, 1991). Far from the interface, in the bulk pore fluid, the potential and net charge density are essentially zero. A measure of the double layer thickness, called the Debye length and commonly denoted K \" 1 , is given by where £/is the permittivity of the liquid, k is Boltzman's constant (1.38xl0~23 J/K), Tis temperature in Kelvin, v is the ionic valence, e is the electronic charge (1.602xl0\"19 C), and n° is the electrolyte concentration in ions/m3. Equation (1.1) shows that the thickness varies inversely with the valence and the square root of the concentration, and directly with the square root of temperature and permittivity. For monovalent electrolytes with concentrations in the range 0.001 to 0.1 mol/litre (typical concentrations for aqueous solutions in soils), the thickness varies from about 10 to 1 nm. These values are much smaller than pore radii in the majority of rocks and soils with the exception of materials like clays (Mitchell, 1993, p. 117; Pride and Morgan, 1991). Electrokinetic effects arise because only part of the electrical double layer is free to move with the liquid relative to the solid. The shear plane separating immobile ions that are tightly bound to the solid from mobile cations in the diffuse layer lies very close to the outer (1.1) / Chapter 1: Introduction 12 limit of the Helmholtz layer. As discussed below, the electric potential at the shear plane, called the zeta potential £ , plays a central role in determining the magnitude of electrokinetic effects. Most crustal rocks and rock minerals have a negative surface charge and zeta potential when in contact with natural pore liquids. The mobile charge in the diffuse layer therefore tends to be positive. Laboratory measurements made by Ishido and Mizutani (1981) and Morgan et al. (1989) indicate that zeta potentials of -10 to -100 mV are typical. Mitchell (1993, p. 270) indicates that values between 0 and -50 mV are common for clays, the lowest magnitudes being associated with high pore water salinities. 1.2.2 The Streaming Potential When liquid flows through a porous medium under the influence of a pressure gradient, charges in the diffuse part of the electrical double layer are displaced in the direction of flow, and an electrical potential difference develops between the ends of the sample. This effect, known as the streaming potential, is the main electrokinetic effect of interest in this study. The classic equation relating the streaming potential to a pressure gradient is derived for the case of flow through a capillary tube as illustrated in Figure 1.3(a). Here, the pressure gradient along the tube Ap/L is provided by the difference in hydraulic head between two reservoirs. Figure 1.3(b) is an enlargement of part of the tube and shows schematically the distribution of excess positive charge (cations) in the diffuse layer near the tube walls, as well as the parabolic velocity profile given by Poiseuille's equation for the flow of a viscous fluid through a circular tube (e.g. Hunter, 1993, p. 105). The amount of charge transported in the liquid depends not only on the concentration of excess cations in the diffuse layer, but also on the thickness of that zone. When flow begins, positively charged liquid exits one end of the tube while neutral electrolyte enters the other. This results in the accumulation of excess positive charge at the outlet, and of excess negative charge (surface charge) at the inlet. A potential difference A y known as the streaming potential therefore develops between the ends 13 (a) High Input Impedance Voltmeter Figure 1.3. (a) Schematic diagram of the type of apparatus used to measure streaming potentials in the laboratory, (b) Cross-section of a capillary tube showing the parabolic velocity profile characteristic of laminar (Poiseuille) flow, and the distribution of excess charge in the electric double layer along the tube walls. The thickness of the double layer is exaggerated for illustrative purposes. Chapter 1: Introduction 14 of the tube. Charge accumulation ceases, and A\\\\t stabilizes when the conduction current generated by the electric field balances the streaming current caused by convection of charge in the moving liquid. The apparatus in Figure 1.3(a) could also be used to observe an electrokinetic effect known as electroosmosis which is the converse of a streaming potential. In electroosmosis, a pressure gradient is induced by the application of an electric field along the capillary. If the two reservoirs were initially at the same level, then the application of an electric field would cause fluid to flow along the tube raising the level of one reservoir and reducing the level of the other. Fluid flow would cease when the back-pressure produced by the difference in hydraulic heads couterbalanced the viscous drag exerted by charged ions moving under the influence of the electric field. Returning to our analysis of the streaming potential, note that the conduction current is proportional to Ay and the streaming current is proportional to Ap . Therefore, under conditions of steady fluid flow (when the two currents must be equal and opposite), Ay must vary linearly with Ap. The constant of proportionality can be determined as outlined below. We begin by writing equations for the streaming current Is and conduction current Ic: where oy is the fluid conductivity, q(r) is the charge density, and v(r) is the fluid velocity given by Poiseuille's equation for laminar flow of a fluid of viscosity r\\ through a capillary of radius a. The assumption of laminar flow is justified provided flow rates remain low enough to prevent turbulence. Equation (1.3) (a statement of Ohm's law) is justified provided the conductance determining the conduction current is dominated by the bulk conductivity of the fluid. In some cases surface conductivity may contribute significantly to the conductance in (1.2) and Ic = 7ta2<7 Ay/ (1.3) L Chapter 1: Introduction 15 which case a correction factor (discussed below) must be added to O/. Using Poisson's equation to substitute for q(r), assuming the double layer is confined to a region near the wall very much thinner than the tube radius, and requiring that Is+Ic = 0 leads to the Helmholtz-Smoluchowski equation for the streaming potential ^ = ^ - (1.4) Ap r}of where C, is the zeta potential, and £/ is the liquid's permittivity. More detail on the derivation of equation (1.4) can be found in many texts (e.g. Hunter, 1993, p.249). Overbeek (1952, p. 204) derived equation (1.4) for a curved capillary of variable cross-section, and thereby showed that it is also applicable to the streaming potential developed by flow through a porous plug such as a soil or rock sample provided the assumptions discussed above are also valid for the tortuous flow paths through the sample (note that Overbeek used cgs units). It is interesting to note that the derivation makes no assumptions about the structure of the double layer except that it is thin compared to the tube and that charge becomes mobile along a shear surface where the potential is t,. Since £ is usually negative in earth materials, the potential gradient A\\|//L is normally in the opposite direction to the pressure gradient Ap/L. It is also significant that the streaming potential as given by equation (1.4) is independent of the capillary radius or pore size. The surface conductivity neglected above represents the enhanced conductivity of the diffuse zone with its excess charge. It tends to be most important at very low pore fluid salinities. If surface conductivity os contributes substantially to the conductance of the sample, then the streaming potentials predicted by equation (1.4) will be too large. A correction factor for the case of a capillary tube can be formulated by noting that the conduction current will be proportional to 7ta2of + 2naas , rather than Ka2o as in equation (1.3). As a result, equation (1.4) becomes Chapter 1: Introduction 16 Ay/ £fC (1.5) V 2os a The presence of the capillary radius a in equation (1.5) shows that pore size does influence the streaming potential when surface conduction is present. This is to be expected as the ratio of surface area to pore volume increases as pores become smaller. Although instructive, equation (1.5) is but only partly successful in accounting for the effect of surface conduction in soils and rocks (Overbeek, 1952, p. 206, 203). Another approach involves reformulating equations (1.2) and (1.3) to model the streaming and conduction currents through a porous plug rather than through a single capillary. In that case, the appropriate conductivity to use in equation (1.3) is the bulk conductivity of the sample o~ B which accounts for conduction through the fluid, along surfaces, and also through the solid grains themselves (the last could be important in the case of metallic ores). Also, equation (1.2) for the streaming current must be multiplied by two factors: (i) the porosity <|) to account for the fact that only a fraction of the plug's cross-sectional area is available for fluid flow, and (ii) an empirical factor c, (0 < c < 1) to account for reduced fluid flow due to tortuosity of the pore network. (Factor c would be 1 if the pore space consisted of an array of capillaries extending straight across the sample.) The resulting expression for the streaming potential across the porous sample is then ^ = c ^ . (1.6) Ap T]OB This equation is the same as that given by Ishido and Mizutani (1981) except they wrote c as a function of tortuosity* (c = 1/x 2 ) , and a B as a function of tortuosity, fluid conductivity, and surface conductivity. The quantity Ay/Ap is sometimes called the streaming potential cross-coupling coefficient. Laboratory measurements by Morgan et al. (1989) of streaming potentials * Tortuosity is commonly defined as the dimensionless ratio % = L'/L where L ' is the shortest distance that a fluid particle must travel on average to traverse a straight-line distance L through the porous medium. Hence, T > 1 by definition. Chapter 1: Introduction 17 generated by a variety of electrolytes flowing through crushed samples of granite, quartzite, and K-feldspar indicated that the coupling coefficient varied nearly linearly with fluid resistivity. Most measurements were made for fluid resistivities in the range of about 5-70 Qm (given as typical resistivities for natural surface waters), and the effect of surface conductivity was relatively small in their samples. For 1:1 electrolytes (such as NaCl), the ratio of the coupling coefficient to fluid resistivity was approximately -4 mV/(atm Qm), while for 2:1 electrolytes (such as CaCh) it was closer to -2 mV/(atm Qm). 1.2.3 Charge Accumulation and the Physical Nature of the Streaming Potential As mentioned above, the physical nature of the streaming potential is the development of regions of excess positive and negative charge. The flow of fluid carrying excess charge in the mobile part of the double layer constitutes a streaming current, but does not, by itself, create a streaming potential. This point can be illustrated by considering the flow of fluid through a donut-shaped ring as shown in Figure 1.4. The fluid flow could be generated by spinning the ring for some time and then suddenly stopping it; the fluid inside would continue to circulate relative to the walls of the ring and would gradually slow down due to viscous drag. By the symmetry of this experiment, the fluid velocity, and hence any electric potential generated by electrokinetic coupling, must be independent of angular position 8. Thus, the streaming potential, defined as A y =Y ( 0 2 ) - V ( 0 I ) (holding r constant) must be zero. Symmetry also dictates that the pressure must be independent of 9. Thus Ap = p(Q2)~ Pfii) is also zero and the linear dependence of Ay on Ap is maintained. However, this example shows that the existence of fluid flow and streaming current is not sufficient to generate a streaming potential. A steady streaming current must encounter a discontinuity (such as the ends of the sample in a laboratory measurement) in order to produce charge separation and an electric field. Figure 1.5 depicts a simple flow geometry for which it is straightforward to calculate Chapter 1: Introduction 18 the charge accumulation associated with a streaming potential. We consider an ideal porous plug in which the pores are straight capillaries extending from one end to the other. According to Gauss' Law, the surface integral of the net electric flux leaving one end of the sample must equal the total enclosed charge Q: §eE-dS=jjjqdV =Q . (1.7) s v(s) Now if the plug's cross-sectional area A is large compared to its length L then we can neglect the contribution of the fringing field to the total flux in the same way as we neglect the fringing field of a parallel plate capacitor. Then equation (1.7) simplifies to Q = A(emlEsol{l-^) + efE^) (1.8) where is porosity, and esol, Esol ,ef, and Ef are the permittivities and electric fields in the solid and fluid components. Since the tangential component of the electric field must be continuous across a boundary, Esol = Ef =E. Furthermore, since E = -A \\ j / / L, equation (1.8) reduces to Q = ^j^JA(esol(l-<$>) + £,<$>) . (1.9) showing that the charge accumulation is proportional to the streaming potential across the sample. We may substitute for A\\|/ using any of equations (1.4) to (1.6) to obtain the accumulated charge in terms of the pressure gradient Ap/L and fundamental material properties. For example, substitution of (1.5) yields Q -Ap • e , C A(£ml{l-Q>) + £f V . J ^ = V - ( X V p ) . (1.22) Chapter 1: Introduction 25 VX-Vp + XV2p. (1.23) Finally, using equation (1.19) to find V 2 p, and substituting the result into equation (1.23) yields The terms on the right hand side of equation (1.24) represent three source terms for conduction currents generated by steady fluid flow. Since conduction currents give rise to voltage drops, the terms in equation (1.24) are also sources for the streaming potential. Streaming potentials are therefore generated (i) by fluid flow perpendicular to a boundary where the coupling coefficient X changes, (ii) by fluid flow perpendicular to a boundary where the hydraulic mobility K changes, and (iii) wherever there are sources or sinks of fluid. Implications for Streaming Potentials Induced by Seismic Waves The above analysis (and indeed all of sections 1.2 and 1.3) has focused on the case where streaming potentials are produced by steady fluid flows. If we assume that equation (1.24) also holds approximately for flows which are slowly time varying, then a fourth source of streaming potentials, in addition to the three given above, can be expected for the case of fluid flow induced by seismic waves. Within the compressional wave, there is relative flow between the solid and fluid and hence V- J / * 0. Thus during the passage of the wave, the third term in equation (1.24) is non-zero and streaming potentials may arise even in homogeneous media. The topic of seismically-induced streaming potentials is discussed in detail in Chapter 5. ^ V - J ^ V X - V p - ^ V ^ j ' - ^ V - j ' . (1.24) Chapter 2 Field Methods 2.1 INTRODUCTION Unlike established geophysical techniques such as seismic refraction and reflection, there are no commercially available instruments, and no standardized field methods for making seismoelectric measurements. These have been developed incrementally by members of the U B C seismoelectric group since the mid-1980's, and the field experiments of this study have contributed significantly to their testing and improvement. Much of the credit for instrumentation design goes to Anton Kepic, and his PhD thesis (Kepic, 1995) is another useful source of information on that topic. The acquisition of seismoelectric data requires carefully designed instrumentation and field methods in order to obtain sufficient signal-to-noise ratios, and minimize acquisition artifacts. The latter point is of particular importance when making novel measurements such as those in this thesis as one must take special care to ensure that experimental procedures do not introduce unforseen transients that could be mistaken for, or obscure, the true signals. This chapter describes the instrumentation I used for making field measurements, the types of noise that affected the data, and field methods used to reduce noise and identify seismoelectric signals. 2.2 INSTRUMENTATION 2.2.1 Overview The main components of the data acquisition system were a seismic source, geophones, grounded dipole sensors, amplifiers, and an eight channel digitizer. A l l instruments were battery powered in order to eliminate electrical noise produced by portable generators. Unless otherwise stated, seismic arrivals were detected at the earth's surface using 26 Chapter 2: Field Methods 27 conventional vertical component geophones (like those used for oil exploration surveys) equipped with a spike that penetrated a few inches into the ground to improve mechanical coupling. They were manufactured by Geo Space Corporation (model GSC 20D) and were critically damped with a resonant frequency of 10 Hz, and a sensitivity of 31.4 V/(m/s). Geophones measuring the horizontal component of ground motion, or having a higher resonant frequency were used on occasion as discussed in Chapter 4. Hydrophones were used to measure seismic arrivals in water-filled boreholes. Electrical responses were measured by grounded dipole receivers consisting of two stainless steel stakes penetrating about 0.3 m into the ground and typically separated by 2 to 10 m. Single-wire leads ran from each stake to a local differential preamplifier having a high input impedance (2 MQ. differential), a gain of 30, and a bandwidth of 2 Hz to 30 kHz. The preamplifier's gain, and its differential output were used to boost the dipole signal and make it more immune to common-mode noise which was inevitably picked up during transmission to the recording site. More details on the grounded dipole and its preamplifier are given in section 2.2.3. Signals from each preamp and geophone were transmitted to the recording site along separate shielded twisted pair cables 20 to 100 m long. At the recording site, typically located in the back of a truck, a bank of passive, 60 Hz notch filters was available for optional notch filtering. Additional gain and bandpass filtering were then applied by passing each channel through a differential amplifier (Tektronix model AM502 modified for battery operation). Finally, the data were digitized and recorded by a RC Electronics 12 bit A /D board mounted in a portable computer. In order to improve signal-to-noise, records from repeated shots at the same shotpoint were sometimes added together (stacked) in real-time using software, written by R.D. Russell to control the A/D board. A l l instruments (including the dipole preamplifiers) were grounded to a stainless steel stake pounded into the earth at the recording site. Chapter 2: Field Methods 28 Each seismoelectric record contained up to 8 channels of data. Typical sample intervals and record lengths were 0.05 and 400 ms respectively. The bandpass filters were normally set with comer frequencies of 0.1 to 10 Hz at the low end and 1000 Hz at the high end. However, the filter slopes (12 dB/octave and 6 dB/octave in the lower and upper stopbands respectively) were quite gentle. As a result, the analog filters were only moderately effective at screening out high frequency noise, and a fine sample interval (0.05 ms) was required to ensure that aliasing was not a problem. Although the records were broader band than necessary, this made it easier to distinguish noise spikes from the lower frequency seismoelectric signals. Following data acquisition, excessive high frequency noise was routinely removed by digital lowpass filtering, and the number of samples per trace was usually decimated by a factor of 2 to 4 to make the data more manageable for further processing and plotting. 2.2.2 Seismic Sources and Triggering The main seismic sources used in this study were a sledgehammer, blasting caps, and a in-hole shotgun or 'buffalo gun'. The sources were selected for their portability, ease-of-use, and efficiency in the generation of seismic compressional waves (F-waves) in differing ground conditions. Stronger, explosive sources would have been useful for experiments on the Fraser Delta where seismoelectric effects were very weak, but this was ruled out due to logistical difficulties. The decision to focus on producing compressional rather than shear waves was made largely because the former are easier to generate, measure, and interpret in field data. Furthermore, as discussed in Chapter 5, we have good reason to expect that electrokinetic effects are induced by compressional waves, while the connection between electrokinetic effects and shear waves is more tenuous. The sledgehammer is a common seismic source for shallow seismic surveys because of its simplicity, portability, and low cost. In addition, responses from multiple hammer blows are Chapter 2: Field Methods 29 easily stacked or averaged in rapid succession to improve signal-to-noise ratios. A sledgehammer source was used extensively for seismoelectric experiments at the Haney test site. Early experiments revealed that the impact of a hammer on an aluminum base plate, such as that commonly used for hammer seismic surveys, generated an undesirable electromagnetic transient. This interference was removed by using a non-metallic base plate made of plastic or hardwood. It was not necessary to replace the steel hammer. On one occasion at Haney, we had the opportunity to compare the performance of the sledgehammer source to that of a lead ball weight drop. We determined that the seismic energy delivered by one swing of the 16 lb. (7.3 kg) sledge was roughly equivalent to that resulting from one drop of the 50 lb. (23 kg) weight from a height of 1 m. The potential energy of the weight before its release (220 Joules) provides an upper limit for the energy delivered by the sledgehammer blow. The buffalo gun source, described by Pullan and MacAulay (1987), generates seismic waves by detonation of a shotgun shell in a shallow borehole up to about 1 m deep. The design used in this study consisted of a steel barrel 90 cm long and 4.9 cm in diameter with a T-shaped handle at the top, chamber for the shell at the bottom, and a heavy firing rod that dropped onto the shell to initiate detonation. Three different chambers were available to accommodate 12, 10, or 8 gauge shells. Blank shells were preferred for reasons of safety, and the entire assembly was intentionally heavy (16 kg) in order to reduce recoil upon firing. The buffalo gun source was well-suited for the experiments I carried out on the Fraser Delta where the surficial soils were generally fine grained and moist, and shot holes were easily produced with a hand-auger. Under those conditions, the seismic pulse generated by a buffalo gun tends to be significandy stronger and broader band than that produced by a sledgehammer blow at the surface. The gun was not so useful at the Haney site where it was difficult to auger holes because of gravel in the soil. We have also used blasting caps in water-filled boreholes to produce seismic sources at depth. Conventional zero-delay caps, or seismic caps are not suitable because the large Chapter 2: Field Methods 30 initiation current can contaminate seismoelectric records at the time of detonation. Instead, we used fuse caps or electrical delay caps. Small (6 gram) explosive boosters, known as stingers, were occasionally slipped over the caps to provide more energy. Borehole damage was an unfortunate consequence of the blasting cap source. The boreholes at Haney were lined with P V C or ABS pipe having inner diameters of 2 to 3 inches. Detonation pressures destroyed the borehole lining around the blasting cap and blocked the borehole at that depth. We attempted to eliminate this problem by enclosing the blasting cap in an open-ended steel cylinder having a diameter slightly less than that of the plastic liner and a length of about 20 cm. That technique was abandoned when we discovered that it produced an electromagnetic transient at the moment of detonation. The same problem would be expected for detonations in a steel-cased borehole. Apparently, one should avoid deforming metal - either by hammer or explosive -when making electrical measurements close to the shotpoint. The RC Electronics digitizer employs a circular data buffer and a flexible trigger algorithm that allows the trigger point to fall anywhere within the record. That is, a record can include both pre-trigger and post-trigger data. This feature was useful for purposes of recording a sample of the ambient noise on each channel prior to the instant of source impact or detonation. Triggering was usually accomplished by use of an accelerometer mounted on the hammer or buffalo gun, or occasionally by a geophone in the ground close to the shotpoint. The accelerometers, made from piezoceramic audio transducers set in epoxy, had resonant frequencies in the 1-5 kHz range and acted as shock detectors.. We sometimes observed minor accelerometer 'crosstalk' on the dipole channels - apparently due to radiation of the accelerometer signal from its cable. This was eliminated by feeding the accelerometer signal directly into a miniature audio transformer (potted with the transducer) in order to produce a differential signal for transmission to the recording site. When blasting caps were used, data acquisition was triggered by the flash of light that accompanies cap detonation; this flash was transmitted to the recording site by a fiber optic cable, one end of which was taped to the Chapter 2: Field Methods 31 blasting cap (Kepic and Russell, 1996). Experience has shown that this method of triggering is reliable and accurate. Furthermore, it does not generate the electromagnetic noise associated with other triggering methods which rely on the termination or production of an electrical current at the time of detonation. 2.2.3 The Grounded Dipole Sensor A l l of the electrical data in this thesis were obtained using grounded dipole sensors to measure the difference in electric potential between two points in the earth. The sensor consisted of two stainless steel electrodes connected by insulated wires to a custom, high input impedance (2 M Q ) differential preamplifier located about midway between the electrodes. The preamplifier's input impedance must be significantly higher than the resistance between the electrodes in order to obtain an accurate measurement. For dipoles 1-10 m long, we have found that typical resistances lie in the range of a few kQ to several tens of kQ depending on the soil type and its moisture content. The first stage of the preamplifier is a passive bandpass filter used to reject D C voltages as well as very low and very high frequencies that lie well outside the range of interest (~ 10 Hz to 1000 Hz). The need for rejection of very high frequencies became apparent during early seismoelectric experiments (e.g. Butler, 1991; Russell et a l , 1992) when we observed that A M radio broadcasts had added significant noise to our records. Although the carrier frequencies (hundreds of kHz) were well outside the range of interest, nonlinearity in the amplifiers rectified and demodulated the radio signals with the result that we recorded their audio-frequency envelope (music and voices). The problem was discovered by listening to the analog output of the preamplifiers with headphones. It was euminated by constructing new preampli-fiers with operational amplifier chips having better high frequency (slew rate) performance, and by limiting the frequencies that could reach the op-amps through passive filtering. A review of the A M demodulation problem and its solution is given by Kepic (1995, Appendix B.8). 32 T-Box Preamplifier BANANA JACK INPUT P R2 499K VWvV R3 3K32 f—f—VWvV-*— 1 RF 47K5 BANANA JACK INPUT N CF 1 ' 100P CF 11 100P XLR CONNECTOR RF 47K5 • * V W W - + - 1 J 7 R2 499K -VvVVV 100N C1 SIGNAL _ V+ E1 - GND - E2 - V-PIN 2 1 3 5 4 10P C2 The choice of the 499K resistor is based on requirement that the lo-pass formed with the 10 pF capacitor have a corner at approx 30 kHz and that this comer change little with various source impedences. If the noise is a problem then replace the resi stor and capacitor with 47Kand 100pF respectively. Nominal Gain is 30 HP Corner for Dipole configuration is 2 Hz and for Long Wire Antenna is 160 Hz LP Corner for both configurations is 33 kHz with 12 dB/Oct rolloff until 1 Mhz Power Supply rails are decoupled by a 100 nH inductor and a 10 uF tantalum capacitor Geophysical Instrumentation Group Title : The Best Preamplifier Drawn By : AK Jan 1992 Ver: 1.3 Figure 2.1 T-box preamplifier circuit schematic (adapted from Kepic, 1995). Chapter 2: Field Methods 3 3 The 'T-box' preamplifier, designed by A . Kepic (1995) with input from other members of the U B C seismoelectric group and Dr. Yves Lamontagne, was used to acquire the vast majority of the data for this thesis. Its differential input-differential output circuit is shown schematically in Figure 2.1. The design is based on the three op-amp instrumentation amplifier which provides high input impedance, and good common mode rejection (Horowitz and Hi l l , 1980, p. 282). (In this case, the Tektronix amplifier at the recording site served as the third op-amp.) The frequency response for the grounded dipole and its preamplifier is determined by splitting the sensor into two sections: (i) the source impedance in combination with the preamplifier's passive filter, and (ii) the op-amp stage of the preamplifier. Transfer functions for the two sections can be determined independently because the op-amp stage has a very high input impedance and does not load the filter network. Figure 2.2 shows an equivalent circuit for section (i). The voltage source VVrepresents the open circuit voltage across the dipole, and Zs represents the source impedance. For a grounded dipole, at the frequencies of interest in this study (up to about 1 kHz), the source impedance is essentially just the resistance between the electrodes R5. However, a parallel source capacitance Cs is allowed for the sake of generality as it does not complicate the analysis. The output of the filter network V is the input to the op-amp stage of the preamp. The transfer function for the first section, obtained by standard techniques of circuit analysis is V jcoR.C, Vs 1+ +Z,)Ci +(/?, +R2)C2]-co2[Rl{R2 +Z,)+R2Z,]C1C2 Now, since Rx » Zs (recall that typical source impedances are of order 10 kQ), and since Cj » C 2 and RX>R2, equation (2.1) is well approximated by V jcoRxCx (2.1) Vs l + ja)RlC1-a>2R1{R2+Zt)C1C2 For high frequencies, this simplifies to a single pole lowpass filter, while at low frequencies it (2.2) Chapter 2: Field Methods 34 becomes a single pole highpass filter. The corner frequencies are given by 1 . 1 ®hP=— and « , „ = -RXCX (2.3) * (R2+Z,)C2 Note that in the T-box design, R2 was intentionally made much larger than the expected values of Zs so that the lowpass comer is approximately given by I/R2C2 and is insensitive to variability in the source impedance. Substituting Rj=lMCl, R2=A99 kQ, C/=100 hF, C2=10 pF, and Zs ~ 10kQ we find that the comer frequencies for the bandpass filter formed by the first stage of the T-box preamp are 1.6 Hz, and 31 kHz. Zs C1 -AWv-R2 2R1 Vs (%) C1 R2 C2 2 V Figure 2.2. Equivalent circuit for a grounded dipole connected to the input (passive filter) stage of the T-box preamp. The transfer function for the op-amp stage of the preamp can be derived in a straightforward fashion using the common, ideal op-amp approximation. Calling the input voltage V, and output V \" , the transfer function is V\" 2Z,+R, 2R, — = —!- = L + 1 . (2.4) V R3 R3(l + juRfCf) If the feedback impedance Z/ was purely resistive, then the transfer function would simply be a constant real-valued gain factor, determined by the ratio However, the feedback capacitance introduces frequency dependence. Examination of equation (2.4) reveals that it has two limiting forms depending on whether co exceeds or is less than a critical frequency Chapter 2: Field Methods 35 coc = (2Rf +R3yRfCfR3 (i.e. 992 kHz). For co »CO,, the second term in equation (2.4) dominates and the transfer function is simply V\"/V = 1 (unity gain). At frequencies well below 992 kHz the first term dominates and the response is that of a single pole lowpass filter with a comer frequency of co = yRfCf (33.5 kHz), and a passband gain of (1 + 2Rf/R3) or 29.6. The total frequency response of the T-box preamp is the product of transfer functions (2.1) and (2.4) and is illustrated in Figure 2.3, along with data points from calibration tests. Nominally, the preamp acts as a bandpass filter with a gain of about 30 in the 2 Hz to 31 kHz passband. The frequencies in Figure 2.3 extend to very high values (10 MHz) only for purposes of showing the change in preamp response near 1 MHz. However, at frequencies beyond about 100 kHz, the response may begin to deviate from that shown as the source capacitance Cs, and stray capacitance between the dipole wires become more important (Kepic, 1995). (Those capacitances were assumed negligible in producing Figure 2.3.) This is not a major concern since the frequencies of interest for this study are below 1 kHz. The intrinsic noise of the T-box preamp is dominated by the thermal (Johnson) noise associated with the two resistors labelled R2 in the passive filter (this can be shown by standard noise analysis techniques for circuits). Since the thermal noise of any resistive element is proportional to the square root of its resistance, a modified preamp design, known as the V -box, with a noise level about four times lower than that of the T-box, was produced by substituting a smaller resistor for R2. The only other difference between the two designs was the value used for C2. Equations 2.1 through 2.5 given above are therefore equally valid for the V-box preamp provided the new values of R2 =10 kQ and C2 =680 pF are used. Actual noise measurements for the two preamps are shown in Figure 2.4. The measurements were made by putting a 10 kQ resistor across the preamp inputs (to represent a typical source impedance) and measuring the rms output voltage from the preamp with a spectrum analyzer. After accounting for the passband gain of 29.6, the noise levels referred to the preamp input are approximately 130 and 33 nV/VHz for the T-box and V-box respectively in the 36 100 10 0.01 0.001 < U tea? >urc 3-merits alevatec Dy nois« i j ; R rt R/rvtavo \\ 3 • 12 dB/octave t--— • Calibation 1 (A.W. Kepic, 1991) O Calibration 2 (K.E. Butler, 1996) b d ctav e 11 nun I M i l ! ! I l l l l l l I III Hi 1.0E-2 1.0E-1 1.0E+0 1.0E+1 1.0E+2 1.0E+3 1.0E+4 1.0E+5 1.0E+6 1.0E+7 Frequency (Hz) 1.0E-2 1.0E-1 1.0E+0 1.0E+1 1.0E+2 1.0E+3 1.0E+4 1.0E+5 1.0E+6 1.0E+7 Frequency (Hz) Figure 2.3. Calculated amplitude and phase response for the T-box preamplifier. The response changes only marginally with typical variations in the source impedance Zs • Here it is assumed that Zs =10 k Q . Measured values of the amplitude response are also shown for two preamp calibrations. Chapter 2: Field Methods 37 10 Hz - 1000 Hz frequency range. The total integrated noise over a 500 Hz bandwidth in that range is therefore approximately 2.9 \\xV rms for the T-box and 0.74 u,V rms for the V-box design. These noise levels are comparable to the smallest seismoelectric signals I was able to detect during field experiments. The sacrifice made for lower noise was less stability in the preamp bandwidth; since R2 is of the same order as typical source impedances, the passive filter lowpass comer given by equation (2.3) can vary substantially with variations in Zs. For example with Zs = lOkQ, the lowpass corner is 12 kHz while for very low and very high source impedances of 1 kQ, and 100 kQ, the comer frequencies are 21 kHz and 2.1 kHz respectively. Given that the main frequencies of interest for this study are below 1 kHz, this variability was not a major concern. The V-box preamp was used for some experiments on the Fraser Delta where signals were very small and the lower noise design was particularly desirable. N •c cr > CD O 20 18 16 14 12 10 - T-box preamp -V-box preamp 0.1 10 100 Frequency (Hz) 1000 10000 Figure 2.4. Measurements of the intrinsic output noise of the T-box and V-box amplifiers in the range between approximately 0.1 and 1000 Hz. The preamp inputs were loaded by a 10 kQ resisitor to simulate a typical source impedance. Chapter 2: Field Methods ' 38 2.3 SOURCES OF NOISE There are many types of electrical noise that can obscure or completely overwhelm seismoelectric signals. These include ambient noise as well as acquisition artifacts which arise because of the measurement process. The latter can often be eliminated or controlled by careful design of instmmentation and field methods. The main sources of noise encountered during field experiments are listed and discussed below: i) Ambient noise - powerline harmonics - V L F and A M radio transmissions - natural atmospheric electricity (sferics) ii) Acquisition Artifacts - leakage, radiation or crosstalk from geophones or other sensors - E M fields associated with the impact or detonation of a seismic source - sensor shaking or electrode contact phenomena During most seismoelectric field experiments conducted by the U B C group, powerline harmonics have been the noise of highest amplitude. It is not uncommon to find that the ambient electric field in the earth at 60 Hz is 5-10 mV/m peak-to-peak. Such amplitudes are 5 to 1000 times larger than the signals I sought to measure. Harmonics of 60 Hz are weaker than the fundamental but several have amplitudes that are typically larger than seismoelectric signals. Some mines in Ontario use 25 Hz power in addition to 60 Hz with the result that both fundamental frequencies and their harmonics can contaminate electrical recordings. In Europe, electric railways can also be a troublesome source of harmonic noise (Strack et al. 1989). I used both analog instmmentation and digital data processing techniques to remove harmonic noise from field data. The former were sometimes required to attenuate noise prior to recording where seismoelectric signals were simply too small to be resolved in the presence of the much larger noise field given the limited dynamic range of the digitizer. The analog Chapter 2: Field Methods 39 methods included a passive 60 Hz notch filter, and noise cancellation by use of dipole arrays as described in section 2.4. However, these techniques alone were not adequate to reduce powerline noise to an acceptable level. Furthermore, notch filters are undesirable because of the amplitude and phase distortion they apply to the signal. As a result, I routinely relied on the data processing techniques described in Chapter 3 to complement or replace analog methods of powerline noise suppression. High-powered V L F (very low frequency) E M transmitters, used for naval communications with aircraft, ships, and submarines, are scattered all over the globe and have carrier frequencies in the range 15-30 kHz. During field experiments carried out in the Vancouver area, the 24.8 kHz signal transmitted from Jim Creek, Washington was a source of interference. Although data from the dipoles were normally lowpass filtered at 1 kHz prior to recording, the Tektronix filter slope (6 dB/octave) was not steep enough to eUminate the V L F signal. Because data were typically recorded with a sample interval of 0.05 ms (10 kHz Nyquist), the V L F interference was aliased during digitization and appeared at 4.8 kHz in the seismoelectric records. This was easily removed by digital lowpass filtering, but we sometimes chose to avoid the noise altogether by scheduling field experiments during the transmitter's weekly maintenance period (Thursdays). It is prudent to be aware of the potential for aliasing when selecting a sample interval. If, for example, a slightly finer sample interval of 0.04 ms had been chosen, the V L F signal would have been aliased down to 200 Hz - right in the frequency band of interest. The susceptibility of some amplifiers to audio-frequency noise produced by demodulation of A M radio signals has been discussed above (section 2.2.3), and has also been noted by Petiau and Dupuis (1980). We overcame this problem by use of a passive filter in the preamplifier. Kaufman and Keller (1981, p. 506) state that some demodulation of RF signals can also occur at the electrodes of a grounded dipole (i.e. at the metal-electrolyte interface) but that effect, if present, was not significant for our records. Chapter 2: Field Methods 40 Natural E M fields in the frequency range 10- 10000 Hz are caused mainly by atmospheric electricity (sferics) and generally occur as pulses associated with lightning strokes all over the world. At frequencies below about 10 Hz, natural fields are generated by interaction between the solar wind (charged particles and fields produced by the sun) and the earth's magnetosphere (Sumner, 1976, p. 135-144). Telluric (earth) currents induced by natural fields constituted a minor source of noise in records I acquired with grounded dipoles. The noise was rarely visible above cultural interference in the data but was sometimes apparent after the removal of powerline harmonics. Because it originated at great distances, I found that it could be at least partially cancelled out by the subtraction of recordings taken by a remote dipole (see section 2.4). Electrical noise caused by coupling (due to leakage, radiation or crosstalk) of geophone signals onto dipole channels is of particular concern because it could be mistaken for a seismoelectric signal. Crosstalk is typically characterized by its identical form on all affected channels. It is caused by induction or by currents leaking from one channel to another within multi-channel instrumentation such as amplifiers, filters, and digitizers. Electromagnetic radiation from geophones and the leakage of geophone currents into the earth differ from crosstalk in that they are received through the dipole sensors. Seismic signals may also radiate from poorly shielded cables. The standard test used by the U B C seismoelectric group to check for radiation, leakage, and crosstalk involves tapping the geophone (or other seismic sensor) after it and all other sensors have been laid out for an experiment. If any signals similar in character to that of the geophone channel appear on the electrical channels, geophone leakage, radiation, or crosstalk is suspected. Additionally, seismoelectric records were often taken twice - once with geophones in place, and a second time with all geophones removed, and the two records were compared to ensure that they were the same. These precautions enabled us to identify some problems with leakage/radiation from geophones and their cables during early field experiments. These problems were corrected by modifications to connectors and cabling, Chapter 2: Field Methods 41 and there is no sign of geophone interference in any of the seismoelectric data presented in this thesis. Radiation from an accelerometer used to trigger data acquisition is evident in some recordings; it appears at time zero and is can be distinguished from seismoelectric signals because it is dominated by frequencies above 1 kHz. As mentioned above (in section 2.2.2), this radiation appeared to come from the cable connecting the accelerometer to the recording site, and was eventually eliminated by use of differential signal transmission. The impact or detonation of a seismic source can produce acquisition artifacts. As mentioned previously, sledgehammer blows and blasting cap detonations generated E M transients when they applied a mechanical shock to metal. This interference was eliminated by removing the metal strike plate or blast chamber from the shotpoint. Other blast-associated E M artifacts have been reported by a number of investigators (e.g. Kepic, 1995; Kepic et al., 1995; Kondrashev, 1980, Sobolev andDemin, 1980; Neishstadt et al., 1972). These were not observed during my experiments - possibly because I used only blasting caps and very small (6 gram) charges - but they are worthy of mention. Kondrashev (1980, p. 52-59) suggests that blast E M may be produced by several mechanisms including the expansion of thermally ionized explosion products, as well as fracturing and piezoelectric effects in the rock irnmediately surrounding the charge. Experiments by Kepic and Russell (1996) indicated that blast E M problems can be exacerbated by any wires (such as electrical detonator leads) leading into the explosive as such wires can re-radiate the blast E M transient. In the experience of the U B C seismoelectric group, blast E M artifacts can be minimized or eliminated by detonating explosives with a thermal fuse cap, using the fiber optic technique discussed in section 2.2.2 to obtain a time break. For more details on blast-associated E M , the reader is referred to Butler (1991, p. 25-26) for a summary of Russian reports and to Kepic (1995). Shaking of a grounded dipole sensor due to the passage of a seismic wave is another potential source of noise. When dipole wires or the electrodes themselves are vibrated, the magnetic flux (from the earth's field) linking the effective loop formed by the dipole can vary Chapter 2: Field Methods 42 and produce an emf by induction. Petiau and Dupuis (1980) stated that noise amplitudes of about 10 LtV were typical for the oscillation of dipole wires in the wind. However, their dipole lengths (about 100 m) were much longer than those used in this study, and movements due to wind tend to be much larger than those caused by a seismic wave. It therefore appears unlikely that seismic shaking of a grounded dipole could produce significant noise by magnetic induction. The possibility that mechanical vibrations could disturb electrochemical potentials at the electrode -soil interface ('contact phenomena') has also been considered by some authors (e.g. Thompson, 1939; Ivanov, 1949) but mainly for the case where a strong current is being forced to flow through the electrodes. My own experiments at the Haney, B C indicate that sensor shaking is not a source of concern at that site (see discussion of Figure 4.6 (seismoelectric response vs offset)). Nonetheless, the sensitivity of the grounded dipole to seismic shaking is not well known, and it is important to keep this in mind when designing and interpreting the results of seismoelectric experiments. 2.4 FIELD METHODS Although there are no standardized procedures for the acquisition of seismoelectric data, there are some principles which guide the choice of source-receiver geometry: (i) Dipoles should be laid out over a wide enough range of offsets so that a simultaneous seismoelectric arrival can be clearly distinguished from the seismic arrivals at each geophone. (ii) It is useful to place a geophone at each dipole for purposes of correlating seismic and electrical arrivals, and identifying possible artifacts associated with dipole shaking. (iii) The mapping of an interface by its seismoelectric response requires data from a series of shotpoints. (iv) A remote reference dipole, or dipole arrays (discussed below) can be useful for suppression of regionally coherent noise such as powerline harmonics and sferics. (v) Conventional seismic refraction or reflection data should also be acquired to provide knowledge of geologic structure and seismic velocities. Chapter 2: Field Methods 43 Of course, the choice of source-receiver geometry also depends on the number of recording channels available. With my eight-channel system, I commonly re-occupied the same shotpoint with different arrays of receivers in order to satisfy all of the principles listed above. Two examples of useful eight-channel layouts are given in Figure 2.5. Layout (a) has dipoles at 6 different distances from the shotpoint, and a geophone at the electrode closest to the shot to determine when seismic shaking might first begin to affect the dipole sensors. The eighth channel is connected to an accelerometer on the seismic source and provides the trigger for data acquisition. A remote reference dipole is located far from the shotpoint beyond the expected range of the seismoelectric signal (typically 60 to 100 m away during my field experiments). It is used to monitor noise from distant sources such as power lines and sferics which tends to be relatively uniform in character over the survey area (although the amplitude may vary depending on the ground resistivity around each dipole). Such regional noise can be partially cancelled by subtracting scaled versions of the remote dipole record from recordings taken simultaneously by the other dipoles. During some experiments, I used the differential amplifiers at the recording site to carry out the subtraction in real time. This process, known as analog mixing, analog balancing, or analog noise cancellation, reduces the noise level prior to digitization and is therefore desirable for real-time quality control and for making the best use of the digitizer's limited dynamic range. However, I had only eight differential amplifiers at the recording site and every one used for analog balancing reduced the number of channels available for independent measurements. As a result, the remote dipole data were often subtracted from other channels digitally after data acquisition as discussed in section 3.3. Chapter 2: Field Methods 44 a) accelerometer remote dipole r©i r®i r®i r®i dipoles r ® i -4tr SP geophone b) accelerometer r©i r©i r®i KS>i r©i SP geophone Figure 2.5. Two useful source-receiver geometries for the acquisition of seismoelectric data with an eight channel recording system. Layout (b) in Figure 2.5 differs from (a) in the arrangement of the six dipoles. They are located on either side of the shotpoint at three different offsets and there is no remote sensor. One advantage of this arrangement is that it allows (partial) investigation of the radial symmetry of the seismoelectric response. As discussed in section 4.2.3, the symmetry of signal polarity is an important characteristic for distinguishing between seismoelectric effects of electrokinetic origin, and those caused simply by the modulation of earth resistivity in the presence of a telluric current. Furthermore, the symmetry can be exploited for purposes of regional noise cancellation; taking the difference between two dipole placed symmetrically about the shotpoint suppresses the regional noise common to the dipoles while enhancing the desired seismoelectric response. Thompson and Gist (1993) used this technique during large scale field experiments and I have also found it useful on occasion. However, because the method yields an average of the signals measured by two dipoles, it is not ideal for studying how the seismoelectric response varies with dipole position. Chapter 2: Field Methods 45 Some additional points should be made regarding the layouts in Figure 2.5. Firsdy, the number of dipole channels has been maximized at the expense of seismic information. Additional shots (or additional channels) with more geophones would be required to provide proper seismic control. Secondly, the optimum shotpoint-dipole offset depends on the depth to the target. Modelling presented in Chapter 5 indicates that the seismoelectric effect from a horizontal interface is strongest at an offset equal to one-half the interface depth. Finally, data from a series of shotpoints are required to map a boundary by its seismoelectric response. The dipoles may remain stationary as the shotpoint is moved away provided the seismoelectric signal remains strong enough to be reliably detected. However, dipoles would have to be moved with the shotpoint in order to maintain the symmetrical array of Figure 2.5(b). Chapter 3 Data Processing Methods 3.1 INTRODUCTION This chapter documents digital processing methods used to improve the signal-to-noise ratio of seismoelectric records. The techniques were chosen for their relative simplicity and for their ability to produce major improvements with minimal signal distortion. The ability to apply these techniques was an important advantage I had over previous investigators who were limited by the analog recording systems of their day. Two general types of processing were used in this study: (i) frequency filtering, and (ii) noise estimation and subtraction. I have grouped the techniques according to the type of noise they were designed to eliminate. Section 3.2 describes bandpass filtering used to suppress all types of noise at frequencies outside the signal band. In section 3.3,1 present an algorithm for automatically scaling and subtracting regional noise recorded by a remote reference. Section 3.4 describes three techniques used exclusively for the suppression of powerline harmonics -notch filtering, sinusoid subtraction and block subtraction. The latter two techniques have proved to be particularly effective and constitute one of the significant contributions of this thesis. They are also described in a paper by Butler and Russell (1993). Table 3.1: Main programs used to process seismoelectric records. Program Function mlpfdec.for frequency domain lowpass filter bpf.for frequency domain bandpass filter iirfilt.for applies an UR time domain filter (e.g. Butterworth) lstrsub.for remote reference subtraction (with least-squares fit) multq3.for remote reference subtraction (gain specified by user) mbsub.for block subtraction ssub.for sinusoid subtraction dnotch9.for frequency domain notch filter 46 Chapter 3: Data Processing Methods 47 Table 3.1 list the names and functions of the main Fortran programs I wrote to process seismoelectric data. More details are given in Appendix A. The programs were compiled for the MS-DOS operating system so that they could be applied on the same computer used to acquire data in the field. Some in-field processing was required for quality control because seismoelectric signals can rarely be identified in the raw data. Other processing techniques, besides those described in this chapter, have been used by myself and others in the U B C seismoelectric group for processing particular data sets. Some examples include (i) cross-correlation of seismic and seismoelectric traces for the identification of signals in noisy data (Russell et al., 1992), (ii) application of the Karhunen-Loeve transform (also known as principal component analysis) for the enhancement of signals that arrive simultaneously across many traces in a multi-channel record, and (iii) calculation of spectrograms to reveal how the frequency content of seismoelectric responses vary with time (Kepic, 1995). These techniques are not discussed further here because they were not required to enhance or interpret the data presented in this thesis. 3.2 LOWPASS AND HlGHPASS FILTERING The analog filters used for data acquisition commonly admitted a much wider range of frequencies than were necessary to record the signal alone. This was due in part to the limited number of cut-off frequencies available on the Tektronix AM502 amplifiers, and the gentle stopband slopes of those filters. But, we also found it beneficial to make broadband recordings in order to minimize signal distortion, and for purposes of differentiating between noise spikes, and narrower band signals. Digital lowpass and highpass filters, with cut-off frequencies and slopes tailored to the specific noise conditions in each data set could then be applied following acquisition. A digital lowpass filter was applied routinely as the first step in processing to attenuate high frequency noise such as V L F signals, high order powerline harmonics, and atmospheric Chapter 3: Data Processing Methods 48 electricity. This improved signal-to-noise, and allowed data to be decimated to a lower sample interval (e.g. from 0.05 ms to 0.2 ms) without fear of aliasing. Highpass filtering was applied only in special cases to attenuate low frequency noise (rarely a problem), or for purposes of spectral balancing (i.e. to reveal the higher frequency components in records dominated by low frequency signals). Both frequency domain, and time domain filters were used. The former were simpler to design for any desired cut-off frequency and slope, while the latter were advantageous for faster in-field processing. The design and characteristics of digital filters is a broad, well documented subject. My purpose here is to describe briefly how the filters were applied, and highlight some of their characteristics. Filtering a time series in the frequency domain involves multiplying its (discrete) Fourier transform by a transfer function and transforming the product back to the time domain. The forward and inverse Fourier transforms are computed most efficiently by FFT algorithms such as the subroutine R E A L F T from the Numerical Recipes library (Press et al., 1989). Filter transfer functions are generally complex functions of frequency commonly written in the form frequency domain filters used in this study are called zero-phase because they have 0(co) = 0. The phase delay ((J)(co)/co) and group delay are therefore zero for all frequencies. In a practical sense, this means that pulses (such as seismoelectric signals) are not delayed by filtering, and their shapes remain similar apart from the attenuation of the high (or low) frequencies outside the passband. The amplitude response chosen for the lowpass frequency domain filter was where coc is the cut-off frequency, and n is the filter order. The response of the highpass filter is obtained by substituting © c/co for co/co c in the above equation. Note that, regardless of the H(a) = \\H(®)\\eMa>) ( where |//(co)| is the amplitude response, and ())(co) is the phase response of the filter. The 3.1) (3.2) Chapter 3: Data Processing Methods 49 filter order, the amplitude at the cutoff frequency is 1/2 or 6 dB down. At frequencies in the stopband, well beyond the cut-off, the amplitude response falls off at a rate of 6n dB/octave. Digital filtering in the time domain was carried out using an DR (infinite impulse response) form of the classic Butterworth analog filter. Given an input time series xk, and filter coefficients [b0,,...bn; ax, a2,...an ], the filtered output series yk was obtained by convolution: yk = b0xk + bxxk_x + ... + bnxk_n -a,yk_i + ... + anyk_n. Equation (3.3) is really the sum of two convolutions - one between an operator bk and the input data series, and a second between operator ak and the output data series. The challenging part of time domain filtering is determining the required filter coefficients. The lowpass or highpass coefficients for a given order n and normalized cut-off frequency (co c/co N where co N is the Nyquist frequency) were obtained using the Butterworth filter design function 'butter' in the mathematical software package 'Matlab'. This function first determines the transfer function for the appropriate Butterworth analog filter prototype, and then applies the bilinear transformation to find the discrete (z-transform) equivalent and filter coefficients. Following this, a Fortran program was used to apply the filter according to equation (3.3). From equation (3.3) it is clear that the filtering of a time series containing M points would require approximately 4nxM operations (multiplications and additions). By comparison, frequency domain filtering is dominated by the time required for two FFTs each involving of order Mlog 2 M operations. In practice, for a raw seismoelectric trace containing 8192 points, 4th order Butterworth filtering was approximately four times faster than frequency domain filtering (requiring 7.5 versus 30 seconds/trace on our 80286-based field computer). This difference in speed can amount to substantial time savings in the field where several traces may require filtering after each shot. The Butterworth filter transfer function can be found in standard texts on filter design (e.g. Rorbaugh, 1993) and need not be given here, but the amplitude response for an n-pole Chapter 3: Data Processing Methods 50 (or order n) Butterworth lowpass is simply N » l = j 1 • (3-4) This response is very similar to that assigned to the frequency domain filter, except that its value is only 3 dB down (i.e. l/V2) at the cut-off frequency. The rate of decay is 6n dB/octave at frequencies well beyond the cut-off, and the highpass version is again obtained by substituting © c/co for co/co c . However the Butterworth filter differs markedly in that its phase response (])(co) is neither zero nor linear. The phase and group delays vary significandy with frequency especially near the cut-off. As a result filtered pulses will tend to appear distorted in shape, and their main energy will appear delayed in time particularly if the pulse contains frequencies in the vicinity of the cut-off. If necessary all phase shifts can be removed by re-filtering the time series in reverse, but in many cases, the signal distortion is small or unimportant. The Butterworth response may even be preferred over that of a zero phase filter because it is causal and matches the response of its well known analog filter prototype. Digital step responses for two zero-phase frequency domain and Butterworth time domain filters are presented for comparison in Figure 3.1. The lowpass designs are 6th order with a normalized cut-off frequency of 0.2, while the highpass are 2nd order with a normalized cut-off of 0.05. For a sample interval of 0.2 ms, these cut-offs would be 500 Hz and 125 Hz, and the window shown in Figure 3.1 would be 20 ms long. Note that the zero phase filters preserve the symmetric nature of the step; the step is spread out both backward and forward in time showing that this type of filter is acausal. The Butterworth responses are not symmetric and a net delay of the step is clearly evident in the lowpass case. These examples illustrate that the apparent arrival time of a pulse can be modified by frequency filtering. It is important to keep this in mind when picking arrival times or 'first breaks' of seismoelectric or seismic signals. As shown in Figure 3.1, the Butterworth filter is superior to the zero phase filter for highpass filtering if the preservation of first breaks is of utmost importance. 51 Input Step Lowpass Filter Step Responses f c / f N = 0 . 2 s l o p e = 3 6 d B / o c t Highpass Filter Step Responses f c / f N = 0 . 0 5 s l o p e = 1 2 d B / o c t 950 970 990 1010 Sample Number 1030 1050 Figure 3.1: Comparison of step responses for two types of digital lowpass and highpass filters. The zero phase filter preserves the symmetry of the step. The Butterworth design emulates the response of a common analog filter. Chapter 3: Data Processing Methods 52 3.3 R E M O T E REFERENCE SUBTRACTION As discussed in Chapter 2, E M noise from distant sources such as power lines and atmospheric electricity (spherics) is typically similar in amplitude and phase over a broad area. Thus, records acquired by a remote reference dipole, far from the shotpoint and beyond the range of the seismoelectric signal, can be taken as an estimate of the noise recorded simultaneously by other sensors closer to the shot. Subtraction of the remote reference from the other sensors often yields significant improvements in signal-to-noise. In my experience, improvements by factor of about 10 (20 dB) are typical. While this process can be done in the field by use of a differential amplifier (analog balancing), it is logistically easier to do it digitally after acquisition. This reduces demands for instrumentation in the field. It also allows a computer to be used to calculate an optimum scaling factor for the remote reference data prior to subtraction. Here, we define the 'optimum' factor as that gain a by which the reference trace h(t) should be multiplied so as to minimize, in a least-squares sense, the difference between it and a seismoelectric record r(r). That is, we seek to minimize over which r(r) contains noise but minimal signal. Setting dQ/da = 0 and solving for the optimum gain a yields The integrals in equation (3.6) are evaluated numerically and the processed record r'(t) is obtained as 2 (3.5) with respect to a. The interval \\tx ,t2], called the design window, should ideally be an interval (3.6) r'(t) = r{t)-ah(t). (3.7) 53 a) Raw Data M A / v W A W A V v W A \\ A M ^ Record Hanf17 Hammer Impact Dipolel Dipole2 1.0 V Dipole3 Remote Reference Dipole 80 160 240 TIME (ms) 320 400 b) Data after subtraction of the remote reference trace Dipolel Dipole2 0.2 V Dipole3 160 240 TIME (ms) 320 400 Figure 3.2: Illustration of noise removal by the subtraction of a remote reference. Frame (a) shows a raw seismoelectric record generated by stacking 20 sledgehammer blows at the Haney field site. The dashed line indicates the moment of hammer impact. Dipoles 1,2 and 3 were located 12, 14, and 16 m N of the shotpoint. The remote reference dipole was 60 m N of the shotpoint. Processed data, obtained by automatic scaling and subtraction of the remote reference is given in frame (b). A seismoelectric signal is now clearly visible on all three dipoles. Note that the amplitude scale has been expanded by a factor of about 5 relative to the upper plot. Chapter 3: Data Processing Methods 54 This technique was used to process many data sets, particularly when sinusoid subtraction alone was unable to remove all powerline noise. An example involving data from the Haney site is given in Figure 3.2. The raw data in the upper panel are time-varying potential differences recorded simultaneously by four dipoles, one of which was a remote reference located 60 m away from the shotpoint. The responses to 20 sledgehammer blows were stacked (averaged) to form this record, and the signal amplitudes (approximately 1 V peak to peak) include an amplifier gain of 1400. The traces are dominated by powerline noise at 60 Hz and its third harmonic (180 Hz). The character (or phase) of the noise is remarkably similar from trace to trace but its amplitude varies with dipole position by a factor of almost two. The results of subtracting the remote reference from each of the other traces, using the 130 ms period before hammer impact to calculate the optimum gain, are shown in the lower plot. The optimum gains were 0.58, 0.51, and 0.78 for dipoles 1, 2, and 3 respectively. The noise level has been reduced greatly and seismoelectric signals are now visible on all three traces. However, some residual noise remains because the harmonics on the remote reference were not perfectly in phase with those on the other traces. These residual harmonics can be suppressed using other techniques discussed below. 3.4 R E M O V A L OF POWERLINE HARMONIC NOISE Powerline harmonic interference was typically 1 to 3 orders of magnitude larger than the electrical signals I sought to measure in the field, and I routinely relied on data processing to complement or replace the analog methods of powerline noise suppression discussed in Chapter 2. Notch filter processing, discussed briefly in section 3.4(c), is a common and robust method for suppressing harmonic noise. However, notch filters also attenuate and distort the signal at frequencies in the vicinity of each notch. In order avoid this problem, I developed two novel techniques for estimating and subtracting harmonic noise in a time series (Butler and Russell, 1993). These techniques, called sinusoid subtraction and block subtraction are discussed below and proved to be the most useful tools for powerline noise removal. Sinusoid Chapter 3: Data Processing Methods 55 subtraction has been used to process almost all of the seismoelectric data presented in this thesis. Both of the harmonic noise subtraction techniques assume that the seismoelectric record r(t) can be represented by the sum of a signal s(t), nonharmonic noise e{t), and harmonic noise p(i) having a fundamental frequency fa: The amplitude, phase, and frequency of all powerline harmonics are assumed to remain constant over the length of the record. This approximation is generally reasonable for record lengths of a few seconds or less. The fundamental frequency of 60 Hz power transmission in Canada and the United States rarely deviates by more than 0.03 Hz from its nominal value and variations in frequency tend to occur slowly over a period of minutes to hours (Adams et al., 1982; British Columbia Hydro, personal communication). Powerline fields may also include weaker nonstationary components such as sidebands and subharmonics generated by motor loads, and transients generated by switching of current loads or rectifiers. Such nonstationary components would form part of the nonharmonic noise e(t). Block subtraction is a simple process that takes advantage of the fact that p(t) is periodic with period 1/f. An interval of r(t) over which s(t) and e(t) are small is taken for the powerline noise estimate. This interval or block of data is then shifted by an integral number of periods m/fo and is subtracted from r{t) to remove the powerline harmonics from a second interval. r{t) = s(t) + e(t) + p(t) (3.8) where p(t) = ^ ck cos(.2izkf0t+$k). (3.9) Sinusoid subtraction involves estimating the amplitude and phase of each of the harmonics contributing to p(t) and subtracting them from the record. Reductions of up to 45 dB in the ambient noise level have been obtained by using this method to subtract up to 25 Chapter 3: Data Processing Methods 56 harmonics of 60 Hz from seismoelectric data. A start along these lines was made by Butler (1991), but the derivation and analysis presented here is much more comprehensive. The method is similar to that proposed by Nyman and Gaiser (1983), but we have assumed that the fundamental frequency / is known. This simplification allows us to cast amplitude and phase estimation as a least-squares minimization problem and carry out a rigorous analysis of the sources of estimation error. In particular, we demonstrate that the error associated with the presence of multiple powerline harmonics can be eliminated by calculating estimates over an interval that is an integer number of cycles of the fundamental frequency. Finally, we discuss the performance of the method when the assumed value offo is inaccurate and suggest approaches to obtaining an improved frequency estimate. Block subtraction and sinusoid subtraction are most effective when p(t) is estimated over a portion of the record where the nonharmonic components (s(t) and e(t)) are small. Under those conditions, the methods are capable of suppressing harmonic interference with minimal distortion of the signal. In contrast, notch filters always attenuate and distort the signal at frequencies in the vicinity of the notch. 3.3.1 Block Subtraction Block subtraction is a simple method for suppressing powerline noise in a record containing an interval [ti, f2] over which the nonharmonic components are negligible. This interval becomes the powerline noise estimate or noise block. Since p(t) has a period of l/fo, the powerline harmonics in any other interval [tj + m/fo, + m/fj, where m is an integer, can be removed by shifting and subtracting the noise block from the record. It is important to note that any signal or nonharmonic noise present in the noise block will appear (inverted in polarity) in the processed interval following subtraction. As a result, the method is useful only if the nonharmonic components in the noise block are negligible compared to the signal in the interval to be processed. Chapter 3: Data Processing Methods 57 The block subtraction procedure is illustrated in Figure 3.3. Trace A shows the potential difference measured across a 5 m grounded dipole receiver during a seismoelectric test. A small explosive charge was detonated at 47 ms in a shallow borehole 16 m away from the dipole. The data were recorded using 0.1 - 1000 Hz band-pass and 60 Hz notch filters and the sample interval is 0.032 ms. Harmonics of 60 Hz clearly dominate the record. The largest harmonics are 540 Hz, 180 Hz, and 300 Hz although there are several others of significant amplitude. The amplitude of the noise is approximately 5 mV representing an average electric field of 1 mV/m across the dipole. Note that, because/ is 60 Hz, the noise pattern repeats itself every 1/60 s (16.7 ms). The first 47 ms of trace A were recorded immediately before the charge was detonated and therefore contain no signal s(t). The interval also contains little nonharmonic noise (an observation confirmed by applying the sinusoid subtraction technique to the region) and therefore represents a good sample of the powerline noise. Trace B in Figure 3.3 shows the noise block shifted to the right by two cycles of 60 Hz (i.e., delayed by 33.3 ms). Subtraction of B from A removes the powerline noise between 33.3 and 80.3 ms and yields trace C. A seismoelectric signal, totally obscured in the original data, is now clearly visible at 52 ms. This figure also illustrates the computational simplicity of block subtraction; all harmonics of 60 Hz were simultaneously removed by shifting and subtracting a sample of recorded noise. In practice, it is rarely practical to apply a time shift ts that is exactly m/fo seconds. The finite sample interval of the record and uncertainty over the exact value of the fundamental frequency will limit the accuracy of the shift. Although noise cancellation will be optimal when ts = m/f0 , good results will be obtained provided that \\ts -m/f0\\ is substantially less than one-half period of the highest frequency harmonic in the data. In some cases it may be worthwhile to resample the data in order to obtain a better match between the desired and available time shifts. An ideal sample rate would be an integer multiple of / . 58 0 20 40 60 80 100 120 TIME (ms) F i g u r e . 3 .3 : I l l u s t r a t i o n o f the b l o c k s u b t r a c t i o n m e t h o d . T r a c e A i s a r a w s e i s m o e l e c t r i c r e c o r d c o n t a m i n a t e d b y p o w e r l i n e h a r m o n i c s . A s i g n a l i n the c e n t r a l p o r t i o n o f the t r a c e has b e e n u n c o v e r e d b y s u b t r a c t i n g a s h i f t e d b l o c k o f h a r m o n i c n o i s e . Chapter 3: Data Processing Methods 59 A block subtraction-like process can also be applied in real-time during data acquisition using a digital delay line of the type used by musicians for echo effects. These instruments accept an analog input, digitize it, delay it by an adjustable period ranging from about 1 ms to 1 s, and convert it back to an analog output. The current and delayed signals can then be mixed using a differential amplifier to subtract one from the other. I experimented with two such delay lines and obtained reasonable harmonic cancellation by specifying delays of length m/fo , but found that the performance was limited by the fact that the frequency responses of the delay lines began to degrade below about 100 Hz (i.e. near the bottom of the frequency range important for music). This limitation could be overcome by a custom-designed delay line, and the real-time technique would be useful for expediting quality control in the field. However, the application of block subtraction during processing was sufficient for my research needs. 3.3.2 Sinusoid Subtraction A least-squares algorithm for amplitude and phase estimation Powerline harmonics occur at the distinct frequencies nfo where n is any positive integer. As a result, they can be suppressed during processing by subtracting from the record, sinusoids of the appropriate frequency, amplitude, and phase. If / is known then the amplitude and phase of each harmonic can be estimated using the least-squares approach described below. The problem is to estimate the amplitude and phase of the nth powerline harmonic present in a record r(t). The harmonic noise p(t) given by equation (3.9) may also be written in the form p(t) = ^ iakcosk(H0t+bksmk(aot (3.10) where coo = 2K f0 is the fundamental (radian) frequency of power transmission and the amplitude cn and phase § n of the nth harmonic are given by Chapter 3: Data Processing Methods 60 cn=Ja2n+bl (3.11) and ^ t a n - ' K / a J . (3.12) A. The objective is to find estimates aR and for the coefficients an and bn which can be used to calculate amplitude and phase. This may be done by choosing dn and bn such that a sinusoid of frequency nf0 fits the record (or a portion of the record) as closely as possible in a least-squares sense. That is, the integral 9„ = £ [ r ( r ) - a n c o s « c o o r - 6 „ s i n n c o o r ] dt is minimized with respect to an and bn. The constant x denotes the duration of the estimation window. Setting — = 0 and — T T - = 0 and solving for an and bn yields the result dan dbn an = p £ r ( r ) ( a 3 c o s « c o 0 r - a 2 sinn(n0t)dt (3.13) bn = p £ r ( f ) ( a ! sin nco 0 f - a 2 cos n(0ot)dt (3.14) where the a's and /? are constants given by f1 2 \" j 1 1 sin2rtCfl„T a, = cos ntatdt = —+ — J o 0 2 4rtco0 rx . , l-cos2«co T a 2 = cos«coorsin«coorrff = 2— J o 4«(o„ a 3 = J o X s i i . 2 , t sin2«co„x sin ncojdt = 2 4«coo and J3 = l / ( a 1 a 3 - a 2 ) . Equations (3.13) and (3.14) define the least-squares estimates of an and bn. Estimates for the amplitude and phase of the nth harmonic are obtained using equations (3.11) and (3.12). Note that if the length x of the estimation window is an integer multiple of the period of the fundamental frequency (x = mlf = 2iim/(x)o , m an integer) then cc2 = 0 and an and bn - m for x = — Chapter 3: Data Processing Methods 61 simplify to the trigonometric Fourier series coefficients for frequency nu>o: a„ = — J o r(r) cos rtcoordf (3-15) 4 =-F r(t) sin ncajdt. (3.16) x J ° When the estimation window covers the whole record, sinusoid subtraction is equivalent to a very narrow notch filter. This filter could be implemented in the frequency domain by taking the FFT of the record, setting the nfo component equal to zero, and taking the inverse FFT. However, the record length would have to be an integer multiple of 1/f in order to ensure that the discrete frequencies in the FFT representation included nf. Error analysis A, Some insight into the accuracy of estimates an and bn may be gained by using r(t) = s(t) + e(t) + ancosnaot + bn smn(00t + ^ a k cosk(0ot+bk sinktaj to substitute for r(t) in equations (3.13) and (3.14). Equation (3.13) becomes an = aN+Eas+eAE+EAP (3.17) where e ° = p j^(f)(a 3 costtC0 o r-a 2 sinn(d0t)dt (3.18) e ° = P Je (0 (a 3 cosrcco 0 f-a 2 s,inn(o0t)dt (3.19) and eap = ^ >^f (akcoskmj+bk sink&j)(a3cosn($j-a2 sinn(00t)dt. (3.20) t=i 0 Ideally, we would like to obtain the result an = an. The terms , 8 ° , and e° represent errors in the estimate an due to the presence of the signal, nonharmonic noise, and powerline harmonics other than harmonic n. Equation (3.14) may be expanded in the same fashion as equation (3.13) to obtain an expression for bn analogous to that for dn: k=bn+zbs+£be+£bp • (3-21) Chapter 3: Data Processing Methods 62 Expressions for £*, ebe, and e* are obtained by substituting (o^ sinnco of-a 2cosnco or) for ( a 3 c o s « c o o r - a 2 sinncoof) in equations (3.18), (3.19), and (3.20) respectively. When x = m/f0 (m being any positive integer), the errors e° , ebs, e\" and £* are equal to the Fourier series coefficients at frequency nfo for the signal and nonharmonic noise contained within the estimation window. These errors may be minimized by calculating estimates over a portion of the record where s(t) and e(t) are absent, small, or have little energy at frequency nf. The best estimation window could be a short interval over which the nonharmonic J o components are particularly small, or a longer interval over which they may be larger but less correlated with the harmonic of interest. The integral expressions for the errors e\"p and ebp associated with the presence of other harmonics in the record are evaluated in Appendix B. Equation (B-2) shows that these errors are functions of the window length x and are, in general, non-zero. However, for windows that cover an integer number of the periods of the fundamental frequency (i.e., for x = m/f0 where m is any positive integer), the numerators of all terms in equation (B-2) are exacdy zero and eap = ebp — 0. This is a consequence of the fact that two harmonic signals cos(fccoo£+(|)) and cos(no)0t+y) are orthogonal only over an integer number of periods of the fundamental. The effect of window length on amplitude and phase estimates can be illustrated by a simple example. Consider a record composed solely of two sinusoids that are the fifth and sixth harmonics of 60 Hz: r(t) = A[sin(27t 300r) - sin(27t 360r)], where A is an arbitrary amplitude factor. The 5th harmonic is to be estimated and subtracted from the record. Since s{t)=e(t)-Q, we can use equations (3.17), (3.21), and (B-2) to calculate a5 and bs analytically as a function of x. The errors in the amplitude and phase estimates calculated from a5 and bs are plotted in Figure 3.4(a) for window lengths ranging from 10 to 100 ms. Figure 3.4(b) shows the amplitude of the residual 300 Hz that would remain after sinusoid subtraction. Chapter 3: Data Processing Methods 63 Perfect estimates of the amplitude and phase are obtained when x = m/f0 = m/60 seconds; at these lengths the amplitude error, phase error, and residual 300 Hz amplitude are zero. It is also clear that the local error maxima decrease as the window length increases. The performance of the sinusoid subtraction technique on real data is illustrated in Figure 3.5. Trace A shows the potential difference recorded across a 10 m grounded dipole during a seismoelectric experiment. A sledgehammer was used to strike the ground at 64 ms at a distance of 5 m from the dipole. The recording bandwidth was 0.1-1000 Hz and the sample interval is 0.128 ms. The record is clearly dominated by 60 Hz noise having an amplitude of approximately 25 mV (2.5 mV/m) and there is littie evidence of any seismoelectric signal. The other traces in Figure 3.5 show the trace after successive sinusoid subtractions of 60 Hz (B), 180 Hz (C), 300 Hz (D), 360, and 420 Hz (E), and 480, 540, 660, 720, 780, and 240 Hz (F). Note that the raw data are plotted at a greatly reduced scale relative to the other traces. A seismoelectric arrival at 70 ms is clearly visible in trace F after the subtraction of eleven powerline harmonics. The noise level has been reduced by about 45 dB relative to trace A. The amplitude and phase of each of the eleven powerline harmonics were estimated over the first 50 ms of the record. The error terms Eap and ebp were small since the estimates were calculated over an integer number of cycles of 60 Hz. Since there is no signal and negligible nonharmonic noise over this interval, the errors eas,eb,eae, and ebe were also small. Excluding the signal from the estimation window ensured that it did not contaminate the harmonic noise estimates. 64 W I N D O W L E N G T H ( m s ) Figure 3.4: (a) Errors in the amplitude and phase estimates for the 300 Hz component in a synthetic record (described in the text) as a function of the length of the estimation window, (b) Amplitude of the residual 300 Hz sinusoid that would remain, following sinusoid subtraction, as a result of the amplitude and phase errors shown in part (a). 65 50 mV 5 m V 5 m V 5 m V 5 m V 5 mV TIME (ms) Figure 3.5: Illustration of the incremental improvements in signal-to-noise obtained by successive subtractions of powerline harmonics from a raw seismoelectric record (trace A). The signal is clearly visible in trace F following the subtraction of eleven harmonics, as outlined in the text. Chapter 3: Data Processing Methods 66 Effect of an error in the fundamental frequency f0 In the preceding analysis we assumed that/o was known exactly. However, as discussed earlier,/o may deviate from its expected nominal value. As a result it is important to consider the effect of an error in frequency on the performance of the sinusoid subtraction process. For the sake of simplicity, suppose the amplitude and phase of a harmonic signal have been accurately determined but the assumed frequency differs by an amount AGO from the true frequency. The residual R(t) that remains after sinusoid subtraction is R(t) = ccos(cor+(]))-ccos[(CQ + Aco)f+<))] . The above expression may also be written in the form If Aco is small then (2co + Aco)/2 = co and R(f) may be regarded as an amplitude modulated sinusoid of frequency co and amplitude 2csin(Acof/2). The amplitude of the residual is zero at t=0 where the estimated harmonic is optimally aligned with the true harmonic. Beats, of amplitude 2c, occur at times t = m%/Ao) where m is any odd integer, and the beat frequency is equal to the frequency error Aco . Suppose one tried to suppress harmonic noise having a frequency of 60.03 Hz by subtracting a 60 Hz sinusoid having the same amplitude and phase. A frequency deviation of 0.03 Hz is large but not unreasonable for 60 Hz power transmission. After subtraction, the amplitude of the residual 60 Hz noise at times t = 0.1, 0.25, and 1.0 s would be 1.9, 4.7, and 19 percent respectively of the original amplitude c. The first beat would occur at approximately 16.7 s. It is clear therefore that the importance of an accurate frequency estimate increases with record length. Also, beating may be more prevalent in the higher frequency harmonics since an error of Acoo in the fundamental results in an error of «Acoo in the nth harmonic. If a record shows signs of beating following sinusoid subtraction an estimate of the frequency error Aco may be obtained from the beat frequency or from the residual (3.22) Chapter 3: Data Processing Methods 67 amplitude 2csin(Acor/2) at a particular time t. Alternatively, one of the frequency estimation algorithms described by Nyman and Gaiser (1983) or Hancke (1990) could be employed as the first step in the sinusoid subtraction process. Our experience with 60 Hz noise indicates that it is often worthwhile to estimate / rather than assume it is equal to its nominal value. 3.3.3 Comparison with Notch Filtering Although sinusoid subtraction was the technique I used most frequently for removing powerline noise, there were a few occasions where notch filtering was required. In these cases the fundamental frequency could not be accurately determined, or the harmonics were not sufficiently stationary to be modelled as sinusoids of constant frequency, amplitude, and phase. This problem was more common during experiments near populated and industrialized areas on the Fraser Delta than it was in more isolated areas such as the U B C Research Forest in Haney. Significant changes in noise amplitude and character (e.g. the relative amplitudes of various harmonics) were sometimes observed every few seconds during particularly unstable periods. A zero phase, frequency domain notch filter program was developed for these situations. Any number of notches could be specified and the amplitude response of the filter decayed gradually to zero about each notch frequency. The rate of decay was determined by an adjustable notch width. For best results, the input data trace r(t) was truncated to a length x that was an integer multiple of l/f0 ; this ensured that the discrete frequencies in the FFT representation of r(r) included f0 and its harmonics. A practical consequence of this was that the number of points in the truncated record was rarely optimal for efficient computation of the FFT. As a result, notch filtering tended to be much slower than sinusoid or block subtraction. Nonetheless, the speed penalty was occasionally justified by superior noise rejection in areas of unstable harmonics. In Figure 3.6 we have used a synthetic record to compare the performance of block Chapter 3: Data Processing Methods 68 subtraction, sinusoid subtraction, and a notch filter. Trace A is the sum of the record's nonharmonic components s(t) and e(t). The signal is a pulse with a dominant frequency of 50 Hz and an arrival time of 85 ms. The nonharmonic noise includes a random component and a noise spike at 15 ms. The complete synthetic record (B) was obtained by adding harmonic 60 Hz interference to A. This interference was suppressed in four different ways to produce the results labelled C-F. A l l traces are plotted at the same scale and have a sample interval of 0.5 ms. Block subtraction (trace C) was carried out by subtracting the first 85 ms of trace B from the data between 50 and 135 ms. This process eliminated the 60 Hz noise but, as expected, increased the level of random noise and shifted the noise spike into the processed interval. The best result, trace D, was produced by sinusoid subtraction using an estimation window that covered the first 66.5 ms of the record. The noise spike and random noise within this window had no significant effect on the harmonic noise estimate. Trace E was also produced by sinusoid subtraction but, in this case, the estimation window covered the whole record (200 ms). The signal is evident but residual 60 Hz noise remains. The presence of the signal within the estimation window introduced errors e\"s and which contaminated the amplitude and phase estimates for the harmonic noise. Because the estimation window covered the whole trace, the sinusoid subtraction process was, in this case, equivalent to filtering by a very narrow notch filter. Trace F resulted from application of the frequency domain notch filter described above. The notch is broader and decays more gently than the one responsible for trace E. As a result, the oscillatory ringing that precedes and follows the signal is less extensive. Nonetheless, the ringing and signal attenuation are undesirable; attenuation of the second peak is particularly obvious. 69 0 40 80 120 160 200 TIME (ms) Figure 3.6: Comparison of the performance of three harmonic noise suppression methods on a synthetic data example. The synthetic record (B) consists of the nonharmonic components shown as trace A, plus 60 Hz noise. The remaining traces are the results of block subtraction (C), sinusoid subtraction with two different estimation windows (D and E), and notch filtering (F). Chapter 3: Data Processing Methods 70 3.3 .4 C o n c l u s i o n s Block subtraction and sinusoid subtraction are two processing techniques that, given a suitable estimation interval, are capable of suppressing multiple powerline harmonics without attenuating or distorting the signal of interest. Block subtraction is simpler and computationally faster because the harmonic noise estimate is recorded rather than calculated. However, it is suitable only for records containing an interval over which the nonharmonic components are negligible compared to the signal of interest. Sinusoid subtraction is more robust and is used routinely to subtract 10 to 25 harmonics of 60 Hz from seismoelectric records. Errors associated with the presence of multiple harmonics are ehminated by using an estimation window of length m/fo where m is any positive integer. Both methods have been used to obtain major improvements in the signal-to-noise ratios of seismoelectric records. Chapter 4 Field Experiments 4.1 INTRODUCTION Field experiments were the focal point of this research. In fact, it was the unexpected discovery of a remarkable seismoelectric effect at Haney, BC, in 1991 that motivated the choice of this thesis topic. U B C researchers recognized that the Haney site offered the opportunity for detailed investigations of a seismoelectric effect under conditions of unprecedented experimental control. As a result, boreholes were drilled, and a series of field experiments were carried out to methodically investigate the site over a period of 3.5 years. The resulting case study (Butler et al. 1994, 1996) clearly demonstrates the potential utility of seismoelectric effects as a geophysical tool, and constitutes one of the main contributions of this research. Other topics were examined largely because they were required to support the field experiments or interpret their results; the data acquisition and processing techniques in the preceding chapters were developed to enable reliable and low-noise measurements, while the theoretical modelling of Chapter 5 was carried out to offer an explanation for the observed signals. The experiments and surveys carried out for this thesis at several sites are listed chronologically in Table 4.1. The near-surface geology at each site was relatively simple and well-known, and the experiments involved attempts to measure seismoelectric responses from known interfaces within unconsolidated sediments at depths of 20 m or less. A representative subset of the results is presented in this chapter. This field effort would not have been possible without the assistance of many other people. In particular, I would like to emphasize the contributions of time and expertise made by the other full-time members of the U B C seismoelectric group - R. Don Russell, Anton Kepic, and Michael Maxwell. The cooperation of the U B C Faculty of Forestry in granting access to the Malcolm Knapp Research Forest at 71 Chapter 4: Field Experiments 72 Haney, B C is also greatly appreciated. The most interesting results were obtained at our Haney test site where two types of seismoelectric signals were observed. The primary response was generated by the arrival of seismic waves at a shallow boundary between road fill and glacial till. Sledgehammer and blasting cap sources positioned up to 7 m away from that interface produced clear seismoelectric conversions which were observed essentially simultaneously by widely separated dipoles at the surface. Secondary seismoelectric signals, generated in the immediate vicinity of each dipole dominated recordings made by dipoles relatively distant from the shot. Our studies show that the primary seismoelectric response exhibits the same symmetry and spatial decay rate as the field due to a vertical electric dipole beneath the shotpoint. They also show that the response cannot be attributed to piezoelectricity or to resistivity modulation in the presence of a telluric current. We infer that seismically-induced electrokinetic effects are a likely mechanism for the seismoelectric conversion. Experiments carried out at Base Borden, Ontario, and on the Fraser Delta near Vancouver, B C are also described in this chapter. At Base Borden, the interfaces of interest were the water-table, and a sand-clay contact at depths of 3 m and 9-10 m respectively. Several different sites were investigated on the Fraser Delta. There, the targets were sand-clay or sand-silt contacts at depths of up to 20 m (see Table 4.1). The data from Base Borden and the Fraser Delta clearly contain seismoelectric signals but they are not generated at boundaries. Rather - like the secondary responses at Haney, and the signals described by Ivanov (1939, 1940) - they are observed as seismic waves pass by each dipole sensor. These signals are of some interest in their own right, and are useful for comparison to the order of magnitude estimates for electrokinetic amplitudes presented in Chapter 5. Any seismoelectric effects generated at interfaces were either too small to be measured or were obscured by the signals generated in the immediate vicinity of each dipole. Stronger seismic sources (or lower noise recordings), and larger arrays of dipoles will be required to detect interfaces at those sites. Chapter 4: Field Experiments 73 Table 4.1. Chronological list of field studies carried out (in whole or in part) for this research. Filename prefixes used in the data archives are also given. Date Site Filename Prefix Comment Feb. 28/91 Haney S H A K E Instrumentation tests, electrode performance/shaking tests Mar. 14/91 Haney R O L L Curiosity piqued by an interesting seismoelectric effect observed during further instrumentation tests. Apr. 11/91 Haney R A T T L E Cross-road traverse yields strong evidence that the seismoelectric effect originates at a subsurface interface. Jan. 18/92 Haney Seismic refraction survey (with undergraduate John Rennie) Feb. 3/92 Haney Resistivity survey Feb. 4/92 Haney H A N C Test pitting, cross-road traverse, and signal polarity tests Apr. 15/92 Haney H A N D Confirmed radial symmetry of the seismoelectric response, and carried out a traverse along the road. July 9/92 CFB Borden, Ont. D N A Measurements of seismoelectric effects at site of a controlled spill of D N A P L contaminants operated by U. of Waterloo. Results inconclusive - limited by site dimensions. July 10/92 July 11/92 CFB Borden, Ont. WT Measurements of seismoelectric effects in a natural sandy aquifer. (See section 4.3 of this thesis for details.) Aug. 20/92 Aug. 24/92 Haney H A N E Drilling of first two boreholes (BH-92A and BH-92B) reveals existence of the fill/till boundary at 1-3 m depth. Apr. 22/93 Highway Ramp, Richmond, B C (Fraser Delta) RP Attempt to measure seismoelectric effect from base of sandy fil l on an abandoned highway ramp at corner of Miller Rd. and Russ Baker Way. Ambiguous results (insufficient data). Apr. 29/93 Haney H A N F Measurement of seismoelectric response vs. dipole offset May 13/93 Haney H A N G Measurement of seismoelectric response vs. shot depth with blasting caps in borehole 92B. May 20/93 Haney 3 component seismic data acquired with Dr. Jim Wright June 17/93 Haney H A N H Seismic velocity measurements with hydrophones in BH-92A June 24/93 Haney HANI Seismoelectric response vs shot depth experiment in BH-92A Nov. 25/93 Haney HANJ Experiment to test for linearity of the seismoelectric response (with undergraduate student Sean Fleming) June 22,30/94 Haney 5000,6000 Seismic refraction survey July 21/94 Haney Drilling of second set of boreholes (BH-94A and BH-94B) July 26/94 Haney H A N K , 7000 Seismic velocity measurements with hydrophones in new boreholes. Sept. 29/94 Oct. 6/94 Haney H A N L Measurements of seismoelectric response vs. dipole offset with shots in the new boreholes. Table continues on next page Chapter 4: Field Experiments 74 Table 4.1 continued Date Site Filename Prefix Comment Oct. 20/94 Oct. 27/94 Haney H A N M Measurements of seismoelectric response vs. shot depth in deep (10m) borehole 94A demonstrate conclusively that conversion occurs at the fill/till boundary. July 24/95 Delta, B C (Fraser Delta) F D A Site of excellent seismic reflection data acquired by U B C student Kevin Jarvis along 57B St. near Delta Port Way. July 26/95 Arthur Laing Bridge Richmond, B C (Fraser Delta) FDB Site of GSC borehole FD90-1. Sand/clay interfaces at 20 m and 29 m depth. Very high electrical conductivities (500 mS/m between 4 and 20 m depth) probably diminished any seismoelectric responses from those boundaries. Aug. 2/95 Sept. 9/95 Sept. 14/95 M O T Tower Site Richmond, B C (Fraser Delta) FDC F D E Site of GSC boreholes FD92-3 and FD94-3. Holocene/Pleistocene boundary at approximately 20 m depth. (See section 4.4.1 of this thesis for more details.) Aug. 20/95 Tssawassen, B C (Fraser Delta) FDD At N W comer of 52nd St. and 16th Ave., near GSC borehole FD86-1 which intersected the boundary between Holocene sands/silts and Pleistocene till at about 8 m depth. 4.2 EXPERIMENTS AT HANEY, B C 4.2.1 Site Description The Haney site, near Vancouver, Canada, is located on an unimproved dirt road that runs along the side of a steep slope in the Malcolm Knapp Research Forest of the University of British Columbia. The road fill, consisting of a permeable, organic-rich soil, overlies a highly impermeable, silty glacial till. We drilled four boreholes at the site, and drill cores show that the fill varies in thickness from 1 to 3 m across the width of the road. The deepest borehole, drilled to a depth of 10.4 m, penetrated 7.7 m of glacial till without encountering any obvious change in lithology. Attempts to reach bedrock were thwarted by difficult drilling conditions (abundant boulders) in the deeper part of the till layer. However, seismic refraction data indicates that bedrock lies at a depth of 10 to 15 m. The P-wave velocities of the fill, till, and bedrock are approximately 250,2100, and 4000 m/s respectively. Resistivity soundings carried out on the road using a Wenner array indicate that electrical conductivities are approximately 0.4 mS/m in the road fill, and of the same order of magnitude in the glacial till. Chapter 4: Field Experiments 75 4.2.2 Results and Discussion Figure 4.1 shows seismic and electrical responses to a single sledgehammer blow on the road at Haney. The uppermost trace is from an accelerometer attached to the head of the hammer, and simply gives the moment of hammer impact. Traces 2 and 3 are signals from the geophones located at the electrodes closest to the shotpoint. The final five traces show the time-varying potential difference across each dipole (i.e. the potential at the northern electrode minus that at the southern electrode). Powerline noise has been removed from the electrical traces by sinusoid subtraction A clear seismoelectric response is visible on the dipole traces. It arrives simultaneously on all five dipole traces 6 ms after hammer impact. Since it precedes the seismic arrivals, the response cannot be attributed to shaking of the dipole electrodes, or to geophone cross-talk. The simultaneous arrival at different dipoles is consistent with the idea that the response is generated at depth and propagates rapidly to the surface as an electromagnetic signal. Note that the polarity of the response is reversed on opposite sides of the shotpoint. This shows that electrodes near the shotpoint initially detected a drop in electric potential relative to the more distant electrodes. The peak to peak magnitude of the response is about 2.4 mV across dipoles D1-D4, and 0.8 mV across the more distant dipole (D5). The early cycles of both the electrical and the geophone signals exhibit similar dominant frequencies (70-100 Hz) but the electrical frequencies appear to be slightly higher. The later, low frequency cycles in the geophone data are surface waves. Figure 4.2 shows the electrical responses measured by dipole D2 as the shotpoint was moved in 0.5 m increments across the road at Haney. The line of shotpoints lay between the 0 and 6 m marks in Figure 4.1(b). The seismoelectric signal arrived earliest (about 4 ms after hammer impact) at the western edge of the road, and the delay gradually increased to about 14 ms as the shotpoint was moved to the east. This suggests that the distance to the seismoelectric source increased as the shotpoint was moved to the east. 76 Figure 4.1. (a) A shot gather showing the seismic ( G l , G2) and electrical (D1-D5) responses to a single sledgehammer blow at the Haney site. The dotted line indicates the moment of sledgehammer impact as determined by the accelerometer trigger ( A C C ) . Powerline harmonic noise has been removed from the dipole traces by sinusoid subtraction, (b) Plan view of the experimental layout on the road (shaded area) showing the shotpoint (SP), two geophones, and five 10 m dipoles. The geophones are located within a few centimeters o f the electrodes closest to the shotpoint (i.e., 3 m from the shotpoint). 77 < CM X m Figure 4.2. Seismoelectric signals generated at 13 shotpoints distributed across the dirt road. Time zero indicates the moment of hammer impact. The cross-section below shows the estimated position of the seismoelectric target (a dipping interface tangent to the arcs), as well as the actual depths to the road fill/glacial till boundary in two boreholes Chapter 4: Field Experiments 78 The interpretation of the first arrival is also shown in Figure 4.2. At each shotpoint, an arc has been drawn with a radius equal to the product of the seismoelectric delay and the P-wave velocity of 250 m/s. Neglecting out-of-plane effects, the feature responsible for the seismoelectric conversion should be tangent to all of the arcs. The seismoelectric data therefore delineate a dipping interface, 1 m deep in the west and 3.5 m deep at the eastern edge of the road. As shown in Figure 4.2, two shallow boreholes confirmed the existence of a dipping boundary between road fill and glacial till. The 40 to 50 cm discrepancy between the actual and estimated depths can be ascribed to uncertainty in the first break picks (1 to 2 ms), and in the P-wave velocity of the road fill. As a further check on the identity of the target, we measured the seismoelectric responses generated by the detonation of blasting caps at various depths in a borehole. Figure 4.3 shows the experiment layout and a typical shot record (the shot at 7.7 m depth). A l l measurements were made at the surface using one uphole geophone and six dipoles located 2, 4, 6, 8, 12, and 16 m north of the borehole. The uppermost trace is from the fiber optic trigger circuit and gives the time of detonation. The second trace shows the seismic arrival at the surface, and remaining traces are electrical responses measured by the dipoles. The responses differ in character but the first arrival occurs simultaneously at all six sensors 2 ms after detonation. This arrival is clear at near offsets but barely visible above the noise at the most distant dipole. The 2 ms delay corresponds to the time required for the seismic wave to travel 5 m up to the base of road fill where the seismoelectric conversion takes place. There is a further delay of 13 ms before the seismic wave reaches the geophone; this is the time required for the seismic wave to travel through 2.7 m of fill. In Figure 4.4 we have plotted the seismoelectric responses measured by dipole D3 as the shotpoint was moved up from the bottom of the borehole. The most striking feature is the abrupt change in signal character that occurs opposite the fill/till interface. Blasting caps detonated below the boundary yielded higher amplitudes and frequencies. We speculate that 79 a) T 2.7 m D1 D2 D3 D4 D5 D6 geophone 2 m road fill < O I x DQ 4r s h o t glacial till N 10 Detonation Time 10 20 30 TIME (ms) 40 50 Figure 4.3. (a) Illustration of the experimental layout used to measure seismoelectric responses to the detonation o f blasting caps in a borehole. One geophone and six 2 m dipoles were deployed on the surface, (b) Response observed for a shot at 7.7 m depth. Chapter 4: Field Experiments 80 this is indicative of better seismic coupling in the dense, competent glacial till than in the loose, highly compressible road fill. The voltage spike at time zero on the trace at 8.7 m was caused by the use of a particularly narrow steel blast chamber. Blast chambers were not used for any of the shallower shots. The main point illustrated by Figure 4.4 is that seismoelectric conversion clearly occurs at the fill/till interface. As expected, the delay between the instant of detonation and the reception of a response was proportional to the distance between the shot and the road fill/glacial till boundary. There was no delay when the shot was located direcdy opposite the boundary. The first arrivals can be fit well by two straight lines that intersect at time zero near the fill/till interface. The slopes of these lines - 200 and 2300 m/s - provide estimates of the seismic velocities in the road fill and glacial till respectively. These are in reasonable agreement with the average P-wave velocities of 250 and 2100 m/s derived from seismic measurements at the site. Given that electrokinetic effects involve the motion of pore water, the state of water saturation in the subsurface is of obvious interest. We have monitored the position of the water table at Haney by measuring the natural water level in boreholes. These measurements ranged from about 0.5 m above, to a few meters below the fill/till boundary depending on the season and the road fill thickness at the point of measurement. Our studies indicate that the seismoelectric effect is present, and originates at the same point (the base of fill) regardless of the water level. For example, the natural water level was 1 m below the fill/till interface during the experiment displayed in Figure 4.4, but the same type of data were acquired in a different borehole the previous year when the water table was 0.35 m above the interface (see Figure 4.5). In both cases shots detonated at the boundary generated immediate electrical responses while those detonated at the water table yielded delayed signals. 81 TIME ( m s ) Figure 4.4. Seismoelectric response vs. shot depth in borehole 94A. These electrical signals were measured by a 2 m dipole on the surface 6 m north of the borehole. Note that the delay between detonation and response depends on the distance to the fil l / t i l l interface at 2.7 m depth. The peak to peak amplitude of each trace is indicated on the right. 82 T I M E (ms) Figure 4.5. Seismoelectric response vs. shot depth in borehole 92B (one of the boreholes shown in Figure 4.2). The signals were measured by a 2m dipole on the surface 6m north of the borehole as shown in the accompanying cross-section. These data; like those in Figure 4.4, confirm that seismoelectric conversion occurs at the lithologic (fill/till) boundary. 2m d i p o l e ^ v B H - 9 2 B 1 j I 2.8m road fill g lacial till Chapter 4: Field Experiments 83 Seismoelectric Amplitude vs. Offset We have also carried out experiments to determine how the seismoelectric response varies with dipole offset. Figure 4.6 shows a comparison of seismic and seismoelectric arrivals at offsets of 2 to 26 m from a shotpoint on the road surface, about 2 m above the fill/till inter-face. Thirteen 2 m dipoles and twelve vertical geophones (40 Hz geophones manufactured by O Y O Geospace) were laid out along the road to the north of the shotpoint as shown on the map. The two data sets were collected separately but with the same sledgehammer source. The seismic data appear to be dominated by surface waves and it is difficult to identify any meaningful arrivals apart from the first breaks. Beyond a distance of 6 m, the breaks are due to a head wave traveling along the interface between road fill and glacial till. The rise in the first breaks at an offset of about 14 m is probably due to localized thinning of the road fill. In contrast, the seismoelectric record exhibits negligible surface wave interference, and at least three coherent arrivals. The primary arrival appears 9 ms after hammer impact - essentially simultaneously on all traces out to about 15 m offset. This is the seismoelectric response produced by the arrival of the seismic wave at the road fill/glacial till boundary below the shotpoint. Beyond 15 m, it is too weak to be clearly identified. However, two secondary sub-parallel and non-simultaneous arrivals are evident. The secondary signals appear in a 40- 50 ms time window beginning 2 to 3 ms before the first seismic arrival at each dipole. Furthermore, their apparent velocities are comparable to those of seismic body waves. This suggests that they may be associated with the arrival of seismic waves beneath each dipole - beginning with the head wave that travels along the fill/till boundary. The two clear secondary events that appear between 30 and 60 ms exhibit moveout that is roughly hyperbolic; we speculate that they might be generated near each dipole by seismic waves that have been reflected from the till/bedrock interface. Unfortunately, any reflections of that type are obscured by surface waves in the seismic record. At this time, the origin of the secondary arrivals is not well understood. 84 Chapter 4: Field Experiments 85 The seismoelectric traces in Figure 4.6 have been normalized for display purposes. As shown in Figure 4.7, the true peak to peak voltages measured across the 2 m dipoles varied by two orders of magnitude from a high of almost 3 mV at close range, to 20 | l V at the farthest offsets. The amplitude of the primary seismoelectric response (the simultaneous first peak) is also plotted for the offset range where it is visible. Beyond an offset of 4 m, both amplitude curves can be approximated by straight lines in the log-log plot. The slopes of these lines indicate that peak to peak amplitudes decay approximately as l/x2, while the amplitude of the primary signal falls off approximately as l/x 4 , x being the dipole offset. The latter rate of decay is the same as would be exhibited by the horizontal component of the electrostatic field from a vertical electric dipole at the fill/till boundary beneath the shotpoint. The source of the primary seismoelectric response can therefore be modeled as a buried vertical dipole for offsets greater than 4 m. 0 5 10 15 20 25 30 1 10 100 Dipole Offset (m) Dipole Offset (m) Figure 4.7. Linear and log-log plots of seismoelectric amplitude vs. offset for the experiment shown in Fig. 4.6. The solid symbols give the overall peak to peak amplitude of each trace, while the open sym-bols give the amplitude of the primary (simultaneous) arrival. The dashed lines in the log-log plot have slopes of -2 and -4. Notes: (i) offset is taken as the distance from the shotpoint to the near electrode of each dipole; (ii) the plots include data points from dipoles at 3, 5, 7 and 9m offset not shown in Fig. 4.6 Linearity of the Seismoelectric Response As mentioned in Chapter 1, seismoelectric effects can be divided into 2 major categories - linear and non-linear. The presence or absence of linearity is important because it Chapter 4: Field Experiments 86 constitutes an important constraint on possible conversion mechanisms. In 1993/94, the linearity of the seismoelectric response at Haney was investigated by Sean Fleming (then an undergraduate student in the Geophysics Program at UBC) with assistance from the author, R.D. Russell, and A.W. Kepic (Fleming, 1994). The experiment involved measuring both seismic and seismoelectric responses to the drop of a 22.7 kg (50 lb.) lead weight from heights of about 10 cm to 2.5 m above the road surface. The weight was suspended from a rope which passed through a pulley tied to the top step of a large step-ladder. As shown in Figure 4.8(a), seismic signals were measured by a vertical geophone buried near the shot while seismoelectric responses were measured by four 2 m dipoles oriented radially about the shotpoint. In order to obtain optimum signal-to-noise ratios (important for the drops from very low heights), we took the difference between two dipoles on opposite sides of the shotpoint. This process doubled the signal and cancelled much of the regional noise due to powerlines and other sources. The resulting 'quadxupole' trace was also subjected to sinusoid subtraction in order to remove the remaining powerline harmonic noise. Figure 4.8(b) shows that the primary seismoelectric response at Haney does indeed vary linearly with the seismic input. Each point in the plot represents data from a separate weight drop. The abscissa give the rms (root-mean-square) particle velocity measured by the geophone over a 120 ms interval following each drop, while the ordinate gives the rms voltage measured by the quadrupole sensor over the same interval. The locus of data points is well-approximated by a straight line having a slope of 1 on the log-log plot. The linearity of the response is not surprising, but this experiment confirms that we can rule out nonlinear conversion mechanisms in our search for the origin of the seismoelectric effect at Haney. Chapter 4: Field Experiments 87 0.01 Figure 4.8. (a) Plan view of the experimental layout used for the linearity test on the road surface at Haney. The geophone was buried at a depth of 0.7 m. (b) Log-log plot of the seismoelectric response (rms quadrupole voltage) vs. the seismic input (rms particle velocity) generated by a series of weight drops from various heights. The straight line has a slope of one (in log-log space), showing that the seismoelectric amplitude is proportional to the seismic amplitude. 4.2.3 The Seismoelectric Conversion Mechanism The preceding experiments have served to identify the seismoelectric target and put some constraints on possible models for the conversion mechanism. First of all, we conclude that the conversion is not a piezoelectric effect. Although the glacial till contains quartz, there is no reason to expect that the quartz-rich grains would have been deposited with the alignment required to produce a measurable piezoelectric response. Furthermore, i f the till was piezoelectric, blasting caps detonated within it should have produced immediate electrical responses regardless of the distance to the fill/till boundary. Resistivity modulation is another mechanism that must be considered. According to this model, the resistivity of a volume of earth varies with stress during the passage of a seismic wave. Electric potentials due to any natural currents flowing in that volume therefore vary as well. To investigate this possibility we used an experiment devised by Ivanov (1939). We measured electrical responses to seismic sources located 3 m away from either end of a dipole Chapter 4: Field Experiments 88 sensor. To first order, we can assume that both shots changed the effective resistance between the dipole electrodes in the same way. Then, given that the ambient telluric currents were expected to be essentially unidirectional beneath the dipole, any telluric potential drop should have varied in the same way (had the same polarity) regardless of whether the shot was to the left or the right. However, the two shots actually yielded responses with opposite polarities, indicating that resistivity modulation is not the relevant mechanism at this site. Figure 4.9 shows the results of a more detailed investigation of signal polarity. Seven dipoles were arranged in a radial pattern about the shotpoint. Each measured the potential at its outer electrode relative to that at its inner electrode. Apart from some early source-generated noise (the dipoles were very close to the shotpoint), the main feature is the seismoelectric arrival at 10 ms. The signal polarities indicate that the response (measured at the surface) begins with a flow of current radially inward towards the shotpoint. In the absence of a current source or sink, this net horizontal flow toward the shotpoint must be balanced by a net vertical flow downward beneath the shotpoint. Again, we cannot envisage any likely scenario by which resistivity modulation could cause uniform telluric currents to change in this fashion. One interesting alternative is suggested by the observation that natural potential differences of up to a few tens of mV can exist between adjacent layers of soil depending on the degree of oxidation or reduction occurring in them (Mitchell, 1993, p. 283). Suppose such a potential difference existed between the organic-rich road fill and the underlying glacial till. If so, there would be a nearly vertical electric field across the dipping fill/till interface*. If a seismic wave, incident from above, were to modulate the resistivity or change the potential distribution in any way, then the vertical electric field would be expected to change in a way * Vertical electric fields could also be produced by fluid rxercolating naturally across the interface (i.e. by a natural streaming potential (Sill, 1983)) but this is considered less likely at Haney because the glacial till appears to be highly impermeable. <8^> \\ 2 m b) TIME (ms) Figure 4.9. (a) Plan view of seven 1.5 m dipoles arranged in a radial pattern to measure the directionality of the seismoelectric response, (b) This record, generated by a hammer blow at the center of the pattern, shows that the seismoelectric signal is approximately radially symmetric at the surface. The accelerometer trace (top) indicates the moment of hammer impact. Chapter 4: Field Experiments 90 that is radially symmetric about the shotpoint - in agreement with the signal symmetry shown in Figure 4.9. This idea is of practical interest because large potential differences between soil layers can cause water to flow through electroosmosis and accumulate at an interface, leading to slope stability problems if the soils lie on a hillside (Mitchell, 1993). A seismoelectric method for mapping the depth to possible failure planes would be a valuable tool for geotechnical investigations. However, the idea is speculative and untested at this point; measurements of electric potential vs. depth in a borehole or trench at Haney would be required to determine if there is a substantial potential difference across the fill/till interface, and a proper theoretical model is required to estimate the signal strengths that might be obtained by this mechanism. Seismically-induced electrokinetic effects appear to offer the best explanation for the conversion mechanism. The rigorous theoretical treatment of this problem is complicated, but a useful qualitative explanation can be given by the simple electrostatic model proposed in Chapter 5. The model accounts for the two types of arrivals present in our field data, and provides physical insight into the nature of the conversion process. 4.2.4 Summary of the Haney Experiments We have measured seismoelectric responses from a sedimentary boundary that are unusual in their clarity and detail. The boundary was mapped by experiments involving multiple shotpoints, and a few stationary electric field receivers. Other experiments, involving one shotpoint and an array of receivers revealed the existence of two types of seismoelectric arrivals. Information on the nature of the conversion mechanism has been obtained by measuring the rate of signal decay with offset, by determining the symmetry of the response, and by confirming its linearity. These studies and other observations show that the response cannot be attributed to piezoelectricity or to the modulation of resistivity in the presence of a uniform telluric current. Chapter 4: Field Experiments 91 The boundary mapped at Haney is an interface between permeable, organic-rich road fill and impermeable, silty glacial till, and seismically-induced streaming potentials are believed to be responsible for the seismoelectric conversion. Given the high signal-to-noise ratio obtained using a relatively weak seismic source, we expect the effective depth of exploration for other suitable interfaces to be much greater than the 3 m probed at this site. 4.3 EXPERIMENTS AT BASE BORDEN On July 10 and 11,1992, seismoelectric measurements were made at a well-characterized site within the University of Waterloo's hydrogeology research station at Canadian Forces Base Borden near Alliston, Ontario. The surficial geology consisted of a clean sand aquifer overlying a silty-clay aquitard at a depth of 9-10 m. The targets of interest were the water table at approximately 3 m depth, and the sand-clay interface. Some relevant physical properties have been measured by other researchers at this site. Conant (1991) measured hydraulic conductivities in core samples and quantified the permeability contrast at the sand-clay boundary. His results indicate that average values for permeability are about 6 x 10\" 1 2m 2 (6 darcy) in the aquifer, and at least 2 to 3 orders of magnitude less in the aquitard. Gilson and Redman (1995) measured electrical conductivity in a borehole with an inductive logging tool (the Geonics EM39). They found conductivities of about 1-2 mS/m in the sands above the water table and 10-12 mS/m in the saturated sands below. Figure 4.10 shows one of the experimental layouts used at this site. The shotpoint (SP3) was located in a trench 1 m deep, and the water table depth, as measured in a well about 30 m away, was 3.4 m. It is expected therefore that the water table lay about 2.4 m below the shotpoint. The six dipoles were oriented as shown because the powerline noise was weakest in that direction. Four vertical geophones were deployed at the southern electrodes for dipoles 1-4 in order to determine when the electric field measurements might be affected by electrode shaking. Chapter 4: Field Experiments 92 Q Q Arrow indicates survey north (not true north). (§)<$)<£)<£) (D 2m N A shotpoint S P 3 T T T T geophones Figure 4.10. Plan view of one layout used to make seismoelectric measurements in sandy sediments at Base Borden in July, 1992. Dipoles D1-D5 were 5 m long; D7 was 3 m long. The polyethylene tank sitting on the surface was being prepared by University of Waterloo geophysicists for another study. It is shown for geographical reference only. Figure 4.11 shows the seismic and electrical arrivals generated in response to a stack of hammer blows at SP3. All traces are mean stacks of 30 hammer blows and the dipole traces have been processed by sinusoid subtraction to suppress powerline noise. The peak-to-peak voltage gradients or electric fields across the dipoles ranged from 50 to 15 LtV/m. The seismoelectric signals are similar to the secondary responses observed at Haney in that their arrival times clearly increase with dipole offset. In fact, the electrical signal appears to travel with the seismic wave at an apparent velocity of 350-400 m/s. It is significant however that the electrical arrival precedes the seismic arrival at each dipole by 2-3 ms. This shows that the seismoelectric response cannot be attributed merely to shaking of the dipole sensors by the seismic wave. It could represent a seismically-induced electrokinetic effect that travels with the seismic wave through the sandy sediments. A conceptual model for this effect is given in the Chapter 5. A comparison of the seismoelectric data in Figure 4.11 to that from Haney (Figure 4.6) reveals two significant differences. First, the peak electric field strength at Borden is a factor of about 30 lower. Secondly, there is no clear simultaneous electrical arrival in the Borden data. This suggests that the water table and sand-clay interface were poor seismoelectric Chapter 4: Field Experiments 94 targets compared to the fill/till interface at Haney. Simultaneous signals may have been generated at these targets but, if so, they were weak and obscured by noise or by the other seismoelectric response that travels with the seismic wave. 4.4 EXPERIMENTS ON THE FRASER DELTA The Fraser River Delta, located immediately south of the city of Vancouver, BC, is the largest delta on the west coast of Canada. It is comprised primarily of Holocene sands, silts and clays deposited during the past 10 000 years over glacial and interglacial (Pleistocene) sediments (Lutemauer et al., 1995). Much of the delta is urban and industrialized, and because the earthquake risk is considered moderately high, the geologic and geotechnical characteristics of the Fraser Delta have been the focus of considerable study by the Geological Survey of Canada (GSC) - particularly since the mid-1980's (e.g. Clague et al., 1991; Hunter et al., 1994; Dallimore et al., 1995). Their methods have included extensive drilling and coring, stratigraphic analysis, high resolution seismic reflection profiling, seismic cone penetrometry, surface and downhole geophysical logging of electrical conductivity, magnetic susceptibility, and natural gamma radiation, and mapping of the depth to bedrock using industry seismic data. The detailed information provided by these studies at sites scattered across the delta constitutes an excellent database for researchers such as myself interested in the development of shallow geophysical methods. On a regional scale, the deltaic Holocene sediments are now known to reach depths of at least 236 m. The bedrock surface lies at an average depth of about 500 m but exhibits considerable topography and plunges to depths in excess of 1000 m in some areas. Seismoelectric experiments were carried out at four different sites on the Fraser Delta during the summer of 1995 (see Table 4.1). The objective at each site was to determine if known interfaces between naturally deposited sands and silts or clays at depths of 5-20 m could be detected by the seismoelectric method. Data were acquired for one to three days at Chapter 4: Field Experiments 95 each site. Seismoelectric signals were measured, but like those at Base Borden, they were associated with the arrival of seismic waves at each dipole sensor. Responses from interfaces were not observed. The characteristics of the Fraser Delta data are well-represented by the example given below from the M O T Radio Tower site in Richmond. FD92-3 (T.D. 41 m) (approximate location) N A 100 m •o CO 2 2 CD 840m to No. 6 Rd. ® + FD94-3 (T.D. 305 m) MOT Radio Tower Compound Location of Dipoles used for seismoelectric survey 800m to No. 7 Rd. Westminster Hwy Figure 4.12. Map of the MOT Radio Tower site on the Fraser Delta in Richmond, BC. The two boreholes lie within approximately 100 m of the location of the seismoelectric survey. 4.4.1 The MOT Tower Site, Richmond This site was chosen for its relatively simple geology and for the excellent geological control provided by two GSC boreholes FD92-3, and FD94-3 drilled to depths of 41m and 305 m. It is located next to a Ministry of Transport (Transport Canada) radio tower, in an agricultural area of Richmond, BC, between No. 6 and No. 7 Roads just north of Westminster Highway. The seismoelectric survey was carried out within about 100 m of the boreholes and the radio tower in a dry ditch beside the gravel drive leading to the tower compound (see Figure 4.12). According to samples from the two boreholes, the generalized geology is as Chapter 4: Field Experiments 96 shown in Figure 4.13. A surficial layer of silt and silty sands extends to a depth of 1-3 m, and is underlain by a thick package of loose, well-sorted, fine to medium-grained sand with occasional clayey-silt interbeds. At 19 m there is a layer of clay and cobbles about 1 m thick, underlain by a more compact sand which extends to a depth of at least 40 m. Borehole FD94-3 penetrated through an additional 265 m of sands, silts, and clays and was terminated at 305 m without encountering bedrock. The clay and cobble layer at 19 m is believed to mark the boundary between the post-glacial Holocene sediments deposited by the Fraser River, and the underlying Pleistocene section. In fact, this site is of special geological and geotechnical interest because the Holocene sediments are anomalously thin (Dallimore et al, 1995). Shear wave refraction data acquired by the GSC indicate that the Holocene/Pleistocene boundary dips gently to the south reaching a depth of approximately 30 m at Westminster Highway. The depth beneath our survey location is therefore probably 20-25 m (pers. comm., Dr. Jim Hunter, GSC). The water table depth is not given in the borehole logs, but a depth of 5-7 m is suggested by seismic refraction data that I acquired at the site. A wide variety of geological, geophysical, and geotechnical logs are available for the boreholes at this site (Hunter et al., 1994; Dallimore et al., 1995). Electrical conductivities obtained from borehole logs are approximately 30 mS/m in the thin silty layer at the surface, and 10 mS/m in the Holocene sands. Conductivity increases below the clay/cobble layer reaching a value of 100 mS/m at 25 m depth and 300 mS/m by 40 m depth. P-wave velocities obtained by a seismic refraction survey were approximately 700 - 800 m/s above the water table, and 1530 m/s in the saturated sands. Borehole measurements of P-wave velocity (only available below the water table) were 1500-1570 m/s in the Holocene sands, and only modestly higher - about 1600 m/s - in the Pleistocene sands at 20 to 40 m depth. However, it is interesting to note that borehole measurements of S-wave velocity almost doubled at the Holocene/Pleistocene boundary changing from about 200 to 400 m/s. This indicates that the rigidity of the Pleistocene sands is significantly greater than it is in the Holocene. 97 MOT Tower Site, Richmond horizontal «-» and vertical t geophones shotpoints 4 ® L 2m dipole 4 ^ o=30mS/rn silt, silty sands 1-3m v p=700m/s o=10mS/m 5-7 m . . 1 water table v p=1530 m/s o=10mS/m / oose sand •20mdepth 1m I v p=1600m/s a=100mS/m compact sand FDELAY3.CDR Figure 4.13. Cross-section o f the near-surface geology at the M O T Tower Site based on information from G S C boreholes FD92-3 and FD94-3 located approximately 100 m away. The experimental layout used to acquire seismoelectric data on September 9, 1995 is shown at the surface; it included three 2 m dipoles (with northern electrodes positive), three vertical geophones, and three horizontal (radial) geophones. Data for 15 different source-receiver offsets were obtained by occupying five different shotpoints wi th the buffalo gun seismic source. A remote dipole (not shown) was used to monitor the regional noise 70 m south o f the northernmost dipole. Chapter 4: Field Experiments 98 The target interface at this site was the boundary between the Holocene sands and the clay/cobble layer at about 20 m depth. The conditions which might contribute to a seismoelectric conversion at this depth include the marked increases in conductivity, and rigidity mentioned above, as well as an expected drop in permeability. Unlike the fill/till * interface at Haney, this boundary is probably a poor reflector of P-waves given that the clay layer is thin compared to typical P-wave wavelengths, and that the velocities of the Holocene and Pleistocene sands are almost the same. Another reason for the choice of this site was the fact that the conductivities above the target interface were low compared to those measured in most of the GSC's other Fraser Delta boreholes (compare the logs for 22 boreholes compiled by Hunter et al., 1994). Since seismically-induced streaming potentials (like the steady flow potentials discussed in Chapter 1) are expected to vary inversely with conductivity, a low conductivity site was considered advantageous. Seismic and seismoelectric arrivals measured as a function of offset from the shotpoint at the M O T Tower Site are compared in Figure 4.14. The seismic source was an 8-gauge buffalo gun shot detonated in a hole approximately 1 m deep. Seismoelectric signals were measured using 2 m dipoles. The seismic arrivals were recorded with both vertical geophones, and horizontal geophones which measured the radial component of particle velocity in the direction of a line running from the shotpoint to the geophone. More details on the geophone specifications are given in Table 4.2. Table 4.2 Specifications of geophones used at the MOT Tower site. Due to a shortage of horizontal geophones, two different models were used. Geophone Resonant;Frec|uency_ Transduction..Constant Vertical component (Geospace GSC 20D) 10 Hz 31 V/(m/s) Horizontal Component 2-20 m offset: Mark Products LIB 22-30 m offset: Mark Products L l 8 Hz 4.5 Hz 39 V/(m/s) 39 V/(m/s) Chapter 4: Field Experiments 9 9 The signal-to-noise ratio for electrical measurements on the Fraser Delta was much worse than that at Haney. The seismoelectric signals themselves were smaller, and the ambient noise was more troublesome. The powerline harmonic noise was unstable, and sudden changes in the relative amplitudes and phases of its harmonic components were observed quite frequently. Such changes diminished the effectiveness of the sinusoid subtraction processing technique which assumes that the frequency, amplitude, and phase of each harmonic remains constant over the length of each trace. It is also possible that the M O T antenna (a 500 W omni-directional aircraft navigation beacon for Vancouver International Airport) contributed to the noise at this site. The antenna broadcasts a repetitive A M signal with a carrier frequency of 266 kHz (pers. comm., Transport Canada Technical Services, Vancouver). Although the dipole preamplifiers were designed to reject such high frequencies and minimize the possibility of A M demodulation, some noise of this type may contaminate the records given the proximity of the transmitting antenna. In order to reduce electrical noise prior to recording, the output of a remote dipole, located approximately 75 m south of the shotpoint, was subtracted from that of the other dipoles using the Tektronix AM502 differential amplifiers at the recording site. This reduced noise levels by a factor of about 10 to 20. Additional noise reduction was obtained by stacking (averaging) up to 3 shots for each dipole trace. The recording passband was 10 Hz - 1 kHz. Following acquisition, the seismoelectric data were digitally filtered with a 6 pole, 120 Hz Butterworth lowpass filter and subjected to sinusoid subtraction to remove as much of the remaining harmonic noise as possible. The final noise level in the processed data (Figure 4.14) is approximately 3 p:V peak to peak. By comparison, the intrinsic noise level of the V-box preamplifiers used here is only 0.35 U.V rms over a 110 Hz bandwidth, and is reduced further by signal stacking. This discrepancy, together with the semi-periodic character of the residual noise, shows that the signal-to-noise ratio of the data is limited by our inability to remove all of the ambient noise in the data - not by the intrinsic noise of the amplifier. 100 u CO CD o Q. o CD a (3 o <-i as n \"c ch OZU Ho OO > 111 o CO n CD oo cz •i—i o CD _ c oo CL o T3 CD a \"cO O o • i—i t : *-( CD o > CD CD O CJ oo • i—i CD oo y ^ CO CD o O CL a T3 VJ 0 E • r H CM S-t c3 O +-* £ o CD O CD U mo CO 'CD 00 50. uV 90. uV 200. uV 100. uV 500. uV 13 £ >> +-> U -13 -a o o c \" '53 OH W • a -s S3 o s o a N « .CO S crJ S a •3 S McrJ '•B ° 3 o Ii * O o o U co 5 3 . 2 ? L 1 H 1 c o O N •60 O S I * . * 3 ! T3 • s « R T3 O u o o 'B d O co 6 a 0) c o g T3 o o c c o a . a o ° u O o o • i—I e \"53 CO o cx . 6 £ 60 & ' - 3 S g o •I 2 § • N ^ c > co ^ > o e CT O CO oj ^ & s O _ o M [> D (U (H t-i » ctp ^ is porosity K is permeability T| is fluid viscosity Chapter 5: Modelling 106 P12 is the inertial coupling term (pi 2 < 0) which depends on the pore geometry Pn = (l-)P,-Pi2 p 2 2 =( j )p / -p 1 2 and ps,Pf are the densities of the solid and fluid components respectively. The terms on the left hand side of equations (5.1) are forces per unit volume. They are balanced by two terms on the right which represent (i) the time rate of change of linear momentum in the solid and fluid (the inertial term), and (ii) a viscous drag force. The assumption that viscous drag varies linearly with the relative flow between solid and fluid (i.e. the assumption of Poiseuille flow) is valid provided the frequency is less than a certain value dependent on the viscosity of the fluid and the size of the pores (Biot, 1956a). In water-saturated coarse-grained rocks or sediments with rather large pores, the assumption may not be justified even at seismic frequencies in the range from 10's to 100's of Hz. Biot (1956b) proposed an alternate form for the viscous drag force in that case, but the low frequency/small pore approximation is used here for the sake of simplicity. The left hand side of equations (5.1) can be expressed in terms of the displacements and Ui by substituting stress-strain relations which Biot (1956a) derived by analysis of the strain energy in a porous medium. Making this substitution and combining equations (5.1) with the analogous equations for motion in the y and z directions yields the following vector wave equations in terms of u and U: N V 2 u + V[(A + N ) V - u + 2V -U] = | ^ ( P l l u + p 1 2 U ) + ^ | - ( u - U ) (5.2) V[GV • u + R V - U] = |^- ( P l 2 u + p 2 2U) - (u - U). A,N,Q, and R are elastic constants with A and N corresponding to Lame's parameters (N is the shear modulus). The physical significance of these constants, and methods for their Chapter 5: Modelling 107 measurement in the laboratory, were elaborated on by Biot and Willis (1957). The vector equations (5.2) represent, in compact form, the complete set of six equations required to solve for the six components of displacement (ux,uy,uz,Ux,Uy, and Uz) accompanying elastic wave propagation in homogeneous porous media. To simplify analysis, they may be separated into two pairs of scalar equations describing shear and compressional wave propagation by applying the curl and gradient operators respectively. (The same procedure is used to separate the vector wave equation for purely elastic media.) Biot (1956a) next sought plane wave solutions to the coupled wave equations. In the case of compressional waves, such solutions are of the form The direction of propagation/displacement can be arbitrarily chosen because the medium is isotropic, but choosing a direction parallel to one of the coordinate axes (e.g. parallel to the x-axis as in equations (5.3)) simplifies the algebra. The wavenumber k = kr + jkt is generally complex, its imaginary component specifying the attenuation due to viscous damping. The wave's phase velocity is given by v = o)jkr. Substitution of equations (5.3) into the coupled wave equations (5.2) and analysis of the resulting system yields solutions for the velocity v, wavenumber k and amplitude ratio U°/u° as functions of frequency*. The details of this process are not reproduced here. By solving the coupled wave equations (5.2), Biot (1956a) revealed that one shear wave, and two types of compressional wave can propagate in homogeneous poroelastic media. As discussed below, it is the compressional waves that are of most interest to us because they * Biot (1956a) and many other authors discuss solutions for v and k only; the amplitude ratio is not important to them. However, calculation of the amplitude ratio is quite straightforward (e.g. see Neev and Yeatts (1989)). (5.3) V = U°Qxp[j(kx-m)]x. Chapter 5: Modelling 108 are accompanied by relative fluid flow that generates charge separation and electric fields through electrokinetic coupling. Of the two compressional waves, the first, or fast, one is the only one normally observed. The fluid and solid motions associated with it are in-phase but generally not equal, so there is a small amount of relative fluid flow and attenuation. The second compressional wave is unique to porous media (i.e. it is not supported in purely elastic solids), and is known as the slow wave because of its very low velocity. The fluid and solid displacements associated with the slow wave are in opposite phase and yield significant relative flow. As a result the slow wave is heavily damped due to viscous dissipation. It attenuates by a factor of exp(-2ft) or approximately 500 for every wavelength travelled. The existence of the slow wave has been verified by laboratory experiments, but it has yet to be definitively observed in the field. Geertsma and Smit (1961) showed that slow waves should be generated by mode conversion whenever a fast compressional wave or a shear wave (at oblique incidence) is reflected from a boundary in porous media. However, because they attenuate so rapidly, and are not observed, slow waves are generally ignored in seismological studies. Given the large relative fluid flows associated with slow waves, they may play a more important role in seismoelectric phenomena. It is appropriate to note that basic Biot theory, as outlined above, does not account for all effects observed in experimental measurements of elastic wave attenuation and velocity. Modifications have been advanced by both Biot himself (1962a,b) and other researchers to account for anisotropy, and dissipation mechanisms other than viscous losses associated with bulk fluid flow. (Other significant loss mechanisms include inelasticity of the skeletal frame (friction between grains) and viscous losses due to local fluid motion near grain contacts (squirt flow)). Both can be accounted for by the use of complex, and frequency dependent bulk and shear moduli for the porous frame as discussed by Biot (1962a,b), Stoll (1989), and Dvorkin et al. (1995).) Recently, Russell (1996) has used a three component model (composed of solid, coupled fluid, and uncoupled fluid) to develop wave equations for Chapter 5: Modelling 109 poroelastic media that are very similar to those obtained by Biot (1956a). However, his analysis has lead him to question Biot's conclusion that the constants pn and Q are of opposite sign. The origin of these discrepancies has not yet been resolved. Despite the qualifications mentioned here, basic Biot theory is useful for modelling seismically-induced electrokinetic effects because it provides a first-order model for macroscopic fluid motion. Using Biot Theory to Generalize the Transport Equations governing Electrokinetic Coupling The second of Biot's equations of motion (5.1) can be re-written 4>(U, ~ K) = ^ | -0 ^ - (pnux + p22Ux) K 0771 where we have taken

VL = (v-w)y0£ + jv-jfdvEdt. (5. i i ) Chapter 5: Modelling 113 The total fluid volume within the element at any time is given by vf = vf-Vf V V 0 Vout where V0f is the volume of fluid in the element prior to the application of any load. Substituting for Vfut using equation (5.11), and dividing by VE yields yf yf cu.™\\\\/E + 6VE VE+8VE VEi V VQ i r Q and, neglecting second-order terms, this equation reduces to yf yf —=- = S- - V - w . (5.12) yE yE Equation (5.12) states that the volume of fluid, or fluid content, in any elementary volume attached to the solid frame fluctuates about an average value V0f /VE, and that the amplitude of the fluctuation depends on the divergence of the relative displacement w. Appropriately, Biot (1962a) called V - w the increment of fluid content. Since the mobile pore fluid carries excess counter-ions, we expect to find accumu-lations of the counter-ions where the fluid content is high and depletions where it is low. These regions will be separated by 1/2 of the seismic wavelength like the peaks and troughs of the compressional wave. As discussed in Chapter 1, the counter-ions are normally positive under most geological conditions, so regions of high and low fluid content should correspond to regions of excess positive and excess negative charge respectively. Of course, as charges begin to separate, electric fields develop which tend to pull them back together; the degree of charge separation is determined by the balance between streaming and conduction currents. Finally, as an aside, note that we can identify the term on the left hand side of equation (5.12) as the porosity 0(x,r) at a particular point and time, and the first term on the right hand side as the average porosity 0O. In the special case of a compressional wave propagating without relative fluid flow, w = 0, and the porosity does not fluctuate. Variations in porosity are normally neglected in studies of the acoustics of porous media because, like the density Chapter 5: Modelling 114 variations that accompany compressional wave propagation in pure elastic solids, they can generally be treated as second order effects. 5.3.2 Overview of the Conceptual Model As discussed above, the compressional wave produces regions of excess positive and negative charge parallel to the seismic wavefronts. A simple electrostatic model is now introduced to show that the electric fields generated by these charge separations are observed (i) when seismic wavefronts pass by a sensor, and (ii) when the spherical symmetry of the charged wavefronts is broken by a boundary or other heterogeneity. The existence of these two types of arrivals is also supported by more rigorous (but complicated) theoretical treatments (e.g. Haartsen and Pride, 1994; Haartsen, 1995). In the modelling that follows, I focus on the case where the electric field receivers, and a point source of compressional waves are located at surface of a (possibly layered) halfspace. This is the geometry that has been used most commonly for making seismoelectric measurements in the field. However, the same two types of seismoelectric arrivals are expected for any source-receiver geometry. The conceptual model is illustrated in Figure 5.1. First consider the case where compressional waves are generated at a point on the surface of a homogeneous halfspace as in Figure 5.1(a). The regions of excess positive and negative charge are approximated by charged spherical shells separated by one-half of the dominant seismic wavelength. The total positive charge on one shell exactly balances the negative charge on the other to preserve overall electroneutrality. For a boundary condition, we require that the normal component of the electric field be zero at the earth's surface. Then, as long as the charged shells remain in a homogeneous medium, it is straightforward to show that the electrostatic field remains confined to the region between them (details are given below). As a result, the contained electric field is not observed until the seismic wave reaches the dipole sensor. This provides a qualitative explanation for seismoelectric signals that travel at the same velocity as the seismic Chapter 5: Modelling 115 wave. Such arrivals dominated the data I acquired at Base Borden, and at sites on the Fraser Delta. The secondary arrivals observed at Haney, and the signals reported by Ivanov (1939, 1940) are also of this type. However, the model does not explain why the onset of these seismoelectric arrivals commonly appears to precede the seismic first breaks by a few milliseconds. The arrival time discrepancy can be explained if we consider the electric fields associated with a critically-refracted compressional wave travelling below a thin surficial layer such as the low-velocity weathered layer often encountered at the earth's surface. More details on that generalization to the homogeneous model are given in the next section (5.3.3). Electric fields can also be observed when the spherical symmetry of the charged regions is broken by an inhomogeneity such as a boundary (Figure 5.1(b)). In principle, the boundary could separate regions with differing permeabilities, elastic or electrokinetic properties (any property that affects the distribution of charge transported by the seismic wave). Here we consider the case of a perfectly reflecting boundary - a reasonable approximation for the case at Haney where loose fill overlies highly competent glacial till. As shown in Figure 5.1(b), the reflected charge distribution can be modelled as the superposition of two simpler distributions: (1) a pair of hemispherical shells in homogeneous media, and (2) a pair of spherical caps joined at the boundary. As already discussed, the field due to (1) is confined to the region between the shells. The electrostatic field generated by (2) however is everywhere non-zero and can be observed simultaneously by widely separated sensors. For the case where the electrical properties are uniform across the boundary, it can be calculated analytically as the sum of the fields from two charged caps as demonstrated below. The horizontal component of this field has radial symmetry at the surface, and for sensor offsets significantly greater than the interface depth, it decays as 1/x4. Furthermore, the arrival time of this signal is equal to the seismic traveltime from shotpoint to boundary. The electrical response generated at a boundary in this way therefore provides an explanation for the primary seismoelectric arrival observed during the Haney field experiments. 116 a) Figure 5.1. Conceptual modelling of two types of seismoelectric arrivals. Through electrokinetic coupling, regions of excess positive and negative charge form parallel to the peaks and troughs of a seismic P-wave. This charge separation produces electric fields which can be observed (a) when the P-wave passes by a receiver, and (b) when the spherical symmetry of the charged regions is broken by a boundary (here a perfect seismic reflector). Chapter 5: Modelling 117 5.3.3 Detailed Analysis of the Conceptual Model The Seismoelectric Effect in Homogeneous Media The characteristics of the conceptual model summarized above are now derived in a more formal manner. In particular, analytic equations are determined for the electric field patterns produced by the charged wavefronts shown in Figure 5.1. The model could be expressed in terms of streaming currents in conductive media rather than charge accumulations in dielectric media. However, I have chosen to focus on the electrostatic formulation because I find the concept of charge accumulation more appealing physically. We begin with the simplest case - the electric field due to compressional seismic waves emanating from a point source at the surface of a homogeneous halfspace. Figure 5.2 shows a diagram useful for the analysis of this problem. Because we are using an electrostatic model, the halfspace is treated as a dielectric with permittivity e but no conductivity. In order to obtain electric field (E) patterns like those that would exist in a conductive halfspace, we require that the normal component of E be zero at the surface. Given this boundary condition, and the spherical symmetry of the charge distribution, it is clear that the electric field in the halfspace must be purely radial. In order to determine the magnitude of E , we can now use Poisson's equation in integral form which relates the electric flux D = eE flowing through a closed surface S to the total enclosed charge Q as follows: e£(27cr 2) = <2 Chapter 5: Modelling 118 —. =*E(r) = Q/(e2Tir2). (5.14) (Note that 2 7 1 ^ is the area of the hemisphere S\".) Now if r < b, then the enclosed charge Q is zero. If r > a we also have Q = 0 because the total charges on the two shells are equal and opposite. However, for b< r < a , Q is equal to the charge on the inner shell Qi. Thus, from equation (5.14) we have E(r) = 0 , ra , ba where the constant coefficients are G = (cos o p ( c o s op 2n + l e is the permittivity q is the surface charge density, and P„(x) represents the Legendre polynomial of degree n evaluated at x. Chapter 5: Modelling 121 Note that the potential \\|/ is a function of r and 0 only, and is therefore symmetric about the z-axis like the charge distribution itself. The infinite series arises because Jeans (1948) represented the charge distribution by an infinite series of spherical harmonics. In practice approximate solutions for y( r ,6) are obtained by truncating the series after a finite number of terms. The number of terms required to obtain a good approximation is greatest at radii r ~ a where a is the radius of curvature of the cap. In the results presented below, the number of terms used to approximate the infinite series is called N. An expression for the electric field could be obtained by taking the gradient of equation (5.16) in spherical coordinates, but for purposes of modelling (and for purposes of comparison to potential differences measured in the field) I found it more convenient to calculate \\|T and then obtain E by computing the gradient numerically. Returning to the specific problem at hand, we need to calculate the potential due to a pair of spherical caps at a depth d as illustrated as the bottom of Figure 5.1(b). The first step is to use equation (5.16) to calculate \\|/z, and \\\\fu for the lower and upper caps individually and then apply the principle of superposition \\ |/(r,e)=\\|/ t(r,9)+ V l ,(r /,e /) . (5.17) The primed coordinates emphasize that the origin for the upper cap is different from that for the lower cap. Secondly, in order to satisfy the boundary condition of En = 0 at the earth's surface, we appeal to the method of images and place a pair of charged caps (image caps) a height d above the surface. The total potential at any point below the surface is then determined by superposition of the individual fields from four charged caps. Along the earth's surface, the calculations are easier because the effect of the image charges is simply to double the potentials that would be measured in that plane if the surface (and image charges) were not present. That is, for points at the earth's surface, the total electric potential is simply double that given by equation (5.17). ChapterS: Modelling 122 Figure 5.4 is a log-log plot of the electric field profiles generated at the surface by the reflection of a compressional wave from a boundary at 5 different depths according to the conceptual model. The field associated with the direct seismic wave passing by the sensors has not been included. The half-wavelength (or distance between positive and negative charged wavefronts) A/2 is taken as 5 m; in a medium such as saturated sand having a typical seismic velocity of 1500 m/s, this would correspond to a dominant frequency of 1500/(2x5) = 150 Hz. In order to ensure that the total charge Qi on the inner wavefront remains equal and opposite to that on the outer wavefront (thereby preserving electroneutrality in the halfspace), I assume that the charge density varies inversely as the wavefront radius squared - i.e. that q(b) = Q,/(2%b2) and q(a) = -Q{j (27ta2). In fact, this may overestimate the rate of decay of q with wavefront expansion because the far-field decay rate for particle displacements associated with a spherically spreading elastic wave decay is only 1/r. The profiles plotted in Figure 5.4 have been calculated for the instant in time when the outer wavefront's radius is given by a = d + A/4. Then, as shown in the accompanying diagram, the maximum height of the charge distribution at the boundary is A/2. The absolute electric field amplitudes depend on the value chosen for Qi and are therefore somewhat arbitrary, as this model provides no means for estimating the amount of charge separation. Here we have assumed Qi = 2 pC which is the approximate value required to obtain modelled amplitudes comparable to those measured in the field at Haney, B C (discussed further below). Figure 5.4 indicates that the peak electric field at the surface decreases by approximately one order of magnitude for every doubling in the depth to the reflector*. This rather severe depth dependence probably represents a worst case scenario because, as discussed above, I have assumed that the charge density q is strongly depth (radius) dependent. * In real conductive media a time-varying electric field would also be affected by conductivity-related attenuation which is not included here. However, at 100 Hz, the skin depth 8 = ^ j2/(cop(j) in a rather conductive medium with o = 100 mS/m, and p = |io, is 160 m which is large compared the depths and offsets of most interest here. 123 1.0E-03 • I I I I I I I I | I I I I I I I 11 l I I I I l l 11 : E increases as x : = = = = = = : : : = E decays as 1/x4 •Q 1.0E-06 • —^-1!— = = = = = = = = = ~ = mnrn^^ - - - - - - - - -~% 1 .OE-07 • = = = = d = z - ^ * — \" \" ^ j = ^ | = = ^ = = = 1.0E-09 1.0E-10-i 1 1 ' 11111) i i i 111 m i i i 11111j i i i 11111| 0.1 1 10 100 1000 Offset x (m) x o-I Figure 5.4. Log-log plots of the electric field profiles generated at the surface by the reflection of a compressional seismic wave from a boundary at 5 different depths. The modelling parameters were X/2 = 5m, a = d+KIA, <2, = 2 pC, and e = to • The number of terms N used to approximate the infinite series expansion for Mf was 100 for depths of 5,10, and 20m, and 200 for depths of 40 and 80 m. Note: At offsets x > 0, the field Ex at the surface is actually in the -x direction; the absolute value of the field is plotted here. Chapter 5: Modelling 124 Despite the uncertainty regarding absolute amplitudes, we can gain some useful information on the expected shape of the electric field profiles - i.e. on their amplitude vs. offset behaviour. Regardless of the boundary depth, the amplitude increases in proportion to the offset x at short distances, and decays rapidly as 1/x4 at long distances from the shotpoint. The position where the electric field is greatest - an optimum offset for making seismoelectric measurements -increases with the depth to the boundary. For depths much greater than the wavelength X, the optimum offset is xop = d/2. The potential due to a single cap, as expressed by equation (5.16) is a multipole expansion. One would expect that at sufficiently large distances from the charge distribution, the lowest order term in the expansion will dominate. For a pair of caps carrying equal and opposite charge, this should be the dipole term (the monopole term being zero because the net charge is zero). It is interesting therefore to compare the electric field of a vertical dipole to that produced by the spherical caps. The electric potential due to an electrostatic dipole at a distance r much greater than the dipole length is / „ s. mcos8 V * M ) = — — (5-18) where m is the dipole moment. If the dipole lies at a depth d below the surface as shown in 1/2 Figure 5.5, then cosG = d/r and r = [x2 + d2) for points along the surface and the horizontal component of the electric field is E =-^- = —L dx dx fmd^ V r 3 J d_ dx md (x2 + d2) 3/2 => £ a = _ * ~ f r ( 5 .19) (x2+d 2f Equation (5.19) shows that at short offsets [x2 «d2), Ex increases proportionally with x while at far offsets {x2 » d2) the field decays as 1/x4. To find the point where Ex reaches its maximum, we set dEx/dx = 0 and obtain the solution xop = d/2. Thus, the electric field Figure 5.5. Diagram of the vertical dipole used to approximate a pair of charged spherical caps. 1.0E-03 1.0E-10 -I 1 1 • • 1 1 0.1 1 10 100 1000 Offset x (m) Figure 5.6. Comparison between the electric fields generated at the surface by a pair of charged spherical caps and a vertical dipole at five different depths. The pair of caps is well approximated by the simple dipole when the depth to the boundary is large compared to the seismic wavelength (see text for details). Chapter 5: Modelling 126 profile due to the vertical dipole exhibits nearly identical characteristics to those described above for the spherical caps. In contrast, it can be shown that the field Ex due to a horizontal dipole decays as l /x 3 for long offsets and is maximum at x = 0. Figure 5.6 shows electric field profiles generated by vertical dipoles at various depths compared to those generated by pairs of caps. The dipole moment m used in the calculations is given by m = l{^){^\\ (5.20) where the factor of 2 accounts for the image dipole (required to satisfy the halfspace boundary condition En = 0), A/2 represents the effective dipole length, and the factor aq/2e is the same as the scaling factor that appears in equation (5.16). Recalling that a = d + A/4, and our assumption that q = -QJ(27ia2) on the outer wavefront, we see that the dipole moment is determined by three variables <2,, d, and A. The important point illustrated by Figure 5.6 is that for depths much greater than the seismic wavelength (d » A), the spherical cap charge distribution is very well approximated by a simple vertical dipole. The match between the two sources at shallow depths might be improved if we were to take into account the true length of the dipole rather than make the short dipole assumption inherent in equation (5.18). The seismoelectric amplitude vs. offset profile measured at Haney, and shown previously in Figure 4.7, provides a real data set for comparison to the spherical cap model. The amplitudes of the seismoelectric signal generated at the fill/till boundary and observed by several 2 m dipoles at the surface are plotted in Figure 5.7 along with the potential differences predicted by the conceptual model. The peak amplitude of the modelled profile was set approximately equal to that of the real data by choosing £); = 2 pC. The important features to compare are the basic shape and decay rates of the two profiles. At offsets greater than 4 m the measured amplitudes decay approximately as 1/r4 in agreement with the model. The maximum amplitude in the real and synthetic profiles occur across dipoles centred at Chapter 5: Modelling 127 approximately 4 m and 2 m respectively. Factors that could account for this discrepancy include possible in-line dip on the fill/till interface, variations in the near-surface resistivity (probably small), and the electric field trapped within the direct wave which is not modelled here. Figure 5.1(b) shows that this field is directed opposite to that from the pair of spherical caps. In this case, the charged wavefronts associated with the direct wave lie between offsets of 1 and 3 m and would therefore reduce the amplitudes measured by 2 m dipoles centred between 0 and 5 m offset. Thus the effect of the direct wave on the synthetic profile would be to shift its peak to a greater offset as observed in the real data. On the whole, the agreement between the decay rates and general shapes of the modelled and real profiles is encouraging. The way in which the Haney data differ most from that predicted by the conceptual model is in fact in the apparent absence of any strong electric field associated with the direct wave. One would expect such a seismoelectric signal to sweep past the dipole sensors at a velocity equal to that of the compressional wave in the fill (a very slow 250 m/s). According to equation (5.15) the electric fields associated with that wave should be approximately one order of magnitude greater than that produced by reflection at the fill/till boundary. The absence of any such strong arrival moving slowly across the dipole array (see Figure 4.6) is puzzling. It is possible that the road fill is so highly attenuative of seismic waves that the direct wave and its electric field diminish rapidly with offset. In any case we must acknowledge that the model presented here provides a partial but not perfect explanation for the seismoelectric effect at Haney. 1.4 0.0001 4 • i ' •— i • •— | 0.1 1 10 100 Offset x (m) Figure 5.7. Comparison of measured seismoelectric amplitudes to those predicted by the charged spherical cap model. The real data are potential differences generated by the arrival of a seismic wave at a fill/till boundary at 2 m depth. They were measured by eleven 2 m dipoles on the surface, centred at offsets of 3 to 15 m from the shotpoint. The parameters used for modelling were d= 2m, a = 3m, A/2 = 2m, <2, = 2 pC, e = to , and N = 7 terms. Chapter 5: Modelling 129 x a •> Figure 5.8. Charge distribution induced by a critically refracted compressional wave. The Seismoelectric Effect associated with a Refraction A third scenario that is amenable to investigation by the concept of charged wavefronts is the case of a critically refracted wave travelling below a thin low velocity surface layer. This case is of practical interest because a layer of low velocity sediments or weathered rock is commonly found at the surface of the earth. For point source at the surface, and a layer of uniform thickness d, the refracted wavefront geometry beyond a certain critical offset x' * is as shown in Figure 5.8. Now the charge densities q\\ and q2 induced above and below the boundary will generally be different due to differences in the zeta potential, conductivities, permeabilities, porosities, and other factors affecting relative fluid flow. The abrupt change in charge density and the lack of spherical symmetry give us reason to doubt that the electric field remains contained between the wavefronts as it does in the case of homogeneous media. We can again employ the model of charged spherical caps to calculate electric field profiles if we make two simplifying assumptions. First we assume that the charge densities induced in the upper layer are negligible. This could be the case if for example, the upper layer was * The critical distance is the offset at which a refracted wave first returns to the surface (e.g. see Telford et al., p. 277). ChapterS: Modelling 130 completely dry or a layer of water such as a lake. Secondly, we assume that the wavefronts in the lower medium can be approximated by hemispherical shells centred at point O on the interface directly below the shotpoint. The validity of this approximation improves with time as the refracted wavefronts travel farther and farther through the second medium. The traveltime required to justify this approximation is small when the layer is thin and the velocity contrast v 2/vj is high. With these two assumptions, the charge distribution reduces to two hemispherical shells centred at point O (see Figure 5.9), and a fringing electric field analagotis to that observed at the edges of a parallel plate capacitor must exist in the surficial layer. If the electrical properties (dielectric constants) of the layer and underlying media are the same then the electric field profile Ex— — dy//dx at the surface can be calculated by superposing the potentials from two charged hemispheres (plus two image hemispheres to satisfy the boundary condition En = 0), and then calculating the derivative of the total potential field numerically as described previously. Equation (5.16) simplifies slighdy because a = 90° for both hemispheres and hence 7Jn(cosa) = Pn(0) = 0 for odd n. Synthetic profiles for Ex calculated in this way using the parameters a = 20 m, X/2 = 5 m, e = e0, and <2< = 2 pC are presented in Figure 5.9 for layer thicknesses between 0 and 5m (0 to 1/2 wavelength). In the case of d = 0, the model reduces to that of the homogeneous halfspace and the electric field is completely confined between the charged wavefronts in agreement with equation (5.15). As the layer thickness increases, a fringing electric field develops and becomes detectable farther and farther ahead of the seismic wave. As a result, the electric field associated with the refracted wave is observed at the surface before the seismic arrival. However, the electric field observed at a given sensor still obtains its maximum value when the charged wavefronts pass directly below it, so seismoelectric arrivals still appear to travel at the seismic wave velocity. This model is useful because it provides a simple explanation for why seismoelectric signals observed in the field often seem to travel at seismic velocities but precede any seismic arrivals by a few milliseconds. Figure 5.9. Electric field profiles (top) generated by critically refracted charged wavefronts travelling below a thin surficial layer. Profiles are shown for values of d ranging from 0 to 5 m (zero to one-half wavelength). Chapter5: Modelling 132 Finally, I point out that the seismoelectric effect due to a refracted wave is essentially a hybrid of the two end-member cases presented in Figure 5.1. In the limit of d = 0, this case reduces to that of the direct wave in homogeneous media. In addition it helps to explain how the seismoelectric effect due to an interface can evolve long after the initial seismic reflection. 5.4 QUANTITATIVE M O D E L S In the preceding section I have shown how a conceptual model can be used to explain the origins of seismically-induced electrokinetic effects and estimate their relative ampltudes. However, a more quantitative model is required to estimate absolute amplitudes. Progress made by others in that direction is summarized below. I pay particular attention to the theory of Neev and Yeatts (1989) showing how it links Biot's wave equations to the transport equations governing streaming potentials which were introduced in Chapter 1. Finally, I calculate order of magnitude estimates for the amplitudes of seismically-induced electrokinetic effects and compare them to the measurements presented in Chapter 4. 5.4.1 Frenkel's Suggestion: Use of the Helmholtz-Smoluchowki Equation The first significant attempt to model seismically-induced electrokinetic effects in a quantitative fashion was made by Frenkel (1944). He suggested that the effect be estimated by calculating the relative fluid flow and using a modified form of the Helmholtz-Smoluchowski equation. The Helmholtz -Smoluchowski equation was presented in Chapter 1 (as equation 1.4); it is the simplest equation relating the streaming potential difference A\\|/ to the pressure difference Ap driving fluid flow across a porous sample. Writing the equation in terms of potential gradients rather than differences gives V \\ | / = - ^ - V p . (5.21) T\\Of For Poiseuille flow through a capillary tube of radius r, the average velocity of the flow is given in terms of the pressure gradient by Chapter 5: Modelling 133 r2 = - — Vp . 87] (5.22) Frenkel assumed that this formula could be applied to flow through porous media by use of an effective pore radius r, and substituted for Vp in equation (5.21) to obtain In the case of fluid flow induced by seismic waves, t> represents the velocity of fluid relative to the porous frame, i.e. v> = -^(U - u). Frenkel solved for the relative velocity induced by compressional waves by developing wave equations similar to Biot's. He substituted his result for \"0 into an equation like (5.23) to obtain an equation for the seismoelectric field. Further work in this direction was carried out by Tome (1975) in an M.Sc. thesis. He used Biot theory to find an expression for the relative velocity, and calculated numerical estimates of the seismoelectric signal strength expected in sandy sediments over a broad range of frequencies. He also provided a summary of the differences between Frenkel's and Biot's derivations of the poroelastic wave equations, and pointed out an inconsistency in Frenkel's approach. Using equation (5.23) to calculate seismically-induced electrokinetic effects is attractive because of its conceptual simplicity. However, two limitations should be noted: (i) It is assumed that the induced electric fields have no effect on the relative fluid flow. This seems reasonable since that effect is second order in the electrokinetic coupling. (ii) The Helmholtz-Smoluchowski equation was developed under the assumption that the net electrical current is zero everywhere in the medium (i.e. J\" = + J', = 0 as discussed in Chapter 1). But, in the case of time-varying flows and currents, such as those induced by a seismic wave, the correct statement of charge conservation is Since time differentiation corresponds to multiplication by y'co in the frequency domain, the (5.23) (5.24) Chapter 5: Modelling 134 error associated with the assumption that J* = 0 is expected to increase with frequency. A theory that overcomes the two limitations listed above, and draws on the general principles of coupled flows discussed in Chapter 1, is presented next. 5.4.2 A quasistatic theory for electrokinetic coupling in the elastic wave equations Neev and Yeatts (1989) proposed to incorporate electrokinetic effects in Biot's theory for elastic wave propagation by adding terms to the equations of motion (5.2) to account for the electric forces exerted on the macroscopic charge densities qj and q, in the solid and fluid phases respectively. Using the quasistatic (low frequency) approximation, these forces are qfVy/ and q^y and the equations of motion become NV2u + V[(A + 7V)V. u + QV • U] = 4r ( P n » + P12U) + — ( u - U) + qsVy/ at K dt (5.25) V[£V. u + * V • U] •= | l (p 1 2 u + p 2 2U) - ^ f - ( u - U) + qfV ¥ . dt K dt In order to solve for the additional field \\|/, a seventh equation is required in addition to the six represented by (5.25). Also, equations are required relating qj and qs to material properties of the porous medium. These requirements were met by introducing the transport equations for coupled flows of fluid and electrical current y =-(K/TI)V/? - Z V y (5.26) Je=-XVp - oVy. (5.27) Note that these are the transport relations discussed in Chapter 1 for the case of steady flows; in using them, Neev and Yeatts ignored the inertial terms which should be added to the pressure gradient as shown in equations (5.5) and (5.6) for wave-induced flow. This does not affect their solution for the electric charge densities because, as outlined below, qj and qs were found by considering a limiting case where the inertial terms are zero. However, it does affect the form of their seventh equation which is based on equation (5.27). I will ignore this problem for the time being and proceed with an overview of the rest of Neev and Yeatts' Chapter 5: Modelling 135 derivation. Neev and Yeatts observed that in the limiting case of steady flow through a solid at rest (i.e. U = u = ii = 0), the second of equations (5.25) reduces to the transport equation for fluid flow (5.26) provided that qf = (T\\§/K)X . Substituting forX using equation (1.18) yields (J)2ez af = \\ with z = c-J 1 e, v x 2 e j (5.28) (5.29) where e is the bulk permittivity of the saturated porous medium ef is the permittivity of the pore fluid and £ , x are the zeta potential and pore tortuosity introduced in Chapter 1. Equation (5.28) shows that the charge density in the fluid is proportional to the zeta potential. For the solid's charge density, it is assumed that qs ~ -qf (and hence that the total charge density q = qf +qs ~0 everywhere). This assumption should be reasonable for purposes of calculating the electric forces <7JV\\|/ and ), X((o), and K ( W ) are well approximated by their dc values. However, it is not clear that the modified driving force (-Vp+co2p /u) will be well approximated by simply -Vp. In the case of fast Chapter 5: Modelling 139 compressional waves for which the solid and fluid motions are in phase, this modification to the direct driving force for appears to account for the fact that part of the pressure gradient is required simply to accelerate the fluid at the same rate as the solid (note that p / i i = co2p /u for the assumed harmonic time dependence); only a portion of the pressure gradient then remains to drive relative flow. Because Neev and Yeatts (1989) did not include this inertial term in the transport equations, their theory may overestimate the size of the electric field generated by electrokinetic coupling. Unfortunately, a comparison of the electrokinetic fields predicted by these two different models is not available as Pride (1994) was wholly concerned with deriving the governing equations and did not present their plane wave solutions. Fully Coupled Seismoelectric waves in Layered Media Each of the quantitative theories described to this point has been concerned with the case of waves travelling through unbounded homogeneous media. In his Ph.D. thesis, Haartsen (1995) built upon the work of Pride (1994) and modelled seismoelectric effects that arise as seismic waves cross boundaries in porous media. The issue of boundary conditions was investigated in detail, and numerical techniques were developed to account for seismic point sources and for a variety of source-receiver geometries of interest in exploration geophysics. Haartsen's results confirm the existence of the two types of seismoelectric effects predicted by the conceptual model (Figure 5.1). His model also indicates that the electromagnetic effects associated with shear waves are primarily magnetic fields. As one test of his modelling theory and algorithms, Haartsen (1995, Chapter 7) calculated synthetic seismic and seismoelectric records for a layered model representing the Haney test site; he then compared his results to records that I had measured there. One of the data sets modelled is the seismoelectric response vs. offset data shown in Figure 4.6. The synthetic records confirm that a seismoelectric conversion is produced at the fill/till boundary, and that the earliest of the secondary arrivals is probably generated by a refracted P-wave Chapter 5: Modelling 140 travelling along the fill/till interface. In order to match the measured amplitude of the primary seismoelectric response, the model required a seismic source equivalent to a 0.25 kg charge of dynamite. This is clearly a stronger source than the sledgehammer impact actually used in the field. The discrepancy may be partly explained by the fact that the boundary was placed at a depth of 3 m in the model while its true depth at this location was only 2 m. Perhaps the most significant difference between the real and synthetic records was the lack of evidence (in the real data) of strong seismoelectric signals associated with compressional waves arriving at the dipoles. The lack of such signals was also pointed out in section 5.3.3 above where I compared this data set to the predictions of my conceptual model. Both Haartsen and I have speculated that these signals may be smaller than expected because P-waves attenuate rapidly in the highly compressible road fill. In any case, the road fill probably differs significandy from the ideal poroeleastic materials handled by Haartsen's model. 5.4.4 Order of Magnitude Estimates for Electrokinetic Signal Strength Order of magnitude estimates for the seismically induced electrokinetic effects contained within a seismic wave can be calculated from the Neev and Yeatts (1989) model outlined in Section 5.4.2. The ratio of the induced electric field to the seismic particle velocity for a compressional fast wave is given in terms of fundamental material properties by equation (5.33). Neev and Yeatts provide estimates of these material properties for a permeable, water-saturated, medium-grained quartz sandstone*. In Appendix D, I provide estimates of the same parameters for the water-saturated sands at the MOT Tower Site on the Fraser Delta in Richmond, BC. Substituting these two sets of material properties into equation (5.33) yields the following ratios: * I have corrected errors in the values given by Neev and Yeatts for Biot's elastic constants Q, R, P, and H. Their values give an impossibly low value of 1400 m/s for the fast wave's characteristic velocity vc while the corrected values give a much more reasonable 3570 m/s. The correction makes the ratio Elu for sandstone a factor of 2.5 times smaller than that obtained using Neev and Yeatts' values. Chapter 5: Modelling 141 (i) for the sandstone, E (co/o) 1.9x10- -7 V / m (5.35) u m/s (ii) for the sand, E (co/a) 8.8x10' ^ V / m (5.36) Because conductivity can vary widely, it is left as a variable in the above expressions. Note that the expected electric field in sands, is a factor of almost 50 greater than that in sandstones. This is partly due to the fact that the factor £ / x 2 used for the sands is 5 times larger than that used by Neev and Yeatts for the sandstone. Ratio (5.36) can now be used to predict seismoelectric amplitudes for comparison to those measured in sandy sediments on the Fraser Delta and at Base Borden. An example of data acquired on the Fraser Delta using the 8-gauge shotgun source is shown in Figure 4.14. The conductivity in the sands at this site is approximately 10 mS/m, and the first seismic/seismoelectric arrivals have a dominant frequency of about 30 Hz. At offsets of about 20 m from the shotpoint the total seismic particle velocity (vector sum of the components measured by the vertical and horizontal geophones) is about lxlO\" 3 m/s and the electric field amplitude is approximately 6 u.V/m. Substituting for o~, CO, and u in equation (5.36) gives an estimate of 170 p:V/m which is a factor of 30 greater than the measured value. Data from the Base Borden site are presented in Figure 4.11. A representitive conductivity for this site is about 6 mS/m (the average of conductivities above and below the water table). The first seismic and seismoelectric arrivals at offsets of 5-8 m from the sledgehammer source have a dominant frequency of roughly 100 Hz and peak to peak amplitudes of about 3xl0\" 3 m/s and 25 u.V/m respectively. The electric field predicted by equation (5.36) is 2.8 mV/m or 110 times higher than the measured value. There are a number of reasons that may explain why the theoretical estimates exceed the measured values in both cases above: (i) The material properties may not be correct; the values of C, and x2 are most uncertain -Chapter 5: Modelling 142 both could be in error by factors of two or three. (ii) The seismoelectric signals may well be associated with a refracted wave in which case the measured fields could be decreased by a factor of three or more depending on the depth to the refracting layer (e.g. water table) as illustrated in Figure 5.9. Furthermore, particle velocities measured by geophones on the surface may not give an accurate measure of the particle velocities inducing the field at depth. (iii) The conductivity assumed for the Fraser Delta site may be too low given that the main sand body is overlain by 2-3 m of silty sediments having a conductivity of about 30 mS/m -three times higher than that of the sands. (iv) As discussed in section 5.4.4,1 have speculated that the transport equations used by Neev and Yeatts (1989) may overestimate the electrokinetic amplitude. Given these sources of uncertainty, the discrepancy between measured and predicted values is not unreasonable. Tome (1975, p. 82) also calculated order of magnitude estimates for the case of a particular water-saturated sand. The assumed conductivity is not explicitly stated. Ffis results indicate that for a dominant frequency and particle velocity comparable to those mentioned above for the Fraser Delta experiment, the expected electric field is of the order of 30 piV/m. This is only a factor of about 5 times larger than the measured field. At the Haney test site, the measured seismoelectric amplitudes range from about 1.3 mV/m peak-to-peak for the primary arrival at near offsets to 10 \\i\\/m peak-to-peak for the secondary arrivals at offsets of about 20 m. The large primary response comes from the fill/till interface and therefore should not be compared directly to the theoretical estimates given above for fields associated with waves in homogeneous media. Furthermore, the highly compressible, partially saturated road fill and underlying glacial till are both rather different from the sandy sediments on which those estimates are based. However, there are two reasons why we might expect relatively strong seismoelectric signals at Haney. First, the conductivities Chapter 5: Modelling 143 of both the road fill, and the underlying till are quite low (approximately 0.4 mS/m in the fill and comparable in the till). Secondly, the seismic particle velocities at the fill/till interface are large because of its proximity to the sledgehammer source. Measurements with a buried geophone for example indicated that typical particle velocities produced at the boundary by a sledgehammer blow on the surface were approximately 1 to 2xl0\"2 m/s - considerably larger than those considered above. Chapter 6 Conclusions Claims regarding the utility, and even the existence of, seismoelectric effects have been treated skeptically since they were first reported in the 1930's. The results contained in this thesis show that, despite the presence of strong ambient noise, seismoelectric effects can be measured in the field, and can be used as a tool for the mapping of shallow boundaries. The well-controlled field experiments carried out at Haney, BC, are not the first to yield promising results, but they are among the most convincing. The primary seismoelectric response at that site was generated as compressional waves crossed a lithologic boundary between permeable road fill and impermeable glacial till. It arrived simultaneously at widely separated dipoles, and was readily exploited to map the fill/till interface from the surface using a low-powered (sledgehammer) seismic source. As the shot-dipole offset increased, the amplitude of the primary response decayed rapidly, and secondary arrivals became evident. They differed from the primary response in that their arrival times increased with dipole offset, and they appeared to be associated with the arrival of seismic waves beneath each dipole sensor. Studies of the site geology, and of the polarity, arrival times, and linearity of the seismoelectric effect at Haney have revealed that it cannot be attributed to piezoelectricity, to the modulation of telluric currents flowing across the survey area, or to non-linear phenomena. The possibility that a vertical electric field, generated by natural oxidation/reduction potentials, exists at the fill/till interface and is modulated by the seismic wave deserves further investigation but is considered highly speculative. The favoured explanation for the seismoelectric conversion is an electrokinetic mechanism. Seismoelectric experiments carried out at sites on the Fraser Delta and at Base Borden, Ontario yielded electrical responses that were similar to the secondary signals at Haney. However, attempts to measure seismoelectric effects from interfaces including the water table, 144 Chapter 6: Conclusions 145 sand/silt, and sand/clay contacts at those sites were not successful. The fill/till interface at Haney was clearly a better target for seismoelectric detection. This may be explained in part by the lower electrical conductivities at Haney, and by the large particle motions we were able to induce at the fill/till interface because of its shallow depth. The interface itself differed from those investigated elsewhere in that it represented not only a large change in permeability, but also a huge contrast in elastic properties. A l l of these differences would be expected to contribute to an unusually large electrokinetic conversion. It is evident that stronger seismic sources, or noise levels lower than 1 pV/m will be required for seismoelectric mapping of interfaces at the Base Borden and Fraser Delta sites. The difficulty of resolving small signals in the presence of strong electrical noise has probably delayed the development of seismoelectric methods. As a member of the U B C seismoelectric research group, the author has contributed to the development of equipment and techniques that can be used with confidence to induce and measure seismoelectric conversions in the earth. Chapters 2 and 3 of this thesis detailed the steps taken to meet challenges in the design of reliable, low-noise instrumentation, in the development of field methods free from acquisition artifacts, and in the elimination of ambient or environmental noise. Powerline harmonics were the strongest type of interference in field records acquired for this study. Noise was suppressed primarily by the subtraction of estimates which were either (i) recorded by a remote dipole, (ii) recorded immediately before the shot, or (iii) calculated by means of a least-squares method of harmonic decomposition. The remote dipole subtraction technique reduced interference from both power lines and spherics, and typically lowered noise levels by a factor of 10 to 50. I suspect that its effectiveness could be enhanced significantly if the remote record were to be match filtered or adaptively filtered to match the noise in the other record as closely as possible prior to subtraction. Block subtraction and sinusoid subtraction constitute two novel techniques for the removal of powerline harmonics from a time series. Unlike notch filters, they are capable of removing multiple harmonics Chapter 6: Conclusions 146 without attenuating or distorting the signal of interest. The sinusoid subtraction technique was used to process all of the seismoelectric data presented in this thesis and yielded signal-to-noise improvements of up to 45 dB (a factor of about 180). At the time this thesis was started, quantitative models for seismically-induced electrokinetic effects had been published only for the case of compressional waves in homogeneous media (Frenkel, 1944; Neev and Yeatts, 1989). There was a need to determine if electrokinetic responses could also originate at boundaries in layered media as suggested by our experimental results at Haney, and by other researchers (e.g. Martner and Sparks, 1959; Migunov, 1987; Thompson, 1990). The conceptual model developed in Chapter 5 predicts that two fundamental types of electrokinetic effects are generated as a result of the charge separation induced by compressional waves in porous media. The first of these is an electric , field that is contained within the seismic pulse and travels at the seismic wave velocity. The second effect is generated when the distribution of charge transported by the seismic wave is altered as the wave crosses a boundary or other inhomogeneity; the resulting electric fields propagate at the electromagnetic wave velocity and are therefore observed essentially simultaneously by widely separated antennas. In principle the boundary may separate regions with different elastic properties, electrical conductivities, permeabilities, or pore fluids since all of these factors influence the magnitude of the charge accumulations that form parallel to the seismic wavefronts. It remains to be seen which of these properties is most important in reality. The potential for detection of high permeability formations and pore fluid contacts is of particular interest as these targets are not easily detected by conventional geophysical techniques. The two types of arrivals predicted by the conceptual model are consistent with those observed during field experiments. Signals of the first type, which travel at seismic wave velocities, were observed during all field experiments. A signal of the second type, originating at a shallow interface and arriving simultaneously at widely separated sensors, was observed at Chapter 6: Conclusions 147 Haney. It exhibited the same radial symmetry, and rate of decay with offset, as the vertical dipole-like field produced by the model. According to the quantitative theory of Neev and Yeatts (1989), the magnitude of the electrokinetic field contained within a fast compressional wave propagating through homogeneous media at low frequencies is given by equation (5.33). The field is proportional to frequency and seismic particle velocity, and to the ratio of several physical properties, some of which also govern the magnitude of streaming potentials generated by steady flows through porous plugs (compare equations (1.6) and (5.33)). Order of magnitude estimates for the electrokinetic fields expected at the Base Borden and Fraser Delta sites were calculated using this result. To my knowledge, this represents the first attempt to test the theory of Neev and Yeatts against experimental data. The predicted amplitudes exceed those measured in the field at Base Borden, and on the Fraser Delta by one to two orders of magnitude. The discrepancy is not unreasonable considering the experimental uncertainties, and uncertainty about the effect of an approximation in the model's derivation. The principal conclusion of this thesis is that seismoelectric effects can be measured and used as a tool to map shallow lithologic boundaries in sedimentary materials. Seismically-induced electrokinetic effects are considered to be the most likely mechanism for the seismoelectric responses presented here. The experimental results from Haney support recent theoretical developments which predict that seismoelectric effects may be used to map the boundaries of permeable formations. However, further research is required to identify the types of boundaries and physical properties most amenable to detection. References Adams, R . K . , Mclntyre, J . M . , and Symonds, F .W. 1982. Characteristics of the eastern interconnection line frequency. IEEE Transactions on Power Apparatus and Systems, PAS-101, 4542-4547. Berryman, J .G. 1980. Long wavelength propagation in composite elastic media, I. spherical inclusions. Journal of the Acoustical Society of America, 68, 1809-1819. Biot, M . A . , 1956a. Theory of propagation of elastic waves in a fluid-saturated porous solid. 1. Low frequency range. Journal of the Acoustical Society of America, 28, 168-178. Biot, M . A . 1956b. Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range. Journal of the Acoustical Society of America, 28, 179-191. Biot, M . A . 1962a. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics, 33, 1482-1498. Biot, M . A . 1962b. Generalized theory of acoustic propagation in porous dissipative media. Journal of the Acoustical Society of America, 34, 1254-1264. Biot, M . A . , and Willis, D . G . 1957. The elastic coefficients of the theory of consolidation. Journal of Applied Mechanics, 24, 594-601. Blake, E . W . 1992. The deforming bed beneath a surge-type glacier: measurement of mechanical and electrical properties. Ph.D. Thesis, Dept. of Geophysics and Astronomy, The University of British Columbia. Blau, L . W . and Statham, L . 1936. Method and apparatus for seismic-electric prospecting. US Patent No. 2 054 067, September 15, 1936. Bourbie, T. , Coussy, O., and Zinszner, B . 1987. Acoustics of Porous Media. Gulf PubHshing Company, Houston. Broding, R . A . , Buchanan, S.D. and Hearn, D.P. 1963, Field experiments on the electroseismic effect. IEEE Transactions on Geoscience Electronics, GE-1 , 23-31. Butler, K . E . 1991. Enhancement and modelling of signals from natural piezoelectric targets. M.Sc. thesis, Dept. of Geophysics and Astronomy, The University of British Columbia. Butler, K . E . , and Russell, R .D. 1993. Subtraction of powerline harmonics from geophysical records. Geophysics, 58, 898-903. Butler, K . E . , Russell, R.D. , Kepic, A . W . , and Maxwell, M . 1994. Mapping of a stratigraphic boundary by its seismoelectric response. Proc. of the 7th Symposium on the Application of Geophysics to Engineering and Environmental Problems (SAGEEP) , 689-699. 148 References Butler, K .E . , Russell, R.D., Kepic, A.W., and Maxwell, M . 1996. Measurement of the seismoelectric response from a shallow boundary. Accepted for publication in Geophysics in 1996. Clague, J.J., Lutemauer, J.L., Pullan, S.E., and Hunter, J.A., 1991. Postglacial deltaic sediments, Fraser River delta, British Columbia. Canadian Journal of Earth Sciences, 28, 1386-1393. Conant, B. Jr. 1991, Field studies of vertical well-screen placement: effects on recovery of stratified contaminants from an unconfined aquifer. M.Sc. Thesis. Dept. of Earth Sciences, University of Waterloo. Dallimore, S.R., Edwardson, K . A . , Hunter, J.A., Clague, J.J., and Luternauer, J.L. 1995. Composite geotechnical logs for two deep boreholes in the Fraser River delta, British Columbia. Geological Survey of Canada Open File Report 3018. de Groot, S.R., and Mazur, P. 1962. Non-equilibrium thermodynamics. North-Holland, Amsterdam. Dvorkin, J., Mavko, G., and Nur, A. 1995. Squirt flow in fully saturated rocks. Geophysics, 60,97-107. Fleming, S.W. 1994. Characterization of the relationship between seismic input and seismoelectric response at the Haney, B C test site. Report on undergraduate research project. Dept. of Geophysics and Astronomy, The University of British Columbia. Frenkel, J. 1944. On the theory of seismic and seismoelectric phenomena in moist soil. Journal of Physics (USSR), 8, 230-241. Freeze, R.A., and Cherry, J.A. 1979. Groundwater. Prentice-Hall, Englewood Cliffs, N . J. Gaskarov, I.V., and Parkhomenko, E.I. 1974. The seismoelectric effect in rocks and the preconditions for its application in geological prospecting work. Izvestiya, Academy of Sciences, USSR, Physics of the Solid Earth (English translation by A.G.U.), no. 1, 71-74. Geertsma, J. and Smit, D.C. 1961. Some aspects of elastic wave propagation in fluid-saturated porous solids. Geophysics, 26, 169-181. Gilson, E., and Redman, D. 1995. Field Evaluation of the PulseEkko Ground Penetrating Radar (GPR) Borehole Antennae Systems. Progress Report for Terrain Sciences Division ofGSC,Feb. 1995. Haartsen, M.W. 1995. Coupled Electromagnetic and Acoustic Wavefield Modeling inPoro-Elastic Media and its applications in Geophysical Exploration. Ph.D. thesis, Massachusetts Institute and Technology. References Haartsen, M . , and Pride, S. 1994. Modeling of coupled electroseismic wave propagation in layered media. 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1155-1158. Hancke, G.P. 1990. The optimal frequency estimation of a noisy sinusoidal signal. IEEE Transactions on Instrumentation and Measurement, 39, 843-846. Hassanzadeh, S. 1991. Acoustic modeling in fluid-saturated porous media. Geophysics, 56, 424-435. Horowitz, P., and Hi l l , W. 1980. The Art of Electronics. Cambridge University Press, Cambridge. Hunter, J.A., Lutemauer, J.L., Roberts, M.C. , Monahan, P.A., and Douma, M . 1994. Borehole geophysics logs, Fraser River Delta. Geological Survey of Canada Open File Report 2841. Hunter, R.J. 1987. Foundations of Colloid Science, Vol. I, Oxford University Press, New York. Hunter, R.J. 1993. Introduction to Modern Colloid Science. Oxford University Press, New York. Ishido, T. and Mizutani, H . 1981. Experimental and theoretical basis of electrokinetic phenomena in rock-water systems and its application to geophysics. Journal of Geophysical Research 86, 1763-1775. Israelachvili, J. 1991. Intermolecular and surface forces. Academic Press, San Diego. Ivanov, A . G . 1939. Effect of electrization of earth layers by elastic waves passing through them. Comptes Rendus (Doklady) de VAcademie des Sciences de I'URSS, 24, No. 1, 42-45, (In English). Ivanov, A . G . 1940. The seismoelectric effect of the second kind. Izvestiya Akademii Nauk SSSR, seriya geograficheskaya i geofizicheskaya. 4, No. 5, p. 699-727, (English translation obtained from US Library of Congress). Ivanov, A . G . 1949. On the seismo-electric effect of the first kind (J) in the region near the electrode. Dokl. Akad. Nauk USSR, 68, No. 1, 53-56, (English translation from National Technical Information Service, U.S. Dept. of Commerce, document TT 6211280). Jeans, J.J. 1948, The Mathematical Theory of Electricity and Magnetism (5th edition). Cambridge University Press, Cambridge. Jou, D., Casas-Vazquez, J., and Lebon, G. 1993. Extended Irreversible Thermodynamics. Springer-Verlag, Berlin. References Kaufman, A .A . , and Keller, G.V. 1981. The Magnetotelluric Sounding Method. Elsevier, Amsterdam Kepic, A.W. 1995. Seismoelectric responses from sulphide orebodies. PhD Thesis, Dept. of Geophysics and Astronomy, University of British Columbia. Kepic, A.W., and Russell, R.D., 1996, Fiber optic time break: Geophysics, 61, 294-298. Kepic, A.W., Maxwell, M . , and Russell, R.D., 1995. Field trials of a seismoelectric method for detecting massive sulfides. Geophysics, 60, 365-373. Kondrashev, S.N. 1980. The Piezoelectric Method of Exploration. Nedra, Moscow (in Russian). An English translation is available at the University of British Columbia. Long, L.T. and Rivers, W.K. 1975. Field measurement of the electroseismic response. Geophysics 40, 233-245. Luternauer, J.L., Harris, J.B., Hunter, J.A., and Finn, W.D.L. 1995. Fraser River delta geology: modeling for seismic hazard assessment. The BC Professional Engineer, 46, No.9, 11-14. Martner, S.T. and Sparks, N.R. 1959. The electroseismic effect. Geophysics, 24, 297-308. Maxwell, M , Russell, R.D., Kepic, A.W., and Butler, K . E . 1992. Electromagnetic responses from seismically excited targets B: non-piezoelectric phenomena. Exploration Geophysics, 23, 201-208. Migunov, N.I. 1987. Electroseismic effect of ore bodies. Izvestiya Academy Sciences USSR, Physics of the Solid Earth (English translation by A.G.U.), 23, 946-951. Mitchell, J.K. 1993. Fundamentals of Soil Behavior. John Wiley and Sons, Toronto. Morgan, F.D., Williams, E.R., and Madden, T.R. 1989. Streaming potential properties of Westerly granite with applications. Journal of Geophysical Research 94, 12449-12461. Neev, J., and Yeatts, F.R. 1989. Electrokinetic effects in fluid-saturated poroelastic media. Physical Review B, 40, 9135-9141. Neishtadt, N . M . , Mazanova, Z.V., Binevich, L .Ya . and Maiko, M.I. 1972. Piezoelectric method of exploration (methodological recommendations). ONTIVITR, Leningrad. Nyman, D.C. and Gaiser, J.E. 1983. Adaptive rejection of high-line contamination. 53 r a\" Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 321-323. Onsager, L . 1931a, Reciprocal relations in irreversible processes I. Physical Review, 2nd Series, 37, 404-426. References Onsager, L . 1931b, Reciprocal relations in irreversible processes II. Physical Review, 2 Series, 38, 2265-2279. Overbeek, J.Th. 1952 Colloid Science, Vol. I, Irreversible Systems, edited by H.R. Kruyt, Elsevier, New York. Parkhomenko, E.I. 1971. Electrification Phenomena in Rocks, (translation by G.V. Keller). Plenum, New York. Parkhomenko, E.I., and Gaskarov, I.V. 1971. Borehole and laboratory studies of the seismoelectric effect of the second kind in rocks. Izvestiya Academy Sciences USSR, Physics of the Solid Earth (English translation by Am. Geophys. Union), 7, 663-666. Petiau, G., and Dupuis, A. 1980. Noise, temperature coefficient, and long time stability of electrodes for telluric observations. Geophysical Prospecting 28, 792-804. Press, H.P., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T. 1989. Numerical Recipes (Fortran version). Cambridge University Press, New York. Pride, S.R. 1994. Governing equations for the coupled electromagnetics and acoustics of porous media. Physical Review B, 50, 15678-15696. Pride, S.R., and Morgan, F.D. 1991. Electrokinetic dissipation induced by seismic waves. Geophysics 56, 914-925. Pullan S.E., and MacAuley, H.A. 1987. An in-hole shotgun source for engineering seismic surveys. Geophysics 52,985-996. Rorbaugh, C.B. 1993. Digital filter designer's handbook. Mcgraw-Hill. Russell, R.D. 1996. Speculations on Biot 'slow' waves. Accepted for publication in Canadian Journal of Exploration Geophysics in 1996. Russell, R.D. and Barker, A.S., Jr. 1991. Seismo-electric exploration: expected signal amplitudes. Geophysical Prospecting 39, 105-118. Russell, R.D., Maxwell, M . , Butler, K.E. , and Kepic, A.W. 1992. Electromagnetic responses from seismically excited targets A: piezoelectric phenomena at Humboldt, Australia. Exploration Geophysics, 23, 281-286. Sill, W.R. 1983. Self-potential modeling from primary flows. Geophysics, 48, 76-86. Sobolev, G.A. and Demin, V . M . 1980. Mechanoelectric Phenomena in the Earth. Nauka, Moscow (in Russian). A partial English translation is available at the University of British Columbia. References Sobolev, G.A. and Demin, V . M . 1980. Mechanoelectrical radiation from orebodies. Doklady AkademiiNauk SSSR, 252 , 1353-1355. Sobolev, G.A., Demin, V . M . , Los' and Maybuk, Yu.Ya. 1982. Study of the electromagnetic radiation of rocks containing semiconductor and piezoelectric minerals. Izvestiya, Academy of Sciences USSR (English translation by AX3.U.) 18, 888-897. Sobolev, G.A., Demin, V . M . , Narod, B.B. and Whaite, P. 1984. Tests of piezoelectric and pulsed-radio methods for quartz vein and base-metal sulfides prospecting at Giant Yellowknife Mine, N.W.T., and Sullivan Mine, Kimberley, Canada. Geophysics 49, 2178-2185. Stern, O. 1924. Ziir theorie der electrolytischen doppelschist, Z. Elektrochem., 30 , 508-516. Stoll, R.D. 1989. Sediment Acoustics, (Vol. 26 of the series Lecture Notes in Earth Sciences). Springer-Verlag, Berlin. Strack, K . M . , Hanstein, T.H., and Eilenz, H.N. 1989. L O T E M data processing for areas with high cultural noise levels. Physics of the Earth and Planetary Interiors, 5 3 , 261-269. Sumner, J.S. 1976. Principles of Induced Polarization for Geophysical Exploration. Elsevier, Amsterdam. Thompson, A . H . 1990. Electroseismic prospecting by detection of an electromagnetic signal produced by dipolar movement. US Patent 4 904 942, Feb. 27, 1990. Thompson, A . H . , and Gist, G.A. 1991. Electroseismic prospecting. 61st Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 425-427. Thompson, A .H . , and Gist, G.A. 1993. Geophysical applications of electrokinetic conversion. The Leading Edge, 12, 1169-1173. Thompson, R.R. 1936. The seismic-electric effect. Geophysics 1, 327-335. Thompson, R.R. 1939. A note on the seismic-electric effect. Geophysics 4, 102-105. Tome, R.F. 1978. Studies on the Seismoelectric effect. M.Sc. Thesis, Department of Physics, University of Toronto. Volarovich, M.P., Sobolev, G.A. and Parkhomenko, E.I. 1962. The piezoelectric effect in pegmatite and quartz veins. Izvestiya, Academy of Sciences USSR, Physics of the Solid Earth (English translation by A.G.U.), p. 145-152. White, D.L. 1962. Amplification of ultrasonic waves in piezoelectric semiconductors. Journal of Applied Physics, 3 3 , 2547-2554. Appendix A Notes on Fortran Programs for Data Processing The programs listed below were developed by the author to process seismoelectric data acquired with the RC Electronics (RCE) digitizer. A l l were written in Fortran and compiled for the MS-DOS operating system using the 'Microsoft Fortan' compiler (version 5.00). Each program reads processing parameters from a separate parameter file. The input/output data file format is the binary, multiplexed format used by RCE. The programs will read data in the offset 2-byte integer format used by RCE, as well as data that have been converted to 4-byte integers or 4-byte real numbers. The output data are written as 4-byte real numbers. SSUB.EXE Removes harmonic noise (such as powerline harmonics) from traces by sinusoid subtraction. See remarks at top of source file (SSUB.FOR) for more details. M B S U B . E X E Removes harmonic noise (such as powerline harmonics) from traces by block subtraction. See remarks at top of source file (MBSUB. FOR) for more details. SNOTCH5.EXE Calculates and applies a zero-phase frequency domain N O T C H filter to traces. Up to 30 different notches may be specified. The number of points in each trace must be a power of two (a requirement imposed by the FFT subroutine) and less than or equal to 32768. In order to handle such long traces, this program reads just one channel at a time, and performs all calculations in single precision (i.e. with REAL*4 numbers). SNOTCH is in fact an acronym for Single Precision Notch. DNOTCH9.EXE This program calculates and applies a zero-phase frequency domain notch. Up to 50 different notches may be specified. The FFT subroutine can handle traces having an arbitrary number of points (not necessarily a power of 2) up to a maximum of 8192. This is particularly useful if you have accurate knowledge of the fundamental frequency f0 you want to remove (e.g. 60 Hz for most powerline noise in North America). In that case the time series can be truncated to the longest length of time t that is an integer multiple of l/fo. Then, the discrete frequencies in the fft domain will include fo and all its harmonics, and those frequencies can be notched most effectively. This program is slower than SNOTCH5.EXE because the FFT subroutine takes an arbitrary number of points (it will be particularly slow if the number of points cannot be expressed as a product of small prime numbers). D N O T C H is an acronym for Double Precision Notch. See remarks at top of source file DNOTCH9.FOR for more details. M L P F D E C . E X E Calculates and applies a zero-phase frequency domain LOWPASS filter to traces. An FFT subroutine is used to convert to and from the frequency domain, and the number the number of samples per trace must be a power of 2. Sample interval can be reduced (decimated) after filtering. See remarks at top of source file (MLPFDEC.FOR) for more details. 154 Appendix A BPF.EXE Calculates and applies a zero-phase frequency domain BANDPASS filter to traces. The program structure is essentially the same as MLPFDEC.FOR and the number of samples per trace must be a power of 2. See remarks at top of source file BPFDEC.FOR for more details. URFILT.EXE Applies an Infinite Impulse Response (IIR) filter to traces. The filter is applied in the time domain. Filter coefficients must be supplied in a separate file having the following format: Line 1: Description of the filter (ASCII) Line 2: Filter order n (INTEGER) Line 3: bl,b2,...b(n+l) (numerator coefficients - REAL*8) Line 4: al,a2,...a(n+l) (denominator coefficients - REAL*8) Butterworth highpass and lowpass coefficients for a given order n and normalized cut-off frequency were obtained using the filter design function 'butter' in the mathematical software package 'Matlab'. See remarks at top of source file URFILT.FOR for more details. FILSTK3.EXE This program is used to stack multiple data files. The total stacks (not mean stacks) of corresponding channels in the files are calculated. The stacked data is written to an output file called FILESTK3.BIN. See remarks at top of source file FILSTK3.FOR for more details. LSTRSUB.EXE This program is used to replace one or more traces in a seismoelectric data file by the difference between that trace and a 'reference' or 'noise' trace (also taken from the file). The amplitude of the reference trace is automatically adjusted so that each difference trace is as small as possible in a least-squares sense. The program is intended for subtracting remote reference or 'noise' records from traces that contain both signal and regional noise. LSTRSUB is an acronym for 'Least Squares Trace Subtraction'. See remarks at top of source file LSTRSUB .FOR for more details. MULTQ3.EXE Like LSTRSUB.EXE, this program is used to make multiple 'quadrupole traces' by subtracting one or more reference traces from the other channels in a file. However, the user must specify how much gain to apply to the reference prior to subtraction. The program is intended for subtracting remote reference or 'noise' records from traces that contain both signal and regional noise. The program can also be used to simply amplify or invert the polarity of individual channels, or to change the number of channels in the output file (e.g. to write only a single channel to the output file). See remarks at top of source file MULTQ3.FOR for more details. Appendix B Mathematical Detail for Sinusoid Subtraction EVALUATION OF ERROR TERMS E\"P AND e* Equation (3.20), which gives the error in the estimate of an due to the presence of other powerline harmonics, may be rewritten in the form E\"P = $ra^ckfcos(kc»j+ 0t+§ k -ya)]dt. ^ k=\\ k*n Evaluating this integral over the length of the estimation window [0,x] yields the result k*n sin((£ + n)coox +<\\>k + y a) -sin()p, + 4> pf = 1990 kg/m 3 (bulk density) P i2 = ( l - ^ P / =-300 kg/m 3 Pn = (l-)P,-Pi2= 1890 kg/m 3 p 2 2 = §pf - p 1 2 = 700 kg/m 3 Note that, according to equation (5.31), the values assigned to the permeability K and inertial coupling term pn have no direct effect on the magnitude of the electrokinetic potential \\|/ associated with the fast compressional wave at low frequencies. They affect the amplitude of the electrokinetic field -V\\ | / (equation (5.33)) in only a minor way through their influence on the wave velocity v. Thus the values used for K and pn were not very important for calculating the estimates of electrokinetic signal strength given in this thesis. However, it is important to appreciate that, in real rocks, K and p i 2 are very dependent upon the tortuosity x which does appear explicitly in the equations for the electrokinetic potential and field. Notes regarding estimates for the physical properties in Table DI Electrical Conductivity a Electrical conductivity was measured directly as a function of depth with a borehole logging tool in boreholes FD92-3 and FD94-3. It was approximately 10 mS/m throughout the Holocene sands at 2- 20 m depth. Porosity (J) The GSC calculated porosity as a function of depth in borehole FD94-3 from measurements of the weight percent water in core samples taken at 1 m depth intervals (Dallimore et al., 1995). This method gives accurate porosities provided that the core is fully saturated, that all water is successfully extracted, and that the density of the solid grains is known. (A density of 2650 kg/m3 - appropriate for quartz sand grains - was assumed.) The porosity was approximately 0.4 throughout the Holocene sands above 20 m depth. Appendix D Shear Modulus [im The shear modulus of the porous frame is assumed to be independent of water saturation (a common assumption that is reasonable for our purposes). The given estimate for \\x.m is based measurements of the shear wave velocity vs in borehole FD94-3 (Dallimore et al., 1995). In the Holocene sands above 20 m depth vs ~ 180 ± 30 m/s. Since vs = ^ [im/p , and p - 2000 kg/m 3, we find u.m - 6.5 x 107 N/m 2 . Bulk modulus of the dry porous frame Km This value cannot be obtained directly from in-situ measurements of P-wave velocity because the Holocene sands were at least partially saturated at all depths. Instead we rely on an empirical relationship between \\im and Km. Stoll (1989, p. 92) states that granular materials at relatively low confining pressures generally have Poisson's ratios •& in the range 0.1 to 0.2; given that Kj\\x.m =[2(l+'Dv)]/[3(l-2'r>)], this means that Km can be expected to lie in the range 0.9u.m ^Km< 1.3p:m. We chose an intermediate value Km = 1.02p:m to obtain the estimate Km = 6.6 x 107 N/m 2 given in Table DI . Note that this estimate, together with those for the other elastic constants and 0, yields a characteristic velocity v c of 1630 m/s for the fast compressional wave (see Table D2). This velocity seems reasonable because it is just a little higher than the measured P-wave velocities of 1500-1570 m/s obtained in the Holocene sands below the water table. Permeability K The permeability of the sand at this site is not known. The estimated value of 5 darcy (5 x 10' 1 2 m2) is a guess based on the fact that the sand is described as predominantly massive, fine to medium grained sand with occasional silt interbeds (Hunter et al., 1994; Dallimore et al, 1995). According to Freeze and Cherry (1979) typical permeabilities in silty sands range from about 10\"2 to 102 darcy, while those in clean sands range from about 10\"1 to 103 darcy. Tortuosity x For a medium composed of spherical grains, x = (l +1 /()) )/2 (Berryman, 1980). This equation, together with (j) = 0.4, gives an estimate of 1.75 (rounded to 1.8) for the tortuosity. Zeta potential C, Laboratory measurements of streaming potentials indicate that zeta potentials of -10 to -100 mV are to be expected in water-saturated geologic materials (Ishido and Mizutani, 1981; Morgan et al., 1989). We have selected an intermediate value, £ = -50 mV, for the estimate in Table D I . This value, together with those given for given for (J), £/, rj, and T, yields a streaming current cross-coupling coefficient X = -§eft,/(n^2) of 4.6 x 10~9 m2/(Vs) which is quite reasonable; for example, measurements of X in 14 different fine-grained sediments are listed by Mitchell (1993, p. 270) and most lie in the range 2 x 10\"9 to 10 x 10\"9 m2/(Vs). "@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "1996-11"@en ; edm:isShownAt "10.14288/1.0087918"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Environmental Sciences"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Seismoelectric effects of electrokinetic origin"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/6270"@en .