MODELING GPR DATA TO INTERPRET POROSITY AND DNAPL SATURATIONS FOR 

CALIBRATION OF A 3-D MULTIPHASE FLOW SIMULATION 

Open-File Report 02-451 

U.S. Department of the Interior 

U.S. Geological Survey 

This report is preliminary and has not been reviewed for conformity with U.S. Geological 

Survey editorial standards nor with the North American Stratigraphic Code. Any use of 

trade names in this publication is for descriptive purposes only and does not imply 

endorsement by the U.S. Government. 

1 
MODELING GPR DATA TO INTERPRET POROSITY AND 

DNAPL SATURATIONS FOR CALIBRATION OF A 3-D 

MULTIPHASE FLOW SIMULATION 

By Kristen W. Sneddon

1

, Michael H. Powers1, Raymond H. Johnson1,2, and Eileen 

P. Poeter2 

______________________________________________________________________________ 

U.S. GEOLOGICAL SURVEY 

Open-File Report 02-451 

Denver, Colorado 

2002 

1 U.S. Geological Survey, Denver CO

2

 Colorado School of Mines, Golden, CO 

2
U.S. DEPARTMENT OF THE INTERIOR 

BRUCE BABBITT, Secretary 

U.S. GEOLOGICAL SURVEY 

Charles G. Groat, Director 

The use of trade, product, industry, or firm names is for descriptive purposes only and does not imply endorsement 

by the U.S. Government. 

For additional information write to: Copies of this report can be 

can be purchased from: 

U.S. Geological Survey U.S. Geological Survey

Box 25046, Mail Stop 413 Branch of Information Services 

Denver Federal Center Box 25286

Denver, CO 50225-0046 Denver, CO 50225-0425

3
MODELING GPR DATA TO INTERPRET POROSITY AND DNAPL 

SATURATIONS FOR CALIBRATION OF A 3-D MULTIPHASE FLOW 

SIMULATION 

By Kristen W. Sneddon, Michael H. Powers, Raymond H. Johnson, and Eileen P. Poeter 

ABSTRACT 

Dense nonaqueous phase liquids (DNAPLs) are a pervasive and persistent category of 

groundwater contamination. In an effort to better understand their unique subsurface behavior, a 

controlled and carefully monitored injection of PCE (perchloroethylene), a typical DNAPL, was 

performed in conjunction with the University of Waterloo at Canadian Forces Base Borden in 1991. Of 

the various geophysical methods used to monitor the migration of injected PCE, the U.S. Geological 

Survey collected 500-MHz ground penetrating radar (GPR) data. These data are used in determining 

calibration parameters for a multiphase flow simulation. 

GPR data were acquired over time on a fixed two-dimensional surficial grid as the DNAPL was 

injected into the subsurface. Emphasis is on the method of determining DNAPL saturation values from 

this time-lapse GPR data set. Interactive full-waveform GPR modeling of regularized field traces resolves 

relative dielectric permittivity versus depth profiles for pre-injection and later-time data. Modeled values 

are end members in recursive calculations of the Bruggeman-Hanai-Sen (BHS) mixing formula, yielding 

interpreted pre-injection porosity and post-injection DNAPL saturation values. The resulting interpreted 

physical properties of porosity and DNAPL saturation of the Borden test cell, defined on a grid spacing of 

50 cm with 1-cm depth resolution, are used as observations for calibration of a 3-D multiphase flow 

simulation. Calculated values of DNAPL saturation in the subsurface at 14 and 22 hours after the start of 

injection, from both the GPR and the multiphase flow modeling, are interpolated volumetrically and 

presented for visual comparison. 

INTRODUCTION 

Due to unique physical properties of dense nonaqueous phase liquids (DNAPLs), characterization 

and remediation of DNAPL contamination is a challenging problem (Pankow and others 1996). When 

introduced into the subsurface, DNAPLs migrate primarily in response to gravity as a discrete liquid 

phase from groundwater, resulting in a flow pattern more tightly coupled to geological heterogeneity than 

to groundwater flow (Sander and others 1992). DNAPLs have relatively low viscosity, low degradability 

and low solubility in water, so they can move to depth quickly to act as persistant contaminant sources for 

large quantities of groundwater. The small amounts of contaminant that do dissolve (parts per billion) are 

orders of magnitude above EPA-established maximum concentration limits for drinking water. Extensive 

research in recent years has been devoted to improved characterization of DNAPL-contaminated sites and 

understanding of the processes governing DNAPL movement. 

In efforts to better understand and predict complex physical processes, numerical simulation 

models are widely used in the earth sciences. Though it is not possible to truly validate simulation models 

of natural systems (Oreskes and others 1994), it is advantageous to strive for an optimized fit of 

simulation predictions to available data observations in anticipation of evaluating and improving our 

hypotheses regarding DNAPL behavior. Using the first 22 hours of time-lapse 500-MHz GPR data 

acquired during the controlled perchloroethylene (PCE) injection at Canadian Forces Base Borden in 

1991, this paper describes a procedure to interpret GPR data for pre-injection porosity and subsequent 

DNAPL saturations for use as calibration parameters for multiphase flow simulation (Figure 1). 

4 
Closeness-of-fit evaluation of the multiphase flow simulation results to the GPR-derived calibration data 

guides the development of our hypotheses regarding processes governing DNAPL flow in the Borden test 

cell. 

The ability of GPR to detect DNAPL in the subsurface and to be used in calculating saturation 

values relies on the availability of time-lapse data including the pre-injection or pre-release condition. For 

a single time snapshot of GPR data, variations in relative dielectric permittivity may be due to variation in 

pre-existing layer porosity or possibly due to the presence of DNAPL in the subsurface. DNAPL cannot 

be recognized directly in GPR data until it is seen to move through inspection of time-lapse data sets. 

When time-lapse data are available, DNAPL saturation and pre-release porosity can be quantified through 

iterative applications of the Bruggeman-Hanai-Sen (BHS) mixing formula (Bruggeman 1935; Hanai 

1968; Sen and others 1981) applied to the dielectric properties interpreted from electromagnetic modeling 

of the GPR data (Figure 1). These values are then available to estimate hydraulic permeability fields 

through iterative inversion of multiphase flow simulations. 

GEOPHYSICAL DETECTION OF PCE IN THE BORDEN CELL 

The Borden Cell 

A controlled and closely monitored injection of perchloroethylene (PCE), a typical DNAPL, was 

undertaken in the Borden sand aquifer northwest of Toronto in 1991 at the Canadian Forces Base Borden. 

Redman (1992) provides a detailed description of the controlled injection into the 9-m x 9-m x 3.3-m 

deep Borden cell, which was isolated on four sides by double steel sheet piling and underlain by a 

particle-sized clay aquitard (Bauman 1989). The natural aquifer consists of fine to medium grained sand 

with silt lenses, deposited in a progradational foreshore environment (Bohla 1986). The cell was 

considered water-saturated by artificially maintaining the water table at 15 cm beneath the smoothed 

ground surface. Approximately 770 liters of PCE were injected under constant head at 60 cm below 

ground surface over a 70-hour period. 

Geophysical Data 

Many technical papers and theses have been published regarding the extensive geophysical 

monitoring of the Borden injection (a few are Olhoeft 1992; Sander and others 1992; Schneider and 

Greenhouse 1992; Redman 1992; Brewster and others 1992; Brewster and Annan 1994; Sander 1994; 

Brewster and others 1995). Figure 2 (after Brewster and others 1995) illustrates the locations of each 

geophysical method employed throughout the cell. Access tubes AT-1 through AT-8 (equispaced on a 

circle of radius 3 m from the central injection point) and AT-9 (located 1 m from the injection) were used 

for neutron, gamma ray, density, and induction logging (Schneider and Greenhouse 1992), with AT-1 

