Published online 25 February 2008
Published in Vadose Zone J 7:140-159 (2008)
DOI: 10.2136/vzj2006.0142
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: GROUND PENETRATING RADAR IN HYDROGEOPHYSICS
GPR Attenuation and Scattering in a Mature Hydrocarbon Spill: A Modeling Study
Nigel J. Cassidy*
Applied and Environmental Geophysics Group, School of Physical and Geographical Sciences, Keele Univ., Staffordshire, ST5 5BG UK
* Corresponding author (n.j.cassidy{at}esci.keele.ac.uk).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 27 September 2006.
 |
ABSTRACT
|
|---|
Subsurface contamination by light nonaqueous-phase liquids (LNAPLs) is a pressing environmental issue with many industrial sites having some degree of near-surface pollution. Noninvasive, geophysical investigation methods, such as ground penetrating radar (GPR), have become increasingly popular, but their use has tended to be restricted to plume "mapping" and recent contamination events only. Mature LNAPL contamination, where the hydrocarbons undergo significant alteration and biodegradation, complicates the electrophysical properties of the subsurface, making GPR data interpretation more difficult. This is particularly true in shallow, coastal environments where LNAPLs form a disaggregated, smeared zone associated with the seasonally changing groundwaters. Therefore, any method that can improve data interpretation has a significant practical benefit in these complex environments. By using dielectric material property analysis, empirical mixing model evaluation, and three-dimensional, finite-difference time-domain (FDTD) modeling, it has been possible to investigate the nature and spectral content of GPR signal attenuation and scattering within the vadose or smeared zone of a mature LNAPL-contaminated site. Using FDTD models to simulate real subsurface conditions, a comparison between the GPR response of "clean" and contaminated environments has been made where the subsurface materials contain different LNAPL–porewater saturations and contaminant geometries. The results show that uniform mixtures, disaggregated scattering mixtures, and stratified pools of contaminated material all exhibit different, yet characterizable, temporal and spectral GPR responses and that the recorded data can provide new insights into the mode of contaminant distribution and saturation as long as the incident field of the GPR system can be determined accurately.
Abbreviations: CRIM, complex refractive index model EM, electromagnetic FDTD, finite difference time domain GPR, ground penetrating radar LNAPL, light nonaqueous phase liquid
 |
INTRODUCTION
|
|---|
Ground penetrating radar has become an increasingly popular noninvasive tool for the investigation and characterization of groundwater contamination by LNAPL contaminants such as petroleum fuels, solvents, and other mobile hydrocarbons (e.g., Marcak and Golebiowski, 2006; Lopes de Castro and Branco, 2003; Bradford, 2003; Burton et al., 2004; Campbell et al., 1996; Benson, 1995; Daniels et al., 1995; DeRyck et al., 1993; Olhoeft, 1992). Despite these successes, there is still some debate on the interpretation of the observed GPR results and how the data can be used to determine the physical characteristics of the contaminating materials (such as obtaining reliable estimates of the LNAPL saturation index). Investigations have shown that recent LNAPL spills can produce increased amplitude reflections, or "bright spots" associated with the sharp interfaces created at the boundary between the high-permittivity, highly attenuating, saturated groundwater and the low-permittivity, low-loss, LNAPL-saturated materials (e.g., Campbell et al., 1996; Benson, 1995). Conversely, it is also possible to observe "dim spots" in the reflection response from the groundwater interface where isolated or thin zones of free-phase LNAPLs are preferentially attenuating or scattering the GPR signal.
Mature spills, however, undergo a significant process of alteration and biodegradation that actually increases the bulk conductivity (and therefore loss) of the LNAPL-contaminated zone. This is particularly true where the LNAPL is in a disaggregated, free-phase form and associated with saturated groundwater under anaerobic conditions (Cozzarelli et al., 2001, 1994; Baedecker and Cozzarelli, 1994). In situ geophysical investigations and laboratory experiments have supported this hypothesis (e.g., Abdel Aal et al., 2006; Atekwana et al., 2004a,b,c, 2002, 2000; Cassidy et al., 2001; Bradford, 2003; Osella et al., 2002; Kim et al., 2000; Sauck et al., 1998) and, in general, show that the presence of microorganisms in the vadose and shallow saturated zones generate elevated levels of biosurfactants and organic and carbonic acids as a byproduct of the physical and chemical alteration of the hydrocarbon contaminants. These mobile acids attack the mineral constituents of the matrix materials (e.g., feldspars, quartz, etc.) that, in turn, release dissolved ions into the pore waters. The combination of these free ions, organic and carbonic acids, and dispersed bacteria increases the electrical conductivity of the interstitial pore waters and, therefore, the bulk attenuation or loss characteristics of the materials.
The complexity of the biodegradation process, plus the natural variation present in the near-surface environment, means that, in practice, a range of GPR responses will be evident at any LNAPL-contaminated site. It is common to find highly attenuating regions (which produce "shadow zones" in the GPR data) being associated with "bright spot" reflections and signal enhancement. This level of complexity is more evident in recent studies (e.g., Lopes de Castro and Branco, 2003) as our understanding of the physical evolution of a LNAPL spills moves away from a simple, four-component vapor-, residual-, free-, and dissolved-phase model to a more complex, mixed-component, multiphase "smeared zone" model (Lowe et al., 1999). In this case, the mobile LNAPL products have a tendency to preferentially propagate along pathways of increased permeability as the water table depth changes seasonally. This process results in the LNAPLs being "smeared" into a zone of disaggregated units that, depending on the contaminant's viscosity and localized capillary pressure, may become immobile under unconfined conditions. This complexity has led to a focus on contaminant plume identification only, with less effort being applied to developing more advanced characterization techniques. Studies by Marcak and Golebiowski (2006), Lopes de Castro and Branco (2003), Carcione et al. (2006) and Carcione and Seriani (2000) have shown, however, that it is possible to obtain important information on the nature and properties of the contaminants (e.g., material property information, volumetric LNAPL saturation, and contaminant distribution and migration pathways, etc.) through the judicial use of advanced data processing methods such as GPR signal and attribute analysis.
A recent GPR signal attenuation and dielectric property study by Cassidy (2007, 2006, 2004b) has identified specific attenuation characteristics associated with a mature, near-surface (<2-m) benzene–toluene–styrene–ethylbenzene LNAPL contamination problem at a coastal hydrocarbon storage and processing facility. As part of a larger research program, the work at this study site is intended to provide a practical evaluation of attenuation-based GPR analysis and processing methods and has revealed that the most significant, identifiable attenuation (i.e., where the attenuation can be wholly attributed to material property losses) is associated with contaminated pore waters in the "central", or major, part of the smeared zone. Figure 1
shows a summary model of the physical, geologic, and GPR signal attenuation characteristics at the site, as determined from laboratory dielectric measurements and GPR signal analysis. Key to this research was the identification of a dominant static or ionic conductivity loss associated with the high-permittivity, contaminated pore waters and enhanced levels of relative signal attenuation (>2–6 dB/m) in the smeared zone only. In contrast, the pure free-phase LNAPLs and the aeolian sands sampled from within the vadose zone were low permittivity and low loss, and showed relatively little signal attenuation. Within the groundwater-saturated zone (beneath the notional water table) the relative attenuation between the "clean" and contaminated areas was not strong enough to be uniquely identifiable using signal-processing techniques alone. The laboratory results indicated, however, that the contaminated waters were more attenuating than the natural waters. Within the smeared zone, where the attenuation was highest, there was a slight but observable frequency dependence to the attenuation spectrum, with the lower frequency data (<300 MHz) showing the greatest degree of loss. This suggests that frequency-related processing methods (such as instantaneous-frequency attribute analysis) may provide additional material property information that can be extracted from the frequency spectrum of the GPR data. Unfortunately, the frequency dependence of the recorded signal spectrum cannot be directly related to the material property losses per se. For any particular survey, subsurface attenuation factors (such as antennae-to-ground coupling loss, scattering losses, geometrical spreading losses, etc.) combine with the less dominant system-related factors (e.g., antennae losses, system efficiency) to give an overall level of attenuation that is site specific (Daniels, 2004). For sites where the lateral variation in geology is gradual, or at scales much larger than the dominant GPR wavelength, the geometrical spreading losses at each data collection point can be considered equivalent. Consequently, any variation in signal amplitude or frequency will be predominantly related to the attenuation or scattering losses associated with the subsurface materials only. What attribute or signal analysis cannot reveal is whether the observed attenuation (and its frequency dependence) is related to the signal scattering, material property attenuation, or a combination of both. As mature LNAPL spills are likely to be highly disaggregated, either as distinct "pools" or "blobs", it is likely that the scattering within the vadose and smeared zones will play a significant role in determining the strength and spectrum of the recorded GPR signal. Therefore, the aim of this work was to utilize three-dimensional, finite-difference time-domain (FDTD) numerical GPR modeling in combination with dielectric property mixing models to evaluate the contribution of both scattering and material attenuation mechanisms to the overall signal loss observed at the study site.

