l

.

r q
.

.r .
an)»
.‘AQ
. .. 1.;
r572,
.31.“... .n .. .
a .xmnim.
:-\45

tan

f. ..u . A
WAiuunAm.

5%

m- "z"

2'.

 

. q
.mw
v‘mflw...
grant.
a? Wu

._
..

(a

y .1‘\.:
.1}: .. E... r

I. , 0., r‘ .
C .‘II. ' ‘ .
- Iv] 1
r

p .

r! EVI. 2:33.1L

 

 

A ‘H- ‘1)

LIBRARIES
MICHIGAN STATE UNIVERSITY
, EAST LANSING, MICH 48824-1048

This is to certify that the
thesis entitled

QUANTIFYING THE SPECIFIC CONDUCTIVITY OF
CONTAMINATED GROUNDWATER USING GROUND
PENETRATING RADAR AT THE FORMER
WURTSMITH AIR FORCE BASE, OSCODA, MICHIGAN

presented by

JOHN DAVID MOSS

has been accepted towards fulfillment
of the requirements for the

MS. degree in Geological Sciences

 

 

WM gig/«aw?

Mafor Professor’s Signature

/2 /8/04/

Date

MSU is an Afi‘irmative ActiorVEqual Opportunity Institution

.—--.-.—o-.—.-.—.—;.._

-—.._-—-u—-_.-

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

(fiyazé {ii @05

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/ClRC/DateDue.p65-p.15

 

QUANTIFYING THE SPECIFIC CONDUCTIVITY OF CONTAMINATED
GROUNDWATER USING GROUND PENETRATING RADAR AT THE
FORMER WURTSMITH AIR FORCE BASE, OSCODA, MICHIGAN
By

John David Moss

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
MASTER OF SCIENCE

Department of Geological Sciences

2004

ABSTRACT
QUANTIFYING THE SPECIFIC CONDUCTIVITY OF CONTAMINATED
GROUNDWATER USING GROUND PENETRATING RADAR AT THE FORMER
WURTSMITH AIR FORCE BASE, OSCODA, MICHIGAN
By
John David Moss

A calibrated model was developed using ground penetrating radar (GPR) to quantify the
specific conductivity of contaminated groundwater at the former Wurtsmith Air Force
Base in Oscoda, Michigan, USA, which transcends previous GPR use to image and/or
delineate plume geometries based on electromagnetic (EM) wave attenuation and the
observation of a ‘shadow zone’. Data processing techniques consisting of: l) editing, 2)
frequency spectrum analysis, 3) band-pass filtering, 4) automatic gain control application,
5) velocity analysis, 6) normal moveout correction, and 7) spherical divergence
correction are utilized to quantitatively assess the recorded GPR signals and media
parameters as they relate to the propagation of EM waves. The amount of attenuation
recorded in each GPR trace is calculated by assessing the exponential decay of the EM
wave with increasing depth. The change in conductivity of the media between individual
GPR traces is calculated using the derivative of the wave equation solved for attenuation
with respect to conductivity. Empirical relationships between the electrical conductivity
of groundwater and the electrical conductivity of saturated sediments, which are
dependent on porosity and tortuosity, are used to calculate the specific conductivity of
groundwater from the effective conductivity calculated from the GPR data. The model is

calibrated to specific conductivity measurements recorded in groundwater samples

collected from the on-site multi-level wells.

ACKNOWLEDGMENTS

The research is funded by a grant from the National Science Foundation Environmental
Geochemistry and Biogeochemistry Program, EAR-9708487. I would like to thank
David W. Hyndman for his guidance, assistance with all aspects of the project, and
persevering through all obstacles, David T. Long for his guidance, assistance with project
development, and field-sample collection, Phanikumar S. Mantha for his guidance and
assistance with data analysis, Ran Bachrach for his assistance with geophysical
interpretations and data analysis, Jennifer T. McGuire for her assistance with project
development and field-sample collection, and Erik W. Smith for his assistance collecting
GPR data and providing site aquifer parameters. I would like to thank Loretta, Cathy,
and Jackie for their administrative guidance and assistance. I would also like to thank the
National Center for Integrated Bioremediation, Research, and Development (NCIBRD)

and EFX Systems, Incorporated for their cooperation.

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................. vi
LIST OF FIGURES....., ..................................................................................................... vii
KEY TO SYMBOLS OR ABBREVIATIONS .................................................................. ix
INTRODUCTION ............................................................................................................... 1
BACKGROUND ................................................................................................................. 4
SITE HISTORY ............................................................................................................... 6
DATA ACQUISITION ....................................................................................................... 9
Groundwater Sampling Procedures and Frequency .................................................. 12
DATA ANALYSIS ........................................................................................................... 18
EDITING GPR DATA .................................................................................................. l9
FREQUENCY SPECTRUM ......................................................................................... 24
FILTER DATA .............................................................................................................. 27
AUTOMATIC GAIN CONTROL (AGC) .................................................................... 28
VELOCITY ANALYSIS ............................................................................................... 34
Velocity Spectrum ....................................................................................................... 39
Results of Velocity Spectrum Analysis .................................................................... 42
X24" Method .............................................................................................................. 46
Results of X2 4" Analysis ......................................................................................... 47
Interval Velocities ....................................................................................................... 49
CONVERTING TIME TO DEPTH IN A GPR TRANSECT ....................................... SI
Comparison of Calculated Depths to Measured Depths ............................................ 51

TABLE OF CONTENTS (continued)

SPHERICAL DIVERGENCE CORRECTION ............................................................. 57
AMPLITUDE AND ATTENUATION ANALYSIS .................................................... 58
SPECIFIC CONDUCTIVITY CALCULATION .......................................................... 65
CONCLUSIONS ............................................................................................................... 75
APPENDIX. MATLAB® CODES .................................................................................. 78
REFERENCES .............................................. 89

Table I

Table 2

Table 3

Table 4

Table 5

Table 6

Table 7

Table 8

Table 9

Table 10

LIST OF TABLES

Multi-level wells and screened depths below grade ........................................ 12
Pulse Ekko 100 system specifications ............................................................. 14
Input values for velocity spectrum code ......................................................... 43
Calculated NMO velocities ............................................................................. 46
Calculated RMS velocities .............................................................................. 49
Calculated interval velocities .......................................................................... 50
Calculated depths to capillary fringe and clay layer ....................................... 51
Differences between measured and calculated depths to the water table

and capillary fringe ......................................................................................... 52
Differences between measured and calculated depths to the clay layer ......... 57
Weighted Average of Measured specific conductivity values ........................ 74

vi

Figure 1
Figure 2
Figure 3
Figure 4
Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10
Figure 11
Figure 12
Figure 13

Figure 14

Figure 15

Figure 16

Figure 17

LIST OF FIGURES

(Note: Images in this thesis are presented in color)

Study site location ............................................................................................. 7
Multi-level wells at PTA-02 site with surveyed GPR transect ....................... ll
GPR data collected using 50 MHz and 100 MHz antennas ............................ 15
GPR transect at FTA-02 site, raw data ............................................................ 20
CMP gather at ETA-02 site, raw data ............................................................. 21

GPR transect at PTA-02 site, raw data excluding bad trace and data
prior to “time zero” ......................................................................................... 22

CMP gather along transect at FTA-OZ site, raw data excluding data
prior to “time zero” ......................................................................................... 23

Representative power spectrum within GPR transect as a function of

frequency ......................................................................................................... 25
Representative power spectrum within CMP gather as a function of

frequency ......................................................................................................... 26
GPR transect at FTA-02 site, filtered .............................................................. 29
CMP gather at FTA-02 site, filtered ............................................................... 3O
GPR transect at PTA-02 site with AGC applied ............................................. 32
CMP gather at PTA-02 site with AGC applied ............................................... 33
GPR transect at PTA-02 site with AGC applied and contoured May

1999 groundwater specific conductance data ................................................. 35
GPR transect at FTA—02 site with outlined reflectors ..................................... 37

CMP gather centered at 175-meters along GPR transect with outlined
reflectors .......................................................................................................... 38

A) Velocity spectrum. B) Contoured image of the absolute value of
the velocity spectrum ...................................................................................... 4]

vii

Figure 18

Figure 19

Figure 20
Figure 21

Figure 22

Figure 23

Figure 24

Figure 25
Figure 26
Figure 27
Figure 28

Figure 29

LIST OF FIGURES (continued)

Points identifying locations of maximum coherency for positive and

negative phases of water table reflector within the velocity spectrum ........... 44
Points identifying locations of maximum coherency for positive and

negative phases of clay layer reflector within the velocity spectrum ............. 45
XZ-T2 analysis for each reflector ..................................................................... 48
Interpolated ground and water table elevations .............................................. 53

Ground elevations (it) at former Wurtsmith Air Force base (Stark,

Cummings, and Twenter, 1983) ...................................................................... 55
Clay surface elevations ( ft.) at former Wurtsmith Air Force base

(Cummings and Twenter, 1986) ..................................................................... 56
GPR transect at ETA-02 site with gain applied to correct for spherical
divergence ....................................................................................................... 59
Amplitudes of water table and clay layer reflections ...................................... 61
Attenuation of EM wave front ........................................................................ 64
Calculated specific conductivity of groundwater ............................................ 71
Interpolated groundwater specific conductivity values ................................... 72
Calculated and measured specific conductivity of groundwater ..................... 73

viii

a
A

AGC
BTEX
CMP

e

EM
EPA

F

FFT

FT A-02
GPR

k

K
LNAPL
p

m

m/s

NCIBRD

NGVD
NMO
ns

KEY TO SYMBOLS OR ABBREVIATIONS

attenuation coefficient

amplitude

automatic gain control

benzene, toluene, ethylbenzene, xylenes
common mid-point

electrical permitivity
electromagnetic

Environmental Protection Agency
formation factor

fast Fourier transform

Fire Training Area-2

ground penetrating radar

wave number

dielectric constant

light non-aqueous phase liquid
magnetic permeability

meters or cementation factor
meters per second

National Center for Integrated Bioremediation Research and Development

National Geodetic Vertical Datum
normal moveout
nanoseconds

angular frequency

porosity

radius

reduction / oxidation

root mean squared
conductivity

effective conductivity
groundwater conductivity
Siemens

signal-to-noise

time or two-way travel time
micro-gram

United States Geological Survey
velocity

interval velocity

stacking velocity
Wurtsmith Air Force Base
distance or offset

CHAPTER I

INTRODUCTION

In areas where groundwater contamination is pervasive it is Often difficult to
identify the extent of plumes because the installation and sampling of monitoring wells is
both invasive and expensive. In addition, the timeframe required for conventional
contaminant plume analysis (i.e., install and sample monitoring wells and analyze the
collected groundwater samples) is often on the order of weeks to months. The research
presented herein provides a new mechanism to evaluate the nature and extent of
groundwater contamination using remote geophysical methods from the ground surface.
Information regarding the location of plume boundaries is nearly instantaneous and, with
minimal data processing, the contaminant plume conditions such as specific conductivity
can be accurately quantified.

This research was conducted as part of an interdisciplinary research project designed to
better understand reduction/oxidation (redox) processes in a contaminant plume undergoing
natural attenuation based on an interdisciplinary study of hydrogeology, geochemistry, and
microbiology. The goals of the study were to (l) quantitatively assess hydrogeologic,
geochemical, and microbiological constraints that influence redox processes in different regions
of a plume; (2) identify sets of reactions that describe the evolution of groundwater from the
water table through the contaminant plume; and (3) develop dynamic flow and transport models
to explore contaminant degradation and redox processes at the site. The overall research effort
offers an improved understanding of the interactions among physical, chemical, and biological
processes that lead to redox zonation in both pristine and contaminated aquifers. It also provides

useful information for the design of appropriate sampling, monitoring, and remediation strategies.

The analyses, results, and conclusions presented herein will supplement the
larger, interdisciplinary study using a novel approach to assess contaminant plume
conditions (i.e., specific conductivity) with hi gh-resolution ground penetrating radar
(GPR). The approach is based on the theory that contaminant plumes with well
established microbial activity and plumes containing such contaminants as petroleum
hydrocarbons, BTEX compounds (benzene, toluene, ethylbenzene, and xylenes), and
chlorinated compounds, will exhibit groundwater with elevated specific conductivities as
a result of ions leaching from the sediments. Organic or carbonic acids generated during
microbial degradation of hydrocarbons dissolve ions such as calcium, iron, magnesium,
manganese, silica, and others from the aquifer solids, which leach into the groundwater
(Sauck, 2000; Atekwana, 1998; Lendvay et. al., 1998; Sauck, 1998; Nash et al., 1997;
McMahon and Chapelle, 1991). As a result, the highly conductive groundwater
significantly attenuates electromagnetic (EM) waves propagated through the subsurface
by the GPR transmitting antennae, and creates what is referred to as a ‘shadow zone’ or
‘mute zone’ in a GPR profile collected through the contaminant plume (Atekwana et. al.,
2000; Bennejo et. al., 1997).

In addition to the natural occurrence of microbial activity and the subsequent
increase in groundwater specific conductivity, anthropogenic ions, most notably sodium,
have been introduced to the shallow aquifer at the study site between June 1998 and June
1999. A pilot-scale pump-and-treat system designed to remediate groundwater
contamination by passing it through a bioreactor was installed on-site by EFX Systems,
Inc. A large amount of sodium was discharged from the bioreactor as part of the

remediated groundwater returned to the study site, which also contributed to the

increased groundwater specific conductivity identified in down-gradient multi-Ievel
wells.

A calibrated site-specific model is the basis for the unique approach to analyzing
GPR data for the purpose Of quantitatively evaluating the specific conductivity of
groundwater within a contaminant plume. Previous studies have been designed to image
and/or delineate plume geometries based on GPR data rather than quantitatively assessing
hydrogeologic and geochemical properties of groundwater in such complex contaminant
systems. Although the approach is presented with data from the Wurtsmith site, it was
developed for use on any site where GPR signal penetrates the contaminant zone and a
reasonable reflector exists below the plume. The approach can be used to quickly and
accurately assess a specific aspect of the groundwater quality (i.e., specific conductivity)
without installing monitoring wells and collecting groundwater samples for laboratory

analysis.

CHAPTER II

BACKGROUND

