See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/236595409

Forward modeling of Ground-Penetrating Radar using digitized outcrop

images and multiple scenarios of water saturation

Article in Water Resources Research · June 2001

DOI: 10.1029/2001WR900015

CITATIONS READS

59 447

4 authors, including:

M. B. Kowalsky Peter Dietrich

Michael B. Kowalsky Consulting Helmholtz-Zentrum für Umweltforschung
103 PUBLICATIONS 3,183 CITATIONS 426 PUBLICATIONS 6,323 CITATIONS

SEE PROFILE SEE PROFILE

Yoram Rubin
University of California, Berkeley
252 PUBLICATIONS 8,651 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Methodologies to measure and characterise fine sediment input to alpine rivers and their effects on health and reproduction of gravel spawning brown trout (Salmo
trutta) View project

AquiferAnalogue View project

All content following this page was uploaded by M. B. Kowalsky on 03 June 2014.

The user has requested enhancement of the downloaded file.

 WATER RESOURCES RESEARCH, VOL. 37, NO. 6, PAGES 1615–1625, JUNE 2001

Forward modeling of ground-penetrating radar data

using digitized outcrop images and multiple
scenarios of water saturation
M. B. Kowalsky1
Department of Civil and Environmental Engineering, University of California, Berkeley, California

P. Dietrich and G. Teutsch

Institute of Applied Geology, University of Tübingen, Tübingen, Germany

Y. Rubin
Department of Civil and Environmental Engineering, University of California, Berkeley, California

Abstract. Simple petrophysical models and a sedimentologically interpreted outcrop

photograph corresponding to the plane of a ground-penetrating radar (GPR) survey are
combined to create models for the simulation of GPR. This makes possible the
comparison of GPR field data, synthetic GPR sections, and a lithology image. On the
basis of this comparison the usefulness of the method for identifying hydrologically
significant lithofacies and the sensitivity of the results to different subsurface conditions
may be investigated. In particular, GPR simulations are performed for an outcrop model
at three states of water saturation: uniformly drained (uniform residual saturation),
nonuniformly saturated, and fully saturated. As predicted by reflection coefficient
calculations, comparison among the synthetic simulations highlights the importance of the
existing pore water distribution in determining the “visibility” of lithologic elements in
GPR sections. Comparisons of the synthetic GPR sections with the field data show overall
agreement, though the occurrence of various observed reflections depends on the presence
and distribution of pore water. Conclusions are also drawn about extending outcrop
analog-derived results to investigations of real (fully saturated) aquifers.

1. Introduction logic parameters, a methodical approach is desirable to assist

in the analysis of GPR data and to evaluate the influence of
Ground-penetrating radar (GPR) surveys are increasingly soil conditions on such data.
used to assist in subsurface characterization. The potential of One tool that can help in this goal is provided by the outcrop
the method for gaining information about subsurface structure analog concept. The introduction of outcrop analogs in the
arises from the apparent correlation between material type and geosciences has offered a glimpse at real subsurface heteroge-
electrical properties, in which contrasts cause reflections of neity and corresponding physical parameters [e.g., Davis et al.,
electromagnetic waves. Successful applications of GPR are as 1997] and has allowed for overall conceptual advancements in
far ranging as agriculture [Freeland et al., 1998], archeology such applications as predicting contaminant transport and the
[Tohge et al., 1998], and hydrological analyses [Greaves et al., effect of absorption kinetics on flow modeling [Klingbeil, 1998]
1996; Rubin et al., 1998; Chen et al., 2001]. While many GPR and using geostatistical methods to build flow models from
investigations have been carried out above the groundwater borehole data [Whittaker and Teutsch, 1999].
table, the sensitivity of GPR to the presence of pore water in Outcrop information has been used to validate geophysical
the unsaturated zone is well known [Asprion, 1998]. In some methods as well. For example, Dietrich et al. [1998] evaluated
cases, partial saturation is advantageous because enhanced tomographic methods using an outcrop model. Rea and Knight
contrasts in electromagnetic parameters result and can even [1998] compared GPR data with a nearby outcrop and went so
allow for the estimation of water content and permeability far as to evaluate the agreement between the geostatistical
[Hubbard et al., 1997]. However, the presence of water can also parameters seen in the outcrop and those seen in the GPR
render poor GPR data quality [Asprion, 1998; Vandenberghe data. Reflections were assumed to occur with changes in ma-
and van Overmeeren, 1999]. Whether for applications aiming to terial type, and since the outcrop was bimodal (it consisted
delineate subsurface structures or aiming to estimate hydro- mainly of two materials), reflections were assumed to corre-
spond to boundaries between the materials. With images of
material boundaries as detected by GPR and as seen in a
1
Also at Earth Sciences Division, Lawrence Berkeley National Lab- digitized photo of the outcrop, variogram analyses were per-
oratory, Berkeley, California.
formed, and the variograms between the respective cases were
Copyright 2001 by the American Geophysical Union. compared. The effect that entrapped water could have on the
Paper number 2001WR900015. results was not addressed.
0043-1397/01/2001WR900015$09.00 Another application of GPR is radar stratigraphy, i.e., the
1615
1616 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES

use of GPR to recognize characteristic radar facies and to 2. Petrophysical Models