View larger version (38K):
[in this window]
[in a new window]
|
FIG. 1. Summary model of the physical, geologic, and ground penetrating radar (GPR) signal attenuation characteristics at the study site consisting of a near-surface "aged" light nonaqueous-phase liquid (LNAPL) that is smeared across the vadose and saturated zones due to seasonal variations in water table depth.
|
|
 |
Site Overview and Data Collection
|
|---|
The study site is a demolished hydrocarbon storage facility with a 40-yr legacy of LNAPL leakage and surface contamination, resulting in the presence of residual, free-phase and dissolved-phase contaminants in the top 1 to 3 m of the subsurface (predominantly styrene, benzene, ethylbenzene, toluene, and diethylbenzene). The subsurface geology consists of laterally uniform, clean, well-sorted, medium-to-coarse-grained aeolian beach sands with a shallow covering of mixed sand, hardcore, and concrete "made ground". The sands have an average porosity of 42% (maximum 45%) and consist of clean quartz grains (87%), lithics (4%), shelly fragments (4%), vegetable matter (2%), and a 3% clay fraction. Sampled water contents range from 3 to 5% (i.e., dry conditions), while groundwater levels were typically within 1 to 2 m of the surface and exhibit strong seasonal variations. At the time of the survey, the whole region had undergone a prolonged period of very low rainfall and the groundwater table was at a level of about
1.5 m below ground level. Known contaminant hotspots produce localized, smeared zones of mixed free- and dissolved-phase materials having layer thicknesses of between 0.1 and 0.7 m at a depth of approximately 1.2 to 1.5 m (sampled from closed boreholes and open trial pits). Trial pit investigations showed that the free-phase LNAPLs tended to be associated with the higher permeability pathways created by natural geologic variations (e.g., dune forests and cross bedding), resulting in disaggregated pools, lenses, and blobs of contaminant that ranged in spatial extent from centimeters to meters. These observed features have led to the smeared–disaggregated model of Fig. 1, where the mixed LNAPL and contaminated pore waters range in saturation from about 10% (of total pore volume) in the vadose zone to full saturation within the center of the smeared zone. In comparison, the "clean" areas of the site show little (or only minor) contamination and invasive sampling has indicated that there is an observable, yet narrow (
0.05–0.10 m) capillary zone associated with the water table horizon.
Ground penetrating radar surveys were undertaken in both the "clean" and contaminated areas of the site using a Sensors and Software PE1000 with 225- and 450-MHz shielded antennae (Sensors and Software, Mississauga, ON, Canada). Common-offset, co-planar, broadside mode reflection surveys were combined with common midpoint (CMP) velocity profiles with all data collected in step mode with 0.05-m spatial increments, 0.2-ns temporal sampling, and 16 pulse stacks per sample point. Processing steps included "de-wowing" to remove low-frequency signal bias, time zero adjustment, and topographic variation correction (where required). The GPR section depth estimates were based on a CMP-derived uniform velocity of 0.1 m/ns with a general error of ± 15%. An attenuation analysis of the recorded data has shown that, in general, the contaminated areas show a 2- to 6-dB level of signal attenuation increase (or enhancement) and that there is a small, but significant, frequency dependence to this attenuation (higher attenuation at lower frequencies). The most significant attenuation contrast occurs within a region that encompasses the very lower part of the vadose zone, the smeared zone, and the very top of the saturated zone.
 |
Contaminant Evaluation
|
|---|
The variable and disaggregated form of the contaminants in the smeared zone, plus the nature of GPR wave propagation in the near surface, means that it is difficult to truly determine whether the GPR signal response is related to scattering from low-permittivity, dispersed LNAPL materials or direct attenuation from the "bulk" lossy nature of the smeared zone as a whole. Attribute methods, whether they are amplitude or frequency based, only operate on the recorded GPR traces and cannot uniquely separate the scattered or reflected response of the subsurface from the incident wave field. This is particularly true in near-surface environments where the target objects are still within the near-field zone of the radiating antenna. In these cases, the resonant effect of the antenna's incident field (seen as a low-amplitude, extended-time "tail" in the recorded air or ground wave pulse) becomes an integral part of the recorded signal. The averaging effect of the attribute analysis method (i.e., the GPR signal is usually averaged across a specific time–space window, e.g., one pulse width and a number of traces) also compounds the problem as finer scale, subwavelength reflection or scattering responses tend to get combined into a overall amplitude (or frequency) value for a given depth or layer thickness.
To understand the GPR signal response at the site more comprehensively, a series of potential contaminant scenarios was devised to replicate the range of subsurface conditions observed throughout the smeared zone. If it is assumed that (i) the vadose zone is dry, sandy, low loss,
1 m thick with a relative permittivity of about
r = 5; (ii) the smeared zone is moderate-to-high loss,
0.5 m thick, and contains low-loss, low-permittivity (
r = 3–4) disaggregated, free-phase LNAPLs mixed with conductive, high-permittivity (
r = 80), high-loss contaminated groundwater; and (iii) the saturated zone contains mobile, contaminated groundwater mixed with clean groundwater, then there are five contaminant scenarios (or contaminant models) that describe the likely range of pore fluid conditions that would occur across the smeared zone (Fig. 2
):
- Uniform mixing of pore fluids. A fully saturated, single horizon containing different free-phase LNAPL and contaminated groundwater pore fluid ratios (note: mixing is at the micro or subpore scale with LNAPL/groundwater ratios of 10:90, 30:70, 50:50, 70:30, and 100:0).
- Small-scale, disaggregated LNAPL units. A single horizon containing centimeter-scale, randomly distributed, disaggregated "blobs" of fully saturated, free-phase LNAPL within a background of fully saturated, uniformly mixed contaminated groundwater. This model is intended to reproduce the macroscale disaggregation of the more viscous, less mobile LNAPLs into distinct, isolated units. As with the uniform mixing scenario above, the range of LNAPL/contaminated groundwater ratios is 10:90, 30:70, 50:50, 70:30, and 100:0.
- Larger scale pools of spatially distinct LNAPL separated into a number of thick lenses. Larger pools of flat-lying, free-phase LNAPL (with some associated centimeter-scale, isolated blobs) that sit within the fully saturated, uniformly mixed contaminated groundwater. The individual pools are laterally extensive, isolated, have lateral scales of between 0.5 and 1 m, are 0.03 to 0.04 m thick, and are randomly distributed into distinct horizons within the smeared zone. The model is intended to represent larger scale pooling of fully saturated free-phase LNAPL along specific, high-permeability pathways. This matches the observed features in some areas of the site where the contamination pathways appear to create sub-meter-scale hot spots of free-phase LNAPL saturation. The ratio of LNAPL/contaminated groundwater in this case is almost 25:75 (i.e., 25% LNAPL).
- Larger scale pools of spatially distinct LNAPL separated into a large number of thin lenses. As above, but with twice the number of thinner pools (0.17 m) and a ratio of LNAPL/contaminated groundwater of
20:80.
- Larger scale pools of spatially distinct LNAPL separated into a small number of thin isolated lenses. As above, but with half the number of thinner pools and a ratio of LNAPL/contaminated groundwater of
12:88.

View larger version (29K):
[in this window]
[in a new window]
|
FIG. 2. Modeling scenarios of mixed, contaminated pore and ground waters and "aged" light nonaqeous-phase liquids (LNAPLs) in the vadose and smeared zones based on the observed contamination characteristics.
|
|
In each case, there is a thin (0.088-m) capillary zone associated with the top of the smeared zone that represents the interface between the vadose zone (low loss) and the attenuating smeared zone. Note that the free-phase LNAPL within the smeared zone is never considered as being 100% pure (i.e., each macroscopic, free-phase LNAPL volume has some component of contaminated pore water associated with it—approximately 5%). This is based on site and laboratory observations where the free-phase LNAPL contaminants always seemed to contain some trapped water within their volume.
Although the geometry of the contaminant pools and blobs are designed to reflect the observed form of the LNAPL contamination, their size is important in terms of the probable range of GPR scattering or reflection responses. For an interface to generate coherent, observable reflection pulses in the recorded GPR trace, its area must be a substantial fraction of the GPR footprint. In other words, the isolated regions of permittivity contrast (i.e., the LNAPL pools) must be physically large enough to reflect a significant amount of coherent energy back to the receiving antenna. This minimum cross-sectional area is related to the antenna directivity (or beam width), signal frequency, depth, material permittivity, and relative antenna–target orientation. For flat-lying, semirough interfaces, it usually approximated to the size of the first Fresnel zone (refer to Reynolds [1997] for a more detailed explanation) and is a function of target depth and signal wavelength at the interface. At scales smaller than the first Fresnel zone, energy will start to be scattered throughout the subsurface, with increasingly incoherent responses being observed in the GPR traces. Figure 3
illustrates the range of Fresnel zone widths calculated for the permittivities and frequencies likely to be encountered within the smeared zone. In general, isolated bodies or pools of LNAPL would have to be between 0.5 and 1 m wide to generate coherent, individual reflection responses, which equates to the approximate size of the modeled pools and lenses. The LNAPL blobs of the small-scale, disaggregated LNAPL units contaminant model (Model 2), however, are much smaller than this (
0.018 m) and the GPR energy will be scattered in a random manner. The exact nature of the scattering response is dependent on the frequency, with lower frequencies having proportionally less energy scattered throughout the volume.

