ATTENUATION ANALYSIS OF GROUND PENETRATING RADAR DATA ACQUIRED OVER A CRUDE OIL SPILL
by
Hydrocarbon spills in the subsurface are a widespread problem and have been extensively studied with ground penetrating radar (GPR). A change in apparent amplitude response has been repeatedly observed over areas of known spills. This research focuses on determining the dominant attenuation mechanism in a multi-offset GPR dataset acquired at a crude oil spill site outside of Bemidji, Minnesota.
A frequency spectrum attribute analysis was performed to test if particular spectra attributes (peak frequency, half-amplitude width, and low and high frequency slopes) would indicate the presence of oil. Several factors were applied to the data to account for non-attenuation related changes in amplitude: radiation pattern, antenna-ground normalization, and geometrical divergence. This analysis exhibited no consistent trends that would indicate the location and presence of the oil if it were not already known.
In the final analysis, attenuation and loss tangent were calculated after additionally taking into account the water table reflection coefficient and the estimated transmitted wavelet’s
amplitude spectrum. The analysis was performed on several, selected reflectors and all exhibited trends that are indicative of conductive losses over the low range of calculated frequencies and scattering losses over the high range. The loss tangents of the contaminated water table reflectors shift upward compared to the uncontaminated water table reflectors’ loss tangents at lower frequencies, indicating an increase in the conductive losses. This increase is supported by specific conductance measurements made on water samples and by a DC resistivity survey. An EM-31 dataset, however, contains less than 1% variation in conductivity over the site. If these data are considered, then the shift upward of the contaminated loss tangent at low frequencies may be due to a decrease in the permittivity caused by the oil replacing the water in the pores. The lack of a similar shift upward at the higher frequencies may indicate a decrease in the amount of scattering from a compressed capillary fringe. To better determine the dominant mechanism causing the observed shift upward in the loss tangents, this same analysis could be performed on a dataset acquired with an antenna of a different frequency, or a multi-polarization GPR survey could be acquired to exploit the depolarization caused by scattering at the capillary fringe.
ABSTRACT... iii
LIST OF FIGURES ... vii
LIST OF TABLES... xi
ACKNOWLEDGMENTS ... xii
INTRODUCTION ...1
Chapter 1 GROUND PENETRATING RADAR ...5
1.1 Theory...5
1.2 Velocity Analysis...8
1.3 Attenuation Analysis...10
1.3.1 Geometrical Divergence ...12
1.3.2 Angular Dependent Gain Functions...13
1.3.3 Amplitude Spectrum of Transmitted Wavelet ...13
1.3.4 Reflection Coefficients ...14
1.3.5 Electromagnetic Loss Mechanisms...17
1.3.6 Scattering Losses ...19
1.3.6.1 Surface Scattering ...20
1.3.6.2 Volume Scattering ...23
1.3.7 Dominant Mechanism Determination...26
Chapter 2 BEMIDJI RESEARCH SITE DESCRIPTION...30
2.1 Site History ...30
2.2 Geologic and Hydrgeologic Setting...31
2.3 Well and Core Data...33
2.4 Previous Geophysical Surveys Performed...41
2.4.1 GPR Surveys...41
2.4.2 Electromagnetic Surveys ...41
Chapter 3 GPR DATA ACQUISITION AND PRE-PROCESSING ...44
3.1 Data Acquisition ...44
3.2 Data Pre-processing ...47
Chapter 4 VELOCITY ANALYSIS ...54
4.1 Semblance Plot and Constant Velocity Stack Analyses ...55
4.2 Constant Velocity Panel Analysis...57
4.3 Complex Trace Attribute Analysis on Stolt Migrated Data...65
4.4 Discussion...69
Chapter 5 ATTENUATION ANALYSIS...73
5.1 Frequency Spectra Attribute Analysis ...74
5.1.1 Non-Attenuation Related Factors Affecting Amplitude ...74
5.1.1.1 Angular Dependent Gain Functions...76
5.1.1.2 Antenna-ground Coupling Normalization ...78
5.1.1.3 Geometrical Divergence ...79
5.1.1.4 Reflection Coefficients and Transmitted Wavelet Spectrum...81
5.1.2 Attribute Analysis Methodology...83
5.1.3 Results...86
5.1.4 Discussion...89
5.2 Attenuation and Loss Tangent Determination ...97
5.2.1 Amplitude Spectrum of Transmitted Wavelet ...97
5.2.1.1 Airwave and Groundwave Separation ...98
5.2.4 Results...114
5.2.5 Discussion...123
Chapter 6 SUMMARY AND CONCLUSIONS...131
REFERENCES CITED...136
Appendix A RADIATION PATTERN ANALYSIS...142
Appendix B MATLAB CODES AND SU SCRIPTS...170
Appendix C RAW AND BINNED GPR DATA ...174
1.1 Plots of reflection coefficients for a perpendicular-polarized electric field ...18 1.2 Surface scattering patterns and a schematic of volume scattering ...21 1.3 Graph of normalized radar cross section for a single spherical particle versus
normalized incident wavelength...25 1.4 Plot of conductive and volume scattering losses in terms of loss tangents as a function of frequency ...28 1.5 Flow chart of the general steps necessary to perform the velocity and attenuation
analyses...29 2.1 1991 aerial photograph of the Bemidji research site with the spill-related
features overlaid...32 2.2 Map of the well positions at the Bemidji site ...34 2.3 Specific conductances and total dissolved solids (TDS) of water samples
collected in June 1987 within the North Oil Pool...37 2.4 Topographic map of the Bemidji site showing the positions of the North Pool
(A-A’) and South Pool (B-B’) transects shown in Figure 2.5 ...38 2.5 Subsurface oil distributions of the North and South Oil Pools ...39 2.6 Three parallel, single antenna, 80 MHz profiles acquired over the spray zone
and South Oil Pool in January 1985 ...42 3.1 1991 aerial photograph with spill site features and positions of the two GPR
profiles acquired in January 2000 overlaid...45 3.2 Photograph of the acquisition of the W602 Line...46 3.3 South Pool and W602 Lines 0.60 m offset profiles...49
drift correction...51
3.6 An expanded pre-processing flow chart of processes performed on data before velocity or attenuation analyses...53
4.1 Expanded velocity analysis flow chart of processes performed on the data ...56
4.2 Sample South Pool Line CMPs ...58
4.3 Sample constant velocity panel from W602 Line with pick velocities ...60
4.4 South Pool Line stacking velocity profile from the constant velocity panel analysis ...61
4.5 W602 Line stacking velocity profile from the constant velocity panel analysis...62
4.6 South Pool Line smoothed stacking velocity profile...63
4.7 W602 Line smoothed stacking velocity profile ...64
4.8 South Pool Line interval velocity profile ...66
4.9 W602 Line interval velocity profile ...67
4.10 Sample instantaneous amplitude plots of a section of the W602 Line 0.60 m offset Stolt migrated data...70
5.1 Expanded flow chart of the steps involved in the frequency spectrum attribute analysis ...75
5.2 Field radiation pattern in relative amplitude versus incidence angle for φ=90°...77
5.3 Stacking velocity profiles and corresponding geometrical divergence functions for the South Pool and W602 Lines ...80
5.5 Schematic of a frequency spectrum...84 5.6 Schematic of a sample trace showing time-windows...87 5.7 South Pool Line full trace attributes that exhibit trends that correspond with the
position of the oil pool...88 5.8 South Pool Line time-windowed trace frequency spectra half-amplitude width
attribute...90 5.9 South Pool Line time-windowed trace frequency spectra low and high
frequency slope attributes ...91 5.10 South Pool Line 0.60 m offset profile showing positions of time windows 3
through 6...94 5.11 Low and high frequency slopes for time windows 3 through 6 calculated from
the South Pool Line 0.60 m offset profile ...95 5.12 A diagram showing the concept of how the infiltration of oil in the subsurface
affects the RDP contrasts...96 5.13 Three binned, un-zeroed profiles from the W602 Line showing traces from
which air and groundwave wavelets extracted...99 5.14 Transmitted wavelet in air of the GSSI 300 MHz antenna and air and
groundwave wavelets extracted from three traces indicated in Figure 5.13 ...100 5.15 Frequency spectra plots that correspond to the recorded time wavelets shown
in Figure 5.14...101 5.16 Time wavelets and corresponding spectra of the estimated transmitted wavelet
into the ground...104 5.17 0.60 m offset profiles showing the positions of the selected reflectors for the
attenuation and loss tangent analysis...110
attenuation and loss tangent analysis...112 5.20 Expanded flow chart of the steps involved in the attenuation and loss tangent
analysis ...115 5.21 Loss tangent and wavelet spectra plot for trace 870 (sample 59-93) of the
W602 Line 0.60 m offset profile ...117 5.22 Loss tangent and wavelet spectra plot for trace 1778 (sample 111-142) of the
South Pool Line 0.60 m offset profile ...118 5.23 Loss tangent and wavelet spectra plot for trace 1778 (sample 161-196) of the
South Pool Line 0.60 m offset profile ...119 5.24 Loss tangent plots for same three traces on square log-log scale...121 5.25 Model of the first four wavelets to arrive at the receiving antenna...122 5.26 Loss tangents for three reflectors from uncontaminated zones of the South
Pool Line 0.60 m offset profile...124 5.27 Loss tangents for three reflectors from contaminated and uncontaminated
zones of the South Pool Line 0.60 m offset profile ...125 5.28 Loss tangents for three reflectors from contaminated and uncontaminated
zones near the water table of the South Pool Line 0.60 m offset profile...126 5.29 Loss tangents for three reflectors from below the water table in contaminated
and uncontaminated zones of the South Pool Line 0.60 m offset profile ...127
1.1 Relative dielectric permittivities and conductivities of common geologic
materials and pore fluids at 100 MHz...7
3.1 Profile acquisition parameters ...47
4.1 Summary of velocity analysis method parameters...71
5.1 Summary of model parameters for estimated transmitted wavelets...106
5.2 Summary of selected reflectors for attenuation and loss tangent analysis ...109
This research was funded by the U.S. Geological Survey’s (USGS) GPR Fundamental Research Project, which is led by Jeff Phillips of the Crustal Imaging and Characterization Team. The geophysical data acquired at the Bemidji research site were funded by the USGS’s Toxic Substances Hydrology Program. My education was also partially funded by the National Science Foundation’s (NSF) Graduate
Scholarship/Research Fund, granted through the Colorado School of Mines, and by a scholarship from the Society of Exploration Geophysicists (SEG) Foundation.
I would like to thank Jared Abraham, Jeff Lucius, and Dr. Mike Powers from the USGS Crustal Imaging and Characterization Team for acquiring the GPR dataset
analyzed in this research. Jeff Lucius has also been particularly helpful in providing additional information on the Bemidji site and previous work performed.
Appreciation is extended to Dr. Powers and Dr. Tom Boyd for serving as my committee members. I am grateful for the ideas and support that Dr. Powers has provided me throughout this research project. I would also like to extend my sincere gratitude to my advisor, Dr. Gary Olhoeft, for his support and guidance throughout my graduate studies.
INTRODUCTION
Ground water contamination by various chemicals is a widespread problem. The list of contaminants includes hundreds of different compounds with physical properties that vary greatly (Lucius et al., 1992). A particular class of contaminants, non-aqueous phase liquids (NAPLs), are immiscible with water and include many organic compounds. NAPLs that are denser than water are known as dense non-aqueous phase liquids
(DNAPLs) and will continue to migrate downward below the water table due to gravity-driven flow. Examples include creosote, coal tar, and chlorinated hydrocarbons such as trichloroethylene (TCE) and tetrachloroethylene (PCE). NAPLs that are less dense than water are called light non-aqueous liquids (LNAPLs) and will float on top of the water table. Examples include gasoline, diesel fuel, and crude oil (Fetter, 1999).
The use of a geophysical technique called ground penetrating radar (GPR) to investigate soils contaminated with LNAPLs has been extensively researched over the past decade. In particular, research involving hydrocarbon spill sites has generated great interest. Results from field investigations involving natural spills and from controlled-spill experiments have produced contradictory and inconsistent results compared to the expected and observed GPR responses.
GPR responds to changes in the electromagnetic properties of the subsurface. The electromagnetic properties that control GPR propagation are dielectric permittivity (ε ), electrical conductivity (σ ), and magnetic permeability (µ). These properties are defined by the response of a material to an applied electric and magnetic field. In general, materials that conduct electricity well cannot support propagation of GPR energy. Within electrically insulating materials where GPR propagation can occur, it is the variations in the electrical and magnetic properties that create a GPR response.
In an attempt to locate and to characterize the many hydrocarbon spills, GPR is often utilized because of the electrical property contrasts between the soil and the different pore fluids (air, water, and LNAPL). The relative dielectric permittivities (RDP) and conductivities change drastically as the water content in the pores changes. Air and LNAPLs both have low conductivities and RDPs (1 and 2-3, respectively) (Lucius et al., 1992). Natural ground water, however, is much more conductive and has an RDP of 80 (Davis and Annan, 1989).
When an LNAPL is spilled at or near the ground surface, it infiltrates vertically downward through the vadose, or unsaturated, zone until it reaches the water table. Because an LNAPL is less dense than water, it accumulates on top of the water table and will eventually migrate laterally with the topography of the water table. Once in the subsurface, these hydrocarbons become a continuous source of ground water
contamination. Although LNAPLs have low solubilities, large volumes of ground water are contaminated by free-phase product slowly dissolving into the surrounding ground water. This dissolved phase can be a result of contact with the pool itself or with the residual contaminant in the vadose zone above the pool. A hydrocarbon with a high vapor pressure can also volatilize and contaminate the vadose zone as a vapor plume (Fetter, 1999; Domenico and Schwartz, 1998).
Based on current research, GPR is generally not effective at directly detecting LNAPLs in the subsurface. Controlled LNAPL spill experiments have yielded conflicting radar responses at the top of the saturated zone. To explain an observed increase in amplitude at the top of the saturated zone under the spill, Campbell et al. (1996) hypothesized that the introduced LNAPL created a sharper water table with a thin capillary fringe by altering the soil-water-oil wettability. In the case where a decrease in GPR amplitude was observed (Redman et al., 1994), a possible explanation was that the LNAPL acts as a transition, or intermediate, layer between the low-RDP, air-filled pores of the unsaturated zone and the high-RDP, water-filled pores of the saturated zone. The smoothly varying RDP gradient that was generated could not effect a GPR response. In
all of these controlled-spill experiments, however, the surveys were taken within months of the spill and do not reflect what is observed in data taken over field sites where often many years have passed between the spill and the GPR survey. Hydrocarbons age with time. They lose volatiles, oxidize, and change the soil-water-hydrocarbon chemistry, while they also change the soil wettability.
In data acquired at “natural” hydrocarbon spill sites, there tend to be “shadow zones” in areas with contaminant pools that are not present in the controlled-spill experiments. These shadow zones are areas with reduced reflectivity or changes in the reflectivity pattern, from an apparently increased attenuation, as compared to surrounding areas. Because LNAPLs are resistive and therefore thought to be conducive to radar signal propagation, several theories have been proposed to explain this apparent attenuation phenomenon. The first published explanation did not deal with the signal attenuation but rather with a lack of electrical property contrast between the surrounding sand and the LNAPL to effect a reflection (Olhoeft, 1986). More recent theories deal directly with the attenuation but differ in their explanations because of the differences in position of the shadow zones in relation to the water table. Most investigations have noted a shadow zone over the plume in the vadose zone and therefore suggest a
relationship between attenuation and vapor plume effects (Daniels et al., 1992; Grumman and Daniels, 1995; Maxwell and Schmok, 1995). The most current theory, however, is based on data that show the shadow zone directly below the water table in the saturated zone. The explanation is a thin conductive zone surrounding the resistive pool that is attenuating the radar signal. This conductive zone is a result of biodegradation of the oil pool that is in contact with the subsurface for an extended period of time (Sauck et al., 1998). A conductive zone could also be created by leaching of salts from the oil or by flushing salts from clays (Olhoeft, 1992), resulting in a conductive “halo” around the advancing front of the contaminant.
An initial motivation of the research in this thesis was to determine the dominant attenuation mechanism that causes these shadow zones in the presence of LNAPLs. A
multi-fold GPR dataset containing shadow zones acquired at a natural crude oil spill site was used for the analysis. The spill occurred in 1979 when an underground pipeline burst near Bemidji, Minnesota. This particular dataset was acquired in 2000 at what is now a well-characterized, highly monitored U.S. Geological Survey (USGS) research site (Delin et al., 1998).
A detailed velocity analysis was performed on the data in an attempt to gain additional information on the location of the contaminants and to get a more correct depth section. The velocity profiles were also used to account for losses in amplitude in the data due to geometrical divergence. Several seismic-style velocity analysis techniques were used, such as semblance plots and constant velocity stacks, constant velocity panels, and complex trace attribute analysis on Stolt migrated data. A frequency spectrum
analysis was performed to determine if trends in specific spectrum attributes correlated with the presence of oil in the subsurface. This analysis was done on the data after accounting for factors that affect the amplitude and frequency dependence of the wavelet that are not related to attenuation due to material losses. These factors include the antenna radiation pattern, changes in antenna-ground coupling, and geometrical divergence. By stepping through the radar equation, which includes all factors that contribute to the power recorded at the receiving antenna, attenuation as a function of frequency was calculated. The dominant attenuation mechanism determination was then performed by analyzing the frequency dependent nature of the attenuation and estimated loss tangent.
CHAPTER 1
GROUND PENETRATING RADAR
1.1 Theory
A GPR system consists of a timing device that controls transmission of an electromagnetic signal, transmitting and receiving antennas, and a receiver that digitizes and stores the returned signal. The timing device ensures that the transmitting antenna radiates pulses at a constant rate. The resulting transmitted pulses from the antenna have a particular frequency bandwidth with a center frequency that is dependent on the size and other properties (shape, type, etc.) of the antenna. When the signal is both
transmitted and received by the same antenna, the system is called monostatic. A more common mode, however, is known as a bistatic system with separate transmitting and receiving antennas. Bistatic antennas can be contained in the same housing, parallel to one another, or each in its own housing. The returned signal that is detected by the receiving antenna is digitized by the system electronics and stored. Each received signal is called a trace, which is a record of the received electric field amplitude variation with respect to time for a given position of the antennas. Traces can be acquired in constant time increments according to the pulse timing electronics, or in constant distance increments according to a wheel rotation as the antennas are moved along a traverse. Each trace is made up of multiple samples. The total recording time, or time range, is the sampling rate in nanoseconds multiplied by the number of samples per trace. In
commercial systems, traces are digitized and recorded using the equivalent-time sampling technique. This sampling method is used with high frequency waveforms and requires an identical, repetitive signal. For GPR data acquisition, source signals are transmitted at a fixed pulse rate, and a single, sequential sample is recorded from each pulse. The result
is a single trace constructed from a number of transmitted pulses equal to the number of samples (Witte, 2002). For example, a trace with 512 samples is a compilation of 512 separate, transmitted signals into the subsurface that are received at slightly later times with each pulse, then accumulated into one equivalent radar trace. This method of digitization therefore assumes that the GPR system is pulled along the ground at a slow enough rate that a single trace accumulated from many transmitted pulses is a
measurement of a relatively small ground surface area.
When a pulse is sent to the transmitting antenna, an electromagnetic wave is radiated into the ground. As the wave propagates, whenever a contrast in
electromagnetic properties is encountered, some of the energy may be reflected back toward the receiving antenna. The propagating energy is sensitive to dielectric
permittivity (ε ), electrical conductivity (σ ), and magnetic permeability (µ). Dielectric permittivity can be defined as a complex material property that is a measure of the material’s ability to both store and dissipate energy as a function of frequency in the presence of an external electric field. The real part describes energy storage as a result of charge polarization, and the imaginary part describes energy lost due to the movement of charges in the process of polarization in response to the applied external electric field (Olhoeft, 1998). It is most often expressed as a dimensionless, relative value, normalized with respect to the dielectric permittivity of free space,
0 r ε ε ε = (1.1)
where εr is the relative dielectric permittivity,
ε is the absolute dielectric permittivity of the medium [farads/m], and ε0 is the dielectric permittivity of free space,
12
8.854 x10− farads/m.
Free space, or air, therefore has a real relative permittivity of one. Table 1.1 gives the relative permittivity values for some typical geologic materials and pore fluids. Sharma (1997) and Reynolds (1997) contain values for additional media.
Table 1.1 Relative dielectric permittivities and conductivities of common geologic materials and pore fluids at 100 MHz (Davis and Annan, 1989; * Ulaby, 2001).
Material εr σ (mS/m) Air 1 0 Distilled water 80 2x10-3 Fresh water 80 0.1 Sea water 72-80* 1000 Dry sand 3-5 0.01 Saturated sand 20-30 0.03-0.3 Limestone 4-8 0.4-1 Shales 5-15 1-100 Silts 5-30 1-100 Clays 5-40 1-300 Granite 4-6 0.01-1 Dry salt 5-6 0.01-1 Ice 3-4 0.01
Electrical conductivity is a measure of a material’s ability to transport charge, or conduct electricity, and so is always a measure of energy loss for wave propagation. The measure of conductivity used in GPR is a real, direct current (DC) electrical conductivity that describes the charge transport and loss at the lowest, limiting frequency (Olhoeft, 1998). These values are quite variable for a given material (Table 1.1) because of the high dependence on water content, quality, and pore interconnection and on the presence of interstitial, mineralogical clay.
Magnetic permeability is similar to dielectric permittivity in that it is a complex material property that is a measure of energy storage and loss. The difference is that these processes are a result of electron spin responses due to an applied external magnetic field (Olhoeft, 1998). Magnetic permeability is also expressed as a relative value with respect to the permeability of free space,
0 r µ µ µ = (1.2)
where µr is the relative magnetic permeability,
µ is the absolute magnetic permeability of the medium [henries/m], and µ0 is the magnetic permeability of free space, 4 x10π −7henries/m.
Most geologic materials have a relative magnetic permeability value around one. Small variances from a value of one greatly change the recorded waveforms. Because it is difficult to independently determine dielectric permittivity and magnetic permeability of a material without measuring soil samples in a lab using a network analyzer, a magnetic permeability of one is often assumed.
The magnitude of the contrast between these electromagnetic properties determines in part how much energy is reflected, or scattered, at the interface upward toward the receiving antenna. A larger contrast therefore results in larger amplitude in the data. Other factors such as antenna geometry and polarization also affect the
magnitude of the recorded amplitudes. The relative direction of change in contrast (high to low or low to high property values) determines the polarity of the reflection. For additional information on GPR, refer to Davis and Annan (1989), Reynolds (1999), Sharma (1997), and Powers (1995).
1.2 Velocity Analysis
Accurately determining a detailed velocity profile of the subsurface is important for ascertaining the dielectric permittivity and magnetic permeability for each of the layers through which the energy propagates and the depth and thickness of these layers. Most GPR surveys are acquired with a fixed source-receiver offset. Velocities can be estimated from these constant offset profiles with a hyperbolic curve fit to a reflection event in the data. A hyperbolic reflection is a result of a point or 2D scatterer as the antennas are dragged across the surface. The shape and extent of the hyperbola in time
are dependent on the propagation velocity, and therefore material properties, between the surface and the scatterer and of the scatterer itself, as well as the transmitter-receiver offset, the diameter of the reflector, and the trace spacing. RDP and magnetic
permeability determine the velocity of the energy propagation according to the following equation when assuming a low loss material:
r r
c v
µ ε
= (1.3)
where v is the velocity of the medium [m/s] and
c is the velocity of electromagnetic wave propagation in air, 3.0 x 108 m/s. As a rule of thumb, a material is considered low loss if the conductivity is below 30 mS/m. A simplification of Eq. 1.3 assumes that µr is equal to one, so the velocity can be estimated by using onlyεr, as in
r
c v
ε
= . (1.4)
Hyperbola fitting, however, gives velocity information only according to the number and position of edge reflectors. If there is only a single pipe, for example, embedded in a layered medium, then only one velocity can be determined. This velocity is then usually assumed to be the average velocity over the whole dataset. If there are multiple hyperbolas at varying depths, then the velocity profile can be better constrained and a more accurate depth section can be calculated when two-way travel-time is
converted to depth using the velocity model.
With multi-offset GPR data, however, more sophisticated methods can be applied to obtain a better velocity profile that varies both laterally and with depth. Seismic-style methods are often employed, because the radar data can be sorted into multi-offset common midpoint (CMP) gathers. As with seismic data when assuming a horizontally layered subsurface, a single CMP is composed of multiple traces of varying transmitter-receiver offsets that are sampling the same point of a given reflector. The travel-time
curve for a single reflection target is hyperbolic in shape when the traces are ordered from near to far offset. Normal moveout (NMO) is the difference in two-way travel-time from the surface to the reflector between a zero-offset trace and a trace at another offset. An NMO-corrected gather has therefore had all non-zero offset traces shifted upward in time to flatten the reflection across all offsets. The velocity that corresponds to the proper NMO correction can be used to determine the velocity of the medium above the reflector (Yilmaz, 1987).
There is, however, one major difference between seismic and GPR velocity profiles. Seismic velocities tend to increase with depth. The velocity of an
electromagnetic wave in the ground is dependent on µr and εr, as shown in Eq. 1.3. As water content increases with depth, εr also increases because a medium containing water-filled pores will have a higher relative permittivity value than that same medium with air-filled pores. Assuming µr is equal to one, GPR velocity profiles therefore tend to decrease with depth. The maximum velocity of m/s occurs at the surface with the direct wave and decreases by at least half at the first reflection. The velocity decreases much more dramatically if the near surface is even partially saturated or is a finer-grained soil with higher water retention.
8 3.0 x 10
1.3 Attenuation Analysis
Propagating energy spreads out as it travels into the ground. Eventually, the amplitude of the propagating energy is weakened to the point where the local reflected energy no longer overcomes the noise level. This decreased, received energy is due to geometrical divergence and to attenuation mechanisms present in the medium through which the energy travels. Geometrical divergence is not a true loss but a spreading out of a finite amount of energy over an increasingly larger area; it is still electromagnetic
energy. Attenuation mechanisms, however, cause energy loss by converting
electromagnetic energy to heat (Keller, 1987) or by scattering the wave in directions that do not return to the receiver.
True signal attenuation can be determined once other factors that affect recorded amplitudes and frequency response are taken into account. The signal power that is returned to the receiving antenna is a function of several factors and is described by the radar equation, 2 1 2 2 0 2 1 1 ( " ') 4 j j n r r t j r j P P G e K G P P α λ π − − = ⎡ ⎤ ⎛ = ⎢ ⎜ ⋅ ⎝ ⎠ ⎣
∏
⎦ ⎞ ⎥ ⎟ r G (1.5) where is the power received [W], Pris the power emitted from the transmitting antenna [W], 0
P
Gt and are the dimensionless, angular dependent gain functions for the transmitting and receiving antennas, respectively,
2 1
( " ')P P⋅ is the geometrical spreading term taken from May and Hron (1978) and described in more detail by Powers (1995) [1/m4] (This is a typo in the text of Powers (1995), but the term as written here is correct.),
is equivalent to the radar cross section of the scatterer [m 1 2 2 1 j j n r j j e α K − − =
∏
2] ,is the complex reflection or transmission coefficient of the j
j
K th segment,
j
α is the attenuation constant [Np/m],
is the radius of the propagating wavefront in the j
j
r th segment [m],
j is the increment number of ray segment, or layer, being calculated,
n is the total number of layers between the receiving and transmitting antennas, and
The term is the exponential power material loss (Powers, 1995). After accounting for antenna effects (angular dependent gain functions and power spectrum of the
transmitted wavelet into the ground), reflection coefficients, and geometrical divergence, all of which are terms that affect the received power, then remaining losses in power are due solely to attenuation. Because all terms of the radar equation except geometrical divergence are frequency dependent in lossy media, attenuation is therefore also a function of frequency. Although the radar equation describes the power recorded at the receiving antenna, this research focuses on comparing amplitudes of the electric field, which is simply proportional to the square root of power.
2 j jr
e−α
1.3.1 Geometrical Divergence
Geometrical divergence is a result of geometrical spreading as the wave
propagates. A finite amount of transmitted energy is spread over an increasingly larger area as the wavefront propagates downward from the source, like being spread over the increasing surface area of an expanding balloon. This spreading mechanism is a function of the distance traveled in the subsurface and can be compensated for if an accurate velocity profile is known. There are several different functions that can be used to account for this energy loss, all of which are again borrowed from the seismic industry, including the May and Hron (1978) term in Eq. 1.5. A simpler factor, however, will be used for this analysis. In general for a homogeneous earth, amplitudes decay inversely with the radius, , of the propagating wavefront. For a layered earth, the non-frequency dependent correction factor applied to the data becomes (Yilmaz, 1987)
r
2( )
SD
G =v t ⋅ (1.6) t
where v is the root mean square (RMS) velocity [m/s] and
1.3.2 Angular Dependent Gain Functions
The angular dependent gain functions are also known as directivity functions or radiation patterns of the antennas. Radiation patterns describe the relative amplitude of the energy output or received by the antenna as a function of angle or direction. It is a three-dimensional, non-circularly symmetric pattern that is a function of the antenna design. The shape of the pattern is also dependent on the properties of the ground on which the antenna rests. Both the transmitting and receiving antennas have radiation patterns that can each be altered by the interaction of the two patterns in a bistatic system. These patterns can be frequency dependent to the extent that the antennas themselves are frequency dependent. Appendix A includes a full description of the theoretical far-field radiation pattern for an infinitesimal dipole and an extension of the theory to approximate a bowtie antenna. A radiation pattern derived from the multi-offset dataset acquired in Bemidji, Minnesota is also compared to these theoretical patterns. This field radiation pattern is used to account for the angular dependent gains. This process is described in more detail in Section 5.1.1.1.
1.3.3 Amplitude Spectrum of Transmitted Wavelet
The power spectrum as a function of angular frequency of the transmitted wavelet ( ) must also be taken into account to study and accurately quantify attenuation. Again, because this analysis focuses on the amplitudes of the data, the transmitted amplitude spectrum is estimated instead.
0
P
This transmitted wavelet, however, can be difficult to determine from most bistatic GPR data. It is not simply the first arrival wavelet in the data, which is the airwave, because the fastest medium for energy travel is air. Even the arrivals seen
directly below the airwave in the near field have not all traveled in the subsurface between the transmitter and receiver. These arrivals are most likely a combination of energy bouncing off the surface back toward the receiver and “refracted” arrivals that traveled along the ground surface and back into the air to the receiver.
Because the transmitted wavelet into the ground is dependent on the properties of the ground on which it rests, it can change as the ground properties change along a given profile. In the case of a homogeneous ground surface or acquisition over snow, however, the transmitted wavelet can be assumed constant along the profile.
One method to estimate the transmitted wavelet is to model either a Ricker wavelet or a known transmitted wavelet in air from the antenna used in the survey with known properties of the ground into the which the energy is transmitted. A second method is possible if there is a large enough transmitter-receiver offset. In this case, the airwave arrival and first real reflection from the subsurface are separated in time. The first strong reflection below the airwave can be used as an estimate for the transmitted wavelet.
1.3.4 Reflection Coefficients
The reflection and transmission coefficients describe the portions of incident energy on an interface that are reflected upward toward the receiver and that are
transmitted past the interface further into the subsurface as a function of incidence angle. Transmission coefficients are neglected here. The following equation is taken from Lehmann (1996) and describes the reflection coefficient, , for a perpendicular-polarized electric field with respect to the plane of incidence:
refl K 2 1 2 2 1 2 1 1 1 1 1 1 1 1 ( cos ) ( cos ( cos ) ( cos r r r refl r r r q i p K q i p 1 1 ) ) r r µ β θ µ µ α θ µ µ β θ µ µ α θ µ − + − = + + + (1.7)
where the subscripts 1 and 2 refer to the upper and lower layers at the boundary, respectively,
α is the imaginary part of the wavenumber known as the attenuation constant [Np/m],
β is the real part of the wavenumber known as the phase constant [rad/m], 1 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 4 2 2 2 2 2 2 1 1 1 2 2 2 2 2 1 1 1 2 2 1 1 2 2 1 ( ) ( )sin ... 2 ( ) sin ( ) ... 1 2 2sin ( )( ) 4 k k q k k k k k k α α θ α θ α θ α α α α ⎡ ⎡ − − − ⎤+ ⎤ ⎣ ⎦ ⎢ ⎥ ⎢ = ⎢ ⎡ + + + − ⎢ ⎢ ⎥ ⎥ ⎢ ⎢⎣ ⎡⎣ − − + ⎤⎦⎥⎦ ⎥ ⎣ ⎦ ⎥ ⎥ ⎤ , (1.8) 2 2 2 1 1sin 1 k k p q α − α θ = , (1.9)
θ is the vertical angle (angle of incidence) [rads], and k is the wavenumber of the medium [1/m].
In lossy, dispersive materials, the wavenumber is ( r' r'')( 'r r'') ( r' r'')( ) o k i i i i c ω µ µ ε ε µ µ σ ε ω = − − − − (1.10)
where ω is angular frequency [rad/s],
εr' and εr'' are the real and imaginary parts of relative dielectric permittivity, respectively,
µr' and µr'' are the real and imaginary parts of relative magnetic permeability, respectively, and
σ is electrical conductivity [S/m] (Powers, 1995).
Both the attenuation and phase constants are functions of angular frequency. The attenuation constant, α , and phase constant, β, can be expressed in terms of complex electromagnetic material properties as
2 2 2 A B A c ω α =⎛ ⎞⎜ ⎟ ⎡⎢ + − ⎤⎥ ⎝ ⎠ ⎢⎣ ⎥⎦ and (1.11) 2 2 2 A B A c ω β =⎛ ⎞⎜ ⎟ ⎡⎢ + + ⎤ ⎝ ⎠ ⎢⎣ ⎥⎦⎥ (1.12) where ' ' ''( '' ) r r r r o A µ ε µ ε σ ωε = − + and (1.13) '' ' '( '' ) r r r r o B µ ε µ ε σ ωε = + + . (1.14)
A common simplification of Eqs. 1.11 and 1.12 involves neglecting the imaginary parts of both εr and µr and assuming that µr is equal to one (no magnetic materials present). Substituting the simplified versions of Eqs. 1.13 and 1.14 into Eq. 1.11 results in the following expression: 2 2 0 1 2 r c ω σ r α ε ωε ⎡ ⎛ ⎞ ⎤ ⎛ ⎞ ⎢ =⎜ ⎟ + ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ε ⎥ − . (1.15)
The result of these simplifying assumptions on Eq. 1.12 is (Powers, 1995) 2 2 0 1 2 r c ω σ r β ε ωε ⎡ ⎛ ⎞ ⎤ ⎛ ⎞ ⎢ =⎜ ⎟ + ⎜ ⎟ ⎢ ⎝ ⎠ ⎣ ⎝ ⎠ ε ⎦ ⎥ + ⎥ α . (1.16) To convert the attenuation constant from Np/m to dB/m, the expression
(1.17) [ / ]
[dB m/ ] 20 logeαNp m 8.686 [Np m/ ]
α = ⋅ ≅ ⋅
can be used (Balanis, 1989). The phase constant can also be defined by ω
β ν
= (1.18) where ν is the phase velocity [m/s] of angular frequency ω (Powers, 1995).
Reflection coefficients are typically neglected in GPR data, because the data are most often single offset and the contrasts at the majority of reflection interfaces are not sufficient to affect the amplitudes as a function of incident angle. Two interfaces,
however, that can be important to multi-offset datasets are the air-ground and water table boundaries. The reflection coefficients of two hypothetical interfaces were calculated to determine their potential effect on the recorded amplitudes. Figure 1.1 shows plots of the calculated reflection coefficients as a function of incidence angle for air-ground and water table interfaces for varying RDP contrasts. These reflection coefficients assume a smooth, specular interface with no frequency dependence. Frequency dependence, however, will be a factor with fine-grained media and rough interfaces.
1.3.5 Electromagnetic Loss Mechanisms
After taking into account the factors mentioned above, the remaining cause for loss of radar signal amplitude in the data as the energy propagates through different media is attenuation. Attenuation is comprised of several different mechanisms. One such mechanism is conductive loss, or the transfer of electromagnetic energy to heat (Keller, 1987). The loss factor, equal to the ratio of conductivity to the product of angular frequency and dielectric permittivity, is often expressed in terms of a loss tangent, 0 tan ' c r σ δ ωε ε = (1.19)
where tanδc is the conductive loss tangent.
The loss conductive loss tangent physically describes the cotangent of the phase angle between J, the current density vector in A/m2, and E, the electric field vector in V/m (Powers, 1995).
Figure 1.1 Plots of reflection coefficients for a perpendicular-polarized electric field with respect to the plane of incidence as a function of incidence angle for two interfaces. Each plot contains three sample media combinations with varying RDP contrasts. The top plot is of the air-ground interface, and the bottom plot is of the water table interface. The subscripts 1 and 2 refer to the upper and lower layers at the boundary, respectively. A larger reflection coefficient corresponds to a larger percentage of incident energy being reflected upward at the interface.
Another attenuation mechanism is dielectric relaxation. It is a result of charges moving from a non-polarized to a polarized state in response to an applied external field (Olhoeft, 1998). In the case of polar water molecules, energy is lost as the electrically neutral, dipolar molecules rotate in place to align themselves with the applied field in a viscous medium (Olhoeft, 1986). The dielectric loss tangent ( tanδd) is expressed as
'' tan ' r d r ε δ ε = (1.20)
and is equal to the cotangent of the phase angle between D, the displacement current vector in C/m2, and E. The magnetic relaxation losses are related to the complex magnetic permeability of the medium,
" tan ' r m r µ δ µ = (1.21)
where tanδm is the magnetic loss tangent.
The magnetic loss tangent is physically equal to the cotangent of the phase angle between B, the magnetic induction vector in teslas or Wb/m2, and H, the magnetic field vector in A/m. The electrical and magnetic loss tangents can be combined into a single
electromagnetic loss tangent ( tanδem) equal to tanδem α
β
= . (1.22)
Again, the total electromagnetic loss tangent is physically equal to the cotangent of the phase angle between E and H (Powers, 1995) and is frequency dependent.
1.3.6 Scattering Losses
While the attenuation constant takes into account conductive and dielectric and magnetic relaxation losses, scattering can also be a major attenuation mechanism.
Scattering includes all reflected energy, that which is recorded at the receiving antenna and that which is reflected in every other direction, so it is not a fully undesirable mechanism. Too much scattering, however, can be a problem when it is not directed toward the receiving antenna or is a result of the interaction between the dominant wavelength and grain size or surface roughness. There exist both surface and volume scattering. Figure 1.2 shows surface scattering patterns for varying degrees of rough surfaces and a schematic of volume scattering. Rough surface and volume scattering cause both frequency and polarization dependencies.
1.3.6.1 Surface Scattering
Surface scattering occurs at an interface between two media where energy is reflected and transmitted. A rougher, or less specular, surface will cause energy to be reflected in multiple directions and less energy to be reflected in the direction of the receiving antenna. When energy encounters a smooth boundary, the reflection is
specular, like a mirror. A specular surface scatters all reflected energy at an angle equal to the angle of incidence according to the Fresnel reflection laws. The energy that is scattered as a result of specular reflection is the coherent scattering component. As the surface becomes rougher, relative to the dominant wavelength of the incoming energy, the coherent component becomes smaller because some of the energy is scattered at other angles. This part of the scattered energy at other angles is the diffuse scattering
component. If the surface becomes rough enough, the reflected energy is composed only of the diffuse component (Ulaby et al., 1982).
Two parameters that define the degree of roughness of a boundary are the
standard deviation of the surface height variation (rms height) and the surface correlation length. Both parameters are statistical descriptions of the variation between the random surface height component (high frequency) and the background component (periodic, low
(a)
(b)
Figure 1.2 (a) Surface scattering patterns for varying surface roughnesses and (b) a schematic of volume scattering (Ulaby et al., 1982).
frequency or flat). Refer to Ulaby et al. (1982) for an explanation of surface correlation length. The discrete standard deviation of the surface height variation, , assuming a one-dimensional roughness of the interface in the direction of the GPR profile, is defined by h SD 1 2 2 2 1 1 ( ) ( ) 1 N h i i SD z N z N = ⎡ ⎛ ⎤ =⎢ ⎜ − − ⎝ ⎠ ⎣
∑
⎦ ⎞ ⎥ ⎟ (1.23) where N is the number of samples,is the height of the i
i
z th sample above a reference elevation [m], and
z is the mean height of the surface, expressed by
1 1 N i i z N = =
∑
z . (1.24)As a rule of thumb for defining a surface as smooth, the Fraunhofer criterion is used. This criterion requires that the surface be in the far field of the antenna, defined by a phase difference of the emitted energy from the center and edge of the antenna that is less than π/8 rads. A surface is considered smooth when the Fraunhofer criterion is satisfied, as defined by (Ulaby et al., 1982)
32cos
h
SD λ
θ
< (1.25)
where λis the dominant, incident wavelength [m].
Quantifying surface, as well as volume, scattering in the subsurface is difficult. To estimate the loss, however, equations describing atmospheric scattering and
attenuation mechanisms of airborne radar due to particulates and clouds are borrowed. Scattering is most often approximated by calculating the radar cross section, σ, equivalent to the product term in Eq. 1.5. The radar cross section is defined by the target area that is numerically equal to the ratio of the energy scattered in the direction of the receiving antenna to the total energy incident on the target (Skolnik, 1970). The bistatic
dependent on the wavelength and standard deviation of the surface roughness (Skolnik, 1970; Ulaby et al., 1982). An analytical expression for surface scattering, however, has not been derived. An example of a possible rough surface is the capillary fringe above the water table.
For calculating attenuation, however, volume scattering will be the dominant scattering mechanism of interest. The scattering caused by the theoretical vapor plume effects would be caused by scattering within the layers and not by the interfaces.
1.3.6.2 Volume Scattering
Volume scattering is also dependent on the wavelength of the propagating energy and how it interacts with the surrounding media. The amount of scattering is a function of the number, density, and scale of inhomogeneities (such as changes in RDP) and grain size distribution, relative to the dominant incoming frequency (Ulaby et al., 1982). The general rule is that this type of scattering becomes important as the size of the
heterogeneities or grain size approach 1/3 the length of the wavelength in the material. As the energy comes in contact with these scatterers, it is reflected in all directions, leaving increasingly less energy to be reflected toward the receiver at the surface.
As described in the previous section, volume scattering is difficult to properly quantify. Again, the equations presented in this section were originally developed for use in airborne radar techniques that must take into account particulate scattering from gases, clouds, and precipitation. Focusing on the theory of scattering by individual spherical particles, the backscattering cross section of a single spherical particle, σs, in m2, is
(
)
(
)
2 5 6 4 1 64 1 s R ε π σ λ ε − = + (1.26)ε is the absolute dielectric permittivity of the scatterer [F/m]. For Eq. 1.26 to be valid, the criterion 2 Rπ λmust be met (Skolnik, 1970).
Summing over the number of scatterers, N, per unit volume, V , extends this single particle scattering equation to volume scattering, as shown by
(
)
(
)
2 5 6 4 1 1 1 1 1 N 1 N s i i i R V V ε π η σ λ ε = = − = = +∑
∑
(1.27)where η is the radar reflectivity [1/m].
A particle size distribution can be substituted for a single particle radius for a heterogeneous soil. The λ−4 (equivalent to
4 2 ω πν ⎛ ⎜ ⎝ ⎠ ⎞
⎟ ) dependence of scattering is the Rayleigh scattering model. This model is used when the scattering particles have smaller diameters than the length of the incident wavelength. The Mie scattering model is used when the particle diameter is equal to or larger than the wavelength. For particles with diameters much larger than the wavelength, the scattering is not frequency dependent. There is, however, a resonant region when the particle diameter and wavelength are about equal that exhibits limited frequency and higher order ω dependence (Skolnik, 1970). Figure 1.3 shows the relationship between the particle size, wavelength, and normalized scattering cross section for both Rayleigh and Mie scattering regimes. Both of these models are for independent spherical particles. These models will most likely break down at higher densities of scatterers and when the particles are in contact with one another. They also assume that no multiple scattering occurs, which of course,
Figure 1.3 Radar cross section, σs, for a single spherical particle with radius, R, and incident wavelength, λ. The horizontal axis is the ratio of the circumference of the scatterer to the wavelength, and the vertical axis is the ratio of the radar cross section to the area of the scatterer. Increasing normalized radar cross section values along the vertical axis correspond to increasing permittivities of the scatterers. The Rayleigh region has a normalized cross section proportional to ω4 and is for particles with diameters smaller than the wavelength. The Mie region has a resonant, higher order dependence on ω and is for particles with diameters greater than the wavelength
(Skolnik, 1970). For GPR frequencies and typical scatterers in the ground, the region of interest is the Rayleigh region.
1.3.7 Dominant Mechanism Determination
Conductive, dielectric and magnetic relaxation, and scattering losses, unlike geometrical divergence, cannot be recovered in the data. The total attenuation, however, can be estimated after taking into account other, non-attenuation factors that affect amplitude in the data. These factors include the radiation patterns of the antennas, changes in antenna-ground coupling as the antennas are pulled along the surface,
geometrical divergence, and the transmitted wavelet into the ground. Once these factors are accounted for in the data, then any amplitude loss observed in the data should only be a result of the true attenuation mechanisms.
The next step is to determine which of the mechanisms is dominant. Based on the current research and theories of what is causing the shadow zones in GPR data acquired over an oil spill, the focus of the research is on scattering (due to vapor plume effects) and conductive (due to biodegradation or salt leaching or flushing) losses and their differing frequency dependences. Theoretically, when plotted as a function of loss tangent versus frequency, conductive and scattering losses can be distinguished by their relative slopes. As can be seen from Eq. 1.19, the conductive loss tangent is inversely proportional to frequency. Attenuation due to scattering, however, is assumed to be proportional to the radar reflectivity, η, (Eq. 1.27) which is directly proportional to frequency to some power. The loss tangent approximation therefore is directly
proportional to the radar reflectivity. The exponent is dependent on the model used to describe the volume scattering. For this study, because the loss tangent physically describes the phase difference between a lossy and a non-lossy medium, the following expression was used to estimate scattering losses:
tanδsc η β
= (1.28)
Figure 1.4 shows a plot of the conductive loss tangents versus frequency for εr = 9 and 25 and σ = 1 and 30 mS/m and the estimated scattering loss tangents for εr = 9 and 25 and number of silt particles per unit volume, N = 10 and 106.
As mentioned previously in Section 1.2, a material is considered low loss if the conductivity is below 30 mS/m. For typical GPR frequencies of 80-1500 MHz, this conductivity value roughly corresponds to the boundary between electromagnetic wave propagation and induction regimes. When the conductive loss tangent is plotted on a log scale, as shown in Figure 1.4, this boundary corresponds to a value of zero (tanδ =1). Negative log values correspond to wave propagation (necessary for the GPR method), and positive values correspond to the inductive regime, where diffusion occurs.
Therefore, in order to determine which of these attenuation mechanisms is dominating in this Bemidji GPR dataset, the trend of the frequency dependence of the loss will be plotted. If the loss relationship is directly proportional to frequency with some slope greater than one, then scattering losses are dominating. If the losses are greater with lower frequencies, however, then conductive losses are dominating. Figure 1.5 summarizes the general steps involved in this research.
Figure 1.4 A plot of conductive and volume scattering losses in terms of loss tangents as a function of frequency. The conductive loss tangent is calculated using Eq. 1.19, and the scattering loss tangent is an approximation using a ratio of the sum of backscattering cross sections of N discrete scatterers per unit volume to the phase constant, β (Eq. 1.28). The conductive loss tangent, plotted in green and blue, is calculated for all four
combinations of εr = 9 and 25 and σ = 1 and 30 mS/m, and the scattering loss tangent, plotted in black and red, is calculated for all combinations of εr = 9 and 25 and N = 10 and 106 for silt-sized grains.
Velocity & Attenuation Analyses Flow
Chart
Figure 1.5
Flow chart of the general steps necessary to perform
the velocity and attenuation analyses.
B
inned dat
a
R
em
oved range gai
n Flatten ed first arriv al tim es o n sel ect Sout h Pool Li ne offset s R e-appl ie d range gai n
Re-zeroed each offset
C onvert ed t o SU form at Im po rted in to VISTA Sem bl ance Pl ot & C onst ant
Velocity Stack Analyses
Secti on 4.1 Co ns tan t Velo city Pan el An alysis Secti on 4.2 St ol t M igrat ion C om pl ex
Trace Attribute Analysis
Section 4.3 C onvert ed t o ASC II for M A TLAB VELOCITY ANALYSIS R adi at ion pat te rn fact ors appl ie d Ant enna-ground coupl ing norm al izat io n Geom
etrical spreading factors applied
Ran sumaxspfx on ful l and tim e-windowed traces
Calculated attenuation as a funct
ion of frequency Frequency Spect ra Attrib ut e An alysis Secti on 5.1 Dom inant At te nuat ion Mech an ism Determ in atio n Secti on 5.2
Calculated loss tangent as a
funct ion of frequency Plo tted resu lts in MATLAB ATTENUATION ANALYSIS Determ ined frequency spectrum o f tran sm itted wav elet C onvert ed t o SU form at Sm oo th ed stack in g v elo cities Fi gure 4.1 Raw d ata Fi gure 3.6 Fi gure 5.1 Figure 5.20
CHAPTER 2
BEMIDJI RESEARCH SITE DESCRIPTION
2.1 Site History
The Bemidji research site is a crude oil spill site located approximately 16 km northwest of Bemidji, Minnesota. The area in which the spill occurred is a remote, forested area used mainly for recreational purposes (Hult, 1984). On August 20, 1979, a subsurface crude oil pipeline burst and spilled 1,700,000 L (about 10,700 barrels) of crude oil, contaminating both the surface and near subsurface. The pipeline that broke is 34 in. in diameter and is the largest of three parallel, closely spaced pipes (Herkelrath, 1999). At the rupture point, the pipe was 1.5 to 3 m below land surface (bls) and about 5 m above the water table. Directly below the break, some of the oil percolated downward through the vadose, or unsaturated, zone and accumulated on top of the water table to form the North Oil Pool. This particular pool is the largest accumulation of oil at the site with a 1.5 m thick layer of oil on the water table (Bennett et al., 1993). Based on
calculations using data taken from sediment cores obtained since 1989, it is estimated that the current total amount of free-phase oil present at this pool is 147,000 L (Herkelrath, 1999).
There is also a surface spray zone covering more than 7,500 m2 to the southwest of the pipeline break. The oil that sprayed onto the surface traveled eastward and then southward along a topographic low on the surface and formed an oil slick on a small, seasonal pond. This migrated oil eventually percolated down to the water table to form a second oil pool, called the South Oil Pool (Delin et al., 1998). It is estimated that 94,000 L of oil is present in this pool. In 1998, the Lakehead Pipe Line Company discovered that a third, smaller oil pool, since named the Middle Oil Pool, had formed in between the
North and South oil pools stretching westward from the oil pipeline (Herkelrath, 1999). Figure 2.1 is a 1991 aerial photograph of the Bemidji site with the spill-related features overlaid on top.
Cleanup efforts at the site began immediately after the spill and continued only for one month (Herkelrath, 1999). This remediation included excavating a large area around the pipeline break to repair it, as well as to recover some of the oil. The free-phase oil was pumped from the surface and from trenches that had been dug. When completed, it was estimated that 400,000 L (about 2,500 barrels) of oil remained in the subsurface (Delin et al., 1998). In 1982, the Office of Hazardous Waste Hydrology of the U.S. Geological Survey (USGS) selected the Bemidji site as a long-term,
interdisciplinary study site. Research began in May 1983. The goal of the research at this site is to better understand the mobilization, fate, and transport of crude oil in both the unsaturated and saturated zones by studying the physical, chemical, and biological processes that govern contamination in the subsurface (Delin, 2003a).
2.2 Geologic and Hydrogeologic Setting
The research site is located on the Bagley glacial outwash plain. This glacial outwash, consisting mostly of sand and gravel, is about 8 m thick and is underlain by other stratified glacial deposits that range in grain size from clay to gravel and extend to a depth of 20-25 m bls. Near the water table, thin, discontinuous silt lenses are interbedded with sand. The majority of the aquifer is quartz with 30% plagioclase and potassium feldspar, 5% dolomite, calcite, and low-Mg calcite, and about 5% heavy minerals including hornblende and ilmenite. Clay minerals, primarily kaolinite, smectite, and chlorite, total less than 1%. The carbonate percentage is higher in the silt fraction as compared to the bulk aquifer matrix (Bennett et al., 1993). Below the glacial deposits is a gray, calcareous, clay-rich, areally extensive till (Hult, 1984).
Figure 2.1 1991 aerial photograph of the Bemidji research site with the spill-related features overlaid on top. The three bold polylines indicate the approximate extents of the oil pools as of August 1998. The gray dashed lines leading to the east from the spray zone indicate the migration pathway of the oil on the surface that eventually formed the South Oil Pool. The square markers show the position of some of the wells at the site (Delin, 2003b).
The water table ranges in depth from 0 to 11 m bls (Smith and Hult, 1993). Near the South Oil Pool, the water table is about 2 m bls (Herkelrath, 1999). The linear flow velocity of groundwater in this area averages 10-50 cm/day toward the northeast (Bennett et al., 1993) and discharges into a small lake 400 m downgradient (Smith and Hult, 1993). Based on water table level measurements, the average hydraulic gradient is about 0.003. The till at the Bemidji site is similar to till elsewhere in Minnesota that has a vertical hydraulic conductivity of 9 x 10-3 cm/day (Hult, 1984), so the till at the site should act as an aquitard. The average annual precipitation from 1948 to 2000 was 613.7 mm (Delin, 2003c).
2.3 Well and Core Data
Figure 2.2 shows the position of the numerous wells that have been installed at the Bemidji research site. The wells provide continuous monitoring of the site, including changes in the water and oil levels over the oil pools. Driller’s logs were also included for some of the wells and provide useful information regarding the near surface geology and the approximate depths and lateral extents of contamination in both the saturated and vadose zones. Some geophysical logs have been performed including hole-to-hole radar.
The North Oil Pool has historically been the most well-documented and well- instrumented of the two main oil pools, as can be seen in Figure 2.2. For this reason, the GPR data in this thesis were acquired over and near the South Oil Pool to have minimum interference from the metallic wells.
8
Figure 2.2 Map showing the position of most of the wells at the Bemidji site. Note that the South Oil Pool is much less instrumented than the North Oil Pool. The W602 well cluster is located east-southeast of the South Oil Pool past the extents of the map (Delin, 2003d).
2.3.1 Oil Biodegradation
Based on samples acquired at two wells in the North Oil Pool, the oil is classified as a light paraffinic crude. Its composition is 58 to 61% saturated hydrocarbons,
(dominated by normal alkanes (C6-32)), 33 to 36% aromatics (benzene, naphthalene, and phenanthrene), 4 to 6% resins (N-, S-, and O-containing species), and 1 to 2%
asphaltenes (Eganhouse et al., 1993).
The contaminant plumes that extend both up and downgradient from the oil pools beyond the extents of the oil bodies are a result of geochemical reactions due to
biodegradation, volatilization, oxidation, and leaching of the oil. The biodegradation results in ground water with a lowered pH and increased levels of dissolved inorganic solutes. Based on the water data acquired, a model for the reaction sequence has been proposed. In the native ground water, the dominant solute is from carbonate dissolution with minor constituents including Na, K, and SiO2 derived from silicate minerals. Under aerobic conditions, an inorganic solute plume consisting of significantly increased levels Ca, Mg, and is formed by the complete biological oxidation of oil. Once this inorganic solute plume reaches the anoxic zone below the oil, reducing conditions exist, and Fe and Mn concentrations increase because of the dissolution of mafic silicate minerals. SiO
-3 HCO
2 concentrations also increase more than ten-fold. As the plume continues to flow downgradient and back into oxic conditions, Fe and SiO2 precipitate. The levels of the other inorganic solutes are then controlled by the amount of mixing with the native ground water (Bennett et al., 1993).
2.3.2 Contaminant Plume Migration
The North Oil Pool’s plume has been extensively studied, and the processes occurring there are assumed to also occur at the South Oil Pool. Based on water samples acquired in 1987, the plume is divided into five zones, each of which is defined by a particular aqueous chemistry. For a more complete description of each zone, refer to Bennett et al. (1993), Eganhouse et al. (1993), and Baedecker et al. (1993b). The North Oil Pool plume, as defined by the five zones, extended more than 100 m downgradient from the center of the oil body and about 5 m below the water table as of 1992. This migration, however, is slower than expected, considering models that used measured ground water flow velocities and conservative transport, and had actually slowed slightly between 1992 and 1995. It is thought that the lower hydraulic conductivities of
discontinuous silt layers impede the plume advancement and that natural biodegradation of the oil is a considerable factor (Baedecker et al., 1993a; Cozzarelli et al., 1999). A plume defined by elevated total dissolved inorganic solutes (TDS) extended more than 150 m downgradient from the center of the north oil body in 1987 (Figure 2.3b). Figure 2.3a is a plot of the specific conductances at 25ºC, measured from water samples
collected in June 1987 from several wells along the North Oil Pool axis. The background TDS and specific conductance levels of the native ground water are typically less than 400 mg/l and between 25 and 35 mS/m, respectively (Bennett et al., 1993). The elevated TDS and specific conductance values are supportive evidence for the conductive losses that were proposed by Sauck et al. (1998) as the cause of shadow zones observed in GPR data acquired over LNAPL spills.
Although the South Oil Pool is smaller, a similar contaminant plume must logically extend from that oil body as well. The size and shape of this plume will be dependent on local variations of the geology (i.e. presence of silt lenses) and distribution of the oil in the subsurface. Figure 2.4 shows the locations of two cross section
(a)
(b)
Figure 2.3 (a) Specific conductances, in mS/m at 25ºC, of water samples collected in June 1987 within the North Oil Pool and (b) total dissolved solids (TDS) levels, in mg/l, of the North Oil Pool in 1987. The vertical lines extending downward from the land surface indicate the sampling intervals for the data presented in (a) and (b) (Lucius, 2000; Bennett et al., 1993).
C’ E’
D’
E D
C
Figure 2.4 Topographic map of the Bemidji site showing the positions of the North Pool (A-A’) and South Pool (B-B’) transects shown in Figure 2.5 (Dillard et al., 1997). The gray, dotted square indicates the approximate extent of the aerial photograph shown in Figure 2.1, and the gray, dashed square indicates the approximate extent of the map in Figure 2.2. The three dotted-dashed lines labeled C-C’, D-D’, and E-E’ indicate the positions of three 80 MHz GPR profiles acquired in 1985 that are shown in Figure 2.6.