Although petroleum hydrocarbons as a pure phase are good insulators and exhibit
low electrical conductivity values, an aquifer contaminated with dissolved phase
petroleum hydrocarbons is often characterized by groundwater with elevated specific
conductivity (Sauck, 2000; Atekwana, 1998; Lendvay et. al., 1998; Nash et al., 1997;
McMahon and Chapelle, 1991). Microbial degradation of hydrocarbons creates organic
or carbonic acids that dissolve ions such as bicarbonate, sulfate, nitrate, iron, manganese,
silica, and others from the aquifer solids. These mobile ions leach into the groundwater
and increase its specific conductivity (Sauck, 2000; Atekwana, 1998; Lendvay et. al.,
1998; Sauck, 1998; Nash et al., 1997; McMahon and Chapelle, 1991). Microbial
degradation of dissolved phase BTEX and chlorinated compounds at the research site
also contributes to the formation of organic or carbonic acids and the subsequent leaching
of ions from the aquifer solids. Organic acids such as benzoic, methylbenzoic,
trimethylbenzoic, toluic, cyclohexanoic, dimethylcyclohexanoic, and carbonic acid have
been identified during studies of petroleum spills and were each produced from aerobic
and anaerobic reactions by microorganism in the subsurface (Sauck, 2000). Bacterial
activity is generally greatest within areas of residual free-phase hydrocarbons near the
base of the vadose zone. It has been suggested that the leachate from this zone is
introduced into the underlying aquifer and/or contaminant plume during recharge events,
which results in an increase in total dissolved solids and a subsequent increase in specific

conductivity (Sauck, 1998).

In addition, anthropogenic sodium was introduced into the shallow aquifer at the
study site between June 1998 and June 1999 during operation of a pilot-scale pump-and-
treat system and bioreactor by EFX Systems, Inc. The sodium discharged from the
bioreactor contributed to the increased groundwater specific conductivity identified in
down-gradient multi-level wells.

The specific conductivity of the groundwater can be measured directly by
analyzing water samples from a well, or, as demonstrated by this current research,
estimated using GPR and subsequent data analysis. As mentioned above, highly
conductive groundwater significantly attenuates EM waves propagating through the
subsurface by the GPR transmitting antennae, and creates a ‘shadow zone’ in a GPR
profile through a contaminant plume (Atekwana et. al., 2000; Bermejo et. al., 1997).
Previous studies have suggested that the attenuation of EM waves is caused by high loss
tangent scattering resulting from small, dispersed concentrations of hydrocarbons in the
capillary fiinge, or from contaminant vapor effects. Others have suggested that the cause
is diffuse reflection resulting from the uneven water table surface due to variations in
pore sizes and capillary forces, petroleum hydrocarbons, and soil mixture (Sauck et. al.,
1998). The current and most recognized theory is that the highly conductive waters
associated with contaminant plumes with well established microbial activity attenuate the
EM waves. The lack of measured reflections is not the result of a subsurface void of
reflectors. Instead, it is caused by a low signal-to-noise (S/N) ratio resulting from energy
scattering and absorption in the highly conductive contaminated groundwater (Yilmaz,
1987). A GPR profile collected along a transect extending from within a contaminant

plume into uncontaminated background waters shows a sharp contrast from a low S/N

ratio to a high S/N ratio, respectively (Bermejo et. al., 1997; Lendvay et. al., 1998; Sauck
et. al., 1998). The use of GPR to locate contaminant plumes Offers several benefits such
as non-intrusive characterization of sites contaminated with hydrocarbons and chlorinated
solvents (Bermejo et. al., 1997), laterally extensive 2-dimensional imaging that can be
acquired with relatively short time in the field, and the ability to delineate plume
boundaries that may change over time due to variations in groundwater flow directions,
which is not possible with stationary monitoring wells. GPR also allows for the

exploration of remote places.

SITE HISTORY

The study site, also referred to as Fire Training Area-2 (PTA-02), is located at the
former Wurtsmith Air Force Base (WAFB) in Oscoda, Michigan, U.S.A., which is
approximately 2 kilometers west of the Lake Huron shoreline in Michigan's Lower
Peninsula (Figure 1). The site was historically used for fire training exercises from 1952
to 1986. During these exercises waste fuel and solvents were ignited in the proximity of
an area currently covered with a concrete pad, and were subsequently extinguished.
Various compounds were used to extinguish the fires, including an aqueous film-forming
foam, a multipurpose dry chemical, potassium bicarbonate based soda, and Halon 121 1.
Incomplete combustion of the waste fuels and solvents resulted in the infiltration of those
liquids into the ground and created a large contaminant plume. In an effort to monitor the
spatial and temporal changes in the distribution of contaminants within the plume
migrating from FT A-02, the National Center for Integrated Bioremediation Research and

Development (NCIBRD) installed a series of multi-level sampling wells. The multi-level

dogwoo— oxm beam A B=wE

 

moc< 33w

 

5:83 .9950

 

 

wells are generally oriented in a line perpendicular to groundwater flow approximately
100 m down gradient from the contaminant source area and were used to collect
groundwater chemistry data as part of this study. The largest identified mass of
contaminants occurs within the capillary fringe, but as a result of leaching, an extensive
plume approximately 50 m wide and 400 m long has developed within the shallow
aquifer. The plume is complex in that it contains both BTEX compounds and chlorinated
solvents (e. g., dichloroethane). Concentrations of BTEX compounds in the plume range
from 20 to 1,000 jig/liter while concentrations of dichloroethane, chloroethane, and vinyl
chloride range from 125 to 98,940 jig/liter. Most of the contamination is associated with

aquifer solids, which have an average total hydrocarbon concentration of 13,650 mg/k g

between 4.5 and 5.7 meters below the land surface (USGS, 1991).

CHAPTER III

DATA ACQUISITION

HISTORICAL AQUIFER DATA

The study site is located within a Shallow sandy aquifer comprised of former
glacial lake sediments. Historical geologic mapping of the region has identified the entire
sequence of geologic units at the site as Mississippian sandstones, siltstones and shales,
overlain by unconsolidated glacial sediments and surficial deposits. The glacial deposits
range in thickness from 30.5-76.2 meters and consist of gravel, sand, silt, and clay. The
surficial deposits include lacustrine sediments (i.e., deltas, beaches, and lakebed sand and
clay), ice-contact sediments such as till, and alluvium, which were deposited about
13,000 years ago when the water within ancestral Lake Huron was at a higher elevation.
Materials from receding glaciers were transported by meltwater and deposited in the area.
Aeolian deposits occur in the eastern part of the former Wurtsmith Air Force Base near
Lake Huron (Stark et. al., 1983; U.S.G.S., 1990). The presence of these surficial deposits
has been verified by on-site well drilling activities and sediment core analyses completed
by members from Michigan State University, NCIBRD, and the United States Geological
Survey (U.S.G.S.).

Previous site characterization efforts included perrneameter tests, slug tests, sieve
analyses, and tracer tests, which demonstrated that the principal aquifer at the Wurtsmith
site is a relatively homogeneous sandy unit approximately 20-meters thick, composed of
sands and gravels that are highly permeable with hydraulic conductivities on the order of

30 meters/day. The homogeneity of the aquifer is later demonstrated by the distinct GPR

reflections at lithologic boundaries and the lack of reflections within homogeneous units.
Penneameter measurements of two core segments collected near multi-level well ML-
307 (Figure 2) yield hydraulic conductivity values of 20 and 41 meters/day. NCIBRD
completed 18 Slug tests on—site and calculated hydraulic conductivity values between 6
and 37 meters/day. They also completed 10 separate Sieve analyses and calculated
hydraulic conductivity values between 15 and 67 meters/day, with an average of 34
meters/day. In addition, a sieve analysis was completed on sediment samples collected
near multi-level wells ML-3 and ML-5 (Figure 2). Using the Hazen method, hydraulic
conductivity values between 15 and 38 meters/day were calculated, with an average Of 26
meters/day. It was noted that there was no correlation between depth and hydraulic
conductivity.

The Shallow aquifer is underlain by at least 30.5 meters of silty clay, which forms
a low permeability aquitard between the surficial aquifer and the underlying deeper
confined aquifer. The water table ranges from approximately 4 to 6 meters below land
surface (183-185 meters above sea level) and varies 0.3 to 0.7 m annually. The water
table gradient ranges from 3 to 5 m/km. Assuming an effective porosity of 30%, the

average groundwater velocity is approximately 0.5 m/day (McGuire et. al., 2000).

GROUNDWATER GEOCHEMICAL DATA
As mentioned above, NCIBRD installed a series of multi-level sampling wells at
this site to document the vertical distribution of contaminants within the on-site
contaminated groundwater plume. The wells were installed as nested groups of 2.5 cm

inner-diameter PVC casings with 0.33 m PVC screens set at various depths within a

10

.8855 Maw c0323 £3, 86 No-<.E 3 £25 ERG—-232 .N can—mi

oVMNMoo omnshoo o

L

3.30:: groom

onnmoo owomhoo coonwoo ovomhoo omonnoo oommhoo owm

F

Koo oomtoo

omwmmr

 

 

 

8:83 920 Mm.
Cameo: moo lIl

 

 

25>) 55232 ._

 

  

 

c2396 RD: I RI
22853320 I I mnI I
I I
new.

. ovwmmr

. 8me

- owwmmr

- 08mm?

. ommmme

, 03mm?

 

08mm?

 

(moron!) BurtIuON

Single borehole. The vertical Spacing between the multi—level screens varied from 0.5 to
2 m. A locked steel standpipe protects the PVC wells at the ground surface. Multi-level
wells ML-3, ML-4, ML-5, ML—7, ML-8, ML-9, and ML-IO are oriented in a line
perpendicular to groundwater flow approximately 100 m down gradient from the primary
contaminant source and were used to collect groundwater chemistry data as part of this

study (Figure 2). Table 1 lists the various screened depths within each multi-level well.

 

 

 

 

 

 

 

 

 

 

Table l. Multi-Ievel wells and screened depths below grade.
Multi—Level Well DCJIII] Below Grade (m)

ML-3 5.43, 5.58, 6.71, 7.63, 8.53, 9.99, 11.46
ML-4 5.33, 5.41, 6.53, 7.45, 8.36

ML-S 5.61, 5.74, 6.90, 7.82, 8.65, 10.30, 11.73
ML-7 6.39, 6.41, 8.22, 9.70, 11.20

ML—8 5.70, 5.84, 6.95, 7.88, 8.78, 10.26, 11.75
ML-9 6.52, 7.61, 8.53, 9.44, 10.91, 12.44
ML-10 6.44, 7.37, 9.75, 11.27

 

 

 

Groundwater Sampling Procedures and Frequency

As part of the overall project, seasonal measurements of hydrogeological,
microbiological, and geochemical parameters were collected at the Site and added to a
dataset started in 1995 by Chapelle and others (Chapelle et. al., 1996). Field personnel
from Michigan State University and/or the USGS collected hydrogeological,
microbiological, and/or geochemical data from the on-site multi-level wells on the
following dates: 1995 (October 24), 1996 (June 13, October 24), 1997 (November 18),
1998 (February 25, April 23, May 20, June 8, August 10, November 12), 1999 (February
19, May 12, June 17, October 6), and 2000 (January 14, June 6, July 10, August 16).
During each sampling event, approximately six well volumes of groundwater was purged

from within each well at a rate of approximately 400 ml/min using a peristaltic pump.

12

Prior to sampling, the pumping rate was decreased and an in-line flow cell was attached
to the tubing to continuously monitor temperature, pH, dissolved oxygen, Eh, and
Specific conductivity. Each parameter value was recorded every 3 to 5 minutes. Once
each parameter stabilized within the flow cell, as defined by the Environmental
Protection Agency’s (EPA’S) low-flow sampling guidance documentation, groundwater
samples were collected (McGuire et. al., 2000). Hydrogen (H2) samples were collected
using the bubble-stripping method (McGuire et. al., 2000, Chapelle and McMahon,
1991). Following stabilization of parameters recorded by the in-line flow cell,
groundwater samples were collected for extensive geochemical analysis including
dissolved sulfide and iron 11 (Fe II) by colorimetric analysis, major cations by flame
atomic adsorption, sulfate (804) by capillary electrophoresis and methane (CH4) by
headspace extraction on a DOC analyzer. When applicable, samples were filtered through
0.45 pm filters. In addition, field blanks were collected during each sampling event
(McGuire et. al., 2000). Groundwater elevation data were collected from three hourly
recording pressure transducers triangulated approximately 300 m apart about the ETA-02
site and individual water level measurements were collected manually from each
available multi-level well (about 20 wells) during each sampling event.

It is important to note that despite the large amount of geochemical data collected
from the Site, the model presented herein only includes specific conductivity data
collected in May 1999 and June 2000. The in-line flow cell was used to measure the
Specific conductivity of the groundwater. The conductivity of a solution, including
groundwater, is highly dependent on temperature. Therefore, conductivity measurements

are often immediately compensated for temperature dependence and reported as specific

13

conductivity so any changes in value regardless of changes in temperature may be noted.
Raw conductivity measurements were compensated to 25°C using the following equation:

Conductivity
+ TC * (T — 25)

 

Specific _ Conductivity(25°C) = 1 (1) (YSI Incorporated)

where TC is the temperature coefficient. A temperature coefficient of 1.91%/°C was
used, which is an approximation for a solution containing pure KCl in water. However,
this temperature coefficient also provides close approximations for calculating the
specific conductivity of groundwater containing other similar salts and ions (YSI

Incorporated).

GPR DATA ACQUISITION
Most GPR data for this study were acquired with a Sensors and Software® Pulse

Ekko 100 system with the Specifications listed in Table 2.

 

 

 

 

 

 

 

Table 2. Pulse Ekko 100 system specifications.

Transmitting and receiving antenna 50 MHz

frequency

Transmitting and receivig voltage 400 volts

Recording scan time 1200 nanoseconds (ns)
Stacking 32-fold

 

 

Data were collected and processed using Pulse Ekko hardware and software and an in-
field laptop personal computer (PC) was used as a graphical interface. A small amount of
data were collected using 100 MHz antennae. Figure 3 shows a comparison of two GPR
transects from the same location using 50 MHz and 100 MHz antennas. Each data set

has been filtered and an automatic gain control has been applied (see below for data

processing details). The two data sets appear very Similar with the following exceptions:

l4

.3235 E2 9: a5 £2 om use. 38:8 5% E0 .m “Ema

  
 

A2225 8:220
our 03 our 8.. 8, om ow cm 0
v
H.
.. 9
.. m
. N w
.. 9
68:5 «no NI: 2: “.2 x
AflouoEv oocaflfi

8N . 03 09. 03. our 8.. . on ow oc on o

" ' a
(spuooas) awn

 

88:3 Ea E: on “.2 x

l) the 100 MHz antennas provide greater resolution of the reflection events, and 2) the
100 MHz antennas do not penetrate as deep into the subsurface due to the preferential
attenuation Of the higher frequencies. Therefore, the 100 MHz antennas were unable to
detect the clay layer at approximately 20 m below the surface. All remaining data
analyses presented herein will utilize data obtained from the 50 MHz antennas.