correlate them with specific depositional environments [Van-
The material properties, which govern electromagnetic wave
denberghe and van Overmeeren, 1999]. Vandenberghe and van propagation, are magnetic permeability and electrical permit-
Overmeeren [1999] systematically investigated the use of GPR tivity and conductivity. The magnetic permeability is approxi-
for the identification of various types of sedimentary struc- mately constant and equal to that of the free space (␮0) for
tures. GPR surveys performed near, and in some cases only a most shallow subsurface materials (i.e., those containing no
few meters away from, cliff faces (outcrops) allowed for a metals). The electrical permittivity and conductivity for such
comparison of cliff face photographs and GPR sections. Ad- materials are, however, functions of, for example, porosity,
ditionally, some forward simulations were performed using water content, and mineral composition [Schön, 1996]. A com-
relative estimates of the reflectivity between modeled struc- monly used form of the electric permittivity is the dielectric
tures to help interpret the data; some large-scale features in constant k (or relative permittivity), defined as the dielectric
the GPR sections were identified in this manner as diffractions permittivity ␧ of the medium normalized by that of free space
from structural discontinuities such as intersecting channels ␧0:
and interfaces at the bottoms of channels. In this study the
effects of pore water were considered at another site as well, k ⫽ ␧/␧ 0 . (1)
where sandy channels in a braided river deposit were probed For estimating the dielectric constant in the present work a
with GPR. In interpreting the data a disturbance in continuity mixture model [Wharton et al., 1980] will be used. The use of
of the water table reflection was attributed to a region of this model is desirable since it is physically rather than empir-
fine-grained particles in a channel fill situated above the water ically based and since it allows for the permittivity to be easily
table. Though not arrived at through a systematic modeling calculated while varying parameters such as the water satura-
approach, it was speculated that particles in this region have tion and porosity. The petrophysical model may be formulated
higher moisture content and lower velocity than those in the (as by Hubbard et al. [1997]) as
region lying directly below.
In fact, water has been observed to enhance the detectability k ⫽ 关共1 ⫺ ␸ 兲V cl 冑 kcl ⫹ 共1 ⫺ ␸兲共1 ⫺ Vcl兲 冑 ks ⫹ Sw␸ 冑 kw
of some materials. For example, Beres et al. [1999, p. 15] rea- ⫹ 共1 ⫺ Sw兲␸ 冑 ka兴2, (2)
soned that large changes in “porosity and water content at
erosional surfaces and boundaries of the open-framework where V cl is the clay content in the mixture, ␸ is the porosity,
gravel produce the most-continuous reflections and correlate S w is the water saturation (the fraction of the pore space filled
best with outcrop data.” They explained the hydrogeologic with water), and k cl , k s , k w , and k a are the dielectric constants
relevance of these elements by noting that their hydraulic con- of the clay, sand grains, water, and air, respectively. The di-
ductivities are 3 or 4 orders of magnitude higher than the electric constant may then be converted through (1) into the
neighboring units and concluded that these elements can be electrical permittivity, which is needed for forward modeling.
mapped with 3-dimensional (3-D) GPR analyses. This conclu- Additionally, in low loss media the corresponding velocity may
sion highlights the need for a systematic, physically based be estimated by the following relationship:
method to determine which lithologic elements can be delin- v ⫽ 1/ 冑 ␮ 0 ␧. (3)
eated with GPR and under what conditions.
Recently, subsurface excavations accompanied GPR sur- For estimation of the electrical conductivity the empirically
veys, yielding outcrop photographs collocated exactly with based model

冋 册
planes of GPR surveys [Beres et al., 1999; Bayer, 2000]. Sedi- ⫺1
␴ eff 1 a
mentological interpretation of the photographs in this case ⫽ (4)
␴w S wn ␸ m
allows for a direct comparison of GPR data with lithology. If
electrical properties can be adequately estimated, the oppor- is chosen, where ␴ w is the electrical conductivity of the pore
tunity exists as well for forward modeling based on the out- fluid and a, m, and n are empirically determined [Archie, 1942;
crop-derived models. Aside from the validation of GPR field Schön, 1996]. Typical values for various subsurface materials
data, forward modeling can allow for a controlled investigation are given by Schön [1996].
of electromagnetic wave sensitivity to different soil types and Attenuation is governed by the electromagnetic parameters
conditions. and is especially sensitive to water saturation and clay content,
In the present study, a methodology is proposed in which a an increase of either causing a decrease in penetration depth
carefully interpreted and digitized outcrop image allows real- for a GPR wave [Saarenketo, 1998]. Illustrating this point,
istic models to be created for the simulation of GPR. First, Vandenberghe and van Overmeeren [1999] described a field
petrophysical models based on porosity and water saturation survey and attributed limited radar penetration depth to the
are adopted in order to estimate electrical properties for each presence of electrically conductive clayey material and a near-
lithological element in the outcrop. Then, the sensitivity of surface water table. For a discussion on the estimation of
GPR surveys to pore water is investigated through forward actual GPR penetration depth, given material properties, and
simulations with models representing three cases of water sat- actual GPR system performance, the reader is directed to
uration: uniformly drained (uniform residual saturation), non- works such as that by Noon et al. [1998].
uniformly saturated, and fully saturated. After comparing the
synthetic simulations, the field data, and the outcrop image,
some conclusions are drawn about the usefulness of the mod- 3. Case Study (Herten Gravel Quarry)
eling approach and of GPR in identifying subsurface structures In the present study, an outcrop site was chosen where
given various soil conditions. geophysical measurements were taken before excavation, and
 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES 1617

Klingbeil, 1998; Klingbeil et al., 1999] giving, for example, po-