through AT-8 additionally used for hole-to-hole borehole radar measurements (Sander and others 1992; 

Olhoeft 1992). In situ resistivity probes were employed using access through R-1 and R-2 (Schneider and 

Greenhouse 1992) and borehole time domain reflectometry (TDR) probes were employed using access 

through T-1 and T-2 (Redman and Annan 1992; Brewster and others 1992). Surface resistivity (Schneider 

and Greenhouse 1992) and multiple-frequency surface GPR (Brewster and others 1995; Sander and others 

1992; Sander 1994; Sander and Olhoeft 1994) were also employed at the site, with the grid of GPR 

survey lines shown in Figure 2. The ability of electrical geophysical methods to detect zones of DNAPL 

results from the distinctly contrasting electrical properties of the DNAPL from the groundwater displaced 

in the pore spaces (Brewster and others 1992; Olhoeft 1992; Annan and others 1991; Redman and others 

1991; Redman and Annan 1992). Additionally, neutron logging was expected to respond to the low 

neutron porosity associated with PCE compared to water; density and gamma ray logging were expected 

to respond to the greater density of PCE compared to water (Schneider and Greenhouse 1992). 

5
Fig. 1. Flowc

hart describing the process of determining subsurface porosity and DNAPL saturation from the GPR data using the Brugg

Hanai-Sen (

B

HS) mixin

g formula (

Bruggeman 1935; Hanai 1968; Sen and others 1981) and application of results for starting model an

calibration para

meters for multip

hase flow si

mulatio

n and in

version. 
Fig. 2. Plan view map of 9-m x 9-m experimental Borden cell showing layout of instrumentation 

(after Brewster and others, 1995). Shaded region indicates GPR and multiphase flow model 

domain. 

7
 Grids of 500-MHz surface ground penetrating radar (GPR) data were among the geophysical data 

acquired by the U. S. Geological Survey in conjunction with the University of Waterloo during the PCE 

injection (Sander and Olhoeft 1994) and are used to interpret pre-injection porosity and subsequent 

DNAPL saturations presented here. Sander and others (1992) describe the 500-MHz data acquisition 

using a GSSI SIR-7 radar system over a 340-hour post-injection period. For each suite of data, there were 

a total of 16 radar survey lines, with 8 lines trending north to south (1E to 8E) spaced 1 m from each 

other, and 8 lines trending west to east (1N to 8N) spaced 1 m from each other (Figure 2). 

Also available are a series of pre-injection cores from the eight access tubes and cores obtained at 

1000 hours post-injection from locations CP-1 through CP-8 (Figure 2). Ground truth for the various 

geophysical data was provided from the post-injection cores. Cores CP-1 through CP-4 were evaluated in 

5-cm intervals and cores CP-5 through CP-8 were evaluated in 8-cm intervals for pre-injection porosity 

and 1000 hour post-injection PCE saturation distribution (Brewster and others, 1995). Pre-injection 

porosity was determined from measures of core volume, dry sample weight, and grain density; postinjection PCE saturation was determined from extraction of core PCE into methanol and inferring 

saturation from spectrophotometer and gas chromatography measurements of the methanol. Information 

from each of the interpreted geophysical data sets shown in Figure 2 was later used in constraining the 

GPR modeling parameters of the 500-MHz GPR data presented here. 

500-MHz GPR BORDEN DATA 

Qualitative information on the migration behavior of PCE in the Borden cell may be ascertained 

from visual inspection of the GPR profiles over time and is provided by Sander and others (1992), Sander 

(1994), Brewster and others (1995), and Greenhouse and others (1993). Comparison of pre- to postinjection images and traces indicates which features are pre-existent and which are due to mobilization or 

pooling of DNAPL. To facilitate this comparison, GPR data were regularized to common spatial, time, 

and gain scales. Data were topographically corrected using available survey data (Olhoeft 1998) and 

accurately positioned using cubic spline interpolation between 1-m marks in the raw data. The data were 

then converted to a standard format before binning and smoothing. A bin increment of 1 cm with a bin 

width of 10 cm was applied to the raw data of originally 1-2 traces/cm density resulting in output files of 

one trace per centimeter (i.e., no new data was created in the process). Binned data were de-wowed using 

a 150-MHz low-cut Hanning-window filter. Each file was then zero-time shifted by first removing the 

field range gain, positioning the zero sample at the top of each trace, then applying a common range gain 

across all parcels of data at all elapsed times. 

The GPR acquisition grid over the Borden cell is shown in Figure 2 and a time-series image of 

the regularized data along line 4E is shown in Figure 3 (0, 14, and 22 hours post-injection). The shaded 

area of Figure 2 defines the surface domain of the multiphase flow model and the outlined area of Figure 

3 defines the vertical domain. Line 4E illustrates that accumulation of PCE appears as a bright reflector 

amplified in the GPR record over time. The most prevalent DNAPL feature within the Borden cell is an 

apparently continuous spreading zone at approximately 30-40 nanoseconds (ns) two-way travel time 

(Figure 3), which is equivalent to approximately 1-m depth assuming an average relative dielectric 

permittivity of 25 for the cell (Olhoeft 1992). This zone bows slightly upward near the injection and 

follows an inferred channel to the south of the injection, trending nearly parallel to the southern boundary 

of the cell and dipping toward the eastern boundary. The upward bowing is most likely due to the rate of 

injection exceeding the rate of horizontal flow, rather than due to a pre-injection geometry. 

Aside from the bright reflection horizon indicating the presence of PCE, there are several other 

characteristics of the GPR data not indicative of incoming PCE. These include the diagonal reflections off 

the sheet pile walls, the reflection off a confirmed till layer seen below 80 ns, the direct first-arrival 

energy recorded just after time zero, and various pre-existing layers that do not fill with PCE (Figure 3). 

These details do not change with time and detract from visual interpretation. One way to enhance DNAPL 

features is to difference the data along each trace for two profiles in time. Difference images were used 

only to qualitatively enhance our subsurface view of DNAPL, delineating relevant features to be 

8 
Fig. 3. Regularized 500 MHz GPR time series of line 4E: pre-injection, 14, and 22 hours post-injection. 

Outlined model area indicates GPR and multiphase flow model domain. DNAPL is implied by increased 

amplitudes over time, as seen at approximately 30 ns between 3.5 and 6.5 m. DNAPL can be seen 

entering an inferred channel located at 40 ns two-way travel time between 6.0 and 6.3 m by 22 hours 

post-injection. The bright reflection at approximately 70 ns two-way travel time between 3N and 2N is a 

confirmed till layer, rather than DNAPL accumulation, as its character does not change appreciably with 

time. Diagonal reflections are airwave arrivals from the above-ground portion of the sheet-pile 

containment walls. The reflector at 80 ns two-way travel time extending across the profile is another 

confirmed till layer. 

9
interpreted in GPR modeling. Results of difference images, however, cannot be analyzed on a wavelet 

scale to quantify DNAPL saturation due to slight positioning uncertainties. 

INTERACTIVE FULL-WAVEFORM GPR MODELING 

Background 

For quantitative information on porosity, DNAPL accumulation thickness, and fluid saturation, 

interactive full-waveform GPR modeling was applied to regularized raw waveform data. Sander (1994) 

describes initial attempts at interactive forward GPR modeling, and provides a framework for the research 

continued by Sneddon and others (2000) and developed further here. The one-dimensional (1-D) fullwaveform GPR modeling software used in this research is a module of GRORADAR version 8.99 

(Olhoeft 1998), which is based upon GPRMODV2 by Powers and Olhoeft (1995). Detailed description 

and consideration of the theoretical development and calculation of the GPR forward model is provided 

by Powers (1995). This software requires input of a layered geological model defined by thickness and 

electromagnetic properties of complex dielectric permittivity, electrical conductivity, and complex 

magnetic permeability. Given an earth model, the program computes the associated GPR response by 

calculating the character of the transmitted and reflected waves as the pulse is propagated through a 