View larger version (15K):
[in this window]
[in a new window]
|
FIG. 3. First Fresnel zone diameters for a range of material relative permittivities and signal frequencies.
|
|
This effect can be understood more readily by considering the Mie scattering response of spherical objects as a function of target size, frequency, and scattered energy (Fig. 4
). Figure 4 illustrates the scattering characteristics of conducting, spherical objects in a uniform medium, where the radar cross-section effectively relates to the degree of back-scattered energy (Annan, 1999). Two target sizes are shown, 0.018 mm and 0.5 m, which represent the smallest LNAPL blob and first Fresnel zone sizes, respectively. For the full range of frequencies associated with the higher resolution antenna (i.e., 100–700 MHz), the scattering response of the 0.018-m objects is highly frequency dependent, with low frequencies scattering less energy than higher frequencies. This would result in a quantifiable frequency dependence to the scattering that enhances the low-frequency components of the propagating wave as it passes through the small-scale disaggregated materials. As the scale of the target increases, the scattering response becomes more uniform, with the different frequencies exhibiting similar characteristics. At 0.5 m and above, the target is large enough to become an individual reflector and produces a coherent response at all frequencies. The practical consequence of this behavior is a probable frequency dependence to the scattering within the smeared zone and changes in the spectrum of the recorded water table reflection.

View larger version (18K):
[in this window]
[in a new window]
|
FIG. 4. Scattering characteristics of 0.018- and 0.5-m-diameter spherical objects at 300, 400, and 500 MHz.
|
|
 |
Analysis Methodology
|
|---|
Laboratory Dielectric Analysis and Mixing Models
The effective permittivity of the site materials: groundwater (clean and contaminated), free-phase LNAPL, and clean, subsurface sands, were evaluated in the laboratory using a dielectric analysis technique incorporating an Agilent Technologies (Santa Clara, CA) 8753ES Automated Vector Network Analyzer and a series of calibrated test cells. The technique is well established and has been used to measure a wide range of subsurface materials including soils (e.g., Shang et al., 1999; Starr et al., 1999; Heimovaara et al., 1996; Olhoeft and Capron, 1993), rocks (e.g., West et al., 2003; Chelidze et al., 1999; Fam and Dusseault, 1998; Taherian et al., 1990), and, more pertinently, hydrocarbon-saturated materials (e.g., Hu and Liu, 2000). In this instance, the method developed by Cassidy (2001) was used, which follows the approach adopted by Baker-Jarvis et al. (1993, 1990) rather than the more commonly used Nicholson and Ross method (Nicholson and Ross, 1970; Nicholson, 1968). The method has the advantage of providing a solution that is irrespective of sample length or size and can be solved easily using automated, computational methods. Because the testing technique utilizes the reflection or transmission characteristics of a propagating electromagnetic (EM) wave (at the same frequencies as the GPR system), the measured effective permittivity results include all the loss effects associated with the analyzed material regardless of their physical attenuation mechanisms (e.g., dipolar losses, ionic or static conductivity losses, Maxwell–Wagner losses, etc. [Von Hippel, 1954]). A calibration routine is performed on known materials (polytetrafluoroethylene and deionized water) to reduce the experimental and instrument errors (which are predominantly systematic), therefore ensuring that the results are repeatable and absolute. The procedure is similar to standard, commercial-based dielectric testing calibration software (Agilent Technologies, 2001) but includes a repeat measurement and averaging phase to reduce the random calibration errors to a minimum. Multiple repeat measurements are automatically taken (typically 20–100 repeats depending on the material) and the results averaged to produce a final, mean effective permittivity reading at each frequency. Errors are typically in the order of ±0.2 to 0.8
r, with the lowest errors associated with the higher permittivity materials and at higher frequencies.
Examples of the measured, complex effective permittivity spectrums of 0.05-m-diameter, 0.04-m-thick representative samples of the site materials (i.e., "clean" and contaminated groundwater, free-phase LNAPL, and clean sands) are provided in Fig. 5
complete with measurement errors. The form of the groundwater's permittivity spectra is indicative of a dominant static or ionic conductivity loss superimposed on the permittivity behavior of free water (Von Hippel, 1954). This is supported by electrical conductivity (EC) and ionic chemistry analyses where the contaminated samples show an almost three-fold increase in EC and enhanced levels of ionic sulfate, Na, K, Mg, and Ca. This is consistent with the notion that the GPR signal attenuation is related to the presence of microbial degradation byproducts in the pore water. Interestingly, the free-phase LNAPL has an almost constant zero value for the imaginary component and suggests that, in general, the pure LNAPL can be considered as insulating or lossless for most practical applications.

View larger version (28K):
[in this window]
[in a new window]
|
FIG. 5. Measured effective permittivity spectra of the light nonaqueous-phase liquid (LNAPL), site sand, and "clean" and contaminated groundwater samples.
|
|
Effective Permittivity Mixing Models
To relate the measured effective permittivity characteristics to the GPR data effectively, it is necessary to combine the individual permittivity properties into a bulk spectrum that represents the true macroscopic behavior of the subsurface. Each material can be considered as macroscopic mixtures of solid (sand matrix), liquid (LNAPL and groundwater), and gaseous (air and vapor phase) components with a range of volumetric percentages. Laboratory measurements of carefully prepared mixtures can be achieved but it is time consuming, costly, and, with four potential mobile components (sand, LNAPL, water, and air), difficult to control physically. Consequently, theoretical or empirical mixing models have been developed that replicate the physical properties of the mixture by assuming that it is a multiphase combination of geometrically simple shapes, layers, or inclusions in a surrounding matrix (e.g., solid spheres in a fluid). There is a range of applicable formulations with slightly different approaches to determining the macroscopic permittivity (e.g., Fiori et al., 2005; Johnson and Poeter, 2005; Carcione et al., 2003; Hu and Liu, 2000; Tsui and Matthews, 1997; Sihvola, 2000, 1999, 1989; Sihvola and Alanen, 1991). Extensive experiments have demonstrated that, for relatively simple materials (sands, rocks, low-clay-content soils, etc.), these models correlate well with the experimental data. Of these, the complex refractive index model (CRIM) has become popular for hydrologic contaminant mixtures, as it is simple to apply, robust, and accurate throughout the GPR frequency range (e.g., Ajo-Franklin et al., 2004; Persson and Berndtsson, 2002; Darayan et al., 1998; Endres and Knight, 1992). Strictly a one-dimensional, layered medium model, CRIM has been shown to be effective for medium-to-coarse-grained, multiphase mixtures involving simple granular materials (e.g., semispherical sand grains, etc.) and moderate-to-low-viscosity fluids. It has the advantage of being a volumetric model that only requires knowledge of a material's permittivities and fractional volume percentage and can be used on both the real and imaginary components of the complex permittivity (Sihvola, 1999; Tsui and Matthews, 1997).
In its general form, the CRIM formula is written as
 | [1] |
where
mixe is the bulk effective permittivity of the mixture, fi is the volume fraction of each ith component, and
i is the permittivity of each ith component. Any number of phases can be included but, in this case, a four-phase model has been chosen, with
w,
g,
m, and
n representing the measured effective permittivities of water, gas and air, matrix, and LNAPL phases, respectively. This formula needs to be modified, however, to account for the measured effective permittivity of the sand samples, which can be considered as a two-phase, predetermined mixture of sand grains and air that is a submixture of the main CRIM model. The measured permittivity of the site sands (
ss) is related to the matrix permittivity of the CRIM model (
m) by
 | [2] |
where
s is the porosity fraction of the measured sands (0.42 or 42% in this instance) and
sg is the permittivity of the vapor and air in the pore spaces. The modified CRIM formula for the four phases now becomes
 | [3] |
where Sw is the water saturation index and Sn the LNAPL saturation index (i.e., percentage of pore space filled with groundwater and LNAPL, respectively). Using complex permittivity variables, this formula can be used to determine both the real and imaginary components of the "bulk" effective permittivity of any specific mixture and will be a more effective and accurate representation of the true material's properties. Figure 6
illustrates both the real and imaginary components of the effective permittivity spectrum for a range of CRIM-based mixtures and their corresponding FDTD material parameter models. Note that although errors bars are not shown for visual clarity, the typical permittivity error is between ±5 and 15%.

