Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The application of seismic techniques to hydrogeological investigations Jarvis, Kevin Donald Gibson 2001

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2001-611191.pdf [ 19.16MB ]
Metadata
JSON: 831-1.0053191.json
JSON-LD: 831-1.0053191-ld.json
RDF/XML (Pretty): 831-1.0053191-rdf.xml
RDF/JSON: 831-1.0053191-rdf.json
Turtle: 831-1.0053191-turtle.txt
N-Triples: 831-1.0053191-rdf-ntriples.txt
Original Record: 831-1.0053191-source.json
Full Text
831-1.0053191-fulltext.txt
Citation
831-1.0053191.ris

Full Text

THE APPLICATION OF SEISMIC TECHNIQUES TO HYDROGEOLOGICAL INVESTIGATIONS By Kevin Donald Gibson Jarvis B.Eng., The University of Saskatchewan, 1986 M . S c , The University of British Columbia, 1995  A THESIS SUBMITTED IN P A R T I A L F U L F I L L M E N T OF THE REQUIREMENTS FOR T H E D E G R E E OF DOCTOR OF PHILOSOPHY in T H E 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  T H E UNIVERSITY OF BRITISH C O L U M B I A April 2001 © Kevin Donald Gibson Jarvis, 2001  In  presenting  degree freely  at  this  the  University  available  copying  of  department publication  for  this or of  thesis  this  of  reference  thesis by  in  for  his thesis  partial  fulfilment  of  the  British  Columbia,  I  agree  and  scholarly  or for  her  of  T h e U n i v e r s i t y o f British Vancouver, Canada  Date  DE-6  (2/88)  Columbia  I  further  purposes  gain  shall  that  agree  may  representatives.  financial  permission.  Department  study.  requirements  It not  be  that  the  Library  permission  granted  is  by  understood be  for  allowed  an  advanced  shall for  the that  without  make  it  extensive  head  of  my  copying  or  my  written  ABSTRACT The objective of this thesis is to demonstrate some new applications of seismic techniques for hydrogeological applications.  A compressional-wave,  surface-based,  reflection seismic  technique is used to map aquifer boundaries within a series of Pleistocene sediments.  near-surface  The interpretation uses both water wells and sequence stratigraphic concepts to  identify the boundaries of new and existing aquifers. The use of the cone penetrometer is an integral part of this thesis.  The seismic cone is demonstrated to be both cost-effective and  reliable for the acquisition of high-quality vertical seismic profile (VSP) data. Other data from the cone, in particular the tip resistance data, are shown to be an integral link for the conversion of shear-wave velocities to values of hydraulic conductivity.  Surface-based,  shear-wave  reflection seismic data are used to image an aquifer contained within Holocene deltaic sediments. A Bayesian inversion of the shear-wave seismic amplitudes (using cone-derived velocities) results in the generation of a two-dimensional profile of shear-wave velocity that is a direct indication of aquifer heterogeneity. Conversion of the velocity to hydraulic conductivity (using a cone-derived relationship) results in the distribution of a key hydrogeologic property within the aquifer. The results from the thesis show significant promise for improving groundwater flow models and providing new techniques for the management and protection of our groundwater resources.  n  TABLE OF CONTENTS ABSTRACT  ii  T A B L E OF CONTENTS  iii  LIST OF FIGURES  vi  ACKNOWLEDGEMENTS  ix  C H A P T E R 1: Introduction  1  C H A P T E R 2: Using Compressional-Wave Seismic Data to Locate Aquifer Boundaries  5  2.1  Introduction  5  2.2  Geologic Setting  6  2.3  Seismic Data Acquisition  8  2.4  Seismic Data Processing  9  2.5  2.4.1  Geometry Assignment  10  2.4.2  Amplitude Recovery and Filtering  11  2.4.3  Deconvolution  12  2.4.4  Elevation Correction  12  2.4.5  Re-Sort to C M P Gathers  13  2.4.6  N M O Removal and Residual Statics Correction  13  2.4.7  First-Break Muting and Trace Stacking  14  2.4.8  Post-Stack Noise Reduction  14  2.4.9  Migration  ••••15  2.4.10 Amplitude Spectra  16  Seismic Interpretation  16  2.5.1  Interpretation ofLine20-206  18  2.5.1.1  Nearby water well records  19  2.5.1.2  Identification of a deep, undrilled aquifer  19  2.5.2  Interpretation of Line 216-28  20  2.5.2.1  20  Nearby water well records iii  2.5.2.2 2.5.3  2.5.4 2.6  Identification of discrete aquifers within a glacial till  21  Interpretation of Line 32-184  21  2.5.3.1  Nearby water well records  21  2.5.3.2  Identification of thinly-bedded sand aquifers  22  2.5.3.3  Identification of a high-quality sand and gravel aquifer... 23  2.5.3.4  Identification of deep, undrilled aquifers  23  Summary of Seismic Interpretations  24  Conclusions  24  C H A P T E R 3: Near-Surface VSP Surveys Using the Seismic Cone Penetrometer  34  3.1  Introduction  34  3.2  Geologic Setting  35  3.3  Cone Penetrometer Data  36  3.3.1  VSP Data Acquisition  37  3.3.2  VSP Data Processing  38  3.3.3  Cone Data Acquisition  41  3.4  SH-Wave Seismic Reflection Data  42  3.4.1  Data Acquisition  42  3.3.2  Data Processing  43  3.5  A Comparison of the Geotechnical Logs and VSP Data  44  3.6  A Comparison of VSP and CDP Data  47  3.7  Conclusions  49  C H A P T E R 4: Aquifer Heterogeneity from SH-Wave Seismic Impedance Inversion  57  4.1  Introduction  57  4.2  Geologic Setting  60  4.3  SH-Wave Seismic Reflection Data  61  4.3.1  Data Acquisition  61  4.3.2  Data Analysis and Processing  62  4.4  Seismic Cone Penetrometer Data  63  4.4.1  Data Acquisition  ,  4.4.2  Interval Velocity Logs  65  4.4.3  VSP Data Analysis and Processing  65  iv  64  4.5  Impedance Inversion of Shear Wave Data  66  4.6  Differentiating Lithology Using the Velocity Model  70  4.6.1  Correcting for Effective Stress in Velocity Data  71  4.6.2  Void Ratio in Sand Aquifer  73  4.7  Conclusions  75  C H A P T E R 5: Aquifer Hydraulic Conductivity Using Cone Penetrometer and Seismic Data .... 85 5.1  Introduction  85  5.2  Geologic Setting  89  5.3  Cone Penetrometer Data  90  5.3.1  Data Acquisition  90  5.3.2  Tip Resistance and Shear-Wave Velocity Relationship  91  5.3.2.1  Removing the effect of effective stress  91  5.3.2.2  Normalized tip resistance and shear-wave velocity plot... 93  5.4  SH-Wave Seismic Reflection Data  95  5.4.1  Data Acquisition  95  5.4.2  Data Processing  .96  5.5  Impedance Inversion of Shear-Wave Data  5.6  Aquifer Hydraulic Conductivity  100  5.6.1  Tip Resistance Model of Aquifer  100  5.6.2  Tip Resistance and Grain Size Relationship  101  5.6.3  Hydraulic Conductivity Model of Aquifer  103  5.7  Conclusions  97  106  C H A P T E R 6: Conclusions and Recommendations for Future Work  115  REFERENCES  117  v  LIST OF FIGURES Figure 2.1:  Figure 2.2:  Figure 2.3:  Location map showing the outline of the Brookswood aquifer and the location of the three seismic lines: 20-206, 216-28 and 32-184  25  Representative raw shot profiles (with 0.04 s A G C ) from lines 20-206, 21628 and 32-184  26  The shot profiles from Figure 2.2 with a bandpass filter applied (and 0.04 s  AGC)  '..J..'  26  Figure 2.4:  Processed seismic data for line 20-206; a) stack b) migration  27  Figure 2.5:  Processed seismic data for line 216-28; a) stack b) migration  28  Figure 2.6:  Processed seismic data for line 32-184; a) stack b) migration  29  Figure 2.7:  Average amplitude spectra of the migrated seismic data; a) 20-206 b) 21628 c) 32-184  30  Figure 2.8:  Interpretation of line 20-206  31  Figure 2.9:  Interpretation of line 216-28  32  Figure 2.10: Figure 3.1:  Interpretation of line 32-184 Location map showing the Kidd2 research site on the Fraser River delta with the grid of SH-wave CDP seismic (dotted lines) and the two VSPs, labeled VSP1 and VSP2  33  50 Figure 3.2:  Figure 3.3:  Data from VSP1. A 30-150 Hz bandpass filter and spherical divergence correction have been applied to the data  51  Data from VSP2. A 30-150 Hz bandpass filter and spherical divergence correction have been applied to the data  51  Figure 3.4:  Interval velocities derived from the VSP first break arrival times  52  Figure 3.5:  Processing flow chart for the cone penetrometer VSP data  52  Figure 3.6:  Processed data from VSP1 (bottom). The upgoing wavefield has been aligned and muted. The stacked data trace is duplicated and is shown above (VET1)  53  Processed data from VSP2 (bottom). The upgoing wavefield has been aligned and muted. The stacked data trace is duplicated and is shown above (VET2)  54  Figure 3.7:  vi  Figure 3.8:  Two representative shot profiles from line 20 after application of gain, trace balancing and record subtraction. Shot profile 1004 is at the south end of the line and shot profile 1252 is at the north end  55  Composite plot showing the relationship of the cone penetrometer friction ratio (FR), tip resistance (Q ), and pore pressure (U2) to the first break shear-wave velocity (V ) and V S P l ' s extracted traces (VET1). The data are plotted in time with the equivalent depths indicated on the left  55  Composite plot showing the relationship of the cone penetrometer friction ratio (FR), tip resistance (Q ), and pore pressure (U2) to the first break shear-wave velocity (V ) and VSP2's extracted traces (VET2). The data are plotted in time with the equivalent depths indicated on the left  56  SH-wave CDP reflection data with the extracted traces from VSP1 and VSP2 (VET1 and VET2) inserted  56  Location map showing the Kidd2 research site on the Fraser River delta with the grid of SH-wave CDP seismic (dotted lines) and the two VSPs, labeled VSP1 and VSP2 that were used for impedance constraints. The seismic line (line 20) used in the inversion has been highlighted  77  Two representative shot profiles from line 20 after application of gain, trace balancing and record subtraction. Shot profile 1004 is at the south end of the line and shot profile 1252 is at the north end  77  The migrated SH-wave seismic reflection line 20 that was used as the input to the seismic impedance inversion  78  Interval velocity logs from the two seismic cone penetrometer holes (VSP1 and VSP2); a) standard b) normalized  79  Figure 4.5:  Input velocities to the seismic impedance inversion  80  Figure 4.6:  A comparison of the velocity logs from the two VSPs [a) VSP2 b) VSP1] to the output from the seismic impedance inversion. The line with diamonds is the velocity log and the solid line is the inversion output  81  Figure 3.9:  c  s  Figure 3.10:  c  s  Figure 3.11:  Figure 4.1:  Figure 4.2:  Figure 4.3:  Figure 4.4:  Figure 4.7:  The output impedance model convolved with the source waveform.  This  data closely matches the input data shown in Figure 4.3  82  Figure 4.8:  Velocities output from the seismic impedance inversion  83  Figure 4.9:  Normalized velocities output from the seismic impedance inversion. The sand aquifer can be clearly identified by the yellow and red colors  83  vn  Figure 4.10:  Lithology/void ratio profile obtained by separating the sand and clay zones using a cutoff velocity of 170 m/s and using Equation 4.7 to obtain void ratio within the sand aquifer  84  Figure 5.1:  Location map showing the Kidd2 research site on the Fraser River delta with the grid of CDP seismic (dotted lines) and the cone holes that were used to constrain the inversion ( K l , K 2 , K9404, K9603). The seismic line (line 20) used in the inversion has been highlighted as well as the location of the borehole (CORE) 107  Figure 5.2:  Plot of the normalized shear-wave velocity ( V ) vs. normalized tip resistance (Q ). The line represents a least squares fit of the data 107 sn  tri  Figure 5.3:  Plots showing a comparison between the tip-resistance-derived shear-wave velocity (solid line) and the cone-measured shear-wave velocity (dashed line with diamonds) for a) K1 and b) K 2 108  Figure 5.4:  The migrated SH-wave seismic reflection line 20 that was used as the input to the seismic impedance inversion. The four cone holes that were used to constrain the impedance inversion are highlighted 109  Figure 5.5:  Normalized velocities input to the seismic impedance inversion  Figure 5.6:  Comparisons of the velocity models input to the impedance inversion (thin lines) to the velocity models output from the impedance inversion (thick lines) at the four cone hole locations: a) K 2 b) K l c) K9404 d) K9603 110  Figure 5.7:  Normalized shear-wave velocities output from the impedance inversion  Figure 5.8:  Normalized tip resistance computed from the normalized velocities in Figure 5.7 Ill  Figure 5.9:  Detailed comparisons of the original tip resistance logs (light solid lines) to the tip resistance logs computed from the normalized velocities output from the impedance inversion (thicker lines with diamonds): a) K 2 b) K l c) K9404d)K9603 112  Figure 5.10:  Plot of grain size (dso) versus normalized tip resistance (Qm)- The line represents a least squares fit of the data 113  Figure 5.11:  Color plot showing the hydraulic conductivities computed from the normalized shear-wave velocities output from the impedance inversion 113  Figure 5.12:  Color plot showing the potential error in the computed hydraulic conductivities. The plot was generated by running two inversions with aquifer velocities of +/-8 m/s with respect to the original inversion. The hydraulic conductivities were then calculated and the difference between the two results is illustrated in the plot. The contour interval is 0.02 114  vm  109  Ill  ACKNOWLEDGEMENTS I would like to express my gratitude to the entire Rock Physics Research Group at U B C for providing ideas and encouragement over the past six years. In particular, my research advisor, Rosemary Knight was a constant source of enthusiasm and support. This research was funded primarily by an N S E R C Collaborative Research Grant to Rosemary Knight with additional funding obtained from an N S E R C Research Grant. I was personally supported by an N S E R C postgraduate scholarship.  I would like to thank Golder  Associates for providing the Geometries seismograph, Professor John Howie who coordinated the usage of the U B C Department of Civil Engineering cone truck, and Mauricio Sacchi for providing the Bayesian impedance inversion algorithm. The data shown in this thesis required the dedicated effort and assistance from a number of field crews. The P-wave seismic reflection crew had to endure the brambles and ditches south of Langley and consisted of: Gwenn Flowers, Andrew Gorman, Kevin Kingdon, Wendi Milner, Sean Rastead, Jane Rea, Paulette Tercier and Dan Woznow from the Department of Earth and Ocean Sciences. The SH-wave seismic reflection crew had to endure an asthma and hay-fever inducing grassy field and consisted of: Ian Blumel, David Butler, Christina Chan, Jay Joyner, Jane Rea, Richard Taylor and Paulette Tercier from the Department of Earth and Ocean Sciences.  The cone penetrometer data was collected with the assistance of: Chris Daniel,  Heraldo Giacheti, Scott Jackson and Harald Schrempp from the Department Engineering.  ix  of Civil  CHAPTER 1 Introduction The development of accurate models of groundwater flow is of primary importance for the management and protection of our groundwater resources.  A n accurate groundwater model  requires two types of information: the location of aquifer boundaries and (at a smaller scale) the distribution of hydro geologic properties within the aquifer. The objective of this thesis is to detail some non-traditional techniques for providing both types of information. Of specific interest in my research are groundwater aquifers found in near-surface sediments. The traditional geologic technique used to locate aquifer boundaries is the drilling of a well. Ideally, a continuous core will be extracted as the well is drilled, that provides unambiguous detail about the vertical variation in the subsurface sediments. Often a core is not extracted (due to the expense) and instead the hole is logged to varying degrees of detail. The simplest logging technique involves monitoring the drilling rate and analyzing the cuttings of the drill head that are expelled by the drill fluid. More complicated techniques involve the use of geophysical logs to infer the types of sediments that surround the well bore (see Freeze and Cherry, 1979). The distribution of hydrogeologic properties within an aquifer can only be determined by having data away from the well bore. Additional wells can be drilled to assess the regional complexity and additional experiments (e.g. pumping tests and tracer tests) can be performed to assess the aquifer properties but the drilling of well bores is an expensive proposition.  The  ability to locate aquifers before any drilling commences and to model groundwater flow with limited subsurface control are desirable objectives.  Geophysical techniques, in particular  geophysical techniques used on the surface provide the means to assess subsurface aquifers with little or no subsurface control.  1  There are a variety of different geophysical techniques capable of locating subsurface aquifers. I have chosen to use the seismic method for aquifer imaging for two key reasons. The first and most important reason relates to the scale of information that can be obtained using the seismic method. The seismic method is an elastic wave-based imaging technique involving the detection of waves after they have passed through a region of the subsurface. The interaction of the waves with subsurface changes in velocity and density provides detailed information at a scale governed by the wavelength of the propagating waves.  The second reason is the proven  ability of seismic methods for providing images of aquifers in the study area (Pullan et al., 1995). In solid materials two different wave types are possible: compressional or P-waves and shear or S-waves.  The waves are generated by a source at one location and sensed by receivers  (referred to as geophones) at another. P-waves are generally easier to generate using a source such as a sledgehammer or a buffalo gun (Pullan and MacAulay, 1987) but S-waves can provide greater detail in saturated sediments (Dobecki, 1993). Pullan et al. (1995) and Harris et al. (1998) demonstrated how a portable engineering seismograph and a small number of geophones could be used to obtain excellent quality P-wave and S-wave seismic reflection data in the Fraser lowland. I will be demonstrating the use of both P-wave and S-wave seismic methods to locate aquifer boundaries and map out variations in the hydrogeologic properties of the aquifers. The challenge in using seismic reflection data to obtain hydrogeologic parameters such as porosity and hydraulic conductivity is the development of a model that relates the seismic parameters of velocity and density to the hydrogeologic parameters.  This relationship is best  developed by having independent subsurface control to make the problem tractable.. This control is often obtained by the geophysical logging of well bores. A n alternative to this expensive technique is the cone penetrometer, which was developed by geotechnical engineers for rapid subsurface characterization (for a recent review see Lunne et al., 1997). The cone penetrometer involves the pushing of an instrumented steel cone into the ground and making measurements as 2  the cone gets progressively deeper. It works extremely well in finer-grained deltaic sediments (Campanella et al., 1983) but has difficulty in coarser-grained gravel or compacted till sediments and cannot be used at all in rock. To use the cone penetrometer data to constrain the analysis of seismic reflection data requires that relationships be established between the fundamental parameters from each technique.  Establishing some of these relationships and using them  specifically for hydrogeological applications is one of the key contributions of this thesis. The focus of this thesis is the application of seismic methods to hydrogeologic investigations. In Chapter 2 I examine the use of compressional-wave (P-wave) seismic data to locate aquifer boundaries in a number of different sedimentary environments surrounding the Brookswood aquifer (Fraser valley, SW British Columbia). In Chapter 3 I investigate the use of the cone penetrometer to acquire vertical seismic profile (VSP) data within a near surface aquifer on the Fraser delta (Kidd2 research site). Chapter 4 extends the application of the cone penetrometer by demonstrating how cone velocity data can be used in the analysis of a surface based seismic reflection profile to map out the heterogeneity of the Kidd2 aquifer.  Chapter 5 further extends  the integration between the cone and seismic reflection techniques by using the cone tip resistance both to develop a subsurface velocity model and to convert an optimized velocity model to a model of hydraulic conductivity. A l l figures are located at the end of each chapter. There is some repetition of figures and introductory material between the chapters of the thesis. The chapters were written as complete research papers and due to the similarity of the location and subject matter the repetition was unavoidable. The ideas and techniques developed in this thesis provide innovative and cost-effective tools for locating and establishing the distribution of hydrogeologic properties of subsurface aquifers. The judicious application of these techniques in groundwater investigations will have a significant impact on the quality of groundwater flow models developed. The models can then  3  serve as a key component in making decisions to protect and manage our groundwater, an increasingly valuable resource.  4  CHAPTER 2 Using Compressional-Wave Seismic Data to Locate Aquifer Boundaries 2.1  Introduction Locating aquifer boundaries without drilling a well is a desirable objective. This can be  accomplished by using surface-based geophysical techniques.  In particular, the seismic  technique, which involves the propagation and measurement of elastic waves in the subsurface, will be demonstrated in this chapter to have significant promise. This chapter summarizes and builds upon my contribution to Pullan et al. (2000) as part of the Fraser Lowland Hydrogeology Project (Ricketts, 2000) and involved the acquisition of three seismic lines close to the Brookswood aquifer, near Langley, B C . The primary objective of the Fraser Lowland Hydrogeology Project (Ricketts, 2000) was to investigate different techniques for the characterization of aquifers and aquitards. The focus was to examine the 3-dimensional geometry, stratigraphic position, sedimentary character, hydraulic properties, recharge characteristics and groundwater types. The study area involved the Fraser River lowland southeast of Vancouver. In order to meet the objectives of the study a wide variety of different data were utilized including water-well records (Woodsworth and Ricketts, 1994; Ricketts, 1995; Ricketts and Dunn, 1995; Makepeace and Ricketts, 2000; Ricketts and Makepeace, 2000), new drilling and water sampling, cone penetrometry (Campanella et al., 1994) and geophysical surveys (Rea et al., 1994; Best et al., 1995; Pullan et al., 1995; Best and Todd, 2000; Rea and Knight, 2000; Roberts et al., 2000). The geophysical surveys used both electromagnetic (EM-47 and GPR) and seismic techniques. The focus of Pullan et al. (2000) was the use of the compressional-wave (P-wave) seismic reflection technique to provide highresolution images of aquifers and aquitards in the Fraser lowland. A total of 14 seismic lines were interpreted to identify regional geological trends and major stratigraphic units. The coarse5  grained and fine-grained sediments were identified without making references to specific aquifers.  M y contribution to Pullan et al. (2000) was the acquisition, processing and basic  interpretation of three of the seismic lines. In this chapter I will fully expand on the seismic processing sequence used on the three lines and provide a detailed interpretation with specific references to nearby water wells and aquifers. The P-wave seismic reflection technique has the potential to produce very high quality subsurface images. If good source and receiver coupling can be achieved and the subsurface sediments are not too attenuative then excellent data can be acquired.  Many researchers have  demonstrated the viability of near-surface P-wave seismic reflection surveys (e.g. Steeples and Miller, 1990; Roberts et al., 1992; Pullan et al., 1994; Bradford et al., 1998; Cardimona et a l , 1998; Biiker et al., 1998).  A consistent observation by researchers is that the highest quality  near-surface seismic data are obtained in areas with surface sediments that are fine-grained and saturated. The Fraser Lowland is characterized by a variety of different surficial sediments most of which are fine-grained and saturated, therefore making it an ideal location to acquire high quality P-wave seismic reflection data. In this chapter I will detail the local geological setting, seismic acquisition parameters and discuss the three lines, which have been acquired and the hydrogeologic information that was obtained from the data. The seismic data presented in this chapter has also been presented and discussed in a Geological Survey of Canada publication (Pullan et al., 2000).  2.2  Geologic Setting  The seismic investigation was carried out on and nearby the Brookswood aquifer, just south of Langley, British Columbia, Canada.  Figure 2.1 shows the outline and location of the  Brookswood aquifer in the Fraser lowland with the location of the three seismic lines highlighted. The Brookswood unconfined aquifer represents the deltaic deposits generated by a  glacial meltwater channel (Roberts et al., 2000). The surficial sediments to the west are part of the Capilano Formation and represent offshore marine sediments deposited at the same time as the Brookswood. The surficial sediments to the east are part of the Fort Langley Formation and consist of marine sediments with greater amounts of coarser-grained material than the Capilano Formation. The Fort Langley sediments are the landward equivalent of the Capilano sediments and there is evidence of compaction by glacial ice (Campanella et al., 1994). The boundary between the two formations is the approximate shoreline at the time of Brookswood deposition. The uppermost sediments of the Capilano Formation have never been compacted by an overlying ice mass and consist primarily of clay. The relatively low surface elevation where Capilano sediments are located places the water table very close to the surface resulting in the full saturation of the near-surface sediments. The Brookswood unconfined aquifer is the most obvious source of groundwater in the area but by no means the only source of groundwater. Underlying the surficial sediments lie a sequence of sedimentary deposits containing a number of confined aquifers. These deeper, "hidden", aquifers are the most challenging to locate and those interested in a reliable source of groundwater will benefit the most using a remote sensing technique such as the seismic method. The deeper sedimentary sequence can be characterized by a series of glacial advances and retreats within the Fraser valley during the Pleistocene (Armstrong, 1984; Clague, 1994). When glacial ice covered the Fraser lowland the surface of the earth was pushed lower to compensate (isostatic adjustment). During interglacial periods the ice melted but the depressed land surface required time to readjust to the reduced loading. At these times the sea level was relatively high with a shoreline far up the Fraser valley. As the surface gradually rebounded isostatically the shoreline advanced down the valley.  The various depositional environments created by the  changing position of the shoreline resulted in a variety of distinct sediment types.  Marine  sediments were typically fine-grained silts and muds. Stream sediments varied from sands to 7  coarse gravels and nearshore sediments were typically sandy in texture.  From a hydrogeologic  standpoint, the sands and gravels are the primary aquifers and the silts and clays from the deeper marine deposits are the primary aquitards.  2.3  Seismic Data Acquisition The common depth point (CDP) compressional-wave (P-wave) seismic data were collected  in June of 1995 using a 24 channel Geometries SmartSeis seismograph, 48 Mark Products vertical geophones (L40A2-100 Hz) attached to CDP spread cables and a Geostuff roll switch. The data sample interval was 0.00025 s and a 1200 Hz high-cut anti-alias filter was applied before analog to digital conversion. A split-spread shooting configuration was used with a geophone and source spacing of 2 m, giving a nominal far offset of 23 m and 1200% subsurface coverage. The seismic source consisted of a buffalo gun (Pullan and MacAulay, 1987) using 12 gauge popper load blanks. A shothole was drilled for every shot using a portable post-hole auger with a 5 cm bit. The shothole depths varied depending on drilling conditions but were typically in the range of 0.7 to 1.0 m deep. A rubber mat with a hole cut out was placed over the shothole to prevent the release of high-velocity sand and gravel during detonation.  A total of 3 lines of  data were acquired with subsurface coverage of approximately 600 m. The location of these seismic lines can be seen in Figure 2.1. A l l 3 surveys were carried out in ditches adjacent to roadways, providing easy access to the sites but also resulting in acquisition delays during times of high traffic volume or the movement of large trucks and equipment on the road.  The ditches often contained fine-grained and  saturated material, providing an excellent medium to firmly plant geophones and to set off wellcoupled shots. The ditch along 32  nd  Ave had a stream running down it (outflow from the  Brookswood aquifer) with a higher stream velocity in the steeply graded, eastern portion of the line. The movement of water against the geophones was undoubtedly a large source of noise. 8  Having the shots detonated under water seemed to provide the best seismic signal on the records so the shots were saturated whenever possible. At the beginning and end of each line the split-spread geometry was not maintained. The active patch of geophones (connected by the roll-box) consisted of the first or last 24 geophones on the line. At the start of a line the first shot was between the first and second geophone and then subsequent shots moved through the spread until it reached the middle of the spread, after which the spread was then rolled to maintain the split-spread geometry. At the end of the line, when rolling the spread would result in non-existent traces, the patch remained stationary as the shot was moved from the middle of the spread to the last geophone.  Since the planting of  geophones is one of the most time-consuming aspects of reflection seismic acquisition, moving the shot through the spread at the beginning and end of the lines is the best way to maintain a reasonable level of subsurface coverage. The elevations of all the geophones were measured by running a level survey after the geophones had been planted.  The survey rod was placed on the top of each geophone being  surveyed. This procedure ensured that geophones planted deep into muddy ditches and those moved to avoid obstacles (e.g. gravel, piping, buried cables) had the best possible elevations assigned to them. The final relative elevations were assumed to be accurate to within 5 cm, which should be suitable to properly image even the highest frequencies contained within the data.  2.4  Seismic Data Processing  The seismic data were processed using a fairly standard processing scheme used in the processing of seismic for the oil and gas industry (Yilmaz, 1987) as well as for shallow seismic research (Steeples and Miller, 1990).  The objective was to provide the highest resolution  possible and to image reflections very close to the surface. A number of characteristics of the  data, some of which were location dependent, limited both the resolution-and shallow imaging capabilities.  A l l of the processing was done using the ITA commercial seismic processing  software accessible on a Sun workstation network installed as part of the University of British Columbia Geophysical Research Processing Facility. The end result of the seismic processing was three seismic sections that imaged different subsurface stratigraphy within and around the Brookswood aquifer and illuminated a series of known and potential subsurface aquifers. Each of the main processing steps is detailed below.  2.4.1  Geometry Assignment  The first step in the processing of the data involved the assignment of line geometry.  Well-  documented field notes were critical in assigning the proper set of geophones to each shot. As with all field surveys the shots did not always go according to plan. Shots were occasionally skipped and the roll box position was not necessarily set properly each time. Misplaced shots were interpreted by examining the raw shot profiles and identifying the shot location on the seismic spread based on the projection of the direct arrivals to zero time. Some example raw shot profiles from each line are shown in Figure 2.2. The only process that was applied to the data was the application of a 0.04 s automatic gain control (AGC). The A G C performs a balance of amplitudes to make all potential reflections visible. The ground roll is clearly visible on all three profiles as a cone of low frequency data extending away from the shot location, covering up the reflection energy. Hyperbolic reflections can be seen on each of the profiles at the farther offsets, outside of the ground-roll energy. The A G C was only used in this instance for display purposes and was not part of the processing flow.  10  2.4.2  Amplitude Recovery and Filtering  The gaining and deterministic filtering of the data were key processes in obtaining the highest possible resolution and preserving the relative amplitudes. The traces were gained using a time-power correction whereby every sample is multiplied by its time raised to a power. The power used was 1.8, which was determined by trying a number of different values and assuming that the deeper and shallower packages of reflections should have a similar range of amplitudes. The objective of the deterministic filtering was to reduce the amplitude of frequencies dominated by ground-roll prior to deconvolution. The deconvolution process attempts to estimate the signal spectrum so reducing the amplitude of the ground-roll frequencies will help the deconvolution algorithm to do a better job (Yilmaz, 1987). The data that were attenuated were generally associated with source-generated noise. This noise consisted of both ground-roll and air-wave energy that were generally comparable in amplitude or significantly higher than the surrounding reflection amplitudes.  The ground-roll  energy tended to be lower in frequency than the reflection signal so the filtering easily reduced the amplitudes of the ground-roll. However, the near-shot ground-roll (ie. within 5 m of the shot location) still possessed significant amplitudes. In addition, the air-wave energy (a linear event with a velocity equal to the velocity of sound in air - approximately 340 m/s), which contained frequencies similar to the reflection frequencies, is most evident near the shot. The same three shot profiles (Figure 2.2) after bandpass filtering are shown in Figure 2.3. The filter parameters used for each line are annotated below the profiles. As with Figure 2.2 an 0.04 s A G C was applied for display purposes only. The ground-roll energy has been greatly reduced in amplitude and the underlying reflection energy is now visible across each of the three profiles. The air-wave energy is apparent on each profile as a series of high-amplitude, linear events extending away from the shot locations. To remove the effects of the non-reflection energy on the final products the near-source noise energy was muted. 11  The mute was as  conservative as possible to prevent removing too much shallow data.  However, the application  of this process is the main reason that reflections shallower than 10 m could not be imaged on the final products.  2.4.3  Deconvolution  A spiking deconvolution was applied to the data to maximize the bandwidth.  The  deconvolution optimally whitened the data and balanced the waveform from one shot to another (see, e.g., Yilmaz, 1987). Single reflections prior to deconvolution were split into a series of reflections with a very compact (broadband) waveform. The background noise also became more evident. This resulted from the increase in amplitude of higher frequencies with a lower signal to noise ratio, thereby causing the noise to be more noticeable. The stacking process eventually reduced the noise (due to its incoherence), so having the noise in the data was not a concern.  When the deconvolution filter was applied to the data the large amplitude changes  (e.g. those associated with noise spikes or mute zone edges) became apparent in the data as high frequency amplitude bursts.  These non-reflection events were muted.  The deconvolution  process also created large changes in overall trace amplitude so a trace balance was applied to the deconvolved, muted data.  2.4.4  Elevation Correction  To compensate for the elevation variations along the lines a simple static correction was applied to the data.  The near-surface velocities appeared to be relatively constant (from an  analysis of the first-break energy), with no obvious layering.  As a result, an adjustment to a  constant datum was made by assuming a constant velocity (approximately 1000 m/s). To avoid violating the hyperbolic shape assumption of the velocity removal algorithm, a datum was  12  chosen that represented an average elevation of the entire seismic line. This type of assumption is valid only if the topography is not too significant (which it wasn't).  2.4.5  Re-Sort to C M P Gathers  From this point onwards, the data were re-sorted into common mid-point (CMP) gathers. Each gather consisted of traces from the midpoints between the shots and receivers (see Yilmaz, 1987). With 1200% data the C M P gathers had on average 12 traces in them. Because the traces came from approximately the same ground position, the velocity variations could now be properly compensated.  2.4.6  N M O Removal and Residual Statics Correction  The normal moveout (NMO) of the data results from the variation of receiver offset from the shots.  N M O is a function of reflection time, receiver offset and subsurface velocity. The  correction for N M O was done by analyzing the C M P gathers over a range of different velocities (different hyperbolic reflector trajectories) and picking the velocity that did the best job of flattening the reflectors. The N M O analysis was accomplished using the ITA semblance analysis routines; providing a graphical interface for picking the appropriate velocities. Due to the small number of traces in each gather, a series of adjacent gathers were stacked together to increase the signal to noise (valid as long as there was minimal velocity variation and structure among the stacked CMPs). The application of residual statics to the data was one of the key processes used to preserve the high frequencies output from the deconvolution. Due to small variations in velocity at each position along the line and the simplified technique used to correct for the elevation variations the traces within the CMPs had small-scale misalignments. If the data were stacked without correcting for these misalignments high frequencies would be lost in the final products 13  (destructive interference is wavelength dependent).  The residual statics algorithm computed  statics (trace shifts) that when applied to the data resulted in the correction of these misalignments.  The ITA algorithm is surface-consistent, meaning the statics computed can be  attributed to either specific shots or receivers (Yilmaz, 1987). The calculation of residual statics is best accomplished when the N M O has been completely removed from the data. However, the compensation for N M O is best done after the residual statics have been removed from the traces. Due to the circular, interconnected nature of this problem an iterative procedure was followed. After the first-pass N M O removal, the residual statics were calculated and then applied to the data before N M O removal.  The N M O was then reevaluated to establish the final velocity  function for the dataset and a final residual statics solution was computed and applied to provide the optimal compensation of the residual statics.  2.4.7  First-Break Muting and Trace Stacking  The final step before stacking involved the muting of the first breaks.  A relatively  conservative approach was followed, to preserve as much data as possible. The N M O corrected C M P gathers were examined and a mute function was picked that eliminated all severely stretched data and non-hyperbolic reflection energy. The C M P gathers were stacked to form the final stacked traces.  2.4.8  Post-Stack Noise Reduction  The interpretation of stacked seismic data is not always possible due to the presence of diffractions, noise and dipping events that are not properly located in space. To improve the interpretation all 3 lines were processed using noise reduction algorithms. The noise reduction was carried out using an f-x prediction algorithm, which removed non-coherent energy from the data. In the f-x domain coherency is manifested by trace-to-trace similarities of amplitude and 14  phase that can be predicted by a relatively short operator. Incoherent noise, which tends to have a random response in f-x space cannot be predicted by the same short operator. Therefore, the fx noise reduction reduced the noise without disturbing the coherent signals. In order to ensure coherent energy was not removed a difference section between the starting and final data was output and showed only random, incoherent energy.  2.4.9  Migration  The final process applied to the data was migration.  The migration process collapses  diffractions and places dipping events in their proper subsurface location (Yilmaz, 1987). A n f-k based migration algorithm was used that performs the migration in the frequency-wavenumber domain. The velocity function chosen for the migration was based on a smoothed, edited version of the stacking velocities. The stacking velocities were chosen from a portion of the line with good signal-to-noise and little geologic dip. A near-surface velocity was assigned based on the types of sediments apparent in the ditches where the surveys were carried out (typically around 1600 m/s). The final migrated sections showed the complete collapse of diffraction energy and no artifacts that could be attributed to an improper velocity function. The final stacked and migrated seismic sections are shown in Figures 2.4, 2.5 and 2.6 corresponding to lines 20-206, 216-28 and 32-184 respectively.  The stacked and migrated  datasets of each line are plotted at the same horizontal scale to allow for comparisons between the two products. The sections are plotted in time with the maximum time varying between 0.20 and 0.23 s. The horizontal scale changes from section to section because the number of traces on each section varies. The scale was chosen to fit the seismic section to the width of the page, without distorting the data. Good reflections are obvious on each of the seismic datasets and will be discussed in greater detail in the interpretation section.  15  2.4.10  Amplitude Spectra  The initial goal in the acquisition of the P-wave seismic data was to obtain as high a resolution as possible. The resolution is a function of wavelength (approximately 1/3 according to Miller et al. (1995)), which is a function of velocity and frequency. In order to get some idea of the final resolution of each of the 3 seismic lines an average amplitude spectrum was computed for each line and is shown in Figure 2.7.  A total of 71 traces were used in the  computation of the spectra over a time zone of 0.01 to 0.15 s from the migrated seismic datasets. Using a cutoff amplitude of-12dB the bandwidth of each of the lines was determined to be: 200700 Hz for line 20-206, 200-600 Hz for line 216-28 and 200-1100 Hz for line 32-184. The dominant frequency is probably the best indication of resolution and was also determined for each line: 400 Hz for line 20-206, 300 Hz for line 216-28 and 500 Hz for line 32-184. This gives a minimum bed resolution for each of the lines (assuming a velocity of 1700 m/s and using the 1/3 wavelength criteria): 1.4 m for line 20-206, 1.9 m for line 216-28 and 1.1 m for line 32-184.  2.5  Seismic Interpretation  The objective of the seismic interpretation was to obtain as much information as possible about the subsurface aquifers apparent on each of the seismic lines. A reliable interpretation of the data requires a regional context that cannot possibly be obtained by the examination of 600 m of seismic data.  As a result, the interpretation relied on the work of Pullan et al. (2000) to  identify geological boundaries based on other seismic lines and subsurface control (water wells and gas exploration wells).  M y interpretation of these lines will build upon the existing  interpretation of Pullan et al. (2000), incorporating nearby water well data whenever possible. The interpretation of potential aquifers on the seismic lines relied on a relatively simple relationship between the reflectivity in the seismic data and the presence of coarse-grained sediments.  It is assumed that the highly reflective zones are associated with coarse-grained 16  sediments and the zones showing very little reflectivity are associated with fine-grained sediments. The environment of deposition created sedimentary features with characteristic patterns that were seen in the seismic data and used in the interpretation.  The fine-grained  sediments were deposited in a marine setting (probably during high sea level, early interglacial periods) and the coarse-grained sediments were deposited in a near-shore or glaciofluvial environment while the sea level was lower and glaciers were actively advancing. The highenergy environments were characterized by incising channels and sloping depositional surfaces as large quantities of glacial sediments were being transported and accommodated.  Due to the  limited view of the subsurface stratigraphy provided by 3 short seismic lines, the model may be oversimplified but it serves as a first-pass interpretation. Differentiating between fine-grained and coarse-grained sediments does not directly identify potential aquifers. Many of the coarse-grained sediments identified in the water well logs are clay-rich and not suitable as aquifers. The best aquifers are in clean sand and gravel zones, typically adjacent to clay-rich coarse-grained sediments. Interpreting coarse-grained sediments based on large-scale seismic reflection patterns requires the additional interpretation of smallscale reflection patterns to identify the aquifer zones within them. As noted by Ricketts, (2000), aquifers can be identified on seismic data by the identification of distinct stratigraphic sequences. Where an unconformity has fine-grained sediments above and coarse-grained sediments below it indicates that a sea level transgression has occurred. If the transgression was associated with an active shoreline depositing and reworking clean sand, then a clean sand aquifer may now exist directly above the unconformity. The final interpreted seismic datasets were converted from time to depth.  With the data  plotted in terms of depth and with the horizontal scale equal to the vertical scale (1:2000) the dipping events are representative of apparent dip, thereby creating the best possible plots for a detailed geological interpretation. The interval velocity functions chosen for the time to depth 17  conversion were computed by applying the Dix equation (Dix, 1955) to the stacking velocity functions used to migrate the data. Depths are estimated to be accurate to within 10% of the true depth down to a depth of 30 m and possibly much less accurate at deeper depths (where the limited offsets resulted in poorly constrained N M O velocities). Water well records were used to assist with the seismic interpretation. These water wells were unequally distributed on each of the seismic lines, with locations roughly corresponding to the locations of houses along the roadways. The water wells have been assigned unique well numbers (links) by Ricketts and Makepeace (2000) and I have used the same numbers to identify the wells (I have also included the street addresses). The descriptions of the subsurface materials varied among the well records because the wells were drilled by a variety of individuals and contractors that employed their own standards for describing the subsurface sediments.  I've  attempted to make the interpretations as uniform as possible and developed a simple grayscale, symbolic graphical scheme to represent the different sediment types.  2.5.1  Interpretation of Line 20-206  The interpreted version of line 20-206 is shown in Figure 2.8. The stacked seismic section was chosen instead of the migrated section for the interpretation because the migration process seemed to destroy some of the prominent reflections (see Figure 2.4a and 2.4b at a time of 0.15 s). The precise reason for this effect is unknown but is probably related to the limited migration aperture of the deeper zones of the dataset due to the short length of the line. As shown by Figure 2.8 (with no vertical exaggeration) the plot has a greater maximum depth than the length of the line making it difficult to properly migrate dipping events greater than 100 m in depth. Line 20-206 was acquired in a region that has been interpreted to correspond to the paleoshoreline at the time the Brookswood delta was deposited (Roberts et al., 2000). The eastern half of the line is onshore and the western half represents distributary channel deposition. 18  In the middle at approximately 50 m is the transition zone, which was imaged very clearly on GPR data (Rea and Knight, 2000). There were no usable data in the upper 20 m of the seismic line, thereby preventing the imaging of the same near-surface transition zone.  The most  prominent reflective zones are at depths of 120 m and 160 m. The zone at 160 m was interpreted by Pullan et al. (2000) as possibly corresponding to the top of the Tertiary bedrock.  This  interpretation is supported by the westward dip of the interface. Water wells approximately 1800 m east and 1600 m north define the Tertiary interface at a depth of 100 m. The dipping interface at 170 m on Line 20-206 would thus correspond to a fall-off from the bedrock high to the east.  2.5.1.1  Nearby water well records  The water wells closest to line 20-206 do not penetrate very deeply so a deeper water well (link 959) was tied into the line and is shown in Figure 2.8. The water well location is scaled properly horizontally, but also lies off-line approximately 150 m to the north. A number of other deeper water wells approximately the same distance from the line also show similar stratigraphy so it is assumed this water well is a good representation of the sediment profile in the upper 40 m. The water-bearing zone in the well is at its base and consists of clean sand overlain by a tight, clay-rich sand zone. Projecting the top of the water-bearing zone onto the seismic seems to correlate with a seismic reflection at a depth of 30 m, which can be traced to a horizontal distance of 70 m on the line. The signal quality in this upper zone is very poor so the base of the water-bearing sand is uncertain but may be represented by the reflection at a depth of 40 m.  2.5.1.2  Identification of a deep, undrilled aquifer  Additional aquifer zones on the line are probably related to the reflective zone at a depth of 120 m, the top of which is likely an unconformity surface similar to the one identified by Ricketts (2000). One of the potential aquifers would thus lie on top of this unconformity surface 19  within clean, shoreline sand that cannot be separately resolved on the seismic line. The reflective zone beneath the unconformity appears to be thicker to the west possibly indicating the presence of a buried channel with another potential aquifer. In general, line 20-206 represents the poorest quality seismic data acquired on the Brookswood. As a result, it is more challenging to define the aquifer zones.  2.5.2  Interpretation of Line 216-28  The interpreted version of line 216-28 is shown in Figure 2.9. This line is east of the Brookswood aquifer with Fort Langley glaciomarine sediments at the surface. The base of the fine-grained sediments is clearly apparent on the seismic line as a reflection at a depth of 40 m and has been interpreted as an unconformity surface.  2.5.2.1  Nearby water well records  There are two water wells (links 1644 and 1612) that are located reasonably close to the line (as shown in Figure 2.9). Both wells are completed in coarse-grained sediments close to the interpreted unconformity surface. The good, water bearing sand and gravel at the base of well 1644 may be an example of a shoreline aquifer but well 1612 shows the sand and gravel zone at the same level as being hard-packed (presumably tight) and produces from a deeper sand and gravel zone. The wells are separated from each other by approximately 100 m, with the seismic line between them, so a definitive interpretation is difficult.  The unconformity surface is  reasonably easy to correlate to nearby water wells (links 1613, 1618 and 1638) and appears to drop-off to the north (keeping in mind that well 1638 is not located properly on the section).  20  2.5.2.2  Identification of discrete aquifers within a glacial till  The dominant reflections on the line are contained within a zone extending from 40 m to 120 m.  Pullan et al. (2000) suggest the reflection at 110 m corresponds to the Tertiary bedrock,  which ties the bedrock identified at the base of water well 1638 (400 m to the north on 216  th  Avenue). The well records show a sequence of glacial till with interbedded sand and gravels beneath the unconformity surface and extending to the final depths. The interbedded sand seems to correspond with the most prospective aquifer zones and may also correspond with the coherent reflective zones within the coarse-grained deposits interpreted on the seismic data. This fact strongly supports the interpretation of these dipping zones as continuous sand aquifers and demonstrates how these aquifers can be mapped to the north and south. There is not enough information in the well records to indicate whether the glacial till falling between the sandy aquifer zones also may have a reasonable hydraulic conductivity. It is definitely noted to be tilllike suggesting the presence of larger grains, possibly gravel, but the intergranular matrix also has a large amount of fine-grained, possibly clay material suggesting a very low permeability.  2.5.3  Interpretation of Line 32-184  The interpreted version of line 32-184 is shown in Figure 2.10. This line is west of the Brookswood aquifer with Capilano glaciomarine sediments at the surface.  The Capilano  sediments most likely extend to the strong reflection at a depth of approximately 55 m (Pullan et al., 2000). The reflection at a depth of 200 m was interpreted by Pullan et al. (2000) as possibly the top of the Tertiary bedrock.  2.5.3.1  Nearby water well records  The water well records on this line are sparse and the ones that are available have not differentiated the sediments. Water well 1120 classifies all sediments as "glacial deposits and 21  alluvium", suggesting a lot of similarity of the sediments with no major lithology changes. There is also no indication as to the depth where potable water is being extracted and I am assuming that a fairly large interval was screened and production is out of a series of discrete sand bodies.  Water well 1130 provides a very descriptive record but is also located  approximately 300 m to the east (off the eastern end of the line). The uppermost sand and gravel in well 1130 is the Brookswood aquifer extending to an elevation of 19 m (a depth of 1 m below the seismic datum). The lower part of the well consists of interbedded clay and water-bearing silt, sand and gravel layers. This same zone seems to be laterally consistent with a series of high amplitude reflections on line 32-184 and are interpreted to correspond to the water-bearing zones of well 1120.  2.5.3.2  Identification of thinly-bedded sand aquifers  Unlike line 216-28 where the near surface sediments were uniform down to their base, on line 32-184 the fine-grained marine sediments are interbedded with a series of thin coarsegrained deposits (as seen in water well 1130). A well-formed channel at a depth of 25 m and a distance of 150 m suggests that these sediments may not have been deposited in a deepwater, marine setting.  The channel is more characteristic of a fluvial environment but might also  represent a submarine or nearshore distributary channel.  The presence of a series of high  amplitude reflections forming a package extending from a depth of 15 to 30 m also supports this hypothesis.  These high amplitude reflections cannot possibly represent uniform fine-grained  marine sediments and my interpretation is they represent a series of nearshore continuous sand bodies probably very close to a paleoshoreline.  22  2.5.3.3  Identification of a high-quality sand and gravel aquifer  Water well 1130 produces water from a coarse water-bearing sand and gravel at its base. This zone correlates with a wedge of material that is directly on top of the interpreted unconformity surface. The wedge appears to infill a channel on the unconformity at a distance of 280 m on the seismic line and the wedge pinches out at a distance of 190 m. This wedge is interpreted as a very high quality aquifer. It has impermeable sediments above it and to the west and possesses composition.  an anomalously low reflective quality suggesting a very homogeneous It has a thickness of close to 10 m and one lateral dimension of possibly as great  as 500 m (assuming the correlation with well 1130). It was most likely deposited as a shoreline deposit during a transgressive sea-level rise. It appears to be separated from the unconfmed Brookswood aquifer, thereby protecting it from near-surface contaminants.  2.5.3.4  Identification of deep, undrilled aquifers  Additional aquifers on line 32-184 are most likely present below the upper unconformity surface. The series of high amplitude reflections with a depth range of 60 to 70 m suggest a series of interbedded sand and gravel deposits (interpreted as Vashon Drift by Pullan et al. (2000)). The entire zone between the two unconformities has excellent aquifer potential. If the large channel was incised and infilled during the same interglacial period, there is a good chance the entire sequence consists of sands, gravel and till. On the other hand, i f the deeper channel was part of a large glacial lake, the lower part of the channel could be filled in with fine-grained, glaciolacustrine silt. The precise interpretation of the data requires additional control. Another geophysical technique such as EM-47 could be used to determine if the channel corresponds with a low conductivity zone (as demonstrated for similar sedimentary sequences by Best and Todd, 2000) or a well must be drilled to determine the true quality of the subsurface aquifer.  23  2.5.4  Summary of Seismic Interpretations  Each seismic line has provided some insight as'to the location of additional aquifers. The lines also show how quickly the stratigraphy changes, with deeply incised channels and sipping depositional surfaces that have no obvious surface expression. Some of the wells (e.g. link 1120) appear to have stopped just short of an excellent aquifer zone. The seismic identifies the major unconformity surfaces and reflective zones most likely to be associated with good aquifers and the future application of seismic in the area would ensure that the best aquifers are located and fully utilized.  2.6  Conclusions The use of P-wave seismic for aquifer delineation seems to be most effective when the  surficial sediments are fine-grained and saturated.  These types of surficial sediments provide  excellent source and receiver coupling resulting in the acquisition of very high-quality seismic reflection data.  From a groundwater resource perspective the areas with fine-grained, possibly  clay at the surface are also the most desirable places to use P-wave reflection seismic; in order to produce groundwater the thickness of the clay zone must be determined. The seismic reflection data not only determines the thickness of the near-surface clays but also provides an indication of the thickness of aquifer-hosting coarse-grained sediments and the continuity and structure of those sediments. On a regional scale, P-wave seismic could be used quite effectively to map the deep aquifers. This would require a large survey effort involving the acquisition of 10s of km of seismic reflection data in a grid pattern to ensure the identification of three-dimensional structure. The work of Pullan et al. (2000) serves as the start of this process in the Fraser lowland. Recognition of the importance of the Fraser Valley groundwater resources by the government would be required to obtain the necessary funding for an expanded seismic program. 24  236 30' 49' 30  237° 00'  237' 30' 49" 30'  British i,, C o l u m b i a Pacific  49' 20' H  h 49' 20'  Ocean \ J W^ington^  49" 10'  Strait of Georgia  h  49' 00' 236' 30  49' 00' 237' 30'  237* 00  o — CN  40th A v e .  • ' .CN  32nd A v e ,•• .|216-28||j  24th A v e .  28th A v e  20th A v e . 16th A v e  8th A v e . Brookswood Aquifer _i  F i g u r e 2.1:  2 km I  Location map showing the outline of the Brookswood aquifer and the location of the three seismic lines: 20-206, 216-28 and 32-184.  25  49' 10'  Figure 2.3: The shot profiles from Figure 2.2 with a bandpass filter applied (and 0.04s AGC). 26  0  25  50  75  100  Distance West of 206th Street (m)  0  25  50  75  100  Distance West of 206th Street (m) Figure 2.4: Processed seismic data for line 20-206; a) stack b) migration 27  0  30  60  90  120  150  Distance South of 28th Avenue (m)  a  Figure 2.5: Processed seismic data for line 216-28; a) stack b) migration. 28  180  0  50  100  150  200  250  300  250  300  Distance East of 184th Street (m)  0  50  100  150  200  Distance East of 184th Street (m) b  Figure 2.6: Processed seismic data for line 32-184; a) stack b) migration. 29  a  500  1000  1500  2000  Frequency (Hz)  f J-io  ^ -  -  \  j  -  V  \  -J <: -20 b  > CD  -  —\  I -30 -35  500  1000  1500  2000  1500  2000  Frequency (Hz)  1000  Frequency (Hz)  Figure 2.7: Average amplitude spectra of the migrated seismic data; a) 20-206 b) 216-28 c) 32-184.  30  31  32  33  1  CHAPTER 3 *  Near-Surface VSP Surveys Using the Seismic Cone Penetrometer  3.1  Introduction Vertical seismic profile (VSP) surveys are a means of improving the use of seismic data in  obtaining information about the subsurface.  In a standard V S P survey, a source at the surface  generates seismic waves and receivers are deployed in an existing borehole.  The receivers,  which can be either downhole hydrophones or clamped geophones, are raised from the bottom of the borehole and stopped at regular intervals to record the seismic data. In most cases, the main benefit in conducting a V S P survey is the improved interpretation of surface common depth point (CDP) seismic data. The velocity control provided by the V S P survey allows improved accuracy in the conversion of time sections to depth sections; multiples and primary reflections can be identified and deconvolution operators designed to eliminate the multiples. In addition, the increased bandwidth of the VSP data, compared to surface reflection data, often makes it possible to resolve smaller scale features. VSP surveys have been extensively used in the petroleum industry (see, e.g., Gal'perin, 1974; Wuenshcel, 1976; Hardage, 1985) and have been shown recently to be a useful technique for near-surface investigations (Skvortsov et al., 1992; Milligan et al., 1997). The Milligen et al. survey was carried out in recent deltaic sediments (similar to the sediments in our investigation) consisting of fine silts interbedded with sand and gravel lenses. Downhole hydrophones were used to acquire V S P data over a depth range of 3 to 28 m. A multi-offset technique was used that resulted in a profile of seismic reflection data extending 9 m away from an existing borehole * A modified version of this chapter is published as: Jarvis, K . D . and Knight, R.J., 2000, Near-surface VSP surveys using the seismic cone penetrometer: Geophysics, 65, 1048-1056. Used with permission. Copyright by the Society of Exploration Geophysicists.  34  and with a depth range of 3 to 50 m. The bandwidth of their VSP data is significantly greater than the surface P-wave data acquired at the same site with usable frequencies as high as 900 Hz. The use of a 24-level hydrophone array greatly reduced the acquisition time. However, tube waves were a significant problem and were only partially eliminated by baffling the array. The development of a near-surface V S P technique, which does not rely on a borehole would eliminate any problems with tube waves and result in higher quality VSP data. The use of VSP surveys for near-surface investigations has been relatively limited. One of the reasons for this limited use is the need for a borehole to acquire the data.  The drilling of a  borehole significantly increases the costs associated with geophysical site characterization. In addition, there are often environmental or safety considerations such that the drilling of a borehole is not allowed. Given the significant potential value of VSP data from a site, we have investigated the acquisition of a VSP survey using a cone penetrometer; this has the advantages of reduced cost, and the cone penetrometer is generally considered to be "minimally invasive". A n added benefit is the simultaneous measurement of other geotechnical parameters such as penetration resistance and friction ratio, which can then be used for stratigraphic correlations. We demonstrate the technique using a single offset S-wave source.  However, the same  procedure is also applicable for the acquisition of P-wave and multi-offset VSPs.  3.2  Geologic Setting This investigation was carried out on the Fraser River delta, southwestern British Columbia,  Canada at the Kidd2 research site. The location is shown in Figure 3.1. This site has been the focus of a number of studies including the Canadian Liquefaction Experiment (CANLEX), which resulted in a large number of cone penetrometer holes and the retrieval of three frozen cores for laboratory measurements (Hofmann, 1997). A number of hydrogeologic experiments have also been carried out at the site by using a pumping well, multi-level sampling wells and 35  zone-specific piezometer installations. Neilson-Welch (1999) used the multi-level sampling wells to obtain groundwater chemistry data, which were combined with measured hydrogeologic parameters to develop a groundwater flow model.  The Geological Survey of Canada has  obtained core, participated in geotechnical studies and performed geophysical surveys at the site (Hunter et al., 1998b) in an ongoing effort to study the liquefaction potential and earthquake site amplification effects of the Fraser Delta. Based on information from these previous studies, we have a simple stfatigraphic model involving the major lithologic units at the site: there exists a near-surface silt layer, 4 to 5 m thick; a sand-dominated unit, 17 to 18 m thick; a clayey-silt layer, 26 to 28 m thick; and a compact Pleistocene till at a depth of approximately 50 m. Due to the cooperation of B C Hydro, Kidd2 has become an ideal field research site where it is possible to carry out both hydrogeological and geophysical investigations. The inset map in Figure 3.1 shows the acquisition grid of SH-wave CDP reflection data at the site with the two cone VSPs falling on one of the lines. The water table is approximately 1 m below the surface with some variation due to tidal influences. Due to the nearby estuary, a salt water wedge is present within the aquifer with a depth range of approximately 12 to 22 m.  3.3  Cone Penetrometer Data The cone penetrometer is a technology developed by geotechnical engineers for in-situ site  characterization. The basic idea is to push a steel rod, with a cone-shaped tip, into the ground while making measurements with sensors mounted close to the tip. The sensors, which typically consist of pressure transducers, inclinometers, thermistors and accelerometers or geophones, are used to measure a variety of parameters such as tip resistance, sleeve friction, temperature, pore pressure and first arrival shear-wave velocity (Campanella et al., 1983, Robertson et al., 1986). The cone is pushed into the ground using hydraulic rams with the pushing force offset by the weight of the truck (typically 5 to 10 tons). The technique works extremely well in areas of 36  alluvial and deltaic deposits consisting of clays, silts or sands. It is not well-suited for use in gravelly sediments or rocks.  3.3.1  V S P Data Acquisition  Two zero-offset shear-wave VSP surveys were conducted using the cone penetrometer truck from the Department of Civil Engineering at the University of British Columbia (UBC). A large sledgehammer struck against the baseplate of the truck was used as a source of SH-waves and a single cone-mounted accelerometer used as a receiver. Two records were collected with hammer blows in one direction, then two records were collected with hammer blows in the opposite direction. The records from impacts in the same direction were stacked and subtracted from each other to enhance the signal to noise ratio of the SH-waves and reduce P-wave interference. The data were recorded using a Nicolet oscilloscope with 15 bit data digitizing and storage capability. The depth levels at which the cone was stopped to record V S P data were chosen as a compromise between field time and spatial aliasing. The truck engine (which must be on to run the hydraulic system) was turned off before every measurement to reduce the noise. Two V S P surveys were recorded. The first survey was conducted using a constant depth interval of 0.5 m, from 3.0 m to 46.5 m, for a total of 88 traces.  The second survey was  performed using variable depth intervals: 0.33 m from 2.0 m to 10.0 m, 0.50 m from 10.5 m to 25.0 m and 1.0 m from 26.0 to 49.0 m, for a total of 80 traces. The variation in depth intervals on the second survey was used as a means of reducing field acquisition time, and obtaining increased coverage of the shallowest reflections. For the VSP experiments no modifications were made to existing instrumentation except to extend the record length. With standard cone penetrometer surveys, first arrivals are recorded to determine interval velocity, and record lengths are commonly less than 0.5 s. For the V S P survey, the record length was extended to 2 s in order to record the later arrivals. Saving data on 37  the Nicolet oscilloscope was a time consuming part of the operation.  A digital acquisition  system, using a PC-based acquisition card and software, would significantly reduce acquisition times. A n important aspect of first arrival measurements is the time and depth accuracy. Time zero, which corresponds to the instant the baseplate is struck, must be known with a very high degree of accuracy. The U B C truck is equipped with an electrical contact trigger, detailed by Robertson et al. (1986).  When the metal head of the sledgehammer contacts the metal baseplate, an  electrical circuit is closed which triggers the oscilloscope. Depth accuracy also determines the accuracy of calculated interval velocities. The cone truck uses 1.000-m-rod lengths, resulting in extremely accurate depths.  The estimated time accuracy is +/- 0.03 ms and the estimate depth  accuracy is +/- 0.001 m.  3.3.2  V S P Data Processing  Data from the first cone V S P survey (referred to as VSP1) are shown in Figure 3.2. Data from the second cone V S P survey (VSP2) are shown in Figure 3.3. The first survey has two regions with a lower signal-to-noise ratio at depths of 18.0 and 31.0 m; these are attributed to sporadic instrumentation problems. The upgoing reflection from the top of the Pleistocene till is at a time of approximately 0.5 s at a depth of 3 m and 0.25 s at a depth of 46.5 m. A downgoing free surface multiple from the Pleistocene till is parallel to the first breaks with a delay of approximately 0.5 s after the first break times. The variation in trace spacing on the second survey (Figure 3.3) is due to the changing depth increment previously discussed. The upgoing Pleistocene till reflection can be observed in Figure 3.3; it is obscured at depths less than 7.0 m and most obvious between 7 and 25 m at a time of approximately 0.4 s. A downgoing free surface multiple from the Pleistocene till is not obvious on this dataset. Upgoing reflections  38  above the Pleistocene till are not obvious in Figure 3.2 and Figure 3.3 due to the high-amplitude downgoing wavefield. Even though both surveys were acquired at the same site, separated by only 34 m, there are considerable differences in data quality. VSP2 has a series of near surface multiples, which are not apparent in the data from VSP1. This is most obvious on the traces in the upper 7 m in Figure 3.3, at times greater than 0.4 s, but the effect is also noticeable in the character of the first arrival, which is more compact on VSP1 (Figure 3.2).  These multiples are attributed to a  variation in the near surface conditions. The surveys were done at different times of the year (VSP1 in the fall and VSP2 in the spring) so the water content of the near surface sediments is likely to have varied. Changes in the saturation conditions of the near-surface have been shown to have significant effects on the reflection character and quality of near-surface reflection data (Jefferson and Steeples, 1995) and undoubtedly would have similar effects on V S P data. The study area also has a heterogeneous near-surface with variations on the scale of meters due to the development of the site for an electrical sub-station.  The fill layer in the upper 2 m varies  throughout the site from hard compacted till to loosely compacted sand, silt and organic material. This likely resulted in significant differences in source and receiver coupling which thus affected the VSP data. The differential arrival times from the direct wave first breaks were used to generate detailed interval velocity with depth functions shown in Figure 3.4.  The determination of interval  velocities using the V S P technique was an important objective of the survey. For most nearsurface studies borehole sonic logs are not available, and the VSP-generated velocity profiles become critical for subsurface depth control. The estimated error in the interval velocities derived from these VSP surveys is 5 percent. The processing of the two VSP datasets followed a generally accepted processing flow (as outlined by Hardage, 1985). The main objective of the processing was the separation of the 39  upgoing and downgoing wavefields, with improvement of signal bandwidth using deconvolution. A processing flow chart is shown in Figure 3.5. Compensating for spherical divergence involves a simple amplitude correction. Interval velocities from Figure 3.4 were used to.compute a predicted amplitude decay determined by the velocity-squared time method [(V T)'](Hardage, 1985), where V is the root-mean-square (RMS) 2  velocity and T is the total traveltime to a particular depth. The slope of the least-squares fit line was 1.55. A n amplitude correction was applied to the data by multiplying each amplitude value byT  155  .  To isolate the downgoing wavefield the singular value decomposition (SVD) technique outlined by Freire and Ulrych (1988) was used. deconvolution and phase adjustment.  The downgoing wavefield was used for  Minimum phase deconvolution operators were designed  based on the SVD-filtered data and applied to both the original data and the SVD-filtered data. The deconvolved SVD-filtered data were used to design match filters that were applied to the deconvolved original data. A zero-phase wavelet (suggested by Berkhout, 1984) with frequency cutoffs of 15 and 150 Hz was chosen for the match. The match filter had the effect of removing trace-to-trace wavelet variations and reshaping the embedded wavelet to a zero phase wavelet, to allow for easy interpretation of the final results. An f-k filter was applied to the deconvolved, match filtered data. In order to balance the trace amplitudes for the f-k filter and avoid edge effects the direct arrivals were muted on all traces and a linear taper  was applied to the top, bottom and sides of the data.  The compact first break on the data after deconvolution and match filtering did not require the application of a significant mute, thereby preserving as many data as possible. A 2-D Fourier transform of the data revealed the f-k spectrum and a pie shaped filter (with edge tapers) was designed to isolate the aligned downgoing wavefield. The reverse transformed data were then  40  subtracted from the input data, to remove the downgoing wavefield without affecting the upgoing wavefield or background noise. In order to flatten the upgoing wavefield, traces were shifted downwards in time by an amount equal to twice the first break arrival times. A top mute was designed to eliminate spurious, high-amplitude f-k-filter artifacts, which tended to be generated near the initial first break times. A bottom mute was designed to remove data dominated by residual multiples and waveform changes due to attenuation (Hardage, 1985). The application of this mute is analogous to the generation of a corridor stack. A final step involved stacking the muted data to create the VSP extracted traces (VETs).  The muted upgoing wavefield and VETs for the two surveys,  referred to as VET1 and VET2, are illustrated in Figure 3.6 and Figure 3.7. In Figure 3.7, data below 25 m were not used to generate the V E T for the second survey because the deeper data were aliased. The f-k spectra (not illustrated) clearly showed that data at a sampling interval of 1.0 m were spatially aliased at frequencies greater than 60 Hz, whereas the maximum usable frequency appeared to be around 130 Hz. Deep data at a depth-sampling interval of 1.0 m were therefore not used to construct VET2. The elimination of the deeper data and the successful generation of VET2 illustrates that the cone does not necessarily have to penetrate the complete sediment column in order to obtain useful VSP information.  3.3.3  Cone Data Acquisition  The interval velocity measurements obtained using the cone penetrometer can only be used for a limited geological interpretation. This is partly due to the coarse depth sampling and also partly due to the error in the measurements, which is close to the total change in velocity associated with the different sediment types. As a result, some other independent information should be used for precise geologic calibration. This information is obtained by the other sensors  41  installed in the cone that provide 2.5 cm sampling and are very sensitive to variations in the sediment types. While the cone was being.pushed from one V S P depth level.to a deeper level the tip resistance, sleeve friction, and pore pressure were measured. The tip resistance is a measure of the penetration force and is low in fine-grained sediments such as clay or silt and high in coarser grained material such as sand (Robertson and Campanella, 1983).  The sleeve friction is a  measure of the sliding friction of the cone. Geotechnical engineers have determined that the percent ratio of sleeve friction to tip resistance (referred to as the friction ratio) is a very sensitive indicator of sediment type. The friction ratio has a low, flat response in sands and a higher reading in silts and clays (Robertson and Campanella, 1983).  Penetration pore pressure is  relatively low in sands (hydrostatic) and higher when clay minerals are present. When the cone is stopped, dissipation of pore pressure with time can be observed. The dissipation rate aids in sediment classification and can be used to interpret the consolidation characteristics of the material (Campanella et al., 1983). Typically, measurements are made at 5 cm intervals and the combination of tip resistance, friction ratio and pore pressure are used to identify distinct stratigraphic sequences. These sequences are associated with interfaces which should correlate with reflections on VSP surveys.  3.4  SH-Wave Seismic Reflection Data  3.4.1  Data Acquisition  The C D P shear-wave seismic data were collected in June of 1996 using a 24 channel Geometries SmartSeis seismograph, 48 Geo Space horizontal geophones (20DM-28 Hz) attached to CDP spread cables and an I/O RLS-100 roll switch. The data sample interval was 0.0005 s and an anti-alias, high-cut filter was applied with a corner frequency of 600 Hz. A n end-on  42  shooting configuration was used with a geophone and source spacing of 1 m, giving a far offset of 24 m and 1200% subsurface coverage. The seismic source consisted of a 35 kg steel block with metal fins attached to the base to increase the source coupling. SH-waves were generated by swinging a 5 lb sledgehammer against the sides of the block. To further increase source coupling the block was stood upon by the person swinging the hammer. The SH-wave technique was utilized to minimize mode conversions and P-wave interference. A n average of 10 hammer blows were stacked for each record to improve the signal to noise ratio and remove the effects of the occasional misdirected blow. Data were recorded with blows in opposite directions to allow for the subtraction of the records to further reduce P-wave interference and provide additional signal to noise enhancement of the SH-waves. A total of five lines of SH-wave seismic data were acquired and processed, covering 420 m. The locations of these lines are shown in Figure 3.1. This paper will only examine line 20, a 100 m north-south line (highlighted in Figure 3.1).  3.4.2  Data Processing  The seismic data were processed using ITA's UNIX based seismic processing system installed on a Sun workstation network. The processing flow consisted of: spherical divergence correction, bandpass filtering, record subtraction, deconvolution, normal moveout removal, residual statics, stacking and migration. Refraction statics were not applied to the data because the elevation at the site varied by less than 1 m. Most of the large velocity variations in the nearsurface should be confined to the upper 2 m, which can be dealt with adequately by the residual statics. A spherical divergence correction was applied using a time scaling function with a power of 2.0. filter.  The data were bandpass filtered (25/50-150/400 Hz) using a minimum phase Butterworth The record subtraction process was preceded by the application of a trace to trace 43  amplitude balance. The surface consistent deconvolution used both shot and receiver gathers, with trace by trace editing required afterwards to remove the noisy trace segments.  The  elevation statics compensated for a topographic variation of less than 1 m across the site. A semblance velocity analysis was used iteratively with surface consistent residual statics. The mute function tended to be quite conservative to ensure that the amplitude variations at shallow depths did not dominate the stack. The data were migrated using an interval velocity model based on the extensive cone data available at the site, which defines the depth of key reflecting horizons. Two shot profiles from both ends of the line are shown in Figure 3.8. represent the data quality after the records have been subtracted.  These profiles  Hyperbolic reflections are  obvious at a time of 0.15 s and a reflection from the top of the Pleistocene till is present at a time of 0.5 s.  Additional reflections cannot be seen at this stage of the processing; the stacking  process is needed for enhancement.  3.5  A Comparison of the Geotechnical Logs and VSP Data We obtained high quality VSP data over a region of the subsurface where we also have cone  penetrometer measurements of geotechnical parameters. The two types of data were compared to determine how each technique responds to changes in lithology and near-surface conditions. Figure 3.9 and Figure 3.10 are composite plots, which compare the geotechnical logs: [friction ratio (FR), tip resistance (Qc) and pore pressure (U2)] with the shear-wave velocity (Vs) and the VSP data (VET1 and VET2). The interpretation of the data involves a number of factors related to the complex nature of near-surface, unconsolidated sediments or soils. One of the dominant influences on the SH-wave velocity and tip resistance (Qc) is the rapid increase in effective stress in the upper 20 m that is essentially a compaction effect. As the sediment load on a uniform sand layer increases with depth, the porosity decreases and the 44  material becomes stiffer, thereby making it harder to penetrate and increasing the shear modulus. In Figure 3.9 the overall trend of increasing tip resistance from 5 M P a at 4 m to 15 M P a at 23 m is an example of the increase in penetration resistance due to compaction. In Figure 3.10 the velocity trend from 150 m/s at a depth of 5 m to 220 m/s at a depth of 22 m is an example of the increasing shear modulus due to compaction. Smaller scale variations within these dominant trends are possibly related to non-uniformities within the layers. In comparing the geotechnical logs to the VSP data, it is important to recognize that the logs are sensitive to a variety of physical properties, which are not necessarily related to seismic velocity or density. As a result, it is possible to observe significant changes in the geotechnical logs with no corresponding seismic reflection at that depth on the VETs. For example, the tip resistance shows large variations due to changes in sand grain size; this variation in grain size does not necessarily cause a variation in seismic velocity or density. Civil engineers (e.g. Lunne et al., 1997) have attempted to model the cone responses in terms of specific soil conditions such as density, clay content and grain size but have had little success. Interpretation of the cone data is therefore accomplished using empirical relationships. Based on repeated measurements at different sites, a crossplot of friction ratio and tip resistance is used to classify soils (Campanella and Robertson, 1983). Another complexity is the fact that the cone does not fully respond to "thin" layers (Campanella and Robertson, 1983). For the standard cone used, the tip resistance will not reach a full response unless the layer is at least 36 cm thick (Campanella and Robertson, 1983). Due to these inherent complexities, we cannot expect there to be a simple relationship between the cone data and the seismic data. The interpretation of the geotechnical logs and the V S P data is based on the stratigraphic model available for the site: a near-surface layer of silt, underlain by a sand-dominated unit, underlain by clayey-silt, with the compact Pleistocene till at the base. A fine-grained silt layer 45  near the surface (upper 4 to 5 m) is identified by the high friction ratio and low tip resistance. The shear-wave velocity increases rapidly in this interval from 60 m/s to 150 m/s (Figure 3.10). This rapid increase in velocity is most likely a compaction effect due to the high gradient in nearsurface effective stress.  There are no V S P data in this depth range because data were not  recorded until the cone was near the base of the layer. The zone from 5 m to approximately 23 m is a sand-dominated unit. In this depth interval, the friction ratio remains low while the tip resistance shows an overall increase with depth. This can be seen in the geotechnical data from both surveys.  A significant variation in the tip  resistance (most sensitive to sand grain size) is observed in both data sets in this depth interval. In this unit we can interpret two fining upward sequences (approximately 5 m to 13 m and 13 m to 22 m) within which the tip resistance decreases with decreasing depth. These two fining upward sequences can be correlated throughout the Fraser Delta (Monahan et al., 1993). The seismic propagation velocity increases within this sand layer. Most of the increase is due to the rapidly changing effective stress. There also appears to be a correlation between the zones with the highest tip resistance and those with slightly higher velocity. Reflection events on the VETs are also identified within the sand unit. In Figure 3.9 VET1 shows 3 prominent reflections between 250 and 300 ms that are related to the lower fining upwards sequence. The base of the sequence at a time of 280 ms correlates with a trough on VET1 and a corresponding decrease in Vs. The transition from a high velocity to a low velocity zone will result in a trough on zero-phase-processed seismic data. The peak on VET2 (Figure 3.10) at 110 ms seems to be related to the top of the coarse component of the uppermost fining upwards sequence. A peak at 190 ms corresponds with an increase in velocity that appears to occur within a finer grained sand unit with a relatively low tip resistance. The base of the lower sequence is also represented by a trough on VET2 at a time of 260 ms.  46  An obvious boundary is evident on all of the geotechnical data at approximately 23 m that is a response due to the change from a sand dominated unit to a silty unit. This silty unit has a significant amount of clay in the pore spaces. The clay component is apparent by examining the tip pore pressure plots (U2), where a large change in pore pressure is evident at a depth of 23 m. This finer grained unit is continuous down to the Pleistocene till (46.5 m for VSP1 and 49.0 m for VSP2). The tip resistance is very low in the clayey-silt unit and shows little variation. The friction-ratio, pore-pressure and shear-wave- velocity plots show subtle variations, suggestive of variations within the silt, that are confirmed to some degree by the presence of relatively small amplitude reflections on VET1 and VET2. The depth of cone refusal was at the top of the Pleistocene till. The VSP data were able to resolve the top of the Pleistocene till and reflections from interfaces within the till. The top of the Pleistocene till represents a change from a low velocity clay to a high velocity till, observed as a peak on VET1 at 495 ms and on VET2 at 505 ms. In VET2 there are a series of peaks and troughs above the Pleistocene till, starting at 460 ms that seem to correspond to a low-velocity zone, seen on the shear-wave velocity and friction-ratio plots, which may be due to a sand or silt layer. The true nature of this layer cannot be confirmed because there is no core at this depth. The reflections within the Pleistocene till can be seen in both VET1 and VET2 down to times of approximately 0.75 s or depths of approximately 80 m. Thus, data were obtained from beyond the penetration depth of the cone; this could be extremely valuable at other sites where there is concern about contamination or hazards beyond a certain depth.  3.6  A Comparison of V S P and C D P Data One of the main reasons for acquiring V S P data was to improve the analysis of CDP  reflection data. In order to accomplish this objective there should be a reasonable correlation 47  between the VETs and the CDP data. Figure 3.11 illustrates the surface SH-wave reflection seismic data with the two VETs inserted at the tie points. The overall reflectivity trends compare favorably between the surface reflection data and the VSP data. The reflectivity is moderate within the sand unit, very low within the clayey-silt and high within the Pleistocene till. Individual reflections also match extremely well. The base of the lower sand sequence (at a time of approximately 270 ms) is obvious on both the CDP data and the VETs; the relative amplitudes are also comparable.  The reflection from the top of the Pleistocene till at a time of  approximately 500 ms also matches from the VETs to the CDP data. The top of the Pleistocene till is evident as well as a number of reflections within the till. A series of smaller amplitude reflections from 300 to 450 ms also seem to match well and are suggestive of velocity and/or density variations within the lower clayey-silt layer. The peak at a time of 200 ms on V E T 2 seems to be the most significant mismatch between the VETs and the CDP data. There is no reflection on the CDP data that compares in relative amplitude.  There is however, a lower amplitude reflection on the CDP data so it probably  represents differences in the amplitude preservation of the two techniques. The CDP data were dominated by significant source generated noise in the upper 300 ms, so the relative amplitudes are probably not well preserved. The V S P technique is not affected by this source generated noise so the amplitude information should be more reliable. The higher frequency content of the VSP data is evident. The VETs show a series of distinct reflections just below the top of the Pleistocene till, while the same reflections on the CDP data are overlapping and complex. As has been demonstrated by other researchers (Milligan et al., 1997; Hardage 1985), this difference in frequency content allows for smaller features to be resolved in the VSP data than in the CDP data.  48  3.7  Conclusions The seismic cone penetrometer was an effective method for obtaining high quality VSP data  at the Kidd2 site. The cone VSPs have a number of advantages compared to conventional borehole V S P surveys including decreased cost, better acoustic coupling and the avoidance of the problems caused by tube waves. The V S P data were acquired along with other cone data, thereby providing a matched set of data for subsurface interpretation. The seismic reflectivity is well understood and limited to variations in velocity and density, while the cone log response is governed by additional parameters such as grain size and effective stress and relies on empirical relationships for interpretation. Subsurface velocities measured in situ are often not available for near-surface studies and are important for the processing of surface CDP seismic data as well as the accurate conversion of time to depth. The oil and gas industry has long recognized that the analysis of VSP data can significantly improve the interpretation of reflection seismic data (Hardage, 1985).  Ambiguities in the  interpretation of reflection data (such as multiples) are reduced and a precise reflectivity model can be generated to serve as the starting point for detailed seismic inversion. The inversion of seismic data based on seismic cone VSP constraints has the potential to become an effective method for characterizing the near-surface.  49  Figure 3.1: Location map showing the Kidd2 research site on the Fraser River delta with the grid of SH-wave CDP seismic (dotted lines) and the two VSPs, labeled VSP1 andVSP2.  50  Time (s)  (Upgoing Reflection)  Multiple)  Figure 3.2: Data from VSP1. A 30-150 Hz bandpass filter and spherical divergence correction have been applied to the data.  Time (s)  Pleistocene Till (Upgoing Reflection)  Figure 3.3: Data from VSP2. A 30-150 Hz bandpass filter and spherical divergence correction have been applied to the data.  51  Figure 3.4: Interval velocities derived from the VSP first break arrival times.  Figure 3.5: Processing flow chart for the cone penetrometer VSP data. 52  Time (s)  Figure 3.6: Processed data from VSP1 (bottom). The upgoing wavefield has been aligned and muted. The stacked data trace is duplicated and is shown above (VET1).  53  Time (s) 0.0 iii  0.1  0.2  0.3  0.4  0.5  0.6  0.7  i i i i i i i i i i i i i i i i i i i i i ' i i i i ' i i i i ' i i i i i i i i i i i i i i i i ' i — i i i 11 i i i i i i i i i  Figure 3.7: Processed data from VSP2 (bottom). The upgoing wavefield has been aligned and muted. The stacked data trace is duplicated and is shown above (VET2).  54  Offset (m) 12  Offset (m) no ,  1  1 2  0.1 0.2 0.3  IH 0-4 0.5 0.6 0.7 Shot Profile 1252  Shot Profile 1004  Figure 3.8: Two representative shot profiles from Line 20 after application of gain, trace balancing and record subtraction. Shot profile 1004 is at the south end of the line and shot profile 1252 is at the north end.  FR(%)  Qc(MPa)  U2 (MPa)  Vs (m/s)  Figure 3.9: Composite plot showing the relationship of the cone penetrometer friction ratio (FR), tip resistance (Qc) and pore pressure (U2) to the first break shear wave velocity (Vs) and VSPl's extracted traces (VET1). The data are plotted in time with the equivalent depths indicated on the left. 55  FR(%)  Qc(MPa)  U2 (MPa)  Vs (m/s)  VET2  Figure 3.10: Composite plot showing the relationship of the cone penetrometer friction ratio (FR), tip resistance (Qc) and pore pressure (U2) to the first break shear wave velocity (Vs)and VSP2's extracted traces (VET2). The data are plotted in time with the equivalent depths indicated on the left.  South  VET2  \  1 0 M  4  H  VET1  North  Figure 3.11: SH-wave CDP reflection data with the extracted traces from VSP1 and VSP2 (VET 1 and VET2) inserted. 56  CHAPTER 4 Aquifer Heterogeneity from SH-Wave Seismic Impedance Inversion 4.1  Introduction A critical part of many groundwater or environmental studies is obtaining quantitative  information about the hydrogeologic properties of groundwater aquifers. Numerous techniques have been developed to provide such information; these techniques vary in terms of the scale and the accuracy of the measurement.  Most of the commonly used methods of subsurface  hydrogeologic testing (e.g. tracer tests, pumping tests, slug tests) require a well bore.  The  drilling of the well bore can be expensive, so there are often too few wells to provide good spatial coverage. In addition, there can be problems associated with the invasive nature of this form of testing; this is of particular concern in characterizing contaminated regions. A n alternate approach is to develop geophysical methods as a means of non-invasive imaging or sampling of the subsurface. Hydrogeologists have recognized the benefits that can be gained by having seismic data from a subsurface region of interest.  A synthetic model study by Copty et al. (1993)  demonstrated that sparsely sampled permeability and pressure data can be combined with densely sampled seismic velocity data to invert for densely sampled permeability; a relationship between P-wave velocity and permeability developed by Marion et al. (1992) was used to make the connection between seismic and hydraulic properties. In general there does not exist a unique relationship between velocity and permeability, so two other studies have avoided this difficulty by developing ways of calibrating and incorporating the seismic information. The study by McKenna and Poeter (1995) used crosshole seismic tomographic data to identify certain ranges in P-wave velocity values with distinct hydrofacies; permeabilities were assigned to the units using measurements on core and from packer tests. Hyndman and Gorelick (1996) 57  also defined distinct zones in the subsurface based on measured seismic velocities, and assigned uniform hydraulic properties to these zones using information from the inversion of tracer test data. The basic concept behind all of these studies is the idea that the subsurface can be divided into regions each of which has a distinct range of seismic velocity. These regions, defined by their seismic properties, can then be assigned hydraulic properties. Regardless of the specific methodology followed in using seismic data for hydrogeologic applications, all studies to date illustrate the importance of two key steps: the first is obtaining an accurate model of the velocity variation in the subsurface;  the second is providing the link between velocity and hydraulic  properties. Obtaining velocity information from seismic data can take two approaches. One approach is to use existing boreholes to carry out crosshole tomographic experiments. requires the placement of cased boreholes in the area to be studied.  This procedure  In addition to being  relatively expensive, cased boreholes pose risks in contaminated aquifers. A second approach uses surface seismic data.  Surface seismic is non-invasive and provides a detailed two-  dimensional image of the subsurface. The focus of this study is the use of surface based shear wave (S-wave) seismic reflection data as a means of describing the heterogeneity of a sand aquifer; specifically I use the seismic data to obtain a void ratio model of the aquifer. While there have been a number of studies illustrating the usefulness of compressional wave (P-wave) reflection data for near-surface applications (e.g. Steeples and Miller, 1990; Biiker et al., 1998), there has been surprisingly little work done on the use of S-wave data, with the exception of the few studies by Hasbrouck (1991), Harris et al. (1996), and Carr et al. (1998). A n obvious advantage of S-wave data in many shallow environments is the improved vertical resolution, which is taken to be 1/3 of a wavelength (Miller et al., 1995); in water-saturated sediments the S-wave velocity is generally much less than the P-wave velocity, resulting in a smaller wavelength. 58  Another significant  difference between P-wave and S-wave techniques is the lack of sensitivity of the shear modulus to fluid saturation. This characteristic of S-wave techniques results in a larger change in S-wave velocity compared to P-wave velocity for a given variation in void ratio or lithology. The first issue, which I addressed is the critical step of obtaining an accurate velocity model from the collected data. Interval velocities can be obtained using the Dix equation (Dix, 1955) applied to the stacking velocities, but this technique assumes horizontal reflectors and is prone to large unconstrained errors. I chose to use an impedance inversion algorithm (e.g. Lines and Treitel, 1984; Tarantola, 1984), which uses the seismic trace amplitudes. In order to obtain usable results from impedance inversion, velocity constraints are required that are typically obtained by velocity logging in a nearby cased borehole. In this study I obtained the.nearsurface velocity information, required for the inversion, from the seismic cone penetrometer. This is an attractive alternative to the use of the more expensive and invasive borehole logging, and is ideally suited for the inversion of shallow surface seismic data. The second important step, i f seismic methods are to be used for hydrogeologic applications, is to convert the velocity model to a "map" of lithology. M y approach involved using information derived from the cone seismic data and from laboratory measurements; the final product being a lithology and void ratio model of the aquifer system. The overall objective in this study was to investigate the optimal way of using shear wave seismic reflection data to estimate hydraulic properties. A novel aspect of the methodology was the use of cone seismic data as an integral part of both the inversion of the surface seismic data and the interpretation of the velocity model.  I found that cone-measured data can contribute  significantly to the use of near-surface seismic data for hydrogeologic applications.  59  4.2  Geologic Setting The location of the investigation is shown in Figure 4.1. The work was carried out on the  Kidd2 research site located adjacent to a B C Hydro electric substation on No. 4 Road, Richmond, British Columbia on the Fraser River delta. The site has been the focus of research by both hydrogeologists (e.g. Wood, 1996 and Neilson-Welch, 1999) and geotechnical engineers (e.g. Hoffmann, 1997). The geotechnical engineers were interested in measuring in situ soil properties using the cone penetrometer and the hydrogeologists were interested in measuring the hydraulic conductivity and modelling groundwater flow within a subsurface aquifer.  The  Geological Survey of Canada has also used the site for extracting core and testing geophysical surveying techniques (e.g. Hunter et al., 1998b) as part of an earthquake site amplification study. The Fraser River delta is composed of sediments deposited during the progradation of the delta during the Holocene. Clague et al. (1983) identifies the typical characteristics of Fraser River deltaic sediments.  The surficial sediments were deposited in a floodplain and upper  intertidal environment and consist of fine-grained silts, with possible peat deposits.  The  underlying sediments were deposited in estuarine river channels in an intertidal environment and consist of sands and silts.  Sands deposited in a delta-slope environment underlie the channel  deposits and tend to be more uniform in texture. Beneath the sand/silt sediments lie clay-rich silt sediments deposited in a marine basin environment.  A compact Pleistocene glacial till lies at  the base of the Holocene sediments. Based on work done by Wood (1996) and Neilson-Welch (1999) I have developed a simple model of the site: silt and clay in the upper 4 to 5 m; sands form the aquifer to a depth of 21 to 23 m; marine clays and silts to a depth of 49 to 52 m; and the Pleistocene till underlying the entire sequence. Core analysis by Wood (1996) revealed a series of fining upwards sequences within the sand aquifer that divides the aquifer into upper and lower zones.  The upper aquifer  represents the uppermost fining upwards sequence and has a depth range from 5 to 13 m. The 60  lower aquifer has a depth range of 13 to 22 m. The upper aquifer has a number of very fine sand and silt layers up to 10 cm in thickness that are probably related to river channel deposition along with a few thin zones (1 to 2 cm) of clay. The lower aquifer tends to be uniformly composed of sand. Hoffmann (1997) measured the void ratio in the frozen cores retrieved predominantly from the lower aquifer (from 12 to 17 m) and obtained values ranging from 0.88 to 1.12.  The  water table is approximately 1 m below the surface at the site with some variation in depth due to tidal influences  4.3  SH-Wave Seismic Reflection Data  4.3.1  Data Acquisition  The common depth point (CDP) shear wave seismic data were collected in June of 1996 using a 24 channel Geometries SmartSeis seismograph, 48 Geo Space horizontal geophones (20DM-28 Hz) attached to CDP spread cables and an I/O RLS-100 roll switch. The data sample interval was 0.0005 s and a 600 Hz anti-alias filter was applied before analog to digital conversion. A n end-on shooting configuration was used with a geophone and source spacing of 1 m, giving a far offset of 24 m and 1200% subsurface coverage. The seismic source consisted of a 35 kg steel block with metal fins attached to the base to increase the source coupling. SHwaves were generated by swinging a 5 lb sledgehammer against the sides of the block. To further increase source coupling the person swinging the hammer stood on the block. The SH-wave technique was utilized to minimize mode conversions and P-wave interference. An average of 10 hammer blows were stacked for each record to improve the signal to noise ratio and remove the effects of the occasional misdirected blow. Data were recorded with blows in opposite directions to allow for the subtraction of the records to further reduce P-wave interference and provide additional signal to noise enhancement of the SH-waves. A total of five  61  lines of SH-wave seismic data were acquired and processed, covering 420 m. The locations of these lines are shown in Figure 4.1. This paper will only examine line 20, a 100 m north-south line (highlighted in Figure 4.1).  4.3.2  Data Analysis and Processing  The seismic data were processed using ITA's UNIX based seismic processing system installed on a Sun workstation network.  The processing flow was the same as that commonly  used in the petroleum industry and consisted of: spherical divergence correction, bandpass filtering, record subtraction, deconvolution, normal moveout removal, residual statics, stacking and migration. Some of the key parameters used in the processing are summarized below. A spherical divergence correction was applied using a time scaling function with a power of 2.0.  The data were bandpass filtered (25/50-150/400 Hz) using a minimum phase Butterworth  filter.  The record subtraction process was preceded by the application of a trace to trace  amplitude balance. The surface consistent deconvolution used both shot and receiver gathers, with trace by trace editing required afterwards to remove the noisy trace segments.  The  elevation statics compensated for a topographic variation of less than 1 m across the site. A semblance velocity analysis was used iteratively with surface consistent residual statics. The mute function tended to be quite conservative to ensure that the amplitude variations at shallow depths did not dominate the stack. The data were migrated using an interval velocity model based on the extensive data from the site that defines the depth of key reflecting horizons. Two shot profiles from both ends of the line are shown in Figure 4.2. These profiles were created after the shot records had been subtracted. Hyperbolic reflections are apparent at a time of approximately 0.15 s. A reflection from the top of the Pleistocene till is present at a time of approximately 0.5 s.  62  The migrated dataset is shown in Figure 4.3.  The reflectivity variations on the seismic  section are directly related to the dominant lithologic units. The upper 0.3 s has relatively high reflectivity due to lithologic variations within the sand aquifer. From 0.3 to 0.5 s the reflectivity is reduced due to a lack of lithologic variation in the underlying marine clay/silt sequence. The reflection off the Pleistocene till is the most dominant coherent reflection at a time of 0.5 s that results from a very large velocity change (Hunter et al., 1998a).  A large channel is obvious  within the sand unit, extending from CDP 1112 to the northern end of the line at a time of approximately 0.15 s.  This channel is most likely a distributary channel formed during the  progradation of the delta.  Other reflections within the aquifer are believed to be caused by  depositional interfaces; i.e. interfaces within the laterally extensive sand unit were noted by Clague et al. (1983) and attributed to delta slope texture variations.  4.4  Seismic Cone Penetrometer Data The cone penetrometer is a technology developed by geotechnical engineers for rapid in situ  site characterization and is considered minimally invasive. A steel rod with a cone shaped tip is pushed into the ground hydraulically while making measurements with sensors mounted close to the tip (Robertson et al., 1986).  The technique works extremely well in deltaic deposits  consisting of sand, silt or clay. In the past, most of the data obtained with the cone penetrometer consisted of standard cone parameters such as tip resistance, sleeve friction, and pore pressure response.  Over the years additional sensors such as geophones or accelerometers have been  added to the cone to allow measurement of seismic arrival times.  These seismic cone  penetrometers are referred to as SCPTU (Robertson et al., 1986). Conventionally, shear wave velocities have been determined from first arrival or first cross-over times of polarized shear waves generated at the ground surface during pauses in penetration in a SCPTU sounding. Jarvis and Knight (2000) [Chapter 3] outlined an extension of the standard SCPTU sounding to 63  carry out vertical seismic profiling (VSP). The VSP survey is a way to obtain high resolution reflection images of the subsurface using a subsurface sensor (typically in a borehole) and surface seismic source. The cone penetrometer is particularly well suited to VSP surveying due to the increased coupling of the sensor, the decreased cost when compared to conventional borehole VSPs, and the reduction of unwanted wave arrivals such as tube waves. In this study the cone V S P data were used in two ways to aid with the inversion of the surface seismic data. The determined velocities were essential for obtaining the initial velocity model of the subsurface that provides the constraints for the inversion of the surface seismic data. In addition the V S P data provides a way to identify the seismic waveform (another key constraint for the inversion). A third use of the VSP data was in the interpretation stage, where the observed variation in velocity with depth was used to remove the effect of effective stress on the velocity model.  4.4.1  Data Acquisition  The cone penetrometer has been used extensively in a number of other studies at the Kidd2 site. Most of the cone holes have involved the acquisition of data associated with standard geotechnical parameters such as tip resistance, sleeve friction and pore pressure. Velocity data have been acquired in very few holes. In this investigation the cone accelerometers were used to obtain near-surface vertical seismic profile (VSP) data, as described in detail in a recent paper by Jarvis and Knight (2000) [Chapter 3]. Two shear wave VSP surveys, referred to as VSP1 and VSP2, were conducted along line 20, located as shown in Figure 4.1, using the cone penetrometer truck from the Department of Civil Engineering at the University of British Columbia. A large sledgehammer against the baseplate of the truck was used as a source of SH-waves and a cone-mounted accelerometer used as a receiver. Bi-directional hammer blows and record subtraction were also used in this experiment 64  to optimize the signal. The data were recorded using a Nicolet oscilloscope with 15 bit data digitizing and storage capability. The depth levels at which the cone was stopped to record VSP data were chosen to minimize field time and spatial aliasing. The truck engine (which must be on to run the hydraulic system) was turned off before every measurement to reduce the noise.  4.4.2  Interval Velocity Logs  The interval velocities are obtained from the first arrivals of the V S P surveys. The interval velocity logs from the V S P surveys are shown in Figure 4.4a. Each V S P trace is recorded at a particular depth and the arrival times are used to calculate the interval velocity between adjacent depth recordings. The accuracy of the interval velocities is dependent on the accuracy of the picking of the arrival time as well as the accuracy of the assigned depths. The original data were recorded at a sample interval of 0.5 ms, which is adequate to record the frequencies of the data (without aliasing) but insufficient for arrival time picking. To increase accuracy the data were sub-sampled down to 0.05 ms. The depths are based on the lengths of rod used by the cone penetrometer, which are machined to very close tolerances and should be accurate to within 1 mm. It should be noted that the first arrivals (and interval velocities) could have been obtained using standard seismic cone measurement techniques, without the requirement of recording complete VSP traces.  4.4.3  VSP - Data Analysis and Processing  In addition to the use of the first arrivals in the VSP data to obtain information about subsurface velocities, the V S P data also provided a way of defining the seismic wavelet in the surface CDP seismic data. The seismic wavelet is one of the required inputs to the impedance inversion algorithm.  The processing of the VSP datasets followed a generally accepted  processing flow (as outlined by Hardage, 1985). The main objective of the processing was the 65  separation of the upgoing and downgoing wavefields, and the improvement of the signal bandwidth using deconvolution. The processing is detailed in Jarvis and Knight (2000) [Chapter 3] and consisted of singular value decomposition, minimum phase deconvolution, match filtering and f-k filtering. The final step involved the stacking of the data to create the VSP extracted traces (VET). For this particular investigation the VETs played a critical role in allowing us to characterize the wavelet of the S-wave surface seismic data.  At a subsurface location sampled by both the  VSP and the CDP data, I determined the phase shift that had to be applied to the migrated CDP data so as to match the known zero phase of the VETs. This phase shift was then applied to all the CDP data resulting in an embedded wavelet that can be assumed to have zero phase.  4.5  Impedance Inversion of Shear-Wave Data The goal of this research was to provide a way to map the heterogeneity within the Kidd2  aquifer. The shear-wave seismic reflection data provide a well-sampled image of the aquifer. The challenge is to convert the amplitude variations of the seismic traces to variations in subsurface properties that demonstrate the heterogeneity. The conversion of seismic amplitudes was accomplished using an impedance inversion algorithm. The algorithm was developed at the University of British Columbia by M . D . Sacchi and T.J. Ulrych. In this section I will outline the impedance inversion algorithm (as detailed by Sacchi and Ulrych, 1996) and show the connection between impedance and shear-wave velocity.  Finally, I will demonstrate the  application to the shear-wave seismic reflection dataset. The seismic impedance inversion is in essence a seismic reflectivity inversion from which the seismic impedance can then be calculated. A Bayesian approach is followed that allows a priori information and associated error to be incorporated into the inversion. The convolutional model of a seismic trace is used, which implies that the seismic trace is equal to the convolution 66  of a reflectivity series and a wavelet, plus random noise (assumed to be Gaussian). Impedance is obtained by integrating the reflectivity series. Constraints on the impedance will have errors associated with them that are also assumed to be Gaussian. These errors are due to inaccuracies or uncertainties in the measurements of impedance. Sacchi and Ulrych (1996) provide the final cost function (J) in matrix form associated with the inversion for reflectivity (x): J = aJ +±1| - ( W x - y) || 2"<r •"" x  2  +11| S - (Cx 2  )  (4.1)  where W is the wavelet matrix, y is the seismic trace vector, o~ is the standard error in the seismic trace amplitudes, C is an integrator matrix, £, is the impedance constraint vector and S is a diagonal matrix of errors associated with the impedance constraints. The noise in the seismic trace and the error in the impedance constraints are the only source of misfit between the model reflectivity series and the data. The Sacchi and Ulrych (1996) algorithm minimizes the cost function subject to standard % statistical measures. The regularization term J is based on the L i 2  x  criterion that enables sparseness to be specified for the model reflectivity series as:  (4.2)  Sparseness in the reflectivity series is a key component of the inversion. A sparse reflectivity series will consist of a limited number of large reflection coefficients that result in a blocky impedance model and are probably the most representative of impedance variations in the earth. The degree of sparseness is controlled by the weighting parameter a.  This weighting parameter  must be solved for in the inversion and involves a computationally intensive search for a value that matches the data misfit to the  statistic (% - number of traces + number of impedance 2  constraints). Due to the computational effort required to obtain a, the value of a was determined  67  for a representative part of the seismic line and then held constant for the inversion of all other traces. The inversion of the seismic data proceeds using a non-linear conjugate gradient algorithm (Sacchi and Ulrych, 1996).  The output reflectivity model relies on maximizing the a posteriori  probability density function ( M A P solution).  This Bayesian approach results in a reflectivity  model that honors the probabilities of the priors and adheres to the sparseness criterion. The final impedance model was obtained by integrating the obtained reflectivity model. Seismic impedance is the product of velocity and density.  To convert from seismic  impedance to shear-wave velocity I have made a constant density assumption. For this particular site, the water table is within one meter of the surface and all sediments are saturated so this is taken to be a reasonable assumption. With a constant density assumption the impedance model is equivalent to a SH-wave velocity model.  The SH-wave velocities determined using the cone  penetrometer and interpolated between cone holes formed the input velocity model and were used as impedance constraints in the inversion. The a posteriori reflectivity model output from the inversion was converted to impedance and resulted in the output velocity model. The extensive data available from the Kidd2 site are the basis of the input velocity model. The data consisted of depth picks and interval velocities from cone penetrometer logs, an interpretation of the reflection times on the surface seismic data, and previous investigations (e.g. Neilson-Welch, 1999; Hofmann, 1997), which provided information on the basic geologic framework and the expected depths and continuity of the different sedimentary units. The input velocity model is shown in Figure 4.5. The model is based on three control points. The first control point is the location of VSP2 at CDP 1080 and the second control point is the location of VSP1 at CDP 1216. At these control points interval velocity data from the seismic cone penetrometer are available from a depth of 3 m down to the top of the Pleistocene till. The velocity of the Pleistocene till had to be estimated because the cone could not penetrate into the 68  Pleistocene till. Based on work by Hunter et al. (1998a), the velocity of the Pleistocene till is at least double the velocity of the overlying deltaic sediments; the model reflects this large increase in velocity. Due to the variations to the north of CDP 1216 a third control point was inserted at the north end of the line. At CDP 1373 the data from a cone hole were available. Unfortunately, there were no velocity data, only tip resistance and sleeve friction logs. In the upper 10 m, where we interpreted there to be the low velocity edge of a channel, velocities were obtained from the measured tip resistance, using a relationship between velocity and tip resistance established with the data from VSP1 and VSP2. Below 10 m the velocities from VSP1 were used. The reflection time (on the CDP data) at the base of the aquifer and the top of the Pleistocene till were used as a guide to position the VSP1 velocities. A comparison of the input and output velocity models is shown in Figure 4.6. The two sets of velocity traces represent the location of the two VSPs, where the input velocity model has been obtained from the interval velocity logs. The dashed line is the input velocity model and the solid line is the output from the inversion. The increase in velocity at the Pleistocene till is readily apparent in the 47 to 50 m depth range. The trends in the velocities are very similar thereby demonstrating the viability of the inversion. The variations in the output model below the top of the Pleistocene till result from the presence of a series of closely spaced reflectors in the seismic dataset, which are likely to be associated with velocity variations within the Pleistocene till.  However, no velocity information is available to properly constrain the  inversion of this zone so the accuracy of the output model in this region is questionable. One test of the success of any inversion is the ability to reproduce the input data.  The  derived velocities can be converted into reflection coefficients and convolved with the source wavelet to create a seismic dataset.  This seismic dataset is shown in Figure 4.7.  69  When  compared with Figure 4.3, there are very few differences  between the two datasets,  demonstrating that the inversion output is a viable solution. The output velocity model from the inversion is shown in Figure 4.8. The velocity variation is similar in form but differs in the detail when compared to the input model. The continuity of the velocities in the input model is due to the interpolation of three control points. The output model does not show this same continuity because the recorded changes in seismic amplitudes have been honored by the inversion. Our working hypothesis is that the variation in seismic velocity of the subsurface can be used to map out heterogeneity in the hydraulic properties of the aquifer.  In Figure 4.8 we have  acquired, through inversion of our collected surface seismic data, a model of seismic velocity. The next step was to transform this into an image showing sand and clay regions. The separation of the sand unit and the use of laboratory measurements made it possible to determine spatial variability in the void ratio of the sand.  4.6  Differentiating Lithology Using the Velocity Model  A desirable objective in using seismic data for hydrogeologic studies is to be able to relate the measured velocity to permeability. however, is not straightforward.  The relationship between permeability and velocity,  The use of laboratory-derived relationships, such as that  presented by Marion et al. (1992), requires key assumptions about the microstructure of the in situ material; e.g. a sand framework with constant intergranular void ratio. Alternatively some means of calibrating field-measured velocity data with field-measured permeability is required; as in the study by McKenna and Poeter (1995). In my study, given insufficient knowledge about the subsurface materials to accurately obtain permeability information, I elected to use the velocity model to determine variation in void ratio.  70  As with most rock physics relationships, the relationship between velocity and void ratio can be better defined i f applied to single lithologic units. In this study I was primarily interested in the sand aquifer. In order to identify the distinct lithologic units in my velocity model I had to first remove the influence of effective stress, which causes a velocity increase with depth and dominated the output model shown in Figure 4.8. Correcting for effective stress enabled me to identify the individual lithologic units, and made it possible to see more clearly the velocity variation within the sand that was related to variation in material properties.  4.6.1  Correcting For Effective Stress in Velocity Data  A complicating factor in the interpretation of any near-surface velocity measurements is the effective stress, which in general acts to increase the shear wave velocity with depth. Removing this effect from the measured velocity is a procedure referred to as stress normalization. The concept of stress normalization was developed by geotechnical engineers to aid in identifying lithologies based on seismic velocity. Robertson et al. (1992) used a formula based on lab studies by Hardin and Drnevich (1972). The normalized velocity V  V =V (P /rj' ) - , 0  sn  s  a  s n  is given by  (4.3)  25  vo  where V is the measured shear wave velocity, P is the atmospheric pressure (~100kPa) and a ' s  a  v0  is the vertical effective stress (which is obtained by multiplying the unit weight of soil by the depth and subtracting the hydrostatic pressure).  It has been suggested that the value of the  "stress exponent" is not a constant of 0.25 but is variable and possibly site-specific. Other researchers (e.g. Tokimatsu et al. 1991; Lee and Campbell, 1985) have found exponents ranging from 0.31-0.46 in their normalization function.  71  Due to the variability of the stress normalization exponent, I used the local velocity data from the interval velocity logs to derive site-specific normalization parameters. Referring to Figure 4.4a, the increase in velocity due to effective stress is most obvious in the upper 10 m. The normalization function, which I chose to use, was determined by examining the velocities within the thick, uniform clay layer interpreted based on CPT data to extend from 22 to 50 m. The clay layer was apparent on the cone penetrometer data with low tip resistance and high friction ratio (Woods, 1996).  The correct normalization function should produce a constant normalized  velocity in a layer of uniform lithology. I found that the normalization function used by Robertson et al. (1992) was unsatisfactory for this data set.  The problem with this particular relationship is the fact that it produces an  unrealistically large normalized shear wave velocity close to the surface, where o"' is very small. v  At the surface, where a '  is zero, V  v0  s n  is predicted to be ~2V . Therefore, a modified form of this S  expression was used to normalize the velocity data:  V n=V [P /(a' +k)] - , 0  S  s  a  (4.4)  3  vo  where k is a constant that provides an offset of the velocity normalization function at the surface. For this study k was taken to be 10 kPa. Lee and Campbell (1985) also recognized the need for a constant (k) to account for the non-linear intersection of shear-wave velocity at the ground surface. The stress exponent was increased to 0.3 partly to compensate for the shifting of the exponential function and ultimately to provide the best normalization function for this site.  The  quality of the normalization function was evaluated by ensuring that the clay zone velocities fell within a relatively constant range and the sand aquifer velocities exhibited neither overcorrection nor under-correction of the stress-effect.  72  The normalized velocity data are shown in Figure 4.4b.  The differentiation of velocity  between the clay dominated and sand dominated units is now apparent.  The clay-dominated  unit (below 22 m) has a constant normalized velocity of approximately 160 m/s.  The sand  aquifer (5 to 22 m) has a higher normalized velocity ranging from 160 to 210 m/s. For purposes of separating the section into "sand" and "clay" we selected 170 m/s as the cutoff; this is the maximum normalized velocity found within the clay-dominated unit. The normalized output model is shown in Figure 4.9. The upper and lower bounds of the aquifer are readily apparent on the color plot, with a scale chosen to make the higher velocity silts and sands appear red and yellow and the low velocity clay-rich material to appear in blue and green. The scale has been chosen to focus on the aquifer velocities; the large increase in velocity at the Pleistocene till is off-scale and shown in white. The base of the sand channel in the upper 10 m is marked by a low velocity zone which probably represents a fine-grained silt or clay zone. The impedance inversion results demonstrate that impedance inversion can be used to map the velocity variations within a near-surface region. The output is highly dependent on the quality of the input and the absolute velocities will still have error associated with them. The relative variations should be valid and are related to the amplitude variations on the input seismic data.  Having removed the stress effect by normalizing the output velocities, normalized  velocities have been used to differentiate lithologies and can now be interpreted to obtain information about void ratio.  4.6.2  Void Ratio in Sand Aquifer  The primary objective of this study was to develop a means of using surface seismic data to delineate the heterogeneity in the properties of a sand-dominated aquifer. The profile of shear wave velocity output from the impedance inversion provides some insight into the variations in 73  the aquifer.  Rather than consider these variations in velocity, they can be transformed to  illustrate variations in void ratio; a parameter which has more relevance for hydrogeologic investigations. Previous work (Wood, 1996) has described the aquifer as an upper aquifer composed primarily of sand and silt with a few very thin clay stringers and a lower aquifer composed of a uniform sand; the absence of any significant amount of clay simplifies the interpretation of the velocity data in terms of void ratio. The relationship between S-wave velocity and void ratio in such a system was studied by Hardin and Richart (1963) through an extensive set of laboratory experiments on dry materials ranging from coarse sands to fine silts with void ratios from 0.7 to 1.3. Their laboratory study showed that velocity depends primarily on void ratio and confining pressure, with very little dependence on grain size, and resulted in the following empirical relationship: V =(104-35e)a s  (4.5)  1/4 o  where e is void ratio and a is confining pressure (in kPa). 0  In order to adapt this relationship for use in my study, I needed to account for the saturated state of the aquifer unit. This was done by adding a term to Equation 4.5 to account for the saturated density, and defining a revised empirical relationship from the best-fit linear relationship between V and e values, generated using the modified form of Equation 4.5. In s  making this density correction, I assumed a density for the solid phase of 2.65 g/cm . The 3  resulting empirical relationship is V =(96-37e) a s  1/4  (4.6)  0  which can also be written as:  e = 2.6-  (4.7) 1/4  -  74  Working with the normalized velocity data (Figure 4.9), the confining pressure is removed as a variable and is fixed at the normalizing pressure of 100 kPa. The velocity values, which ranged from 170 m/s to 205 m/s within the sand aquifer were converted to void ratio values, using Equation 4.7. The calculated void ratios ranged from 0.85 to 1.15, in excellent agreement with the range of 0.88 to 1.12 determined from the C A N L E X experiment (Hoffman, 1997). The final plot of lithology and void ratio is shown in Figure 4.10. The upper and lower aquifers can be identified in the figure by the narrow transition zone from low to high void ratio at a depth of approximately 13 m.  The upper aquifer shows both vertical and horizontal  variations in void ratio that most likely represent the wider distribution of sediments expected in a higher energy near-shore environment.  The lower aquifer shows a relatively uniform  distribution of void ratio, with void ratio decreasing towards the bottom of the aquifer.  4.7  Conclusions  Using seismic velocities to quantify hydraulic properties within a near-surface aquifer has a number of challenges, which must be addressed. The first is to obtain reliable velocity values. For near-surface studies standard velocity measurement techniques (borehole based sonic logging tools) are not always available or affordable so non-traditional techniques such as the cone penetrometer and surface based reflection seismic must be used wherever possible. The reflection seismic technique is advantageous because the result is a very dense subsurface coverage and the cone penetrometer has the advantage of rapid deployment, reduced cost, and the simultaneous acquisition of other cone logs for stratigraphic interpretation. Employing a Bayesian inversion technique to convert the seismic amplitude variations to velocity variations is a robust method that honors the probabilities of the priors and adheres to a geologically reasonable sparseness criterion. 75  The second challenge involves the separation of lithology using the seismic velocity data. The near-surface velocities are dominated by the rapidly increasing effective stress and this effect must be removed from the velocities. The normalization of the velocities accomplishes this task.  The analytic function used to normalize the data will be dependent on the local  sediments and their depositional history so local velocity data are critical to deriving a sitespecific normalization function.  Once the velocities have been normalized a separation of  lithologies may be possible depending on the types of sediments under investigation.  The  availability of core and independent stratigraphic information from the site enables the sediments to be associated with ranges of velocity. As long as the ranges do not overlap the lithologies can be separated. The final conversion to hydraulic properties can be achieved on a lithology specific basis. Using empirical or derived relationships each lithologic unit should be treated separately.  The  end result is a spatial map of hydraulic properties that provides direct insight on the subsurface heterogeneity.  76  236° 30'  237' 00'  236° 30'  237* 00'  237' 30'  237° 30'  Figure 4.1: Location map showing the Kidd2 research site on the Fraser River delta with the grid of SH-wave CDP seismic (dotted lines) and the two VSPs, labeled VSP1 and VSP2 that were used for impedance constraints. The seismic line (line 20) used in the inversion has been highlighted. Offset (m)  Offset (m)  Shot Profile 1004  Shot Profile 1252  Figure 4.2: Two representative shot profiles from line 20 after application of gain, trace balancing and record subtraction. Shot profile 1004 is at the south end of the line and shot profile 1252 is at the north end. 77  78  fr  r S  •  <  >— Q if) CM  O O CM  O LO ••-  • O O •<-  O ID  O  (s/iu) U S A  o  o  o  O CO  o  LO CM  o  O CN  o  LO  T—  o  O  T—  (s/iu) S A  o  LO  o  o  O O CO  Q LO CN  O O CM  CD LO •<-  O O i-  (s/iu) S A  79  O LO  O  to  Figure 4.5: Input velocities to the seismic impedance inversion.  80  Figure 4.6: A comparison of the velocity logs from the two VSPs [a) VSP2 b) VSP1 ] to the output from the seismic impedance inversion. The line with diamonds is the velocity log and the solid line is the inversion output.  81  82  CDP Figure 4.8: Velocities output from the seismic impedance inversion.  1050  1100  1150  1200  1250  1300  1350  CDP  Figure 4.9: Normalized velocities output from the seismic impedance inversion. The sand aquifer can be clearly identified by the yellow and red colors. 83  1050  1100  1150  1200  1250  1300  1350  CDP  Figure 4.10: Lithology/void ratio profile obtained by separating the sand and clay zones using a cutoff velocity of 170 m/s and using Equation 4.7 to obtain void ratio within the sand aquifer.  84  CHAPTER 5 Aquifer Hydraulic Conductivity Using Cone Penetrometer and Seismic Data 5.1  Introduction A n objective of hydrogeological investigations is the prediction of fluid flow within  subsurface aquifers.  This objective is achieved by creating a model of the subsurface that is  constrained by a variety of different data. The data are obtained by slug tests, pumping tests and core sample analysis coupled with geologic mapping, hydrostratigraphic unit classification and tracer analysis. The measurement scale, by necessity is always much smaller than the aquifer scale.  The accuracy of the results will depend on how representative the measurements are of  the entire aquifer.  The more complex the aquifer is, the more challenging it will be to detail the  complexities within the aquifer.  Geophysical methods provide a way to detail the aquifer  complexities using non-invasive, surface-based measurements. There are different types of geophysical methods capable of imaging the earth. M y focus has been on the seismic method, which involves the propagation of seismic waves. A wave technique has the potential to provide excellent subsurface detail because the measured response can be attributed to specific regions of the subsurface.  Seismic waves respond to changes in the  mechanical properties of the earth (ie. compressional and shear modulus and density) and allow for the determination of these properties in great detail.  Seismic waves are generated at the  surface using some type of source and as they encounter discontinuities in the mechanical properties of the subsurface they are reflected back to the surface where they are measured by a series of receivers (known as geophones).  The amplitude of the reflected response is due to the  magnitude of the change in mechanical properties. A n inversion of the seismic amplitudes can thus be used to map out the property variations.  85  The seismic method can employ two different types of waves: shear and compressional. Compressional waves are directly analogous to sound waves or oscillating pressure pulses, while shear waves involve an oscillating shearing motion and as a result can only propagate through solids. The velocity of the waves is directly related to the compressional or shear modulus and inversely related to the density. The resolving ability of the seismic method is a function of the length of the propagating seismic pulse or wavelength.  A generally accepted rule is 1/3 of a  wavelength for the smallest bed that can be imaged (Miller et al., 1995). In saturated sediments shear-waves will have a smaller wavelength than compressional-waves and thereby provide the best resolution (Dobecki, 1993).  This is one of the key reasons I have chosen to use a shear-  wave seismic method for our investigation. The source and receiver configurations determine the type of data that are recorded and control the type of information that can be obtained.  A surface reflection seismic survey  involves having both source and receivers on the surface typically in a straight-line arrangement. The recorded data consists of a series of records showing the variation in amplitudes of reflected waves. A tomographic survey records the times of the first arriving waves, which are then used to invert for velocity.  One common tomographic survey is also called a refraction survey with  the sources and receivers both on the surface. boreholes, which improves the  The sources and receivers can also be placed in  "spread" of the propagating waves and when the data are  inverted will provide a better constrained velocity model. I have chosen to use the surface based reflection seismic survey due to its non-invasive nature and reduced cost compared to the drilling of boreholes and deployment of borehole based equipment. Hydrogeologists have recognized the value of having seismic properties for groundwater modeling.  In particular, the velocity of the seismic waves has been demonstrated to be a  valuable tool.  Research by Copty et al. (1993) shows that a profile of velocity variations can  provide important constraints for the inversion of hydrogeological properties such as 86  permeability.  Some type of relationship between seismic velocity and permeability must be  used. Copty et al. (1993) used laboratory measurements by Marion et al. (1992) that demonstrate the non-unique variation of velocity with permeability given changes in clay content within a uniformly graded sand. The non-uniqueness can also be handled to some degree by having other types of data available such as tracer data (e.g. Hyndman and Gorelick, 1996). obtaining the velocity profiles is only addressed by some of the researchers.  Actually  One technique  discussed by McKenna and Poeter (1995) uses seismic tomographic data from adjacent boreholes to provide the velocity data.  In Chapter 4 I have detailed how surface seismic data  can be used to obtain velocity information. Chapter 4 demonstrates how shear-wave velocities were obtained from the amplitudes of a surface based shear-wave reflection seismic survey.  A Bayesian inversion technique was  utilized which used the known velocities obtained by a seismic cone penetrometer as input to a starting model and resulted in an optimized determination of shear-wave velocities along the entire seismic profile.  The shear-wave velocities were then converted to porosity based on  laboratory measurements.  This chapter is an extension of Chapter 4. The Bayesian technique is  once again used to provide a profile of shear-wave velocity. In addition, the cone penetrometer tip resistance is integrated as a way to provide a connection with both grain size and hydraulic conductivity. The cone penetrometer is a commonly used technique by geotechnical engineers for rapid subsurface characterization (see Lunne et al., 1997). It involves the pushing of a steel rod into the ground and making measurements as the rod is pushed using sensors mounted close to the conical tip. One of the most common measurements is tip resistance. The tip resistance involves measuring the force required to push the rod into the ground. As the rod is pushed the sediments have to move away from the tip and are compressed into the surrounding sediments. The tip resistance will be closely related to both the shear modulus (from grains sliding against each 87  other) and the compressional modulus (from grains being compressed into surrounding sediments).  In addition, there will be a grain size dependence because the smaller grains are  easier to push aside than larger grains.  This grain size dependence is apparent when sediment  classification schemes relating tip resistance to sediment type (as illustrated by Lunne et al., 1997) are examined.  The local stress conditions (both horizontal and vertical) also will be  significant because this stress must be overcome in order to move the grains.  Due to all the  factors that influence tip resistance a simple, globally applicable relationship cannot be expected between tip resistance and either grain size or shear-wave velocity. However, it is possible to develop a site-specific, empirical relationship between tip resistance and both shear-wave velocity and grain size. The reliable determination of grain size is the key to obtaining values of hydraulic conductivity. The grain size distribution will influence both the size and the tortuosity of the pore spaces, parameters that have a direct effect on how easily fluids can move through a porous medium.  This is recognized empirically by the Hazen formula (see Freeze and Cherry, 1979),  which relates hydraulic conductivity (K) to an effective grain size (dio): K=A[d ] .  (5.1)  2  10  The value of dio is taken from a sieve analysis of the sands from an aquifer.  A sieve analysis  provides the data for grain size distribution curves and dio is the point on the curves where 10% (by weight) of the sand is finer and 90% of the sand is coarser. The constant, A , is equal to 100 x 10" m/mm s. The Hazen formula seems to work best in uncemented sand aquifers with very 4  2  little clay content. Wood (1996) has demonstrated how a sieve analysis of the sands from core taken at the Kidd2 site resulted in reliable estimates of the hydraulic conductivity (using the Hazen formula). The sieve analysis data from Wood (1996) were used in this chapter for the determination of a relationship between tip resistance and grain size.  88  In this chapter I have detailed the development of a site-specific relationship between shearwave velocity and tip resistance and between tip resistance and grain size. As a result, the tip resistance was used to develop the constraining velocity model for the inversion of the shearwave seismic data and the refined velocity model output from the inversion was then converted to a model of tip resistance.  The ultimate objective of obtaining a model of hydraulic  conductivity was achieved by converting the tip resistance to grain size and then to hydraulic conductivity.  5.2  Geologic Setting This research was also carried out at the Kidd2 hydrogeology research site. Figure 5.1 shows  the location of the site and the relative locations of the seismic lines, cone penetrometer holes and the core hole used in the investigation. In terms of hydrogeology the sediments can be classified as a three layer confined aquifer system. The aquitard above the aquifer consists of fine-grained silty clay and peat. The aquifer consists of fine to medium sand with sub-angular to sub-rounded grains that have been reworked by wave and tidal action during deposition (Wood, 1996).  The aquitard at the base of the  aquifer consists of fine-grained clayey-silt that was most likely deposited in a deep-water marine setting. Underlying the entire sequence is the Pleistocene till, which does not influence the hydrogeology within the aquifer but has a large change in mechanical properties and can readily be identified on the seismic data. Neilson-Welch (1999) summarizes the measurements of hydraulic conductivity that have been measured at the site using a variety of techniques. The hydraulic conductivity of the lower aquifer ranges from 1 x 1CT m/s to 4 x 10" m/s and the upper aquifer has a similar range but 4  4  with values as low as 2 xlO" m/s in the siltier regions (typically at depths less than 8 m). The 5  89  clay-dominated aquitard between the aquifer and the Pleistocene till has a hydraulic conductivity of 1 x 10~ m/s. 8  5.3  Cone Penetrometer Data  The cone penetrometer is an effective tool for the examination and characterization of aquifers.  The technique is fast and relatively inexpensive and provides detailed information for  the classification of subsurface sediments.  The cone was developed by geotechnical engineers  for in situ site characterization. Over the years it has proven to have numerous applications to hydrogeology including permeability estimation and groundwater sampling (see Lunne et al., 1997).  The technique relies on having the ability to push a steel rod hydraulically into the  sediments so it can only be used where this can be achieved; fine-grained deltaic sediments are ideal whereas gravelly deposits can be problematic. I have used the cone penetrometer for the acquisition of velocity and tip resistance data.  5.3.1  Data Acquisition  The cone penetrometer has been used extensively in a number of other studies at the Kidd2 site. The cone holes with tip resistance data and locations close to the seismic data are shown in Figure 5.1 as K9404, K9505 and K9603. Velocity and tip resistance data have been acquired in only two holes ( K l and K2 in Figure 5.1). The shear-wave velocity surveys ( K l and K2), were conducted along line 20, located as shown in Figure 5.1, using the cone penetrometer truck from the Department of Civil Engineering at the University of British Columbia. A large sledgehammer against the baseplate of the truck was used as a source of shear-waves and the receiver was mounted close to the conical tip. The data were recorded using a Nicolet oscilloscope with 15 bit data digitizing and storage capability. The velocities were obtained by determining the arrival times of the seismic 90  waves at adjacent depth levels. A ratio of the depth interval to the difference in arrival time gives the velocity over that interval.  5.3.2  Tip Resistance and Shear-Wave Velocity Relationship  A critical input to the impedance inversion algorithm is a reliable set of shear-wave velocities. These velocities can be obtained directly using the seismic cone penetrometer but the acquisition of good-quality shear-wave velocity data can be time consuming and expensive compared to the acquisition of cone logs such as tip resistance.  At a particular site, there may be  a historical set of tip resistance logs but no corresponding velocity logs, making it advantageous to use the existing cone information. The tip resistance logs can be used to obtain velocities i f a relationship between the two parameters is known. Other researchers (Olsen, 1994; Lunne et a l , 1997) have also recognized the benefits of having a conversion function between tip resistance and shear-wave velocity. I have developed a site-specific relationship based on available cone data.  The development of a relationship between shear-wave velocity and tip resistance is  complicated by the influence of effective stress, which must first be removed.  5.3.2.1  Removing the effect of effective stress  In the near-surface the effective stress increases very rapidly in the upper 20 m.  The  increasing effective stress will affect the measurement of any parameter sensitive to horizontal or vertical stresses.  For the cone penetrometer, this means the tip resistance and shear-wave  velocities will both be affected.  This was recognized by geotechnical engineers (Robertson et  al., 1992; Olsen, 1994) who were interested in identifying distinct lithologies based on the cone log responses.  91  The series of normalization functions recommended by Robertson et al. (1992) to remove the stress effect from tip resistance (q ) and shear-wave velocity (V ) are: t  s  Qtn=(qt-o"vo)/ o-' o,  (5.2)  V  V =V [P /o-' f .  (5.3)  25  sn  s  a  vo  The normalized tip reistance (Q ) is a function of the vertical stress (a ) and the vertical tn  vo  effective stress (a' ). The normalized shear-wave velocity ( V ) is a function of the atmospheric vo  sn  pressure (P ) and the vertical effective stress. The problem with these functions, in particular the a  normalization of q , is that different lithologies show the stress effect in different ways. Equation t  5.2 works well for clay rich sediments but in sands the shallower tip resistances are overcorrected. A uniform sand layer would show a decreasing normalized tip resistance with depth. This was recognized by Olsen (1994) who recommended evaluating the stress exponent for every different lithology. However, any type of correction that is lithology dependent requires some type of interpretation in advance, which defeats the automated, non-biased approach favored by geotechnical engineers (Lunne et al., 1997).  M y goal was to create a simple  normalization function that properly compensated for the effective stress at our site, and in particular, that properly dealt with the stress effect within the sand aquifer, which is the primary zone of interest. The normalization functions I used were: Qtn=(qt-avo)(Pa/a vo)°' ,  (5.4)  V =V [P /(a' +k)] ,  (5.5)  ,  5  a3  sn  s  a  vo  where P is the atmospheric pressure (~100kPa) and k is a constant of lOkPa that compensates a  for the non-linear nature of the shear-wave velocity profile at the earth's surface.  I first  presented the velocity normalization equation in Chapter 4. Geotechnical engineers (see Lunne et al., 1997) define the relative density of sands using a very similar relationship to Equation 5.4 92  and they also use a normalization exponent of 0.5.  Defining a constant relative density with  depth is directly analogous to the effective stress compensation.  With the effective stress  removed from the tip resistance and shear-wave velocity the data were plotted to determine a relationship between the two parameters. This relationship enabled me to transform between shear-wave velocity and tip resistance, a key step in the determination of the hydraulic conductivity model.  5.3.2.2  Normalized tip resistance and shear-wave velocity plot  The goal of this research was to estimate hydraulic conductivities from cone penetrometer and seismic data.  M y approach involved the development of a relationship between tip  resistance and shear-wave velocity.  Both tip resistance and shear-wave velocity have been  acquired from the same cone penetrometer holes in two locations ( K l and K2) within the same aquifer system. The influence of effective stress has been removed from the tip resistance and shear-wave velocity curves so that the responses could be attributed to variations within the sand aquifer. The next step was to plot the normalized values to determine a relationship between the two parameters. Before the data were plotted, the vertical sampling of velocity and tip resistance had to be equalized.  The tip resistance down a cone hole is measured at 2.5 cm intervals whereas the  shear-wave velocity sampling varied from 0.33 m to 2.0 m.  The shear-wave velocity was  assumed to be an average of the shear-wave velocities over the range of the measurement interval. The tip resistances were averaged over the same intervals. The plot of normalized shear-wave velocity vs normalized tip resistance is shown in Figure 5.2. A straight line fit to the data is also shown in Figure 5.2 resulting in the equation:  93  V =2.82Q +158. sn  (5.6)  tn  The standard error in the fit was 8.4 m/s.  Other researchers (e.g. Lunne et al., 1997) have  characterized the relationship between normalized tip resistance (Qm) and shear-wave velocity (Vs„)  as Qtn=[V /Y] ,  (5.7)  4  s n  where Y represents a constant that is site-specific and depends on the sand grain characteristics such as sand compressibility, age and degree of cementation (Lunne et al., 1997). not easily support this non-linear relationship that applies only to sands.  M y data do  The velocities  computed using this relationship for a clay layer would be completely unrealistic (too low). M y intent was to use a relationship that defines the velocities within sand and clay because both will be required to develop the velocity model. I decided that the linear relationship was the simplest and most appropriate for the tip resistance and velocity data from the Kidd2 site. The standard error in the linear fit was 8.4 m/s, which attributes all of the error in the fit to the shear-wave velocities. This represents a small relative error in the shear-wave velocities of approximately 5%. The ability to determine shear-wave velocities from tip resistance to within 5% indicates that the linear relationship is an excellent way to estimate shear-wave velocities from tip resistance measurements at this site. To look in detail at the difference between measured and predicted velocity data the tip resistance data from the two cone holes ( K l and K2) were transformed to velocity and plotted along with the cone determined shear-wave velocities. Figure 5.3 is a plot of the two velocity profiles.  The figure demonstrates a very close match between the estimated and actual  velocities. The derived velocities within the clay layer from a depth of 22 to 50 m show none of the variations obvious with the measured velocity data; a result of the lack of sensitivity of the tip resistance to clay dominated sediments. 94  The linear relationship between normalized tip resistance and normalized shear-wave velocity developed for this site has proven to be a useful way to obtain reliable shear-wave velocities.  The development of a relationship has the added benefit of providing a way to  convert shear-wave velocity to cone tip resistance, a relationship that will be critical in the determination of hydraulic conductivity. The next step was to look at using the derived shearwave velocities to constrain the inversion of surface seismic data.  5.4  SH-Wave Seismic Reflection Data The common depth point (CDP) seismic reflection method is a non-invasive geophysical  technique that has the potential to supply a significant amount of subsurface information. The technique is based on sending seismic waves into the ground and sensing the vibrations on the surface as they reflect off of subsurface changes in seismic velocity and density.  The CDP  technique relies on having a large number of sensors (geophones) laid out on the surface for every source activation. The overlap between subsequent recordings provides redundancy in the subsurface sampling that is later used in the processing to enhance the signal compared to the noise.  5.4.1  Data Acquisition  The data were acquired using a 24 channel Geometries SmartSeis seismograph, which is capable of recording data from 24 geophones simultaneously.  In order to reduce the time  required to collect the data a roll box was interfaced with the seismograph that allowed for the placement of 48 geophones (Geo Space 20DM - 28 Hz) on the ground and the recording of data from 24 geophones.  The time sample interval was 0.0005 s and a 600 Hz high-cut filter was  95  applied to the analog data before digitization to prevent aliasing. Shear-wave CDP reflection data were collected; the shear-waves being generated by a horizontal shearing of the ground surface. The source of the shear-waves was a 35 kg steel block that was impacted on the side. Horizontal geophones placed along the survey lines sensed the horizontal motions associated with the propagating shear-waves. A total of 5 lines of CDP shear-wave data were recorded in June 1996 at the Kidd2 research site. Line 20 is the only line examined in this paper, primarily because there are two cone holes with velocity data on the line and two other cone holes with tip resistance data on the line. The velocity data were used in Chapter 4 to constrain an impedance inversion of line 20. M y intent was to build on this earlier research by demonstrating how tip resistance could be used for two key operations: obtaining the shear-wave velocity constraints for the impedance inversion and providing a way to convert the inverted velocities to hydraulic conductivity.  5.4.2  Data Processing  The seismic data processing has been discussed in both Chapter 3 and Chapter 4. generally accepted processing flow was followed (e.g. Yilmaz, 1987). seismic section is shown in Figure 5.4.  A  The final migrated  The dataset consists of 192 traces with a displayed  recording time of 0.6 s. Each trace in the seismic section can be thought of as being the result of a seismic experiment with both source and a single receiver in the same location. The amplitude variations along each trace correspond to a combination of background noise and changes in velocity or density in the subsurface. The size of the amplitude variation is indicative of the magnitude of the subsurface variation.  96  The location of the aquifer is from a time of 0.08 to 0.28 s. Within this time zone there are a number of reflective surfaces; some are short and horizontal others are long and dipping. The time zone from 0.28 to 0.48 s represents the clay-dominated zone.  There are few lithologic  variations expected within this layer and as a result there are very few reflections visible on the seismic data. The top of the Pleistocene till is represented by the high amplitude reflection at a time of 0.48 s.  The velocity contrast between the clay and the till is significant therefore  resulting in a large amplitude reflection. Looking at a migrated seismic dataset (such as Figure 5.4), it is possible to qualitatively assess the location and structure of the different lithologic units.  However, my goal was to  obtain quantitative information about the sand aquifer. The seismic amplitudes are due solely to changes in velocity and density so this was the only information I could hope to extract from the data. A seismic impedance inversion was used to achieve this goal.  5.5  Impedance Inversion of Shear-Wave Data The methodology used for the seismic impedance inversion is the same as that used in  Chapter 4.  The impedance inversion uses a Bayesian approach that allows a priori information  and associated errors to be incorporated into the inversion. The Bayesian inversion algorithm was developed at the University of British Columbia by M . D . Sacchi and T J . Ulrych and is detailed in Sacchi and Ulrych (1996).  Seismic impedance is the product of velocity and density.  I have assumed the density is relatively constant, so the seismic impedance model is simply a velocity model. For a shear-wave seismic inversion the velocities that are relevant are shearwave velocities.  The shear-wave velocities for this inversion were derived from cone  penetrometer tip resistance logs.  97  I have not directly used the known velocity data (from K l and K2) in the development of the input model. Ignoring the velocity data has three main purposes. M y technique illustrates how tip resistance alone can successfully be used to develop a velocity model. At a particular site the velocity and tip resistance relationship can be developed at one location and then applied to the tip resistance logs falling close to the seismic lines.  Secondly, I will eventually be converting  the output velocity model back to tip resistance as a way of determining hydraulic conductivity. A n input model developed strictly from tip resistance will be laterally consistent. The third purpose is related to the detail provided by tip resistance compared to velocity logs.  The  sampling of tip resistance (typically 2.5 cm) is much finer than the sampling of velocities (typically 0.5 to 1.0 m) resulting in a very detailed velocity model that accurately locates significant impedance boundaries. The normalized input velocity model is shown in Figure 5.5. The model is based on four control points that are labeled on the migrated seismic data in Figure 5.4. The starting point of the model is the tip resistance data from each of the control points. Both K l and K 2 have tip resistance data down to the top of the Pleistocene till.  Using Equation 5.6 the normalized  velocity for each of these control points was computed. The velocity of the Pleistocene till had to be estimated because the cone could not penetrate into the Pleistocene till. Based on work by Hunter et al. (1998a), the velocity of the Pleistocene till is at least double the velocity of the overlying deltaic sediments; the model reflects this large increase in velocity. The control point associated with K9603 had tip resistance data that ended at a depth of 25 m, within the clay layer beneath the aquifer. The cone tip resistance down to the Pleistocene till had to be synthesized.  The synthesizing of the cone logs was necessary because the  interpolation algorithm used to create the input model required full-length velocity logs. The tip resistance of the clay does not change very much because the lithology is uniform (as with K l and K2) so material with a similarly low tip resistance was added to the bottom to a depth 98  corresponding to the known reflection time of the Pleistocene till on the shear-wave seismic data. The tip resistance values were then converted to normalized velocities using Equation 5.6. The velocities of the Pleistocene till are the same as the Pleistocene till velocities used with K l and K2. The control point associated with K9404 had tip resistance data that ended at a depth of 21 m (within the aquifer but close to the base), thereby requiring the addition of synthesized tip resistance data for the aquifer, underlying clay and Pleistocene till. The tip resistance values from the nearby K 2 cone hole were averaged over the missing depth interval (21 to 22 m). A constant 0.8 m block of aquifer material (with this average tip resistance) was then added to the bottom of the K9404 tip resistance log. Finally, the underlying aquitard thickness was estimated and an appropriate amount of low tip resistance values were added to the base of the synthesized log to tie the known reflection time of the Pleistocene till. A comparison of the input and output velocity models is shown in Figure 5.6. The four sets of velocity traces represent the location of the four cone holes, where the input velocity model was derived from the tip resistance logs. The line with triangles is the input velocity model and the solid line is the output from the inversion. The trends in the velocities are very similar thereby demonstrating the viability of the inversion. The output velocity model from the inversion is shown in Figure 5.7. The velocity variation is similar in form but differs in the detail when compared to the input model. The continuity of the velocities in the input model is due to the interpolation of four control points. The inversion process uses both the input velocity model as well as the input seismic amplitudes. The output model does not show this same continuity because the changes in seismic amplitudes have been honored by the inversion. M y working hypothesis is that the variation in seismic velocity of the subsurface can be used to map out the hydraulic conductivity of the aquifer. In Figure 5.7 I have obtained, through 99  inversion of our collected surface seismic data, a model of seismic velocity. The next step was to use the relationship between tip resistance and velocity and between tip resistance and grain size to finally establish a model of hydraulic conductivity.  5.6  Aquifer Hydraulic Conductivity A desirable objective in using seismic data for hydrogeologic studies is to be able to relate  the measured velocity to hydraulic conductivity.  The relationship  conductivity and velocity, however, is not straightforward.  between hydraulic  The use of laboratory-derived  relationships, such as that presented by Marion et al. (1992), requires key assumptions about the microstructure of the in situ material; e.g. a sand framework with constant intergranular void ratio. Alternatively some means of calibrating field-measured velocity data with field-measured hydraulic conductivity is required; as in the study by McKenna and Poeter (1995).  In this  chapter, I have used the relationship developed between tip resistance and velocity to convert the velocity model to a tip resistance model. The grain size distribution curves and nearby cone tip resistance data were used to relate tip resistance to grain size resulting in a model of sand grain size variations. Finally, using the Hazen formula (Freeze and Cherry, 1979) the sand grain size model was converted to hydraulic conductivity.  5.6.1  Tip Resistance Model of Aquifer  M y objective was to obtain a model of hydraulic conductivity using the Hazen formula which relates hydraulic conductivity and grain size. In order to accomplish this task I required a model in a form that could be readily and uniquely converted. The tip resistance showed significant promise due to previous implied relationships between grain size and tip resistance (see Lunne et al., 1997).  I have already developed a relationship between shear-wave velocity and tip  100  resistance, so the conversion of velocity to tip resistance was relatively straightforward. The tip resistance is given by the rearrangement of Equation 5.6, so that Q =3.542V -559.7. tn  (5.8)  sn  The application of Equation 5.8 to the normalized shear-wave velocity model of Figure 5.6 resulted in the tip resistance model shown in Figure 5.8. The simple linear relationship that gave us the tip resistance suggested that there would be very few differences in the appearance of the tip resistance plot compared to the normalized velocity plot (Figure 5.6); zones of high shearwave velocity mapped directly to zones of high tip resistance. As expected, the clay zones mapped out with a very low tip resistance. The variations within the clay zone are due to the variations introduced by the inversion program in order to honor the amplitude variations of the input seismic data (Figure 5.4). To demonstrate the quality of the derived tip resistance data I compared the derived tip resistance with the original four tip resistance profiles that I started with. The comparison of the four profiles is shown in Figure 5.9.  Both K l and K 2 exhibit close matches between the tip  resistance derived from the impedance model and the original cone tip resistance. The low tip resistance core of the aquifer is matched in both cases. The largest errors appear at the base of the aquifer. K9603 and K9404 are not as well matched, but the fit is still relatively good. The impedance inversion has provided an effective technique for the interpolation of the cone tip resistance data.  The final step involved the conversion of tip resistance to grain size, a  parameter that is readily related to hydraulic conductivity.  5.6.2  Tip Resistance and Grain Size Relationship  It is generally accepted that the cone penetrometer tip resistance is directly related to grain size variations (see Lunne et al., 1997). This relationship was expected to be very complex and due to many factors including the grain roughness and the stress history of a particular deposit. 101  Recognizing the potential complexity I decided to establish a site-specific relationship similar to the way I analyzed tip resistance and shear-wave velocity. At the Kidd2 site I had both core (labeled as C O R E in Figure 5.1) with measured grain size distribution curves (Wood, 1996) and a nearby cone penetrometer hole (labeled as K9505 in Figure 5.1) with tip resistance.  From  these data I developed the tip resistance and grain size relationship. The dsn grain size versus tip resistance data are illustrated in Figure 5.10. D  5 0  represents the  point on a grain size distribution curve where 50% of the sand (by weight) is finer and 50%> is coarser (also referred to as the median grain size). I used a log-log plot because grain size variations are typically log-normally distributed. The grain sizes were based on representative samples taken from distinct lithologic sections of core that were either 2.5 feet (0.76 m) or 5 feet (1.52 m) in length (standard core barrel lengths). In order to match the tip resistance data to the measured grain sizes, the tip resistances were averaged over the same depth intervals. The error in grain size is given by half the sieve size increment and the error in tip resistance is given by the standard deviation of the averaged values.  A least squares line was fit to the data and is  plotted in Figure 5.10. The equation of the line is:  Lo o[d5o]=0.56Log [Qtn]-l.l gl  (5.9)  10  with dso in mm and Q  t n  in MPa. Equation 5.9 provides a simple way to determine the grain size  distribution within the aquifer. I used the d n values from the Wood (1996) grain size distribution curves to develop the 5  relationship between normalized tip resistance and grain size. The dso values were used because the dio values were not determined by Wood (1996) for the finer grained sediments (finer sieves were required).  However, the Hazen constant (A) assumes that dio values are used.  An  examination of the dio and dso grain size distribution curves revealed that the d o sand grain size 5  102  was approximately double the dio. As a result, a constant of 25 x 10" m/mm s was used in the 4  2  Hazen formula to make the d.50 grain size valid for the computation of hydraulic conductivity.  5.6.3  Hydraulic Conductivity Model of Aquifer  The primary objective of this study was to develop a means of using surface seismic data to delineate the hydraulic conductivity of a sand-dominated aquifer.  The profile of shear-wave  velocity output from the impedance inversion provided some insight into the variations in the aquifer.  Rather than consider these variations in velocity, they have been transformed to tip  resistance, from which the grain size variations were computed using Equation 5.9. The final step involved the computation of a hydraulic conductivity model. The conversion between sand grain size and hydraulic conductivity uses the Hazen formula (Equation 5.1, see Freeze and Cherry, 1979). Wood (1996) demonstrated the effectiveness of using the Hazen formula at the Kidd2 site to convert sand grain size to hydraulic conductivity. The results compared favorably with values of hydraulic conductivity measured using a variety of other techniques such as slug tests and pump tests. The Hazen formula and the conversion function between tip resistance and dso grain size were combined into a single relationship:  Lo o[K]=1.13Log [Qtn]-5.948 gl  (5.10)  10  Equation 5.10 converted the tip resistance model directly to a logarithmic hydraulic conductivity model. The final hydraulic conductivity model is shown in Figure 5.11. A cutoff of 5 M P a was chosen as a lower bound for the validity of Equation 5.10.  Zones of the aquifer with a  normalized tip resistance less than 5 M P a were assumed to be clay-dominated and assigned a hydraulic conductivity of 1 x 10" m/s. 8  The aquifer hydraulic conductivities range  from  1.0 x 10" m/s to 3.0 x 10" m/s. These values are in excellent agreement with the values of 4  4  103  hydraulic conductivity determined by Wood (1996) as well as values from other independent measurements (see Neilson-Welch, 1999) not dependent on a grain size relationship. The general pattern of hydraulic conductivity variations in Figure 5.11 matches the expected variations extremely well. The highest hydraulic conductivities are expected at the base of the two fining upwards sequences that make up the aquifer (i.e. within the coarser grained sands). The band of lower hydraulic conductivity at a depth of 14 m defines the top of the lower sequence and the base of the upper sequence. The magnitude of the variation (as seen in Figure 5.11) between the high and low conductivity zones within each sequence is approximately 0.1 logio m/s. This is a very reasonable value that matches the hydraulic conductivity variation measured by Wood (1996). In order to quantify the potential error in the derived hydraulic conductivities two additional inversions were done. In the first inversion the input aquifer velocities were increased uniformly by 8 m/s, corresponding to the standard error in the fit of normalized shear-wave velocity to normalized tip resistance (Figure 5.2). In the second inversion the input aquifer velocities were decreased uniformly by 8 m/s. The outputs from these two inversions were then converted to hydraulic conductivity using the same assumptions and equations used for the model in Figure 5.11. The two hydraulic conductivity models were then subtracted to approximate the potential error. The potential error plot is shown in Figure 5.12. Most of the potential error in hydraulic conductivity is within the sand aquifer. The errors in hydraulic conductivity were expected to fall within the aquifer zone because the only input velocities that differed in the two inversions were within the aquifer. A n exception is the spreading out of the error at the base of the aquifer (into the underlying aquitard), which is a result of the low amplitude and coherency associated with the reflection at the base of the aquifer. The error within the aquifer has a maximum value of 0.14 logio m/s and a typical value of 0.10 login m/s. The error is largest within the center of 104  the aquifer and tapers to zero at the top and base. The largest errors occur within narrow, horizontal bands at a depth of approximately 10 m. These bands correspond to high amplitude reflections on the input seismic data (Figure 5.4) and are associated with features that do not have cone penetrometer velocity or tip resistance control. The general error pattern shows the largest error furthest away from constraining velocity variations (i.e. the top and bottom of the aquifer). This results from one of the limitations of seismic inversion. The seismic amplitudes alone do not contain the low frequency data required to identify larger scale changes in velocity with depth. The inversion relies on obtaining this data from the input, constraining velocity model. The inversion is most valuable for identifying smaller-scale relative variations in velocity and hydraulic conductivity such as the variation between the top and bottom of a fining upwards sequence. The inversion of the seismic reflection data has been proven to be a powerful tool for mapping variations in hydraulic conductivity.  The absolute values match other independent  measurements of hydraulic conductivity extremely well. The relative variations also agree with the known sedimentary textures and measured values.  The inversion relies on having good  velocity control. Significant variations in amplitude on the input seismic reflection data require velocity control using borehole based or cone penetrometer measurements. The conversion of velocity to hydraulic conductivity based on empirical measurements will have error associated with each step. Ideally, these steps should be minimized so that velocity is directly related to in situ measurements of hydraulic conductivity. When in situ hydraulic conductivity data are not available then this research demonstrates how a viable map of hydraulic conductivity can still be obtained.  105  5.7  Conclusions The cone penetrometer is a valuable tool for taking seismic velocities determined from near-  surface seismic reflection data and converting them to values of hydraulic conductivity. The cone penetrometer has the advantage of rapid deployment, reduced cost and the simultaneous acquisition of a suite of logs for stratigraphic interpretation and the development of wellconstrained site-specific relationships. The challenge involves identifying a link between conductivity.  seismic velocity and hydraulic  This link must be associated with rock properties such as grain size that have a  reasonably well-understood connection with hydraulic conductivity.  Cone penetrometer  measurements such as tip resistance are not always going to have an easily identified relationship with velocity.  A normalization procedure, which removes the effect of effective stress, is a  critical part of developing a well-constrained relationship.  The use of a relationship derived  from the work of other researchers or even from a different site may not be applicable to the site being investigated.  Certainly, a relationship from another site should first be tested using local  data, or ideally, a new site-specific relationship should be developed to obtain the bestconstrained values of hydraulic conductivity.  106  Figure 5.1: Location map showing the Kidd2 research site on the Fraser River delta with the grid of CDP seismic (dotted lines) and the cone holes that were used to constrain the inversion ( K l , K2, K9404, K9603). The seismic line (line 20) used in the inversion has been highlighted as well as the location of the borehole with core (CORE).  240  220  200  S  c  180  160  •  •  •  ~  • 4  • *  • • • •  i  •• •  V  •  R =0.71 2  140  120 10  12  14  16  18  Qm (MPa)  Figure 5.2: Plot of normalized shear-wave velocity (V ) vs normalized tip resistance (Q ). The line represents a least squares fit of the data. sn  107  tn  300  250  T  •  200 <* :<  a  •  Jf  150  100  50  10  15  20  25  30  35  40  45  50  Depth (m)  300  250  t *  200  •  •  *  •  •  •  >§, 150  100 •  50  0  5  10  15  20  25  30  35  40  45  50  Depth (m) Figure 5.3: Plots showing a comparison between the tip-resistance-derived shear-wave velocity (solid line) and the cone-measured shear-wave velocity (dashed line with diamonds) for a) K l andb) K2. 108  K2 l O l O I  1 O 5 0 1  1 0 9 0 I  K9404 1 1 3 0 1  CDP 1 1 7 0  !  K1 1 2 1 0  !  K9603 1 2 5 0  I  1 2 9 0  !  1 3 3 0  1 3 7 0  I.  Figure 5.4: The migrated SH-wave seismic reflection line 20 that was used as the input to the seismic impedance inversion. The four cone holes that were used to constrain the impedance inversion are highlighted.  CDP  Figure 5.5: Normalized velocities input to the seismic impedance inversion. 109  o  o  CO  a  in  o •4—»  o  o  m co  m CO  & o  oo  o  o  o  co  £ cn  co CN  o  CU  Q  O  Q  CN  o M > -a <D rt5 o  CN  ON cn CD  _ o  (3 rt  o ^ o  o  o  o  o  o  o  o  o  o  o  m m — i c \ r — m c O r t O \ t ^ i n C N C N C N - H — i — . r t ^ H  o  o  o  o  o  o  o  o  o  o  'cn I-H CD  o  •S .2  • n r o - — 0 \ t ^ i n c o > - < o \ r ^ i n C N c N C N ^ r t r t r t r t  (s/ui) S A  . .  on  CD o3  (S/Ul) S A  *0  CD CD —I  o  ^ CD CD  m  rt rt  o  •S £  o o  cn  3 O 3 <4H PH CD  CO  O cn  CN  rt  2  <L>  E .3 CD  in cN  o  00  w  > c  Q  CD O rt3 • — TJ3  O co  CO  >  C S ""' O CD 00  o  'C 1,3  g  is  g rt, M*5 I/) CU iO O O O O O O O O O O i n n > — i O \ r ^ i n n — ' O s t e i n CN CN CN " rt ! rt  O O O O O O O O O O O i n c o ^ o s r ^ i n m ^ o \ t ^ i n  rt  (s/iu) S A  (s/iu) S A  110  s  OJD  1050  1100  1150  1200  1250  1300  1350  CDP  Figure 5.7: Normalized shear-wave velocities output from the impedance inversion.  1050  1100  1150  1200  1250  1300  1350  CDP  Figure 5.8: Normalized tip resistance computed from the normalized velocities in Figure 5.7. Ill  -1.2  0.2  0.4  0.6  0.8  1.2  1.4  Log, [Qta] 0  Figure 5.10: Plot of grain size (d ) versus normalized tip resistance ( Q J . The line represents a least squares fit of the data. 50  1050  1100  1150  1200  1250  1300  1350  CDP  Figure 5.11: Color plot showing the hydraulic conductivities computed from the normalized shear wave velocities output from the impedance inversion. The contour interval is0.05. 113  1050  1100  1150  1200  1250 1300 1350  CDP  Figure 5.12: Color plot showing the potential error in the computed hydraulic conductivities. The plot was generated by running two inversions with aquifer velocities of+/-8 m/s with respect to the original inversion. The hydraulic conductivities were then calculated and the difference between the two results is illustrated in the plot. The contour interval is 0.02.  114  CHAPTER 6 Conclusions and Recommendations for Future Work  This thesis has clearly identified the benefits of seismic reflection techniques for hydrogeological applications. P-wave and S-wave seismic reflection techniques were used for locating aquifer boundaries and the S-wave method was used for mapping the distribution of hydrogeologic properties. The P-wave technique works best when the near-surface sediments are fine-grained and saturated.  The S-wave technique also works best when near-surface  sediments are saturated but additionally requires that a high-velocity zone (relative to the underlying sediments) exists at the surface. The cone penetrometer has been demonstrated to have significant applications for hydrogeological investigations.  The cone can be used for V S P surveys to provide high-  resolution images of aquifer boundaries. The cone VSPs also play a critical role in the utilization of surface seismic reflection data for hydrogeologic applications. The seismic reflection data alone can only provide information on variations in velocity and density, providing some idea about subsurface heterogeneity. Integrating the cone penetrometer (e.g. tip resistance) data with the seismic data provides the ability to obtain additional hydrogeologic parameters such as hydraulic conductivity.  The basis of linking the cone penetrometer to seismic data and  hydrogeologic properties is the generation of site-specific relationships between the parameters. Site-specific relationships are required because the local stress history and sediment types cannot be easily constrained by globally derived relationships. The inversion of seismic reflection data is a powerful way to change seismic amplitudes into variations of subsurface velocity. The Bayesian technique used in this research is extremely robust and simple to use. The velocity variations output from the inversion can be used both to  115  quantify heterogeneity and to provide data for developing relationships with geotechnical parameters such as tip resistance and hydrogeologic properties such as hydraulic conductivity. Compensating for effective stress variations is critical when working with seismic and geotechnical parameters in the near-surface.  The very rapid change in effective stress has been  shown to dominate the response of both shear-wave velocity and cone tip resistance. The effect of the effective stress must first be removed before variations can be attributed primarily to changes in lithology.  The separation of lithologies is critical i f the data are to be used for  hydrogeologic applications because the changes in lithology (particularly between sand and clay) are the most important factor influencing the subsurface hydraulic conductivities. M y investigations in using seismic reflection data for hydrogeologic applications will hopefully serve as the starting point for additional research.  I have shown how a two  dimensional profile of hydraulic conductivity variations can be obtained from surface reflection seismic data.  To determine whether the results make a difference for the development of  hydrogeologic models a test should be performed where model results are computed both with and without the use of the seismic data. The inversions I have completed in the thesis have only utilized SH-wave seismic reflection data. The same procedure also applies to P-wave reflection data and could be a valuable tool where excellent quality P-wave data can be acquired (e.g. the Brookswood aquifer). Finally, the relationships I have developed to convert shear-wave velocity to grain size and hydraulic conductivity are based on a relatively small amount of data. By integrating additional cone penetrometer data such as the friction ratio, pore pressure and resistivity data (from a greater number of cone holes), and additional hydrogeologic data such as direct measurements of hydraulic conductivity and more finely sampled grain-size data, a better constrained series of site-specific relationships can be developed. Well-constrained relationships will result in the development of more accurate groundwater models.  116  REFERENCES Armstrong, J.E., 1984, Environmental and engineering application of the surficial geology of the Fraser lowland, British Columbia: Geological Survey of Canada, Paper 83-23. Berkhout, A . J . , 1984, Seismic resolution: resolving power of acoustical echo techniques: Handbook of Geophysical Exploration, Vol. 12, Geophysical Press. Best, M . E . and Todd, B.J., 2000, Electromagnetic mapping of groundwater aquifers in the Fraser lowland; in Mapping, Geophysics, and Groundwater Modelling in Aquifer Delineation, Fraser Lowland and Delta, British Columbia, ed. B.D. Ricketts: Geological Survey of Canada, Bulletin 552, 27-47. Best, M . E . , Todd, B J . , and O'Leary, D., 1995, Groundwater mapping using time-domain electromagnetics: examples from the Fraser Valley, B . C . ; in Current Research 1995-A: Geological Survey of Canada, 19-27. Bradford, J.H., Sawyer, D.S., Zelt, C.A., and Oldow, J.S., 1998, Imaging a shallow aquifer in temperate glacial sediments using seismic reflection profiling with D M O processing: Geophysics, 63, 1248-1256. Bilker, F., Green, A . G . , and Horstmeyer, H . , 1998, Shallow 3-D seismic reflection surveying: Data acquisition and preliminary processing strategies: Geophysics, 63, 1434-1450. Campanella, R.G., Davies, M . , Boyd, T., Everhard, J., Roy, D., Tomlinson, S., Jackson, S., Schrempp, H . , and Ricketts, B.D., 1994, In-situ testing for the characterization of aquifers: demonstration project: Geological Survey of Canada, Open File 2940. Campanella, R.G., Robertson, P.K. and Gillespie, D., 1983, Cone penetration testing in deltaic soils: Can. Geo tech. Journal, 20, 23-35. Cardimona, S.J., Clement, W.P., and Kadinsky-Cade, K., 1998, Seismic reflection and groundpenetrating radar imaging of a shallow aquifer: Geophysics, 63, 1310-1317. Carr, B.J., Hajnal, Z., and Prugger, A . , 1998, Shear-wave studies in glacial till: Geophysics, 63, 1273-1284. Clague, J.J., 1994, Quaternary stratigraphy and history of south-coastal British Columbia; in Geology and Geological Hazards of the Vancouver Region, Southwestern British Columbia, ed. J.W.H. Monger: Geological Survey of Canada, Bulletin 481, 181-192. Clague, J.L., Luternauer, J.L, and Hebda, R.J., 1983, Sedimentary environments and postglacial history of the Fraser Delta and lower Fraser Valley, British Columbia: Can. J. Earth Sci., 20, 1314-1326. Copty, N . , Rubin, Y . , and Mavko, G., 1993, Geophysical-hydrological identification of field permeabilities through Bayesian updating: Water Resources Research, 29, 2813-2825. 117  Dix, C.H., 1955, Seismic velocities from surface measurements: Geophysics, 20, 68-86. Dobecki, T.L., 1993, High resolution in saturated sediments - a case for shear wave reflection: Expanded Abstracts, SAGEEP 1993 Conference Proceedings, 319-333. Freeze, R.A. and Cherry, J.A., 1979, Groundwater: Prentice-Hall Inc. Freire, S.L.M., and Ulrych, T.J., 1988, Application of singular value decomposition to vertical seismic profiling: Geophysics, 53, 778-785. Gal' perin, E.I., 1974, Vertical seismic profiling: Soc. Expl. Geophys. Gouveia, W.P. and Scales, J.A., 1998, Bayesian seismic waveform inversion: parameter estimation and uncertainty analysis: J. of Geoph. Res., 1 0 3 ( B 2 ) , 2759-2779. Hardage, B.A., 1985, Vertical seismic profiling Part A : Principles, second enlarged edition: Geophysical Press. Hardin, B.O., and Drnevich, V.P., 1972, Shear modulus and damping in soils, design equations and curves: J. of the Soil Mech. and Found. Div., A S C E , 98, 667-692. Hardin, B.O. and Richart, F.E., 1963, Elastic wave velocities in granular soils: J. of the Soil Mech. And Found. Div., A S C E , 89, 33-65. Harris, J.B., Hunter, J.A., Pullan, S.E., Burns, R.A., and Good, R.L., 1996, Shallow shear wave seismic reflection profiling in the Fraser River delta, British Columbia: 66 Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 865-868. th  Hasbrouck, W.P., 1991, Four shallow-depth, shear-wave feasibility studies: Geophysics, 56, 1875-1885. Hofinann, B.A., 1997, In situ ground freezing to obtain undisturbed samples of loose sand for liquefaction assessment: Ph.D. Thesis, Univ. of Alberta. Hunter, J.A.M., Burns, R.A., Good, R.L., and Pelletier, C.F., 1998a, A compilation of shear wave velocities and borehole geophysics logs in unconsolidated sediments of the Fraser River Delta: Geological Survey of Canada, Open File D3622. Hunter, J.A., Pullan, S.E., Burns, R.A., Good, R.L, Harris, J.B., Pugin, A . , Skvortsov, A . , and Goriainov, N . N . , 1998b, Downhole seismic logging for high-resolution reflection surveying in unconsolidated sediments: Geophysics, 63, 1371-1384. Hyndman, D.W. and Gorelick, S.M., 1996, Estimating lithologic and transport properties in three dimensions using seismic and tracer data: the Kesterson aquifer: Water Resources Research, 32, 2659-2670. Jarvis, K . D . and Knight, R.J., 2000, Near-surface VSP surveys using the seismic cone penetrometer: Geophysics, 65, 1048-1056.  118  Jefferson, R.D. and Steeples, D.W., 1995, Effects of short-term variations in near-surface moisture content on shallow seismic data: 65 Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 419-421. th  Lee, M . M . , and Campbell, K.W., 1985, Relationship between shear wave velocity and depth of overburden: Proceedings of Measurements and Use of Shear Wave Velocities, A S C E Fall Convention, 63-76. Lines, L.R. and Treitel, S., 1984, A review of least square inversion and its application to geophysical problems: Geophysical Prospecting, 32, 159-186. Lunne, T., Robertson, P.K. and Powell, J.J.M., 1997, Cone penetration testing in geotechnical practice: Blackie Academic and Professional. Makepeace, A.J. and Ricketts, B.D., 2000, Aquifer mapping and database management using a geographic information system, Fraser Lowland; in Mapping, Geophysics, and Groundwater Modelling in Aquifer Delineation, Fraser Lowland and Delta, British Columbia, ed. B.D. Ricketts: Geological Survey of Canada, Bulletin 552, 17-26. Marion, D., Nur, A., Y i n , FL, and Han, D., 1992, Compressional velocity and porosity in sandclay mixtures, Geophysics, 57, 554-563. McKenna, S.A., and Poeter, E.P., 1995, Field example of data fusion in site characterization: Water Resources Research, 31, 3229-3240. Miller, R.D., Anderson, N.L., Feldman, H.R., and Franseen, E.K., 1995, Vertical resolution of a seismic survey in stratigraphic sequences less than 100 m deep in southeastern Kansas: Geophysics, 60, 423-430. . Milligan, P.A., Rector HI, J.W., and Bainer, R.W., 1997, Hydrophone V S P imaging at a shallow site: Geophysics, 62, 842-852. Monahan, P.A., Luternauer, J.L., and Barrie, J.V., 1993, A delta plain sheet sand in the Fraser River delta, British Columbia, Canada: Quaternary International, 20, 27-38. Neilson-Welch, L . A . , 1999, Saline water intrusion from the Fraser River estuary: a hydrogeological investigation using field chemical data and a density-dependent groundwater flow model: M.Sc. Thesis, Univ. of British Columbia. Olsen, R.S., 1994, Normalization and prediction of geotechnical properties using the cone penetrometer test (CPT): Ph.D. Thesis, Univ. of Calif., Berkeley. Pullan, S.E., Good, R.L., Jarvis, K . , Roberts, M . C . , and Vanderburgh, S., 2000, Application of shallow seismic-reflection techniques to subsurface structural mapping, Fraser lowland; in Mapping, Geophysics, and Groundwater Modelling in Aquifer Delineation, Fraser Lowland and Delta, British Columbia, ed. B.D. Ricketts: Geological Survey of Canada, Bulletin 552, 49-74.  119  Pullan, S.E., Good, R.L., and Ricketts, B.D., 1995, Preliminary results from a shallow seismic reflection survey, Lower Fraser Valley Hydrogeology Project, British Columbia; in Current Research 1995-A: Geological Survey of Canada, 11-18. Pullan, S.E. and MacAulay, H.A., A n in-hole shotgun source for engineering seismic surveys: Geophysics, 52, 985-996. Pullan, S.E., Pugin, A . , Dyke, L.D., Hunter, J.A., Pilon, J.A., Todd, B.J., Allen, V.S., and Barnett, P J . , 1994, Shallow geophysics in a hydrogeological investigation of the Oak Ridges moraine, Ontario: Expanded Abstracts, SAGEEP 1994 Conference Proceedings, 143-161. Rea, J . M . and Knight, R.J., 2000, Characterization of the Brookswood aquifer using groundpenetrating radar; in Mapping, Geophysics, and Groundwater Modelling in Aquifer Delineation, Fraser Lowland and Delta, British Columbia, ed. B.D. Ricketts: Geological Survey of Canada, Bulletin 552, 75-93. Rea, J.M., Knight, R.J., and Ricketts, B.D., 1994, Ground-penetrating radar survey of the Brookswood aquifer, Fraser Valley, British Columbia; in Current Research 1994-A: Geological Survey of Canada, 211-216. Ricketts, B.D., 1995, Progress report and field activities of the Fraser Valley Hydrogeology Project, British Columbia; in Current Research 1995-A, Geological Survey of Canada, 1-5. Ricketts, B.D., 2000, Overview of the Fraser lowland hydrogeology project; in Mapping, Geophysics, and Groundwater Modelling in Aquifer Delineation, Fraser Lowland and Delta, British Columbia, ed. B.D. Ricketts: Geological Survey of Canada, Bulletin 552, 5-15. Ricketts, B.D. and Dunn, D., 1995, The groundwater database, Fraser Valley, British Columbia; in Current Research 1995-A; Geological Survey of Canada, 7-10. Ricketts, B.D. and Makepeace, A J . , 2000, Aquifer delineation, Fraser lowland and delta, British Columbia; mapping, geophysics and groundwater modelling: Geological Survey of Canada, Open File D3828. Roberts, M . C . , Pullan, S.E., and Hunter, J.A., 1992, Applications of land-based high resolution seismic reflection analysis to Quaternary and geomorphic research: Quaternary Science Reviews, 11, 557-568. Roberts, M . C . , Vanderburgh, S., and Jol, H . , 2000, Radar facies and geomorphology of the seepage face of the Brookswood aquifer, Fraser lowland; in Mapping, Geophysics, and Groundwater Modelling in Aquifer Delineation, Fraser Lowland and Delta, British Columbia, ed. B.D. Ricketts: Geological Survey of Canada, Bulletin 552, 95-102. Robertson, P.K. and Campanella, R.G., 1983, Interpretation of cone penetration tests. Part I: Sand: Can. Geotech. Journal, 20, 718-733. Robertson, P.K., Campanella, R.G., Gillespie, D. and Rice, A., 1986, Seismic CPT to measure in situ shear wave velocity: Journal of Geotech. Eng., 112, 791-803.  120  Robertson, P.K., Woeller, D.J., and Finn, W.D.L., 1992, Seismic cone penetration test for evauluating liquefaction potential under cyclic loading: Can. Geotech. Journal, 29, 686-695. Sacchi, M . D . and Ulrych, T.J., 1996, Bayesian priors for sparse inversion: 23rd Ann. Mtg., Canadian Soc. Expl. Geophys., Abstracts, 66-67. Skvortsov, A . G . , Hunter, J.A., Gorianinov, N . N . , Burns, R.A., Tsarov, A . M . , and Pullan, S.E., 1992, High-resolution shear-wave reflection technique for permafrost engineering applications: new results from Siberia: 62 Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 382-384. nd  Steeples, D.W. and Miller, R.D., 1990, Seismic reflection methods applied to engineering, environmental and groundwater problems, in Ward, S.H. Ed, Geotechnical and Environmental Geophysics, Vol. 1: Soc. of Expl. Geophys., 1-29. Tarantola, A., 1984, Linearized inversion of seismic reflection data: Geophysics Prospecting, 51, 332-346. Tokimatsu, K . , Kuywayama, S., and Tamura, S., 1991, Liquefaction potential evauluation based on Rayleigh wave investigation and its comparison with field behaviour: Proceedings, Second Int. Conf. on Recent Advances in Geotech. Earthquake Eng. and Soil Dyn., V o l 1, 357-364. Wood, E., 1996, Comparison of hydraulic conductivity test methods for the Kidd2 site in Richmond, B.C.: B.Sc. Thesis, Univ. of British Columbia. Woodsworth, G.J. and Ricketts, B.D., 1994, A digital database for groundwater data, Fraser Valley, British Columbia; in Current Research 1994-A: Geological Survey of Canada, 207210. Wuenschel, P.C., 1976, The vertical array in reflection seismology - some experimental studies: Geophysics, 41, 219-232. Yilmaz, O., 1987, Seismic Data Processing, Soc. of Expl. Geophys.  121  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0053191/manifest

Comment

Related Items