geological model of flat, horizontal layers of specified thickness and associated electromagnetic 

properties. Solutions are nonunique, but may be constrained by the user with prior knowledge of the 

medium and/or subsurface geometry. 

Model Parameters 

Table 1 summarizes parameters used in the GPR modeling that were based on prior information 

from previous research of the Borden cell. Information summarized in Table 1 was used as direct input to 

the interactive full-waveform GPR models and in subsequent calculations of pre-injection porosity and 

post-injection DNAPL saturation. GPR models were further constrained by previous studies of DNAPL 

pooling history available from multiple-frequency surface GPR results (Brewster and others 1995; Sander 

and others 1992; Sander 1994), surface and borehole resistivity (Schneider and Greenhouse 1992), 

borehole time domain reflectometry (Brewster and others 1992) at location TDR-1 (Figure 2), and 

neutron logging (Schneider and Greenhouse 1992) at location AT-8 (Figure 2). Hole-to-hole radar data 

available from Olhoeft (1998) were also evaluated for use in this study to constrain pooling history. 

Figure 2 shows which available data lie within the GPR model domain. 

Earth properties are interpreted from the field data by interactively changing the geological model 

and associated electromagnetic properties until the computed GPR response to this model matches the 

field data trace. Because equivalent responses may result from a variety of input models, the interpreter 

must use available information to create a model that is consistent with the likely geology and material 

properties. In the Borden cell, time lapse images (Figure 3) suggest a geologic model of material property 

and thickness variation where visible pooling of PCE is observed between approximately 0.70 and 1.17 m 

(Table 2). 

As a result of the distinct contrast between the relative dielectric permittivities and electrical 

conductivities of water and PCE, the bulk electrical properties of the water-saturated Borden aquifer will 

be altered by PCE injection (Lucius and others 1992). As the more resistive and lower dielectric 

permittivity PCE displaces water in the pore spaces, the bulk conductivity and relative dielectric 

permittivity of the cell are also expected to decrease (Table 1); no change is expected in the magnetic 

properties. To approximate this effect in the interpreted post-injection GPR models, electrical 

conductivity values were made to decrease incrementally in a range of 14.0 to 0.25 mS/m as the modeled 

bulk relative dielectric permittivity decreases from 25 to 9. 

10 
Table 1. Measured and assumed parameters used in GPR modeling 

PARAMETER VALUE REFERENCE 

Dry Borden Sand: 

Relative dielectric permittivity 4.5 Sander and others 1992 

Conductivity (mS/m) < 1 Sander 1994 

Borden Water: 

Relative dielectric permittivity 80 Olhoeft 1992 

Conductivity (mS/m) 50 Redman and others 1991 

Water-saturated Borden Cell: 

Average relative dielectric permittivity 25 Olhoeft 1992 

Modeled relative dielectric permittivity 19 - 27 Current research 

Conductivity (mS/m) 14 Schneider and Greenhouse 1992 

Average porosity (fractional %) 40 Brewster and others 1995 

Detailed porosity (fractional %) 31 - 43 Current research 

PCE: 

Relative dielectric permittivity 2.3 Lucius and others 1992 

Conductivity (mS/m) < 1 Redman 1992 

Water-saturated Borden Sand with PCE: 

Relative dielectric permittivity 3.4 - 29.7 Redman 1992 

Conductivity (mS/m) 0 - 14.5 Redman 1992 

Layering 

Thin, low permeability stringers 

controlling PCE accumulation 

cm-scale 

vertically 

Schneider and Greenhouse 1992 

(borehole resistivity; gamma 

ray, density, and induction 

logging); 

Brewster and others 1992 

(borehole TDR); 

Olhoeft 1992 

(hole-to-hole borehole radar) 

Current research 

(visual inspection of cores) 

Additional Assumptions 

Magnetic permeability (real part) 1 

Magnetic permeability (imaginary) 0 For all materials 

Dielectric permittivity (imaginary) 0 

No frequency dependence 

Horizontal Layering 

Table 2: Conceptual model of pre-injection cell conditions 

Layer 

Dielectric 

Permittivity Depth to Bottom (m) Relative Porosity 

Fractional 

Porosity 

(%) 

I 25 0.63 – 1.00 Average 40 

II 19 – 23 0.70 – 1.05 Less than average 31 – 37 

III 22 – 27 0.74 – 1.10 Greater than average 35 – 41 

IV 20 – 23 0.89 – 1.17 Less than average 33 – 37 

V 23 – 25 2.00 Average 37 – 40

11
Assumptions 

The forward modeling program assumes infinitely horizontal layer boundaries of homogeneous, 

isotropic material in the far-field region of the antenna. The simplified 1-D, flat-layer models will not 

account for amplitude focusing effects resulting from 2-D or 3-D features with non-horizontal boundaries 

(such as a channel), nor for out-of-plane reflections (Olhoeft 1994) from the sidewalls, but it does include 

multiples within a layer. Further assumptions are summarized in Table 

1. 

In addition, the GPR modeler must assume a starting wavelet. 

The forward modeling code computes the response of this wavelet 

originating from the antenna as it interacts with the user-defined 

geoelectric model of the earth. The true starting wavelet is not known 

because it is a property of the antenna and its near-field conditions 

(including ground coupling), properties that vary within each GPR 

profile. In this research, a Ricker wavelet with a center frequency of 

273 MHz (approximating the output of a 500-MHz antenna after 

loading due to ground coupling) was used as the starting pulse for the 

modeling. For the particular 500-MHz antenna used in data acquisition 

at Borden, the Ricker wavelet is expected to resemble the true starting 

pulse (Duke 1990), with polarity of the wavelet dependent on its 

reflection behavior through the material half-space. Though it was not 

the intention of this research to characterize the 500-MHz antenna 

used, similar models calculated with the varying starting wavelets 

(Figure 4) are presented in Figure 5 to demonstrate the effects of the 

wavelet selection on modeling results. For the single trace modeled 

with the modified airwave measured by Duke (1990), a geologically 

reasonable alternative model is produced. This wavelet has not yet 

been used for an extensive remodeling of the cell, but models shown in 

Figure 5 provide an estimate of model uncertainty resulting from 

starting-wavelet uncertainty. A significant result of varying the starting 

wavelet is seen in the decreased number of thin layers necessary to 

create a match to the field trace using the modified air wavelet. To 

match asymmetry in the data, constructive and destructive interference 

is simulated with the introduction of thin layers whereas the modified 

air wavelet has an inherent asymmetry. 

Interpretation 

While the true geologic layering within the cell is complex and 

locally variable, a GPR model that matches the field data is created 

using the simplest geology reasonable. In an effort to elucidate the 

actual zone of PCE accumulation without modeling more layers than 

necessary, each trace to be modeled was first windowed according to 

Fig. 4. Comparison of the 

Ricker starting wavelet 

used in GRORADARTM 

version 8.99 (Olhoeft 

1999) to the modified 

airwave (after Duke 

1990). Polarity of the 

wavelet is dependent on 

its reflection behavior 

through a material halfspace. 

amplitude changes evident in the first 22 hours post-injection. Figure 6 shows an example of the 

windowing process for a field trace located at the intersection of lines 4E and 3N at times 0, 14, and 22 

hours post-injection. The field traces do not appreciably differ for times preceding 30 ns two-way travel 

time and variations noted at times later than 75 ns are due mainly to radio-frequency interference 

(Canadian Forces Base Borden is Canada’s main communications base). In the first 22 hours at this 

location, the main zone of accumulation lies between 32 and 64 ns, with greatest amplitude variation 

between 32 and 44 ns. The domain of the multiphase flow simulation is defined to a depth of 2.5 m 

12
Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

Fig. 5. Comparison of interactive GPR forward modeling results with different starting wavelets. On the left 

are results from use of a Ricker wavelet. On the right are results from use of a modified airwave (after Duke 

1990) measured with the 500-MHz antenna used at Borden.