GPR data were acquired along numerous transects across the FTA-02 site using
the Pulse Ekko 100 equipment. Many transects were oriented in an east-west direction,
which is perpendicular to the groundwater flow direction. Therefore, those GPR
transects provide a cross-sectional “view” of the subsurface, and subsequently a cross-
sectional “view” through the on—site contaminant plume. Additional transects were
oriented north-south, which is parallel to the groundwater flow direction. The data
analysis presented herein focuses on the GPR survey completed along a 200 m transect
located immediately up gradient of the series of multi-level monitoring wells installed by
NCIBRD (Figure 2). The location was selected based on its proximity to the transect of
wells that lie perpendicular to groundwater flow direction and the abundant geochemical
database associated with these wells as a result of temporal groundwater sampling.

During data acquisition along this transect, which is referred to as a common shot
gather, the transmitting and receiving antennas were separated by 2-meters. Between
points of data acquisition (a.k.a. shots) each antenna was moved 0.5-meters in the same
direction along a measuring tape extended along the ground between the start and finish
points of the GPR transect. The 2—meter spacing between antennas was maintained by
tying a 2-meter length of rope between the antennas and drawing the string tight at each

Shot location. In addition, select points along the transect were surveyed to accurately

l6

locate the transect in relation to the fire training area and the on-Site wells and also
establish ground elevations relative to the National Geodetic Vertical Datum of 1929
(NGVD, Figure 2). The antennas were manually moved and data were collected in 0.5-
meter increments along the measuring tape. The antennas were placed in firm contact
with the ground surface at each trace location to maximize coupling, although inherent
variability in surface elevations caused differences in coupling at different points along
the transect.

In addition to completing the ZOO-meter common shot gather, a common mid-
point (CMP) Shot gather was collected to analyze the velocity of the soils. The CMP was
centered at a point corresponding to l75-meters along the ZOO-meter transect, which is
indicated by a yellow star in Figure 2. During completion of the CMP shot gather, the
initial antenna separation was 2-meters and at each successive time step the antenna
offset was increased by 0.5-meters (i.e., each antenna was moved 0.25-meters away from
the center point). Once again, a measuring tape was used for accurate antenna spacing
and the antennas were placed in firm contact with the ground surface at each shot
location to maximize coupling. Analytical methods for interpreting the GPR data and the

subsequent results are presented below.

17

CHAPTER IV

DATA ANALYSIS

Previous studies designed by others to image and/or delineate plume geometries
based on GPR data collected at the PTA-02 site and other sites often only present the
graphical results of the survey. The observation of a ‘shadow zone’ is interpreted as
indicating the presence of contaminated groundwater. However, further work to quantify
the attenuation of the EM wave front and calculate aquifer and/or contaminated
groundwater parameters was not completed as part of the previous studies completed by
others. The data analysis techniques and the model created for quantitatively assessing
the attenuation of the EM wave front as it passes through contaminated groundwater and
calculating the specific conductivity of the groundwater are presented herein.

The GPR data was processed to identify the physical features (e.g. unsaturated
material, saturated material, change from sand to clay) and physical properties (e. g.
velocity of EM wave propagating through the material) of the underlying sediments and
groundwater, with the ultimate goal of estimating the specific conductivity of the
groundwater based on the level of measured signal attenuation. Conventional data
processing sequences as described in Yilmaz (1987 ) and Sheriff and Geldart (1995) were
followed to quantitatively assess individual parameters associated with the data without
losing Significant information. In general, data processing consisted of the following
sequence: 1) editing, 2) frequency spectrum analysis, 3) band-pass filtering to remove
low and high frequency noise, 4) applying an automatic gain control (AGC), 5) velocity

analysis, 6) N MO correction, and 7) correcting for spherical divergence. Information

18

obtained from the conventional data processing techniques was then used to quantify the
amount of signal attenuation and calculate the Specific conductivity of groundwater
within a contaminant plume at the FT A-02 Site.

All processing of acquired GPR data was completed using Matlab®, which has
efficient processing capabilities and data imaging techniques that make it a useful tool for
analyzing acquired GPR data. Matlab® is capable of treating entire GPR data sets as a
single data matrix, which can then be mathematically analyzed using matrix functions,
linear algebra, Fourier transforms, etc. Large data sets are easily managed within
Matlab®. The raw data from the common shot gather and CMP gather displayed using

Matlab® are shown in Figures 4 and 5, respectively (see Appendix).

EDITING GPR DATA

The first step in processing both the common shot gather and CMP data sets was
to remove all data prior to the arrival of the airwave because it contains no useful
information for the data analysis presented herein (Figures 6 and 7). The “time-zero”,
which corresponds to the initiation of data acquisition, was recorded by the Pulse Ekko
software once the GPR survey began. In addition, any “bad” data traces were removed.
More specifically, “bad” data traces are generally created by any number of factors such
as, but not limited to, poor coupling of the antennas to the ground, antennas placed
directly on a foreign object (i.e., rock, metal, etc.), or malfunction of the equipment. Any
“bad” traces generated within the data sets presented herein (1 bad trace out of 476 traces
for the 50 MHz data sets and 13 out of 680 traces for the 100 MHz data set) were

removed by averaging the two adjacent traces (see Appendix). Averaging the adjacent

l9

GPR Transact

O
O
,_

8 8 8 8
V

N o") In
(,,oraz suit}, or .roud papiooai map) sdais auru

20

 

600

120 140 160 180 200
Distance (meters)

100

80

Figure 4. GPR transect at ETA-02 site, raw data.

as in Se 8-5m a 55% ~20 .m an

@295 8:0
mm on an on 9 2 m

o o
o o
v IO

LI.

0
O
0')