rosity and hydraulic conductivity values and geochemical pa-
rameters. Using the sedimentologically interpreted image, the
spatial distribution of the lithological elements was combined
with representative porosity values (that is, each lithologic unit
is assumed to have uniform properties throughout the entire
model and is assigned a single value) to yield the porosity
distribution shown in Figure 2b.
Though a more complete sedimentological description of
the outcrop image is available from Bayer [2000], the major
zones representing separate sedimentary processes are delin-
eated and labeled from 1 to 6 and described briefly. Zones 1,
4, and 6 are mostly composed of sand- and stone-rich compo-
nent-supported gravel. In zones 1 and 2 some thin sequentially
graded deposits occur with thin and discontinuous open frame-
work layers positioned horizontally in zone 1 and angled to
trough shaped in the right side of zone 2. (The higher porosity
wedge-shaped element in the left side of zone 2, which is
pinched out in the middle of the cross section, is a well-sorted,
well-rounded sand-gravel formation.) On average, the hydrau-
lic conductivity in zone 2 is lower than in the other zones.
Whereas zones 1, 4, and 6 represent typical accretionary ele-
ments, zones 3 and 5 contain mostly cut-and-fill sequences, in
which the highly conductive open framework gravel units oc-
cur. These sequences are formed by the deposition of graded
Figure 1. Location of the Herten gravel quarry field site. material, alternating between sand-gravel mixtures with low
porosity and permeability and open framework gravel with
high porosity and permeability. The delineation of zones is
a detailed photograph of the resulting outcrop was taken after somewhat arbitrary. For example, the boundary between zones
excavation. With representative hydrological parameters for 2 and 3 on the right-hand side of the outcrop photo is not
each lithologic element (measured in the laboratory) the dig- clearly defined.
itized outcrop image and the aforementioned petrophysical Ranging from 6.0 ⫻ 10⫺7 to 1.0 m/s, the range in hydraulic
models are used to construct models for GPR simulation; this conductivity values represented by the various elements is
is performed for three different cases of water saturation. large; the largest hydraulic conductivities correspond to the
open framework gravel. However, it is important to note that
3.1. Description of Field Site and Measurements the relevance of various hydrological elements depends on the
The field site is at a gravel quarry in the city of Herten, goals of a hydrogeological investigation. The connectivity of
situated at the southwest border of Germany (Figure 1). The units with high hydraulic conductivity could be important in
sedimentary deposits in this region were formed in a braided the scenario in which the first arrival of contaminants is im-
river environment and consist mainly of layers of poorly to portant. However, for planning the site remediation strategy of
well-sorted sand and gravel with no silt or clay. a reactive contaminant, for example, it is conceivable that in
During August and September 1999, a series of GPR surveys addition to or instead of the high-conductivity zones, the iden-
accompanied excavation of the gravel quarry at Herten. The tification of materials with high adsorption capacity, such as
excavation was performed to a depth of ⬃9 m (never reaching sand, is sought [Kleineidam et al., 1999]. Whether the use of
the groundwater table) in such a way that the exposed face of geophysical methods improves characterization in these cases
the previously GPR-surveyed region could be photographed; is an issue which may be considered through outcrop modeling.
high-resolution photographs were taken every 1–2 m, yielding However, before such issues can be addressed for GPR meth-
six parallel images in a span of 10 m, the dimensions of each ods, it is necessary to first gain an understanding about the
image being 16 m in length by ⬃7 m in depth. Inspection of the detectability of various hydrological elements for different sub-
parallel images indicates that the variation in the third dimen- surface conditions.
sion (perpendicular to the outcrop slices) is significant but
relatively gradual. It is also worth noting that the advancing 3.2. Water Distribution Scenarios
outcrop surface was somewhat uneven as a result of the un- The effects of entrapped water on GPR are not always
stable poorly consolidated sediments; ensuing distortion of the considered in geophysical surveys conducted in the unsatur-
image may have caused some inaccuracy in the sedimentologi- ated zone. Nevertheless, downward infiltration of surface wa-
cal mapping of the outcrop. ter and fluctuations in the groundwater table leave entrapped
Using sediment size and texture information along with con- pore water, the amount of which and the uniformity of which
sideration of the sedimentological processes as constraints, the is a function of time, mean grain diameter, and pore distribu-
outcrop photographs were then carefully interpreted to yield tion [Bear, 1988].
maps of lithology [Bayer, 2000]. In the present study, one pro- In evaluating the aforementioned field data it is important to
file from the work of Bayer [2000] is chosen, the photograph of consider that it rained some days before the GPR measure-
which is shown in Figure 2a. For each representative lithologi- ments. Retained water was therefore suspected to be present
cal unit, measurements were performed in the laboratory [e.g., in the subsurface during the measurements and to have influ-
1618 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES

Figure 2. Outcrop modeling procedure. (a) A sedimentologically interpreted [Bayer, 2000] and digitized
photograph from an outcrop at the Herten field site. (b) Representative porosity distribution. This image is
used to construct models with various water saturation distributions (S w ) for GPR forward modeling. (c)
Nonuniformly saturated model, obtained by assuming that the open framework gravel (shown in black) is at
residual saturation (S w ⫽ 0.08) and that the remaining materials (shown in white) have additional retained
water (S w ⫽ 0.17). The units on the axes are meters.
 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES 1619

Figure 3. Variation of (a) dielectric constant, (b) velocity, and (c) electric conductivity with porosity for
various values of water saturation S w (in increments of 0.1) as calculated with petrophysical models. Note that
in Figure 3c the lowest curve corresponds to S w ⫽ 0.01 (⫽0); the remaining curves in Figure 3c proceed with
S w ⫽ 0.1, 0.2, 䡠 䡠 䡠 , 1.0. On the basis of the porosity estimates for each lithologic element the values for the
elements used in the models for forward simulation are also plotted (crosses, uniformly drained; circles,
nonuniformly saturated; triangles, fully saturated).

enced the GPR data [Bayer, 2000]. Given the complex se- showed negligible amounts of clay. Therefore, in estimating
quence of lithologic units with highly contrasting permeabili- the electrical permittivity and conductivity the volume of clay
ties, the distribution of water was thought to be highly V cl is set to zero. The individual k values are set to 6.9 for sand
heterogeneous. Since investigating the general response of (the value for quartz as measured in the laboratory by Knoll
GPR to soil conditions is the goal of this study, rather than and Knight [1994]) and to 80 and 1 for water and air, respec-
trying to reproduce exactly and explain every feature of the tively.
Herten GPR field data, a simple approach is used to put forth To illustrate the influence of soil porosity and saturation on
some reasonable water distributions. The three scenarios of the electromagnetic parameters, the dielectric constant and
water distribution chosen for forward modeling are as follows: velocity are calculated as a function of porosity for various
1. The first model contains so-called drained elements at a water saturation values using (1)–(3) and are shown in Figures
uniform residual saturation (S w ⫽ 0.08) and will be referred 3a and 3b. On the basis of the representative porosity mea-
to as the drained model. surements and chosen S w values the distributions of values
2. The second model is of nonuniform saturation with used for the outcrop elements in the three synthetic models are
open framework gravel modeled as drained and at residual shown as well.
saturation (S w ⫽ 0.08) and with the remaining elements Electrical conductivity is calculated with (4) and plotted as a
modeled at higher water content (S w ⫽ 0.17); the motivation function of porosity for varying values of S w , as shown in
for this is based on the relationship between mean grain di- Figure 3c. The parameters a and m are assigned values of 0.88
ameter and retained water [e.g., Bear, 1988]. The simplified and 1.37, respectively (average values for unconsolidated
approach to assigning water distribution is intended to yield a
sand), and n is set equal to 2 [Schön, 1996]. Site-specific mea-
model with plausible contrasts in water saturation rather than
surements could improve the accuracy of these values and
one with the most accurate water distribution possible (a more
insure that these models best represent the materials at the
accurate water distribution might be obtained through flow
Herten site. However, for the basic understanding of the in-
simulations). The distribution of water saturation for this
fluence of water saturation on the visibility of different litho-
model is shown in Figure 2c. The black regions are those
logical units such high accuracy of the parameters is unneces-
modeled as completely drained (at residual saturation), while
sary. The electrical conductivity of the pore fluid is taken to be
those with higher water content (less drained) are shown in
white. The first and second models are intended to represent 0.4 mS/cm, a typical value for the investigated site.
potential conditions in the unsaturated zone. The relative change in electrical conductivity with increasing
3. The third model contains all fully saturated elements water saturation is seen to be much larger than that in velocity.
(S w ⫽ 1) and represents an aquifer below the groundwater The nonlinear response of ␴ to S w is evident in (4), in which S w
table. This case is included to help determine whether conclu- is raised to the power of n. The values for the elements in each
sions regarding the detectability of hydrologically relevant tar- of the three synthetic models are shown as well in Figures 3a
gets with GPR at an unsaturated outcrop site can be extended and 3b.
to GPR surveys in the saturated zone (i.e., in a real aquifer). In considering the plots shown in Figure 3, one expects only
small reflections to occur in the drained model, with no units in
3.3. Calculation of the Electromagnetic Parameters particular causing dominant reflections. Since the dielectric
A systematic approach for estimating the electromagnetic constant and electrical conductivity values slightly increase
parameters on the basis of the porosity values shown in Figure with porosity, one also expects the order of elements to help
2b and the S w distributions (as described in section 3.2) can be determine the relative reflectivity of the elements (that is,
achieved through the use of the petrophysical models shown in interfaces between elements with a larger porosity contrast will
(1)–(4). Sieve analyses of soil samples from sites geologically have larger differences in electromagnetic parameters than
analogous to the site which is modeled in the present study those with a smaller contrast in porosity). In the nonuniformly
1620 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES

Figure 4. Velocity distributions calculated using the introduced petrophysical models and porosity mea-
surements for the (a) uniformly drained, (b) nonuniformly saturated (partially drained), and (c) fully saturated
models. Because of the large decrease in velocity for fully saturated materials, the gray scale for Figure 4c is
different than the scale for Figures 4a and 4b.

saturated model the largest contrast should exist between the 3.4. Reflection Coefficients
drained open framework gravel layers and the undrained ele- To help determine which material interfaces correspond to
ments. Among the undrained elements the contrasts in elec- contrasts in electrical properties high enough to create signif-
tromagnetic parameters are comparatively small. For the fully icant GPR reflections, the reflection coefficients are calculated
saturated model, there should be an overall increase in reflec-
for each synthetic model. As an approximation, normal inci-
tivity between all units as compared to the drained model.
dence of the vertically traveling wave front is assumed in the
Depending on the spatial order of occurrence, some reflections
calculations, although the wave fronts clearly intersect litho-
might dominate in this case as well.
logic elements at oblique angles. A general indication of high-
The velocity distributions resulting from the above petro-
physical considerations are shown for each model in Figure 4. reflectivity regions is nonetheless expected.
Inspection of these distributions shows interfaces between The normal reflection coefficient R n between two materials
lithofacies with significant velocity contrasts that are antici- is calculated by
pated to cause GPR reflections. However, the magnitudes of
reflections are functions of additional factors and are better K2 ⫺ K1
Rn ⫽ , (5)
quantified through the calculation of the reflection coefficients. K2 ⫹ K1
 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES 1621

Figure 5. Reflection coefficient distributions calculated assuming normal incidence for each model and a
frequency of 100 MHz for the (a) uniformly drained, (b) nonuniformly saturated, and (c) fully saturated
models. Four regions are indicated with ovals to highlight similarities and differences.

where K is the complex propagation constant for each material reflection coefficient (as shown by shades of gray) for all mod-
and is a function of ␴, ␧, ␮, and radar wave frequency [e.g., els. Regions 2 and 3 contain open framework gravel elements
Turner and Siggins, 1994]. To account for resolution in a GPR that should cause relatively large reflections only in the second
survey being typically around a quarter wavelength, the reflec- model. Thus strong reflections within the zones containing
tion coefficients were calculated sequentially from the top of open framework gravel are not expected without a contrast in
the model downward, and regions positioned less than a quar- water saturation (as in the second model) between the open
ter wavelength below a reflector were assigned reflection co- framework gravel and the surrounding materials within the
efficient values of zero; that is, it is assumed that no reflections zones. Additionally, region 4 shows the bottom border of the
would be generated which are distinct from the reflection gen- sand-gravel wedge formation and is seen to be reflective in the
erated by the overlying reflector. The average wavelengths are uniformly drained and fully saturated models but not in the
112, 107, and 73 cm for the uniformly drained, nonuniformly nonuniformly saturated model.
saturated, and fully saturated models, respectively. These examples already show the importance of water sat-
The distributions of R n values for each model are shown in uration in terms of the reflectivity of the different lithologic
Figure 5. From these images it is apparent that regions of high units. In addition, they demonstrate that the visibility of a
reflectivity vary between the models. Four numbered regions reflector depends on (1) its reflectivity relative to that of
are indicated with ovals to help anticipate and explain poten- nearby interfaces or nearby regions with contrasting water con-
tial differences in the forward simulations. Region 1 indicates tent (i.e., a “weak” reflector might be “masked” by a nearby
an interface between zones that is expected to contain a high stronger reflector) and (2) the resolution of a survey; increas-
1622 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES

Figure 6. GPR simulation geometry. To simulate the GPR survey done at the Herten Site, a 2-D staggered
grid finite difference code (fourth order in space, second order in time) was developed which requires as input
the electromagnetic wave propagation parameters. A Ricker wavelet source with a center frequency of 100
MHz is used to simulate a surface survey with source and receiver points located ⬃1 m apart and with traces
collected about every quarter-wavelength. The discretization in space and time is set to 5.7 cm and 0.2 ns,
respectively, and the air-ground interface is not modeled. The snapshot of a radiating wave is superposed over
the model space for illustrative purposes.

ing the frequency of a survey increases resolution, but the potential importance of including frequency-dependent behav-
benefit of increased frequency can be offset by the correspond- ior in the simulation of GPR, such as that of the potential
ing increase in attenuation and therefore decreased penetra- bound-water relaxation mechanism occurring in materials at
tion depth [e.g., Noon et al., 1998]. The simulation of surveys at low saturation. They present a relatively straightforward nu-
additional frequencies is left for future work. merical procedure to implement such processes. In the present
Some of the issues raised might be resolved by refining the study, determination of frequency dependence for the electri-
method for calculating reflection coefficient distribution within cal properties of the lithological units in their varying degrees
highly heterogeneous environments. However, not all wave of water saturation is not attempted. Instead, the values for the
phenomena, such as 2-D effects and complex wave interfer- electrical permittivity and conductivity are assumed constant
ence patterns, are easily predicted from reflection coefficient with frequency and are calculated for simplicity from the
images alone. Effects such as these may be further explored petrophysical models shown in (1), (2), and (4) with the mea-
through the simulation of GPR, which is applied in section 3.5. sured porosity values and chosen S w distributions. These cal-
culated parameter fields, along with the value of the magnetic
3.5. Simulation of GPR permeability, which is assumed constant and equal to that of
In the present study, a 2-D staggered grid finite difference free space (␮0), are used as input for the finite difference
time domain solution (second order in time and fourth order in procedure. The effect of frequency-dependent wave phenom-
space) of the electromagnetic wave equation is used to com- ena on such GPR modeling of an outcrop image is left for
pute synthetic waveforms (see Levander [1988] for a descrip- future investigation.
tion of how to implement such finite difference schemes). To To minimize wave reflection by the boundaries back into the
simulate the GPR reflection survey performed at the Herten model space, adsorbing boundaries are implemented [e.g.,
site, the source/receiver configuration corresponding to that of Casper and Kung, 1996]. Although slight reflections from the
the field measurements is used, the source and receiver points boundaries to the left and right of the source remain, they are
being located ⬃1 m apart. The air-ground interface is not subtracted from the simulated waveforms (after being calcu-
modeled to simplify analysis; instead, a Ricker wavelet source lated through an additional simulation). See Figure 6 for a
with a center frequency of 100 MHz (chosen to approximately depiction of the model geometry; a simulated wave field is
match the frequency seen in the field data) is placed in the superposed over the model space as a visual aid.
upper layer of the model. Approximately four traces per wave-
length are simulated along the survey line since finer resolution
is typically not expected with closer spacing. Since the code is 4. Results
2-D, 3-D effects (reflections arriving from reflectors not situ- The GPR field data [from Bayer, 2000] corresponding to the
ated in the modeled 2-D slice) are not modeled and are not outcrop plane modeled along with the simulated GPR sections
anticipated to be significant since, as noted in section 3.1, the are shown in Plates 1a–1d. The average velocities of the
variation in the direction parallel to the modeled profile is drained, the nonuniformly saturated, and the fully saturated
gradual. models were used to convert the synthetic traces, recorded as
Many subsurface materials have been shown to have fre- functions of time, into functions of depth assuming vertical
quency-dependent electrical properties [Turner and Siggins, travel paths. To accomplish the same for the field data, the
1994]. Furthermore, Bergmann et al. [1998] demonstrate the average velocity of the nonuniformly saturated model was as-
 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES 1623

Plate 1. Comparison of (a) field data and (b) uniformly drained, (c) nonuniformly saturated, and (d) fully
saturated simulations. The average velocities of the models were used to convert reflection times to depth for the
simulations, and the average velocity for the partially drained model was used to do the same for the field data.
Region 1 highlights a reflection seen in all models, regardless of saturation. Regions 2 and 3 show reflections due
to the drained open framework gravel dominant in the field data and the nonuniformly saturated model but not
in the others. Region 4 corresponds to a dominant reflection seen in the uniformly drained and fully saturated
simulation (but not in the nonuniformly saturated simulation); a dominant reflection in the field data in this region
is not clearly seen.
1624 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES

sumed. In addition, an energy decay function was used to data is improved when the presence of pore water is included
correct for overall attenuation in both the field data and the in the forward simulations, and (3) the synthetic GPR section
simulated waveforms, and for the simulated waveforms a cor- shown for the fully saturated (real aquifer analog) model dif-
rection for cylindrical divergence was additionally performed. fers substantially from the GPR sections of the unsaturated
Since loss of angled features (which corresponded to real struc- simulation models.
ture) in the GPR images was observed with migration, migra-
tion was not performed on the field or simulated data. Al-
though the resulting reflector images can include reflections 5. Summary and Conclusions
arriving from 2-D paths (paths other than directly down from An approach has been described which allows for a better
the source to a reflector and directly back up to the receiver), understanding of what can be seen with GPR. Using simple
the overall image of the reflections corresponds very well to the petrophysical relationships, porosity estimates, and some water
lithologic units from the outcrop image (compare Plate 1 with distribution scenarios, the calculation of reflection coefficients
Figure 2b). for and the GPR modeling of an outcrop analog yield images
The simulated images also show strong agreement with the that predict main reflections seen in field data. Simulation
reflection coefficient images (Figure 5), the differences be- through a drained model represents what should be seen if the
tween simulations arising where predicted. In some cases, re- assumption of a “dry” analog were correct. Simulation through
flections from lithologic units are seen regardless of water a partially drained model represents what might more realis-
saturation, whereas in other cases the visibility of reflections tically be seen in the unsaturated zone, where pore water
depends on water saturation. Specific examples from the sim- remains entrapped in a heterogeneous distribution. In reality,
ulation results will next be described. A dominant reflection is the extent of entrapped water and the degree of its spatial
seen in region 1 for all simulations (i.e., regardless of satura- heterogeneity are determined by many factors, including the
tion). This results from the contrast in porosity and therefore permeability distribution, which is related to the pore size
in electromagnetic parameters (regardless of saturation) be- distribution and, ultimately, to the sedimentary environment
tween one of the open framework gravel units (with ␸ ⫽ 0.23) responsible for sediment deposition. The simulated GPR sec-
and a matrix-supported unit (with ␸ ⫽ 0.13). However, the tion for the fully saturated model is that which is expected in a
open framework units on the right side of zone 2 (with ␸ ⫽ real aquifer (below the groundwater table). On the basis of the
0.23 or 0.26), for example, are surrounded by a component- comparison of the simulations it is possible to see if structures
supported gravel with similar bulk porosity (␸ ⫽ 0.22); as identified in the unsaturated zone should be similarly identifi-
explained in the reflection coefficients discussion (see above), able in the saturated zone. In fact, the simulated GPR sections
an increased reflection might be expected to occur in this case, obtained for the unsaturated zone (the nonuniformly saturated
with increasing water saturation contrast between the ele- model in particular) and for the fully saturated zone show
ments. This is observed in the simulations; because of the substantial differences. This suggests that further consider-
drained open framework gravel, regions 2 and 3 show (1) ation is required in order to extend conclusions drawn about
strong reflections in the nonuniformly saturated model (where which lithologic units are visible at a relatively shallow outcrop
the surrounding units are at a higher water content) and (2) analog site to those which should be visible much deeper, in a
weak reflections in the remaining simulations, where there is region within the saturated zone, for which no outcrop analog
no contrast in water saturation. Region 4, in contrast, corre- is available for study.
sponds to a dominant reflection seen in the uniformly drained As shown with the synthetic examples in this study, the
simulation and seen, though less easily, in the fully saturated detection of subsurface structures depends on the degree and
simulation, though not clearly seen in the nonuniformly satu- distribution of subsurface water saturation as well as on the
rated simulation. physical properties of the materials present (including porosity
The synthetic GPR sections show overall agreement with the and clay content). In reality, additional factors influence data
field data (Plate 1a) but show, individually, distinct differences. quality such as 3-D effects and unidentifiable sources of noise;
When compared to the uniformly drained simulation, improve- the presented modeling is likely a best-case scenario. However,
ment is seen in the nonuniformly saturated simulation, in it is shown that even with noise-free (simulated) GPR sections,
which the reflection in region 2 is emphasized as it is in the successful detection of subsurface targets depends on the pres-
field data. In the uniformly drained simulation the reflection in ence and distribution of pore water.
region 4 is instead emphasized. Seen in region 3, reflections off As noted, compared to the simulated image for the drained
of the internal features (open framework gravel units) within model, the modeling of retained water improved the agree-
zone 5 are also emphasized in the nonuniformly saturated ment between the synthetic waveforms and the GPR field data
simulation and the field data. However, in the region beneath in some regions. This suggests that water was indeed present in
oval 3 in the field data (the right side of zone 3 in Figure 2b), a heterogeneous distribution and influenced the GPR field
a strong reflection is seen which is not seen in any of the measurements, confirming a hypothesis of Bayer [2000]. It
simulations; the lithological element in this location is nowhere might be further derived from these results that a suitable
else present in the outcrop profile and contains a high propor- improvement for shallow subsurface characterization lies in
tion of fines which are conducive to large amounts of retained the collection of multiple GPR surveys at the same site under
water. Therefore this reflection could be due to a saturation different soil moisture conditions (e.g., during different sea-
contrast that was present in the field but not modeled with the sons or before and after rain storms).
simple assumptions of water saturation used for the present Concerning the preparation of a field survey, the prescribed
work. approach could be very useful if outcrop analogs are available;
Overall consideration of Plate 1 suggests that (1) water sat- analog images might be used together with forward modeling
uration and distribution affect the visibility of individual lith- or reflection coefficient estimation to help predict what sub-
ologic units, (2) modeling of some regions of the GPR field surface conditions are required or if it is even possible to
 KOWALSKY ET AL.: MODELING OF GPR DATA USING OUTCROP IMAGES 1625

adequately delineate hydrologically relevant targets. The scope Dietrich, P., T. Fechner, J. Whittaker, and G. Teutsch, An integrated
of the outcrop analog concept should be further extended with hydrogeophysical approach to subsurface characterization, in
Groundwater Quality: Remediation and Protection, edited by M. Her-
the intent being (1) to determine the hydrologically relevant bert and K. Kovar, IAHS Publ., 250, 513–520, 1998.
features for different sedimentary environments, (2) to deter- Ezzedine, S., Y. Rubin, and J. Chen, Bayesian method for hydrological
mine which of these features may be identified using GPR site characterization using borehole and geophysical survey data:
methods (and how best to use the data to do so), and (3) to Theory and application to the Lawrence Livermore National Labo-
examine more thoroughly the possibility of extending the re- ratory Superfund site, Water Resour. Res., 35(9), 2671–2683, 1999.
Freeland, R. S., R. E. Yoder, and J. T. Ammons, Mapping shallow
sults from an outcrop analog investigation done in the unsat- underground features that influence site-specific agricultural pro-
urated zone to the characterization of a real aquifer. Going duction, J. Appl. Geophys., 40, 19 –27, 1998.
one step beyond, the methodology could be applied to assist in Greaves, R. J., D. P. Lesmes, J. M. Lee, and M. N. Toksoz, Velocity
the estimation of water saturation. For example, in addition to variations and water content estimated from multi-offset, ground-
penetrating radar, Geophysics, 61, 683– 695, 1996.
using GPR for delineating the structure of subsurface lithol- Hubbard, S. S., Y. Rubin, and E. Major, Ground-penetrating-radar-
ogy, the presented methodology allows for an integrated ap- assisted saturation and permeability estimation in bimodal systems,
plication of geophysical and flow modeling using models de- Water Resour. Res., 33(5), 971–990, 1997.
rived from field data and information from an outcrop. Kleineidam, S., H. Ruegner, B. Ligouis, and P. Grathwohl, Organic
In reality, perhaps the delineation of individual lithological matter facies and equilibrium sorption of phenanthrene, Environ.
Sci. Technol., 33, 1637–1644, 1999.
elements is not always possible. In this case, a more practical Klingbeil, R., Outcrop analog studies: Implications for groundwater
use of GPR data might be, first, for the identification of the flow and contaminant transport in heterogeneous glaciofluvial qua-
sedimentary environment and then in combination with mul- ternary deposits, dissertation, pp. 36 – 43, 47–52, Geowiss. Fak.,
tiple types of data through geostatistics-based procedures [e.g., Univ. Tuebingen, Tuebingen, Germany, 1998.
Klingbeil, R., S. Kleineidam, U. Asprion, T. Aigner, and G. Teutsch,
Ezzedine et al., 1999; Chen et al., 2001]. The combination of Relating lithofacies to hydrofacies: Outcrop-based hydrogeological
some GPR information (that with higher confidence) and characterisation of Quarternary gravel deposits, Sediment. Geol.,
borehole information, for example, along with information 129, 299 –310, 1999.
known about the sedimentological environment (from an out- Knoll, M., and R. Knight, Relationships between dielectric and hydro-
crop analog or otherwise), may be an optimal way to incorpo- logic properties of sand-clay mixtures, paper presented at Fifth In-
ternational Conference on Ground Penetrating Radar, Environ. and
rate as much good information as possible, including that be- Eng. Geophys. Soc., Kitchener, Ont., Canada, 1994.
low the resolution of GPR or borehole interpolation alone, Levander, A. R., Fourth-order finite-difference P-SV seismograms,
into site characterization. Geophysics, 53, 1425–1436, 1988.
Noon, D. A., G. F. Stickley, and D. Longstaff, A frequency-
independent characterisation of GPR penetration and resolution
Acknowledgments. The authors thank Evert Slob and the anony- performance, J. Appl. Geophys., 40, 127–137, 1998.
mous reviewers for their thorough reviews of the manuscript and their Rea, J., and R. Knight, Geostatistical analysis of ground-penetrating
insightful comments. This research is part of the special research radar data: A means of describing spatial variation in the subsurface,
program (SFB) 275, TP C3, and was supported by NSF grant EAR Water Resour. Res., 34(3), 329 –339, 1998.
9628306. The first author would also like to acknowledge the support Rubin, Y., S. S. Hubbard, A. Wilson, and M. Cushey, Aquifer charac-
by the DAAD (German Academic Exchange Program). terization, in The Handbook of Groundwater Engineering, edited by
J. W. Delleur, pp. 10.1–10.68, CRC Press, Boca Raton, Fla., 1998.
Saarenketo, T., Electrical properties of water in clay and silty soils,
References J. Appl. Geophys., 40, 73– 88, 1998.
Schön, J. H., Physical Properties of Rocks: Fundamentals and Principles
Archie, G. E., The electrical resistivity log as an aid in determining of Petrophysics, pp. 379 – 478, Pergamon, New York, 1996.
some reservoir characteristics, Trans. Am. Inst. Min. Metall. Pet. Tohge, M., F. Karube, M. Kobayashi, A. Tanaka, and K. Ishii, The use
Eng., 146, 54 – 62, 1942. of ground-penetrating radar to map an ancient village buried by
Asprion, U., Ground-penetrating radar (GPR) analysis in aquifer- volcanic eruptions, J. Appl. Geophys., 40, 49 –58, 1998.
sedimentology: Case studies, with an emphasis on glacial systems of Turner, G., and A. F. Siggins, Constant Q attenuation of subsurface
SW Germany, Tuebinger Geowiss. Arb., A43, Univ. Tuebingen, Tu- radar pulses, Geophysics, 59, 1192–1200, 1994.
ebingen, Germany, 1998. Vandenberghe, J., and R. A. van Overmeeren, Ground-penetrating
Bayer, P., Aquifer-Anolog-Studie in grobklastischen “braided river” radar images of selected fluvial deposits in the Netherlands, Sedi-
Ablagerungen: Sedimentaere/hydrogeologische Wandkartierung ment. Geol., 128, 245–270, 1999.
und Kalibrierung von Georadarmessungen, Diplomkartierung, Wharton, R. P., G. A. Hazen, R. N. Rau, and D. L. Best, Electromag-
Univ. Tuebingen, Tuebingen, Germany, 2000. netic propagation logging: Advances in technique and interpreta-
Bear, J., Dynamics of Fluids in Porous Media, pp. 485– 486, Dover, tion, SPE, 9267, Am. Inst. of Min. Metall. and Pet. Eng., New York,
Mineola, New York, 1988. 1980.
Beres, M., P. Huggenberger, A. Green, and H. Horstmeyer, Using two- Whittaker, J., and G. Teutsch, Numerical simulation of subsurface
and three-dimensional georadar methods to characterize glacioflu- characterization methods: Application to a natural aquifer ana-
vial architecture, Sediment. Geol., 129, 1–24, 1999. logue, Adv. Water Resour., 22(8), 819 – 829, 1999.
Bergmann, T., J. O. A. Robertsson, and K. Holliger, Finite-difference
modeling of electromagnetic wave propagation in dispersive and P. Dietrich and G. Teutsch, Institute of Applied Geology, University
attenuating media, Geophysics, 63, 856 – 867, 1998. of Tuebingen, Sigwartstrausse 10, 72076 Tubingen, Germany.
Casper, D. A., and K.-J. S. Kung, Simulation of ground-penetrating (peter.dietrich@uni-tuebingen.de; georg.teutsch@uni-tuebingen.de)
radar waves in a 2-D soil model, Geophysics, 61, 1034 –1049, 1996. M. B. Kowalsky and Y. Rubin, Department of Civil and Environ-
Chen, J., S. Hubbard, and Y. Rubin, Estimating the hydraulic conduc- mental Engineering, University of California, 435 Davis Hall, Berke-
tivity at the south oyster site from geophysical tomographic data ley, CA 94720. (MBKowalsky@lbl.gov; rubin@ce.berkeley.edu)
using Bayesian techniques based on the normal linear regression
model, Water Resour. Res., in press, 2001.
Davis, J. M., J. L. Wilson, F. M. Phillips, and M. B. Gotkowitz, Rela-
tionship between fluvial bounding surfaces and the permeability (Received August 28, 2000; revised January 10, 2001;
correlation structure, Water Resour. Res., 33(8), 1843–1854, 1997. accepted January 11, 2001.)
 1626

View publication stats