13
Fig. 6. Demonstration of windowing process using example from field traces at intersection of line 

4E and 3N: pre-injection, 14 hours, and 22 hours post-injection. Little or no change in character of 

the trace before approximately 32 ns indicates little or no DNAPL accumulation in that zone. Clear 

DNAPL accumulation is evidenced in the model window by increases in amplitude and phase shifts. 

Radio-frequency interference at times later than 75 ns precludes modeling in that zone. 

14
(corresponding to approximately 42 ns assuming an average radar velocity of 0.06 m/ns), though there is 

evidence of another DNAPL accumulation around 55 ns. 

The time-lapse images in Figure 3 suggest that a change in geology at approximately 30 to 45 ns 

is the main control on PCE flow in the first 22 hours post-injection. The 1-D field scans at the intersection 

of line 4E and 3N shown in Figure 6 also indicate that the largest changes occur within this same time 

window corresponding to depths between 

approximately 0.9 and 1.4 m (assuming an average 

cell relative dielectric permittivity of 25). Core data 

illustrate the complexity of the geologic layering on a 

centimeter scale, but confirm a broad change occurring 

at depths between approximately 0.6 and 1.5 m 

(Figure 7). The challenge in modeling GPR field data 

for use in multiphase flow simulation is creating a 

simple geologic model that is complex enough to 

represent the main controls on PCE flow, which in this 

case appear to be most significant in the 30 to 45 ns 

time window. 

Within the selected time window, the GPR 

modeler must prioritize matching criteria regarding 

phase and amplitude. In modeling presented here, 

greater priority was assigned to matching phase of 

peaks, troughs, and zero crossings than was assigned 

to matching amplitude. At 500-MHz frequency in a 

medium of average dielectric permittivity 25, 

decreases in relative dielectric permittivity and 

increases in conductivity from PCE infiltration affect 

GPR velocity (evident as phase and zero crossing) 

more than they affect GPR attenuation (evident as 

amplitude variation) within the assumptions of the 

forward modeling program. 

As the interpreter models each trace along a 

GPR profile, care is taken to maintain geologic 

consistency. The 3-D geological model resulting from 

consideration of all the independent, 1-D models 

created on a 0.5-m grid over the test cell approximates 

the subsurface. This 3-D model is related to the true 

subsurface only in that the 1-D flat-layer computed 

responses to the model, given the estimated and 

idealized starting pulse shape, match the acquired field 

data, over a limited window of two-way travel time, at 

every 0.5-m grid point. Any additional information 

used to constrain the geological model, including 

measured values of material properties and geological 

understanding of the depositional environment, 

increases the likelihood that the model is an adequate 

approximation of the true subsurface. 

Though a three-layer model adequately 

described the portion of the cell modeled by Sander 

Broad zone of 

PCE accumulation 

Fig. 7. Pre-injection core from access tube 

AT-8 illustrating cm-scale complexity of 

geologic layering with more broad variation 

between approximately 0.6 to 1.5-m depth 

where PCE tended to accumulate in the 

Borden cell. (Feature at 0.88 m is a break in 

the core.) 

(1994), further modeling throughout the cell indicated the necessity for at least a five-layer conceptual 

model for a good match with broad consistency. Figure 8 shows an example of constructing a five-layer 

model using a top down approach and adding a layer at a time for the field trace windowed in Figure 6. 

15 
Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

Relative real

permittivity

Depth (m)

a b

c d

Fig. 8. Example of building a model from two layers (7a) to three layers (7b) to four layers (7c) to 

five layers (7d) for the point along 4E intersecting with lines 3N at 22 hours post-injection. Depth 

of modeled layers is determined from the windowing process shown in Figure 5.

As expected, the amplitude of the field trace increases with the introduction of PCE. The zone 

preceding 30 ns is held constant at a relative dielectric permittivity of 25, the background for the cell, 

based on the assumption that complex cm-scale layering and material property variation can be averaged 

as a single m-scale layer above and below the main zone of PCE accumulation. Further layers may also 

be indicated in the zone after 45 ns, but were not modeled at this time and were assumed to remain

constant at an average relative permittivity of 25 within the first 22 hours. Using a Ricker starting

wavelet, a minimum of five layers is needed to reasonably model the main zone of accumulation. 

The GPR model of dielectric permittivity versus location, depth, and time, was determined for a 

5-m x 5-m x 2-m subvolume of the test cell with traces modeled at a spacing of 0.5 m. Figure 9 shows

these interpolated models for the pre-injection condition, 14 hours post-injection and 22 hours postinjection. Pre- and post-injection dielectric permittivity values may then be converted to pre-injection 

porosity and post-injection PCE saturation values as described below. 

16
Pre-injection 14 hours 22 hours 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

Relative Dielectric 

P ii i 

Fig. 9. The three columns show sliced views of the same subvolume of the Borden cell as DNAPL was

injected. The color scale represents relative dielectric permittivity variations determined from interactive

forward modeling of the GPR data. The column on the left shows the pre-injection condition, where all

changes are due to variations in porosity. The center column is modeled from GPR data acquired at 14

hours, and the right column is modeled from data acquired at 22 hours. 

17
DETERMINING POROSITY AND DNAPL SATURATION VALUES 

The relationship of porosity to the complex dielectric permittivity of two materials and their 

mixture may be described by the Bruggeman-Hanai-Sen (BHS) mixing formula: 

* * * * C 

Φ = (εmatrix − εcomp ) (ε fluid / εcomp ) 

* * (εmatrix − ε fluid )

where: 

Φ = fractional porosity 

ε*matrix = complex relative dielectric permittivity of the matrix 

ε*fluid = complex relative dielectric permittivity of the fluid 

ε*comp = complex relative dielectric permittivity of the composite mixture 

C = shape factor (1/3 for spherical grains). 

This formula is an effective-medium-theory mixing method, confirmed for use in this application 

by Kutrubes (1986). Other mixing methods are often used, including empirically based formulas (Topp 

and others 1980), and simple volumetric averaging methods (Sihvola 1999). In this case, the BHS 

formula is used because it can be iterated to estimate values of DNAPL saturation in multiphase mixtures 

of matrix, water, and DNAPL (Sneddon and others, 2000). The BHS is a volumetric mixing formula 

describing a physical mixture assuming no interaction between the matrix and the fluid and no scattering 

(e.g., long wavelength compared to pore and particle sizes). These assumptions are valid for the quartz 

silica sand of the Borden aquifer, which does not react with water or PCE. In addition, PCE does not react 

with water. 

The first step in applying the BHS formula to modeled GPR parameters (Figure 1) is to determine 

the porosity of the pre-injection, water-saturated Borden sand for soil horizons of the cell where variations 

in the pre-injection radar data are assumed to be due to variations in porosity. If the sand is assumed to be 

100% water saturated, and a relative dielectric permittivity of 25 is associated with the average measured 

porosity of 40% for the cell (Table 1), then variations in the relative dielectric permittivity of the preinjection volume, attributed to variations in water content, must be the result of changes in porosity. In 

this step, ε*matrix refers to the relative permittivity (approximately 4.5 [Table 1]) of the dry Borden sand 

(ε*s); ε*fluid is the relative permittivity (approximately 80 with a conductivity of 50 mS/m) of the Borden 

water (ε*w); and ε*comp refers to the relative permittivity (variable from approximately 19 to 27) of the 

sand-water composite (ε*s-water) modeled from the pre-injection GPR data at each horizon. Assuming 

spherical matrix particles and full saturation of the aquifer, fractional porosity of the pre-injection watersaturated sand is calculated for every GPR model location (Figure 1, curve A-B). Table 2 summarizes the 

results of these pre-injection fractional porosity calculations, and Figure 10 shows the pre-injection 