View larger version (40K):
[in this window]
[in a new window]
|
FIG. 6. Effective (relative) permittivity spectrum of Complex Refractive Index Model (CRIM)-based, light nonaqueous-phase liquid (LNAPL)–groundwater mixtures and their generalized finite-difference, time-domain (FDTD) parameter models.
|
|
Finite-Difference Time-Domain Modeling
For complex electromagnetic wave propagation problems, such as GPR, the FDTD technique has become one of the most popular modeling tools, particularly with the recent increase in accessible and inexpensive computational resources. In comparison to integral methods, the technique has specific advantages (Cassidy, 2001), namely:- It uses the direct, time-domain implementation of Maxwell's electromagnetic field equations in three-dimensional space.
- It can be solved with explicit, closed form equations using matrix, parallel, or incremental methods and is robust, accurate, and relatively straightforward to implement.
- It provides a total-field solution (i.e., fully three-dimensional) that operates on both the electric (E) and magnetic (H) field vectors in both time and space.
- It can be applied to broadband, narrow-band, or harmonic time-domain problems without having to change the computation scheme or the inherent properties of the model.
- It incorporates a wide range of conductive, lossy, and dispersive materials without the need to alter the mathematical description of the scheme.
- It can include arbitrary, three-dimensional subsurface geometries, complex material features, and sophisticated antenna designs through the use of different grid schemes and layouts.
The flexibility of the technique has led to a number of different FDTD formulations, usually designed around specific applications (Taflove, 1995, 1998); however, each scheme tends to share common characteristics. The subsurface volume is usually subdivided into a three-dimensional, orthogonal grid of individual "field cells" ranging between 1/8 and 1/15 of the dominant signal wavelength. Within each individual cell, E and H are described by discrete field components (Ex, Ey, Ez, Hx, Hy, Hz) that are uniformly staggered in space and time. For orthogonal grids, the most common form is the Yee Cell (Yee, 1966), with magnetic and electric field components staggered at half-cell intervals. With the use of a finite-difference approximation to the differential form of Maxwell's electromagnetic field equations, it is possible to calculate the electric field at any point in space (and time) from a knowledge of its neighboring magnetic fields, and vice versa. If an incremental, time-varying electric source is applied to one or more of the field components in the volume (as with the GPR antenna signal) then, for a specific time increment, the total three-dimensional EM field is calculated by computationally stepping through each of the alternating E and H field calculations of every individual cell in sequence. This procedure is then repeated for the next time increment, etc., resulting in the time-dependent propagation of the EM wave throughout the whole volume. To model the subsurface physically, individual material properties are assigned to each cell (i.e., a permittivity, magnetic permeability, and conductivity) and geometrically complex structures are constructed by grouping together cells of equal properties.
For our study site examples, the propagating GPR wave was modeled by a total-field, explicit, robust scheme (Cassidy, 2001) that is similar to the formulations of Giannopoulos (2005), Bergmann et al. (1998), and Teixeria et al. (1998). The scheme provides accurate models of GPR wave velocity and dispersion without the need for fine spatial sampling. It includes a wide range of materials, each having frequency-dependent permittivity characteristics and contains realistic antenna designs including shields, efficient signal damping, and accurate source signal descriptions. "Memory variables" are used to determine the dispersive, frequency-dependent effect of the materials and the scheme is an adaptation of the visco-elastic modeling approach of Blanch et al. (1994). Full details on the mathematical and computational properties of the scheme can be found in Cassidy, (2005, 2004a, 2001) and Bergmann et al. (1998), along with the appropriate stability conditions, error estimates, and temporal and spatial sampling criteria, etc. The complex, effective permittivity spectrum of each material is modeled by a combined, complex permittivity–conductivity relaxation function that is based on the visco-acoustic and electromagnetic analogy of linear components (Carcione and Cavallini, 1995). In this instance, the complex permittivity spectrum
*(
) is described by a superposition of individual Deybe electric field and electric flux density relaxation times (
El and
Dl, respectively) combined with a static permittivity,
s (Bergmann et al., 1998):
 | [4] |
The complex conductivity
*(
) is described by a static conductivity
s and a conductivity relaxation time, 
:
 | [5] |
In any particular material, there are different overlapping relaxation mechanisms, each resulting from the influence of pore fluids, grain interfaces, surface charges, etc. The combination of these individual, yet independent, relaxation times allows the frequency-dependent EM response of any material to be fully described, however complex. For the study site's material mixtures, the bulk effective permittivity spectrum is represented by a simplified relaxation response of one or two different permittivity relaxation mechanisms and a single conductivity relaxation. The relaxation parameters are optimized to fit the general form of the real and imaginary effective permittivity components across the whole of the site-specific frequency range (i.e., 100–700 MHz). Representative examples are provided in Fig. 6 along with their associated FDTD modeling parameters in Table 1
.
View this table:
[in this window]
[in a new window]
|
TABLE 1. Finite-difference, time-domain modeling, effective permittivity relaxation parameters for the contaminated water and light nonaqueous-phase liquid (LNAPL) mixtures in Fig. 6.
|
|
 |
Numerical Modeling Examples
|
|---|
Although both 225- and 450-MHz surveys were collected at the site, only the higher frequency data (i.e., from the 450-MHz antenna) have been modeled in the following examples. As the modeled spectrum of this antenna is relatively broad (i.e., a bandwidth of
100–700 MHz), it allows us to investigate the shallow (1.2–1.5-m) vadose and smeared zone to a high degree of spatial resolution yet retain the lower frequency information of the 225-MHz data. To model each of the potential contaminant scenarios identified above, a base (or background condition) three-dimensional FDTD model was built that replicates the clean conditions of the site (i.e., uniform, unsaturated dry sands with 5% water content and a clean groundwater table at 1.48 m). Superimposed onto this base model is each of the contaminant scenarios, with the smeared zone being an identical thickness in each case. The FDTD computational domain is built from a volume of 56 (x) by 36 (y) by 110 (z) equal-sided, orthogonal cells of 0.0177-m width. Two grid cells per face are used for the boundary calculations and the axis of the centrally located antenna pair is orientated in the y direction. This model geometry is shown in schematic form in Fig. 7
, along with the corresponding contaminated, smeared zone model. For each modeling scenario, a single trace, relating to the higher frequency 450-MHz antenna, is simulated over a two-way travel time of 56 ns, with the antenna and base model conditions remaining the same. All the relevant modeling parameters are provided in Table 2
(where 
= 0 in each case) and the materials are considered as having nonmagnetic properties (i.e., the losses are associated with the permittivity and conductivity only). In each case, only the properties of the smeared zone are changed, with the following models representing each of the contaminant scenarios:- Uniform mixing of pore fluids: 0.31-m-thick smeared zone, uniformly mixed LNAPL/contaminated pore water ratios of 30:70, 50:50; 70:30, and 0:100
- Small-scale, disaggregated LNAPL units: 0.31-m-thick smeared zone with randomly distributed LNAPL blobs within each layer of the model. Volumetric LNAPL/contaminated groundwater ratios are 30:70, 50:50, and 70:30, with the final model having full, free-phase LNAPL saturation (see examples in Fig. 8)
- Larger scale pools of spatially distinct LNAPL separated into a number of thick lenses: 0.31-m-thick smeared zone with
0.5-m-diameter, 0.0354-m-thick lenses randomly distributed into four, evenly spaced horizons (Fig. 9). Model includes a 5% (by vol.) amount of isolated 0.018-m LNAPL blobs within the contaminated horizons
- Larger-scale pools of LNAPL separated into a large number of thin lenses: As above, but with twice the number of equally spaced thinner pools (0.017 m thick)
- Larger-scale pools of LNAPL separated into a small number of thin isolated lenses: As above, but with half the number of thinner pools (not equally spaced)

View larger version (36K):
[in this window]
[in a new window]
|
FIG. 7. Finite-difference, time-domain (FDTD) computational model for the clean and contaminated areas of the study site. Numbers in square brackets relate to the x–y–z coordinates of the individual field cells. The effective permittivity relaxation parameters for the FDTD models are provided in Table 2.
|
|
View this table:
[in this window]
[in a new window]
|
TABLE 2. Finite-difference, time-domain modeling, effective permittivity relaxation parameters for the site models.
|
|

View larger version (67K):
[in this window]
[in a new window]
|
FIG. 8. Distribution of disaggregated light nonaqueous-phase liquid (LNAPL) "blobs" (black squares and gray cubes) within the main computational domain (boundary cells not shown). Individual layer examples taken from a single layer within the (A) 10%, (B) 30%, (C) 50%, and (D) 70% LNAPL volumetric saturation models.
|
|