("0192 awn“ or Joud papiooai crap) sdars our

 

.6500 n=20

21

.283 05?. 2 corn 8% was 89: can mics—88 See 38 .86 Stat a 8853 ~30 O .o Bani

A2395 oocafifi
OON ow r 00—. ovw ON _. CO r 00 cm ov ON 0

(spuooas) emu

 

88th and s2 x

22

...Eon 0:5: 9 8E Sec mam—ciao Sec 38 .BE NO.<.E “a 8855 wee—e 55¢» #20 .5 Emmi

<-
(spuooas) auILL

t")

 

3:30 n=20

23

data traces preserved the spatial representation of the data that would be lost by simply

deleting a “bad” trace.

FREQUENCY SPECTRUM

A frequency Spectrum (i.e., range of frequency components within the data) for
the GPR data collected from the ETA-02 site was generated and analyzed to visually
identify the very low and very high frequency components of the signal that degrade the
data quality (i.e., noise). Understanding the frequency components of GPR data is
necessary to filter the low- and high-frequency noise. In general, data processing can
often be more easily implemented in the frequency-domain, thus the GPR signals are
analyzed by converting the data from the time-domain to the frequency-domain using a
forward Fast Fourier Transform (FFT). Each GPR data trace collected at a particular
location is a time-dependent signal that can be uniquely and completely described as a
summation of sinusoidal responses to a range of frequencies. Conversely, the synthesis
of frequencies allows us to transform the data from the frequency—domain into the time-
domain, which is accomplished mathematically using the inverse FFT. The significance
of these Operations is that no GPR data is lost or altered by changing between the time-
and frequency-domains (Yilmaz, 1987).

The frequency spectrum of an individual GPR data trace was created for the GPR
transect data and the CMP data set (Figures 8 and 9, respectively) to analyze the

frequency components of the data. The frequency Spectrum was created by transforming

24

xenon—ooh mo cove—Eu a 8 885.: ~30 5.23 8.58% 532. guacamoaom .w 0.5mm"—

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

oor x 3:262“.
m6 m m6 v m.” m ON N m... P m.o o
_ d _l q 1 a d a q q
I o
r _.
T N
ml II n
r L V
r I n
_ _ _ _ _ L _ _ _ L
82803 god :or x
.2 x are:
m we» v 0n n 0N N mé _. ad 0 m6-
J _ a . a _ a . i J v:
I f J N.
I l l 4 o
I : C I N
_ _ _ _ _ p _ _ b _ V

tomcat ”50 Eat 8a: “2.“. .9. x

25

..>o:o=coc No conga a we 65....» E20 5.23 55.50% 830m Basswoaom .o 93mm"—

 

 

 

 

 

 

mor x 5:032”.
0m n we a On n ON N m9 r md o
J1 1 1.1 q _ _ _ _ _ _
o
I _.
I N
I n
I . v
I - m
I I w
_ _ — L L IF Ll _ _ _
E2809. 830a :9 x
.2 x are:
m me c On m 0N N me P nd 0 md-
al Al a a _ _ _ _ _ _ VI
I I NI
I I o
I A N
bl _ _ _ _ _ _ _ _ _ V

 

 

 

553 $20 E0...— 39: Vet .3 x

26

the data trace from the time-domain into the frequency-domain by using an FFT and then
multiplying the trace by its complex conjugate to obtain the real portion of the data (see
Appendix). The frequency Spectrum shows the individual component frequencies that
comprise the entire data trace. The low and high frequency noise was then removed (i.e.,

filtered) from the data, as discussed below.

FILTER DATA

As mentioned above, no information is lost when data is transformed from the
time—domain into the frequency-domain and vice versa using a FFT. Therefore, a data set
can be transformed into the frequency-domain, edited to remove very low and very high
frequency noise that degrades the data quality, and returned to the time—domain. The
resultant data will be identical to the original data with the exception of the applied
frequency modifications. A band-pass filter works in this manner and transforms the data
into the frequency-domain, removes the unwanted frequency components of the signal,
and returns the data to the original time-domain. A band-pass filter is often used because
it can remove the common occurring low- and hi gh-frequency noise present within data,
thus specifying a frequency bandwidth. This process is referred to as defining a boxcar

amplitude spectrum:

A(t)={1,fl<f<t2

A(f) = {0, elsewhere (Yilmaz, 1987)

27

where f = frequency, fl = minimum frequency, and f2 = maximum frequency. For the
common shot gather transect completed at the FTA-02 Site, the low- and high-frequency
cut-off values selected were 19 and 250 MHz respectively, and for the CMP gather the
values were 13 and 250 MHz respectively (see Appendix). For each data set, the low-
and hi gh-frequency cut-off values were manually selected from the frequency spectrums
shown in Figures 8 and 9. In the common Shot gather and CMP gather data sets, “real”
data from reflections and/or direct wave arrival contains frequencies greater than 19 and
13 MHz, respectively. Therefore, lower frequency components are filtered from each
data set. 250 MHz was selected as the high-frequency cut-off value because no
frequency components greater than 250 MHz are present in either data set. Figures 10
and 11 show the common Shot gather and CMP gather data sets, respectively, without the

low- and high-frequency noise that was removed using the appropriate band-pass filter.

AUTOMATIC GAIN CONTROL (AGC)

An automatic gain control (AGC) is a scaling method applied to GPR data in
which the scaling function is commonly derived from the data (Yilmaz, 1987). An AGC
is often applied to GPR data for display purposes because it enhances weak reflection
zones that could otherwise go unseen. However, it is important to note that an AGC-type
gain destroys the original signal character and was thus only used for visualization.

The GPR data collected at the ETA-02 site was visually enhanced using a root mean
squared (RMS) amplitude AGC gain function, which computes the RMS amplitude
within a specified time “window” on a data trace. The gain function squares the

amplitude of each sample within the Specified time “window” and then takes the square

28

8N

 

our

682:. .36 NOISE 3 Boom—ab ”EU .2 BamE

@295 «83.0.5
oer 9... 0N.. 8.. on om

«03:2... ~50

0?

ON

(spuooas) sum

29

623:. .9.“ 8-5... a 55% .36 .: 2%....

AflofiEv «onto
mm on N ON m .. o _. m

V

(‘0

(spuooas) auru

 

.0230 $20 N.ov x

30

root of the mean value, which is considered the RMS amplitude within that time
“window”. The RMS amplitude from each “window” along the trace comprises the gain
function, which is defined for each individual data trace. The RMS amplitude values
from each “window” are then interpolated to create a complete gain function equal in size
to the original data trace. Dividing the original data set by the gain function, which is
equivalent to normalizing the data by the RMS amplitude, completes the application of
the AGC.

It is important to select the appropriate time “window” when an AGC is applied to
GPR data. If the window is too small, each normalization factor used by the AGC (i.e.,
RMS of the data within the window) will be approximately equal to the data within the
window. Therefore, the applied gain (i.e., normalization factor) will be too high within
each window. As a result, variability in signal and reflection amplitudes between each
window will be lost and all reflection amplitudes will appear the same. If the “window”
is too large, the effectiveness of the AGC process is lessened. As the size of the time
“window” increases, the gain function approaches a constant value (i.e., the RMS of the
entire data set). If the applied gain function were a constant value, each data point would
be divided by the same amount and there would be no net gain or loss amongst the data
for visual enhancement.

The common shot gather data set was averaged over a window of 12 traces, while
the CMP gather was averaged over a window of 10 traces (Ergures 12 and 13,
respectively), (see Appendix). The number of traces averaged as part of the AGC was
selected based on visual preferences and the ability to ‘see’ later events in the GPR data.

The AGC was applied to strengthen reflections occurring near the bottom of the GPR

31

  
 

9......“ OO< as, as 8-5... a .822. EU .N. as»...

EofiEV 35%;...

our oi. ON? OS. on on , ov

«03:8... mQO

V

('0

(spuooas) emu

32

CMP Gather

x 10‘7

llw

(spuooas emu

   

33

10

Offset (meters)

Figure 13. CMP gather at PTA-02 site with AGC applied.

transect and create a visually enhanced and revealing data set. Note that the AGC was
applied to the filtered GPR data set. The ‘shadow zone’, corresponding to the location of
the contaminant plume, can be easily seen GPR transect analyzed using an AGC. The
‘shadow zone’ also corresponds to the location of elevated specific conductance
measured May 10-12, 1999 in groundwater from multi-level wells ML-3, ML-4, ML-5,

and ML—8 (Figure 14).

VELOCITY ANALYSIS

The GPR data, as acquired in the field, records the amplitudes of various
reflections as a function of two-way travel time (seconds). In order for this data to be
more applicable for modeling, calculating EM wave attenuation, and calculating the
specific conductivity of groundwater, it is necessary to convert the two-way travel time
between the ground surface and individual reflections into depth. This transformation
into real space allows for a direct comparison between physical measurements of depth to
the water table, depth to the clay layer, thickness of reflective layers, and distances of
wave propagation to those calculated from the GPR data. The basic principles for
analyzing seismic data were used to analyze GPR data and quantitatively assess the
velocity of the EM waves through the unsaturated and saturated sediments.

A quantitative assessment of GPR signal (i.e., electromagnetic waves) attenuation
due to contaminated groundwater must be based on an accurate velocity analysis of the
subsurface media because the amount of attenuation is a function of distance traveled by

the EM wave. The velocities of the EM waves were estimated by analyzing GPR data in

34

See 8:80:38 ocean Emacs...» 8m— »«2 95.9.80 E... 3.3% UO< 5.3 8% No-<.E a 88:8. MED .3 23m."—

Agzv Og‘a
OON 00—. 8.. .Otr ON? 8r .00 8 0* ON 0

    

 

«009.9% mac NOr x

 

\. + + /
2.2%.. 3525.5 .3... 9d,. 32 32
958% 8326365

35

the time-domain. The velocity estimates resulted in proper time-tO-depth conversions,
based on a comparison between GPR derived depths to the water table and clay layer and
those measured in the field. The velocity estimates also result in proper time-to-distance
conversions for the propagating EM wave, which is required to calculate the attenuation
of the wave as it passes through the subsurface media.

The first subsurface layer encountered at the ETA-02 site was the unsaturated
zone. The physical and/or dielectric property differences between the unsaturated and
saturated soils caused a reflection at the capillary fringe above the water table, which was
recorded in both the common shot gather and CMP GPR surveys. It is important to note
that the EM wave reflection off the top of the capillary fringe results in a slight difference
between depth to the actual water table measured in the multi-level wells on-site and the
same depth calculated from the GPR data. Subsequent reflections observed in the GPR
data set are the result of differences in the physical and/or dielectric properties of various
stratigraphic units within the saturated subsurface sediments. The two reflectors at the
PTA-02 site shown in Figure 15 are characterized as: l) the capillary fringe, and 2) a clay
unit at the base of the shallow aquifer.

A general understanding of the different types of velocity that can be estimated is
required prior to performing the velocity analysis. When dealing with a horizontal layer
comprised of a single-velocity media, the curve generated by plotting two-way signal
travel time as a function of transnritter/receiver offset is a hyperbola (Figure 16). The
difference in travel time measured at a given offset and travel time measured at zero-
offset is referred to as normal moveout (NMO). The NMO velocity is therefore defined

as the velocity required to correct for normal moveout and is also equal to the velocity of

36

  
 

2902.2 85:8 5.3 a... 8-5... a 3855 Eu .2 2%.”.

.9295 35%....
89 .01. ON... 8.. . 8 8 ow

«nomad; maO

(spuooas) sum

.93..— >20

28.... .255

37

CMP Gather

Water Table

(spuooas) sun

38

 

Clay Layer

Offset (meters)

Figure 16. CMP gather centered at 175-meters along GPR transect with outlined reflectors.

the media above the reflector. If several horizontal layers are present with distinct
horizontal reflectors at each layer interface, the curve generated by plotting signal travel
time as a function of Offset is still a hyperbola. It is important to note that the NMO
velocity is calculated from the hyperbolic data set using only small offsets relative to
depth. The general rule of thumb is that offsets should be kept less than the depth to a
particular reflector. In addition, when only small offsets are considered, the NMO
velocity for the horizontally layered media is equal to the RMS velocity of the layer
being considered.

Another velocity type that is often calculated is the stacking velocity. The NMO
velocity utilizes small offsets only, while the stacking velocity is equivalent to the
hyperbola that best fits the entire data set. Due to the inherent similarities, NMO velocity
and stacking velocity are Often considered equivalent. Once again, this holds true at
small offsets. The interval velocity, which is also referred to as the Dix velocity, is the

average velocity within a layer or stratigraphic unit between two reflectors.

Velocity Spectrum

A common technique for analyzing any of the above-mentioned velocity types is
to compute a velocity spectrum. The idea behind a velocity spectrum is to display some
type of signal coherency as a function of the NMO velocity selected to “flatten” a
reflector (i.e., correct for normal moveout) plotted as a hyperbola in a CMP data set.
Signal coherency is calculated by summing the amplitudes within the data set along each

two-way travel time. When the correct NMO velocity is selected and the hyperbola is

39

“flattened”, the signal coherency is greatest because the positive and negative phases of
the sinusoidal EM wave are aligned.

The signal coherency of the CMP data collected at FT A-02 was determined for a
wide range of NMO velocity values. A velocity spectrum was created by plotting the
coherency values at each two-way travel time (i.e., the sum of amplitudes at each offset
distance) as a function of NMO velocity (Figure 17). Peak values corresponding to a
coherent signal are shown as bright spots in the velocity spectrum. The type of media
and presence of unsaturated and saturated sediments at the PTA-02 site constrained the
range of velocity values used in the velocity spectrum analysis (i.e., x-axis on Figure 17).
As mentioned above, if only the portion of the CMP data set comprised of small
transmitter/receiver offsets relative to reflector depths are used, then the NMO velocity is
equivalent to the RMS and stacking velocities (Yilmaz, 1987).

A Matlab® code was created (see Appendix) to perform velocity analyses on CMP
gathers and generate a velocity spectrum. It is a universal code that can be applied to any
data set as long as several basic assumptions are met (e.g. proper coupling of the
antennae with the ground surface, subsurface reflectors are horizontal, etc.). The code
works, as described above, by using a range of user~supplied NMO velocities to
deterrrrine the “best-fit” velocity values to “flatten” the hyperbolic representation of an
individual reflector. The necessary inputs into the code are the CMP data set (referred to
as agx by the code), a vector of values corresponding to the incremental
transmitter/receiver offset values (offset), time of the first data arrival (t1), sampling time
interval (dt), minimum and maximum NMO velocity to be analyzed (vnrin and vmax,

respectively), increment to increase vmin by during each iteration of the calculation (dv),

40

.EEooam .0629. on. we 2...? 83an o... «o owns: @5800 E 3.58% b.8.u> 2 .2 25mm.

.2. 83.0%.. .2. 823..
sis:ls..... 392:2...3

 

 

      
     

.. III
I (W'IIII'FII lIrllIl I‘IIIIIWIIII III

glib. I... II. .I.n......i

 

 

 

 

 

     

 

 

 

 

 

 

a. . w
s . I'NII
if I III
a. . I I. Ill . I I m
.11. W. J
2. l I
n. I, J
m IIIIIIIIII - M
S . IA 5
a . .l l r u
( 9. I m(
IIV.II.IIIHrIIH. :hIIIAIIUII . .. r - ON I .
Err tut-“Mg.” “Mural II U I E‘I‘IN
Irll blur..- .. a. l .I-. - .1 III“... a
. r.I.~ IF
2
. . on s. I I
spasm gaugesoasme eats a .2. Sea. 2.8.... 2 .2.

41

and the mutestretch value (mutestrch). The mutestretch value truncates data
corresponding to large transmitter/receiver offset values in relation to the depth to the
reflector. The mutestretch value is also used to delete (mute) data within the shallow
events that were stretched and subjected to frequency distortions during normal moveout
velocity correction (Yilmaz, 1987). Note that the code works best using a CMP data set
that has been visually enhanced by an AGC. By using a data set with an AGC applied to
it, reflector strengths are enhanced and signal coherency is more prominent. The velocity
analysis does not utilize the actual amplitude values of the reflections as part of its
calculations. The velocity analysis only considers the coherency of signals. Therefore,
using a CMP data set that has been visually enhanced by an AGC is appropriate. The
output parameters generated by the code are NMO velocities (equivalent to RMS and
stacking velocities) for the media above each reflector, interval velocities of the media
between each reflector (calculated using the Dix equation), the two-way travel times
corresponding to the initiation of each interval velocity, and the thickness of each layer.
The interval velocities and their corresponding times were used to convert two-way travel

times recorded in the GPR data into depths.

Results of Velocity Spectrum Analysis
For the analysis of the GPR data from the PTA-02 site, the input values into the

code are listed in Table 3.

42

 

Table 3. Input values for velocity spectrum code.

 

 

 

 

 

 

 

 

 

 

 

Input Abbreviation Value(s)

CMP data with AGC applied agx Data from FI‘A-OZ site
A vector of values corresponding to the offset 2 to 39 m, at 0.5 m
incremental transmitter/receiver offset values increments

Time of the first data arrival t1 0 seconds
Sampling time interval (it 1.2 x 10'9 seconds
Minimum NMO velocitlto be analyzed vmin 8 x 107 m/s
Maximum N MO velocity to be analyzed vmax 1.394 x 108 m/s
increment to increase vmin by during each dv 0.6 x 105 m/s
iteration of the calculation

Mutestretch value mutestrch 35

 

 

The coherency image generated by the code is shown in Figure 17. The NMO velocities
for the media above each of the two reflectors shown in Figures 15 and 16 were then
determined from the velocity spectrum. Figures 18 and 19 show close-up images of each
reflector’s coherency within the velocity spectrum. In each case, points are shown on the
figures that identify maximum coherency values for the positive and negative phases of
the reflected electromagnetic wave. The NMO velocities corresponding to the maximum
coherency values for the positive and negative phases of the individual reflector were
averaged to best approximate the velocity for the zero-phase portion of the reflected
electromagnetic wave. In addition, the two-way travel times corresponding to the
maximum coherency values for the positive and negative phases of the individual
reflector were averaged to best approximate the two-way travel time for the zero-phase
portion of the reflected electromagnetic wave.

Based on this analysis, the NMO velocities for the material above each of the

identified reflectors were calculated (Table 4).

43

 

.8320on .062? 05 £53, 860%. 038
.5225 mo meme—E gunmen can 928m .8 3:80:00 825.88 he 20232 9.3352 850m .w_ BamE
.2 x 3.5. 23.5
v.2 on; one vné mm; m; an; ma... em;

x

 

(395) mm

 

 

£2.88 a_8_o> we. x

.Esboomm .0829, 05 55.3 580:2 com”:
.86 Mo 3923 gamma: was 9.282 .8 20:20:00 E2832: “o «cowa8_ @5252 $50..— .2 QEME

.2 x 3.2. £825
S 3 m m... o n. a 3 o

 

o..-

(095) emu

on. .

 

 

 

£88on b_oo_o> ..o_. x

45

 

Table 4. Calculated NMO velocities

 

 

 

 

Reflector/ Material Above Reflector NMO Velocity of Material
Capillary fringe / Unsaturated Zone 1.35 x 1(Trm/s

Clay layer/ Saturated and Unsaturated 8.07 x 107 m/s

Zone

 

 

The following relationship was used to compare the unsaturated velocity value with that

obtained in laboratory settings:

3*108

J?

V:

 

(2) (Annan, 2001)

where V = velocity (m/s), 3 *108 = speed of light in a vacuum (m/s), and K = dielectric
constant. Dry soil has a dielectric constant of 4. Therefore, the velocity of an EM wave
in dry soil is 1.5 x 108 m/s, which is similar to the value of 1.35 x 108 m/s estimated at the
PTA-02 site. Differences in velocity calculated in the laboratory and in the field are due
to slight differences in the electrical properties of different soil types and residual liquids
and/or vapors within a small fraction of the pore volumes. The velocity of the EM wave
is expected to decrease rapidly when it encounters saturated soils because the dielectric
constant of water is 80, which is much greater than 4, the dielectric constant of soil.

(Annan, 2001).

Xz-T2 Method
A second technique used to estimate the velocity of the subsurface media from a

CMP gather is the XZ-T2 method, where x is equal to the transmitter/receiver offset

46

 

distance and t is equal to the two-way travel time of the reflected electromagnetic wave.
The two-way travel time is a function of the velocity of the media above the reflector, the
depth to the reflector, and the offset distance between the transmitter and receiver. This

method is based on the relationship:

I2 = 7):,— (3) (Sheriff and Geldart, 1995)

V5 +10
where Vs equals the stacking velocity and to equals the time of the initial recorded signal
(Sheriff and Geldart, 1995). If we plot t2 as a function of x2, and the velocity is constant,
the result is a straight line with a slope of V,2 and a y-intercept of toz. Sheriff and Geldart
(1995) and Yilmaz (1987) noted that when the reflectors are horizontal and offset

distances are small when compared to the depths to the reflectors, the stacking velocity

calculated using the above relationship is nearly equivalent to the RMS velocity.

Results of X2 -T2 Analysis

The two-way travel times and offset distances measured from the reflectors
shown in Figure 16 were analyzed using the XZ-T2 method to determine the velocity of
the horizontal subsurface layers (Figure 20). Note that the velocities within each layer
are inferred to be constant and the layers are essentially horizontal based on the high R2
values ranging between 0.9929 and 0.9999. The calculated RMS velocities for each
identified reflector are listed in Table 5. The velocities within each layer, calculated

using the velocity spectrum and the Xz-T2 method, are very similar. Each velocity

47

deuce—me :03 8a max—«5 «Wax .on BamE

 

83.0" Na
2.mn+x2.w~ u»

$3.: u “a
3.2 +5.36 a»

 

92!...» .28. 5 «3.5
cm: comp omN— ooow 0mm com own 0

 

 

 

 

o
v Tum
n Twp
Sauce... .595
.osofim 30'.
n 7mm. _.
2:23
833.3. .525 n TNN
23 afisl
n Tmmfi
Bantam >20 III M Tun
n Tumd
2%... 55.5 lol
M Tu...

 

(ponnbs ‘spuoaofl emu

 

 

nos—Er «fl 2! 323:3—
oE—Qfloah .32.:

48

 

Table 5. Calculated RMS velocities

 

 

 

 

Reflector / Material Above Reflector RMS Velocity of Material
Capillary fringe/ Unsaturated Zone 1.33 x 108 m/s

Clay layer/ Saturated and Unsaturated 7.84 x 10%3

Zone

 

 

analysis method has similar estimated velocity values for the subsurface media.

All assumptions about the subsurface materials, data acquisition methods, and
data analysis methods described below are satisfied at the PTA-02 site. Proper antennae
coupling with the ground surface was maintained at all times. The selected mutestretch
value eliminated data collected from large offsets relative to reflector depths. Also,
compiled well-log data, core analyses, and common shot gather analyses all indicate

nearly horizontal layering within the shallow subsurface at the PTA-02 site.

Interval Velocities

The interval velocity (Vi) is defined as the average velocity over some interval of
the electromagnetic wave travel path. At the PTA-02 site, the interval velocity of the
media between the capillary fringe reflector and the clay layer reflector was calculated

using the Dix Equation:

2 2
___Vn :l‘tri‘-"/rl—l*tl
t —t_l

II II

V.2

”‘1 (4) (Sheriff and Geldart, I995)

 

where VII and tn are the RMS velocity and zero—offset two-way travel time corresponding

to the nth reflector, respectively (Sheriff and Geldart, I995). The Dix Equation is very

49

 

unstable, which has been documented by Sheriff and Geldart, 1995 and Yilmaz, 1987.
Slight deviations in RMS velocities often result in exaggerated differences in calculated
interval velocities. As a comparison to the interval velocity values obtained from the Dix
Equation, a weighted velocity average based on the travel-time within each layer was
used to calculate the interval velocity. Table 6 summarizes the interval velocities

calculated using each method.

 

Table 6. Calculated interval velocities.

 

 

 

Analytical Method
Velocity Spectrum XTuTY Method
(m/sec) (m/sec.)
Medra above caprllary fringe (i.e., 1.353 x 103 1.329 x 103

unsaturated media)

 

Interval velocity of saturated media
below capillary fringe and above clay 6.70 x 107 6.40 x 107
layer calculated using Dix Equation

 

 

Interval velocity of saturated media
below capillary fiinge and above clay
layer calculated using weighted
average

7.13x107 6.87x107

 

 

 

An estimated interval velocity value of 6.78 x 107 meters/second for the saturated zone
will be used for firrther data analysis, which is equivalent to the average of the four
interval velocity values presented in Table 6. Equation (2) was again used to compare the
estimated saturated velocity value with those published in literature. Wet soil typically
has a dielectric constant of approximately 25 (Annan, 2001). Therefore, the velocity of
an EM wave passing through wet soil is approximately 6 x 107 m/s, which is similar to
the value of 6.78 x 107 m/s estimated at the ETA-02 site. Differences between laboratory
and site-specific estimates are due to slight differences in the electrical properties of

different soil types and the presence of different ions in the groundwater (Annan, 2001 ).

50

 

CONVERTING TIME TO DEPTH IN A GPR TRANSECT

To estimate depths from measured reflection times, an accurate velocity analysis
is required. As noted above, accurate velocity calculations were made to properly
convert time to depth within the GPR data. The depths equate to EM wave propagation
distances, which are necessary in calculating the attenuation of the wave as it passes
through the subsurface media. The RMS velocities within the media above the capillary
fringe and above the clay layer were calculated using the velocity spectrum analysis and
the Xz-T2 method. The interval velocity, also referred to as the average velocity of the
media between two reflectors, was calculated for the media between the capillary fringe
and the clay layer. Using these calculated velocity values, the depths to the capillary

fringe and clay layer identified at the PTA-02 site were calculated (Table 7).

 

 

 

 

 

 

 

 

Table 7. Calculated depths to capillary fringe and clay layer

Analytical Method2 2
Depth from Ground Surface Velocitif‘pectrum Meir-rill (m)
Depth to Capillary Fringe 5.56 5.58
Depth to Clay Layer (based on two-way travel
time to clay reflector and estimated NMO 22.66 21.91
velocity)
Depth to Clay Layer (based on depth to
capillary fiinge + depth to clay using interval 21.81 21.69
velocity)

 

 

 

Comparison of Calculated Depths to Measured Depths

Depth-to-water was recorded from the multi-level well network at the PTA-02
site on May 12, 1999, two days prior to completing the GPR survey (Figure 2). No
precipitation event occurred between May 12 and 14, 1999. Therefore, the change in the

water table elevation during that time is considered negligible. The depth-to-water

5|

(below ground surface) measured in multi-level well ML—9, which is the closest well to
the center of the CMP gather, was 5.67 meters. Multi-level well ML—9 is located 36.5
meters southwest of the CMP gather’s central location. The location, and ground surface
and easing elevations were surveyed for each well, and referenced to NGVD. Therefore,
the groundwater elevation at ML—9 is calculated by subtracting the depth-to—water
(measured from top of easing) from the surveyed casing elevation. Surveyed ground
surface elevations and calculated water table elevations were interpolated between known
points along the GPR transect and are shown in Figure 21. Note that the water table
elevations were extrapolated to the southwest of MW-lO and to the northeast of ML-9.
Since there is no well located at the center of the CMP gather, the extrapolated water
table elevation at this location is the best approximation. The error associated with this
approximation is minimal because the maximum difference in groundwater elevations
measured between any two multi-level wells during this sampling event and located
along the GPR transect is 0.27 meters. Subtracting the extrapolated groundwater
elevation from the surveyed ground surface elevation at the center of the CMP gather
results in a depth-to-water estimate of 5.66 meters. Table 8 lists the differences between

the calculated depth to the capillary fringe and the estimated water table depth.

 

Table 8. Differences between measured and calculated depths to the water table
and capillary fringe.

 

 

 

 

 

 

 

Estimated Depth to Calculated Depth to Calculated Depth to
Water Table (m) Capillary Fringe Using the Capillirry2 Fringe Using
Velocity Spectrum (m) the X -T Method (m)
5.66 5.56 5.58
Difference Between
Measured and Calculated 0.10 0.08
Depth (m)

 

 

 

52

 

.8335? 038 533 c5 95on 836985 .3” Esmfi

23.2.38an
8? cm: o: ofi 8F 8 8 8 ow o

 

 

A _

 

mm _.
03.3 .38?

IN
Q
(srerotu) uonamla

53

1mg

\\.\I1/./ll\l r
\

llllll/ \\

, .

._..\. 353% 335.5 . omw

 

_ _ _ _ _ F _ _ _ FQF

 

mco_~a>o_m ”.38. .235 was 9.590 nos—0902.

As shown above, each velocity analysis method has calculated the depth-to-water very
closely to the estimated value. As previously noted, the electromagnetic wave actually
reflects off the top of the capillary fringe resulting in a slight discrepancy (i.e., shallower
estimate) between depth to the actual water table and the same depth recorded by GPR
(Barber and Morey, 1994). Therefore, the differences between the estimated depth to
water table (from extrapolated direct water table measurements) and the calculated depth
to the capillary fiinge (from GPR data) are expected. Annan and Chua (1992) noted that
the water table is a highly gradational boundary encountered by an electromagnetic wave.
Capillary action can smear the water distribution between the unsaturated and saturated
zones as a function of grain sizes (Annan and Chua, 1992). The typical capillary rise
observed in medium to coarse sands, which are present at the site, is 0.25 to 0.15 meters,
respectively (Fetter, 1994). These capillary rise estimates are very close to the
differences between the estimated depth to water table and the calculated depth to the
capillary fringe listed above.

The clay layer identified in the GPR transect at the FTA-02 site ranges from 24.4
meters below grade in the northern part of the former Wurtsmith Air Force base to 9.1
meters below grade at the west edge of the former base, with an average depth of 19.8
meters below grade. In general, the deepest occurrences of the clay layer are in the
northeast portion of the former base (Stark et. al., 1983; Nash et. al., 1997). The exact
depth to the clay layer has not been measured at the location of the GPR transect and
CMP gather. However, calculating the difference between the ground elevation (Figure
22) and the clay surface elevation (Figure 23) in the southwest portion of the former

WAFB at the PTA-02 site estimates the depth to clay as approximately 19.81 meters (65

54

  

Amwa— .eoEoBH can anemia—=0 .fiva 83 8.5"— :< 5:353 583m 3 St 80930:. 9:55 .NN Bamfi

nip-.88. n u . o a. 3 a)...

.L..|_.J|_l. it...» £21.g.191
3.5; Eggsgaioozlaroool

§k<2¢4§u
g.

.8...

o

bag:

55

.Gwfl .5225. can magi—:59 omen 88m .:< sugar—=3 5.58 “a Cb gouge—o course. .35 .mm BawE

m0<0 mOCOu ¢.<

 

 

 

 

 

2:30.524: 8 >¢<0238 lrl

#540 if...

I Q ”Eda—u
80$ Zghowa no
zO..-.<0°.. wh<X§O¢tt< 3"

ZQ+<2<Jnxm

 

 

OCWFW3°i¥ N p O
O “'
ON...’ a p O L a D - 0313’ 800“."0
W .38 3!
.ONO" colt . .1
\uk .0 r/ L I8 :8;
IL / LI- - .
1.0 . I v I: I3 :2:
9
O O .
:95: m. _
0:3
9 .
l I
.<
7
W A
o
9
393 m. a g
0? cow-duo“:

0
MO

 

d9}

04’

\
\o"

 

 

3an

50-000

56

feet). Table 9 lists the differences between calculated and estimated depths to the clay

layer.

 

Table 9. Differences between measured and calculated depths to the clay layer.

 

 

 

 

 

. *Calculated Depth to *Calculated Depth to Clay
5:211:13:ng $3533) Clay Layer Using the Layer Using the Xz-T2
Velocity Spectrum (m) Method (mL
19.81 22.66 21.91
Difference Between
Estimated and Calculated *2.85 *2.10
Depths to Clay (m)

 

 

 

*Note that the depth to the clay layer was calculated using velocity values obtained from
the velocity spectrum and the )(2-T2 method and the two-way travel time from the surface
to the reflector. Depth calculations using the calculated interval velocities were very
similar. The degree of accuracy to calculating the depth to the clay layer from the

acquired GPR data is unknown without an exact depth measurement.

SPHERICAL DIVERGENCE CORRECTION

The energy passed through the electronics of the transmitter and induced into the
subsurface as a propagating wave front is not equal to the energy recorded in the
receiving antennae. Electrical losses and spherical spreading of the wave front attenuate
the energy propagating through the ground. In addition, some of the energy that reaches
the target reflector is scattered by the reflector. The reflected energy is once again
subject to attenuation by electrical losses and spherical spreading as it propagates through
the media prior to reaching the receiving antennae (Annan and Chua, 1992). As such, in

order to quantify the amount of attenuation an EM wave experiences due to electrical

57

 

losses as it propagates through the subsurface media, which is a function of the
conductivity of the media, energy loss from spherical spreading must also be quantified.
The wave field generated by the transmitting antennae is conceptually thought of
as a point source that generates a spherical wave field. As the spherical wave field
spreads outward from the point source in a homogeneous media, the energy density
decays proportionately to (1/r2), where r is the radius of the wave front, which
corresponds directly to the propagation depth of the wave (Yilmaz, 1987). As such, a
time-variant gain was applied to the filtered GPR transect data to correct for spherical
divergence. Amplitude is proportional to the square root of the energy. Therefore, the
recorded energy at each point within the transect was multiplied by its corresponding
depth (depth calculation discussed above). The GPR transect corrected for spherical
divergence is shown in Figure 24 (see Appendix). The spherical divergence correction
resulted in a GPR data set with slightly enhanced reflection amplitudes at depth. The

gain applied to the GPR data is relatively small due to the shallow features at the site.

AMPLITUDE AND ATTENUATION ANALYSIS
The novel approach to quantitatively assessing the attenuation of the transmitted
GPR signal and subsequently calculating the specific conductivity of the groundwater is
to analyze amplitude changes with distance. As mentioned above, previous studies by
others have not analyzed GPR data to quantitatively assess hydrogeologic and
geochemical properties of groundwater. Instead, GPR has been used to simply image
and/or delineate plume geometries. Researchers have established typical attenuation

values for a range of common materials by using empirical relationships between the

58

.ouaomhozu Eoconam H8 “cocoa 9 teams Ea» 53» 8mm «CANE an 8853 ~50 .vm 0.59.;

A2395 8:820
SN 09. 09. 03. our 00—. on om ow on o

V
(spuoaas) cum

6")

 

sneak Eo so? x

59

electrical properties of the materials and EM waves (Sensors & Software, 1996).
However, these typical attenuation values are not useful when a complex contaminated
system is encountered and many of the factors contributing to attenuation may be
unknown. As such, an approach has been developed to quantify the amount of
attenuation directly from the GPR data.

The decrease in amplitude of a propagating EM wave front in rocks and
sediments due to attenuation (i.e., absorption of the wave front by electrically conductive

materials) is exponential with distance and can be represented as:

A = Aoe'a“ (5) (Sheriff and Geldart, 1995)

where A is the amplitude of the initial wave front, A0, afier traveling a fixed distance, x
and experiencing attenuation coefficient, a (Sheriff and Geldart, 1995). By measuring
the reduction in amplitude of the transmitted GPR signal, the attenuation coefficient (a)
of the shallow aquifer material at the PTA-02 site was quantified.

The decrease in amplitude between the water table and clay layer reflections
observed at the PTA-02 site was analyzed to quantify the amount of attenuation. Figure
25 shows the amplitudes of the water table and clay layer reflections. This method of
calculating attenuation as a function of distance has the advantage of utilizing well-
defined reflection events that are easily identified and analyzed in the GPR data set.
However, this method is limited to conditions where reflectors are identified and does not
include amplitude reduction independent of reflection amplitudes, which may have

resulted in the boundary effects discussed below.

60

.mcoLLoocoL 8.3L .36 EB 038 583 Le mucBLLLLELLL .3 am

 

 

 

 

 

 

 

 

 

 

 

 

 

 

LmLoLoELoocaLmLO
com 02 09 o: om? 8., 8 om 3 cm o S
L a L L L 1 L d L o W
no”...
_ _ . I; L»: m
L. L: L. LLLL _ L . - .L w
_L a. ‘L _ W
- -Bm.
m
- -N w
m.
L L L L L L L r L . W...
man

93 L66 .6 83.95 LS x

68.058535
con 02 09 03 our 9: on om 9. on o S
L A L L L 4 L L L O n
w
w
l [N u
m
r .4 W
w
r b ‘- ..m m
D:< VLF; in ...<4..1 m
Ly 5;) : >> ' 1 ‘ ILQ W
LLL. LLLLLLLL LLLLLL LLLLL L... .L .L.
F r L L L L m.
ore

 

 

 

can» Lass Lo 82.3.5 L2 x

61

The amplitude of each reflection was calculated by summing the absolute
amplitude values along one complete wavelength of the reflected signal (see Appendix).
The amplitude of each reflection was also calculated using the maximum absolute
amplitude value along one complete wavelength of the reflected signal. The relative
reflection amplitudes differ depending on the calculation method, but the energy loss
(i.e., attenuation) between reflections is nearly identical. Therefore, it is believed that
either method can be used. Also note that the amplitude of the water table reflection is
significantly reduced between approximately 160-meters and 244-meters within the GPR
transect, which corresponds to the location of the identified groundwater contaminant
plume. This observation recorded in the GPR data may be caused by the presence of
light non-aqueous phase liquid (LNAPL) contaminants on the water table surface and/or
a significant concentration of residual and/or vapor phase contaminants within the vadose
zone. These conditions do not meet the assumptions of the model presented herein.
Therefore, the amplitude of the water table reflection between approximately l60-meters
and 244-meters within the GPR transect will not be used to quantify the attenuation.

Trend lines were fit to the water table and clay layer amplitude curves (Figure
25). The trend lines smooth slight amplitude differences between each GPR trace that
may have been caused by differences in antennae coupling and incident angles. Trend
lines were fit to the western and eastern portions of the water table amplitude curve,
excluding the area between l60-meters and 244-meters within the GPR transect (Figure
25). The trend lines for the clay layer amplitude curve and the western portion of the

water table amplitude curve were best fit using a fourth-order polynomial. The trend line

62

for the eastern portion of the water table amplitude curve was best-fit using a first-order
polynomial (see Appendix).
By solving the exponential decay relationship presented in Equation 5 for a, the

attenuation of the EM wave front is defined as:

moi)

 

a:

(6)

—X

Therefore, the attenuation of the EM wave front within each GPR trace between the

water table ( A0 , trend line) and clay layer (A , trend line) reflections are quantified using

Equation 6. The distance (x) is calculated from the GPR data as one half the difference in
two-way travel times between the water table and the clay layer, multiplied by the
velocity of the saturated media (see Appendix; Figure 26).

As shown in Figure 26, the increased attenuation values near the western and
more notably eastern ends of the GPR transect appears to be an artifact of the best-fit
trend line for the clay layer reflection amplitude curve and apparent reduced clay layer
reflection amplitudes. The trend line is a function of a fourth-order polynomial and
decreases in value near the western and eastern edges of the GPR transect. As such, the
calculated attenuation values in these areas have been artificially elevated. The
reductions in clay layer reflection amplitudes in the western and eastern portion of the
transect were unexpected and the cause is unknown. A different method of analyzing

amplitude loss may be required to improve the attenuation quantification in these areas.

63

com

09.
L

cow
L

.LLLoLL 2&3 Em Lo =0ua==§< .cm Baum

03.
#

LEoLoEL 00535
ON _. 00 L. ow
L q

om

ov

on

 

 

L

L L L

 

cud

and

mud

n .

no
N_ N
o o
(”J/8P) UODWUOW

a)
“l
o

‘9.
o

5.0

mmd

 

 

D>a>> 2W *0 COBGDCDLL<

mmd

SPECIFIC CONDUCTIVITY CALCULATION

Elevated specific conductivity was measured in groundwater samples collected
within the contaminant plume (discussed later and detailed in Table 10) and, as
previously discussed, a ‘shadow zone’ was observed in a GPR transect through the
plume. The ‘shadow zone’ is caused by EM wave attenuation, which can be calculated
from the GPR data within each trace of the GPR transect. As described below, empirical
relationships between attenuation and conductivity parameters enable the calculation of
the specific conductivity of the groundwater at the PTA-02 site.

Gueguen and Palciauskus, 1994 indicated that the properties of an
electromagnetic wave can be described by using solutions to Maxwell’s equations in
three dimensions. However, these solutions are mathematically very complex.
Therefore, to investigate the important features of the electromagnetic wave, they
simplified the mathematics by applying the solution for a harmonic plane wave to
Maxwell’s Equations for electromagnetic waves. This is a well recognized and standard
approximation for analyzing electromagnetic waves. The results define the equation for

the electromagnetic wave propagation and attenuation in a medium, such that:

 

k = \[(e,ua)2 +1100) = km, + ik (7) (Gueguen and Palciauskus, 1994)

imaginary

where k is the complex wave number, e is the electrical permitivity (F arads/meter), p is
the magnetic permeability (Newton/Amperez), ? is the angular frequency (MHz), and s

is the conductivity of the media (Siemens/meter). The real part of It describes the

65

propagation of the electric field and the imaginary part of k describes the amount of

attenuation with distance. The imaginary part of k, or attenuation (a) is defined as:

a = (0* Jig—l— * J‘/1+(—0;)2 -l (8) (Gueguen and Palciauskus, 1994)
£60

The attenuation of the EM wave is therefore a function of conductivity and other

 

electrical properties of the media. As described above, the attenuation of the wave can be
calculated from the GPR data by evaluating the exponential decay of reflection
amplitudes. The direct relationship between attenuation and conductivity can be

evaluated as the derivative of Equation 8 with respect to 0:

LEE
6a_ 0' 2

53—
flrul)? —1*,/1+(£)2
£60 £0)

If we assume that changes in conductivity are greater than any other changes in the x-

 

 

(9)

dimension and are the only cause of differing amounts of attenuation, then the following

relationships are also true:

80 A0

66

A0:
= Orr—1 + ‘53-— (l l)
/
Go

0'

II

where n is the individual trace number within the GPR transect. As such, the change in
conductivity (A o) of the media within each GPR trace can be calculated using the change
in attenuation (A a) between traces and the derivative of the equation that defines the
attenuation (a a / a o). A a between individual traces within the GPR transect can be
calculated directly from the attenuation values shown in Figure 26. A o is calculated
using the relationships above. However, an initial conductivity value is required to
calculate A o. The initial conductivity value used for the western end of the transect is
180 micro-Siemens/centimeter, which is equal to the background specific conductivity
values measured by others at the site (Bermejo et. al., 1997; Nash et. al., 1997; Sauck et.
al., 1998; Sauck, 2000) and approximately equal to the average specific conductivity
value of 188 micro-Siemens/meter measured at four depths within multi-level well ML-
10.

The conductivity calculated from the relationships presented above should be
referred to as the effective conductivity (om) because the propagating EM wave front is
subject to absorptive losses from the conductivity of both the sediment grains and the
pore waters. Because there is such a high contrast between the conductivity of the
silicate sediments ((38 z 10’14 to 10'‘0 Siemens/meter) and ion-rich fresh water solutions
occupying the pore volumes (oW 2' 10'1 to l Siemens/meter), it is understood that most
electrical currents flow through the pore waters. Therefore, the effective conductivity of

saturated sediments is equal to the specific conductivity of the fluid (0,.) reduced by a

67

formation factor (F), which is dependent on pore-size parameters, including porosity (F )

and tortuosity:

a; (12), (Gueguen and Palciauskus, I994).

 

0",” z

The microstructure details required to directly calculate the formation factor are not
known for the sediments at the FTA-02 site. As such, Archie’s Law was used to
determine the relationship between the formation factor and porosity. Archie’s Law

defines the formation factor as:

F z ¢“'" (13), (Gueguen and Palciauskus, 1994)

where m is referred to as a cementation factor, which increases as tortuosity of the pore
volumes increases (Backe et. al., 2001). As mentioned above, the approximate porosity
(F) of the sediments is 0.3, which falls within the range of typical porosity values for
unconsolidated sand deposits (Freeze and Cherry, 1979) and can be expected of glacial
deposits with various grain sizes including lacustrine sediments, ice-contact sediments,
and alluvium present at the PTA—02 site. A typical range of values for the cementation
factor for unconsolidated quartz sand is 1.3 to 1.5 (Backe et. al., 2001; Bigalke, 2000;
Boving and Grathwohl, 2001; Fittennan and Deszcz-Pan, 2001; Hwang et. al., 2004;
Kwader, 1986) and a typical range of values for the cementation factor for geologic
materials is 1.5 to 2.5 (Gueguen and Palciauskus, 1994). However, no universal

cementation factor exists for individual sand types because m can vary simply by the

68

orientation and shape of the pores (Boving and Grathwohl, 2001). In addition, modeling
by others has calculated a cementation factor of 3.2 in unconsolidated sands (Backe et.
al., 2001 ). The model results for the calculated groundwater specific conductivities were
calibrated to a weighted average of the measured groundwater specific conductivities,
described below, using a formation factor of approximately 16.9. Therefore, according to
Archie’s Law and if we assume a porosity of 30%, the cementation factor is equal to 2.35
and the effective conductivity of the media and background waters at the FTA-02 site is
equal to approximately 10.6 micro-Siemens/centimeter. A cementation factor of 2.35 is
slightly higher than values typically used for unconsolidated quartz sands. However,
2.35 is a realistic value for typical geologic materials. In addition, the varying redox
environments within the contaminant plume may facilitate the production of various
coatings on the individual sand grains, such as metal-sulfides and/or biomass from the
active microbial communities. These materials would result in decreased pore volumes,
increased tortuosity, and increased cementation factors.

It is important to note that the site specific porosity and cementation factor have
not been independently quantified for the PTA-02 site. Therefore, according to Archie’s
Law, there are an infinite number of porosity and cementation factor combinations that
can result in a formation factor of 16.9. The relationship is highly sensitive to changes in
the cementation factor. If porosity is as low as 15% then the cementation factor required
to maintain the same formation factor is approximately 1.5. The model is calibrated
using a formation factor of approximately 16.9 and additional research is required to

quantify the site specific porosity and cementation factor.

69

The calculated specific conductivity of the groundwater at the PTA-02 site is
shown in Figure 27 (see Appendix). The highest electrically conductive groundwater is
located near the center of the GPR transect, which directly corresponds to the location of
the contaminant plume. Once again, the increased specific conductivity values calculated
near the western and, more notably, eastern ends of the GPR transect appear to be an
artifact of low clay layer reflection amplitudes in these areas. Therefore, the calculated
attenuation in these areas is over-estimated because there is no suspected contamination
there.

The attenuation along each GPR trace is calculated from the water table reflection
to the clay layer reflection, and represents the total attenuation encountered within the
saturated media above the clay layer. As such, the estimated groundwater specific
conductivity values shown in Figure 27 are representative of the average specific
conductivity values along each trace of the GPR transect within the saturated media
above the clay layer. Therefore, the calculated specific conductivity values were
calibrated (by adjusting the formation factor, as discussed above) to a weighted average
of the measured specific conductivity values (Figure 28, Table 10).

The model for calculating the specific conductivity of groundwater using GPR
has been calibrated to the direct measurements recorded within multi-level wells ML-3,
ML—4, ML—S, ML—7, MIL-8, ML—9, and ML-lO at the PTA-02 site between June 5 and
June 10, 2000, which are shown as asterisks on Figure 29. The values measured at each
sampling depth within each multi-level well have been averaged using a weighted
average to represent the composite groundwater the EM wave front passed through in the

GPR transect. The weighted averages were calculated based on depths defined by the

70

oom

.LoLwBLELLBw Lo 335288 #:8an 322830 R ELLMLnL

$3225 8:830

90..
L

on

om

ov

ON

 

 

02.
L

09
L

o:
L

L

our
L

L

L

L

L

h

l
O
C)
('7

l
o
o
v
(ma/Sn) A“Allill'lpuoi'D

102

1 com

com

.. oom

 

oom

 

3.023396 Lo 53:03:00 050on 35330

71

02
Ba
0 .n
omm
8m
one
cow
own
8m
0mm
8N
om L
Esme

biggoo
058mm

 

don—g b35328 058% 58395on 6220985 .wm 85mm

942

wNLLz

mHLSL

mudz

L142

 

03L >32

5-42

OLE—2

 

Doom 0:3.

72

.8333:on Lo b38258 2.38% 33308 was BEL—GEO 6N 2&5

 

 

 

 

L288: 853.0
com our 09 03. our 8.. on oo 9. on 00
- 37/ 1
\ x/ 09
“Z /
#1 K. HgQK/J

I \fx 1 com
«a. m
u + 1 Don P
\l\f<f% flag m...
x / m
I 517 + 1 n
\ + ,K/ + 00¢ W.
\ + u (

D
m 1 8m
32 @2 n $5 Mao +
55 89 $2: pang
I 1 com
55 coon 2:; .. #822
LuzowouL a
u H _ H L _ _ _ _ r 1 con

 

 

 

 

$332390 Lo 33:26:00 05% 8.5322 Luau @28630

73

water table elevation and the midpoints between sample depths. The GPR data used to
create the model presented herein was collected approximately one year earlier in May
1999. However, groundwater specific conductivity was only measured in wells ML-3,
ML-4, ML-S, and ML-8 in 1999 (Figure 28). The June 5 to June 10, 2000 data set of
measured specific conductivity values is the most comprehensive data set for the site.
The measured values from each sampling event are similar, with the exceptions of ML-4

and ML-S, as shown in Table 10.

 

 

 

 

 

 

 

 

 

 

Table 10. Weighted Average of Measured specific
conductivi values.
Specific Specific
Well Conductivity Conductivity
(uS/cm, May 1999) (uS/cm, June 2000)

ML-3 447.6 400.4
ML-4 316.7 530.0
ML-S 672.5 453.2
ML-7 Not Measured 298.4
ML-8 469.5 388.1
ML-9 Not Measured 427.3
ML- 1 0 Not Measured 190.7

 

 

 

 

The measured specific conductivity of the groundwater was more variable amongst the
multi-level wells in 1999 than in 2000. This may be the result of slightly drier conditions
in 1999. Water table elevations measured in the multi-levels wells on a nearly bi-
monthly basis were approximately 0.8 meters lower in 1999 than in 2000, and the
sampling event followed a 3-month period of low precipitation and declining water table
elevations. The drier conditions may have resulted in areas of differential biodegradation
and resultant specific conductivity measurements. In contrast, the 2000 sampling event

occurred afier a 4-month increase in water table elevations, indicative of greater

74

precipitation volumes. The increased precipitation may have leached ions into the
groundwater and created a slightly more homogeneous contaminant plume with regards

to specific conductivity measurements.

CONCLUSIONS

The research presented herein provides a new mechanism to evaluate the nature
and extent of groundwater contamination using high-resolution GPR methods from the
ground surface. The sequence of data processing techniques produced a calibrated site-
specific model that quantitatively evaluates the specific conductivity of groundwater
within a contaminant plume. The model is based on data from the study site and
documented investigations at other sites that show contaminant plumes with well
established microbial activities and undergoing biodegradation will exhibit groundwater
with elevated specific conductivities as a result of ions leaching from the sediments. This
differs from contaminant plumes that have not undergone significant biodegradation
and/or have LN APL at the water table, which is highly resistive and enhances GPR
reflections. The increased groundwater specific conductivity identified on-site is also the
result of increased sodium mass discharged from a bioreactor as part of a pilot-scale
pump-and-treat system installed at the study site. The model increases the usefulness of
GPR from simply a tool used to image and delineate a contaminant plume to a non-
intrusive tool for estimating the specific conductivity of groundwater.

By altering the site-specific parameters, the model can be used on a number of
sites, given general knowledge of the aquifer parameters, sufficient GPR signal

penetration, and limited naturally occurring conductive materials such as clays and

75

naturally conductive groundwater. The model enables the user to quickly and accurately
evaluate a specific aspect of the groundwater quality (i.e., specific conductivity) without
installing monitoring wells and collecting groundwater samples for laboratory analysis.
The presence of elevated groundwater specific conductivity values above background, as
calculated by the model presented above, is consistent with the location of groundwater

contaminant plumes.

76

APPENDIX

77

APPENDIX. MATLAB® CODES

%%°/o°/o°/o°/o°/o°/o%°/o%%%%°/o%%%%%%%%%%%%%%%%°/o%%%%%%°/o°/o°/o°/o°/o°/o°/o
% 200 m long shot gather along PZZ transect line
% WSPZD = raw data

% WSPZD] = raw data with first arrival at zero
% WSPZfiltx = dewowed data

% WSPZagx = AGC gain on data
fid=fopen('WSPZl_50.SGY','r','l');
[WSPZD,H]=readsegy(fid);

fclose(fid);

WSPZdt=H(l).dt*le-12;
[WSPZnt,WSPan]=size(WSPZD);

%%%%%%%%%%%%%%%%%%%%%%%%°/o°/o%%%%%°/o%%%%%%%%%%%%%
%Removal of bad traces

btrl =(WSPZD(:,256)+WSPZD(:,258))/2;

WSPZD=[(WSPZD(:,1:256)) btrl (WSPZD(:,258:WSPan))];

%%%%%%%%%%%%%%%%%%%°/o%%%%%%%%%%%%%%%%%°/o°/o%°/o°/o°/o%°/o
%Set data to time zero. Time zero is obtained from *.hd file

WSPZDl =WSPZD(44:WSPZnt, z);

[WSPZntl ,WSPan]=size(WSPZDl );

WSPZtl=( l :WSPZntl)*WSPZdt;

WSPle =(O:O.5 :(WSPan*O.5));

%%°/o%%%%%%%%%%%%%%%%%%%%°/o°/o°/o%°/o%%%%%%%%%%°/o°/o°/o°/o°/o%%
%Frequency Spectrum

trace=WSPZDl (:,l );

subplot(2,l,l), plot(WSPZtl,trace), axis([-O.5e-7 5e-7 -4e4 4e4])

title('First trace from GPR transect'), xlabel('Time(s)')
frecLspec=(fft(trace)).*(conj(ffi(trace)));%Frequency Spectrum

df=l/(length(WSPZtl )*WSPZdt);%Delta frequency for x-axis

N=df"‘length(WSPZtl);

freq=0zdf:N-l;

subplot(2,l ,2), plot( freq,freq_spec)

axis([O 5.5e8 -O.5ell 5.5el 1])

title('Power spectrum'), xlabel('Frequency (Hz)')
%%°/o%%%%%%%%%%%%°/o%%%°/o%%%%%%%%%%%%%°/o%%%°/o°/o°/o°/o°/o°/o°/o%
%Filter data

function filtx=bandpass(xx,dt,flow,filigh,n,ps)
%%°/o%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%°/o°/o°/o°/o%%°/o°/o%
% function filtx=bandpass(xx,dt,flow,fhigh,n,ps)

%

% This function filters an array of data xx

°/o with sampling rate of dt with a bandpass filter

% xx-data size[nt,nx]

78

% fiow- lowfrequency (min frequency, Hz)

% fhigh- high frequency (max frequency, Hz)

% n - number of points in filter (optional input)

% ps- phase=0 for zero phase input, (optional input)

% the filter coeficients are calculated with firl ,

% with n point filter. default n=128, if It not specified on input.

% The filtering is done with ffifilt, this introduces a delay as it

% is not a zero-phase filter. delay is ~ O.5*n time points. for zero phase
% filter insert the optional parameter ps=O.

% Written by Ran Bachrach
%%%%%%%%%%%°/o%%%%%%%%%%%%%°/o°/o%%%%%%%%°/o°/o°/o°/o°/o°/o°/o%%°/o

[nt,nx]=size(xx);

if nargin<5, n=l 28; end;

if nargin<6 ps=l; end;

fN=0.5/dt; %nyquist frequency

fl =flow/fN;

f2=fhighlfN;

df=l/(dt*nt);

B=firl (n,[fl ,f2]);

if ps==

filtx=ffifilt(B,xx); %introduces linear phase delay
elseif ps==0

filtx=filtfilt(B,l,xx); %zero phase filter
end

%%°/o%%%°/o%%%%°/o%%%%%%%%%%%%%°/o%%%°/o°/o%%%%%%%%%%%%%%
WSPZfiltx=bandpass(WSPZDl ,WSPZdt, l 9e6,250e6, 1 28,0);
%%%%%%%%%%%%%%%%%°/o%%%%°/o°/o°/o°/o°/o%%°/o%%%%%°/o°/o°/o°A)°/o%°/o°/o"/o°/o

%Apply AGC to data

function agx=agc(xx,t,dt,opt);

% function agx=AGC(xx,t,dt,opt);

% This function calculates the smooth AGC display of the data.

% xx - data matrix size(nt,nx) (Can be shot gather or cdp or just one trace)

% t - time length of the age operator

% dt - sampling rate, or t=#of samples to be averaged dt=l.

% opt- normalization option. opt=2 gives RMS; opt=l normalize by sum of abs values. Defualt
is RMS

%

% see Yilmaz for details.
%%%%%%%%%%%%%%%%%%%%%°/o°/o°/o%%°/o%°/o°/o%%°/o°A)°/o°/o%%°/o°/o%°/o%%°/o
%Written by R. Bachrach

%Modified by T. Mukerji

if nargin==3 opt=2; end

ng=fix(t/dt);

[nt,ny]=size(xx);

for i=1 :fix((nt/ng))
v=[ng*(i-l )+l :ng*i];

79

if opt=2
M(i,:)=sqrt((sum(xx(v,:)."2))/ng);
elseif opt==
M(i,:)=(sum(abs(xx(v,:))))/ng;
end
end
xl=[l; [ng/2:ng:ng*fix(nt/ng)]'; nt]; Ml=[M(l,:); M; M(end,:)]; xl=[l :nt]';
Ml=interpl(x l ,M l ,xI);
agx=xx./Ml;

%%%%%%%%%%%%°/o%%%%%%%%%%%%%°/o%%%°/o°/o°/o%°/o°/o°/o%°/o%°/o%°/o°/o%
WSPZagx=agc(WSPZfiltx, l 2,WSPZdt);

%%%%%%°/o%%°/o%°/o°/o°/o%°/o%%%%%%%%%%°/o%%%%%%%%%%%%%%%%%%
%Velocity Analysis

function
[spectrum,vel,Vstk,tstk,Vdix,v,lyrthk]=velocityspectrum2(agx,offset,tl ,dt,vmin,vmax,dv,mutestr
ch)

°/o fimction
[spectrum,vel,Vdix,v,lyrthk]=velocityspectrum2(agx,offset,tl ,dt,vmin,vmax,dv,mutestrch);
%

% velocityspectrum calculates the velocity spectrum of a CDP gather
% It uses NMO for calculating constant velocity NMO

% Then, the user picks the velocity function using interactive mode.
% Also uses t2dvf to convert the velocity function to depth

%

% spectrum - Velocity spectrum pannel

% vel - axis of velocity: vmin:dv:vmax

% Vdix - interval velocities

% v - velocity vector (Iength(tl ))

°/o lyrthk - thickness(m) of layers differentiated by velocities.

% vector size is dependent on the number of velocities selected

% during formation of the velocity spectrum

% agx- CDP gather (size [nt,nx])

% offset- distance vector of antenna offset (size ([ l ,nx]))

% tl- time vector (size ([nt,l]))

% dt- sampling interval

% vmin,vmax,dv- velocity values for NMO input

% mutestrch- optional parameter for strech mute. Default = 35%

% see NMO for further details

% written by John MossJuly 2| , 2000

% Modified by Ran Bachrach July 25, 2000
%%%%%%%%%°/o%°/o%°/o°/o%%%%°/o%%"/o%°/o°/o°/o%°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o
if nargin<8;mutestrch=35;end;

temp=2; %used later for 'while' statement

[nt,nx]=size(agx);

%°/o%°/o°/o%%°/o%%%°/o%°/o%%°/o%°/o%°/o%%°/o°/o°/o°/o%%%

figure( 1 ), subplot( l ,2,1 ), imagesc(offset,tl ,agx)

80

inpopt=questdlg('Would you like to remove data to preserve small
offsets?','velocityspectrum','No');
switch inpopt
case 'Yes'
defans= {'10'}; prompt: {'Select x-intercept (integer).'};
dlgtitle='Velocityspectrum';
answer=inputdlg(prompt,dlgtitle, l ,defans);
xint=str2num(answer{ l });
%Create line that truncatcs agx to preserve small offsets
line=offset.*(max(t l )/xint);
for j=l :length(offset);
if offset(j)=xint;
temp=j;
end
end
slope=(max(tl )-t l ( l ))/(temp-l );
line=l :temp;
line(l)=tl(l);
for j=2ztemp;
line(j)=tl(l)+(slope*(i-l));
end
subplot(l ,2, 1 ),hold on,plot(offset(l :length(line)),line,'r')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Creates new agx data set with large offset data deleted
for i=1 :length(line);
for j=l :length(tl );
if offset(i)>xint;
agx(:,i)=0;
end
if tl(j)<line(i);
agx(j,i)=0;
end
end
end
if length(line)<length(offset); %if xintercept<max offset, length of line < length offset
temp=length(offset)-length(line);
for j=l :temp;
agx(:,nx-j+l F0;
end
end
%%%%%%%%%%%%%%%%%%%%%%%°/o%%
subplot(l ,2,2),imagesc(offset,tl ,agx)
case 'No',end
%°/o%%%%%%%%%%%%%%%%°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o°/o%°/o°/o%%%%°/o%%%%%%%
vel=vmin:dv:vmax;nv=length(vel);
[nt,nx]=size(agx);spectrum=zeros(nt,nv);
for j=1 :nv;
[ZZ]=NMO(agx,offset,0,dt,vel(j),0,mutestrch);

%°/o%%%%%%%%%%%%°/o°/o%%°/o%%%%°/o°/o%%°/o°/o%°/o°/o°/o°/o%%°/o%%°/o°/o°/o°/o%%
function [22] =NMO(Tr,offset,tO,dt,v,flag,strch);

81

%function [ZZ] =NMO(Tr,offset,tO.dt,v.flag,strch);

%

% Input Tr: data to NMO (if flag=0) or spray (if flag=l );
% [nt,nx]=size(T r)

% v: velocity scalar for constant velocity NMO

% or vector(size [lznt]') for vriable velocity NMO .

% offset: Vector with offset values. size [I :nx].

% tO- two way travel time till lst live point (use 0 if no delay)

%

% dt: sampling rate

% flag to NMO correct=0 (flatten hyperbula)

% flag to spray (or un-NMO) =1

% strch: Optional Mute strech limit (%):

%

% 22: output data after NMO or un-NMO

% NMO uses the function nmoOv.

%

% NOTE: due to NMO stretch the amplitudes are not consistant from input to output.
% best use AGC after NMO correction. See Jon's BEI for details.
%

%

% Modified by Ran Bachrach, 4-12-2000.

[nt,nx]=size(Tr);

if length(v)==l ;V=v+0*[ l :nt]';else;V=v;end;
h=waitbar(0,'NMO: Please wait...');

for ix=l :length(offset)

%ix

tr=Tr(:,ix);
ZZ(:,ix)=nmoLv(tr,V,t0,dt,offset(ix),flag,strch);
waitbar(ix/length(offset),h);

end

close(h);

function zz=nmoLv(Tr,V,tO,dt,x.flag,strch,opt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function zz=nmoO(Tr,V,tO,dt,x,flag,);

%

%This takes an trace and NMO correct it or spary a zero offset trace into an hyperbula
%

% Tr: Trace at distance x

% t0- two way travel time at x=O till first live point(use 0 if ther is nodelay)

% If you don'nt know what to use, use [lznt]'*dt.

% V- velocity (size(Tr))

% dt- Sampling rate

% zz- nmo'ed trace

% flag- 0 for nmo correction (make hyperbula flat), 1 for spraying trace

% strch -optional parameter NMO streach (%). above this streach coef. the data

% is muted. default is set to no strech mute.

82

% opt- large offset correction:

% Note: Due to NMO streach the amplitude are not consistant from input/output

% Best to use age after nmo correction. See Jon's BEI for details
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if nargin<7;strch=0;end

n=length(T r);
zz=Tr*O;

iz=[l :n]';
z=t0+dt*(iz-l );
%size((x./V(iz))),size(z)
tx=sqn(z."2+(xJV(iz))/‘2);
%detx-tO;

if strch:—=0

assrv=tx*0+l;

else

assr=txj(z+dt/le3); %alowable sample streatch ratio-for NMO strech limit
%

assrv=(assr-l)<=(strch/100);% mute above the strch value

end

itl=l+round((tx-tO)/dt);% Round to nearest neighbor
it: (it1<n).*itl+1;
%
if flag=1%undo NMO
zz(it)=zz(it)+Tr(iz);
else % Do NMO<=>Hyperbula_flat
%
zz(iz)=zz(iz)+Tr(it).*assrv(it):
end

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

spectrum(:,j)=sum(ZZ,2);
end
figure, subplot( 1,4,1 ), imagesc(offset,t l ,agx),grid
title('AGC CMP'), xlabel('Offset (m)'), ylabel('I‘ime (sec)')
subplot( 1 ,4,2), imagesc(vel,t l ,spectrum),contour(vel,tl ,abs(spectrum),15)
colorbar,xlabel('Velocity (rn/s)'),grid
set(gca,'ydir','reverse'),title('Velocity Spectrum')
subplot( 1 ,4,3),plot(max(abs(spectrum),[],2),t l ,'+-'),grid,set(gca,'ydir','reverse')
axis([min(max(abs(spectrum)))-10 max(max(abs(spectrum)))+10 min(tl) max(t1)])
%axis([min(templ)-IO max(templ)+lO min(tl) le-7])
title('Max power'),xlabel('Power')
temp2=tl;
for j=l :length(t l );

for i=1 :length(vel);

if spectrum(j,i)==max(spectrum(j.:)); %if spectrum valuezmax value of spectrum row

83

temp2(i)=vel(i); %identifies velocity value at max spectrum value along each row
end
end
end
subplot( l ,4,4),plot(temp2,tl ,'+-'),grid,set(gca,'ydir','reverse')
axis([min(vel)-(dv*10) max(vel)+(dv*10) min(tl) max(t1)])
%axis([min(vel)-(dv*10) max(vel)+(dv*lO) min(tl) le-7])
title('Max velocity'),xlabel('Velocity (m/s)')

inpopt=questdlg('Choose Stacking Velocity?','velocityspectrum','No');
switch inpopt
case 'Yes'
while temp>l ,
defans={'2'}; prompt={'How many velocity values do you need (integer)?'};
dlgtitle='Velocityspectrum';
answer=inputdlg(prompt,dlgtitle, 1 ,defans);
nmv=str2num(answer{ l });

hdlg=msgbox('Use Mouse for selecting your specified number of points to create a velocity
function.',...
'velocityspectrum'); waitfor(hdlg);
figure(2),subplot( l ,4,4), [Vstk,tstk]=ginput(nmv);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on, plot(Vstk,tstk,'r')
v=l :length(tl); v=v.*0; Vdix=l :nmv; Vdix=Vdix.*0; Vdix( l )=Vstk( l );
for i=22nmv;
Vdix(i)=sqrt(((Vstk(i)’\2*tstk(i))-(Vstk(i- l )"2*tstk(i-l )))/(tstk(i)-tstk(i- l )));
end
for j=1:nmv;
for i=1 :round(tstk(j)/dt);
if v(i)==0
v(i)=Vdix(i);
end, end, end
for j=round(tstk(nmv)/dt):length(v);
v(j)=Vdix(nmv);
end
figure(3),hold off, plot(v,tl ,‘r'), set(gca,'ydir','reverse')
hold on, plot(Vstk,tstk)

inpopt=questdlg('Would you like to make a new velocity spectrum?','velocityspectrum','No');
switch inpopt
case 'Yes'
temp=2;
case 'No'
temp=0;
end
end

case 'No'

Vstk=[];tstk=[];
return

84

end

figure,plot(v,t l ,'r’), set(gca,'ydir’,'reverse')
hold on, plot(Vstk,tstk)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%°/o°/o%°/o°/o°/o%°/o%%%%%%°/o°/o%
% Correct for spherical divergence (v*t)
WSPZfiltx=zeros(size(WSPZfiltx, l ),size(WSPZfiItx,2));
for j=l :size(WSPZfiltx,2);
for k=l :wtable_zero_notopo(j);
WSPZfiltx(k,j)=WSPZfiltx(k,j )*((WSPZt l (k)/2)*(l .353e8));
end
for k=wtable_zero_notopo(i)+l :size(WSPZfiltx, 1 );
WSPZfiltx(k,j )=WSPZfiltx(k,i)*(((wtable_zero_notopo_t(j )/2)*(1 .353e8))+(((WSPZt l (k)-
wtable_zero_notopo_t(j))/2)*(6.78e7)));
end
end

%%%%%%%%%%%%%%%%%°/o%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate amplitude of reflections from sum of reflection energy
for j=l :401;
ampwt_sum(i)=sum(abs(WSPZfiltx6_2(wtable__zero_notopo(j )-6 :wtable_zero_notopo(j )+7,j )));
ampcl_sum(i )=sum(abs(WSPZfiltx6_2(clayreflector_zero_notopo(j )-
5:clayreflector_zero_notopo(j)+5,j)));
end
figure
subplot(2,l ,l), plot(WSPle ,ampwt_sum), title('Amplitude of Water Table')
xlabel('Distance (meters)'), ylabel('Summation of Reflection Amplitude')
subplot(2,l,2), plot(WSPle ,ampcl_sum), title('Amplitude of Clay Layer')
xlabel('Distance (meters)'), ylabel('Summation of Reflection Amplitude')

% Fit polynomials to plots of amplitude values

P=polyfit((l : l 60),ampwt__sum(l : l 60),4);
poly_ampwt_west_sum=polyval(P,(l : l 60));
P=polyfit((245:401),ampwt_sum(245:401),l );
poly_ampwt_east_sum=polyval(P,(245 :401 ));

P=polyfit(1 :401 ,ampcl_sum,4);

poly_ampcl_sum=polyval(P,(l :401 ));

subplot(2,l ,1), hold on, plot(WSPle(l :160),poly_ampwt_west_sum,'r')
plot(WSPle (245:401 ),poly_ampwt_east_sum,'r')

subplot(2,l ,2), hold on, plot(WSPle,poly_ampcl_sum,'r')

% Calculate attenuation using polyfit of Water table to polyfit of clay layer
for j=l : l 60;
alpha_west_sum(j )=(log(poly__ampcl_sum(j )/poly_ampwt_west_sum(j )))/...
(-1 *(6.78e7*((clayreflector_zero_notopo_t(j )-wtable_zero_notopo_t(j ))/2)));
end
forj=l : 157;
alpha__east_sum(i)=(log(poly_ampcl_sum(j+244)/poly_ampwt_east_sum(j )))/...
(-l *(6.78e7*((clayreflector_zero_notopo__t(j +244)-wtable_zero_notopo_t(j+244))/2)));
end

85

figure

plot(WSPZx l (l : l 60),alpha_west_sum)

hold on

plot(WSPle(245:401),alpha_east_sum), title('Attenuation of EM Wave')
xlabel('Distance (meters)'), ylabel('Attenuation (dB/m)')

% Calculate delta alpha across the transect. Delta alpha is used later to calculate delta sigma
delta_alpha_west_sum=zeros( l ,160);
delta_alpha_west_sum(l )=0;
for j=2: l 60;
delta_alpha_west_sumfi )=alpha_west_sum(j )-alpha_west_sum( l );
end
delta_alpha_east_sum=zeros(1 , l 57);
delta_alpha_east_sum( 1 )=0;
forj=2z 157;
delta_alpha_east_sum(j )=alpha_east_sum(j)-alpha_east_sum( l );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Then run 'sigma_fwd' to calculate delta_sigma and sigma using the calculated delta_alpha
value(s)
%%%%%%%%%%%%%%°/o%%%%%%°/o%%%%%%%%%%%%%%%%%%%%%%%
[sigma__west_sum,delta_sigma_west__sum]=sigma_fwd(50,6.7867,delta_alpha_west_sum,0.0l 8);
[sigma_east_sum,delta_sigma_east_sum]=sigma_fwd(50,6.78e7,delta_alpha_east_sum,0.0423);
sigma_west_sum=sigma_west_sum.*((1e6/l 00)*(0.3"(-2.35)));
sigma_east_sum=sigma_east_sum. *(( 1 e6/ 1 00)*(0.3"(-2.35)));

figure, plot(WSPZx l ( l : l 60),sigma_west_sum)

hold on, plot(WSPZx] (245:401 ),sigma_east_sum)

title('Calculated Specific Conductivity of Groundwater')

xlabel('Distance (meters)'), ylabel('Conductivity (uS/cm)')

axis([0 200 0 700])

figure, plot(WSPZx l (l : 1 60),sigma_west_sum)

hold on, plot(WSPZx] (245 :40 l ),sigma_east_sum)

title('Calculated and Measured Specific Conductivity of Groundwater')

xlabel('Distance (meters)'), ylabel('Conductivity (uS/cm)')

axis([O 200 0 705])

hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [sigma,delta_sigma]=sigma_fwd(f,v,delta_alpha,ini_cond);
%function [sigma,delta_sigma]=sigma_fwd(f,v,delta_alpha);

% f = frequency (MHz)

% v = velocity (m/s)

% delta_alpha = vector obtained from alpha_eq.m

% sigma = conductivity (S/m)

% delta_sigma = change in conductivity (S/m)

°/o ini_cond = initial (background) conductivity (S/m)

mu=4e-7*pi;

omega=2 *pi*t"‘ 1 e6;

epsO=l/(((3e8)"2)*mu);

kapa=(3e8/v)"2;

86

eps=kapa*epsO;

sigma( l )=ini_cond/(0.3"(-2.35));
for j=l :length(delta__alpha)-l ;
deriv_kim(i)=((sqrt((eps*mu)/2))*sigma(j))/...
(2*eps"2 *omega*(sqrt((sqrt( l +((sigma(j )/(eps*omega))"2)))-l ))*...
(SqrtU+((Sigma(i)/(¢PS*0mega))"2))));
sigma(jH )=sigma( l )+(delta_alpha(j+l )/deri v_kim(j ));
end

delta_sigma( 1 )=0;
for j=l :length(sigma)-l;

delta_sigma(j+l )=sigma(j+ l )—sigma( l );
end

%%%%%%%%%%%%%%%%%°/o%%%%%°/o%°/o%%%%%%%%%%%%%%%%%%%
% June 2000 data

% Weighted avg of sigma measured

plot(WSPle (100), l 90.707,'r*')

plot(WSPle(l28),298.388,'r*')

plot(WSPZx l ( l 60),530.028,'r*')

plot(WSPZx 1 ( l 87),400.353,'r*')

plot(WSPZx] (216),453.245,'r*')

plot(WSPZx](246),388.075,'r*')

plot(WSPZx] (278),427.342,'r*')

% May 1999 Data

% Weighted avg of sigma measured
plot(WSPZx] (160),3 l 6.71 ,‘rs')
plot(WSPle ( l 87),447.568,'rs')
plot(WSPle(216),672.484,'rs')
plot(WSPZx](246),469.482,'rs')

87

REFERENCES

88

REFERENCES

Annan, A. P. and Chua, L. T., 1992. Ground Penetrating Radar Performance Predictions.
Ground Penetrating Radar, edited by Pilon, J. Geologic Survey of Canada, Paper 90—4,
pp. 5-13.

Annan, A. P., September 2001. Ground Penetrating Radar Workshop Notes, Sensors &
Sofiware Inc.

Atekwana, E. A., Sauck, W. A., and Werkema, D. D., 2000. Investigations of
Geoelectrical Signatures at a Hydrocarbon Contaminated Site. Journal of Applied
Geophysics, Vol. 44, pp. 167-180.

Atekwana, E. A., Sauck, W. A., and Werkema, D. D., 1988. Characterization of a
Complex Refinery Groundwater Contamination Plume Using Multiple Geoelectric
Methods. Proceedings of the Symposium on the Application of Geophysics to
Engineering and Environmental Problems (SA GEEP 1998), pp. 675-683.

Bachrach, R. various Matlab® codes. Michigan State University.

Backe, K. R., Lile, O. B., and Lyomov, S. K., 2001. Characterizing Curing Cement
Slurries by Electrical Conductivity. SPE Drilling & Completion, December, pp. 201-207.

Barber, W. B. and Morey, R., 1994. Radar Detection of Thin Layers of Hydrocarbon
Contamination. Proceedings of the 511. International Conference on Ground Penetrating

Radar, Vol. 1, Waterloo Ventre for Groundwater Research, Waterloo, Ontario, pp. 1215-
1228.

Bermejo, J. L., Sauck, W. A., and Atekwana, E. A., 1997. Geophysical Discovery of a
New LNAPL Plume at the Former Wurtsmith AFB, Oscoda, Michigan. Groundwater
Monitoring & Remediation. Vol. 17, pp. 131-137.

Bigalke, J ., 2000. A Study Concerning the Conductivity of Porous Rock. Physics and
Chemistry of the Earth, Vol. 25, No. 2, pp. 189-194.

89

Boving, T. B. and Grathwohl, P., 2001. Tracer Diffusion Coefficients in Sedimentary
Rocks: Correlation to Porosity and Hydraulic Conductivity. Journal of Contaminant
Hydrology, Vol. 53, pp. 85-100.

Burger, H. R., 1992. Exploration Geophysics of the Shallow Subsurface: New Jersey,
Prentice Hall, pp. 144-164.

Chapelle F.H., Haack, S.K., Adriaens, P., Henry, M.A., and Bradley, P.M., 1996.
Comparison of Eh and H2 Measurements for Delineating Redox Processes in a
Contaminated Aquifer. Environmental Science and Technology, 30: 3565-3569.

Cummings, T. R. and Twenter, F. R., 1986. Assessment of Ground-Water Contamination
at Wurtsmith Air Force Base, Michigan, 1982-85. US. Geological Survey, Water
Resources Investigation Report 86-4188.

Fetter, C. W., 1994. Applied Hydrogeology: Prentice Hall, Upper Saddle River, New
Jersey.

Fitterman, D. V. and Deszcz-Pan, M., 2001. Saltwater Intrusion in Everglades National
Park, Florida Measured by Airborne Electromagnetic Surveys. First International
Conference on Saltwater Intrusion and Coastal Aquifers Monitoring, Modeling, and
Management. Essaouira, Morocco, April 23-25.

Freeze, R. A. and Cherry, J. A., 1979. Groundwater. Prentice Hall, Englewood Cliffs,
New Jersey.

Gueguen, Y. and Palciauskus, V., 1994. Introduction to the Physics of Rocks. Princeton
University Press, Princeton, New Jersey.

Hwang, S., Shin, J., Park, I., and Lee, S., 2004. Assessment of Seawater Intrusion Using
Geophysical Well Logging and Electrical Soundings in 3 Coastal Aquifer, Youngkwang-
gun, Korea. Exploration Geophysics, Vol. 35, pp. 99-104.

Kwader, T., 1986. The Use of Geophysical Logs for Determining Formation Water
Quality. Ground Water, Vol. 24, No. 1, pp. 11-15.

90

Lendvay, J. M., Sauck, W. A., McCormick, M. L., Barcelona, M. J ., Kampbell, D. H.,
Wilson, J. T., and Adriens, P., 1998. Geophysical Characterization, Redox Zonation, and

Contaminant Distribution at a Groundwater/Surface Water Interface. Water Resources
Research. Vol. 34, No. 12, pp. 3545-3559.

McGuire, J. T., Smith, E. W., Long, D. T., Hyndman, D. W., Haack, S. K., Klug, M. J .,
and Velbel, M. A., September 2000. Temporal Variation in Parameters Reflecting

Terminal-Electron-Accepting Processes in an Aquifer Contaminated With Waste Fuel
and Chlorinated Solvents. Chemical Geology, Vol. 169, pp. 471-485.

McMahon, P. B. and Chapelle, F. H., 1991. Microbial Production of Organic Acids in
Aquitard Sediments and its Role in Aquifer Geochemistry. Nature. Vol. 349, pp. 233-
235.

Mukerji, T. Various Matlab® codes. Stanford University.

Nash, M. S., Atekwana, E., and Sauck, W. A., 1997. Geophysical Investigation of
Anomalous Conductivity at a Hydrocarbon Contaminated Site. Proceedings of the
Symposium on the Application of Geophysical Engineering to Environmental Problems
(SA GEEP ’97), Environmental and Engineering Geophysical Society, Englewood, CO,
pp. 675-683.

Reynolds, J. M., 1997. An Introduction to Applied and Environmental Geophysics.
Chichester, England, John Wiley & Sons Ltd., pp. 355-357.

Sauck, W. A., 2000. A Model For The Resistivity Structure of LNAPL Plumes and Their
Environs in Sandy Sediments. Journal of Applied Geophysics, Vol. 44, pp. 151-165.

Sauck, W. A., 1998. A Conceptual Model For The Geoelectrical Response Of LNAPL
Plumes in Granular Sediments. Proceedings of the Symposium on the Application of
Geophysics to Engineering and Environmental Problems (SA GEEP '98), Chicago, IL,
pp. 805-817.

Sauck, W. A., Atekwana, E. A., and Nash, M. S., 1998. High Conductivities Associated
with an LNAPL Plume Imaged by Integrated Geophysical Techniques. Journal of
Environmental and Engineering Geophysics. Vol. 2, pp. 203-212.

91

Sensors & Software, 1996. Pulse Ekko User’s Guide, Version 1.2.

Sheriff, R. E. and Geldart, L. P., 1995. Exploration Seismology. Cambridge University
Press.

Stark, J. R., Cummings, T. R., and Twenter, F. R., 1983. Ground-Water Contamination
at Wurtsmith Air Force Base, Michigan. US. Geological Survey, Water Resources
Investigation Report 83-4002.

United States Geological Survey, Water Resources Division, Lansing, Michigan,
September 1990. Installation Restoration Program, Phase II — Confirmation /
Quantification Stage 2. Final Report for Wurtsmith Air Force Base, Michigan:
Hydrogeology Near Wurtsmith Air Force Base, Michigan.

Yilmaz, 0., 1987. Seismic Data Processing: Society of Exploration Geophysics, Tulsa,
OK.

YSI Incorporated. YSI 6-Series Environmental Monitoring Systems. pp. 171-172.

92

1111111111111111