porosity distribution in the modeled area of the cell (complete porosity calculation results available in 

Sneddon, 2001). Depth to bottom of a given layer varies according to location within the cell. GPRinterpreted porosity values are not direct input to the flow model, but are used by flow modelers to 

estimate the initial pre-injection permeability field through iterative forward modeling and optimization of 

parameters. 

Once the injection of DNAPL begins, and assuming there is no air in the system, the aquifer 

becomes a three-phase system of sand, water, and DNAPL. Assuming that fluid type, movement, and 

replacement does not change the soil porosity, that sand is always the matrix, and that porosity is always 

fully fluid saturated, another iteration of the BHS is applied to determine ε*comp, which now refers to the 

relative permittivity of the sand-DNAPL composite (ε*s-DNAPL) (Figure 1, curve B-C). At a given porosity 

determined from the pre-injection BHS water-sand curve (Figure 1, curve A-B), ε*s-water and ε*s-DNAPL 

define the end members of ε*matrix and ε*fluid, respectively, for additional iterations of the BHS (Figure 1, 

18 
Fractional P

orosity (%) 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Fig. 10. Modeled variations in porosity of Borden sand within a subvolume of the cell. This model is 

calculated using the BHS formula applied to relative dielectric permittivity values interpreted from interactive 

forward modeling of the pre-injection GPR data. 

19 
curve D-E). The new ε*comp of the sand-water-DNAPL composite (ε*s-water-DNAPL) is determined from 

modeling the GPR response as DNAPL is injected into the cell and displaces water. Variations in the 

values of the post-injection composite ε*s-water-DNAPL are now assumed to be due to the introduction of 

DNAPL by displacement of water rather than to variations in porosity and solving for (1-φ) is then 

equivalent to solving for percent DNAPL saturation (Figure 1, final BHS curve E-D). 

Using the process described above, DNAPL saturations throughout a subvolume of the Borden 

cell were calculated from the GPR data acquired at 14 and 22 hours after the start of injection (complete 

saturation calculation results available in Sneddon, 2001). Vertical profile views of these PCE saturation 

models at the location of line 4E are shown in Figure 11a. This view allows direct comparison with the 

regularized GPR profile images of Figure 3. These GPR-interpreted PCE saturation values are provided 

as calibration data for 3-D, multiphase flow modeling. 

3-D MULTIPHASE FLOW MODELING 

Multiphase flow modeling has its origins in the oil industry (overviews provided by Crichlow 

1977; Peaceman 1977), but has rapidly found use in the field of environmentally related modeling of 

nonaqueous phase liquids. Three-dimensional, multiphase contaminant transport modeling programs have 

been developed over the past decade (Faust 1985; Faust and others 1989; Letniowski and Forsyth 1991; 

Huyakorn and others 1994; and Gerhard 1995) to address this application. Because DNAPLs form longterm and complex contaminant source zones, accurate determination of location and prediction of 

DNAPL movement are critical to successful remediation efforts. 

Calibrating multiphase flow models to match field core data is difficult due to 

uncertainties associated with the properties and continuity of subsurface units that may occur 

between locations of core samples, the collection of which often risks remobilization of DNAPL. 

For this reason, use of GPR data is an attractive alternative because it is noninvasive and 

provides a detailed continuous 3-D view of the subsurface. The dynamics of multiphase flow 

through a matrix of material with heterogeneous porosity, hydraulic permeability, and chemistry 

are, however, difficult to accurately simulate; and the effects of assumptions and simplifications 

are not easily determined, as discussed by Zheng and Bennett (1995, page 375). Initial work 

using the Borden 500-MHz GPR data set to calibrate a 3-D, multiphase flow model is described 

in Sander (1994). The computer code used in that work was abandoned after several attempts to 

execute it led to inconsistent and unsuccessful results. An alternative commercial code was 

considered before we were presented with an opportunity to use an available research code. The 

program 3D2PHASE, described by Gerhard (1995), was offered with support from the author. 

3D2PHASE is an extension of a two-dimensional multiphase flow modeling program developed 

by Kueper and Frind (1991a,b), incorporating capillary-pressure hysteresis effects and 

formulated in terms of wetting phase pressure and saturation. Details of the program can be 

found in Gerhard (1995). 

To optimize hydraulic permeabilities, the forward modeling program 3D2PHASE was 

executed from the universal inversion code (UCODE) developed by Poeter and Hill (1998). 

UCODE uses nonlinear regression to iteratively update the input parameters of the forward 

modeling program to optimize the match between the simulated fluid flow and the GPRinterpreted PCE saturation values at 14 and 22 hours. If the GPR-interpreted saturation values 

are accurate representations of the true subsurface values and the multiphase flow model is an 

appropriate representation of the multiphase flow process, then this inversion would result in 

parameter values approaching true subsurface properties. 

In the inverse problem, distribution of the dependent variable (PCE saturation in this case) is the 

most well known aspect of the system while the distributions of the independent variable (in this case, 

20 
pre-injection porosity and its associated intrinsic permeability) are less well known. Calibration is the 

process of manipulating the independent variable values to best match the simulated distribution of the 

dependent variable (fluid-flow simulated saturations) to the observation data (GPR-interpreted 

saturations). Calibration is an iterative process at multiple levels. The fluid-flow modeler iteratively 

establishes improved model frameworks (e.g., numbers and geometry of layers; formulation of the 

pressure/saturation relationship) within which a nonlinear regression algorithm iteratively varies 

parameter values to estimate the optimal values to minimize the residual between simulated saturations 

and GPR-interpreted saturations. In some instances, the residuals are biased in space or time, or may be 

nonrandomly distributed. Such a situation suggests that the underlying model framework is not 

representative of field conditions; this may be due to errors in the development of the fluid-flow model 

framework, due to assumption violations in use of the fluid-flow simulation code, or due to errors in the 

GPR interpretation. In such a case, the questionable points are re-examined; the GPR data may be 

remodeled for some points, the framework is redefined, and the fluid-flow inversion is repeated. 

Multiphase flow modeling results presented here result from inversion of a multiphase flow simulation 

based on a five layer GPR-interpreted geologic model. Varying fluid-flow geologic models and 

formulation of the pressure/saturation relationship will affect estimated multiphase flow parameters, 

resulting in varied simulated saturation results. 

21 

0 

50 

100 

150 

200 

250 

0 

50 

100 

150 

200 

250 

6N 5N 4N 

6N 5N 4N 3N 2N 

0 

50 

100 

150 

200 

250 

6N 5N 4N 3N 2N 

0 

50 

100 

150 

200 

250 

0 

50 

100 

150 

200 

250 

De

pth (cm) De

pth (cm) De

pth (cm) 

Fractional Porosity (%) DNAPL Saturation (%) 

(a) 

Fig. 11. On the left (11a), a vertical slice along the location of line 4E is shown of the GPR- interpreted pre-injection 

porosity at top; DNAPL saturation at 14 hours (middle) and 22 hours (bottom). On the right (11b), the same vertical 

slice is shown of the simulated multiphase flow model saturation values at 14 hours (top) and 22 hours (bottom). 

(b) 
The final GPR-interpreted porosity model (Figure 10), and the final GPR-interpreted PCE 

saturation models for 14 and 22 hours (Figures 11, 12, and 13), are instrumental in defining an initial 

geologic model and calibration parameters for the multiphase flow simulation (Figure 1). As these 

available values were utilized in this research, they were not used as direct input to the multiphase flow 

model. 3D2PHASE assumes a constant porosity model with varying hydraulic permeability while GPR 

interpretation provides a starting model varying in porosity and does not describe variations in terms of 

hydraulic permeability. The challenge for the fluid flow modeler is to develop a relationship between 

variations in porosity and hydraulic permeability that will result in a simulation match to the saturation 

data. Additionally, the fluid flow model is defined with cell sizes of 15 cm laterally and 5 cm vertically 

while the GPR data was interpreted at spacings of 50 cm laterally and 1 cm vertically. As a result, the 