View larger version (35K):
[in this window]
[in a new window]
|
FIG. 9. (A–D) Distribution of thick light nonaqueous-phase liquid (LNAPL) "pools" (black areas) within the main computational domain (boundary cells not shown). Each pool is two cells thick (0.0354 m) and separated by a two-cell horizon containing no free-phase LNAPL. (E) Three-dimensional distribution of thin LNAPL "lenses" within the main computational domain of Model 5—"a small number of thin isolated lenses" (boundary cells not shown).
|
|
In addition to the contaminant scenario models above, a separate "incident field" model (simulating the GPR response in the vadose zone only) was generated to help characterize the nature of the recorded near-field pulse and includes the air and ground waves generated by the GPR system. This model is identical to the base clean area model but without the saturated groundwater.
 |
Numerical Modeling Results
|
|---|
The results of the modeling are shown in Fig. 10 to 13

, with Fig. 10 illustrating the recorded trace characteristics for the incident field and base clean area models and Fig. 11 to 13
the uniform mixing, disaggregated scattering, and large pool models, respectively. There are some interesting aspects to the incident field and base clean area models that provide an insight into the complicated nature of GPR signal characterization in the near surface. The recorded trace from the incident field model shows the characteristic dual-peak, near-field ground and air wave signature of a well-coupled, shielded antenna with what appears to be little "ringing" in the recorded signal. The spectrum is smooth, almost symmetrical, and centered around a value that is close to the antennae's quoted frequency. Because the air–ground interface of the model is effectively flat and smooth (with respect to the signal wavelength), this model represents almost ideal antenna-to-ground coupling conditions with no discernable air gap between the antennae and surface. If we consider a time frame that represents the likely two-way travel time of signals emanating from the smeared zone (
17–39 ns), however, it is clear that the incident field still contains a significant amount of energy when compared with the amplitude of the water table reflection from our base model. When compared with the main ground and air wave signature, the frequency spectrum is narrower, with a "down-shifted" central peak at
200 MHz. These features are consistent with the effect of late-stage mutual coupling between the transmitting and receiving antennae, which produces a decreasing-amplitude, low-frequency tail to the incident signal. Unfortunately, the base model trace contains both the low-frequency tail of the incident field and the reflected energy from the water table interface. As the recorded signal is a superposition of the two propagating EM fields within the capture area of the receiving antenna, it is impossible to separate the two without an accurate knowledge of the incident signal. This is possible with our models, as the incident field can be simulated separately, but is very difficult to achieve with real data. A range of different processing methods can be used in an attempt to extract the reflected response from the background (e.g., background removal, predictive deconvolution, etc.) but the success of these methods can never be guaranteed. Consequently, the frequency spectrum of the real near-surface signals will always contain a component of the incident field's spectrum and, therefore, attribute analysis methods may be misrepresentative of the true reflection and scattering response. The influence of the incident field on the base model's spectrum is shown in the central spectrum plot of Fig. 10, where the profile is split into two components: a narrow, low-frequency component, which is similar in amplitude and shape to the incident field spectrum, and a broader, higher frequency component that falls between
300 and 500 MHz. It is reasonable to assume that this higher frequency component represents the true spectral characteristics of the water table reflector, as significant dispersion (through the loss of higher frequencies) will not occur in the overlying, low-loss sands.

View larger version (27K):
[in this window]
[in a new window]
|
FIG. 10. Finite-difference, time-domain modeling results for the "incident field" and "base clean area" models. Upper plots show the spectrum and recorded trace of the whole "incident field" model, whereas the central and lower plots show the spectrum and recorded trace from the smeared-zone region only (i.e., the near-field response has been omitted). In the central two plots, the incident and base fields are shown together while the lower two plots illustrate the comparison between the incident field and the 70:30 contaminated porewater/light nonaqueous-phase liquid (LNAPL) ratio, uniform mixing model.
|
|

View larger version (32K):
[in this window]
[in a new window]
|
FIG. 11. Finite-difference, time-domain modeling results for the "uniform mixing" models (compared with the "base clean model") with varying ratios of contaminated porewater/light nonaqueous-phase liquid (LNAPL). In each case, the data are from the smeared zone and have had the incident field component removed from the spectral response (i.e., the illustrated spectrum represents the frequency-dependent characteristics of the smeared zone's materials only).
|
|

View larger version (32K):
[in this window]
[in a new window]
|
FIG. 12. Finite-difference, time-domain modeling results for the "disaggregated scattering" models (compared with the "base clean model") with varying percentages of free-phase light nonaqueous-phase liquid (LNAPL) "blobs." In each case, the data are from the smeared zone only and have had the incident field component removed from the spectral response.
|
|

View larger version (32K):
[in this window]
[in a new window]
|
FIG. 13. Finite-difference, time-domain modeling results for the "large pool" models (compared with the "base clean model") with varying geometrical arrangements. In each case, the data are from the smeared zone only and have had the incident field component removed from the spectral response. Note that the 100% light nonaqueous-phase liquid (LNAPL) results represent a single layer of pure, free-phase LNAPL having the same thickness as the smeared zone.
|
|
Without being able to extract the reflected and scattered field data, incorrect interpretations could easily be made on the attenuation characteristics of the smeared zone's materials. This can be seen in the lower spectral plot of Fig. 10, where the recorded trace from the incident field model is compared with the 70:30 contaminated pore water/LNAPL, uniform mixing model. Taken "as read," the spectrum of signals from the smeared zone would appear to be skewed toward lower frequencies with a peak at approximately 200 MHz, which implies that the materials are preferentially attenuating higher frequencies. When the incident field spectrum is subtracted from the mixing model spectrum (as illustrated in the 70:30 model of Fig. 11), however, the frequency-dependent characteristics of the materials become clearer, with a broader response and high-frequency spectral components. Consequently, all the spectra shown in Fig. 11 to 13
have had the incident field response subtracted from the data, including the base models used as the clean area comparison. The upper plot in each figure is the 100% LNAPL saturation model and is intended to act as a comparative control (i.e., the temporal and spectral response of each modeling scenario should be the same in each case).
Uniform Mixing Models
When compared with the base model, the contaminated models (Fig. 11) show a rising level of general signal attenuation (i.e., the temporal and spectral signal amplitudes are both reducing) with increasing amounts of contaminated pore water. Interestingly, the 30:70 and 50:50 ratio models show comparable temporal characteristics with similar peak and mean pulse amplitudes and travel times for individual reflections within the smeared zone. If standard, temporal-based attribute analysis methods were used on these traces, there would be little identifiable difference between the clean and contaminated data and, therefore, no indication of the increased levels of contaminated pore water. Only at a contaminated pore water/LNAPL ratio of 70:30 does the temporal amplitude reduce enough to be potentially identifiable, with a reduction in mean and peak amplitude of approximately 50%. This equates to an attenuation increase or enhancement of approximately 6 dB, which is similar to the maximum level of enhancement observed at the study site. In contrast, the form of the modeled signal spectrum changes significantly with increasing contaminated pore water/LNAPL ratios (>50:50), with the higher frequency components being preferentially reduced and the profile showing a distinct, low-frequency skew. This loss of high frequencies would, at first, appear inconsistent with the laboratory analysis and mixing model results for the contaminated mixtures, which show the lower frequencies (<300 MHz) as having an increased imaginary effective permittivity (and therefore a greater degree of loss). The actual GPR signal attenuation (in decibels per meter), however, is not just dependent on the imaginary component of the effective permittivity, but is, instead, a combined function of the loss tangent (i.e., the ratio of imaginary and real components) plus the frequency (Daniels, 1996). Therefore, individual materials that exhibit an apparent low-frequency bias to the effective permittivity spectrum can actually produce an increased attenuation effect at higher frequencies despite the form of the imaginary component.
Disaggregated LNAPL Scattering Units or Blobs Models
The spectral and temporal response of the scattering models (Fig. 12) is noticeably different than the equivalent, uniform mixing models (i.e., 30, 50, and 70% volumetric LNAPL saturation), and, in general, show a higher degree of temporal mismatch to the base model. Instead of seeing a proportional reduction in overall signal amplitude (relative to the base model) as the percentage of LNAPL drops, we have a low-amplitude pair (50 and 30% volumetric saturation) and a corresponding high-amplitude pair (10 and 70%). The temporal and spectral form of the low-amplitude pair is reasonably similar, with an approximate reduction in peak and mean signal amplitude of 50 to 60% (or 6–8-dB attenuation increase or enhancement) and a relatively broad, symmetrical spectrum. In contrast, the high-amplitude, low volumetric saturation model (10%) has a weak water table reflection and a spectrum skewed toward high frequencies, while the 70%, high volumetric saturation model shows a strong water table reflection and broader spectrum. These responses reflect the differential nature of scattering within the smeared zone, where, at low saturations, the LNAPL blobs are too isolated and dispersed to have any significant influence on the general attenuation properties. In this instance, the high permittivity of the volumetrically dominating groundwater means that the smeared zone is acting like a single water-saturated layer, with a large permittivity contrast at its upper surface and low contrast at its base. Therefore, a higher proportion of the incident signal is reflected back from the top of the zone (with a relatively high-frequency component) and less transmitted into its middle. As the volumetric saturation increases, the blobs form larger coalesced units and account for a greater component of the smeared zone's overall permittivity. As a result, less energy is reflected off the top of the zone and more passes down into the contaminated materials, where it is scattered and attenuated in a uniform manner. In this instance, it would appear that the combination of both mechanisms results in higher degrees of overall signal attenuation (when compared with the uniform mixing model) and a broader, symmetrical signal spectrum. As the LNAPL saturation increases above 70%, the low-loss nature of free-phase materials starts to dominate the signal response, with less overall attenuation occurring, a reduced level of scattering, and the production of a strong water table reflection that is greater than that of the base model. This would create a bright spot in the GPR section (with a relative signal amplitude increase of
2–3 dB) that would be strong enough to be identified by attribute analysis methods.
Large Pools and Lenses Models
In this example (Fig. 13), the form and amplitude of the temporal and spectral responses is generally similar to the high-saturation-percentage model above. The "thick pools" model produces a very strong, high-amplitude reflection response (
5 dB above the base model), with individual pulses that are regular, narrow, and similar in size. This temporal form is reflected in the spectral response, with a high-frequency skewed profile that is an order of magnitude larger than the corresponding base model. These features are consistent with the response of multiple, coherently reflecting bodies acting almost as individual layers within the smeared zone. Because there is (i) a large permittivity contrast between the LNAPL pools and the background groundwater-saturated sands and (ii) a wide separation between individual pools, a significant amount of constructively interfering energy is reflected back to the surface, producing the regular, high-amplitude response seen in the recorded trace. Relatively little energy is transmitted down into the saturated zone, and, as a result, the GPR section will show very bright reflections in the smeared zone and a shadow, or dim, region beneath it. The "small number of thin lenses" model exhibits similar characteristics, but with reduced temporal and spectral amplitudes, so it is reasonable to assume that the same physical processes are operating, albeit in a weakened manner. Only the "large number of thin lenses" model generates a significant reduction in peak and mean signal amplitude, with a general level of relative attenuation increase or enhancement on the order of 4 to 5 dB. The spectral response, however, is still broad and exhibits a distinct, if reduced, high-frequency bias. Given the large number of thin pools, and the small gaps that exist between them, it is reasonable to expect a considerable amount of overlap and "reverberation" occurring in the reflected signals, with the likelihood of significant destructive interference in the recorded signal. In some ways, this is similar to the physical effect of scattering in the resonance region (see Fig. 4) but with a more predictable and coherent frequency response. With a volumetric LNAPL saturation percentage of approximately 25%, this model is similar to the 70:30 ratio uniform mixing model and 30% LNAPL saturation, disaggregated scattering model and their temporal responses are similar. The spectral response is slightly different (i.e., skewed to higher frequencies), however, with a peak that matches the incident pulse frequency. This suggests that material attenuation and scattering effects are not dominating the response of this model and that the signals are being reflected back to the surface relatively unmodified.
 |
