UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Ground penetrating radar applications in hydrology Rea, Jane Mary Anastasia 1996-12-31

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

Item Metadata

Download

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

Full Text

GROUND PENETRATING RADAR APPLICATIONS IN H Y D R O L O G Y  by Jane Mary Anastasia Rea  B.Sc, Queen's University, 1987 M.Sc, Queen's University, 1992  A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E O F D O C T O R O F PHILOSOPHY in The Faculty of Graduate Studies Department of Earth and Ocean Sciences  We accept this thesis as conforming to the required standard  T H E UNIVERSITY O F BRITISH COLUMBIA DECEMBER 1996 ©Jane  Mary A n a s t a s i a Rea, 1996  In  presenting  degree at the  this  thesis  in  University of  partial  fulfilment  of  of  department  requirements  British Columbia, I agree that the  freely available for reference and study. I further copying  the  or  understood  or  her  representatives.  it  agree that permission for extensive granted  his  an advanced  Library shall make  this thesis for scholarly purposes may be by  for  It  is  by the that  head of copying  my or  publication of this thesis for financial gain shall not be allowed without my written permission.  Department The University of British Columbia Vancouver, Canada ,  DE-6 (2/88)  ABSTRACT  The goal of this thesis is to develop techniques for using ground penetrating radar (GPR) to characterize aquifers for groundwater modelling. GPR is a shallow geophysical technique which can be used to image up to 30 m into the subsurface. This technique is sensitive to many of the same parameters that affect the hydraulic properties of geologic material. In addition it is a rapid, relatively inexpensive means of obtaining high resolution images of the subsurface. For these reasons GPR is well suited for providing the information needed to constrain ground water models. I have developed three different techniques for analysing and interpreting GPR data for ground water modelling. The first is a method of obtaining the correlation structure of the subsurface through geostatistical analysis of GPR data. The second technique involves analysing the attenuation of radar datatoestimate the electrical conductivity of the geologic material being imaged. The final technique is a method of dividing the radar data into separate units, called radar architectural elements, which can be used as building blocks for a ground water model. This technique has been tested by applying it to a shallow unconfiried aquifer in south-western British Columbia. The research presented in this thesis leads to a new integrated approach to GPR interpretation in which the same data set can be used to provide different types of information to constrain ground water models.  Table of Contents Abstract  1  1  Table of Contents  1  1  List of Figures  1  iv  List of Tables  vii  Acknowledgements  viii  Chapter 1: Introduction  1  Chapter 2: Geostatistical analysis of ground penetrating radar data: a means of describing spatial variation in the subsurface Introduction Geostatistical analysis of data Description of the gravel pit study Radar length scales in architectural elements Variograms from 3D elements Conclusions  3 3 4 5 14 19 21  Chapter 3: Obtaining Estimates of Electrical Conductivity from the Attenuation of Ground Penetrating Radar Data Introduction Analysis of radar attenuation Field experiment: Comparison of GPR and DC resistivity results GPR results over two layers Conclusion  23 23 24 29 39 41  Chapter 4: Determination of Hydraulic Boundaries using Ground Penetrating Radar Introduction The interpretation of GPR data Study to investigate the top of the saturated zone Cliff face study The Brookswood aquifer study Discussion Conclusion  43 43 44 48 57 64 76 78  Chapter 5: Summary and Conclusions  79  References  82  Appendix  86  iii  List of Figures Figure 2.1a: Two common semivariogram models are used in this study: the spherical model and the exponential model. Figure 2.1b: A semivariogram showing the hole effect. The hole effect occurs when there is a natural cyclicity in the data set. Figure 2.2: Photograph of the studied cliff face in the Kirkpatrick sand and gravel pit. Two facies are visible. The light areas are composed of medium to coarse sand; the dark areas are fine sand and silt. Figure 2.3: Binary image of the cliff face photograph shown in Figure 2.2. Black represents the fine facies and white represents the coarse facies. Figure 2.4: 2D semivariogram map of the binary data set shown in Figure 3. The blue colors represent low variogram values; the red colors represent high values. The direction of maximum correlation corresponds to the major axis of the blue band or ellipse. Figure 2.5: Semivariogram of the binary image of the cliff face photograph calculated in the direction of maximum correlation of 4 degrees. The semivariogram is modelled using a nested spherical model with two ranges: 2.1m and 0.5 m. Figure 2.6: Radar section collected along the top of the cliff with 100 MHz antennas. Figure 2.7: 2D semivariogram map of the radar data presented in Figures 2.6. The blue colours represent low variogram values; the red colours represent high values. The direction of maximum correlation corresponds to the major axis of the blue band or ellipse. Figure 2.8: Semivariogram of the 100 MHz radar data calculated in the direction of maximum correlation, 4.0 degrees. The variogram is modelled using a nested spherical model with ranges of 2.0 m and 0.8 m. Figure 2.9: Semivariogram of the cliff face photograph calculated in the vertical direction. The oscillation in the variogram is due to the hole effect. Figure 2.10: Semivariogram of the cliff face radar data calculated in the vertical direction. Figure 2.11: Radar survey showing the dipping radar element, dipl, which is associated with dipping foreset beds. Figure 2.12: Semivariograms of the dipping radar element: a) dipl is modelled with a spherical model with a range of 1.4 m, b) dip2 is modelled using an exponential model with a range of 2.1 m, and c) dip3 is modelled using an exponential model with a range of 1.55 m. Figure 2.13: Radar survey showing the short sub-horizontal element, shortl, which is associated with topset sedimentary elements with short correlation lengths Figure 2.14: Semivariograms of the short sub-horizontal radar elements: a) shortl is modelled with a spherical model with a range of 1.3 m, b) short2 is modelled using a spherical model with a range of 1.0 m, and c) short3 is modelled using an exponential model with range of 0.7 m. Figure 2.15: Radar survey showing the long sub-horizontal element, longl, which is associated with topset sedimentary features with long correlation lengths. iv  Figure 2.16: Semivariograms of the long sub-horizontal radar elements: a) longl is modelled using an exponential model with a range of 2.4 m, b) long2 is modelled using a nested spherical model with ranges of 4.8 m and 1.5 m, and c) long3 is modelled using an exponential model with a range of 4.4 m. Figure 2.17: 3D radar survey collected with 100 MHz antennas over dipping foreset beds. Figure 2.18: Semivariogram calculated in the down dip direction of the 3D radar survey. The variogram is modelled using a nested spherical model with 2 ranges: 4.0 m and 1.2 m. Figure 3.1: Block diagram of radiated and returned power for a radar system, (after Annan and Davis, 1977). Figure 3.2: GPR data collected across a boundary between sand and gravel and a clay-rich till with a) 100 MHz antennas and b) 200 MHz antennas. Figure 3.3: Instantaneous amplitude plots of two traces from the pinchout site. One trace was collected over sand, the other over clay rich till. Figure 3.4: A plot of the function F for two traces from the pinchout site, one over sand and one over clay till. The data were collected using a) 100 MHz antennas and 2) 200 MHz antennas. Figure 3.5: Results of an inversion of DC resistivity data collected over the pinchout. The resistive regions are shown in blue, the conductive regions in red. Conductivities are given as log(S/m). Figure 3.6: Radar data collected over a 2-layer earth. The sediments in the upper layer are composed of sand, silty sand and gravel. The lower layer is a blue-gray soft plastic clay. The clay layer dips at an angle of about 25 degrees and outcrops at one end of the survey line. Figure 3.7: A plot of the function F for a trace from the radar data shown in Figure 6. There are two distinct slopes to the function F: one corresponds to the sand and one to the clay. Figure 4.1: A schematic of the PulseEKKO IV GPR unit. Figure 4.2: GPR data from Site 1. Figure 4.3: GPR data from Site 2. Figure 4.4: a) GPR data from Site 1 with SEC gain applied, b) GPR datafromSite 2 with SEC gain applied. Figure 4.5: a) GPR datafromSite 1 with no gain applied, b) GPR data from Site 2 with no gain applied. Figure 4.6: a) Instantaneous phase plot of GPR datafromSite 2. b) Instantaneousfrequencyplot of GPR data from Site 2. Figure 4.7: a) Instantaneous phase plot of GPR datafromSite 1. b) Instantaneousfrequencyplot of GPR data from Site 1. Figure 4.8: a) SEC gained GPR data collected at site 2. b) SEC gained datafromsite 2 collected 2 weeks later, c) subtraction of data in Figures 4.8a and 4.8b. Figure 4.9: a) Photograph of Cliff 1. b) A GPR survey over Cliff 1. c) Interpretive diagram of Cliff 1. v,  Figure 4.10: a) Photograph of Cliff 2. b) GPR survey over Cliff 2. c) Interpretive diagram of Cliff 2. Figure 4.11: a) Photograph of Cliff 3. Horizontal scale is 20 m; vertical scale is 8 m. b) A 100 MHz GPR survey over Cliff 3. Station spacing is 0.2 m. c) A 200 MHz GPR survey over Cliff 3. Station spacing is 0.2 m. Figure 4.12a: The Fraser Lowland in south-western British Columbia showing the location of the Brookswood aquifer. Figure 4.12b: Hydrostratigraphic units of the Brookswood aquifer, (from Halstead, 1986) Figure 4.13: Examples of radar elements from the Brookswood data set: a) Dipping radar element b) channel radar element with interpretive diagram c) horizontal radar element d) horizontal radar element e) attenuation element with interpretive diagram. Figure 4.14: Map of the Brookswood aquifer showing the location of radar sections Figure 4.15: a) GPR survey of Line 24208 and b) interpretation of Line 24208 showing the radar AEs. The white region represents a horizontal element, the light grey is an attenuation element and the dark grey is a channel element Figure 4.16 a) GPR data from line 20204 and b) interpretive diagram showing radar AEs. The light grey indicates an attenuation element, the dark grey indicates a channel element and the white indicates a horizontal element Figure 4.17: SEC gain of a section of Line 20204 showing the reflection from the TSZ. Figure 4.18: a) GPR survey in Stokes pit. b) Interpretation of GPR line from Stokes pit showing radar AEs. Figure 4.19: a) 3D GPR survey collected in the stokes pit and b) interpretive diagram. Figure 4.20: a) and b) GPR surveys collected along 192 Avenue. nd  vi  List of Tables  Table 3.1: Average attenuation constants. Table 3.2: Conductivity results from GPR attenuation analysis (G=0) and DC resistivity survey. Table 3.3: Range in G and a values. Table 4.1: CMP survey results  vii  ACKNOWLEDGEMENTS I am greatly indebted to my supervisor Rosemary Knight who provided excellent advice, support and enthusiasm throughout the past four years. I have been very lucky to work with the Rock Physics Group who created an enthusiastic and pleasant working environment. Dave Butler and Paulette Tercier in particular provided valuable advice and support both technical and emotional. I would like to thank all the friends who helped me out and made the whole experience much more enjoyable. Colin and Shawn in particular provided much needed encouragement, coffee breaks and other life support. To my family I owe a great debt for being the best, and most decidedly biased, support network I could wish for throughout my university career. Most importantly I thank Greg for his patience and support (even from half way around the world) which kept me mostly sane and relatively stable. Finally I owe a debt to the Snapper who provided the motivation and drive needed to finish this thesis. This research was conducted over a period of four years with funding in the first two years from various sources: Geological Survey of Canada, U.S. Environmental Protection Agency, Environment Canada's Environmental Innovation Progra, with funding from the EIP Program, Canadian Space Agency and Energy, Mines and Resources. The last two years of this study were funded by the U.S. Air Force Office of Scientific Research under grant number F49620-95-1-0166.  viii  CHAPTER 1: Introduction The distribution of hydraulic properties in the subsurface is one of the primary controls on groundwater flow. The key to groundwater modelling is therefore an accurate definition and simulation of this subsurface heterogeneity.  Recent studies have attempted to simulate the heterogeneity of the  subsurface through geostatistical modelling (Sudicky, 1986). It has been shown however, that these models do not account for heterogeneity involving discrete geologic boundaries (Webb and Anderson, 1996), a critical issue in that highly conductive geologic units are a major control on groundwater flow. An alternative approach has been proposed which uses discrete structures, based on geologic units, as the main building blocks of a hydraulic model (Anderson, 1989, Webb and Anderson, 1990; Webb and Anderson, 1996; Fogg, 1996; Scheibe and Freyberg, 1990). Average values for hydraulic conductivity are assigned to each discrete unit and the small scale heterogeneities are described through geostatistics. Initial results show a much improved model of groundwater flow (Webb and Anderson, 1996). A critical aspect of this approach is obtaining the initial data with which to constrain these models. Two main types of data are needed: firstly, the geometry of the main hydrogeologic units which are used as the basic building blocks for the model; secondly, a measure of the internal properties of the units, including both the overall hydrogeologic properties and the heterogeneity of these properties within the unit. The two types of data are closely related as determining internal properties can assist in defining the location of important hydrogeologic boundaries. At present these data are most often obtained from outcrop studies, driller's logs and traditional hydraulic techniques such as pump and slug tests. These data are expensive to obtain however, and therefore are often sparsely collected. An additional issue is the sample size of these techniques, which is often so small with respect to the total aquifer volume, that the true heterogeneity of the aquifer cannot be determined. One solution to the problem of data acquisition for groundwater modelling is through the use of ground penetrating radar (GPR). GPR was first used to investigate ice thickness in glaciers. Since then improved technology has expanded the applications of GPR to include detection of buried objects and high resolution imaging of the subsurface for geotechnical and groundwater applications. In a GPR survey, electromagnetic energy in thefrequencyrange between 1 and 1000 MHz is transmitted into the ground.  The behaviour of the electromagnetic waves is affected by the electrical properties of the  subsurface- the electrical conductivity and the dielectric constant.  The electrical conductivity of the  subsurface attenuates the GPR signal so that penetration can be very poor through an electrically conductive medium. Changes in the dielectric properties of the subsurface cause reflections of energy 1  which are detected at a receiving antenna on the surface. The result of a radar survey, the radar image, is a map of reflectors marking interfaces across which there are changes in dielectric properties. The objective of this thesis was to develop specific methods for extracting information from GPR data for use in groundwater modelling. We approach this problem by focusing on the two types of data needed in groundwater modelling: the large scale structure and the internal properties.  The techniques for obtaining these data  are presented in the form of stand alone papers, with each chapter in the thesis representing a paper. The first two papers, presented in Chapter 2 and Chapter 3, investigate techniques for extracting subsurface properties from GPR data. Chapter 2 introduces a method for obtaining the correlation structure of the subsurface from the geostatistical analysis of radar data.  In Chapter 3 a method is developed to  determine the electrical conductivity of subsurface sediments through the analysis of GPR decay rates. As electrical conductivity has been linked in numerous studies to hydraulic conductivity, it provides important constraints on the hydraulic properties of the subsurface. In Chapter 4 the use of GPR to determine large scale hydrogeologic units is investigated. In this chapter a framework for using GPR as a reconnaissance technique to map large scale hydrogeologic units is developed and tested. One of the interesting aspects of this study is that the data were not collected over an ideal site, with optimum coverage and radar penetration. It therefore constitutes a real world problem, in which the data are not ideal. In Chapter 5 the conclusions of each paper are summarised and some recommendations for future research are presented. A location map identifying the GPR surveys and photographs included in this thesis is presented in the Appendix.  2  CHAPTER 2: Geostatistical Analysis of Ground Penetrating Radar Data: A Means of Describing Spatial Variation in the Subsurface INTRODUCTION The modelling of groundwater flow and contaminant transport requires, as a starting point, a realistic description of the hydraulic properties of the subsurface. While it is often possible through direct sampling to obtain estimates of hydraulic properties at specific locations, the challenge is to adequately account for the spatial variation that is present in heterogeneous geologic systems. One way of characterising this spatial variation is in terms of the continuity of units with the same, or similar, hydraulic properties. In modelling of flow and transport, it is generally the continuity of units with high and/or low hydraulic conductivity that is of most interest. A geostatistical approach can be used to determine the spatial correlation of the hydraulic parameter, such as hydraulic conductivity, and this taken as a basis for modelling of the heterogeneous system. In field studies, the data which are analysed to determine the spatial variation or correlation structure of hydraulic properties are usually obtained from cores (e.g. Sudicky, 1986) or well logs (e.g. Ritzi et al., 1994). A recurring concern with the analysis of these forms of data is that only a small volume of the subsurface is actually sampled. To overcome this limitation, we have investigated the use of data obtained from a geophysical technique, ground penetrating radar (GPR), which can provide highresolution, continuous 3D coverage of the subsurface. In this study we use geostatistical techniques to quantify the correlation length of the radar reflectors through semivariogram modelling. We suggest that the correlation structure seen in the radar data is representative of the correlation structure in the hydraulic properties of the subsurface. In this study we build upon the work of Olhoeft et al. (1994) who carried out preliminary studies of geostatistical analysis of GPR data. Olhoeft applied a semivariogram function to GPR data collected over a mixed fluvial-aeolian environment. His results show patterns of variability in correlation length scales, which he suggests can be used to locate the position of changes between statistical regimes in the subsurface. Two types of radar data sets, collected at sites in south-western British Columbia, were used in this study. The first type of data set is one for which we can directly examine the section being imaged with GPR. Data were collected across the top of a cliff face in the Kirkpatrick sand and gravel pit, providing us with an opportunity to assess the relationship between the spatial variation seen in the radar data and the spatial variation seen visually in the cliff face. The other data sets were collected over the  3  Brookswood aquifer, which is composed of glacio-deltaic sediments.  We conducted a geostatistical  analysis of these GPR sections as a means of characterising the correlation structure of the different sedimentary units.  GEOSTATISTICAL ANALYSIS OF DATA Geostatistics is the area of study concerned with understanding and modelling spatial variability. The underlying assumption in using geostatistics to characterise heterogeneity in geologic systems is that properties in the earth are not random, but have some spatial continuity or are correlated over some distance. This means that a parameter measured at one location provides information about parameter values at other locations.  In this study, we have used the standard semivariogram as a way of  summarising information about the spatial relationship between data values. The geostatistical software used to construct and model the semivariograms is GSLIB (Deutsch and Journel, 1992). The experimental semivariogram is described by the following equation (Isaaks and Srivastava, 1989): (1) where h is the lag, or separation vector, between two data points, vj and VJ; and N is the number of data pairs used in each summation. The semivariogram function increases with increasing lag, levelling off at a sill value when the lag reaches the point at which the data are no longer correlated. This lag is referred to as the range and is taken to be an indication of the characteristic length scale of the data set. It is the range of the GPR reflections that we are most interested in obtaining in this study. In most cases a semivariogram is directional, that is, the vector h lies in a specified direction. This direction can be set to any direction of interest, but is generally that of maximum correlation. The direction of maximum correlation can be found in several ways. In this study we have made use of the XGAM program by Chu (1993). XGAM produces a 2D semivariogram map of the data set, which is a compilation of semivariograms calculated in all directions with the centre of the map correspondingtoa lag of zero. The map is colour coded such that blue colours represent low semivariogram values and reds represent high values. The result is an ellipse in the centre of the map. The major axis of the ellipse indicates the most continuous direction. An experimental semivariogram obtained using equation 1 would result in values for y only at 4  certain lags, which are determined by the sampling interval of the data. Semivariogram models are used to provide an analytic description of the experimental semivariogram so that y can be calculated for any lag h. This is needed in order to use the semivariogram for simulation algorithms and to provide a better estimate of the true correlation length of the data set. Two of the more common models, the ones used in our study, are the spherical and exponential models, examples of which are shown in Figure 2.1a.  -spherical model - exponential model  20  30  40  50  Lag (m)  Figure 21 .a:  Two common semivariogram models are used in this study: the spherical model and the exponential model.  These models are described by the following equations (Isaaks and Srivastava, 1989):  Spherical model: h (h y(h) = c- 1.5—.5 a Va  , ifh>a  (2)  ifh<a  =c Exponential model: t \ y(h) = c- [ 1-exp f - —^ 3h  v  aJ  (3)  where c is the sill, and a is the range, which is equivalent to the conelation length of the data set. Some semivariograms require a combination of several ranges to model the semivariogram shape correctly. As suggested by Isaaks and Srivastava (1989), we have modelled our semivariograms with the minimum number of ranges required to produce a reasonable model. In some geologic systems there exists a natural cyclicity, which introduces an oscillation in the semivariogram; the wavelength of the oscillation is approximately equal to that of the cyclicity in the data. This is refened to as the hole effect (Isaaks and Srivastava, 1989) and is modelled with  y(h) = c • {l - [(e~ ) • cos{2nh / X)]}  (4)  h/&  where 5 is a factor to model damping of the hole effect amplitude with increasing lag, and A. is the hole effect wavelength (Journel and Froidevaux, 1982). An example of a model semivariogram corresponding to a pure hole effect is shown in Figure 2.1b. In this example X is 5 m , 8 is 10 m, and c is 0.55 x , where 2  x is the units of data being used in the semivariogram equation.  ; - 0 . 5 -0.4 -0.3 -0.2 0.1  0  |  -I  1  1  1  1  0  10  20  30  40  50  Lag (m)  Figure 2.1b: A semivariogram showing the hole effect. The hole effect occurs when there is a natural cyclicity in the data set. An underlying assumption in the calculation of semivariograms is that the data conform to a requirement called stationarity.  Stationarity means that any subset of the data will have the same  statistical description as any other subset. This is of concern in most geostatistical problems as the earth is not inherently stationary, with different geologic units displaying very different spatial variability. For the purposes of our study, this implies that we must separate, and analyse individually, data obtained from different sedimentary units. In our geostatistical analysis of the GPR data, our objective is to determine the correlation length of the radar reflectors.  The data values we use to construct the experimental semivariogram are the  amplitude values recorded in the radar traces. We must therefore first correct for a violation of the stationarity assumption introduced by the GPR method due to the decay of amplitude down a radar trace. Uncorrected, the data at early times, having much higher amplitudes, would dominate the semivariogram calculation. We first correct the data for any time shifts due to elevation changes or problems with the electronics, then use an automatic gain control with a short window length to make the amplitudes of all the reflectors approximately equal. The high amplitude air and ground arrivals at the start of each trace are removed by simply cropping the data. A further processing step, migration of the data, might be 6  needed to recover the true correlation length in areas with steeply dipping reflectors; migration was found to have little effect on the analysis of the data from our study areas. As a final step before calculating a semivariogram, the radar data must be converted from a time section to a depth section so that the lag is a physical length scale. This is done using the velocity calculated from a common midpoint radar survey conducted at the site under investigation.  DESCRIPTION OF T H E G R A V E L PIT STUDY One of the issues which must be addressed in the proposed use of radar data as a means of characterising the heterogeneity of the subsurface, is the link between the correlation structure obtained from radar data and that of the hydraulic properties. As with most geophysical surveys, radar surveys do not provide the parameter of interest directly. They do not measure hydraulic properties but rather delineate changes in dielectric properties of the subsurface.  There is, however, a fundamental link  between dielectric and hydraulic properties, since both are related to the sedimentary processes and materials which formed the region of the subsurface being investigated.  The basic parameters which  describe a sedimentary unit: grain size, porosity, composition and packing, control to a large extent the dielectric properties of the unit. The same parameters, weighted differently, also control the hydraulic properties of the unit. Our working hypothesis is that the links between sedimentary and hydrogeologic properties, and between sedimentary and dielectric properties, can be exploited to infer an underlying relationship between the correlation structure of the radar reflectors and that of the hydrogeology. As a means of testing this idea we conducted a study of a cliff face in the Kirkpatrick sand and gravel pit, where we were able to image the section with GPR and actually see the sedimentary sequence which we imaged. The cliff consists of a single sedimentary unit: a sand sheet composed of parallel, slightly dipping beds of two different lithologies (Tribe, 1994). One is a coarse to medium grained sand, 2 to 20 cm in thickness. The other is a silt and fine sand layer which can be up to 5 cm thick. In order to quantify the correlation lengths of the lithologic units, we relied on geostatistical analysis of a black and white photograph of the cliff face, shown in Figure 2.2.  The technique of calculating  semivariograms from photographs has been used by other groups, particularly in the field of remote sensing (Woodcock et al, 1988; Atkinson, 1993; Lacaze et al, 1994). In converting this photograph into a format compatible with geostatistical analysis, we have followed the principals of photogeology which state that different lithologies on a photograph can be distinguished by photographic tone, texture patterns and brightness (Drury, 1987). When the photograph was taken, the cliff face was smooth and clear, and 7  the sun was shining directly on it. This minimised the effect of shadows which make it difficult to rely on the greyscale level as a means of separating the image into facies. Image analysis methods were used to  20 m  Figure 2.2: Photograph of the studied cliff face in the Kirkpatrick sand and gravel pit. Two facies are visible. The light areas are composed of medium to coarse sand; the dark areas are fine sand and silt. <  20 m  •  Figure 2.3: Binary image of the cliff face photograph shown in Figure 2.2. Black represents the fine facies and white represents the coarse facies.  8  quantify the distinction between the two cliff face lithologies through thresholding and smc»othing procedures. In the thresholding procedure, any pixel with a greyscale level below the threshold value was changed to white, and any above the threshold was changed to black. The result was a conversion from a grayscale image to a black and white image. The threshold level was adjusted until the black and white regions corresponded closely to the fine and coarse facies visible in the photograph.  The smctothing  procedure corrected for any irregularities at the facies boundaries that were introduced by the thresholding procedure. Thefinalresult, shown in Figure 2.3, is a binary image with a pixel size of 0.27 m x 0.27 m, in which the black represents the fine facies, and the white represents the coarse facies. The binary image of the cliff face photograph was converted to GSLIB format for analysis. In this analysis, we were interested in the maximum correlation direction which corresponds to the dip direction of the bedding planes. There is some variation in this direction throughout the cliff face due to the natural variation found in all sedimentary systems. To quantify our dip estimate, we used a 2D semivariogram map created with the XGAM program. In the semivariogram map shown in Figure 2.4, the blue regions show low semivariogram values and the red regions show high values at, or near, the sill.  North  South  Figure 2.4: 2D semivariogram map of the binary data set shown in Figure 3. The blue colors represent low variogram values; the red colors represent high values. The direction of maximum correlation corresponds to the major axis of the blue band or ellipse. 9  The major axis of the large blue ellipse, or band, indicates the direction of maximum connectivity. In this case, the final result for the direction of maximum correlation is 4 degrees from the horizontal dipping to the north. The experimental semivariogram calculated in the direction of maximum correlation is shown in Figure 2.5. The model of this semivariogram, shown by the solid line in Figure 2.5, is a nested spherical model with two ranges: 2.1 m and 0.5 m. There is excellent agreement between the semivariogram model and the experimental semivariogram, which have a correlation coefficient of 0.99.  The high  degree of correlation is made possible by the large number of data points in the photograph which create a very smooth experimental semivariogram. The two ranges mean that there are two spatial correlation lengths present in this data set. The shorter one, 0.5 m, may be an artefact caused by the variations in dip direction of bedding planes as described above. If the orientation of the lag vector is such that it cuts across these bedding planes, a short correlation length will result. A semivariogram was also calculated using the greyscale image of the cliff directly; it was modelled using a single spherical model with a range of 1.9 m.  0.25  Lag (m)  Figure 2.5: Semivariogram of the binary image of the cliff face photograph calculated in the direction of maximum correlation of 4 degrees. The semivariogram is modelled using a nested spherical model with two ranges: 2.1 m and 0.5 m. Radar data were collected along the top of the cliff with 100 MHz antennas, using a station spacing of 20.0 cm and a 1.0 m offset between antennas. A common midpoint survey was conducted along the cliff and the data analysed to obtain a velocity of 0.1 m/ns. The radar section is 20.0 m long and is displayed as a depth section in Figure 2.6. The radar section shows a series of parallel, slightly dipping reflectors which can be classified as a single sedimentary unit. When we compare the radar data with the photograph we see that the overall structure of the interbedded sands is clearly imaged. The 10  scour surface, however, which can be clearly seen in the photograph running from the top left to the middle right, cannot be seen in the radar image. This is due to the fact that the sediments on either side of the surface are identical, and the surface itself is too small to be resolved in the radar data. D i s t a n c e (m)  Figure 2.6: Radar section collected along the top of the cliff with 100 MHz antennas.  The data were processed as described in the previous section to account for decay in amplitude down the trace and to ensure stationarity. Again, the maximum direction of correlation was found using a 2D semivariogram map produced using the X G A M program. The semivariogram map is shown in Figure 2.7. Note that when analysing radar data the axes of the semivariogram map are not equal because the spatial sampling interval in the horizontal and vertical directions is different.  In the cliff face radar  example, a pixel in the horizontal direction represents a distance of 0.2 m, the station spacing for the radar survey; in the vertical direction it represents a distance of 0.08 m, a value related to the temporal sampling interval of the pulseEKKO radar system and the velocity of the electromagnetic waves in the medium. The apparent angle must therefore be corrected to account for this distortion. The maximum correlation direction, after applying this correction, was determined to be 4.0 degrees, dipping to the north. This is identical to the maximum correlation direction determined for the photograph. The semivariogram calculated in this direction for the 100 MHz data is shown in Figure 2.8. Modelling of the semivariogram resulted in a nested spherical model with ranges of 2.0 m and 0.8 m which compares favourably with the ranges of 2.1 m and 0.5 m calculated from the cliff face photographic image. The experimental semivariogram and the model semivariogram have a correlation coefficient of 0.99. This excellent agreement is due to the large number of data points used to calculate the semivariogram. The erratic variation in the semivariogram shape, seen in most studies, is diminished 11  erratic variation in the semivariogram shape, seen in most studies, is diminished by the averaging effect of the 10 to 50 thousand data points used to calculate the radar semivariogram. North  South maximum correlation direction  4P  Figure 2.7: 2D semivariogram map of the radar data presented in Figures 2.6. The blue colours represent low variogram values; the red colours represent high values. The direction of maximum correlation corresponds to the major axis of the blue band or ellipse.  semivariogram values semivariogram model  Lag(m)  Figure 2.8: Semivariogram of the 100 MHz radar data calculated in the direction of maximum correlation, 4.0 degrees. The variogram is modelled using a nested spherical model with ranges of 2.0 m and 0.8 m  0.18 0.16  semivariogram values - semivariogram model  0.14 0.12 - 0.1 0 08 0 06 0 04 0 02 0 0.6  0.8  Lag (m)  Figure 2.9: Semivariogram of the cliff face photograph calculated in the vertical direction. The oscillation in the variogram is due to the hole effect. 12  between photographic data and GPR data. This suggests that GPR can provide a means of characterising spatial heterogeneity. The other direction which is of interest in characterising hydrogeologic properties is the vertical direction. The experimental semivariogram in the vertical direction constructed using the binary image from the cliff face photograph is shown in Figure 2.9.  Note the oscillation in the  semivariogram due to a hole effect. To model the semivariogram from the cliff face photograph, we used a combination of the hole effect model and a spherical model. For the hole effect model, the wavelength is 28 cm and the damping factor, 8, is 1.0 m. We interpret the wavelength of 28 cm as corresponding to the distance we see between repeating packages of beds; we define a single package as composed of a bed of the coarse sand and a bed of the silt and fine sand. We interpret the range of 0.07 m obtained for the spherical model as an indication of the correlation length within the packages. The semivariogram from the 100 MHz radar data is shown in Figure 2.10.  Once again, the  variogram exhibits a hole effect. In this case, the semivariogram was modelled with the hole effect model alone with a damping factor, 8, of 1.0 and a wavelength of 1.5 m. This wavelength, 1.5 m, is much longer than the wavelength of 28 cm calculatedfromthe cliff face photograph and is most likely affected by a source of repetition in the radar data itself.  0  2  4  a  8  10  Lag (m)  Figure 21 .0:  Semivariogram of the cliff face radar data calculated in the vertical direction.  A radar trace is composed of an energy wavelet convolved with changes in the dielectric properties of the subsurface. In the PulseEKKO radar system, the radar wavelet has two positive peaks, separated by a negative peak. Therefore, in a radar trace one expects two positive peaks from each reflecting surface, thereby producing a cyclicity in the radar trace related to the wavelength of the radar wavelet. For a radar wavelet with a centre frequency of 100MHz, the wavelength is approximately 1.0 m, which is less than the wavelength determined for the hole effect. We conclude that there is some 13  additive relationship between the geologic structure and the radar wavelet length scale which results in the semivariogram shown in Figure 2.10. The complicated behaviour of the semivariogram means that we cannot use radar data to characterise spatial variation in the vertical direction. A solution may be to deconvolve the radar data before calculating the semivariogram. However, as deconvolution of radar data is a topic of current research, this approach has not been attempted.  We therefore limit ourselves to calculating radar  semivariograms only in directions for which the results are unaffected by the radar wavelet. This means we are restricted to directions corresponding to the orientation of reflections in the radar data. This, in general, will be parallel to sedimentary boundaries.  RADAR LENGTH SCALES IN ARCHITECTURAL ELEMENTS Radar data were collected over the Brookswood aquifer in south-western British Columbia. The Brookswood is a shallow, unconfined aquifer composed of glacio-deltaic sediments deposited about 11,000 years BP at the mouth of the Campbell river during the retreat of the Sumas glacier (Ricketts and 2 Jackson, 1994). The sediments referred to as the Brookswood aquifer cover an area of 30 km (Ricketts and Jackson, 1994). Following the framework of Miall (1985), these sediments can be divided into packages, each referred to as an architectural element, that are defined based on geometry, lithofacies and bounding characteristics.  The Brookswood aquifer includes a number of distinct sedimentary  architectural elements characteristic of a glacio-deltaic system.  We expect each element to have a  different correlation structure related to the sediments and processes which formed it.  In the same  manner, radar data can be divided into elements with similar reflection characteristics which are related to the architectural elements being imaged. We have used geostatistical analysis of GPR data to determine the correlation length of several of these elements. Twelve km of radar data were collected over the Brookswood aquifer.  The survey was  concentrated along two lines, one running north-south, and the other east-west, so that radar images of a number of the sedimentary features in the aquifer were obtained. We have focused our analysis on two key sedimentary units imaged in this data set: dipping foreset beds and sub-horizontal topsets. Large sand and gravel dipping foreset beds are the most prominent and visually distinctive architectural element in the Brookswood, and can be seen in numerous gravel pits in the area. Foresets are dipping beds of fining upwards sediments which form when a river or stream empties into a body of 14  standing water. The dips of these beds can be up to 35 degrees. Foreset structures can occur over a large area as the delta progrades into the body of water. The dipping foreset beds are easily identified in the radar data by their characteristic reflection geometry of parallel dipping reflectors, which commonly mark the boundaries between fining upwards sequences in the foreset. Corresponding topsets are also typical of a deltaic environment and can consist of a variety of features including fluvial type deposits and wave or tidally reworked deposits.  Sediments in the  Brookswood are predominantly fluvial facies from braided channels. This includes crossbedded sands, channel features and sand lenses (Ricketts and Jackson, 1994). These sediments are heterogeneous and are expected to occur over a range of length scales, from less than a meter up to 10's of meters. The topsets, being a complex and heterogeneous mix of several architectural elements, have been subdivided for the purposes of this study based on correlation length. All topset reflectors appear in the radar data as sub-horizontal reflectors and are separated into elements exhibiting shorter correlation lengths and those with longer correlation length scales. An important issue which must be considered when calculating geostatistical correlation lengths from 2D data sets is the 3D nature of most geological features. The orientation of the survey with respect to the underlying geologic structure can greatly affect the semivariogram results. In some circumstances, it is important to determine the correlation length in a particular direction, irrespective of the orientation of the underlying geology. Often however, the most useful direction is optimally oriented with respect to the geologic structure. The examples which follow are geostatistical analyses of 2D radar lines. The GPR data were processed using the procedure outlined previously. Radar elements were then isolated from the data set as a whole, and semivariograms calculated for the individual elements using the radar amplitude data. With no prior knowledge of the paleoflow direction at any of the locations, we were unable to select the optimal orientation for our radar lines.  We therefore present these data simply as examples of the  technique of obtaining correlation lengthsfromradar data using geostatistical analysis. We calculate the range for the maximum correlation direction of each particular data set but do not claim that this is the maximum possible for the sedimentary unit being imaged. To address the issue of determining the true maximum correlation length of a sedimentary feature, we present a geostatistical analysis of a 3D data set from which the paleoflow direction can be inferred.  15  Distance (m) 10  20  30  40  CM  ii  Q5 Figure 2.11:  Radar survey showing the dipping radar element, dipl, which is associated with dipping foreset beds.  a)  b) Semivariogram of dip2 Semivariogram of dipl 0.5 0.4  semivariogram values semivariogram model  semivariogram values -semivariogram model  0.3 0.2 0.1 0  Lag (m)  Lag (m)  Semivariogram of dip3  0.9-0 8 --  •  semivariogram values semivariogram model  Lag (m)  Figure 2.12: Semivariograms of the dipping radar element: a) dipl is modelled with a spherical model with a range of 1.4 m, b) dip2 is modelled using an exponential model with a range of 2.1 m, and c) dip3 is modelled using an exponential model with a range of 1.55 m.  Dipping reflectors In Figure 2.11 is shown an example of a package of dipping foreset reflectors, dipl, which is the radar image of the foresets. This and two other examples of the dipping reflectors, dip2 and dip3, were analysed. All of these sections are oriented in an east-west direction. The directions of maximum 16  analysed. All of these sections are oriented in an east-west direction. The directions of maximum correlation for dipl, dip2, and dip3 were found, using the XGAM program, to be 20 degrees to the west, 19 degrees to the west and 17 degrees to the west respectively. The experimental semivariograms were calculated in these directions and are shown with the model results in Figures 2.12a-c. Modelling of the semivariogram for dipl resulted in a spherical model with a range of 1.4 m. For dip2, the semivariogram was modelled by an exponential model with a range of 2.1 m. The semivariogram for dip3 was modelled by an exponential model with a range of 1.55 m. In all cases, there is excellent agreement between model and data with correlation coefficients of 0.99.  Distance (m)  Figure 2.13: Radar survey showing the short sub-horizontal element, shortl, which is associated with topset sedimentary elements with short correlation lengths Short sub-horizontal reflectors The short sub-horizontal reflectors in this glacio-deltaic system are interpreted to be crossbedded sands and gravels from a braided outwash stream. Figure 2.13 is a radar section showing these short subhorizontal reflectors. Examples are presented here of semivariograms calculated from three occurrences of the short sub-horizontal reflectors, referred to as shortl, short2 and short3.  The directions of  maximum correlation for shortl (oriented east-west), short2 (oriented east-west) and short3 (oriented north-south) were found using XGAM to be at dips of 0.5 degrees to the west, 2.0 degrees to the west and 0.0 degrees respectively. The semivariograms for these data, calculated in the maximum correlation directions, are shown in Figures 2.14a-c. The semivariogram for shortl is modelled using an exponential model with a range of 1.3 m. Short2 is modelled using a spherical model with a range of 1.0 m. Short3 is modelled with an exponential model with a range of 0.7 m. Once again, there is a correlation coefficient of greater than 0.99 between experimental and modelled semivariograms for all three examples  17  Semivariogram of shortl  semivariogram values -semivariogram model  0.9 0.6 0.7 0.6 -0.5 0.4 0.3 0.2 0.1 0  1  2  3  Semivariogram of short2 1 0.9 0.8 0.7 0.6 -0.5 0.4 0.3 0.2 0.1 0  semivariogram values . semivariogram model  4  1  2  3  4  5  Lag (m) Lag (m)  c) Semivariogram of short3  semivariogram values • semivariogram model  1  2  3  4  Lag (m)  Figure 2.14: Semivariograms of the short sub-horizontal radar elements: a) shortl is modelled with a spherical model with a range of 1.3 m, b) short2 is modelled using a spherical model with a range of 1.0 m, and c) short3 is modelled using an exponential model with range of 0.7 m. Long sub-horizontal reflectors The long sub-horizontal reflectors in the Brookswood data set are interpreted as sand sheets. An example of a radar section showing these long, sub-horizontal reflectors is given in Figure 2.15. Three examples of this element were analysed: longl, long2 and long3, all oriented east-west. The direction of maximum correlation for longl, long2 and long3 was found using X G A M to be at a dip of 3.0 degrees to the east, 4.0 degrees to the east and 3.0 degrees to the west respectively. The semivariograms calculated in the maximum correlation direction are shown in Figure 2.16a-c. Longl and long3 are modelled using an exponential model with ranges of 2.4 m and 4.4 m respectively, while long2 is modelled using a nested spherical model with two ranges: 1.5 m and 4.8 m. The correlation coefficient between experimental and model semivariogram in all cases is 0.99. As with the cliff face example, the shorter correlation length in long2 may be due to a slight misalignment of the semivariogram direction with the dip direction of some bedding planes.  18  Distance (m) 10  30  20  40  PIT  50 ..  70  60 1  rt*)..*.^.  .. • T H . ' .  _  ....  ^  _  Figure 2.15: Radar survey showing the long sub-horizontal element, longl, which is associated with topset sedimentary features with long correlation lengths. a)  b) Semivariogram for longl  c)  Semivariogram for long2  Semivariogram for long3  Lag (m)  Figure 2.16: Semivariograms of the long sub-horizontal radar elements: a) longl is modelled using an exponential model with a range of 2.4 m, b) long2 is modelled using a nested spherical model with ranges of 4.8 m and 1.5 m, and c) long3 is modelled using an exponential model with a range of 4.4 m. VARIOGRAMS F R O M 3D RADAR ELEMENTS In order to investigate the effect of survey line orientation on the variogram results, we collected a 3D radar set at the same site as the dipping reflectors described above. The data were collected at half meter intervals over a grid 25.0 m by 50.0 m. The result of this survey is shown in Figure 2.17. For this example, we are interested in the length scale of the foreset beds in the down dip, or paleoflow, direction. We choose this direction to provide a comparison with the ranges calculated from the 2D radar data of dipping reflectors, which were also calculated in the down dip direction. 19  50 m  16m  Figure 2.17: 3D radar survey collected with 100 MHz antennas over dipping foreset beds.  GSLIB has algorithms which calculate semivariograms from a 3D data set. Therefore, the spatial length scales can be determined in much the same way as the 2D data. However, the large number of data made the use of these algorithms much too computer intensive. To avoid this problem, we used a 3D imaging package to determine the strike of the dipping beds.  Knowing the strike direction, we can  judiciously extract 2D slices from the 3D data set to determine correlation lengths in critical directions. As an example, we present a 2D data set extracted in a direction perpendicular to strike, that is, in the down dip direction: south-west, or 135 degrees west of north. Using X G A M the maximum correlation direction of 30 degrees was found.  The semivariogram from these data, shown in Figure 2.18, is  modelled using a spherical model with two ranges: 4.0 m and 1.2 m; the nested range of 1.2 m is believed, again, to be an artefact caused by slightly irregular bedding directions. The length scale of 4.0 m is considerably longer than the ranges calculated from the 2D radar image of the same dipping 20  reflectors.  >-0.5  • 1  0  1  1  semivariogram values semivariogram model  1  2  3  4  5  6  7  8  9  10  Lag (m) Figure 2.18: Semivariogram calculated in the down dip direction of the 3D radar survey. The variogram is modelled using a nested spherical model with 2 ranges: 4.0 m and 1.2 m. This example points out an important issue in any geostatistical analysis of a 2D data set and highlights the potential value of using GPR to characterise the subsurface.  The correlation length  calculated from 2D data may be significantly different from the maximum length scale of the architectural element.  With a 3D GPR survey, the orientation of the underlying geology can be  determined and the data analysed along relevant lines such as along strike or perpendicular to strike.  CONCLUSIONS We conclude that the geostatistical analysis of GPR data provides a means of characterising the spatial heterogeneity of the subsurface. The technique that we have developed involves the calculation and modelling of semivariograms using the amplitudes contained in the radar data. A critical assumption in applying this method is that the correlation structure seen in radar data is representative of the correlation structure of the subsurface. The results of the cliff face study suggest that this is a reasonable assumption. The semivariogram obtained from a digital photograph of the face showed excellent agreement with the semivariogram from the radar data, in determining both the maximum correlation direction and the correlation length. In areas of complex geologic structure, we found that this method can be applied by first recognising and separating the various sedimentary units, and then determining the correlation structure of each unit. By separating the different units we were able to ensure that the stationarity requirement was met in construction of the semivariograms. 21  In the glacio-deltaic environment, we obtained  correlation lengths ranging from 0.6 m to 4.8 m, for the various units seen in 2D radar sections. The final example of analysis of a 3D radar data set highlighted a significant advantage in the use of radar data for characterising the structure of the subsurface. By collecting 3D data we were able to identify the direction of interest and determine the correlation length in that direction. In this example we chose to find the paleoflow direction and found a correlation length of 4 m. In this study we have obtained excellent results in the geostatistical analysis of radar data. We reiterate that our key assumption is that the correlation structure that we determine for the radar image is representative of the hydraulic structure of the subsurface. We conclude, based on the examples in this study, that this is a valid assumption and that radar data can be used in this way to quantify the spatial heterogeneity of the subsurface.  22  C H A P T E R 3: Obtaining Estimates of Electrical Conductivity from the Attenuation of Ground Penetrating Radar Data  INTRODUCTION  Characterisation of the subsurface for the purpose of modelling groundwater flow or contaminant transport can be described as involving two steps. The first step is to identify and locate the large scale hydrogeologic boundaries. Such boundaries include the water table, and the interfaces between lithologic units with different hydrogeologic properties. Locating these boundaries provides a geometric model of the hydrogeologic structure of the subsurface at a scale that can range from meters to 10's or 100's of meters.  The second step is to determine hydrogeologic properties within this geometric framework.  These two steps of defining the large scale structure and assigning local properties are closely related as determining hydrogeologic properties can assist in defining the location of important hydrogeologic boundaries. A technique that has been widely used to locate hydrogeologic boundaries is ground penetrating radar (GPR). Velhdis et al. (1990) and Rea and Knight (1995) used GPR data to locate and monitor the position of the water table. Van Overmeeren (1994) and Beres and Haeni (1991) used GPR to map the aquifer/aquitard interface.  Jol and Smith (1992), Baker (1991) and Pratt and Miall (1993) have  conducted a number of studies where lithologic units are identified within GPR profiles. There is also the potential use of GPR data to address the second step in characterisation - that is for determining properties. The GPR response to a volume of the subsurface is governed by the spatial variation in both the dielectric constant and electrical conductivity. In order to obtain information about hydrogeological properties from GPR data, we must first analyse the GPR data to determine the values of these two properties, and then relate these measured physical properties to hydrogeologic properties. This would transform a GPR image of the subsurface into a "map" of hydrogeologic properties. Most of the focus to date has been on the relationships between dielectric constant values, obtained from GPR data, and hydrogeologic properties such as water content, porosity and hydraulic conductivity (Knoll and Knight, 1994; Knoll et al., 1995; Greaves et al., 1996). In general however, these relationships are nonunique, and empirical relationships or additional sources of information are needed to link the dielectric constant to the hydrogeologic property of interest. We suggest that extracting information about electrical conductivity from GPR data would improve our ability to use GPR data to obtain a hydrogeologic model of the subsurface. Quantifying electrical conductivity would provide valuable complementary information in relating dielectric constant 23  values to hydrogeologic properties. In addition, the measurement of electrical conductivity has been shown to be a useful indicator of hydraulic conductivity in groundwater studies (Mazac, 1990; Mazac et al., 1988; Kosinski and Kelly, 1981; Urish, 1981; Kelly, 1977). In these studies links are made between electrical and hydraulic conductivity for different geological environments.  Determining electrical  conductivity could therefore be useful both in estimating local properties and in refining models of the large-scale hydrogeologic structure of the subsurface. The effect of electrical conductivity is seen in GPR data as the attenuation in the amplitude of the electromagnetic wave as it propagates through the earth. Analysis of the attenuation recorded in the GPR trace therefore provides a measure of the electrical conductivity. There has been little attempt to use this feature of the GPR response quantitatively. An exception is a study by Livelybrooks (1994) who used a forward modelling approach to match the attenuation in GPR data in order to estimate electrical conductivity in the wall rocks of a mine. The objective of our study is to develop a procedure for analysing the attenuation recorded in GPR data to estimate the electrical conductivity of the sediments being imaged. To develop and assess this procedure, we conducted two field experiments.  In the first we  collected and analysed GPR data across a boundary between a sand and gravel unit and a clay-rich till. Electrical conductivity estimates were obtained from these data and compared with results from a direct current (DC) resistivity survey conducted at the same site. The second field experiment was conducted at a location where a layer of sand and gravel overlies a layer of clay. This allowed us to examine the case where conductivity changes significantly and abruptly with depth. The electrical conductivities obtained from this analysis were compared with results from a cone penetrometer (CPT) survey conducted at the site. ANALYSIS OF THE ATTENUATION OF GPR DATA In a GPR survey electromagnetic energy is transmitted into the ground, reflected back to the surface at interfaces in the subsurface and recorded at the receiving antenna. A radar trace is a record of the time variation in amplitude of the received signal. There is a variation in amplitude along the trace which is determined by the contrast in the dielectric properties across each reflecting interface. In addition there is an overall decay in amplitude due to loss of power in the electromagnetic wave. This power loss is caused by factors related to the design and efficiency of the radar system, the travel path of the wave, and the presence and geometry of reflecting interfaces. There is also power loss associated with the intrinsic attenuation of the geological materials; this is determined primarily by the electrical  24  conductivity. By analysing the decay rates of the amplitude of the radar signal we can therefore extract information about the electrical conductivity of the subsurface. Figure 3.1, modified from Annan and Davis (1977), summarises all the potential causes of power loss from a radar pulse starting at the transmitting antenna, travelling through the subsurface, being reflected at a boundary and finally being recorded at the receiving antenna. Referring to Figure 3.1, there is the initial power at the source, P , which undergoes a power loss associated with the antenna S  efficiency, <| . The power radiated downwards into the ground is further affected by the antenna gain, Tx  G , which is a measure of the directivity of the antenna. The power reaching a reflector has been further Tx  reduced due to the intrinsic attenuation of the medium through which the wave has travelled. This decay occurs exponentially as e~  2aL  , where a is the attenuation coefficient of the medium and L is depth; for  this specific example shown in Figure 3.1, L is the depth to the reflecting boundary. The attenuation coefficient is commonly given in units of nepers/m, where a "neper" is a unitless parameter. Spherical spreading of the wavefront decreases the power radiating directly downwards by a factor of l/(47tL ). 2  When the electromagnetic wave reaches a boundary, some power is reflected back to the surface. The amount of reflected power is governed by the cross-sectional area of the reflecting interface, x, and by the back-scatter gain, g.  Before reaching the receiver, the power is reduced once again by spherical  spreading losses and the intrinsic attenuation of the medium. At the receiver antenna, the power is affected by the antenna effective area A, the antenna gain G^, and the antenna efficiency c^. These factors can be summarised by the following expression for the received power  (Annan and Chua,  1992): £TX '^RX  G  T x - G A-g-x-e" Rx  PRX  16TI L 2  4  The power P ^ can be given as  where  is the voltage at the receiver antenna is the impedance of the receiver antenna.  The effective antenna area is defined as GR A=V G /4n f 2  X  where v is the propagation velocity f is the frequency of the signal.  25  2  RX  4  P<s  (3.1)  POWER FROM  POWER AT RECEIVER E L E C T F RONICS  SOURCE  Ps  PRX=^ =  RADIATED P O W E R  Rx  'P7  RECEIVED POWER  P =  P. = ^ T x - P s  AG P Rx  7  P O W E R RADIATED DOWNWARD  6  P O W E R DENSITY A T RECEIVER „-2aL  P =G 2  T x  -P,  ^6=^-  P O W E R DENSITY REACHING REFLECTOR  P  T  4itL  5  2  POWER DIRECTED BACK T O RECEIVER  -2aL  P =gP4  4TCI/  5  POWER SCATTERED BY R E F L E C T O R  P =X 4  Where: £ = transmitter antenna efficiency ^.RX = receiver antenna efficiency a = attenuation coefficient of medium X = scattering cross sectional area A = receiver antenna effective area TX  g = backscatter gain of reflector G = transmitter antenna gain G,^ = receiver antenna gain L = depth Tx  Figure 3.1: Block diagram of radiated and returned power for a radar system, (after Annan and Davis, 1977). In this figure P denotes power and P denotes power density. In order to extract information about the electrical conductivity of the geological material through which the electromagnetic wave travels, we need to isolate the loss associated exclusively with the attenuation constant, a. If we consider a plane electromagnetic wave in a homogeneous medium, a can be defined in terms of the electrical conductivity a, the dielectric permittivity 8, the magnetic permittivity u and the frequency of the electromagnetic wave f, through the equation:  26  a:  -(27if)  2 Ll  8  (27xf) L eJl + [ ^ r 2  +  l  N  (3.2)  2  This equation assumes that the electromagnetic wave contains only a single frequency component. Although the transmitted GPR signal is a pulse comprising a range of frequencies, we assume that f is equal to the centre frequency, and use this value in applying Equation 3.2. The dielectric permittivity can be determined using the well-known approximation e=e c /v , where c is the speed of light in a vacuum, 2  2  0  3.00xl0 m/s, s is the permittivity of free space, 8.85xl0" and v is the velocity of the electromagnetic 8  12  0  waves in the medium. The velocity can be obtained from analysis of a common midpoint (CMP) survey collected at the site. Having set up Equation 3.1, we could solve for s directly. However the parameters P , Z^, cj , s  %RX> GTX  a n  Tx  d G ^ are in general difficult to evaluate. We therefore need to re-formulate the problem to  avoid the explicit evaluation of these parameters. This can be accomplished by rewriting Equation 3.1 as: F = ctL + c  (3.3)  where f  F=  1  17.38 1  1.738  -10 /og(  ) - 40 • log L + 10 log + I0log(g) + \0log(x) vf y 2  PSZRX^TX^RXGTXGRX  _l  *V  64TI  3  and all logarithms are to the base 10. The constant c represents all the GPR system parameters for which the value is unknown. These parameters can be considered constant for all surveys. The procedure is to evaluate the function F for all observed data points in the GPR trace. The average attenuation constant can be found by determining the slope of the line L versus F. By taking this approach the constant c need not be evaluated since it does not affect the slope. The data at earliest times in a radar trace are dominated by the air and ground arrivals. The amplitude decay in these data cannot be described by Equation 3.3, which considers all received power to have been reflected. These data will therefore not be used in determining the attenuation constant. At later times, the data in a radar trace are dominated by noise. These data will also be omitted from the analysis. In Equation 3.3, the received voltage, the depth, the velocity and thefrequencyare all determined directly from the radar data or from the survey parameters. The scattering cross-sectional area, and the backscatter gain, must be estimated. The method of obtaining each term is described below. 27  Parameters obtained from the GPR survey In the PulseEKKO IV GPR system used in this study, radar data are stored as integers which correspond to the measured voltages V (t). As we are interested in the total power of the signal at an m  instant in time, we convert the measured voltages to instantaneous amplitudes which measure reflectivity strength  (Yilmaz, r  v  Rx  2  1987).  Instantaneous  amplitude  is  given  by  the  expression  "j 1/2  2  (t)= V(t) + v '(t) m  (Yilmaz, 1987), where V(t), the imaginary part of the trace, is obtained by applying a Hilbert transform to the real or measured part of the trace V (t). m  In our analysis we are interested in the overall attenuation of the amplitude. As it is more convenient to deal with units of depth, the sampling times t are converted to depths using the velocity of the electromagnetic waves. The velocity of the waves is obtained from analysis of common midpoint (CMP) data collected at the same site as the radar section. The depth corresponding to each data point is determined assuming a constant using L=v(t/2), where t is the recorded 2-way travel time. Scattering cross-sectional area,velocity, % The scattering cross-sectional area can be estimated following the procedure used by Annan and Davis (1977) in which the reflecting boundary is classified as either a smooth or rough planar spectral scatterer  using the Rayleigh criterion.  The cross-sectional area is equivalent to TCL for a smooth 2  scatterer, and is equivalent to 7tLv/2f for a rough scatterer (Annan and Davis, 1977). In this analysis we do not consider point scatterers since in this case scattering cross-sectional area is a complicated function requiring specific knowledge about the scatterer (Annan and Davis, 1977). The Rayleigh criterion states that a surface can be considered smooth if the average height of irregularities, h, conforms to the  X inequality h <  , where X is the wavelength of the electromagnetic wave and 0 is the angle of  SsinQ incidence off the surface. In GPR surveys, the angle of incidence is assumed to be 90 degrees due to the small separation between receiving and transmitting antennas. The wavelength is estimated by using the central frequency of the antenna and the velocity of the electromagnetic wave in the geological material.  Backscatter gain, g The backscatter gain is the least well constrained of the parameters in Equation 3.3. It arises due to the partitioning of energy at an interface as some energy is transmitted through the interface and some is reflected back to the surface. In Equation 3.3, the backscatter gain is considered for a single boundary. 28  For our purposes we need to account for loss from a series of boundaries. We require a function G(L) which describes the cumulative loss from backscatter gain with increasing depth. This loss is governed by the amount of energy reflected back to the surface at each boundary, and by the number of boundaries. The ratio of the amplitude of the reflected wave to the amplitude of the incident wave is given by the reflection coefficient R, defined as:  (3.4)  where s*,, e* are the complex dielectric permittivities of the upper and lower layers respectively. We 2  consider only the highfrequencycase in which the real part of the dielectric permittivity is much greater than the imaginary part, and so neglect the imaginary part of the dielectric permittivity. The level of reflected power is proportional to R , and the transmitted power to (1-R) . If there 2  2  are n boundaries per meter and each boundary has a reflection coefficient R, the power loss of a pulse travelling from the surface to a depth L and back is R (l-R) . 2  4nL  Multiple reflections are not  considered in this analysis. The equation for the energy partitioning loss G, at a depth L, is  G(L) = 20log(R) + 40L • n • log(\ - Rj.  (3.5)  To determine R, some estimates must be made for the dielectric contrasts in the subsurface. The real part of the dielectric permittivities can be estimated from velocity information obtainedfromCMP surveys, through the approximation E^CVV ,  (3.6)  2  where c is the speed of light in a vacuum, 3.00 xlO m/s and e is the permittivity of free space, 8.85xl0" 8  12  0  F/m. Alternatively, if more is known about the material being imaged, dielectric permittivities can be estimated using an effective medium theory or simple mixing law such as the complex refractive index model (Wharton, et al., 1980).  FIELD EXPERIMENT: COMPARISON OF GPR AND DC RESISTIVITY RESULTS A field experiment was conducted to compare conductivity estimates obtained from GPR data with those obtained from a direct current (DC) resistivity survey. The survey line chosen for this test crosses a boundary between sand and gravel and a clay-rich till. The data therefore include examples from materials with very different electrical conductivities. Both the GPR survey and the (DC) resistivity  29  survey were collected along the same line. The surveys were conducted from east to west along a 100 m transect.  Analysis of GPR Data The GPR data were collected with both 100 MHz and 200 MHz antennas. The use of the two different sets of antennas allowed us to investigate the possibility of a frequency dependence in the observed attenuation. We used a station spacing of 0.2 m and an antenna separation of 0.8 m. The first step in the attenuation analysis was the determination of the velocities of the electromagnetic waves. To do this we analysed the hyperbolic reflections from a common midpoint (CMP) survey. In the sand and gravel, the velocity was found to be 0.065 m/ns for the 100 MHz data. In the till region penetration was so poor that CMP analysis was not feasible. We therefore assumed that the velocity in the till was the same as that in the sand and gravel, as is recommended by Annan and Cosway (1992).  We used the same velocity for both 100 and 200 MHz data.  All depths were therefore  determined using the velocity of 0.065 m/ns and the times recorded in the radar trace. The 100 MHz and 200 MHz GPR data displayed as depth sections are shown in Figure 3.2a and Figure 3.2b. The data were processed to correct for time zero drift and signal saturation. For display purposes an automatic gain correction (AGC) was applied. In both sections penetration was very poor in the eastern 30 m. Reflections in this region are not seen at depths below about 3 m. The penetration increased to the west reaching a maximum depth of about 15 m. Numerous horizontal to slightly dipping reflections are seen in the western portion of the section. There are two drill holes within 100 m of this line. Well number WW092G.007.2.3.1.030* is located at 2027 204A Street near the western end of the radar survey line and consists of sand and gravel, followed by 1 m of sandy gravel with boulders and 13 m of sand and gravel. Well number WW092G.007.2.3.1.041* located at 20633 20 Ave. is just east of th  the survey line. Data from the well indicate that the top 7 m are brown sandy clay followed by 2 m of till and 12 m of stony clay. The water table is inferred to occur at a depth between 3 m and 4 m, due to the continuous reflector which occurs at this depth in both the 100 and 200 MHz data. This is marked in Figure 3.2a and Figure 3.2b. The difference in attenuation in the data from the sand/gravel and the clay-rich till regions can be  * well numbers assigned by the Groundwater Section of the Ministry of Environment and Parks, Province of British Columbia. 30  <3 <D  ©  w  ©  (ui) qjdarj  31  e c CO N  o o  e  CM  3  -o a;  a  C-J  cn <u t3  6X1  e  US (ui) ludarj clearly seen in the instantaneous amplitude plots of the 100 MHz radar traces, examples of which are shown in Figure 3.3. In this figure, the data affected by the air and ground arrivals, and the data which are dominated by noise are not plotted, as these regions were not used in the attenuation estimate. Both traces have similar maximum amplitudes occurring at about 15 ns, but the amplitudes in the clay trace 32  decay much more rapidly, so that by 75 ns the trace is pure noise.  The sand trace, however, has  amplitude peaks well above the noise level past 150 ns. This obvious difference in attenuation was quantified to determine the conductivities of the different sediments. The next step in the analysis required estimates for the cross-sectional area of the reflecting interface and backscatter gain. For the cross-sectional area, we chose between the rough and smooth approximations based upon the Rayleigh criterion. The critical height calculated for the velocity of 0.065 m/ns is 81 mm for 100 MHz data and 40 mm for 200 MHz data. We assumed that the height of surface irregularities on the interface was approximately equal to the radius of the sediments being imaged. At this site the maximum grain size recorded in the drill records is gravel, which according to the USDA scale has a radius between 1 and 40 mm. We therefore assumed that the smooth approximation was appropriate for this site. 1  50000 i 45000 -f  0  50  100  150  200  250  300  350  Time (ns)  Figure 3.3: Instantaneous amplitude plots of two traces from the pinchout site. One trace was collected over sand, at 100 m along the survey, and the other over clay rich till at 20 m along the survey. In order to estimate the backscatter gain, we needed to assign values to the number of boundaries and to the reflection coefficients of these boundaries Without any information about the sedimentary structure, composition etc., such estimates cannot be obtained. We therefore decided to use a backscatter gain of one. While this is obviously not the case, setting the gain to one will yield the upper bound, on possible values for electrical conductivity. This is because all loss due to energy splitting at interfaces is attributed to conductive loss. Further study and investigation may enable us to obtain better estimates of the backscatter gain. 33  Having accounted for all of the factors in the attenuation equation we evaluated the function F for each of 10 traces in the sand and clay-rich till regions for both the 100 and 200 MHz data sets. These traces are located at 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 82, 84, 86, 88, 90, 92, 94, 96, 98 and 100 m along  .  6  J  Depth (m) Figure 3.4a: A plot of the function F for two traces from the pinchout site, one over sand and one over clay till. The data were collected using 100 MHz antennas.  Depth (m) Figure 3.4b: A plot of the function F for two traces from the pinchout site, one over sand and one over clay till. The data were collected using 200 MHz antennas. the survey. Two examples of F versus Lfromthe 100 MHz data, one for the sand and gravel region and one for the clay-rich till region, are shown in Figure 3.4a. Examplesfromthe 200 MHz data are shown in Figure 3.4b. The attenuation constants are equated to the slope of these plots and are found using a least 34  squares technique. The attenuation factor was determined in this way for each of the 10 traces in both the sand and gravel and clay-rich till regions and average values obtained.  A summary of the average  attenuation constants and standard deviations is given in Table 3.1. Table 3.1: Average attenuation constants. sedimentary unit frequency average a standard deviation (MHz) (nepers/m) (nepers/m) sand/gravel 100 0.14 .02 sand/gravel 200 0.18 .04 clay-rich till 100 1.2 0.2 clay-rich till 200 2.3 0.3 The average attenuation constant for the sand was found to be 0.14 nepers/m with a standard deviation of 0.02 for the 100 MHz data, and 0.18 nepers/m with a standard deviation of 0.04 for the 200 MHz data. For the clay-rich till, the average attenuation constant was found to be 1.2 nepers/m with a standard deviation of 0.2 nepers/m for the 100 MHz data, and 2.3 nepers/m with a standard deviation of 0.3 nepers/m for the 200 MHz data. The conductivity was determined from the attenuation constant using Equation 3.1.  The  dielectric permittivity, determined using Equation 3.4 with a measured velocity of 0.065 m/ns, was found to be 21e for the sand and gravel. As there was no velocity information for the clay-rich till, we used the 0  same dielectric permittivity as that for the sand and gravel. The resulting conductivities for the sand and gravel were 3.5 mS/m for the 100 MHz data and 4.3 mS/m for the 200 MHz data. The corresponding results for the clay-rich till were 27 mS/m and 51 mS/m. An interesting observation from these results is that the 200 MHz data produce higher conductivity estimates than the 100 MHz data, for both the sand and gravel unit and the clay-rich unit. We propose two possible causes for this difference. The first is that electrical conductivity is known to be afrequencydependent property, with conductivity observed to increase with increasing frequency. Thus the conductivities of the material being imaged should differ when measured at 100 and 200 MHz, as is the case in our data. A second possible cause of the variation in conductivity withfrequencymay be the difference in the wavelength of the 100 MHz and 200 MHz signals.  The 100 MHz signal has a  wavelength of approximately 1 m, roughly twice that of the 200 MHz signal and so averages over a larger volume of material than the 200 MHz signal.  We therefore expect different results from the two  frequencies, given that the subsurface is likely to be heterogeneous on the same scale as the wavelengths of the measurement. It is not possible, from the data available in this study, to differentiate between these effects. In fact the observedfrequencydependence is most likely a combination of these two effects.  35  Results from the DC Resistivity Survey A DC resistivity survey was conducted along the same line as the radar surveys using a Phoenix Model 1PT1B transmitter and a Huntec Model M-4 receiver. A dipole-dipole geometry was used, with a five metre spread. There were 10 stations with readings taken at 7 offsets from each transmitter location, giving a total of 70 data points. The data from this survey were inverted using a DC resistivity inversion package, DCINV2D (Oldenburg and Li, 1994) to produce a map of conductivity in the subsurface shown in Figure 3.5. The inversion procedure aims to establish the simplest possible subsurface model which honours the DC resistivity data collected in the field. The range of conductivities obtained from the inversion represents average conductivities for each lm x lm grid block in the subsurface model. The result of the inversion reveals a conductive layer, shown in red and yellow, between about 0 and 40 m along the survey line. This  o  —  o  r o u t f - c n » ^ i C D s D « -  o  o  o  o  o  o  o  Distance (m)  o  o  -3.17  o  - . o R 3  4  -3.63 ™ Figure 3.5: Results of an inversion of DC resistivity data collected over the pinchout. The resistive regions are shown in blue, the conductive regions in red. Conductivities are given as log(S/m). layer, interpreted as the clay-rich till, dips to the west at a steep angle and has a conductivity which ranges from about 9.5 to 47 mS/m. The resistive zone, interpreted as the sand and gravel, is shown in blue and has conductivities which rangefromabout 0.2 to 1.6 mS/m.  36  The DC resistivity results compare well with the estimates obtained from the radar data. All results are summarised in Table 3.2. The sand and gravel conductivities from the radar data are slightly higher than those from the DC resistivity survey. This is an expected result since in our analysis of the GPR data we have calculated the maximum conductivity. The clay-rich till conductivities fall within the range found in the DC resistivity survey. An additional comparison of results can be obtained from time domain electromagnetic (TDEM) soundings conducted along 20 Avenue (Best, et al. 1995). A sounding at the comer of 208 Street and th  th  20 Avenue resulted in values of 25 to 50 mS/m for the top 5 m of sediments. This location is 500 m east th  of the clay-rich till imaged in Figure 3.2a. However, a GPR survey conducted in this area showed the same character of poor reflections as in the eastern part of the GPR survey in Figure 3.2a. We therefore infer that the TDEM sounding measured the same sedimentary unit as the GPR survey shown in Figure 3.2, and the DC resistivity survey shown in Figure 3.5. Another sounding at the corner of 204A Street and 20 Avenue resulted in conductivities of 0.1 to 5.0 mS/m for the top 20 m of sediments. th  This  location is about 100 m west of the end of the GPR line shown in Figure 3.2. We again infer that the sediments measured in this TDEM sounding are from the same unit as those imaged in the western section of the GPR survey shown in Figure 3.2.  These results also compare favourably with those  obtained from the attenuation analysis of GPR data. Table 3.2 : Conductivity results from GPR decay rate analysis (G=0) DC resistivity survey and TDEM surveys. sedimentary unit sand/gravel clay till  GPR a mS/m (100 MHz) 3.5 27  Results a mS/m (200 MHz) 4.3 51  DC resistivity results a mS/m  TDEM results a mS/m  0.2-1.6 9.5 - 47  5 33  Discussion One of the issues which must be addressed is the effect of the estimates of the scattering crosssectional area and the backscatter gain on the conductivity results. In order to test the sensitivity of the conductivity calculations to these parameters, we varied both the type of reflector (rough versus smooth) and the estimate of backscatter gain in our calculations. For the sand/gravel, the rough approximation gave conductivities of 2.5 mS/m for both the 100 MHz data and the 200 MHz data, slightly lower values than those obtained using the smooth  37  approximation. The conductivity estimates for the clay-rich till using the rough approximation were 22 mS/m for the 100 MHz data and 45 mS/m for the 200 MHz data. These values are again only slightly lower than the conductivity estimates using the smooth approximation. It is evident that the choice of a rough or smooth approximations will not drastically affect the results. The smooth approximation will result in our obtaining a greater electrical conductivity than would be obtained with the rough approximation; this is because any possible loss due to scattering off a rough interface is attributed to conductive loss. As we have elected to set backscatter gain to zero, which attributes all loss due to this parameter to conductive loss, we decide, for consistency, to use the smooth approximation. The combination of zero backscatter gain and the smooth approximation will ensure that we are obtaining an upper bound on electrical conductivity. By far the most ill-constrained parameter in this method is the estimation of G which quantifies the losses from energy splitting at interfaces. One bound on the conductivity can be obtained by setting G=0, as we have done for the previous calculations. This will result in our obtaining the maximum conductivity for the data set since all the power loss associated with splitting at interfaces is attributed to conductive losses. An upper bound on G can be found by doing the opposite: setting the conductivity to zero, thereby attributing all conductive loss to the backscatter gain. Using the technique of setting first G and then s to zero, we can establish bounds on both of these values. These bounds are calculated and summarised in Table 3.3 for both the 100 and 200 MHz data. Table 3.3 : Range in G and a values. sedimentary unit  range in G dB/m (100 MHz)  range in G dB/m (200 MHz)  sand/gravel clay till  0-2.5 0-21  0-3.1 0-21  range in a mS/m (100 MHz) 0-3.5 0-27  range in a ms/m (200 MHz) 0-4.3 0-51  For the 100 MHz data the maximum value for G is 2.5 dB/m for the sand and gravel and 21 dB/m for the till. For the 200 MHz data the maximum G values are slightly higher: 3.1 dB/m for the sand and gravel and 21 dB/m for the till. To gain some insight into the difficulty in constraining G, let us consider the type of subsurface variation that could produce G values of 0, 3.1 and 21 dB/m. Recall that G is determined by both the number of boundaries and the contrast in dielectric properties across the boundary. A value for G of zero, means of course that there are no reflecting surfaces at all in the subsurface. A backscatter gain of 2.5 could correspond to an infinite number of combinations of number of boundaries and dielectric contrasts. For example, a sedimentary structure consisting of two boundaries  38  with porosity contrasts of 8% results in a backscatter gain of 2.5 dB/m as does a structure consisting of 16 boundaries with porosity contrasts of 1%. Similarly, a backscatter gain of 21 dB/m could be the result of a structure with 4 boundaries and porosity contrasts of 17% or, equally likely, it could be the result of a structure with 40 boundaries with porosity contrasts of 3.5%. We therefore see that it is difficult to obtain a physical picture of the environment corresponding to a given backscatter gain, since an infinite number of sedimentary structures could give the same result. For this reason we do not attempt a sedimentary interpretation of the backscatter gain results. The bounds established for the electrical conductivity of the sand and gravel are 0.0 to 3.5 mS/m for the 100 MHz data and 0.0 to 4.3 mS/m for the 200 MHz data. For the clay-rich till the resulting bounds on conductivity are much wider: 0.0 to 27 mS/m for the 100 MHz data and 0.0 to 51 mS/m for the 200 MHz data. The key to tightening the bounds on electrical conductivity is to constrain the estimate for backscatter gain, G. To accomplish this, information is needed about the sedimentary structure being imaged. Unfortunately, such detailed information is only available through direct sediment sampling. It may be possible, however, to obtain some general insight on the expected range in backscatter gain in different sedimentary environments through cliff face studies to examine the relationship between dielectric constant and sedimentary boundaries.  As we are proposing this technique as a means of  obtaining rapid estimates of electrical conductivity over a large area, we feel that the most useful approach is to assume that G is zero, and calculate the maximum electrical conductivity for the sediments.  GPR RESULTS OVER TWO LAYERS The attenuation analysis technique was further tested in a more complex setting: a sand and gravel pit with a subsurface structure composed of two distinct layers. Trenching data from the pit show that the sediments in the upper layer are composed of sand, silty sand and gravel and the lower layer is a blue-grey soft plastic clay. A 100 MHz radar survey was collected from west to east with the result shown in Figure 3.6. The data consist of a series of steeply dipping reflectors, with penetration depth varying from about 15 m at the western edge of the survey to 3 m on the eastern end. Excavation has lowered the pit floor to the level of the water table. We therefore assumed that the area imaged in the radar survey is fully saturated. The velocity of the electromagnetic waves was determined to be 0.063 m/ns from a CMP survey conducted in the pit.  39  An example of the F versus L plot from one radar trace at this site is displayed in Figure 3.7. Note that the slope from about 1.2 to 3.5 m is markedly less than that between 3.5 and 4.2 m. Therefore without any supporting sedimentary data we would expect there to be two different sedimentary units, each with a different electrical conductivity, in the region imaged in this radar trace. In this case, however, we do have sedimentary information provided by trenching and drilling data which indicates that the two distinct attenuations correspond to the sand and gravel and to the thick layer of massive clay which underlies it. D i s t a n c e (m)  W 0  10  20  30  E 40  50  Figure 3.6: Radar data collected over a 2-layer earth. The sediments in the upper layer are composed of sand, silty sand and gravel. The lower layer is a blue-gray soft plastic clay. The clay layer dips at an angle of about 25 degrees and outcrops at one end of the survey line. -3  i  1  1  1  1  1  -3.5 --  -4 - -  - 4 . 5 --  -5 --  .5.5  J  1  Depth (m)  Figure 3.7: A plot of the function F for a trace at 45 m along the radar survey shown in Figure 3.6. There are two distinct slopes to the function F: one corresponds to the sand (between 1 m and 3.5 m) and one to the clay (between 3.5 m and 4.2 m).  40  In total ten traces were extracted from the data set at equally spaced intervals. In our calculation of F, we assumed a smooth approximation for the cross-sectional area because the largest grain size noted in the trenching and drilling data was gravel sized. Again we assumed the backscatter gain was zero in order to obtain the upper bounds on the conductivity estimates. Two attenuation constants were obtained from each trace, one for the portion of the trace interpreted as corresponding to the sand beds, and one from the region interpreted corresponding to the clay-rich till. The attenuation constants from the sand region of each trace were averaged to obtain 0.09 nepers/m with a standard deviation of 0.05 nepers/m. Similarly the attenuation constants from the claytill region were averaged to obtain an attenuation constant of 1.7 nepers/m with a standard deviation of 0.4 nepers/m. The resulting conductivity estimate for the sand and gravel is 2.2 mS/m and that for the clay is 42 mS/m. The sand and gravel result compares well with the conductivities of 2.8 to 5.0 mS/m obtained from a cone penetrometer survey at the site (Campanella et al, 1994). Laboratory measurements of the electrical conductivity of these sediments were not conducted. However the expected range of electrical conductivities of similar soil types can be obtained from published results from other studies. For example Telford et al. (1976) suggest a range of 10 to 1000 mS/m for clays, and 1.25 to 100 mS/m for alluvium and sands.  CONCLUSION In this study we have introduced a method for obtaining estimates of electrical conductivity from analysis of the attenuation of GPR data. We have applied this technique to a radar data set collected over a lithologic boundary with 100 and 200 MHz antennas. The conductivities obtained for the sand and gravel material were 3.8 mS/m and 4.6 mS/m for the 100 and 200 MHz data. The corresponding results for the clay-rich till were 27 mS/m and 51 mS/m. In these calculations we assumed no loss from backscatter gain, and therefore obtained estimates of the maximum conductivity of the geologic material. These estimates compare very well with DC resistivity data collected along the same line as the GPR survey. A second test of the decay rate technique was carried out using data obtained over a site with two distinct layers: a clay layer overlain by sands and gravels. The attenuation analysis resulted in two distinct conductivities of 42 and 2.2 mS/m for the clay and sand respectively. The largest source of error in this method arises from the estimation of backscatter gain. Setting this value to zero provides an upper limit on the electrical conductivity. Improvements to this estimation 41  can be made with further laboratory and field studies to quantify dielectric changes across sedimentary boundaries. From these results we conclude that reasonable estimates of electrical conductivity can be obtained from the analysis of GPR attenuation. The results obtained using this method are most accurate in more resistive regions where tighter bounds can be placed on the electrical conductivity estimates. Using analysis of GPR attenuation, conductivity data can be obtained rapidly and efficiently over a large area. We suggest that this method can be used to provide data for improving groundwater models.  42  CHAPTER 4: Determination of Hydraulic Boundaries Using Ground Penetrating Radar  INTRODUCTION The first step in developing a groundwater model is to define the boundaries that mark changes in hydraulic properties within the region. These boundaries include the boundary between aquifer and aquitard materials, the boundary between the saturated and unsaturated zones and boundaries within the aquifer that mark changes in hydraulic properties. At present the data needed to identify and map these boundaries are mainly lithologic information obtained from drill holes, and direct measurement of hydraulic properties obtained from pump and slug tests. These methods of data acquisition require large commitments in time and money, yet still may provide only limited data coverage due to the large areas involved. The result is that groundwater models are often based on very simplified structures. An example of this type of simplification is the investigation of the Fraser Lowlands in southwestern British Columbia (Halstead, 1986).  The region was divided into what were called  "hydrostratigraphic" units defined as regions with similar hydraulic properties based on the lithology, permeability and porosity (Halstead, 1986) using information obtained from drill data. The resulting units are on the order of 40 to 100 m thick and cover areas of up to several 100 km . They are useful for 2  delineating large scale variations in hydraulic properties, however, this gross structure is not sufficient to fully characterise the region, particularly for small scale problems such as the prevention and remediation of groundwater contamination. A solution to the problem of data acquisition may lie in a geophysical method such as ground penetrating radar (GPR), which has the capability of imaging electrical changes in the subsurface rapidly, at relatively low cost and with a high resolution. GPR has been used in numerous groundwater studies. These include mapping of the boundary between the saturated and unsaturated zones (Greenhouse et al., 1987; Vellidis, 1990; van Overmeeren, 1994), determining depth to bedrock (Davis and Annan, 1989), mapping clay/sand boundaries (Young, 1995), imaging sedimentary units (Greenhouse et al., 1987; Beres and Haeni, 1991; Jol and Smith, 1991; Rea and Knight, 1994; vanOvermeeren, 1994; Huggenberger, 1995) and determining aquifer heterogeneity (Olhoeft, 1991; Gawthorpe, et al., 1993; Rea and Knight, 1996). Our goal in this study is to define and test methods for using GPR to identify large scale hydraulic boundaries within a region for use in groundwater modelling. The top of the saturated zone (TSZ) may be the most straightforward boundary to identify in GPR data as it is roughly horizontal and 43  also marks a large contrast in properties. The other hydraulic boundaries have more complex geometries and property contrasts. We therefore identify these boundaries through the use of hydrogeologic units. We define hydrogeologic units as sedimentary units that can be considered homogeneous with respect to their hydraulic properties. Hydraulic boundaries are therefore defined as the interface between two hydrogeologic units.  The concept of hydrogeologic units is based on the assumption that the  distribution of subsurface hydraulic properties is closely related to sedimentary architecture, fabric and texture, which in turn reflects the sedimentary depositional environment. Anderson (1989) was one of the first to suggest this approach for groundwater and contaminant transport modelling. • She adapted concepts from sedimentary facies modelling to identify hydrogeologic units in glacial and glaciofluvial environments. The hydrogeologic units are then used to develop a regional groundwater model. This approach is particularly suited to GPR since the high resolution data obtained in GPR surveys can be used both to identify sedimentary features and to obtain information about the internal properties of the sedimentary features. In this paper we present aframeworkfor using GPR data to constrain groundwater models through the identification of the TSZ and hydrogeologic units. Our goal is to use GPR to improve the resolution represented by the hydrostratigraphic units defined in the area. In order to improve our understanding of GPR interpretation, we have conducted two preliminary studies. The first is a study to investigate the GPR response to the TSZ. More specifically we study different processing techniques that can be used to enhance the TSZ reflection in a radar survey. The second preliminary study was a series of cliff face studies that addressed specific issues in the interpretation of GPR data to identify hydrogeologic units. The interpretation techniques developed using the TSZ and cliff face studies were tested on a GPR data set collected over the Brookswood aquifer in south-western British Columbia.  THE INTERPRETATION OF GPR DATA In order to identify hydrogeologic units and reflections from the saturated zone in GPR data, we must first assess the capabilities of GPR to outline these boundaries in the subsurface. This involves a discussion of the basic theory of GPR followed by an attempt to relate sedimentary and hydraulic parameters to the parameters recorded in GPR data. In a GPR survey, electromagnetic waves are transmitted into the ground by an antenna placed on the surface. The response of the electromagnetic waves is governed primarily by two properties of the subsurface material: the dielectric constant and the electrical conductivity. Changes in the dielectric constant cause energy to be reflected back to the 44  surface where it is recorded by a receiving antenna. The electrical conductivity of the earth causes the amplitude of the electromagnetic wave to be attenuated. A radar survey therefore provides an image of the subsurface, with reflections in the radar section marking changes in dielectric constant, and attenuation providing an indication of the electrical conductivity. In order to obtain a more complete image of the ground, a 3D radar survey can be collected. In general this consists of a series of parallel radar sections, which together form a 3D image of the subsurface, showing changes in electrical properties. By electrical properties we mean both the dielectric constant and the electrical conductivity. An image of changes in electrical properties, such as that produced in a GPR survey, is not itself useful for groundwater modelling. Changes in electrical properties are, however, related to changes in both sedimentary and hydraulic parameters, which can be used to identify the saturated zone boundary and hydrogeologic units.  The Top of the Saturated Zone (TSZ) We expect that the TSZ will create a strong reflection in the radar survey due to the large contrast in dielectric constant between water (K=80) and earth materials (K=4 to 12). This reflection can sometimes be identified as a horizontal reflector cutting through dipping reflectors. Therefore a strong horizontal reflector in a radar survey is often interpreted as being caused by the TSZ. However this "rule of thumb" often fails to provide the correct interpretation. In areas of poorly drained soil, for example, the TSZ may not be observed in the radar data due to the large capillary fringe created above the water table. This causes an indistinct TSZ and a gradation in dielectric constant rather than an abrupt change. Thus there is no distinct TSZ reflection. Another possible misinterpretation arises in areas where the stratigraphy is horizontal. Two formations of strongly contrasting dielectric constants can cause a strong reflection that appears similar to that which one expects from a TSZ. We investigate this issue in a study focusing on processing techniques that can be particularly useful in enhancing the TSZ reflection.  Hydrogeologic Units In order to identify the hydrogeologic units in GPR data we must assess the capabilities of GPR to outline hydraulic and sedimentary boundaries in the subsurface.  This involves establishing  relationships between sedimentary and hydraulic parameters and the parameters recorded in GPR data. The usefulness of determining sedimentary parameters is based upon the assumption that a sedimentary 45  unit, having been deposited by a single process, will be composed of grains with similar size, packing and orientation. Since these parameters are related to hydraulic properties, it is reasonable to assume that each sedimentary unit will possess characteristic hydraulic properties, and so can serve as a useful large scale building block for the hydrogeologic model. Hydraulic properties, while assumed to be the same throughout a given unit, may be anisotropic due to bedding structure and grain orientation within the unit. Recent work has improved methods of interpreting radar data in sedimentological terms (Jol and Smith, 1992; Stephens, 1994; Olsen and Andreasen, 1995). Identification of sedimentary features in radar data requires establishing a relationship between sedimentary properties and electrical properties. For this purpose it is useful to recognise that sedimentary classification is based on five fundamental properties of the sedimentary material: grain composition, grain size, grain shape, grain orientation and grain packing. Electrical properties are determined by the composition and geometry of the solid and liquid components of the earth materials. As the composition and geometry of the solid component are clearly related to the five sedimentary properties, it is reasonable to assume that a change in sedimentary properties will correspond to a change in electrical properties, and thus will be observed in the GPR data. A change in sedimentary properties may be manifested in the radar data in two ways. If the sedimentary change corresponds to a change in dielectric constant, a reflection will occur in the radar data; if the sedimentary change corresponds to a change in electrical conductivity, the attenuation of the GPR data will be affected. An important issue in the relationship between sedimentary and dielectric properties is the effect of the liquid component of the material. Water, the most common liquid in the subsurface, plays no part in sedimentological classification, yet due to the high dielectric contrast between water (80) and common minerals (4 to 12), it has a significant role in determining the dielectric properties of a material. The presence of water can therefore cause a radar reflection where no sedimentary boundary exists. The most obvious example of this is the TSZ, which, though clearly not a sedimentary boundary, can produce a very strong radar reflection (Vellidis et al., 1990). In other instances, however, the presence of water can enhance reflections from sedimentary boundaries. In the saturated zone the volumetric water content is equivalent to porosity.  A sedimentary boundary defined by a porosity contrast will therefore  correspond to a dielectric contrast, and so cause a radar reflection. In the unsaturated zone, the situation is slightly more complicated. Water tends to remain preferentially in smaller pores due to higher capillary forces. Therefore, in the unsaturated zone, sedimentary boundaries marked by a variation in 46  pore size will cause radar reflections. The same sedimentary boundary in the saturated zone may not create such a large reflection if porosities are similar. The electrical conductivity of the subsurface may also be affected by the presence of water. This introduces similar complications in interpreting attenuation in radar data as discussed in the previous paragraph with reference to the interpretation of radar reflections. relationship exists between electrical conductivity  However, a more important  and clay-content, which, in a fresh water  environment, is the dominant factor determining electrical conductivity. In addition to being an indicator of the amount of clay present, electrical conductivity has been shown in numerous studies to have a strong inverse relationship with hydraulic conductivity in aquifer environments (Urish, 1981; Mazac, 1990). The attenuation of the amplitudes in a GPR section can therefore be used to map interfaces between the resistive sands and gravels of an aquifer unit and the conductive clay-rich sediments of an aquitard unit.  Classification of radar data There are numerous ways of approaching the classification of GPR data. One of the most common approaches to GPR classification is through radar facies analysis (Beres and Haeni, 1991; Bristow, 1995). Radar facies are defined as patterns of reflectors with similar characteristics such as frequency and amplitude (Gawthorpe et al., 1993). These patterns are then related to a depositional sequence and lithology. An alternative classification scheme is radar sequence stratigraphy (Gawthorpe et al., 1993; Bristow, 1995). In a similar way to seismic sequence stratigraphy this method relies on identification of reflector terminations which are interpreted to be non-depositional or erosional events. Identification of such terminations in radar data can be a very ambiguous process due to noise and interference patterns from the TSZ or from reflectors located off the line of the survey.  The benefit of radar sequence  stratigraphy is that it identifies boundaries of units which may also mark hydrogeologic boundaries. In this study we have adopted a classification scheme called radar architectural element analysis (Stephens, 1994), which is an adaptation of a sedimentary classification technique called "architectural element analysis" (Miall, 1985). An architectural element (AE) is a 2D or 3D sedimentary package described by its bounding characteristics (e.g., erosional or gradational), external geometry, scale, and internal characteristics. Following the definition of a sedimentary AE, a radar AE can be defined as a 2D or 3D package of radar reflections described by its scale and geometry, or orientation, as well as the 47  general character of reflections within the element. This concept is well suited to our goal of defining hydrogeologic unit in the subsurface because both the sedimentary and hydraulic properties of a hydrogeologic units contribute to the definition of a radar AE. A radar AE is linked to the sedimentary structure through the scale, geometry and orientation of the reflections and to the hydraulic properties through the overall character (attenuation).  We suggest, therefore, that radar AEs be used in a  groundwater model as hydrogeologic units. There are several major differences between radar and sedimentary AEs which must be clearly understood. First, radar surveys image changes in sedimentary parameters, and do not give actual values of these parameters. Therefore two quite different sedimentary features, similar only in their internal structures, could produce similar results in a radar survey, and therefore be classified as the same radar element. In addition, some sedimentary changes cannot be imaged at all due to the very small difference in electrical properties between sedimentary units. A radar element could therefore contain within it several sedimentary elements. Following a similar argument, it is possible that a single sedimentary AE could contain within it several radar AEs.  STUDY TO INVESTIGATE T H E TOP OF T H E SATURATED ZONE  A key problem in identifying the TSZ is determining whether a strong horizontal reflector corresponds to the TSZ or to a sedimentary boundary. The goal of the TSZ study was to determine whether GPR could be used to positively identify the TSZ in GPR sections and also to determine which data processing techniques can be used to enhance the TSZ response.  In this section we present two  radar sections in which a strong horizontal reflection exists. In one case (site 1) this reflection is not due to the TSZ, and in the other (site 2) it is.  Figure 4.1: A schematic of the PulseEKKO IV GPR unit. 48  Data were collected using a PulseEKKO IV GPR unit with 100 and 200 MHz antennas. This system, shown in Figure 4.1, consists of a transmitting antenna, a receiving antenna and a control box linked by fibre optic cables.  Acquisition is controlled by a computer which allows the data to be  displayed as they are collected. Site 1: Environment Canada Test Site The first site is at the Environment Canada test site at 500 Clearbrook Road in Matsqui, British Columbia. A GPR line 20 m long was collected at this site with 100 MHz antennas. In the data from site 1, a continuous, horizontal reflector is clearly visible at a depth of 4 m along the length of the survey (Figure 4.2). Other horizontal and slightly dipping reflectors are visible to a depth of about 6 m. Distance (m)  Figure 4.2: GPR data from Site 1. The lithological data from the wells indicate a sand and gravel composition from 2 m to 36 m depth. Variations over this length are due to changes in grain size. Data from a piezometer located within 5 m of the radar line indicated that the TSZ was at a depth of 18 m. Therefore, the continuous horizontal reflection at a depth of 4 m is not attributable to the TSZ.  Site 2: Stokes Gravel Pit The second site is the Stokes gravel pit at 24  th  and 192 in Langley, British Columbia. Data at nd  the Stokes pit, shown in Figure 4.3, were collected along a 20 m long line on a cliff 2 meters above the pit floor. A strong reflector is visible at a depth of 2.5 m along the length of the survey. Shorter horizontal reflectors are visible throughout the rest of the survey, down to a depth of about 5.5 m. A second survey of this line was collected 2 weeks later in order to compare the TSZ response. 49  Distance (m)  Figure 4.3: GPR data from Site 2. The cliff is composed of coarse sand to cobbles. Previous GPR surveys in this pit have identified large sand and gravel foreset beds in this area along with a conformable clay layer. The gravel pit has been excavated down to the level of the TSZ. The TSZ level is therefore well known as the level of ponding in trenches in the pit floor. The reflection at a depth of 2.5 m is therefore due to the top of the saturated zone.  Processing Techniques to Identify the TSZ Our goal in this section is to investigate processing techniques to differentiate between the "true" TSZ reflection at Site 2 and the false T S Z reflection at Site 1. In our analysis of these data sets we investigated four techniques for enhancing and identifying the T S Z reflections. These are: gaining the data, instantaneous phase and frequency, subtraction of data lines and C M P analysis.  Gaining of the Data The type of gain applied to a data set can dramatically change the way a data set is interpreted. A n A G C gain (automatic gain control) amplifies all reflections to the maximum possible level in a given window. Therefore a data set to which A G C gain has been applied has no variation between "strong" and "weak" reflectors. A reflection from a TSZ may have the same appearance as that from an interface between two materials which have only a small dielectric contrast.  A spreading and exponential  compensation (SEC) gain removes the effect of spherical spreading and ohmic losses. For SEC gain, the data processor supplies an estimated attenuation from ohmic losses and velocity of the electromagnetic waves. 50  The data in Figures 4.2 and 4.3 have had an A G C gain applied. Figures 4.4a and 4.4b show the data from site 1 and site 2 with an SEC gain applied. At site 2, the SEC gain enhances the 2.5 m reflection and clearly marks this is a significant reflection. The only other features in this plot are two "reflections" in the top half meter. These are not actual reflections but are direct arrivals through the air and ground from the transmitting to the receiving antenna. The SEC gain applied to the site 1 data does not show a similar differentiation between reflections.  There is a reflection at a depth of 4 m that is  continuous over the first 12 m of the survey, but no other significant features are apparent. Distance (m)  5  S  10=0 + illllillilllllllllllilllllllllliilllllllllllllllllllll Figure 44 .b:  G P R data from Site 2 with SEC gain applied.  Figures 4.5a and 4.5b show the ungained data. The reflection at 2.5 m is still faintly visible in the Site 2 data even without amplification. The horizontal reflection from site 1 seen in the A G C and SEC gained sections is not visible.  51  Distance (m)  Figure 4.5a: GPR data from Site 1 with no gain applied. Distance (m) IS _  IS •  IS CM  | 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 I I  §5.  0  1 0. 0  +  I  Figure 4.5b: GPR datafromSite 2 with no gain applied. The type of gain applied to the data is critical to interpretation. Using an SEC gain may be the best way to discriminate between TSZ and geologic reflections.  The ungained sections may also be a useful  visualisation technique, if the TSZ is not too deep.  Instantaneous Frequency and Phase The software package EKKO Tools can calculate the instantaneous phase andfrequencyplots for radar traces. Briefly, the instantaneous phase of the radar data tends to emphasise continuous reflection, even those which have a lot of interference obscuring them. Therefore it may be a means of extracting the continuous TSZ reflection when geologic reflections are obscuring it. The instantaneous frequency of the radar response will change at interfaces such as the TSZ (Davis and Annan, 1989) and therefore is also investigated as a means of enhancing the TSZ reflection.  52  Figure 4.6a shows the instantaneous phase plot of data from the Stokes gravel pit. The strong continuous feature at 2.5 m depth is caused by the TSZ. The continuous feature at 1 m depth is the ground arrival. Figure 4.6b is an instantaneous frequency plot of the Stokes data. There is a distinct, continuous line of spikes at about 2.5 m depth which correspond to the TSZ depth. The first two spikes at a depth of 0.5 and 1.0 m are caused by the air and ground arrivals. Figure 4.7a is an instantaneous phase plot of a line from Site 1. There is a strong continuous feature at a depth of 6 m which looks similar to the feature in Figure 4.6a that was cause by the TSZ. However, the instantaneous frequency plot in Figure 4.7b does not show the sharp spikes observed at the TSZ depth in Figure 4.6b. Distance (m) 3 [ I  CM 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 II  Figure 4.6b: Instantaneous frequency of GPR data from Site 2.  53  I I I I I I [|  Distance (m)  Figure 4.7a: Instantaneous phase plot of GPR datafromSite 1.  Distance (m)  Figure 4.7b: Instantaneousfrequencyplot of GPR datafromSite 1. Instantaneous phase andfrequencyplots used together may be used as a discriminator between TSZ reflections and geological reflections. However, work still must be done to confirm the validity of this method. Subtraction of Data Lines One way to investigate variations in the TSZ is to subtract data sets collected on the same line over a period of time, as we have done in this survey. Assuming that the TSZ has changed, and that the surveys have been conducted at exactly the same position, the reflections from geologic features will cancel out, leaving only the change in TSZ reflection. We have attempted this technique with our data.  54  Figure 4.8c shows the difference between two data sets collected at the Stokes pit two weeks apart (Figures 4.8a and b). There is a strong reflector visible in the first half meter. This is believed to be due to a change in the character of the ground arrival due to changes in the surficial water content. There is a parallel set of reflectors at a depth of about 2.5 m. We interpret this as being due to the change in position of the TSZ reflection.  Distance (m)  2  $  Distance (m)  2  $  Figure 4.8: a) SEC gained GPR data collected at site 2. b) SEC gained data from site 2 collected 2 weeks later, c) subtraction of data in Figures 4.8a and 4.8b. The difficulty with the subtraction method is the necessity to collect the data lines at exactly the same spot. Slight variations in antenna positioning can vary the reflected signal enough to render this 55  method useless. However if the antennas are position correctly, the subtraction of two data sets should give not only changes in the TSZ, but also changes in water content in the unsaturated zone.  CMP analysis A common midpoint (CMP) survey may be carried out and analysed to determine the average velocity of radar waves through the subsurface. Such surveys are time consuming but essential if a time to depth conversion is required. The velocity in water saturated sand is 0.06 m/ns, quite different from that in dry sand, which is 1.5 m/ns. Therefore the velocity determined from a reflection above the TSZ should give a higher velocity than one from below the TSZ. So theoretically a TSZ could be identified by comparing the velocities above and below the reflector. The results of the CMP survey from site 1 and site 2 are summarised in Table 4.1.  Velocities  were calculated from both the direct ground arrivals and from hyperbolic reflection using the procedure outlined in Annan and Cosway (1992). The velocities obtained from site 1 are between 0.0543 and 0.091 m/ns; velocities from site 2 are between 0.0732 and 0.1005 m/ns.  Table 4.1: CMP survey results Site  Date  1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2  17/01/94 17/01/94 17/01/94 17/01/94 17/01/94 09/03/94 09/03/94 09/03/94 02/02/94 02/02/94 02/02/94 14/02/94 14/02/94 14/02/94 14/02/94 14/02/94 09/03/94 09/03/94 23/03/94  Frequenc y 100 100 200 200 100 100 100 200 100 200 200 100 100 100 200 200 100 200 100  Depth of reflector (m)  Velocity (m/ns)  1.51 2.57 GROUND 2.59 2.98 1.17 4.61 1.33 1.96 GROUND 1.92 GROUND 2.06 3.61 1.95 3.59 2.18 2.36 1.93  0.075 0.069 0.091 0.076 0.054 0.078 0.069 0.083 0.085 0.082 0.085 0.091 0.087 0.083 0.082 0.073 0.093 0.101 0.077  Theoretically CMP surveys should be an ideal method of confirming the TSZ location. However, there are many drawbacks. In a wet area, after a rainfall for example, the "unsaturated" zone 56  may have enough water in it to make its velocity similar to that in the "saturated" zone. Phase shifts in the data can provide a challenge in choosing the hyperbolic reflection points in a CMP analysis which may therefore produce misleading results. CMP analysis becomes more difficult for deeper reflectors as phase shifts and other factors interfere with the ability to "pick" the reflector. The results of the CMP data are inconclusive. Our data do show a decrease in velocity in site 2 above and below the TSZ.  This is an expected result because velocity slows with increasing water  content. However the velocities vary between CMPs carried out at the same time, so the accuracy of the CMP data is in question. In addition, velocities obtained from site 1 also decrease with depth by about the same magnitude as those in site 2.  CLIFF F A C E STUDY In order to test the concept of radar architectural element analysis we conducted a series of cliff face experiments which allowed direct comparison between the radar data and the sediments being imaged. The cliff face study involved three cliff faces in the Kirkpatrick gravel pit in south-western British Columbia. The sediments in the pit form part of a deltaic sequence deposited approximately 13,000 years ago as part of the Fort Langley Formation (Tribe, 1994). The pit has very good exposures of a variety of sedimentary features, and so provided a good opportunity for visual comparison of radar data with the sediments being imaged. Before an interpretation was carried out, the radar data were processed to create a clear and undistorted image. The basic processing steps included a correction for time zero drift and a "dewow" process to correct for signal saturation. An elevation correction was applied to adjust for topographic variations. The final step in processing the radar data was the application of a time variable gain to the data. In this study we used an automatic gain control (AGC) that amplifies all reflectors so that they have equal strength. This gain function is useful for visualising all reflectors to determine their overall structure and trend.  Data The first cliff, oriented from north to south, is shown in the photo-mosaic in Figure 4.9a. Note that there is some distortion of the photo-mosaic due to the variation in the angle at which the photographs were taken. The northern section of the cliff consists of horizontally bedded sand sheets of two interlayered facies: one is a coarse to medium sand, the other is a fine sand with silt and some clay. 57  Although not identifiable in the photographs, well-defined climbing ripple structures can be seen in the cliff, preserved by overlying mud and silt drapes (Tribe, 1994). Fluid escape structures on the order of centimetres are evident throughout the cliff. At a distance of 39 m along the survey line, the sand sheets become severely faulted and deformed. Overlying the deformed sand sheets is a poorly sorted, weakly stratified matrix rich conglomerate (Tribe, 1994).  The sand sheets are interpreted as subaqueous  turbidite beds and the conglomerate unit is interpreted as a debris flow (Tribe, 1994). The three main features seen in the visual inspection of the cliff face can also be distinguished in the radar survey, shown in Figure 4.9b, and are identified as three separate radar architectural elements. An interpretive schematic of the radar survey is shown in Figure 4.9c. The first radar AE is identified by the parallel, slightly dipping reflections which occur between about 0 and 25 m along the radar survey. This element corresponds to the turbidite beds. A second radar AE is defined by the series of dipping and trough shaped reflectors occurring at a depth of 2 to 8 m from 30 to 70 m along the survey. These reflectors correspond to the deformed turbidite beds. The third radar AE can be identified in the upper 3 m of the radar survey due to the termination of reflectors from the sand sheets, and to the different orientation of reflectors from those in the surrounding data. This element corresponds to the debris flow overlying the deformed turbidite beds. Small features such as fluid escape and ripple structures are not resolved in the radar data, and therefore do not play a part in the identification of the radar AEs. Similarly, it is not possible to discern the vertical fractures in the deformed sediments (element 2), nor to discern that this feature is a deformed equivalent of the turbidite beds further to the north (element 1). A second cliff face oriented in an east-west direction in the north end of the pit is shown in Figure 4.10a. Large sand and gravel clinoforms with dips of between 5 and 35 degrees (Tribe, 1994) are exposed in this cliff.  The clinoforms consist of normally graded packages which are interpreted as  foreset beds. Radar data were collected along the cliff top using 100 MHz antennas with the results shown in Figure 4.10b. The radar data consist predominantly of dipping reflections with several isolated concave-up reflections.  The concave-up reflections contain within them faintly discernible dipping  reflections. The large dipping reflectors and the concave-up reflections with internal dipping reflections are identified as separate radar AEs. The dipping reflections correspond well to the orientation and variability seen in the foreset beds in the photograph. However, the concave-up elements do not appear to have corresponding sedimentary features in the cliff face. The shape of this element is suggestive of a channel,  which  is  an  expected  sedimentary 58  feature  in  a  deltaic  environment.  These  Figure 4.10: a) Photograph of Cliff2 b) A GPR survey over Cliff2. c) Interpretive diagram of Cliff showing the different radar AEs.  60  inferred channels which cause the concave up reflections may be obscured by the weathering of the cliff face, and so were simply not observed in the sedimentological analysis of the cliff. An interpretive diagram of the radar image is in Figure 4.10c.  Figure 4.11a: Photograph of Cliff 3. Horizontal scale is 20 m; vertical scale is 8 m. The third cliff face (Figure 4.11a) is oriented in the north-south direction and shows the same turbidite beds as in the first cliff described above. An additional feature of note is a large scour surface which can be seen running from the upper left to the middle right corner. Sedimentary structures and facies on either side of this scour surface are identical.  0  Distance (m)  10  20  Distance (m) 0  10  20  Figure 4.11c: A 200 MHz GPR survey over ClifB. Station spacing is 0.2 m. Both 100 MHz and 200 MHz radar data were collected over this cliff. The results are shown in Figures l i b and 11c. Once again the orientation of the turbidite beds is imaged in the radar data as slightly dipping reflections. These reflections are identified as the same architectural element as in the first cliff face. Note that the number of reflections per meter is not the same in the 100 and 200 MHz data sets. Also note that the scour surface, which is clearly visible in the photograph, is not visible in the radar data.  Discussion From these cliff face studies we'were able to identify a variety of radar architectural elements. Not all aspects of the definition of a radar A E were used in each case. For example, the boundaries between elements 1 and 3 in the first cliff are indistinct and so the external shape is not part of the identification criteria; nor can the scale of these elements be fully defined, as the radar survey images only a portion of the complete unit. We therefore base our identification mainly on the orientation and pattern of the radar reflections. An important observation from the study is that the same radar element may be caused by different sedimentary features. For example, the concave-up reflectors in the first and second cliff face surveys were caused by very different sedimentary features: deformed turbidite beds and channels. On close inspection, however, subtle differences are observable. In the first cliff, the reflections within the 62  concave-up element are also concave-up and conformable with the outer reflection. In the second cliff, however, the internal reflections are dipping and not conformable with the outer concave up elements. It is possible therefore to subdivide the concave-up element into two elements based upon the nature of the internal reflections. The cliff face results also clearly show that not all sedimentary boundaries can be identified in radar surveys. Sedimentary boundaries that are particularly clear in the radar record are those that mark a change in internal structure (e.g. bedding orientation), such as the contrast between element 1 and 3 in the first cliff, or a change in composition, such as between element 1 and 2 in the first cliff. An example of a sedimentary boundary that cannot be identified in the radar record occurs in the third cliff face. In this cliff the sedimentary structure and material on either side of the scour is identical. The scour surface is therefore invisible in the radar record even though it is clearly seen in a visual inspection of the cliff. The ability to resolve fine scale horizontal bedding is related to the frequency of antennas and to the bedding thickness. The resolution limit is defined by the Rayleigh limit of A/4, where A, is the wavelength. Energy reflected from bedding surfaces with a smaller separation distance than this limit will arrive at the surface at intervals so closely spaced that resolution of the different reflections is impossible. An illustration of this point is the radar surveys over the third cliff face, collected with 100 and 200 MHz antennas. The number of reflections is different in the 100 and 200 MHz radar sections and neither corresponds to the number of sedimentary boundaries present in the cliff face. This is a good illustration of the fact that the spacing between reflections in radar sections will not in general correspond to the number of sedimentary boundaries present. There are several major differences between radar and sedimentary AEs which must be clearly understood. First, radar surveys image changes in sedimentary parameters, and do not give actual values of these parameters. Therefore two quite different sedimentary features, similar only in their internal structures, could produce similar results in a radar survey, and therefore be classified as the same radar element. In addition, some sedimentary changes cannot be imaged at all due to the very small difference in electrical properties between sedimentary units. A radar element could therefore contain within it several sedimentary elements. Following a similar argument, it is possible that a single sedimentary AE could contain within it several radar AEs.  63  T H E BROOKSWOOD AQUIFER STUDY Having investigated the identification of the TSZ and radar AEs in environments with good ground truthing available, we proceeded to apply the interpretation techniques to 12 km of GPR data collected over the Brookswood aquifer. The Brookswood aquifer is a shallow, unconfined aquifer in the Fraser Lowland area of south-western British Columbia (see map, Figure 4.12a). It is composed of sediments of the Sumas Drift Formation deposited during the last glacial retreat 11,500 to 11,000 years B.P. The meltwaters from the Sumas Ice deposited their sediment in the sea at this point, forming a large fan-type delta. The sand and gravel of these deltaic sediments form the material of the Brookswood aquifer. 236'  237'  238'  i  ,  i  239'  l_  cn"  Figure 4.12a: The Fraser Lowland in south-western British Columbia showing the location of the Brookswood aquifer. The sedimentary architectural elements found in and around the Brookswood aquifer reflect the depositional environment in which the aquifer was formed. Deposition of the Brookswood sediments has been attributed to a river-dominated delta (Ricketts and Jackson, 1994). Evidence for this are the numerous fluvial facies found exposed in gravel pits in the area.  These consist mainly of coarse  sediments in braided channel facies. Some evidence has been found of estuarine accretionary foresets which may suggest a tidal influence (Ricketts and Jackson, 1994). However, no fossil evidence has been found to suggest major marine reworking of sediments (Ricketts and Jackson, 1994). It is therefore expected that many fluvial elements of the kind described by Miall (1985), are present in the Brookswood aquifer: channel, lateral accretion macroforms, downstream accretion macroforms,  64  sediment gravity flow, laminated sand sheets, gravel bars and bedforms, sandy bedforms and overbank fines. In addition to thefluvialfacies expected in a river-dominated deltaic environment, one also expects to find foreset beds, formed where rivers and streams enter the deeper water of the ocean. Foreset beds occur at dips of 10 to 45 degrees and consist of sands and/or gravel in a repeating sequence of fining upwards sediments. A final element which is found in a glaciofluvial environment is glacial till. Till is defined as poorly sorted unstratified sediments deposited directly from glaciers or ice sheets (Eyles and Eyles, 1992). Tills have a wide range of grain size and composition, often including large amounts of clay. Sediments surrounding the Brookswood aquifer are important as they mark the limit of the aquifer. These include the lodgement and flow tills of the Fort Langley Formation in the east and glaciomarine clays of the Capilano sediments in the west.  ,  !  32nd Ave.  / , ]  { |  j  ^  '• il'  !l 19  » , *  i .  1  ...  ?!  (O  "i*  '  2nd  fkf  %  f  CO  [i -""  »  ;  l  24h Ave »•  .-"  •  ' -.  1  ' l  . f  C  jj. •••• ...  Figure 4.12b: Hydrostratigraphic units of the Brookswood area (from Halstead, 1986). The yellow unit corresponds to sediments of the Brookswood aquifer itself. The blue, red and brown units are aquitard units  There are four hydrostratigraphic units that have been defined for the region in and around the Brookswood aquifer positioned as shown in the fence diagram in Figure 4.12b. Each hydrostratigraphic  65  unit identified by Halstead (1986) may be composed of one or more of the different sedimentary AEs discussed above. The sediments of the Brookswood aquifer itself are classified as hydrostratigraphic unit C, shown in yellow in Figure 4.12b, described as "glaciofluvial sand and gravel, deltas, locally clay covering deltas" (Halstead, 1986). The Brookswood is underlain by Unit E, "mostly marine sediments interbedded with estuarine and fluvial deposits" and shown in brown in the fence diagram. The aquitards surrounding the Brookswood are classified as unit A, defined as "clay, stony clay and silty clays" (Halstead, 1986), unit B, defined as "glacio-marine stony clays" and a small area of Unit D, defined as "till and diamictons" on the eastern boundary. Unit A is shown in light blue, Unit B in dark blue, and unit D in red in the fence diagram in Figure 4.12b. Our goal in this study is to use GPR data to improve upon the interpretation shown in the fence diagram by mapping hydraulic features at a smaller scale. These features include the TSZ and hydrogeologic units.  Data Acquisition Twelve kilometres of radar data were collected over the Brookswood aquifer using a PulseEKKO IV GPR unit with 100 MHz antennas. Most radar lines were collected along the shoulder of roads, although an intensive survey was collected in a gravel pit in the area. A 3D survey was also collected in this gravel pit. Common midpoint (CMP) surveys were collected to determine the velocity of the radar waves and enable calculation of depth to reflectors. Data were processed using the same techniques used for the cliff face study.  Data and Interpretation Radar AEs in the Brookswood Our initial analysis of the Brookswood data set involved identifying and defining radar AEs which were later used to subdivided the radar sections into distinct units. Given the environment which deposited the Brookswood sediments, we expected to find several different radar AEs in the GPR data. The foreset beds, characteristic of a deltaic environment are expected to give rise to packages of dipping reflections corresponding to the orientation of the foreset beds. An example of this type of reflection package, the dipping radar element, is shown in Figure 4.13a. Another significant feature in the Brookswood region is a channel which can act as high permeability conduit for groundwater or contaminant flow. We expect to see this in the radar record as channel shaped reflections with internal reflections which indicate the orientation of sediments deposited 66  within the channel.  An example of the channel radar element is shown in Figure 4.13b. An  interpretative diagram below the radar data indicates the location of the channel in the shaded region.  c)  Distance (m)  Figure 4.13: Examples of radar elements from the Brookswood data set: a) Dipping radar element b) channel radar element with interpretive diagram c) horizontal radar element d) horizontal radar element e) attenuation element with interpretive diagram. A number of other river and delta based processes will produce sedimentary features with roughly horizontal boundaries at a spectrum of scales. Some examples are sheet sands, chute and pool  67  structures, etc. We expected to see these in the radar record as horizontal and slightly dipping reflections present at many scales. Two examples of horizontal radar elements are seen in Figure 4.13c and 13d. The final radar AE we expected to find in this environment is associated with the clays and clayrich tills. Due to their high electrical conductivity, we expected these sediments to correspond with zones of high attenuation, that is, zones with few or no reflections. An example of the high attenuation radar element is shown Figure 4.13e. In the interpretive diagram below the data, the location of the attenuation element is shown as a shaded region. To summarise, we see four different radar AEs in the Brookswood data set: dipping elements, channel elements, horizontal elements and attenuation elements. The Brookswood data set contains numerous examples of these elements in a variety of forms and associations. As an illustration of how these elements occur in the Brookswood data set, we show four 2D radar lines and one 3D section from the Brookswood data set. A map of the Brookswood aquifer showing the location of the radar lines is shown in Figure 4.14. Drill data are used to provide ground truthing for the radar interpretations. The well numbers correspond to numbers assigned by the Groundwater Section of the Ministry of Environment and Parks, Province of British Columbia.  r  /  '  48th Ave  \  1 : Site of 3-D survey ... Lines of radar survey — Major Roads ' ' \ Approx. limits of V J Brookswood Aquifer 0  / I  r* /  192nd St.  a u  3 / <> / •  1  I  36th Ave.  i  03  ?  1  24th Ave.  EC  3  4  5  '  i  H  IO  &  ro o o  th St  \  2  i  i  28th Ave.  1  Scale (kilometers)  i  \l6th Ave. \  \ ; •  /I  \• \I. x  Figure 4.14: Map of the Brookswood aquifer showing the location of radar sections  Line 24208 Line 24208 is a 400 m long GPR survey collected along 208 Street south of 24 Avenue. We th  th  chose this location to investigate a boundary between aquifer and aquitard hydrostratigraphic units, an 68  aquifer "pinchout"- shown in Halstead's fence diagram (Figure 4.12b). The result of this survey, Figure 4.15a, shows good penetration of the radar from 1 to 340 m along the survey line.  At this point  penetration decreases rapidly until by 400 m there are reflections only in the top 2 or 3 m. Numerous types of reflections occur in this radar section. In the top 5 m between 0 and 300 m there are very long horizontal reflections. In the southern section between 320 and 400 m there are dipping reflections from 0 to about 5 m in depth. Below about 5 m the reflections are of different lengths: a few long reflectors with  short  horizontal  or  slightly  concave  downwards  reflections.  Data  from  well  WW092G.007.2.4.3.019 located 2445 208 Street at the northern end of line 24208 show that the th  sediments are sand and gravel for the first 10 m, followed by 7 m of soft grey clay and gravel. The southern end of the line is composed of boulder till down to a depth of 18 m as indicated by the data from well WW092G.007.2.4.3.016 located at 1941 208 Street. th  The interpretation of line 24208 is shown in Figure 4.15b. Three different radar AEs have been identified in this line. The most prominent is the zone of attenuation between 350 and 400 m along the survey. This we interpret to correspond to the hydrostratigraphic boundary indicated in Halstead's fence a  )  0  100  Distance (m) 200  300  400  Figure 4.15: a) GPR survey of Line 24208 and b) interpretation of Line 24208 showing the radar AEs. The white region represents a horizontal element, the light grey is an attenuation element and the dark grey is a channel element diagram.  In addition, however, we can identify other heterogeneities within the region which is  classified simply as "Unit C" in Figure 4.12.  For example, the dipping reflections adjacent to the  69  "pinchout" are identified as a channel element and the remaining section is considered a single element of horizontal reflections. The reflection from the TSZ was identified by applying an SEC gain to the data. From this visualisation of the data, we interpreted the reflection at a depth of 2 m to correspond to the TSZ. This reflection is marked in the interpretive diagram in Figure 4.15b.  Line 20204 Line 20204 was collected from east to west along 20 Avenue between 204 and 208 Streets in th  th  th  Langley, British Columbia. We chose this line to investigate another aquifer "pinchout" which was inferred to exist in the area based on the fence diagrams of Halstead. We successfully located this pinchout in the data at a distance of 250 m west of 208 Street. A portion of the 500 m line, collected th  directly over the pinchout, is shown in Figure 4.16a. There are minimal reflections below a depth of about 5 m from 0 to 30 m along the survey. At this point radar penetration increases until 50 m along the survey line where reflections are seen down to a depth of 15 m. There is a strong horizontal reflector at a depth of 3 m which is continuous from 50 m to 100 m. Most other reflections in the section are on the order of 5 to 10 m long and are horizontal to slightly concave downwards. An exception occurs between 25 and 65 m between 0 and 7 m depth. Visible in this region is a long, concave up reflection which bounds a region of reflections dipping to the east.  There are two drill holes within 100 m of this line.  Well number  WW092G.007.2.3.1.030 is located at 2027 204A Street near the western end of the radar survey line and consists of sand and gravel, followed by 1 m of sandy gravel with boulders and 13 m of sand and gravel. Well number WW092G.007.2.3.1.041 located at 20633 20 Ave. is just east of the survey line. Data th  from the well indicate that the top 7 m are brown sandy clay followed by 2 m of till and 12 m of stony clay. The interpretation of line 20204 is shown in Figure 4.16b.  There is an element of high  attenuation in this section between 0 and 30 m. We interpret this as corresponding to Halstead's Unit C. Two additional radar AEs can be identified in this line. The concave up reflection between 35 and 65. m marks the boundary of a channel element. The dipping reflections within this element are caused by the sediments deposited when the channel was infilled by lateral accretion. identified as a single element of horizontal reflections.  70  The remaining data are  The reflection at a depth of 3 m is interpreted as being caused by the boundary between saturated and unsaturated sediments. An SEC gained portion of the line is shown in Figure 4.17. This clearly indicates the strength of this reflection compared with others in the data.  72  Stokes Pit The Stokes pit is a municipal gravel pit in Surrey, British Columbia, situated at the corner of 192  nd  Street and 24 Avenue. The Stokes pit has been active for 40 years and at present the top 6 meters th  have been excavated. The Stokes pit lies completely within Halstead's Unit C, which extends down to a  111  !a "2 4 5  4 3  s* rt e cS rt  T 3  S3  U  £s .s u  o a  w g  6 " 2 *  P is  fr ^ (3  G  S3  D • £ *  C  4 3  13  w  «> .2 ^o 23 + S c oo J  4 3  s  rt S3  1)  C IU  6  rt  o o  S3  a a.  Z 3  ^ 3  T3  & K rt  8 rt JT^  qjj  S  4 3  M  a ~ -2 ^£ a OS 4>  S  wo  73  depth of about 50 m in this area. The goal of our investigation in this area was to determine whether GPR could identify significant heterogeneities within this area, which previously was identified as a single hydrostratigraphic unit. Trenching and drilling data in the area provide sedimentological data. A hole drilled near the eastern end of the Stokes pit indicates that the top 3 m is a well-compacted medium sand, followed by coarse sands and gravels. A trench close to the western end of the pit uncovered a 1 m thick layer of soft blue-grey clay, underlain by sands and gravels. The water level is about 1 m below the bottom of the pit. Line S, shown in Figure 4.18a, is a 100 m GPR section collected with 100 M H z antennas at 0.5 m spacing. The interpretation of Line S is shown in Figure 4.18b. We identified three different radar AEs in the data. The dipping reflections are classified as dipping radar AEs. There are three examples of this element in the radar section. The clearest is between 0 and 60 m along the survey and the other two occur between 60 and 115 m and between 80 and 150m. The small package of reflections dipping to the right between the 2  nd  and 3 dipping elements is interpreted as a channel element. The abrupt end to rd  penetration beneath the westernmost dipping element marks the boundary of the third radar A E : a zone of high attenuation element. In the Stokes pit the TSZ is within 1 m of the surface. The air and ground arrivals obscure the top meter of reflections and we therefore do not expect to see a reflection from the saturated/unsaturated interface.  Figure 4.19: a) 3D GPR survey collected in the stokes pit and b) interpretive diagram. 74  In order to investigate the true strike and dip of the dipping radar reflections seen in line S we conducted a 3D survey over the westernmost dipping radar element shown in Line S. The survey is 50 m x 25 m with a radar trace every 0.5 m. In this survey (Figure 4.19) we can clearly see the orientation of the dipping reflections. The strike of the dipping reflections is to the north-west, with the reflections dipping to the west at about 25 degrees. The reflection marking the end of radar penetration can also be clearly seen in this 3D section, a) Distance (m)  Line 192 Line 192 was collected from north to south along 192 Street in Langley, British Columbia. In nd  contrast to the radar sections shown above, these data are difficult to interpret and divide into radar AEs. Several segments from this line are shown in Figures 4.20a and 4.20b to demonstrate some of the variation in quality which can result in a reconnaissance survey such as the one conducted over the Brookswood aquifer. Drill data do not indicate the presence of clay or other highly conductive material along this line; for example at 80 m along the survey, well WW 092G.007.4.1.1.031 consists of sand and gravel from the surface to a depth of 10.8 m. However, despite the lack of obviously conductive material, radar penetration is very poor. Those segments of the line which do allow penetration show a confusion of reflections which are difficult to classify. The cause of the poor quality may be due to surficial sediments of high conductivity, either naturally deposited, or added in the road building process. 75  Another cause of the poor data in the section shown in Figure 4.20a may be the presence of underground utilities running parallel to the road just underneath the survey line.  DISCUSSION This study has shown that the identification of radar AEs in GPR data is not always a straightforward procedure and may result in considerable ambiguity in the interpretation. This is not surprising given that in the real world, sedimentary units do not fall into discrete categories but rather form a continuum of grain sizes, orientation, packing etc. This makes any imposed rigid classification scheme problematic. However, despite these difficulties and the ambiguities in interpretation, we suggest that radar AEs are useful in some areas as a means of identifying and describing heterogeneities in the subsurface. Once identified in the GPR data, we suggest that radar AEs may be used as hydrogeologic units for use in a groundwater model. This application of radar AEs is supported by the fact that many features of radar AEs correspond to significant hydrogeologic features. The aquifer/aquitard interface, for example, is both a  major hydrogeologic boundary and a radar boundary.  Similarly, major  sedimentary boundaries have been shown above to correspond to radar boundaries and are also associated with changes in hydraulic properties (Fogg, 1986). The sedimentary boundaries which are less likely to appear in the radar record, i.e. those with small changes in porosity, grain composition, etc., are also less likely to have significant hydraulic contrasts. The internal structure of radar AEs is also significant. The orientations of reflectors within radar AEs are related to the orientation of internal sedimentary boundaries, which play a large role in controlling anisotropy in hydraulic properties (Webb and Anderson, 1996).  In addition, the continuity  of reflectors within radar AEs may contain information about the continuity of hydrogeologic features. Radar AEs can therefore be used as a representation of the structure of hydrogeologic units. In order to be used as hydrogeologic units, however, the hydraulic properties within radar AEs must be determined. For this we propose several alternatives. Ideally the hydraulic properties of each type of radar AE can be measured directly using traditional methods such as a pump or slug test. More practically, these properties can be infened either from lithological data obtained from drill holes or from electrical conductivity information obtained from analysis of attenuation in the GPR data. In these cases, the hydraulic conductivity may only be identified as "good" or "poor" based upon empirical  76  relationships between clay content, electrical conductivity and hydraulic conductivity (Urish, 1981; Mazac etal., 1985). It is clear that the hydrogeologic properties inferred from radar AEs are more qualitative than quantitative. The hydrogeologic properties obtained can only be defined as having "good" or "poor" permeability and as having "long" or "short" continuity. Such imprecise data are of benefit, however, to modellers who presently lack information to constrain values assigned to parameters over large unsampled regions of the subsurface.  Further research may enable more accurate determination of  sedimentary and hydrogeologic parameters from the radar response and so provide better interpretations of radar data. One of the critical issues in this interpretation is the 3D nature of hydrogeologic units. Most radar surveys are collected in 2D and therefore may result in an incorrect interpretation. The length of horizontal reflections, for example, is highly dependent upon the relative orientation of the radar line and the strike direction of the sedimentary feature being imaged.  This emphasises the usefulness of  collecting 3D radar whenever possible in order to provide a more accurate interpretation.  CONCLUSION In this study we have attempted to use a classification scheme based on radar architectural element analysis to improve the hydrogeologic model available for the Brookswood aquifer. We acknowledge the inherent difficulties in imposing a rigid classification scheme but believe it can be beneficial as a first step to understanding the larger scale spatial heterogeneity in the subsurface. In this context, ground penetrating radar (GPR) can be a useful technique for providing data for hydrogeologic modelling in the form of defining hydrogeologic units in the subsurface. GPR may be used in this way as an economical tool to provide data which, in concert with well and/or outcrop data, can improve aquifer characterisation. The basic framework can be applied at both the small (contaminant transport modelling ) and large (regional aquifer modelling) scale.  At the scale of most contaminant studies (10's to 100's of  meters) a completely deterministic approach using GPR may be feasible. GPR can be used to image in detail the 3D stratigraphy of the region and can map exactly the location of no-flow boundaries. At a larger, regional scale, lines of radar data crossing the aquifer can pinpoint critical features in the subsurface and identify sedimentary units which can be used as clues to the overall sedimentary  77  environment. It is likely, given the variability in data quality, that GPR will be found most useful small scale problems at selected sites.  78  C H A P T E R 5: Summary and Conclusions  In this thesis I have developed three different techniques for using ground penetrating radar (GPR) for characterising shallow unconfined aquifers. I present an overall framework for interpreting GPR as radar architectural elements (AEs) which can then be used as basic building blocks in a hydrogeologic model. In order to extract information about the internal properties of these elements, I have developed two new techniques for analysing radar data: a geostatistical analysis to determine correlation structure of radar reflections; and a decay rate analysis to estimate the electrical conductivity of the sediments being imaged in the radar survey. In Chapter 2,1 presented a technique for using geostatistics to determine the correlation structure of reflections in radar surveys. The technique involves the calculation and modelling of semivariograms using the amplitudes contained in the radar data. The main assumption in this analysis is that the correlation structure of radar reflections is related to the correlation structure of the sediments being imaged. In order to test this assumption a cliff face calibration study was conducted. In this study a photograph of a cliff face was analysed using geostatistics to determine the correlation structure of the sedimentary beds. A GPR survey was collected along the cliff to image the sediments in the photograph. The semivariogram obtained from the photograph showed excellent agreement with the semivariogram from the radar data, in determining both the maximum correlation direction and the correlation length. We conclude, based on the examples in this study, that radar data can be used in this way to quantify the spatial heterogeneity of the subsurface. Chapter 3 introduced a technique for estimating the electrical conductivity of the subsurface through analysis of the attenuation of radar data. In fresh water environments electrical conductivity is closely related to the clay-content of the earth which is in turn related to hydraulic conductivity. This property is therefore useful both as an indicator of the lithology of sediments and as a means of putting bounds on the hydraulic conductivity of the sediments.  Electrical conductivity is a major factor  contributing to the attenuation of electromagnetic waves. I have developed a method of accounting for other attenuating factors in a radar survey, to isolate those effects directly caused by electrical conductivity. This technique was tested using several GPR data sets. From these results we conclude that reasonable estimates of electrical conductivity can be obtained from the analysis of GPR attenuation. The results obtained using this method are most accurate in more resistive regions where tighter bounds can be placed on the electrical conductivity estimates. Using analysis of GPR attenuation, conductivity data can be obtained rapidly and efficiently over a large area. 79  In Chapter 4, I discuss the use of GPR for determining hydraulic boundaries in the subsurface, including the interface between the saturated and unsaturated zone, and hydrogeologic boundaries such as the aquifer/aquitard interface.  In this analysis I use a radar classification scheme called radar  architectural element analysis which divides the radar data into packages of reflections with similar characteristics. This technique was investigated using a series of cliff face studies which allow direct comparison between radar data and the sediments being imaged. A specific set of experiments were conducted to develop techniques for identifying the reflection from the saturated zone boundary in the radar data. This study identified SEC gaining of the data as the best ways of differentiating between reflections from geologic boundaries and saturated boundaries. We applied the techniques developed in the cliff face and saturated zone studies to a GPR data set collected over the Brookswood Aquifer, a shallow, unconfined aquifer in south-western British Columbia. The top of the saturated zone was identified in most of the Brookswood data. However, we encountered mixed success in the identification of radar.AEs in the data.  At some sites, where  penetration was good, clear patterns of reflections could be identified. In other regions, however, the radar penetration was poor and the pattern of reflections was difficult to distinguish. Another critical issue is that the orientation of the radar line with respect to the underlying geologic structure can have a significant effect on the pattern of reflections in the GPR data. We therefore concluded that GPR is best suited to smaller scale problems at selected sites which allow good radar penetration. In this way the capability of GPR for high resolution 3D imaging of the subsurface can be used to best advantage.  Future Work The work presented in this thesis has laid the ground work for numerous follow-up studies. The use of radar architectural element analysis needs to be tested at other sites, particularly those in different depositional environments. The aspect of this work which is in most need of further investigation is the nature of correlation between radar reflections and sedimentary or hydrogeologic boundaries. Extensive cliff face studies are essential for a better understanding of this relationship, since they allow direct comparison between the radar data and the sediments being imaged. One of the drawbacks to cliff face studies is that the sediments are always unsaturated, whereas in most applications of this technique, the sediments being imaged will be fully saturated. As water plays a vital role in the determination of dielectric constant it is therefore a critical parameter in radar calibration studies. One solution may lie in conducting a series of laboratory experiments to determine the dielectric properties of saturated 80  sediments. There are serious problems inherent in this solution also, however, mainly regarding issues of scale. The technique of geostatistical analysis of radar data must undergo further examination. Forward modelling studies are suggested to investigate the effects of antenna frequency on the results as well as to answer some remaining questions as to the robustness of this technique. The decay rate analysis technique would also benefit from further investigations through forward modelling. In particular the effect of energy splitting at sedimentary boundaries needs to be assessed. Again cliff face studies would be advisable, perhaps making use of in situ sampling techniques to determine dielectric contrasts across sedimentary boundaries. To summarise, this thesis presents several new approaches to the interpretation of radar data. Using these techniques data can be provided to enhance aquifer characterisation and lead to the creation of more accurate models for both groundwater and contaminant flow problems.  81  References Anderson, M.P., 1989, Hydrogeologic facies models to delineate large-scale spatial trends in glacial and glaciofluvial sediments: Geological Society of America Bulletin, 101, 501-511. Annan, A.P., and Chua, L.T., 1992, Ground penetrating radar performance predictions, in Pilon, J. Ed., Ground penetrating radar, Paper 90-4: Geological Survey of Canada, 5-13. Annan, A.P., and Cosway, S.W., 1992, Ground Penetrating Radar Survey Design: Symposium on the Application of Geophysics to Engineering and Environmental Problems, Proceedings, 329-352. Annan, A.P., and Davis, J.L., 1977, Radar range analysis for geologic materials: GSC Paper 77-1B, Report of Activities, 117-124. Annan, A.P., and Davis, J.L., 1989, Ground Penetrating Radar for high resolution mapping of soil and rock stratigraphy: Geophysical Prospecting, 37, 531-551. Atkinson, P.M., 1993, The effect of spatial resolution on the experimental variogram of airborne MSS imagery: International Journal of Remote Sensing, 14, 1005-1011. Baker, P.L., 1991, Response of ground penetrating radar to bounding surfaces and lithofacies variations in sand barrier sequences: Exploration Geophysics, 22, 19-22. Best, M.E., and Todd, B.J., 1996, Fraser Valley Hydrogeology Poject: time-domain EM surveys August 4 to August 18, 1995: supplement to Geological Survey of Canada Open File Report 30950. Beres, M., and Haeni, F.P., 1991, Application of ground-penetrating-radar methods in hydrogeologic studies: Ground Water, 29, 375-386. Campanella, R.G., Davies, M., Boyd, T., Everard, 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 Report 2940. Chu, J., 1993, X G A M User's Manual: Stanford Centre for Reservoir Forecasting. Davis, J.L., and Annan, A.P., 1989, Ground penetrating radar for high-resolution mapping of soil and rock stratigraphy: Geophysical Prospecting, 37, 5, 531-551. Deutsch, C.V., and Journel, A.G., 1992, GSLIB Geostatistical Software Library and User's Guide: Oxford University Press. Drury, S.A., 1987, Image Interpretation in Geology: Allen and Unwin. Eyles, N., and Eyles, C.H., 1992, Glacial Depositional Systems, in Walker, R. G. and James, N. P., Ed., Facies Models: Response to sea level change: Geological Association of Canada, 73-100. Fogg, G.E., 1986, Groundwater flow and sand body interconnectedness in a thick multiple-aquifer system: Water Resources Research, 22, 679-694. Gawthorpe, R.L., Collier, R.E., Alexander, J., Bridge, J.S., and Leeder, M.R., 1993, Ground penetrating radar: application to sandbody geometry and heterogeneity studies, in North, CP. and Prosser, D.J., Ed., Characterization of Fluvial and Aeolian Reservoirs: Geological Society Special Publications, 73, 421 432.  82  Greaves, R.J., Lesmes, D.P., Lee, J.M., and Toksoz, M.N., Velocity variations and water content estimated from multi-offset, ground-penetrating radar: Geophysics, 61, 683-695. Greenhouse, J.P., Barker, J.F., Cosgrave, T.M., and Davis, J.L., 1987, Shallow stratigraphic reflections from ground penetrating radar: Department of Energy, Mines and Resources, Research Agreement 232. Huggenberger, P., Meier, E., and Beres, M., 1994, Three-dimensional geometry of fluvial gravel deposits from GPR reflection patterns; a comparison of results of three different antenna frequencies: Fifth International Conference on Ground Penetrating Radar, Proceedings, 805-816. Isaaks, E.H., and Srivastava, R . M . , 1989, An Introduction to Applied Geostatistics: Oxford University Press. Jol, H.M., and Smith, D.G., 1992, GPR results used to infer depositional processes of coastal spits in large lakes: Fourth International Conference on Ground Penetrating Radar, Proceedings, 169-177. Jol, H.M., and Smith, D.G., 1991, Ground penetrating radar of northern lacustrine deltas: Canadian Journal of Earth Science, 28, 1939-1947. Journel, A.G., and Froidevaux, R., 1982, Anisotropic hole-effect modeling: Mathematical Geology, 14, 217-239. Keller, G.V., 1989, Electrical conductivity and dielectric constant of minerals and dry rocks, in CRC practical handbook of physical properties of rocks and minerals, ed. R.S.Carmichael: CRC Press, Inc., 361-369. Kelly, W.E., 1977, Geoelectric sounding for estimating aquifer hydraulic conductivity: Ground Water, 15, 420-424. Kosinski, W.K., and Kelly, W.E., 1981, Geoelectric soundings for predicting aquifer properties: Ground Water, 19, 163-171. Lacaze, B., Rambal, S., and T. Winkel, 1994, Identifying spatial patterns of Mediterranean landscapes from geostatistical analysis of remotely-sensed data: International Journal of Remote Sensing, 15, 24372450. Livelybrooks, D., and Fullagar, P.K., 1994, FDTD2D+, A finite difference, time-domain radar modelling program for two dimensional structures: Symposium on the Application of Geophysics to Engineering and Environmental Problems, Proceedings, 87-100. Mazac, O., Cislerova, M., Kelly, W.E., Landa, I., and Venhodova, D., 1990, Determination of hydraulic conductivities by surface geoelectric methods, in Ward, S., Ed., Geotechnical and Environmental Geophysics: Society of Exploration Geophysics, 2, 125. Mazac, O., Cislerova, M., and Vogal, T., 1988, Application of geophysical methods in describing spatial variability of saturated hydraulic conductivity in the zone of aeration: Journal of Hydrology, 103, 117126. Mazac, O., Kelly, W.E., and Landa, I., 1985, A hydrogeophysical model for relations between electrical and hydraulic properties of aquifers: Journal of Hydrology, 79, 1-19. Miall, A.D., 1985, Architectural element analysis: A new method of facies analysis applied to fluvial deposits: Earth-Science Reviews, 22, 261-308.  83  Oldenburg, D.W., and Li, Y., 1994, Inversion of polarization data: Geophysics, 59, 1327-1341. Olhoeft, G.R., 1994, Quantitative statistical description of subsurface heterogeneities with ground penetrating radar at Bemidji, Minnesota: United States Geological Survey Water Investigations Report, 650-653. Olhoeft, G.R., 1974, The electrical properties of permafrost: Ph.D. Thesis, the University of Toronto. Olsen, H., Andreasen, F., 1995, Sedimentology and ground-penetrating radar characteristics of a Pleistocene sandur deposit: Sedimentary Geology, 99, 1-15. Pratt, B.R., Miall, A.D., 1993, Anatomy of a bioclastic grainstone megashoal (Middle Silurian, southern Ontario) revealed by ground penetrating radar: Geology, 21, 223-226. Rea, J and Knight, R, 1996, Geostatistical analysis of ground penetrating radar data: a means of describing spatial variation in the subsurface, submitted to Water Resources Research. Rea, J.M., and Knight, R.J., 1995, The use of ground penetrating radar for aquifer characterization: an example from south-western British Columbia: Symposium on the Application of Geophysics to Engineering and Environmental Problems, Proceedings, 279-288. Rea, J.M., Knight, R.J., and Ricketts, B.D., 1994, Ground penetrating radar, Brookswood aquifer, Lower Fraser Valley, British Columbia: Geological Survey of Canada Open File Report 2821. Ricketts, B., and Jackson, S., 1994, An overview of the Vancouver-Fraser Valley hydrogeology project, southern B. C : Geological Survey of Canada Current Research, 94-A, 201-206. Ritzi, R.W., Jayne, D.F., Zahradnik, A.J., Field, A.A., and Fogg, G.E., 1994, Geostatistical modeling of heterogeneity in glaciofluvial, buried-valley aquifers: Ground Water, 32, 666-674. Stephens, M., 1994, Architectural element analysis within the Kayenta Formation (Lower Jurassic) using ground-penetrating radar and sedimentological profiling, southwestern Colorado: Sedimentary Geology, 90,179-211. Sudicky, E.A., 1986, A natural gradient experiment on solute transport in a sand aquifer: spatial variability of hydraulic conductivity and its role in the dispersion process: Water Resources Research, 22(13), 2069-2082. Tribe, S., 1994, Geometry and sedimentology of a Quaternary delta in Kirkpatrick pit, Fraser Valley, British Columbia: B. Sc. thesis, University of British Columbia. Urish, D.W., 1981, Electrical resistivity-hydraulic conductivity relationships in glacial outwash aquifers: Water Resources Research, 17, 1401-1408. van Overmeeren, R.A., 1994, Georadar for Hydrogeology: First Break, 12, 401-408. Vellidis, G., Smith, M.C., and Thomas, D.L., 1990, Detecting wetting front movement in a sandy soil with ground penetrating radar: Transactions of the ASAE, 33, 186. Ward, S.H., and Hohmann, G.W., 1987, Electromagnetic theory for geophysical applications, in Nabighian, M.N., Ed., Electromagnetic Methods in Applied Geophysics: Society of Exploration Geophysics, 1, 131. 84  Webb, E.K., and Anderson, M.P., 1996, Simulation of preferential flow in three-dimensional, heterogeneous conductivity field with realistic internal architecture: Water Resources Research, 32, 533545. Wharton, R.P., Hazen, G.A., Rau, R.N., and Best, D.L., 1980, Electromagnetic propagation logging: advances in technique and interpretation: Soc. Pet. Eng. of AIME, Paper 9267. Woodcock, C.E., Strahler, A.H., and D.L.B. Jupp, 1988, The use of variograms in remote sensing: II. Real digital images: Remote Sensing of the Environment, 25, pp. 349-379. Yilmaz, O., 1987, Seismic Data Processing: Society of Exploration Geophysics. Young, R., 1995, 3D GPR imaging of a sand-gravel/clay boundary at Hill Air Force Base, Utah: Symposium on the Application of Geophysics to Engineering and Environmental Problems, Proceedings, 905-912.  85  APPENDIX Map showing location of radar lines and photographs discussed in this thesis. Thesis figure numbers corresponding to radar sections and photographs are shown in italics. Detailed maps of survey lines are shown on pages 87 and 88. 235'  49'  236*  237*  H  48-  86  238'  239  40A Avenue  ^s^j^OOm  Figure 4.20;  36th Avenue  32nd Avenue  800m  28th Avenue Figure 4.20b  00 m  24th Avenue  Stokes Pit  20th Avenue  87  24th Avenue  Stokes Pit Site 2 of water table survey (Figures 4.3,4.4a, 4.4b, 4.5b, 4.6, and 4.8)  / 300rJ  Figures 3.6 and 4.18  I J J Q m ^ Z — 3D GPR survey (Figures 2.17, 4.19) ^00m  20th Avenue  24th Avenue  J^7oTT 400 m  Aquifer pinchout (Figures 3.2a, 3.2b, 3.5, 4.16 and 4.17) 20th Avenue  1 100 m  470 m  ' Figure 4.15  

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-0052320/manifest

Comment

Related Items