fluid flow modeler re-interprets the geologic model provided by the GPR interpretation from a 1-cm 

depth scale to a 5-cm depth scale and from a 50-cm to a 15-cm lateral scale, requiring averaging and 

interpolation of the original data. Simulation of the PCE injection, using 3D2PHASE with the optimal 

values determined by the inversion, yields values of PCE saturation for every 15-cm x 5-cm cell at each 

specified time (Figures 11, 12, and 13). 

Limitations 

The inherent limitation of any model simulation of a natural system is that the representation is 

nonunique. Oreskes and others (1994) discuss some of the reasons for this and the implications of this 

premise on considering results of numerical models. Fluid-flow models are nonunique because input 

parameters, such as the number of layers, their geometry, hydraulic permeability, and porosity, are 

incompletely known due to sparsity of data, uncertainty associated with their measurement or 

interpretation (in the case of GPR), and/or loss of resolution due to averaging or interpolation of input 

parameters. We can attempt to minimize some of these limitations by providing GPR-interpreted values 

on a finer lateral grid (and simulating multiphase flow on a finer vertical grid), and by constraining GPR 

modeling parameters with as much other available data as possible. 

RESULTS 

Visual comparison of GPR-interpreted and multiphase flow-simulated saturations requires output 

data to be presented in proportional volumetric formats. Multiphase flow simulation output is defined in 

5-cm x 5-cm cubes within the 5-m x 5-m model domain to a depth of 2.5 m. GPR-interpreted 

permeability, porosity, and saturation models are defined along discrete traces spaced 0.5 m apart within 

the 5-m x 5-m model domain with centimeter depth resolution to 2.5 m. In preparation for visual 

comparison, each set of GPR-interpreted results was first interpolated along profile views with a quadratic 

Kriging routine (Cressie 1991) using an anisotropy ratio and search ellipse reflecting greater layer 

continuity in the horizontal east-west dimension over the north-south or vertical dimensions (Brewster 

and others 1995). Common depth points, in centimeter increments, for all profiles within a grouping of 

GPR-interpreted results were then interpolated between with a similar quadratic Kriging routine assuming 

no horizontal anisotropy among common depth points. On completion of interpolating GPR-interpreted 

data sets between traces and among common depth points, results are defined in 1-cm x 1-cm cubes. The 

5-cm x 5-cm resolution multiphase flow simulation data sets and the 1-cm x 1-cm resolution GPRinterpreted data sets are then each interpolated identically within the volumetric visualization routine. 

Figure 11 compares the GPR-interpreted saturation values and the output of the fluid-flow 

modeling along a vertical profile of the test volume at a location coincident with line 4E of the GPR field 

data. The match between the two saturation profiles could be improved, but is good in terms of general 

rate of migration, extent, and location of PCE accumulation. A volumetric comparison is presented in 

Figures 12 and 13 for 14 and 22 hours post-injection respectively. As seen in the profile view, there are 

strong similarities in general rate of flow, extent, and location of PCE accumulation as well as lateral and 

vertical spreading. The inferred channel shows well in the deepest of the horizontal slices of both models 

22 
GPR-derived Simulated 

saturation saturation 

150 

DNAPL Saturation 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

250 

100 

150 

Dep 200 

th (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

250 

100 

150 

Dep 200 

th (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

Fig. 12. The column on the left shows the GPR- interpreted DNAPL saturation in the modeled subvolume at 14

hours post-injection. The column on the right shows the DNAPL saturation in the same subvolume results of

multiphase flow modeling at 14 hours. 

23
GPR-derived Simulated 

DNAPL Saturation 

saturation saturation

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

Depth (cm) 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

100 

150 

200 

250 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

7E 6E 5E 4E 3E 2N 3N 4N 5N 6N 

0 

50 

00 

150 

200 

250 

1 

Fig. 13. The column on the left shows the GPR- interpreted DNAPL saturation in the modeled subvolume at 

22 hours post-injection. The column on the right shows the DNAPL saturation in the same subvolume results 

of multiphase flow modeling at 22 hours. 

24 
with PCE beginning to fill this linear feature at 22 hours. In the GPR- interpreted model, however, more 

PCE has filled the channel at that time. 

In general, the GPR- interpreted model appears to show PCE spreading laterally and vertically at 

a slightly faster rate than simulated by the multiphase flow model, suggesting that the processes 

governing DNAPL flow and accumulation in the Borden cell seen in the GPR data are not yet fully 

represented by the multiphase flow model. This may be a function of the translation from GPR-available 

parameters and geometry to multiphase flow parameters, or due to limitations or assumption violations of 

the multiphase flow modeling code in this particular application to the Borden cell. Another possibility is 

that due to amplitude focusing in the channel and 3-D effects in the data, 1-D modeling of GPR data does 

not yield sufficiently accurate calibration parameters for the fluid flow modeling to simulate; it does, 

however, provide more detailed information on location, accumulation, and relative saturation of DNAPL 

in the subsurface than is possibly available from core or well data alone. 

DISCUSSION 

All data, modeling, and simulation efforts presented here are attempts to better understand and 

characterize DNAPL behavior in an incompletely known subsurface. There are assumptions and 

limitations inherent in every step of the process from data acquisition through presentation of interpolated 

model results. Improvements to every step of the process are possible. The current effort facilitates 

understanding of the limitations while constraining the assumptions with independently available 

information. The GPR modeling process relies on many assumptions, as does the multiphase fluid flow 

simulation process. In spite of these, GPR interpretation produces a solution that fits the data and provides 

geologic and material property detail that is simplified in the multiphase flow simulation process to result 

in a reasonable resemblance of the simulated flow to GPR observations. 

Efficiently providing a useful GPR-derived calibration data set for successful multiphase flow 

simulation requires a high level of communication between the two disciplines. Ideally, this dialogue 

commences prior to GPR data collection and continues throughout the interpretation process. In the data 

collection phase, a knowledgeable multiphase flow modeler may provide insights into the scale 

(horizontal extent and depth resolution) as well as density of GPR lines desirable to effectively sample the 

suspected zone of contamination and areas favored for calibration of the multiphase flow simulation. 

Upon viewing the collected data, the GPR practitioner may provide additional insights into actual zones 

of significant DNAPL accumulation, which may have been otherwise undersampled in calibration. 

Working together, the multiphase flow modeler and the GPR modeler may determine the most efficient 

density and location of GPR traces to be modeled and interpreted for porosity and DNAPL saturation 

values sufficient for successful calibration of the utilized multiphase flow simulation code. Also working 

together, the flow modeler and the GPR modeler may determine the optimized depth window and 

resolution to be modeled (i.e., it is not efficient for the GPR modeler to resolve 15 fine-scale layers if the 

multiphase simulation code is optimized for 5 coarse-scale layers). Once multiphase flow simulation and 

inversion is performed, examination of residuals and biases in the data may elucidate areas that warrant 

further GPR interpretation. Ideally, the dialogue continues until the process is complete. 

This type of dialogue determined the model domain and led to the simplified five layer 

conceptual model of the Borden cell used in research presented here. In this case, simplification of the 

GPR resolution is both a benefit and a detriment to the outcome of the multiphase flow simulation as it is 

unlikely that the geometric arrangement of earth properties is simple at this scale. In previous DNAPL 

infiltration experiments conducted in the Borden aquifer, intricate interfingering of PCE along zones of 

variable permeability was observed (Poulson and Kueper 1992; Kueper and others 1992). 

This type of dialogue also led to the aspiration for a more robust GPR data interpretation, 

available through existing geophysical technology. A cursory review of the vast literature on the use of 

seismic data to delineate oil fields suggests alternate methods (Kettle 1985; Horvath 1985). One preferred 

method to produce the desired GPR data set for porosity and saturation interpretation consists of: 1) 

carefully processing the data to remove the variable near-field effects on the starting pulse; 2) correcting 