Implications for the Interpretation of Site Data
|
|---|
In general, the characteristics of each family of models can be summarized as follows.- Uniform mixing: Overall, there is an increasing level of temporal and spectral attenuation and a distinct loss of higher frequencies with increasing contaminated pore water saturation. Relative signal attenuation increase or enhancement reaches identifiable levels only above 50% contaminated pore water saturation (i.e., relative attenuation is >2–3 dB), with the water table reflection showing increasing travel times and decreasing signal amplitudes with increasing pore water saturation.
- Disaggregated LNAPL scattering units: The temporal and spectral attenuation is greatest at moderate LNAPL saturations (30–50%), with identifiable levels of relative signal attenuation increase or enhancement. The spectral response is generally broad, with a symmetrical form around a central frequency of
400 MHz. The water table reflection shows increasing travel times and decreasing amplitude with decreasing LNAPL saturation.
- Large pools and lenses: There is a significant increase in temporal and spectral amplitudes for the widely spaced thick and thin layer models, while the multiple, thin layers model shows a decrease in signal amplitudes and identifiable levels of relative signal attenuation increase or enhancement. Overall, the spectral responses are broad, with a high-frequency skew toward
450 MHz. The water table reflector is only distinct in the multiple, thin layers model.
When compared with the base clean model, the temporal characteristics of the uniform mixing models can be considered as fairly predictable, with amplitudes decreasing and travel times increasing at higher contaminated pore water saturations (i.e., attenuation and velocity is increasing within the smeared zone). Similarly, the scattering response of the disaggregated LNAPL models is as expected, with the greatest degree of attenuation occurring when the scattering objects form a high-density volume of larger, yet still distinct, individual units. The relative attenuation increase in both cases is comparable to the site data from the smeared zone (
2–6 dB), which suggests that either (or both) mechanisms are occurring and that the pore water and groundwater saturation levels are likely to be on the order of 30 to 50%. This is consistent with the findings of Cassidy (2006, 2007), who estimated contaminated groundwater saturation levels at the site to be
40% (through an attenuation analysis of the GPR data and direct laboratory measurements). One slight surprise is the degree of relative signal increase associated with the thick or thin, widely spaced pools models when compared with the multiple, narrowly spaced pools models. Although the individual pools are approximately the same size (
0.5 m in diameter and slightly less than the size of the first Fresnel zone), the thick-pools model has peak amplitudes that are approximately 8 to 9 dB stronger than the multiple, thin, narrowly spaced pool model. In fact, the thick-pools model produces pulse amplitudes that are greater than the 100% LNAPL (or single thick layer) model and would, therefore, produce a very strong instantaneous or mean amplitude result if attribute analysis methods were used on the data. These apparently anomalous characteristics illustrate the influence of finer scale lenses, or layers, on the recorded signal, where constructive and destructive interference can dominate the temporal response. This can lead to significantly strong bright and dim spots being present in the data and attribute analysis results that are unrelated to the direct effect of material attenuation and scattering.
In contrast to the temporal characteristics, the spectral response of each family is more regular, with a high-frequency skew to the pool models, a loss of high frequencies in the uniform mixing models, and a generally symmetrical profile in the scattering models. This would imply that the frequency-dependent characteristics could be indicative of the mode, or nature, of the signal response for a specific LNAPL or pore water and groundwater saturation level. In our idealized modeling case, where the incident field spectrum can be accurately subtracted from the recorded signal, this is achievable, and an example is shown in Fig. 14
where the pool, uniform mixing, and disaggregated scattering models are shown for comparable LNAPL saturations of 25 to 30%. The temporal form of the modeled signals is very similar, with comparable pulse amplitudes and travel times, but the spectral response is distinct enough to be specifically identifiable. This is based, however, on the accurate subtraction of the incident field from the recorded data, which, in practice, is difficult to achieve with real, near-surface GPR data. Consequently, the slight frequency dependence observed in the study site's GPR data (i.e., a greater degree of attenuation enhancement at lower frequencies, <300 MHz) cannot be uniquely attributed to any specific attenuation mode (Cassidy, 2007). Nevertheless, this high-frequency bias to the recorded data suggests that the LNAPLs are forming spatially distinct, yet isolated, pools within the smeared zone, which is consistent with the field-based observations.

View larger version (32K):
[in this window]
[in a new window]
|
FIG. 14. Comparison between the finite-difference, time-domain modeling results for 25 to 30% light nonaqueous-phase liquid (LNAPL) saturation, "large number of thin pools" model, 70:30 contaminated porewater/LNAPL ratio, "uniform mixing" model, and the 30% LNAPL "disaggregated scattering" model (compared with the "base clean model"). In each case, the data are from the smeared zone only and have had the incident field component removed from the spectral response.
|
|
 |
Conclusions
|
|---|
Through the use of dielectric material property analysis, CRIM-based mixing model evaluation, and FDTD modeling, it has been possible to evaluate the form and nature of GPR signal attenuation within the near-surface (<2-m), contaminated vadose and smeared zone of a mature LNAPL-contaminated site. A comparison between "clean" and contaminated FDTD models having different LNAPL/pore water or groundwater saturations and contaminant geometries has shown that uniform mixtures, disaggregated LNAPL, scattering mixtures, and mixed LNAPL pools or lenses exhibit different temporal and spectral responses across the near-surface GPR frequency range of approximately 100 to 700 MHz. In general, the temporal responses show predicable characteristics, with signal attenuation increase (or enhancement) being related to increased levels of lossy, contaminated pore water and groundwater or variations in signal scattering associated with changes in the form, or geometry, of the LNAPL units. In contrast, the spectral responses show a distinct mixture-related pattern, with uniform mixtures showing relative loss at higher frequencies, scattering mixtures having symmetrical spectral characteristics, and pools or lenses exhibiting a significant high-frequency enhancement. This implies that the spectral characteristics can be used to provide important, additional information on the mode and nature of the attenuating mechanisms as long as the "true" GPR response of the contaminating materials can be determined. For deeper applications (i.e., where the two-way travel time and amplitude of the contaminant-related signals is greater than the tail of the incident signal), near-field and background effect removal should be relatively straightforward. For near-surface applications, however, accurately describing the incident field characteristics can be difficult, particularly in heterogeneous environments, and more advanced processing methods (such as wavelet analysis, predictive correlation, etc.) may be required before useful spectral information can be extracted from the GPR data. This is also likely to be the case where significant, subwavelength-scale natural heterogeneity occurs in the subsurface that is unrelated to the contaminant properties or geometries. If the effective permittivity contrast of these natural variations is much greater than that of the contaminant or background, then, in reality, it will be almost impossible to separate the contaminant signal response from the natural scattering (or clutter) effects of the subsurface. Despite these drawbacks, the modeling results from this study appear to be consistent with the real GPR data and show that, with suitable incident field descriptions, attenuation increase or enhancement can be related to the saturation index (and possibly contaminant form), albeit in an approximate way.
 |
