"Science, Faculty of"@en . "Earth, Ocean and Atmospheric Sciences, Department of"@en . "DSpace"@en . "UBCV"@en . "Rea, Jane Mary Anastasia"@en . "2009-04-02T23:58:01Z"@en . "1996"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "The goal of this thesis is to develop techniques for using ground penetrating radar (GPR) to\r\ncharacterize aquifers for groundwater modelling. GPR is a shallow geophysical technique which can be\r\nused to image up to 30 m into the subsurface. This technique is sensitive to many of the same parameters\r\nthat affect the hydraulic properties of geologic material. In addition it is a rapid, relatively inexpensive\r\nmeans of obtaining high resolution images of the subsurface. For these reasons GPR is well suited for\r\nproviding the information needed to constrain ground water models.\r\nI have developed three different techniques for analysing and interpreting GPR data for ground\r\nwater modelling. The first is a method of obtaining the correlation structure of the subsurface through\r\ngeostatistical analysis of GPR data. The second technique involves analysing the attenuation of radar\r\ndata to estimate the electrical conductivity of the geologic material being imaged. The final technique is a\r\nmethod of dividing the radar data into separate units, called radar architectural elements, which can be\r\nused as building blocks for a ground water model. This technique has been tested by applying it to a\r\nshallow unconfiried aquifer in south-western British Columbia. The research presented in this thesis\r\nleads to a new integrated approach to GPR interpretation in which the same data set can be used to\r\nprovide different types of information to constrain ground water models."@en . "https://circle.library.ubc.ca/rest/handle/2429/6749?expand=metadata"@en . "12195642 bytes"@en . "application/pdf"@en . "G R O U N D P E N E T R A T I N G 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 FULFILLMENT O F T H E REQUIREMENTS FOR T H E D E G R E E OF 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 \u00C2\u00A9 J a n e Mary Anastasia Rea, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department 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 data to estimate 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 1 List of Figures iv List of Tables v i i Acknowledgements v i i i Chapter 1: Introduction 1 Chapter 2: Geostatistical analysis of ground penetrating radar data: a means of describing spatial variation in the subsurface 3 Introduction 3 Geostatistical analysis of data 4 Description of the gravel pit study 5 Radar length scales in architectural elements 14 Variograms from 3D elements 19 Conclusions 21 Chapter 3: Obtaining Estimates of Electrical Conductivity from the Attenuation of Ground Penetrating Radar Data 23 Introduction 23 Analysis of radar attenuation 24 Field experiment: Comparison of GPR and DC resistivity results 29 GPR results over two layers 39 Conclusion 41 Chapter 4: Determination of Hydraulic Boundaries using Ground Penetrating Radar 43 Introduction 43 The interpretation of GPR data 44 Study to investigate the top of the saturated zone 48 Cliff face study 57 The Brookswood aquifer study 64 Discussion 76 Conclusion 78 Chapter 5: Summary and Conclusions 79 References 82 Appendix 86 i i i 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. i v 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 data from Site 2 with SEC gain applied. Figure 4.5: a) GPR data from Site 1 with no gain applied, b) GPR data from Site 2 with no gain applied. Figure 4.6: a) Instantaneous phase plot of GPR data from Site 2. b) Instantaneous frequency plot of GPR data from Site 2. Figure 4.7: a) Instantaneous phase plot of GPR data from Site 1. b) Instantaneous frequency plot of GPR data from Site 1. 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. 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 192nd Avenue. v i 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 v i i 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. v i i i 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 the frequency range 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 high-resolution, 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): 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 corresponding to a 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 (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 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 2.1a: 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: y(h) = c-h (h 1.5\u00E2\u0080\u0094.5 -a Va = c , i f h>a i f h < a (2) Exponential model: t \ [ f 3h^ y(h) = c- 1-exp - \u00E2\u0080\u0094 v a J (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 \u00E2\u0080\u00A2 {l - [(e~h/&) \u00E2\u0080\u00A2 cos{2nh / X)]} (4) 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 x2, where 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 THE GRAVEL 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 \u00E2\u0080\u00A2 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\u00C2\u00BBothing 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. The final result, 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 XGAM 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 0.14 0.12 - 0.1 0 08 0 06 0 04 0 02 0 semivariogram values - semivariogram model 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 calculated from the 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 2.10: 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 lengths from radar 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 CM ii 10 Distance (m) 20 30 40 Q 5 Figure 2.11: Radar survey showing the dipping radar element, dipl, which is associated with dipping foreset beds. a) b) 0.5 0.4 0.3 0.2 0.1 0 Semivariogram of dipl semivariogram values -semivariogram model Lag (m) Semivariogram of dip2 semivariogram values semivariogram model Lag (m) Semivariogram of dip3 0.9-- \u00E2\u0080\u00A2 semivariogram values 0 8 - - 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 sub-horizontal 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 of short2 0.9 0.6 0.7 0.6 -0.5 0.4 0.3 0.2 0.1 0 c) semivariogram values -semivariogram model 1 2 3 4 Lag (m) Semivariogram of short3 semivariogram values \u00E2\u0080\u00A2 semivariogram model 1 2 3 4 Lag (m) 1 0.9 0.8 0.7 0.6 -0.5 0.4 0.3 0.2 0.1 0 semivariogram values . semivariogram model 1 2 3 4 5 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 XGAM 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 PIT 10 20 Distance (m) 30 40 50 60 70 . . 1 rt*)..*.^. . . \u00E2\u0080\u00A2 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) c) Semivariogram for longl Semivariogram for long3 Lag (m) b) Semivariogram for long2 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 FROM 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 XGAM 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 \u00E2\u0080\u00A2 semivariogram values 1 1 1 semivariogram model 0 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. In the glacio-deltaic environment, we obtained 21 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 non-unique, 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 S , which undergoes a power loss associated with the antenna efficiency, <|Tx. The power radiated downwards into the ground is further affected by the antenna gain, G T x , which is a measure of the directivity of the antenna. The power reaching a reflector has been further reduced due to the intrinsic attenuation of the medium through which the wave has travelled. This decay occurs exponentially as e~2 a L, 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/(47tL2). 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): \u00C2\u00A3 T X ' ^ R X G T x - G R x A-g-x-e\" 4 PRX The power P ^ can be given as 16TI2L4 P< s (3.1) where is the voltage at the receiver antenna is the impedance of the receiver antenna. The effective antenna area is defined as GR X A=V 2 G R X /4 n f2 where v is the propagation velocity f is the frequency of the signal. 25 P O W E R F R O M S O U R C E Ps RADIATED P O W E R P. = ^Tx-Ps P O W E R RADIATED D O W N W A R D P 2 = G T x - P , P O W E R DENSITY R E A C H I N G R E F L E C T O R - 2 a L 4TCI/ P O W E R A T E L E C T F PRX=^ R E C E I V E R RONICS = Rx 'P7 R E C E I V E D P O W E R P7 = AGRxP6 P O W E R DENSITY A T R E C E I V E R \u00E2\u0080\u009E - 2 a L ^ 6 = ^ - T P 5 4itL2 P O W E R D I R E C T E D B A C K T O R E C E I V E R P 5 = g P 4 Where: \u00C2\u00A3 T X = transmitter antenna efficiency ^ .RX = receiver antenna efficiency a = attenuation coefficient of medium X = scattering cross sectional area A = receiver antenna effective area P O W E R S C A T T E R E D BY R E F L E C T O R P 4 =X g = backscatter gain of reflector G T x = transmitter antenna gain G,^ = receiver antenna gain L = depth 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 L l 8 + (27xf) 2L leJl + [ ^ r N 2 (3.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=e0c2/v2, where c is the speed of light in a vacuum, 3.00xl08 m/s, s0 is the permittivity of free space, 8.85xl0\"12 and v is the velocity of the electromagnetic 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 s, Z^, cjTx, %RX> GTX a n 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= 1 17.38 f -10 /og( ) - 40 \u00E2\u0080\u00A2 log L + 10 log v f 2 y + I0log(g) + \0log(x) 1 _l P S Z R X ^ T X ^ R X G T X G R X 1.738 *V 64TI3 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 the frequency are 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 Vm(t). As we are interested in the total power of the signal at an instant in time, we convert the measured voltages to instantaneous amplitudes which measure reflectivity strength (Yilmaz, 1987). Instantaneous amplitude is given by the expression r 2 2 \"j 1/2 vRx (t)= V(t) + vm'(t) (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 Vm(t). 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 velocity, using L=v(t/2), where t is the recorded 2-way travel time. Scattering cross-sectional area, % 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 TCL2 for a smooth 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*2 are the complex dielectric permittivities of the upper and lower layers respectively. We consider only the high frequency case 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 R2, and the transmitted power to (1-R)2. If there 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 R2(l-R)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 \u00E2\u0080\u00A2 n \u00E2\u0080\u00A2 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 obtained from CMP surveys, through the approximation E ^ C V V 2 , (3.6) where c is the speed of light in a vacuum, 3.00 xlO8 m/s and e0 is the permittivity of free space, 8.85xl0\"12 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 20th Ave. is just east of 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 \u00E2\u0080\u00A2 r* / i i 192nd St. i I ' 0 1 2 3 4 5 Scale (kilometers) 192nd St. 36th Ave. i i 24th Ave. 28th Ave. H 03 ? 1 IO & \ EC ro o o \l6th Ave. th St \ \ ; \ \ \u00E2\u0080\u00A2 / I \u00E2\u0080\u00A2 x I. 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 208th Street south of 24th Avenue. We 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 208th Street at the northern end of line 24208 show that the 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 208th Street. 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 Distance (m) a) 0 100 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 20th Avenue between 204th and 208th Streets in 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 208th Street. A portion of the 500 m line, collected 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 20th Ave. is just east of 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 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. The remaining data are identified as a single element of horizontal reflections. 70 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 192nd Street and 24th Avenue. The Stokes pit has been active for 40 years and at present the top 6 meters have been excavated. The Stokes pit lies completely within Halstead's Unit C, which extends down to a 111 G !a \"2 4 5 4 3 s * rt e T 3 cS rt U S 3 \u00C2\u00A3 s .s u o a w g 6 \" 2 * P is fr ^ (3 D S 3 \u00E2\u0080\u00A2 \u00C2\u00A3 * 4 3 C 13 w \u00C2\u00AB> .2 ^ 2 o 3 +J c oo S s S 3 rt 4 3 1) C I U 6 S 3 a a. T 3 rt o o Z 3 ^ 3 8 - S rt JT^ & K rt qjj 4 3 M a ~ -2 ^ \u00C2\u00A3 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 n d and 3 r d dipping elements is interpreted as a channel element. The abrupt end to 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 192nd Street in Langley, British Columbia. In 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, XGAM 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, 2437-2450. 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, 117-126. 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, 533-545. 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. 2 3 5 ' 4 9 ' H 236* 237* 2 3 8 ' 239 48-86 Figure 4.20; 40A Avenue ^ s ^ j ^ O O m Figure 4.20b 24th Avenue 36th Avenue 32nd Avenue 800m 28th Avenue 00 m 20th Avenue Stokes Pit 87 24th Avenue 300rJ 20th Avenue Stokes Pit Site 2 of water table survey /(Figures 4.3,4.4a, 4.4b, 4.5b, 4.6, and 4.8) Figures 3.6 and 4.18 I J J Q m ^ Z \u00E2\u0080\u0094 3D GPR survey (Figures 2.17, 4.19) ^00m 24th Avenue 400 m J^ 7oTT ' Figure 4.15 Aquifer pinchout (Figures 3.2a, 3.2b, 3.5, 4.16 and 4.17) 20th Avenue 1 100 m 470 m "@en . "Thesis/Dissertation"@en . "1997-05"@en . "10.14288/1.0052320"@en . "eng"@en . "Earth and Ocean Sciences"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en . "Graduate"@en . "Ground penetrating radar applications in hydrology"@en . "Text"@en . "http://hdl.handle.net/2429/6749"@en .