25 
for variable lateral positioning, timing, and gain; and 3) migrating the data in a manner that best preserves 

variations in geology within the final waveforms. For data presented here, steps 1) and 2) have already 

been performed, as has the preliminary velocity analysis required to perform migration. Migration is 

velocity modeling in 2-D or 3-D, similar to what was performed on this data in 1-D. The Borden data set 

requires a full 3-D migration to best interpret the final waveform (Olhoeft 1994). Processing and 

migration is not trivial and must include all geological information from well control and other sources. 

Oil industry case histories, however, indicate that migration results can be accurate enough to extract 

detailed geology from the data (Hardage and others 1994; Connelly and others 1991; Morse and others 

1987; Dahm and Graebner 1982). The information produced by processing and migration could include 

values of porosity and DNAPL saturation on a grid with the density of the GPR data, and not the 

relatively sparse density of the significant interactive forward modeling effort presented here. Moreover, 

migration is the only way to remove focusing effects of the subsurface geometry from the waveform 

amplitudes. Without this, the amplitude variations due to geometrical effects are falsely interpreted as 

material property contrasts. 

Additionally, the dialogue between multiphase modeler and GPR modeler led to the aspiration for 

improvements in the 3-D multiphase flow modeling effort. Though 3D2PHASE remains our best choice, 

in its current form it requires a uniform porosity throughout the model (allowing for variable 

permeability). The GPR interpretation provides a variable porosity model which, due to the constraints of 

3D2PHASE, are not used as direct input to the model, but are used only indirectly to estimate the starting 

model permeability field. Input of variable porosity would produce a different distribution of DNAPL 

saturation for the same volume of DNAPL. The differences in saturation would influence the relative 

permeability and thus the final distribution. As another consideration, 3D2PHASE is based on the 

Brooks-Corey formulation for capillary pressure/saturation relationship (Brooks and Corey 1966). 

Alternative capillary pressure/saturation models have been shown to produce substantially different 

DNAPL spreading patterns (Rathfelder and Abriola 1998) and may produce a better fit to the GPRinterpreted saturations. 

CONCLUSIONS 

Though there is need for improvement of GPR interpretation and multiphase flow modeling for 

practical applications, this research has enhanced our ability to assess and monitor DNAPL contaminated 

sites. One-dimensional full-waveform GPR modeling coupled with recursive applications of the BHS 

mixing formula can quantify porosity and DNAPL saturation. These DNAPL saturations can be used as 

observations to invert multiphase flow models. 

There are advantages in using GPR to detect DNAPLs, but the GPR user must also be aware of 

its limitations. Such limitations include: the need to repetitively collect GPR data during DNAPL 

movement, difficulties in the quantitative interpretation of DNAPL saturation, and the interference of 

adverse site conditions such as high mineralogical clay content with GPR data collection. In spite of 

these, the GPR- interpreted DNAPL saturations are of value to the fluid-flow modelers. While the 

limitations are undesirable, GPR data provide substantially more spatially dense information about 

DNAPL movement and accumulation in the subsurface than alternative invasive methods, without the 

additional risk of contaminant remobilization. 

Combining the skills of GPR practitioners with those of multiphase flow modelers can improve 

the results obtained by both disciplines. As consumers of GPR surveys, flow modelers demand highly 

detailed and accurate subsurface information, challenging GPR geophysicists to improve both data 

acquisition and interpretation. This research contributes to the continuing effort to improve both GPRinterpreted information and multiphase flow algorithms. 

26
ACKNOWLEDGEMENT AND DISCLAIMER 

The U.S. Environmental Protection Agency through interagency agreement #DW14936568-01-0 

to the USGS has funded the research described in this paper, but it has not been subjected to the Agency 

review. Therefore, it does not necessarily reflect the views of the Agency. Mention of trade names or 

commercial products does not constitute endorsement or recommendation for use. 

We would like to acknowledge the work of Gary Olhoeft in initiating this project and helping to 

see it to completion. Louise Pellerin, Kathy Sander, Terrie Rowley, Jeff Lucius, and Vic Labson have also 

played active roles on this project at various times, and we acknowledge and appreciate their efforts. 

We especially acknowledge the enduring support and patience of Aldo Mazzella who made this 

work possible. He is much appreciated. 

REFERENCES 

Annan, A. P., Bauman, P., Greenhouse, J. P., and Redman, J. D., 1991, Geophysics and DNAPLs: 

Ground Water Management No. 5, Fifth National Outdoor Action Conference on Aquifer 

Restoration, Ground Water Monitoring and Geophysical Methods, Nat. Wat. Well. Assn., 963-

977. 

Bauman, P., 1989, A detailed geophysical investigation of a shallow sandy aquifer: M.Sc. thesis, Univ. of 

Waterloo. 

Bohla, J., Jr., 1986, A sedimentological investigation of a progradational foreshore sequence: C.F.B. 

Borden: M.Sc. thesis, Univ. of Waterloo. 

Brewster, M. L., and Annan, A. P., 1994, Ground-penetrating radar monitoring of a controlled DNAPL 

release: 200 MHz radar: Geophysics, 59, 1211-1221. 

Brewster, M. L., Annan, A. P., Greenhouse, J. P., Kueper, B. H., Olhoeft, G. R., Redman, J. D. and 

Sander, K. A., 1995, Observed migration of a controlled DNAPL release by geophysical 

methods: Ground Water, 33, 977-987. 

Brewster, M. L., Redman, J. D., and Annan, A. P., 1992, Monitoring a controlled injection of 

perchloroethylene in a sandy aquifer with ground penetrating radar and time domain 

reflectometry: Symposium on the Application of Geophysics to Environmental and Engineering 

Problems, Env. and Eng. Geophys. Soc., 611-618. 

Brooks, L. J. and A. T. Corey, 1966, Properties of porous media affecting fluid flow. Journal of Irrigation 

and Drainage Division, ASCE, 92, 61-88. 

Bruggeman, D. A. G., 1935, Berechnung Verschiedener Physikalischer Konstanten von Hetarogenen 

Substanzen: Ann. Phys. Lpz., 24, 636-679. 

Connelly, D. L., Ferris, B. J. and Trembly, L. D., 1991, Northwestern Williston Basin case histories with 

3-D seismic data: Geophysics, 56, 1849-1874. 

Cressie, N. A. C., 1991, Statistics for spatial data: John Wiley and Sons, Inc. 

Crichlow, H. B., 1977, Modern reservoir engineering – a simulation approach: Prentice-Hall Inc. 

Dahm, C. G. and Graebner, R. J., 1982, Field development with three-dimensional seismic methods in 

The Gulf of Thailand - A case-history: Geophysics, 47, 149-176. 

Duke, S., 1990, Calibration of ground penetrating radar and calculation of attenuation and dielectric 

permittivity versus depth. M.Sc. thesis, Colorado School of Mines. 

Faust, C. R., 1985, Transport of immiscible fluids within and below the unsaturated zone: a numerical 

model: Water Resources Research, 21, 587-596. 

Faust, C. R., Guswa, J. H., and Mercer, J. W., 1989, Simulation of three-dimensional flow of immiscible 

fluids within and below the unsaturated zone: a numerical model: Water Resources Research, 25, 

2449-2464. 

Gerhard, J. I., B. H. Kueper, and G. R. Hecox, 1998, The Influence of Waterflood Design on the 

Recovery of Mobile DNAPLs. Ground Water, 36(2), 283-292. 

27 
Greenhouse, J., Brewster, M., Schneider, G., Redman, D., Annan, P., Olhoeft, G., Lucius, J. E., Sander, 

K., Mazzella, A., 1993, Geophysics and Solvents: The Borden experiment: The Leading Edge, 

12, No. 4, 261-267. 

Hanai, T., 1968, Electrical properties of emulsions, in Sherman, P., Ed., Emulsion science: Academic 