ACKNOWLEDGMENTS
|
|---|
I would like to thank Emma Pell for her work on the contaminant chemistry and Fiona Donaghy for her help in proofreading the manuscript. This research has been funded by NERC Grant NER/2002/00100 and NERC Geophysical Equipment Pool Loan no. 730.
 |
REFERENCES
|
|---|
- Abdel Aal, G.Z., L. Slater, and E.A. Atekwana. 2006. Induced-polarization measurements on unconsolidated sediments from a site of active hydrocarbon biodegradation. Geophysics 71(2):H13–H24.[CrossRef]
- Agilent Technologies. 2001. 85070 dielectric probe kit software vD1.00. Agilent Technologies, Santa Clara, CA.
- Ajo-Franklin, J.B., J.T. Geller, and J.M. Harris. 2004. The dielectric properties of granular media saturated with DNAPL/water mixtures. Geophys. Res. Lett. 31:L17501, doi:10.1092/2004GL020672.[CrossRef]
- Annan, A.P. 1999. Ground penetrating radar: Workshop notes, June 1999. Sensors and Software, Mississauga, ON.
- Atekwana, E.A., E. Atekwana, F.D. Legall, and R.V. Krishnamurthy. 2004a. Field evidence for geophysical detection of subsurface zones of enhanced microbial activity. Geophys. Res. Lett. 31:L23603, doi:10.1029/2004GL02157.[CrossRef]
- Atekwana, E.A., E.A. Atekwana, R.S. Rowe, D.D. Werkema, and F.D. Legall. 2004b. Total dissolved solids in ground water and its relationship to bulk conductivity of soils contaminated with hydrocarbons. J. Appl. Geophys. 56:281–294.[CrossRef]
- Atekwana, E.A., W.A. Sauck, G.Z.A. Aal, and D.D. Werkema, Jr. 2002. Geophysical investigations of vadose zone conductivity anomalies at a hydrocarbon contaminated site: Implications for the assessment of intrinsic bioremediation. J. Environ. Eng. Geophys. 7:103–110.
- Atekwana, E.A., W.A. Sauck, and D.D. Werkema. 2000. Investigations of geoelectrical signatures at a hydrocarbon contaminated site. J. Appl. Geophys. 44:167–180.[CrossRef]
- Atekwana, E.A., D.D. Werkema, Jr., J.W. Duris, S. Rossbach, E.A. Atekwana, W.A. Sauck, D.P. Cassidy, J. Means, and F.D. Legall. 2004c. In-situ apparent conductivity measurements and microbial population distribution at a hydrocarbon contaminated site. Geophysics 69:56–63.[CrossRef]
- Baedecker, M.J., and I.M. Cozzarelli. 1994. Biogeochemical processes and migration of aqueous constituents in ground water contaminated with crude oil. p. 69–79. In A.R. Dutton (ed.) Toxic substances and the hydrologic sciences. Am. Inst. of Hydrol., Austin, TX.
- Baker-Jarvis, J., M.D. Janezic, J.H. Grosvenor, and R.G. Geyer. 1993. Transmission/reflection and short-circuit line methods for measuring permittivity and permeability. NIST Tech. Note 1355 (revised). Natl. Inst. of Standards and Technol., Boulder, CO.
- Baker-Jarvis, J., E.J. Vanzura, and W.A. Kissick. 1990. Improved technique for determining complex permittivity with the transmission/reflection method. IEEE Trans. Microwave Theory Tech. 38:1097–1103.
- Benson, A.K. 1995. Applications of ground penetrating radar in assessing some geological hazards: Examples of groundwater contamination, faults, cavities. J. Appl. Geophys. 33:177–195.[CrossRef]
- Bergmann, T., J.O.A. Robertsson, and K. Holliger. 1998. Finite-difference modeling of electromagnetic wave propagation in dispersive and attenuating media. Geophysics 63:856–867.[CrossRef]
- Blanch, J.O., J.O.A. Robertsson, and W.W. Symes. 1994. Viscoelastic finite-difference modeling. Geophysics 59:1444–1456.[CrossRef][Web of Science]
- Bradford, J.H. 2003. GPR offset-dependent reflectivity analysis for characterization of a high-conductivity LNAPL plume. p. 238–252. In Proc. Symp. on the Application of Geophysics to Engineering and Environmental Problems (SAGEEP), San Antonio, TX. 6–10 Apr. 2003. Environ. Eng. Geophys. Soc., Denver, CO.
- Burton, B.L., G.R. Olhoeft, and M.H. Powers. 2004. Frequency spectral analysis of GPR data over a crude oil spill. p. 267–270. In E.C. Slob et al. (ed.) Proc. Int. Conf. on Ground Penetrating Radar, 10th, Delft. 21–24 June 2004. Delft Univ. of Technol., Delft, the Netherlands.
- Campbell, D.L., J.E. Lucius, K.J. Ellefsen, and M. Deszcz-Pan. 1996. Monitoring of a controlled LNAPL spill using ground penetrating radar. p. 511–517. In Proc. Symp. on the Application of Geophysics to Engineering and Environmental Problems (SAGEEP), Keystone, CO. 28 Apr.–1 May 1996. Environ. Eng. Geophys. Soc., Denver, CO.
- Carcione, J.M., and F. Cavallini. 1995. On the acoustic–electromagnetic analogy. Wave Motion 21:149–162.[CrossRef][Web of Science]
- Carcione, J.M., D. Gei, M. Botelho, A. Osella, and M. del la Vega. 2006. Fresnel reflection coefficients for GPR–AVA analysis and detection of seawater and NAPL contaminants. Near Surface Geophys. 4:253–263.
- Carcione, J.M., and G. Seriani. 2000. An electromagnetic modelling tool for the detection of hydrocarbons in the subsoil. Geophys. Prospect. 48:231–256.[CrossRef]
- Carcione, J.M., G. Seriani, and D. Gei. 2003. Acoustic and electromagnetic properties of soils saturated with salt water and NAPL. J. Appl. Geophys. 52:177–191.[Medline]
- Cassidy, D.P., D.D. Werkema, W.A. Sauck, E.A. Atekwana, S. Rossbach, and J. Duris. 2001. The effects of LNAPL biodegradation products on electrical conductivity measurements. J. Environ. Eng. Geophys. 6:47–53.
- Cassidy, N. 2004a. Application of numerical modelling for the interpretation of near-surface ground penetrating radar. p. 47–59. In D.J. Daniels (ed.) Ground penetrating radar. 2nd ed. Radar, Sonar, Navig. and Avionics Ser. 15. Inst. of Electrical Eng., London.
- Cassidy, N.J. 2001. The application of mathematical modelling in the interpretation of ground penetrating radar data. Ph.D. thesis. Keele Univ., Keele, UK.
- Cassidy, N.J. 2004b. Dielectric properties of free-phase hydrocarbon contamination: Implications for GPR investigation. p. 551–554. In Proc. Int. Conf. on Ground-Penetrating Radar (GPR2004), 10th, Delft, the Netherlands. 21–24 June 2004. Inst. Electrical Electron. Eng., New York.
- Cassidy, N.J. 2005. The practical application of numerical modelling for the advanced interpretation of ground-penetrating radar. p. 105–113. In Proc. Int. Workshop on Advanced Ground Penetrating Radar, 3rd, Delft, the Netherlands. 2–3 May 2005. Inst. Electrical Electron. Eng., New York.
- Cassidy, N.J. 2006. GPR signal dispersion and attenuation of LNAPL hydrocarbon contamination. p. 1–8. In Proc. Int. Conf. on Ground Penetrating Radar (GPR2006), 11th, Columbus, OH. 19–22 June 2006. Inst. Electrical Electron. Eng., New York.
- Cassidy, N.J. 2007. Evaluating LNAPL contamination using GPR signal attenuation analysis and dielectric property measurements: Practical implications for hydrological studies. J. Contam. Hydrol. (in press), doi:10.1016/j.jconhyd.2007.05.002.
- Chelidze, T.L., Y. Gueguen, and C. Ruffet. 1999. Electrical spectroscopy of porous rocks: A review: II. Experimental results and interpretation. Geophys. J. Int. 137:16–34.[CrossRef]
- Cozzarelli, I.M., M.J. Baedecker, R.P. Eganhouse, and D.F. Goerlitz. 1994. The geochemical evolution of low-molecular-weight organic acids derived from the degradation of petroleum contaminants in groundwater. Geochim. Cosmochim. Acta 58:863–877.[CrossRef][Web of Science]
- Cozzarelli, I.M., B.A. Bekins, M.J. Baedecker, G.R. Aiken, R.P. Eganhouse, and M.E. Tuccillo. 2001. Progression of natural attenuation processes at a crude-oil spill site: I. Geochemical evolution of the plume. J. Contam. Hydrol. 53:369–385.[CrossRef][Web of Science][Medline]
- Daniels, D.J. 1996. Surface penetrating radar. Radar, Sonar, Navig. and Avionics Ser. 6. Inst. of Electrical Eng., London.
- Daniels, D.J. (ed.). 2004. Ground penetrating radar. 2nd ed. Radar, Sonar, Navig. and Avionics Ser. 15. Inst. of Electrical Eng., London.
- Daniels, J.J., R. Roberts, and M. Vendl. 1995. Ground penetrating radar for the detection of liquid contaminants. J. Appl. Geophys. 33:195–207.[CrossRef]
- Darayan, S., C. Liu, L.C. Shen, and D. Shattuck. 1998. Measurement of electrical properties of soil. Geophys. Prospect. 46:477–488.
- DeRyck, S.M., J.D. Redman, and A.P. Annan. 1993. Geophysical monitoring of a controlled kerosene spill. p. 5–19. In Proc. Symp. on the Application of Geophysics to Engineering and Environmental Problems (SAGEEP), San Diego, CA. 18–21 Apr. 1993. Environ. Eng. Geophys. Soc., Denver, CO.
- Endres, A., and R. Knight. 1992. A theoretical treatment of the microscopic fluid distribution on the dielectric properties of partially saturated rocks. Geophys. Prospect. 37:531–551.[CrossRef]
- Fam, M.A., and M.B. Dusseault. 1998. High-frequency complex permittivity of shales (0.02–1.30 GHz). Can. Geotech. J. 35:524–531.[CrossRef]
- Fiori, A., A. Benedetto, and M. Romanelli. 2005. Application of the effective medium approximation for determining water contents through GPR in coarse-grained soil materials. Geophys. Res. Lett. 32:L09404, doi:10.1092/2005GL022555.[CrossRef]
- Giannopoulos, A. 2005. GprMax: A ground penetrating radar simulation tool. Available at www.gprmax.org (verified 24 July 2007). Univ. of Edinburgh, Edinburgh, UK.
- Heimovaara, T.L., E.J.G. de Winter, W.K.P. van Loon, and D.C. Esveld. 1996. Frequency-dependent dielectric permittivity measurements from 0 to 1 GHz: Time-domain reflectometry measurements compared with frequency domain network analyzer measurements. Water Resour. Res. 32:3603–3610.[CrossRef]
- Hu, K., and C.R. Liu. 2000. Theoretical study of the dielectric constant in porous sandstone saturated with hydrocarbon and water. IEEE Trans. Geosci. Remote Sens. 38:1328–1336.[CrossRef]
- Johnson, R.H., and E. Poeter. 2005. Iterative use of the Bruggeman–Hanai–Sen mixing model to determine water saturations in sand. Geophysics 70:K33–K38.[CrossRef]
- Kim, C., J.J. Daniels, E.D. Guy, S.J. Radzevicius, and J. Holt. 2000. Residual hydrocarbons in a water-saturated medium: A detection strategy using ground penetrating radar. Environ. Geosci. 7:169–176.[CrossRef]
- Lopes de Castro, D., and R.M.G.C. Branco. 2003. 4-D ground penetrating radar monitoring of a hydrocarbon leakage site in Fortaleza (Brazil) during its remediation process: A case history. J. Appl. Phys. 54:127–144.
- Lowe, D.F., O.L. Carrol, and C.H. Ward (ed.) 1999. Surfactants and cosolvents for NAPL remediation: A technologies practice manual. Lewis, Boca Raton, FL.
- Marcak, H., and T. Golebiowski. 2006. A stochastic interpretation of GPR data concerning the location of hydrocarbon plumes. Near Surface Geophys. 4:163–176.
- Nicholson, A.M. 1968. Broad-band microwave transmission characteristics from a single measure of the transient response. IEEE Trans. Instrum. Meas. 17:395–402.
- Nicholson, A.M., and G.F. Ross. 1970. Measurement of the intrinsic properties of materials by time domain techniques. IEEE Trans. Instrum. Meas. 19:377–382.
- Olhoeft, G.R. 1992. Geophysical detection of hydrocarbon and organic chemical contamination. p. 587–595. In R.S. Bell (ed.) Proc. Symp. on the Application of Geophysics to Engineering, and Environmental Problems (SAGEEP), Oakbrook, IL. Soc. Explor. Geophys., Golden, CO.
- Olhoeft, G.R., and D.E. Capron. 1993. Laboratory measurements of the radio frequency electrical and magnetic properties of soils near Yuma, Arizona. Open File Rep. 93-701. USGS, Washington, DC.
- Osella, A., M. de la Vega, and E. Lascano. 2002. Characterization of a contaminant plume due to a hydrocarbon spill using geoelectrical methods. J. Environ. Eng. Geophys. 7:78–87.
- Persson, M., and R. Berndtsson. 2002. Measuring nonaqueous phase liquid saturation in soil using time domain reflectometry. Water Resour. Res. 38(5):1064, doi:10.1092/2001WR000523.[CrossRef]
- Reynolds, J.M. 1997. An introduction to applied and environmental geophysics. John Wiley & Sons, Chichester, UK.
- Sauck, W.A., E.A. Atekwana, and M.S. Nash. 1998. High conductivities associated with an LNAPL plume imaged by integrated geophysical techniques. J. Environ. Eng. Geophys. 2:203–212.
- Shang, J.Q., R.K. Rowe, J.A. Umana, and J.W. Scholte. 1999. A complex permittivity measurement system for undisturbed/compacted soils. Geotech. Test. J. 22:168–178.
- Sihvola, A.H. 1989. Self-consistency aspects of dielectric mixing theories. IEEE Trans. Geosci. Remote Sens. 27:403–415.[CrossRef]
- Sihvola, A.H. 1999. Electromagnetic mixing formulas and applications. Electromagnetic Waves Ser. 47. Inst. of Electrical Eng., London.
- Sihvola, A.H. 2000. Mixing rules with complex dielectric coefficients. Subsurf. Sens. Technol. Appl. 1:393–415.[CrossRef]
- Sihvola, A.H., and E. Alanen. 1991. Studies of mixing formulae in the complex plane. IEEE Trans. Geosci. Remote Sens. 29:679–687.[CrossRef]
- Starr, G.C., B. Lowery, E.T. Cooley, and G.L. Hart. 1999. Soil water content determination using network analyzer reflectometry methods. Soil Sci. Soc. Am. J. 63:285–289.[Abstract/Free Full Text]
- Taflove, A. 1995. Computational electrodynamics: The finite-difference time-domain method. Artech House, London.
- Taflove, A. 1998. Advances in computational electrodynamics: The finite-difference time-domain method. Artech House, London.
- Taherian, M.R., W.E. Keyon, and K.A. Safinya. 1990. Measurement of dielectric response of water-saturated rocks. Geophysics 55:1530–1541.[CrossRef][Web of Science]
- Teixeria, F.L., W.C. Chew, M. Straka, M.L. Oristaglio, and T. Wang. 1998. Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous and conductive soils. IEEE Trans. Geosci. Remote Sens. 36:1928–1937.[CrossRef]
- Tsui, F., and S.L. Matthews. 1997. Analytical modelling of the dielectric properties of concrete for subsurface radar applications. Constr. Build. Mater. 11:149–161.[CrossRef]
- Von Hippel, A.R. 1954. Dielectrics and waves. 2nd ed. Artech House, Boston, MA.
- West, J.L., K. Handley, Y. Huang, and M. Polkar. 2003. Radar frequency dielectric dispersion in sandstone: Implications for determination of moisture and clay content. Water Resour. Res. 39(2):1026, doi:10.1029/2001WR000923.[CrossRef]
- Yee, K.S. 1966. Numerical solution of initial boundary value problems involving Maxwell's equation in isotropic media. IEEE Trans. Antennas Propag. 14:302–307.[CrossRef]
This article has been cited by other articles:

|
 |

|
 |
 
S. Lambot, A. Binley, E. Slob, and S. Hubbard
Ground Penetrating Radar in Hydrogeophysics
Vadose Zone J.,
February 25, 2008;
7(1):
137 - 139.
[Full Text]
[PDF]
|
 |
|