Press. 

Hardage, B. A., Levey, R. A., Pendleton, V., Simmons, J. and Edson, R., 1994, 3-D seismic case history 

evaluating fluvially deposited thin-bed reservoirs in a gas-producing property: Geophysics, 59, 

1650-1665. 

Horvath, P. S., 1985, The effectiveness of offshore three-dimensional seismic surveys - Case-histories: 

Geophysics, 50, 2411-2430. 

Huyakorn, P. S., Panday, S., and Wu, Y. S., 1994, A three-dimensional multiphase flow model for 

assessing NAPL contamination in porous and fractured media: I. Formulation: Journal of 

Contaminant Hydrology, 16, 109-130. 

Kettle, R. W., 1985, Fifty years of case-histories: Geophysics, 50, 2431-2442. 

Kueper, B. H. and Frind, E. O., 1991a, Two-phase flow in heterogeneous porous media: 1. model 

development: Water Resources Research, 27, 1049-1057. 

Kueper, B. H. and Frind, E. O., 1991b, Two-phase flow in heterogeneous porous media: 2. model 

application: Water Resources Research, 27, 1059-1070. 

Kueper, B. H., Feenstra, S., Rivett, M. O., Cherry, J. A., 1992, A series of controlled field experiments to 

study DNAPL behavior: implications for site remediation: HAZMAT International, Hazardous 

Materials and Environmental Management Conference. 

Kutrubes, D. L., 1986, Dielectric permittivity measurements of soils saturated with hazardous fluids: M. 

Sc. thesis, Colorado School of Mines. 

Letniowski, F. W., and Forsyth, P. A., 1991, A control volume finite element method for three 

dimensional NAPL contamination: Internat. Jour. of Numer. Methods in Fluids, 13, 955-970. 

Lucius, J. E., Olhoeft, G. R., Hill, P. L., and Duke, S. K., 1992, Properties and hazards of 108 selected 

substances – 1992 edition: U. S. Geological Survey Open File Report 92-527. 

Morse, V. C., Johnson, J. H., Crittenden, J. L. and Anderson, T. D., 1987, Techniques for seismic 

delineation of complex structure - A case-history: Geophysics, 52, 802-809. 

Olhoeft, G. R., 1992, Geophysical detection of hydrocarbon and organic chemical contamination: 

Symposium on the Application of Geophysics to Environmental and Engineering Problems, Env. 

and Eng. Geophys. Soc., 587-595. 

Olhoeft, G. R., 1994, Modeling out-of-plane scattering effects: in Proc. of the Fifth Int'l. Conf. on Ground 

Penetrating Radar, 133-144. 

Olhoeft, G. R., 1998, GRORADARTM: Acquisition, processing, modeling, and display of ground 

penetrating radar data: ver. 4.0, software distributed on CD-ROM at Seventh Internat. Conf. on 

Ground Penetrating Radar, (http://www.g-p-r.com). 

Oreskes, N., Shrader-Frechette, K., Belitz, K., 1994, Verification, validation and confirmation of 

numerical models in the earth sciences: Science, 263, 641-646. 

Pankow, J. F., Feenstra, S., Cherry, J. A., Ryan, M. C., 1996, Dense chlorinated solvents in groundwater, 

background and history of the problem in Pankow, J. F., Cherry, J. A., in Eds., Dense Chlorinated 

Solvents and other DNAPLs in Groundwater; History, Behavior, and Remediation: Waterloo 

Press, 1-52. 

Peaceman, D. W., 1977, Fundamentals of numerical reservoir simulation: Elsevier Science Publ. Co., Inc. 

Poeter, E. P., and Hill, M. C., 1998, Documentation of UCODE, a computer code for universal inverse 

modeling: U.S. Geological Survey water-resources investigations report 98-4080. 

Poulsen, M. M., and Kueper, B. H., 1992, A field experiment to study the behaviour of 

tetrachloroethylene in unsaturated porous media: Environmental Science and Technology, 26, 

889-895. 

Powers, M. H., 1995, Dispersive ground penetrating radar modeling in 2D: Ph.D. thesis, Colorado School 

of Mines. 

28 
Powers, M.H. and Olhoeft, G. R., 1995, GPRMODV2: one-dimensional full waveform forward modeling 

of dispersive ground penetrating radar data, version 2.0: U.S. Geological Survey Open File 

Report 95-58, 41p. + floppy diskette.Rathfelder, K., and Abriola, L. M., 1998, The influence of 

capillarity in numerical modeling of organic liquid redistribution in two-phase systems: Advances 

in Water Resources, 21, 159-170. 

Rathfelder, K., and Abriola, L. M., 1998, The influence of capillarity in numerical modeling of 

organic liquid redistribution in two-phase systems, Adv. Water Resour., 21(2), 159-170. 

Redman, J. D., 1992, Geophysics and the Solvents-in-Groundwater Program: Symposium on the 

Application of Geophysics to Environmental and Engineering Problems, Env. and Eng. Geophys. 

Soc., 375-382. 

Redman, J. D., Kueper, B.H., and Annan, A. P., 1991, Dielectric stratigraphy of a DNAPL spill and 

implications for detection with ground penetrating radar: Ground Water Management No. 5, Fifth 

National Outdoor Action Conference on Aquifer Restoration, Ground Water Monitoring and 

Geophysical Methods, Nat. Wat. Well. Assn., 1017-1030. 

Redman, J. D., and Annan, A. P., 1992, Dielectric permittivity monitoring in a sandy aquifer following 

the controlled release of DNAPL: Fourth Internat. Conf. on Ground Penetrating Radar, 191-195. 

Sander, K. A., Olhoeft, G. R., and Lucius, J. E., 1992, Surface and borehole radar monitoring of a 

DNAPL spill in 3D versus frequency, look angle and time: Symposium on the Application of 

Geophysics to Environmental and Engineering Problems, Env. and Eng. Geophys. Soc., 455-469. 

Sander, K. A., 1994, Characterization of DNAPL Movement in Saturated Porous media Using Ground 

Penetrating Radar: M.Sc. thesis, Colorado School of Mines. 

Sander, K. A., and Olhoeft, G. R., 1994, 500-MHz ground penetrating radar data collected during an 

intentional spill of tetrachloroethylene at Canadian Forces Base Borden in 1991: U.S.G.S. Digital 

Data Series DDS-25, CD-ROM. 

Schneider, G. W., and Greenhouse, J. P., 1992, Geophysical detection of perchloroethylene in a sandy 

aquifer using resistivity and nuclear logging techniques: Symposium on the Application of 

Geophysics to Environmental and Engineering Problems, Env. and Eng. Geophys. Soc., 619-628. 

Sen, P. N., Scala, C., and Cohen, M. H., 1981, A self-similar model for sedimentary rocks with 

application to the dielectric constant of fused glass beads: Geophysics, 46, 781-795. 

Sihvola, A. H., 1999, Electromagnetic mixing formulas and applications: IEEE electromagnetic waves 

series 47: Institution of Electrical Engineers, London. 

Sneddon, K. W., Olhoeft, G. R., and Powers, M. H., 2000, Determining and mapping DNAPL saturation 

values from noninvasive GPR measurements: Symposium on the Application of Geophysics to 

Environmental and Engineering Problems, Env. and Eng. Geophys. Soc., 293-303. 

Sneddon, K. W., 2001, Determination of porosity and DNAPL saturation values from modeling of 

noninvasive ground penetrating radar data: M.Sc. thesis, Colorado School of Mines. 

Topp, G. C., Davis, J. L., and Annan, A. P., 1980, Electromagnetic determination of soil water content: 

Measurements in coaxial transmission lines: Water Resources Research 16, 574-582. 

Zheng, C. and Bennett, G. D., 1995, Applied contaminant transport modeling: theory and practice: Van 

Nostrand Reinhold, New York. 

29