|USFDC Home | USF Electronic Theses and Dissertations||| RSS|
This item is only available as the following downloads:
xml version 1.0 encoding UTF-8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd
leader nam 22 Ka 4500
controlfield tag 007 cr-bnu---uuuuu
008 s2010 flu s 000 0 eng d
datafield ind1 8 ind2 024
subfield code a E14-SFE0003489
Optimizing the imaging of multiple frequency gpr datasets using composite radargrams :
b an example from santa rosa island, florida
h [electronic resource] /
by Stuart Bancroft.
[Tampa, Fla] :
University of South Florida,
Title from PDF of title page.
Document formatted into pages; contains X pages.
Thesis (M.S.)--University of South Florida, 2010.
Includes bibliographical references.
Text (Electronic thesis) in PDF format.
Mode of access: World Wide Web.
System requirements: World Wide Web browser and PDF reader.
ABSTRACT: Acquiring GPR data at multiple frequencies is useful because higher-frequency profiles have better spatial resolution, although they suffer from reduced depth penetration. Lower-frequencies can generally resolve to greater depths, but at the cost of spatial resolution. For concise presentation of GPR data, it would be useful to combine the best features of each profile into a composite radargram. This study explores effective ways to present GPR data acquired at multiple frequencies. An example is shown from a survey of hurricane overwash deposits from Santa Rosa Island, Florida. The methodology used to create a composite radargram is dependent on which of two goals the composite radargram is designed to achieve. These goals are broadening the spectral bandwidth of GPR data to increase the effectiveness of deconvolution and enhancing the resolution and depth of GPR data by plotting high-frequency data at early two-way travel times, low-frequency data at late two-way travel times, and using filters to smoothly transition from high-frequency to lower-frequency data. The steps towards creating a composite radargram include: 1) applying standard processing to nominal frequency data sets, 2) creating spatially coincident data sets, 3) equalizing the amplitude spectra among each nominal frequency data set, and 4) summing nominal frequency data sets together. Spectral bandwidth broadening is achieved by applying optical spectral whitening and summing nominal frequency data sets using a single ramped. Deconvolving this composite radargram did not show the same success observed by Booth et al. (2009). Enhancing the resolution and depth of GPR data can be achieved by applying amplitude envelope equalization (AEE) and summation using double ramped filters. AEE calculates the coefficients required to make equivalent average amplitude envelopes for GPR data that has been gained with automatic gain control Double ramped filters suppress low-frequency energy for two-way travel times when a higher-frequency data set has adequate signal strength and higher frequency energy for two-way travel times when higher- frequency energy exhibits significant attenuation. A composite radargram built with AEE and double ramped filters achieves the goal enhancing resolution and depth of GPR data. Shallow reflections are interpreted as dune and hurricane overwash stratigraphy.
Advisor: Sarah Kruse, Ph D.
Ground penetrating radar
Multiple frequency antennae
t USF Electronic Theses and Dissertations.
i Optimizing the Imaging of Multiple Frequency GPR Datasets Using Composite Radargrams: An Example From Santa Rosa Island, Florida by Stuart W. Bancroft A thesis submitted in partial fulfillment of the requirements for the degree of Master o f Science Department of Geology College of Arts and Sciences University of South Florida Major Professor: Sarah Kruse, Ph.D. Diana Roman, Ph.D. Mark Stewart, Ph.D. Date of Approval: April 2, 2010 Keywords: ground penetrating radar, multiple freque ncy antennae, hurricane overwash, coastal stratigraphy, barrier island Copyright 2010, Stuart W. Bancroft
ii Acknowledgements I thank Dr. Sarah Kruse for embracing and guiding my creativity towards developing this research, as well as assisting me with computer programming. Invaluable field assistance was provided by Dr. Ping Wang, Rip Kirby, and Mark Horwitz of the USF Coastal Research Laboratory. Rip handled the security clearance required to perform the GPR surveys on Eglin Air Forc e Base. Mark assisted with data acquisition and geologic interpretations. Dr. Wang assisted with data acquisition and guided the trenching and coring process used for this research. I would also like to thank my wife, family, and colleagues for their posit ive encouragement.
iii Table of Contents List of Tables ii List of Figures iii Abstract ix Introduction 1 Chapter One: Ground Penetrating Radar 3 Wave Propagation 3 Frequency Dependent Properties 7 GPR Acquisition 10 GPR in Coastal Environments 11 Chapter Two: Field Example 14 Geologic Setting 14 Acquisition 17 Basic Data Processing 24 Chapter Three: Composite Radargrams 26 Introduction 26 Spatial Coincidence 29 Spectral Bandwidth Broadening 37 Summation Techniques 37 Amplitude Balancing 44 Deconvolution 52 Resolution and Penetration Enhancement 55 Double Ramped Summation 55 Amplitude Envelope Equalization 68 Subjectively Chosen Values 72 Chapter Four: Geologic Interpretation 75 Velocity Analysis 75 Field Verification 77 3D Data Composites 83 Chapter Five: Conclusions 92 References 95 i
iv List of Tables Table 1. Electromagnetic properties of geologic materials common to coastal environments. 12 Table 2. Summary of the basic processing steps applied to GPR data in ReflexW. 25 25 Table 3. A subset of a 100 MHz time slice before jitter correction. 32 Table 4. A subset of a 100 MHz time slice after jitter cor rection. 33 Table 5. This table shows the inaccuracy of the odomet er wheel calibration setting for each survey. 36 Table 6. The period for each antenna center freque ncy is calculated to determine the ramp length chosen to be 15 wave periods. 6 2 Table 7. Each average envelope value is divided by the h ighest average envelope value, resulting in a multiplier value for each nominal frequency data set. 71 Table 8. Chosen values compared to calculated values for composite radargram with double summation and AEE balancing. 7 2 ii
v List of Figures Figure 1. Components of GPR system. 3 Figure 2. Remote station where acquired data can be viewed in real time in a stationary position. 4 Figure 3. The relationship between frequency and vel ocity for materials of varying conductivities. Adapted from Davis and Annan (1989). 7 Figure 4. The relationship between frequency and attenuation for material s of varying conductivities. 8 Figure 5. Flow diagram that accounts for the various ca uses of signal loss within a GPR system. Adapted from Annan and Davis (1977). 9 Figure 6. Map of the state of Florida with inset of Santa Ros a Island. 15 Figure 7. Google Earth image of inset for Figure 6 with symbol over study site. 17 Figure 8. Google Earth image of study site, survey grid drawn in yellow. 17 Figure 9. Survey site looking north (landward). 18 Figure 10. Acquisition with 250 MHz antenna. 20 Figure 11. Diagram of survey site showing location of surveys and field verification sites. 23 Figure 12. Using unscaled amplitudes results in a com pos ite dominated by low frequency energy. 27 Figure 13. Applying time shifts to scaled amplitude synth etic wavelets results in a composite wavelet that more closely resembles a delta function, the ideal wavelet for deconvolution 28 Figure 14. Time slice of 500 MHz data at 14.7 ns, before ji tter corrections were applied. 34 Figure 15. Time slice of 500 MHz data at 14.7 ns, after jitte r corrections were applied. 35 iii
vi Figure 16. Composi te grid (see Figure 11) showing the lo cation of 2D profile used for summation. 38 Figure 17. Composite radargram create by brute summation. 39 Figure 18. 100 MHz radargram. 40 Figure 19. 800 MHz radargram. 40 Figure 20. Composite radargram created by coarse muting. 41 Figure 21. From left to right, 500 MHz, 250 MHz, and 100 MHz radargrams. 42 Figure 22. Composite radargram created by ramped summation. 43 Figure 23. DFAE for synthetic wavelets. 45 Figure 24. OSW least squares algorithm in matrix notation 46 Figure 25. From top to bottom, single trace of 100 MHz da ta, S transform of trace, normalized S transform. 47 Fi gure 26. From top to bottom, single trace of 250 MH z data, S transfo rm of trace, normalized S transform. 48 Figure 27. From top to bottom, single trace of 500 MH z data, S transform of trace, normalized S transform. 49 Figure 28. From top to bottom, single trace of 800 MHz data, S transform of trace, normalized S transform. 5 0 Figure 29. From top to bottom, single trace o f composite data (single ramped summation with OSW balancing), S transform, normal ized S transform. 51 Figure 30. Nominal frequency radargrams after OSW mu ltipliers have been applied. 52 Figure 31. The sum of radargrams from Figure 30. 53 Figure 32. Deconvolved radargram from Figure 31. 54 Figure 33. Double summation ramps using a two pole Butterworth function. 56 Figure 34. Composite radargram produced with double ramped summation 57 iv
vii Figure 35. Double summation ramps using a six pole Butterworth function. 58 Figure 36. Composite radargram produced with double ramped summation. 59 Figure 37. Double summation ramps using linear functions. 60 Figure 38. Composite radargram produced with double ramped summation. 61 Fi gure 39. A comparison between composite radarg rams produced using double linear ramps with lengths of 2 and 15 periods of outgoing frequency. 62 Figure 40. The log amplitude envelope of ungained data plotted against time. 64 Figure 4 1. 250 MHz radargram with black solid line at time when suppression begins. 65 Figure 42. 500 MHz radargram (AGC gained) with black so lid line at time when suppression begins 66 Figure 43. 800 MHz radargram (AGC gained) with bl ac k solid line at time when suppression begins 67 Figure 44. The amplitude envelope for AGC gained da ta and average envelope value plotted against the time period when the nominal frequency data is not suppressed within the composite. 7 0 Figure 45. Composite radargram built with double r amped summation and AEE amplitude balancing. 71 Figure 46. Composite radargram built with calculated values from Table 8. 73 Figure 47. Composite radargram built with ch osen values from Table 8. 74 Figure 48. CMP gather with a best fitting parabola corre sponding to a velocity of 0.10 m/ns. 76 Figure 49. 800 MHz radargram with a fitted di ffraction hype rbola corresponding to an EM wa ve velocity of 0.12 m/ns. 7 7 Figure 50. Core SR 1. Location shown in Figure 11. 78 Figure 51. Trench 2 looking west. 79 iv v
viii Figure 52. Trench 1 looking north (landward). 80 Figure 53. Trench 1 looking east. 81 Figure 54. Composite radargram with structures with interpre ted radar facies outlined in red. 82 Figure 55. Composite radargram of profile 1 meter west of eastern boundary. 83 Figure 56. Compos ite radargram of profile 2 meters west of eastern boundary. 84 Figure 57. Composite radargram of profile 3 meters west of eastern boundary. 84 Figure 58. Composite radargram of profile 4 meters west of eastern boundary. 85 Figur e 59. Composite radargram of profile 5 meters west of eastern boundary. 86 Figure 60. Depth slice at 0.2 meters. 87 Figure 61. Depth slice at 1.5 meters. 88 Figure 62. Depth slice at 2 meters. 89 Figu re 63. Depth slice at 4 meters. 90 Figure 64. Depth slice at 8 meters. 91 vi
ix Optimizing the Imaging of Multiple Frequency GPR Datasets using Composit e Radargrams: An Example from Santa Rosa Island, Florida Stuart W. Bancroft ABSTRACT Acquiring GPR data at multiple frequencies is useful because higher frequency profiles have better spatial resolution, although they suffer from reduced depth penetratio n. Lower frequencies can generally resolve to greater depths, but at the cost of spatial resolution. For concise presentation of GPR data, it would be useful to combine the best features of each profile into a composite radargram. This study explores effec tive ways to present GPR data acquired at multiple frequencies. An example is shown from a survey of hurricane overwash deposits from Santa Rosa Island, Florida. The methodology used to create a composite radargram is dependent on which of two goals the c omposite radargram is designed to achieve. These goals are broadening the spectral bandwidth of GPR data to increase the effectiveness of deconvolution and enhancing the resolution and depth of GPR data by plotting high frequency data at early two way trav el times, low frequency data at late two way travel times, and using filters to smoothly transition from high frequency to lower frequency data. The steps towards creating a composite radargram include: 1) applying standard processing to nominal frequenc y data sets, 2) creating spatially coincident data sets, 3) equalizing the amplitude vii
x spectra among each nominal frequency data set, and 4) summing nominal frequency data sets together. Spectral bandwidth broadening is achieved by applying optical spectra l whitening and summing nominal frequency data sets using a single ramped. Deconvolving this composite radargram did not show the same success observed by Booth et al. (2009). Enhancing the resolution and depth of GPR data can be achieved by applying ampl itude envelope equalization (AEE) and summation using double ramped filters. AEE calculates the coefficients required to make equivalent average amplitude envelopes for GPR data that has been gained with automatic gain control Double ramped filters suppr ess low frequency energy for two way travel times when a higher frequency data set has adequate signal strength and higher frequency energy for two way travel times when higher frequency energy exhibits significant attenuation. A composite radargram built with AEE and double ramped filters achieves the goal enhancing resolution and depth of GPR data. Shallow reflections are interpreted as dune and hurricane overwash stratigraphy. viii
2 Introduction Ground penetrating radar (GPR) is a near surface geoph ysical tool that utilizes high frequency radio waves to detect subsurface changes in electromagnetic properties. Field verification can correlate these changes to interpreted lithologic and water content changes (among others) that may provide geologically useful information. The high frequency radio pulse emitted from commercial GPR antennae typically has a frequency bandwidth of one octave with a center frequency ranging from 100 1000 MHz. The reason for manufacturing antennae with different center freque ncies is that the spatial resolution and depth penetration are highly frequency dependent. Shallow structures may be imaged with high spatial resolution with a high frequency GPR (500 1000 MHz). Deeper structures can only be sensed by lower frequency GPR a ntennae (100 500 MHz), but with lower spatial resolution. Therefore any nominal frequency GPR data set provides an incomplete view of the subsurface in terms of either spatial resolution or penetration depth. This limitation of GPR can be reduced by incorp orating multiple ground penetrating radars with different center frequencies into a single survey. An increased understanding of the subsurface may be achieved by interpreting the multiple data sets jointly. To my knowledge, the joint interpretation of m ultiple frequency data sets has only been done three ways: 1) side by side analysis of multiple spatially coincident profiles of different frequencies, 2) simple summations, and 3) combining profiles with the goal of maximizing bandwidth to increase the ef fectiveness of wavelet deconvolution (Booth et al., 2009). I am not currently aware of studies that explore different methods of combining multiple frequency data sets so that shallow and deep structures may be 1
2 clearly displayed. In this thesis, such metho ds will be developed, evaluated, and compared to the deconvolution based method using a GPR survey of Santa Rosa Island, Florida.
3 Chapter 1: Ground Penetrating Radar Wave Propagation Ground penetrating radar (GPR) is a non destructive geophysical tool that uses high frequency radio waves to detect changes in electromagnetic properties of the shallow subsurface. A GPR system contains a transmitting and a receiving antenna that may be shielded or unshielded, a computer that stores the radar data, and fi ber optic cables that transmit data from the radar to the computer, shown in Figures 1 and 2. An electromagnetic pulse is generated by the transmitting antenna and propagates as a wave through the subsurface. A portion of the electromagnetic energy from t he wave is reflected or scattered towards the receiving antenna when the propagating wave comes into contact with a point or interface that represents a contrast in electric permittivity (Baker et al. 2007). The receiver records the two way travel time an d amplitude of the reflected energy wave. odometer wheel shielded antennas fiber optic cables computer Figure 1. Components of GPR system.
4 Figure 2. Remote station where acquired data can be viewed in real time in a stationary position. Maxwells equations provide the foundation fo r the physical basis of EM wave propagation. They are as follows: x H = J + ( D/ t) (1), x E = B/ t (2), D = q (3), and B = 0 (4). In these equations, H is the magnetic field intensity, J is the electric current density vector, D is the ele ctric displacement vector, E is the electric field strength vector, B is the magnetic flux density vector, q is the electric charge density, and t is time. Amperes
5 law, Equation 1, states that a magnetic field can be induced from a time varying electric field or a current loop. Faradays law, Equation 2, states that an electric field can be induced from a time varying magnetic field. Gausss law, Equation 3, states that the electric field flux is equal to the charge responsible for the field. The last o f Maxwells equations, Equation 4, states that magnetic monopoles have never been observed. Manipulating Maxwells equations results in the derivation of the electromagnetic wave equation. This equation describes how electromagnetic waves propagate throu gh dielectric media. The electromagnetic wave equation can be expressed in terms of the magnetic or electric field as, 2 E= ( 2 E/ t 2 ) g ( E/ t) (5). Equation 5 is expressed in terms of the electric field, where E is the electric field strength vector, is electric permittivity, is magnetic permeability, g is electrical conductivity, and t is time. This equation can describe the propagation of the magnetic component of EM waves by substituting H the magnetic field intensity, for E The elect romagnetic wave equation reveals two fundamental characteristics of EM wave propagation. The first two terms, ( 2 E) and ( ( 2 E/ t 2 )), indicate that the second spatial derivative of E is proportional to the second temporal derivative of E The solutio ns to differential equations in this form are oscillatory functions, revealing that electromagnetic energy propagates as a wave. The third term, g ( E/ t) is the differential form of an exponential decay function. This indicates that the amplitude env elope of an electromagnetic pulse decays exponentially with time as a result of the conductivity of the subsurface medium.
6 It is important to understand the causes of EM wave reflections picked up by the receiving antenna. Waves propagate at a velocity de termined by the properties of the medium through which the wave is traveling. When a wave enters a new medium, the velocity of the wave changes, and its energy is partitioned at the boundary. A portion of the waves energy is reflected at an angle equal t o the angle of incidence. Another portion is transmitted into the new medium as refracted energy. For EM waves emitted by GPR, the main physical property that controls wave velocity is electric permittivity. The contrast in electric permittivity between t wo materials determines the amount of energy reflected towards the receiving antenna. Electric permittivity is a measure of how well a material can store electromagnetic energy and polarize when subjected to an external electromagnetic field (Baker et al., 2007). Materials with a measurable electric permittivity are called dielectric materials. Since most natural materials are dielectric, EM waves can propagate through them for a significant number of wave cycles. However, this does not make GPR necessari ly suitable for all dielectric materials. Penetration depth is discussed further below. Workers commonly use a materials permittivity, relative to the permittivity of air, 0 to describe the materials dielectric capability. Relative permittivity is defined as r = / 0 Davis and Annan (1989) and Daniels et al. (1995) determine the relative permittivity for various geologic materials. Most common minerals have relative permittivities of 2 4. Water has a relative permittivity of 80. For this reason, wa ter content is the principal determinant of bulk permittivity of porous media.
7 Frequency Dependent Properties There has always been a need to provide subsurface geophysical evaluations using methods that provide deep penetration and high spatial resolutio n. GPR at least partially satisfies these requirements and additionally provides results rapidly and economically (Davis and Annan 1989). GPR antennae are available in several frequency bands with a center frequency ranging from 2 1200 MHz. The range of available antenna frequencies arises from the frequency dependent nature of radar wave propagation. Figure 3, adapted from Davis and Annan (1989), shows how wave velocity is dependent on antenna frequency and the conductivity of the wave medium. An import ant observation of this relationship is that the velocity is relatively constant for typical conductivities of most geologic materials (less than 100 mS/m) and typical center frequencies for GPR antennae (10 1000 MHz). Figure 3. The relationship betw een frequency and velocity for materials of varying conductivities Adapted from Davis and Annan (1989). The penetration depth of a radar signal is dependent on several factors including radar system performance, ground attenuation, and the contrast in d ielectric permeability at subsurface boundaries (Davis and Annan 1989). Of those properties, ground
8 attenuation is most dependent on antenna frequency. This figure below adapted from Davis and Annan (1989) shows the relationship between frequency and atten uation at different conductivities for a medium with a dielectric constant of 4. Clearly, attenuation increases significantly with frequency and conductivity. Figure 4. The relationship between frequency and attenuation for materials of varying conducti vities. Adapted from Davis and Annan (1989). Approximating the actual radar penetration depth involves the other non frequency dependent factors and can be approximated using Q = T R G T G R g c e 4 L ___________________________ __ 64 3 f 2 L 4 A description of these variables is found in Figure 5, adapted from Annan and Davis (1977). Q is a quality factor or ratio of the transmitter signal amplitude to minimum receiver sensitivity. The role of frequency in this equation can be traced to the amount of power received by the antenna. The receiver antenna effective area is inversely proportional to frequency. (6).
9 Figure 5. Flow diagram that accounts for the various causes of signal loss within a GPR system. Adapted from Annan a nd Davis (1977). The amount of information that can be extracted from GPR data is dependent on spatial resolution, the ability to distinguish separate reflections or pulses. When two adjacent pulses occur close together in time or space, the two pulses ca n only be distinguished as separate pulses if their wavelength is narrow. The minimum thickness of a reflective layer that can be detected by radar is proportional to the wavelength of the incident wave. As wavelength and frequency are related to each oth er as reciprocals, it is inherent that resolution is related to antenna frequency. The range of available radar antenna frequencies creates a dilemma termed by Davis and Annan (1989) as range resolution trade off. This compromise describes the ability o f higher frequency antennae to have fine spatial resolution at the cost of less depth penetration, and lower frequency antennae provide much deeper penetration at the cost of compromising spatial resolution. In practice, spatial resolution is compromised b efore depth penetration because it is better
10 to image a feature with poor resolution than to not see it at all due to rapid signal attenuation. GPR Acquisition GPR acquisition is in many ways analogous to seismic reflection methods. Different modes of ac quisition are determined by the arrangement of the transmitting and receiving antennae. GPR acquisition is usually done in common offset mode, i.e., the separation between the two antennae remains constant. Common offset profiling reveals the two way trav el time of reflective surfaces versus the distance from a fixed point. GPR reflections are more geologically useful when shown in terms of depth instead of two way travel time. This conversion requires knowledge of the subsurface velocity structure. Anot her type of acquisition mode, common midpoint (CMP) soundings, can be implemented to determine radar velocity. CMP sounding works by increasing the spacing between the transmitting and receiving antennae such that the midpoint between the two antennae rema ins constant. The radargram that results from CMP soundings shows reflections from the air wave, direct wave, and subsurface reflections. The air and direct waves appear as lines whose slope corresponds to the radar velocities of those materials. The subsu rface reflections appear as hyperbolas whose shape corresponds to the root mean square velocity of the subsurface from the depth of the reflection to the surface. Unlike seismic velocities, radar velocities are especially time sensitive due to their depen dency on moisture content. The velocities obtained from a CMP survey can only be applied to a common offset survey if the moisture content has not changed during the elapsed time between surveys. Velocities can also be determined with common offset acquis ition in two ways: 1) fitting hyperbolas to diffractions and 2) correlating the two
11 way travel time of a reflection to the distance of the reflective surface determined through field verification. GPR in Coastal Environments In recent years, GPR has had an immense effect on coastal sedimentology. GPR is commonly favored over other geophysical methods when combined with coring and auger methods for subsurface stratigraphic studies. GPR acquisition provides continuous, high resolution spatial coverage, time and cost effective surveys, and is non invasive. GPR is a proven method in determining stratigraphic architecture, sand body geometry, and correlation and quantification of sedimentary structures (e.g. Bristow and Jol, 2003). Aeolian coastal environments fluvial and alluvial fans, and glacial environments are among the sedimentary depositional environments that can be studied using GPR. Leatherman (1987) was one of the first to evaluate the potential applicability of GPR in coastal settings. Prior to his study, coring and augering were the common methods used to study coastal geomorphology. These methods are limited by a penetration depth of less than two meters and non continuous spatial coverage. Seismic reflection and refraction have shown little promise in coastal settings (Leatherman, 1987). GPR is generally suitable for coastal environments comprised mainly of electrically resistive quartz sand and a thick freshwater lens. Areas with saltwater or significant amounts of silt or clay create an el ectrically conductive setting are less suitable for GPR surveys. Davis and Annan (1989) determined the relative permittivity and EM velocity for various geologic materials including those commonly found in coastal settings, which are summarized in Table 1
12 Table 1. Electromagnetic properties of geologic materials common to coastal environments. Modified from Davis and Annan (1989). Material Relative Permittivity EM velocity (m/ns) Conductivity (mS/m) Freshwater 80 0.03 0.5 Saltwater 80 0.0 3 3000 Sand, dry 3 5 0.12 0.17 0.01 Sand, wet 20 30 0.05 0.09 0.1 1 Although freshwater and saltwater have equivalent relative permittivities and EM velocities, the high conductivity of saltwater causes rapid signal attenuation. Therefore the saltwater interface is easily identifiable in GPR images not as a reflection but as a sudden reduction in signal amplitude. After the pioneering study done by Leatherman (1987) that advocated the use of GPR in coastal stratigraphic studies, Jol (1996) evaluated th e potential use of GPR studies in various coastal setting in Florida, Georgia, Texas, Oregon, and Washington. Jol (1996) found that GPR penetration depths were greatest in clastic environments, but even very shelly sediment and coarse gravel allowed for re asonable penetration. Silt and brackish water were the two main materials responsible for less than excellent penetration: for example, a penetration depth of only 2.75 m was achieved with a 100 MHz antenna due to brackish water on Galveston Island (Jol, 1996). Correlating core or trench data to GPR data allows for determination of the continuous spatial distribution of
13 the water table and saltwater interface. The character of GPR reflections (dip angle and continuity) and ground truthing can be used to determine the geometry of sedimentary structures such as paleo beach face surfaces of prograding and retrograding shorelines, preserved sand bars, tidal inlet fill, and washover fans. Numerous studies of GPR and its application to sedimentology and strati graphy can be found in (Bristow and Jol, 2003) and (Baker and Jol, 2007).
14 Chapter 2: Field Example Geologic Setting A multiple frequency test data set was acquired on Santa Rosa Island, a 95 km long barrier island located along the Florida Panhandle from Destin to Pensacola Pass, shown in Figure 6 (Stone 2004). The site was chosen because it offers stratigraphic heterogeneity on multiple scales and penetration depths across the spectrum of GPR frequencies. The stratigraphy is of interest t o coastal geologists investigating erosional and hurricane depositional characteristics of overwash deposits. Santa Rosa Island is bounded to the north by Pensacola Bay and Choctawatchee Bay, two bodies of water joined by Santa Rosa Sound, and by the Gulf of Mexico to the south. Santa Rosa Island is a Holocene barrier island composed almost entirely of fine grained and well sorted quartz sand. Sediment transported by the Apalachicola River makes its way to Santa Rosa Island via westward littoral drift (Sto ne 1996). Santa Rosa Island has a subtropical climate with an average annual precipitation of 152 cm. In addition to the occurrence of tropical cyclones, the summer months are characterized by heavy precipitation from thunderstorms and light southerly win ds. Cold fronts bring precipitation and strong persistent northerly winds throughout fall, winter, and spring (Miller 2001). Typical wave heights for Santa Rosa Island are approximately 0.7 m. The average tidal range is 0.43 m, but can vary from 0.15 to 0.61 during equatorial and tropic phases (Stone 2004).
15 Figure 6. Map of the state of Florida with inset of Santa Rosa Island. Modified from Preis (2009). Historically, Santa Rosa Island, especially within the vicinity of Navarre Beach, has had a ver y stable shoreline due to balanced rates of offshore and longshore sediment transport. Within recent years (since 1995) Santa Rosa Island has been affected by several major hurricanes including Opal (1995) and Georges (1998), and more recently, Ivan (2004 ) and Dennis (2005). These storm events have had a significant effect on the morphology of the barrier island and have greatly disrupted the equilibrium of sediment transport for the island. The low elevation and relief of Santa Rosa Island make it espec ially susceptible to overwash during tropical cyclone events. Wang and Horwitz (2006) give detailed accounts of the sedimentological impact of hurricanes Ivan and Dennis on Santa Rosa Island. The nature of overwash deposits is dependent on whether or not the barrier island is breached during a storm. Overwash events associated with
16 breaching deposit sediment in the intertidal zone. The overwash deposits may form a flood tidal delta depending on how well the breach connects the ocean to the back barrier b ay. The overwash deposits found on Santa Rosa Island associated with Ivan and Dennis did not breach the island. For such cases, the relative values of the dune field elevation and the elevation of the landward extent of swash determine whether overwash o r inundation occurs. Sedimentary structures associated with the non breaching overwash events include severely truncated dunes overlain by an approximately horizontal overwash platform, and slightly landward dipping beds. Study Site The GPR surveys were conducted at Eglin Air Force Base at coordinates 30.39053 o N, 86.78873 o W (Figures 7 and 8). Previous studies (Horwitz et al., 2006) indicate that the survey site was likely near the transition zone from the dune field and back barrier marsh prior to Hurr icane Ivans landfall in 2004. Post landfall, the survey site was inundated by 3.5 meters of storm surge that deposited a significant amount of sediment over truncated dunes (Horwitz et al, 2008).
17 Figure 7. Google Earth image of inset for Figure 6 with symbol over study site. N 200 m Figure 8. Google Earth image of study site, survey grid drawn in black. Acquisition The GPR surveys were conducted on Eglin Air Force Base April 24 25, 2009. The skies were clear and sun ny; no precipitation occurred during the survey period. Therefore, there is no reason to suspect significant changes in the subsurface velocity profile within the survey period. An optimal survey site for building GPR composite N 10 km
18 radargrams should have geolo gically significant features located at various depths. A reconnaissance survey with the 250 MHz antenna was used to find such an area that clearly exhibited the expected sedimentary structures associated with hurricane overwash. The chosen area is a rec tangular plot 22 m long, north to south, and 10 m across, east to west. This site was also optimal for a GPR survey because there were no physical obstructions on the surface, no sudden changes in topography, and no anthropogenic features such as buried ut ilities or roads. Figure 9. Survey site looking north (landward). There is no significant topography within the survey site. The area surveyed with all antennae (composite grid) is shown in red.
19 For high resolution 3D GPR surveys, it is crucial to have a systematic method for acquiring several 2D lines of data spaced closely together. This survey required the effort of three people at a time. One person was responsible for pulling the antenna in a straight line and carryin g a backpack containing the radar control unit, laptop computer, and appropriate batteries. String was placed along the ground over 25 centimeter intervals within the grid. Another person was responsible for adjusting the position of the string during acqu isition. The third person monitored the quality and progress of data in real time via a remote desktop computer that was wirelessly connected to the laptop on the backpack, and kept detailed field notes on the survey. Since the line spacing for the 500 M Hz and 800 MHz surveys was 5 cm, it was not practical to adjust the guiding string after every 2D line was acquired. Instead, two plastic rods with different colored strips of electrical tape placed every 5 cm dangling down so that the strips and the guid ing string can be used together to guide each line. For all 3D surveys, the odd numbered 2D lines were acquired traveling north while the even numbered lines were acquired traveling south.
20 Guiding string placed parallel to this survey boundary colored strips of tape Figure 10. Acquisition with 250 MHz an tenna. The colored strips of electrical tape spaced 5 cm apart help guide the survey lines. For the first line, the red tape is flush with the survey tape or parallel string placed in multiples of 20 cm left of the labeled survey line; the next line uses t he orange tape as a guide.
21 Acquisition with the 250 MHz antenna was done on the afternoon of April 24 th The antenna was positioned in the southeast corner of the grid. The first of ninety seven lines was acquired with the antennas centered 37 cm wes t of the eastern boundary of the survey grid. The lines ran north to south for 22 m. The spacing between lines was 10 cm, and the trigger interval was 5 cm. The survey distance was measured with the odometer wheel that attaches directly to the radar anten na. Although the odometer wheel was calibrated to the surface conditions of the survey site, there were clearly inaccuracies associated with odometer distance measurements. The survey length measured with survey tape was 22 m. However, the typical survey length measured by the odometer wheel was 22.75 m. Data acquisition with the 100 MHz, 500 MHz, and 800 MHz antennae was completed on the day of April 25 th The antenna setup, odometer calibration, trigger interval, and path of acquisition were all the s ame as in the 250 MHz survey. The only different survey parameter was a line spacing of 20 cm for the 100 MHz survey, resulting in a total of 48 lines. The 500 MHz and 800 MHz antennae acquired data simultaneously with a line spacing of 5 cm. The small s ize of these antennae and the capability of control unit to connect two antennae at once make multiple antenna acquisition very feasible. At the starting and end positions for each survey line, the center of the 500 MHz antenna was flush with the northern and southern boundaries of the grid, and the center of the 800 MHz antenna was 41.3 cm behind. The 800 MHz antenna was secured so that it could move up and down, in order to maximize direct contact with the ground, but not left or right so that the simult aneously acquired data sets would not have spatial coincidence issues. Through review of the results of the 250 MHz survey, we found that
22 it was not necessary to survey the entire 22 by 10 meter grid because the geophysical targets of interest could be ad equately imaged within a subset of the original survey grid. This sub grid is a 10 m, north to south, by 5 m, east to west, rectangle that is the northeast corner of the original grid, as shown in Figure 9. The survey had the following parameters: 5 cm li ne spacing, 2.5 cm trigger interval, and a survey length measured as 10 meters by survey tape and 10.32 meters by the odometer wheel. Although diffraction hyperbolae were present within the data, a common midpoint survey was done to determine the subsurf ace velocity profile of the survey area. The CMP survey was done on the southern boundary of the original grid with the common midpoint being 5.4 meters west of the eastern grid boundary. The first measurement was taken with 100 MHz and 250 MHz antennae touching each other so that the source to receiver distance was 95 cm. The next nine measurements were taken with the source to receiver distance increased by multiples of 10 cm.
23 10 m 5 m Figure 11. Diagram of survey site showing location of surveys and field verification sites. Survey grid for 100 MHz and 250 MHz GPR 10 x 22 meters Survey grid (in red) for 800 MHz, 500 MHz, 250 MHz, and 100 MHz GPR (composite grid) 5 x 10 meters Location of Trench 1 Location of Trench 2 Location of pound core (SR 1) Location of CMP survey
24 Basic Data Processing The process of building GPR composite radargrams consists of two phases of processing: 1) standard processing of single frequency data sets independently using the so ftware package ReflexW (Sandemier Software) and 2) spatially aligning, resampling, amplitude scaling, and summing to produce a composite dataset using algorithms created in the computer program MATLAB. The datasets for each antenna frequency require separ ate directories in ReflexW. The Mala GPR system used in this study saves each 2D profile as a radan file. Profiles are brought into ReflexW by importing radan files using the original file name. Upon completion of the import process, the imported profil es are placed into separate folders within the project directory on the basis of being an odd numbered profile or an even numbered profile. The processing flow for one profile can be saved and applied to the other 2D profiles that make up a 3D survey. Eve n and odd numbered profiles are separated into different folders because slightly different processing flows will be applied to each folder. Processing starts by opening any one 2D profile for a particular survey. The processing steps for each frequency dependent dataset are shown in Table 2. Static corrections involve clipping data prior to the first arrival of reflected energy. The dewow filter subtracts the very low frequencies from the data. AGC (automatic gain control) aids in bringing out weaker ref lections to the same amplitude as stronger reflections within a certain time window. A time cut was used to clip data after a subjectively identified two way travel time beyond which noise exceeded signal
25 Table 2. Summary of the basic processing steps a pplied to GPR data in ReflexW. Process 100 MHz 250 MHz 500 MHz 800 MHz Static Correction 17.75 ns 17.75 ns 17.75 ns 17.75 ns Dewow Filter window length=10ns window length=4ns window length=2ns window length=1.25ns AGC Gain window length=5ns window leng th=5ns window length=2ns window length=1.5ns time cut 180 ns N/A N/A N/A Bandpass Filter (cos taper) corner frequencies 10 50 200 240 MHz corner frequencies 10 100 1000 1100 MHz corner frequencies 50 100 1000 1100 MHz corner frequencies 10 300 1500 2000 MHz After fully processing a single radargram, the processing flow is saved and batch processing can be used to apply the same processing steps to the other radargrams within the single frequency dataset. The processing flow is first applied to the files in the folder where even numbered profiles are located. Since every profile was acquired in the opposite direction as the one before, either all even or odd numbered profiles must be flipped in the x direction. All odd profiles were fl ipped in the x direction so that all profiles appear to be acquired north to south by creating a separate processing flow. The new processing flow that includes the profile flip is applied to the files in the folder for odd numbered profiles. At this point fully processed 2D profiles are placed in the folder PROCDATA. Creating a 3D dataset from several 2D profiles is done with ReflexW in 3D Interpretation Mode. A 3D dataset for each single frequency survey is created and then exported in the form of an ASC II 3 column file. The 3D datasets are viewed by scrolling through the x, y, and z planes of the data cube. These files can then be read into MATLAB where the radargrams are further processed into composite radargrams.
26 Chapter 3: Composite Radargrams Introduction A composite radargram is composed of radargrams from single frequency antennae summed together into a single radargram. Integrating data from both high and low frequency antennae partially resolves the comp romise of range resolution trade off. Although composite radargrams do not provide high spatial resolution at greater depths, shallow and deep structures may be imaged simultaneously with a resolution determined by the highest antenna frequency capable o f imaging that particular structure. Although composite radargrams provide a remarkably more complete view of the subsurface, relatively few workers have attempted to create composite radargrams, and even fewer have been successful. Dougherty et al. (1994) are unable to optimize the contribution of high frequency data in a composite radargram. The relative amplitudes of higher frequency radar signals are smaller than lower frequency radar signals at comparable depths. Summing these signals without any ampli tude scaling will result in little improvement of spatial resolution. Thompson (2000) is unable to justify the increased processing time relative to the marginal improvement made to the composite radargram. A more theoretical approach (Noble 2001) determi ned that amplitude scaling and time shifting the peaks of wavelet pulses to maximize constructive interference should produce a greatly improved composite radargram.
27 Figure 12. Using unscaled amplitudes results in a composite dominated by low frequency energy. Adapted from Strategies for Building a Composite Ground Penetrating Radar (GPR) Profile by Adam Booth.
28 Figure 13. Applying time shifts to scaled amplitude synthetic wavelets results in a composite wavelet that more closely resembles a delta function, the ideal wavelet for deconvolution. Adapted from Strategies for Building a Composite Ground Penetrating Radar (GPR) Profile by Adam Booth. Booth et al. (2009) thoroughly treats the subject of composite radargrams from both the theoret ical basis used by Noble (2001) and a practical perspective using an actual case study. Composite radargrams can also be used to improve the effect of deconvolution of GPR data. Deconvolution is applied more often to seismic data than GPR data because of the limited spectral bandwidth of the GPR wavelet. The summation procedure used by Booth et al. (2009) is focused on broadening the spectral bandwidth of GPR data from one octave to two octaves. Regardless of the motive for creating composite radargrams, t he fundamental steps in creating a composite radargram remain the same. 1) Basic processing steps (static correction, dewow filter, gain, spectral filter, time cut) -4 -2 0 2 4 6 8 10 12 0 10 20 30 40 50 Time (ns) Amplitude 200MHz wavelet; scaled by factor of 1; timeshift from 0ns = +15ns 100MHz wavelet; scaled by factor of 0.25; timeshift from 0ns= +10ns 50MHz wavelet; scaled by factor of 0.0625; timeshift from 0ns = +0ns Composite wavelet
29 2) Spatial coincidence (representing each single frequency dataset with a common set of spatial arrays through interpolation) 3) Amplitude scaling 4) Applying time shifts as in Figure 10 (not deemed as necessary by Booth et al. (2009)) 5) Summation There are multiple methods used to evaluate the procedures used to create composite radargrams. When composit e radargrams are used to optimize the features shown by each individual antenna, the ability to capture as many subsurface reflections at optimal resolution provides the main basis for evaluation. When composite radargrams are used to widen the spectral b andwidth to provide more resolving power for deconvolution, time varying frequency analysis can be used to determine whether the spectral bandwidth reaches the goal set by Booth et al. (2009) of two octaves. Additionally, summation should not result in the introduction of random elements due to interference of out of phase radar signals. Although the theoretical treatment of composite radargrams involves time shifting pulses to eliminate out of phase signals, in practice applying dynamic time shifts is oft en difficult and not essential in creating high quality composites (Booth et al. 2009). Spatial Coincidence ReflexW offers a quick and easy way to apply standard processing techniques to GPR data (dewow, mute, filter, gain, migration, etc). More advanced or customized processing techniques require the use of programs such as MATLAB that allow users to
30 transform GPR data into matrices, vectors, or time series. Data in this form can be manipulated with various techniques of digital signal processing and ma trix algebra. Creating composite radargrams involves summing some form of the matrices that represent the amplitude of energy received by the GPR antenna. A prerequisite to summing is that these matrices must be the same size and spatially coincident wi th each other. Although 2 D GPR images technically only have one spatial dimension, x; time (t) can also be thought of as a spatial dimension because t can be converted to depth using subsurface velocity. The steps taken to achieve spatial coincidence inc lude 1) shift the position of the 800 MHz antenna, 2) correct for 2D profiles that were too long, 3) fix odometer calibration, 4) extract a subset of data for 100 MHZ and 250 MHz data sets, and 5) re interpolate data to the sampling interval of the 800 MHz data set. 1) Shift the position of the 800 MHz antenna One of the advances in GPR technology that makes multiple antenna surveys feasible is the capability of the control unit to record data from two different antennae simultaneously. The Eglin AFB survey i mplemented this time saving advantage through simultaneous acquisition with the 500 MHz and 800 MHz antennae. The 800 MHz antenna was trailing behind the 500 MHz antenna so that the center of each antenna was separated by a distance of 41.3 cm (the 100 MH z and 250 MHz antennae were used separately due to their bulkiness). Although simultaneous acquisition cut survey times down by several hours (approximately 3 hours for a 3 D survey built from about 100 2 D profiles), extra processing is required to make t he distance arrays for each data set to be spatially coincident. The antenna offset (41.3 cm) must be subtracted from the distance array of the trailing antenna. The data array must be manipulated by cropping the
31 appropriate number of columns that corres pond to the space behind the starting position of the other three antennae. 2) Correct for 2 D profiles that were too long Great care was taken to ensure that each 2 D profile was acquired strictly within the boundaries of the survey. However, the element o f human effort introduces the inevitable reality that some profiles may be a couple centimeters too long or too short. Even errors of such a small magnitude produce a jittery visual effect on time slices within the 3 D datasets. In an effort to minimiz e this jitter, processing steps were taken to shift and crop certain 2 D profiles. After each 2 D profile was acquired, the antenna was carefully repositioned to the starting point of the next profile while the computer was saving the last acquired prof ile and creating a new file for the next acquisition. This survey procedure assumes that the starting position for each profile is correct and that the spatial error is associated with the ending position. Another relevant survey procedure is that acquis ition occurred in a mow the lawn type pattern. The first profile and every other one after that was acquired south to north, and the second line and every other one after was acquired north to south. During processing with ReflexW the odd numbered prof iles were flipped so that each profile would run north to south. Also, the distance array (x) measures distance from the northern boundary of the survey. This distance array is calculated by ReflexW from the length of the longest 2 D profile. Initially, zeros are used to pad shorter profiles so that a matrix with no absent elements can represent a time slice. Each profile is placed so that it starts at the northern survey boundary, x=0. The first step in the procedure is to determine the correct survey length in terms of array elements. As an example, for the 100 MHz dataset, the initial survey consisted of 456
32 elements, in the x direction, and the correct survey length consists of 455 elements (equivalent to the survey length measured by the odometer wheel, 22.7 m). Since the assumption is that each even numbered profile starts at x=0, these profiles need not be shifted. Instead, these profiles are either cropped or padded with extra zeros so that they have the correct survey length of 455 elements. The next step in the procedure is to shift the odd numbered profiles so that they all end at element position 455. This is done by locating the last non zero element in each profile. The difference between this element number and 455 is the amount that each profile is shifted. Looping this procedure for each time slice results in a jitter corrected data cube. Table 3. A subset of a 100 MHz time slice before jitter corrections. The southern grid boundary is represented by row 455. 1 2 3 4 5 444 31820 31354 32175 29959 30533 445 31534 31369 32035 30048 30660 446 31937 31327 32008 29954 30348 447 32462 31456 31070 29200 30861 448 31447 30833 30858 29669 30299 449 32442 0 30871 28838 30598 450 31243 0 31862 2911 1 30866 451 32060 0 31226 29207 29993 452 32447 0 31189 30846 29537 453 31488 0 31221 0 30427 454 31208 0 31422 0 30170 455 0 0 31535 0 0 456 0 0 31727 0 0
33 Table 4. A subset of a 100 MHz time slice after jitter correcti ons. The southern grid boundary is represented by row 455. 1 2 3 4 5 444 31820 31097 32175 30922 30533 445 31534 31847 32035 31436 30660 446 31937 31161 32008 30749 30348 447 32462 31715 31070 29959 30861 448 31447 31079 30858 30048 30299 449 32442 31721 30871 29954 30598 450 31243 30727 31862 29200 30866 451 32060 31354 31226 29669 29993 452 32447 31369 31189 28838 29537 453 31488 31327 31221 29111 30427 454 31208 31456 31422 29207 30170 455 0 30833 31535 30846 0 456 Unfortunately, this procedure proved to be ineffective in improving the image quality of time slices. The jitter effect was actually increased for the higher frequency data sets, 500 and 800 MHz. This can be show n by comparison of Figures 15 and 16. The jitter is likely due to acquiring profiles of incorrect length and instances when the odometer did not have direct contact with the surface. Accounting for the odometer wheel losing direct contact requires rescal ing the profiles at necessary locations. Poor knowledge of these locations do not allow for the jitter effect to be fully addressed.
34 Figure 14. Time slice of 500 MHz data at 14.7 ns, before jitter corrections were applied.
35 Figure 15. Time sli ce of 500 MHz data at 14.7 ns, after jitter corrections were applied. Clearly, improved algorithms are necessary for clearer time slices. This work is beyond the scope of this thesis. 3) Fix odometer calibration The odometer wheel must be calibrated to the surface conditions of the survey site. Although the odometer wheel was calibrated prior to the survey, the calibration settings did not work perfectly. This issue was corrected by creating new rescaled distance arrays that started at 0 and went to the c orrect survey length with the appropriate linear spacing that would preserve the number of elements in the original distance array. Table 5 shows that the odometer wheel overestimated the survey length.
36 Table 5. This table shows the inaccuracy of the odo meter wheel calibration setting for each survey. Note that the 500 MHz and 800 MHz surveys were acquired simultaneously. Antenna Frequency (MHz) Survey Length measured by Odometer Wheel (m) Actual Survey Length (m) 100 22.7 22 250 22.75 22 500 10.32 10 800 10.32 10 4) Extracting the subset of the survey space common to each individual survey The dimensions of the surveys done with the 100 MHz and 250 MHz antennae were 22 m, north to south, and 10 m, east to west. A subset of this survey was selected fo r the 500 MHz and 800 MHz survey because 1) using higher frequency antennae requires an especially small spacing between 2D profiles to avoid spatial aliasing (0.05 m) and 2) using a smaller survey space did not compromise the geologic value of the data. The northeast corner of the original survey measuring 10 m, north to south, and 5 m, east to west was used. The 100 MHz and 250 MHz datasets were cropped to be spatially coincident with the 500 MHz and 800 MHz datasets. Cropping in the x direction was d one in MATLAB to retain only the 10 northernmost meters of the original surveys. Cropping in the y direction was done in ReflexW. 5) Interpolation Even when the datasets for each antenna represent the same space, the matrix sizes are still different due t o the different line spacing between 2D profiles, the different sampling rates for each antenna, and the different trigger interval for each
37 antenna. These issues can be overcome by interpolating all arrays so that each dataset can be represented by a com mon set of spatial and time arrays. Creating these common arrays is based on the premise that the time increment should be based off the data set with the highest sampling rate. The other data sets that use lower sampling rates are re interpolated to have a time increment equivalent to the highest sampling rate. The spatial arrays for the 800 MHz dataset were chosen to represent the common spatial arrays because this frequency antenna survey required smaller line spacing (0.05 m) and trigger interval of 0.0 25 m. The time arrays have various sampling intervals and maximum values. The different maximum values of each time array are due to the different penetration ability of each antenna. Acquisition parameters were set to stop acquisition beyond the expecte d depth of penetration. The common time array was created by using the sampling interval of the 800 MHz dataset and the maximum value of the 100 MHz dataset. Once the common time and space arrays have been formed, the amplitude arrays can be interpolated to match the common arrays. Although the amplitude arrays are composed of time series that resemble periodic functions much more than a linear function, linear interpolation is sufficient to use on non linear time series for sufficiently small intervals. Spectral Bandwidth Broadening Summation Techniques Building composite radargrams involves taking spatially coincident radargrams and summing them together so that the properties of each single frequency radargram are optimized. The three summation techniques that have been developed by previous workers can be described as 1) brute summation (Figure 17), 2) coarse mute (Figure 20),
38 and 3) ramped mute (Figure 22). These three techniques are used to broaden the spect ral bandwidth of GPR data. Brute summation, used by Dougherty et al. (1994), involves summing nominal frequency radargrams with no modifications. Unless otherwise noted, all 2D GPR radargrams represent data taken along transect shown in Figure 16. Each GP R acquired a profile directly along this transect. Figure 16. Composite grid (see Figure 11) showing the location of 2D profile used for summation. Position along the transect is measured as the distance from the northern terminus of profile. The profile length is 9.587 m due to the trailing position of the 800 MHz antenna.
39 Figure 17. Composite radargram create by brute summation. The dominance of low frequency (Figure 18) energy minimizes the contribution of high frequency energy (Figure 19) as obser ved by Dougherty et al. (1994).
40 Figure 18. 100 MHz radargram. The direct wave arrival and horizontal reflections directly below can be clearly identified in Figure 17. ` Figure 19. 800 MHz radargram. Higher frequency energy provides enhanced spatial re solution not well displayed in Figure 17.
41 Coarse mute summation (Figure 20) addresses the issue of low frequency direct waves overwhelming the presence of high frequency energy by zeroing the low frequency direct waves. However, coarse muting introduces an undesirable visual attribute to composite radargrams. This visual attribute can be described as an abrupt transition from zeroed data to fully scaled data. abrupt transition Figure 20. Composite radargram created by coarse muting. Muting low frequency energy improves spatial resolution from 0 20 ns. The abrupt introduction of 100 MHz energy can be clearly identified at 20 ns.
42 Ramped summation is the solution to eliminating an abrupt transition. In addition to zeroing out low frequency dir ect waves, a transition zone is established that gradually introduces lower frequency data into the composite diagram (Figure 21). The ramp from 0% wave energy to 100% wave energy can be represented by a linear function or smoother filters such as the Butt erworth function (used only for resolution and depth enhancement). The idea of using a Butterworth function arises because the ramp acts as a time domain filter. Butterworth filters may be favored over trapezoidal filters, the result of using a linear ramp because they are less likely to introduce artifacts due to having no sharp corners (Figures 33 and 35). A composite radargram built with linear ramped summation is shown in Figure 22. Figure 21. From left to right, 500 MHz, 250 MHz, and 100 MHz radar grams. The direct wave of each radargram is muted, and early arriving energy is gradually introduced with a linear ramp as used by Booth et al. (2009) to reduce the abrupt transition associated with coarse muting.
43 smoother transition Figure 22. Com posite radargram created by ramped summation. The enhanced spatial resolution of high frequency energy is better preserved, and the transition to lower frequency energy is less abrupt compared to Figure 20. Similar radargrams are used by Booth et al. (2009 ) for deconvolution.
44 Amplitude Balancing The overwhelming presence of low frequency data when summed with higher frequency data as described by Dougherty et al. (1994) confirms this relationship. Therefore some form of amplitude scaling is required to ba lance the contribution of each frequency component within a composite radargram. The two methods of amplitude balancing described here are 1) dominant frequency amplitude equalization (DFAE), and 2) optical spectral whitening (OSW). Only OSW is applied to the data. Dominant frequency amplitude equalization, DFAE, uses the maximum amplitude spectrum value associated with each nominal frequency data set to calculate balancing multipliers (Figure 23). Not considering a broader portion of the amplitude spectru m may be problematic if the spectrum is asymmetric, because the useful energy may not be associated with the dominant frequency.
45 A a. b. Figure 23. DFAE for synthetic wavelets. A) Synthetic GPR pulse created from a Berlage wavelet. B) The normalized amplitude spectra of synthetic GPR pulses. The maximum value for each spectrum is used as a DFAE multiplier. Adapted from Booth et al. (2009). The objective of OSW is to produce a frequency spectrum over d iscrete time windows that resemble white light, equal amplitude values for all frequencies. OSW produces an array of time varying multipliers for each nominal frequency data set that will produce a composite radargram whose spectrum will match some optimu m spectrum, usually a white spectrum, through least squares analysis algorithm (Figure 24). A time varying Fourier transformation of each nominal frequency data set must be performed prior to applying the least squares algorithm.
46 Figure 24. OSW least sq uares algorithm in matrix notation. The amplitude (A) multiplied by an OSW weight (w) results in some desired spectrum (s). The unknown quantity w can be expressed as, w = (A T A) 1 A T s. Booth et al. (2009) uses OSW to make deconvolution more effective for G PR data. Deconvolution is not usually applied to GPR data because the spectral bandwidth of a nominal frequency data set is approximately one octave. OSW allows for a composite radargram to have a spectral bandwidth of nearly two octaves with an equal con tribution of each frequency within the bandwidth. The effectiveness of the OSW method can be evaluated using the S transform, a MATLAB algorithm that calculates the time varying amplitude spectrogram of a trace. The S transform for a nominal frequency trac e should expectedly have a narrow bandwidth no more than one octave centered at the antennas central frequency (Figures 25 28). The S transform for an OSW composite (Figure 29) should have a bandwidth of nearly two octaves (spanning from the lowest to hig hest center frequencies used).
47 Figure 25. From top to bottom, single trace of 100 MHz data, S transform of trace, normalized S transform. high amplitude low amplitude
48 Figure 26. From top to bottom single trace of 250 MHz data, S transform of trace, normalized S transform. high amplitude low amplitude
49 Figure 27. From top to bottom, single trace of 500 MHz data, S transform of trace, normalized S transform. high ampl itude low amplitude
50 Figure 28. From top to bottom, single trace of 800 MHz data, S transform of trace, normalized S transform. high amplitude low amplitude
51 Figure 29. From top to bottom, single trace of composite data (single ramped summation with OSW balancing), S transform, n ormalized S transform. When the motivation for producing a composite radargram is to optimize the depth and penetration of nominal frequency data sets, as opposed to broadening the spectral bandwidth, OSW may not always be a desirable method for amplitud e balancing. Since each OSW multiplier is determined from the amplitude spectra of each frequency, the method is limited by the penetration capability of the highest frequency antenna. The composite radargram built by Booth et al. (2009) terminates at 80 ns and focuses on improving resolution within the overlapping time period of all nominal frequency data high amplitude low amplitude
52 sets. Figure 29 shows that the OSW method was able to increase the spectral bandwidth from 0 60 ns. Deconvolution Booth et al. (2009) applies deconvolu tion to a composite radargram built with single ramped summation and balanced with OSW. Applying OSW, following the method described by Booth et al. (2009), to nominal frequency radargrams over long time segments results in Figure 30. The sum (single ramp ed summation) of these nominal frequency radargrams results in the composite radargram used for deconvolution (Figure 31). Figure 30. Nominal frequency radargrams after OSW multipliers have been applied. Clockwise from top left, 800 MHz, 500 MHz, 100 MHz, and 250 MHz.
53 Figure 31. The sum of radargrams from Figure 30.
54 Figure 32. Deconvolved radargram from Figure 31. Deconvolution parameters required by ReflexW: autocorr start 0, auto corr end 12, filter length 2, white noise 5%.
55 Resolution and Penetration Enhancement Booth et al. (2009) uses a linear ramped summation method to gradually introduce lower frequency energy. Since the goal of Booth et al. (2009) is to broaden the spectral bandwidth of GPR radargrams, no attempt is made to suppress higher frequency data as it is attenuating. For the purposes of creating a composite radargram that utilizes high frequency energy at shallow depth and lower frequency data at greater depth, I propose the use of an additional ramp that gr adually suppresses high frequency data by the same amount that lower frequency data is being introduced, i.e., double ramped summation. Double Ramped Summation There are three parameters of the double ramp summation filter that may vary 1) the function use d to create the ramps (e.g. Butterworth (Figures 33 and 35) or linear (Figure 37)), 2) the length of time over which ramping takes place, and 3) the time at which a particular nominal frequency data set has attenuated sufficiently to start suppressing it.
56 Figure 33. Double summation ramps using a two pole Butterworth function. A composite radargram created by double ramped summation using four nominal frequency data sets requires three summation zones. The ramps that introduce energy are created with a two pole Butterworth function, and the suppressing ramp is set to 1 (Butterworth function) so that the two ramp values sum to 1 at any given time. A) The first set of ramps suppresses 800 MHz energy and introduces 500 MHz energy B) The second set of ramps suppresses 500 MHz energy and introduces 250 MHz energy. C) The final set of ramps suppresses 250 MHz energy and introduces 100 MHz energy. The radargram produced from these double summation ramps is shown in Figure 34. a. b. c.
57 F igure 34. Composite radargram produced with double ramped summation. The ramps are created with a 2 pole Butterworth function.
58 Figure 35. Double summation ramps using a six pole Butterworth function Using six poles increases the functions range to one. However, the effective length of the ramps is decreased. a) The first set of ramps suppresses 800 MHz energy and introduces 500 MHz energy. b) The second set of ramps suppresses 500 MHz energy and i ntroduces 250 MHz energy. c) The final set of ramps suppresses 250 MHz energy and introduces 100 MHz energy. The radargram produced from these double summation ramps is shown in Figure 36. a. b. c.
59 Figure 36. Composite radargram produced with double ramped summa tion. The ramps are created with a six pole Butterworth function. The visible transitions result from the reduced summation zone of the six pole Butterworth function.
60 Figure 37. Double summation ramps using linear functions. a) The first set of ramps suppresses 800 MHz energy and introduces 500 MHz energy. b) The second set of ramps suppresses 500 MHz energy and introduces 250 MHz energy. c) The final set of ramps su ppresses 250 MHz energy and introduces 100 MHz energy. The radargram produced from these double summation ramps is shown in Figure 38. a. b. c.
61 Figure 38. Composite radargram produced with double ramped summation. The ramps are created with linear functions. A more objective approach to determining ramp length is to measure it in terms of the period of the suppressed frequency (Table 6). The ramp length is also guided by the principle that the length should be long enough in order to create a non abrupt trans ition from higher to lower frequency energy. When this condition has been met, the ramp length should not be further extended because spectral bandwidth broadening absent of deconvolution does not improve the composite radargram in terms of geologic inform ation. Figures 39a and 39b show composite radargrams with a ramp length that is too short and one that better conveys radar stratigraphy at the transition zone, respectively.
62 Table 6. The period for each antenna center frequency is calculated to determi ne the ramp length chosen to be 15 wave periods. Frequency (MHz) period (ns) ramp length (ns) 100 10 150 250 4 60 500 2 30 800 1.25 18.75 a. b. Figure 39. A comparison between composite radargrams produ ced using double linear ramps with lengths of 2 and 15 periods of outgoing frequency. a) 2 periods, b) 15 periods.
63 Selecting a time at which to begin suppressing data can be achieved in a more objective fashion by examining the amplitude envelope. The exp onential decay curve associated with the amplitude envelope of electromagnetic waves in lossy media is obscured when working with AGC gained data. An appropriate amplitude envelope can be created with data that only has been minimally processed (dewow filt er and static correction). The envelope is calculated as the absolute value of the Hilbert transformation of a single trace. The envelopes for each trace in a profile are then averaged to create a single envelope to be evaluated. The general guideline for choosing the time at which suppression starts is that taking the log of the amplitude envelope (Figure 40). This results in a curve whose minimum value correlates well with the point in time when the received energy is too attenuated to provide useful info rmation, shown in Figures 41 43.
64 250 MHz 100.7 ns 500 MHz 60.1 ns 800 MHz 40.9 ns Figure 40. The log amplitude envelope of ungained data plotted against time. The time value corresponding the minimum value of t he curve can be used as the time for starting suppression. a) 250 MHz b) 500 MHz c) 800 MHz a. c. b. c.
65 Figure 41. 250 MHz (AGC gained) radargram with black solid line at time when suppression begins.
66 Figure 42. 500 MHz radargram (AGC gained) with black solid l ine at time when suppression begins
67 Figure 43. 800 MHz radargram (AGC gained) with black solid line at time when suppression begins
68 Amplitude Envelope Equalization AEE involves calculating multipliers that will result in approximately equal ampl itude envelopes for each nominal frequency data set after AGC gain has been applied (Figure 44). Each AEE multiplier is calculated as the ratio of the average envelope value of the lowest frequency data set and the average envelope value of the nominal fre quency data set for which the AEE multiplier is to be applied (Table 7). Since an amplitude envelope is calculated from a single trace of a radargram, the envelopes for each trace are averaged to create a single amplitude envelope for AEE. The only part of the envelope considered for AEE corresponds to the period of time that is not suppressed in the composite radargram, i.e. the amplitude envelopes for the lower frequency direct waves are not considered. A composite radargram balanced by AEE multipliers is shown in Figure 45.
69 average envelope value = 27096 average envelope value = 19039 a. b.
70 average envelope value = 17092 average envelope value = 16864 Figure 44. The amplitude envelope for AGC gained data and average envelope value plotted against the time period w hen the nominal frequency data is not suppressed within the composite. a) 100 MHz, b) 250 MHz, c) 500 MHz, d) 800 MHz. c. d.
71 Table 7. Each average envelope value is divided by the highest average envelope value, resulting in a multiplier value for each nominal frequency dataset. Frequency (MHz) Average envelope value Multiplier value 100 27096 1.00 250 19039 1.42 500 17092 1.59 800 16864 1.61 Figure 45. Composite radargram built with double ramped summation and AEE amplitude balancing.
72 Subjectiv ely Chosen Values Another way to choose the values used in summation and amplitude balancing, e.g. the length of the summation ramp, the starting point of the summation ramp(s) (time to begin suppressing higher frequency data) is to determine them subject ively through visual inspection of each nominal frequency data set. Since the overall goal of radargram composites is image improvement, the calculated values and subjectively chosen values should be similar. Table 8 shows a comparison of the two sets of v alues when the subjective values were chosen with no prior knowledge of the calculated values. Table 8. Chosen values compared to calculated values for composite radargram with double summation and AEE balancing. Chosen value Calculated value 1st ramps start 14 ns 40.9 ns 1st ramps end 50 ns 59.7 ns 2nd ramps start 60 ns 60.1 ns 2nd ramps end 90 ns 90.1 ns 3rd ramps start 110 ns 100.7 ns 3rd ramps end 150 ns 160.7 ns end of data 180 ns 180 ns 100 MHz coeff. 1 1 250 MHz coeff. 1.25 1.42 500 MHz c oeff. 1.575 1.59 800 MHz coeff. 1.75 1.61 Comparison of chosen values to the calculated values reveals that most of the listed parameters can be determined automatically and still give results similar to what looks best to an interpreter. The largest d ifference between a chosen and calculated value is for the start time for the first set of ramps, 14 ns chosen versus 40.9 ns calculated. Composite radargrams built with calculated and chosen values are shown in Figures 46 and 47, respectively.
73 Figure 46 Composite radargram built with calculated values from Table 8.
74 Figure 47. Composite radargram built with subjectively chosen values from Table 8. The earlier introduction of 500 MHz energy shows more detailed reflections from 30 50 ns.
75 Chapter 4: Geologic Interpretation Velocity Analysis Until this point all GPR radargrams have had a vertical axis label of time (ns). This time represents the two way travel time taken for energy to be emitted from the transmitting antenna, reflect ed in the subsurface, and reach the receiving antenna of the GPR. Two way travel time must be converted to depth in order to correlate reflections with ground truthing methods, e.g. core, auger, or trench. The subsurface velocity of GPR waves can be de termined from a common midpoint (CMP) gather (Figure 48). A CMP survey involves separating the transmitting and receiving antennae by fixed distances while maintaining a constant midpoint between the two antennae. CMP surveys assume incident energy is refl ecting off the same subsurface point(s) each time. The change in two way travel time as a function of source receiver distance can be fitted to a parabolic function that represents a unique velocity value.
76 Figure 48. CMP gather with a best fitting para bola corresponding to a velocity of 0.10 m/ns. Although the CMP survey did not indicate the presence of more than one velocity layer, it is reasonable to suspect that there is a velocity change corresponding to the change from unsaturated sediment to sat urated sediment. The diffraction hyperbola shown in Figure 49 indicates that the EM wave velocity for shallow depths is 0.12 m/ns. It is likely that this velocity decreases with depth as the sediment becomes fully saturated. There were no visible diffra ctions to confirm this. Therefore, the velocity value of 0.10 m/ns was used as an average value for the time to depth conversion.
77 Figure 49. 800 MHz radargram with a fitted diffraction hyperbola corresponding to an EM wave velocity of 0.12 m/ns. Field Verification A composite radargram is only useful if it can provide a more efficient image for accurate interpretation. The deconvolved composite radargram built by Booth et al. (2009) shows a stratigraphic discontinuity previously unidentifiable. Howev er, there was no ground truthing that could confirm the observed discontinuity. Although the composite radargrams developed for the Santa Rosa Island data set image both shallow and deep structures, the field verification methods used, pound core (Figure 4 9) and trenches (Figures 50 52), were unable to penetrate more than 1.5 meters below the surface, corresponding to approximately 30 ns of two way travel time. Direct observations of the water table, pre Ivan dune surface, and washover deposits allow for th e uppermost 30 ns of the composite radargram to be geologically interpreted (Figure 53).
78 heavy mineral layer algal mats water table peat layer Figure 50. Core SR 1. Location shown in Figure 11. heavy mineral layer algal mats water table peat layer
79 Figure 51. Trench 2 looking west. The water table coincides with a layer of pe at interpreted as the pre Ivan surface, 1.2 m below land surface.
80 Figure 52. Trench 1 looking north (landward). Exposing the trench to the sea breeze exposes numerous thin layers of sand with varying resistances to erosion.
81 Figure 53. Trench 1 loo king east. Preserved dune vegetation within overwash deposits.
82 1 2 3 4 Figure 54. Composite radargram with structures with interpreted radar facies outlined in red. 1) landward dipping foreset beds observed in Trench 1, 2) outlin e of antecedent dune indicated by preserved sea grass in Trench 1, 3) Water table reflection and pre Ivan surface at 1.2 m observed in Trench 2, 4) Radar facies boundary (combination of horizontal and subhorizontal reflections above and mostly horizontal r eflections below. This reflection may correlate with overwash deposits from Hurricane Opal (Horwitz, personal communication 2010).
83 3D Data Composites Although expensive in terms of acquisition and computation, three dimensional surveys are requ ired to observe the geometry of sedimentary structures imaged by GPR. This can be done either by comparing multiple 2D profiles (Figures 54 58), as seen in all GPR radargrams up to this point, or in time slices (Figures 59 63). The composite radargrams sh own in this section use double summation ramps with subjectively chosen ramp lengths, cutoff times, and amplitude balancing multipliers. Figure 55. Composite radargram of profile 1 meter west of eastern boundary.
84 Figure 56. Composite radargram of profi le 2 meters west of eastern boundary. Figure 57. Composite radargram of profile 3 meters west of eastern boundary.
85 Figure 58. Composite radargram of profile 4 meters west of eastern boundary.
86 Figure 59. Composite radargram of profile 2 meters west of eastern boundary.
87 Figure 60. Depth slice at 0.2 meters. The antecedent dune (at bottom) is easily distinguished from overwash deposits (top).
88 Figure 61. Depth slice at 1.5 meters. In addition to dune and overwash deposits, the pre Ivan back barr ier marsh becomes visible in top left (supported by core data).
89 Figure 62. Depth slice at 2 meters.
90 Figure 63. Depth slice at 4 meters.
91 Figure 64. Depth slice at 8 meters. Three dimensional variations may be more difficult to observe with multi ple 2D profiles. However, time slices show some of this variation more obviously. Where depth slices can be supported by ground truthings, the geometry of observed structures can be imaged very well. Below 2 meters, different radar facies are still observa ble. however geologic interpretations of them can only be speculative.
92 Chapter 5: Conclusions The process of building composite radargrams varies significantly depending on the ultimate goal of the composite. These goals can be summarized as 1) broadeni ng the spectral bandwidth to increase the effectiveness of deconvolution, and 2) merging multiple frequency data sets together to optimize the spatial resolution of high frequency energy and the deep penetration of low frequency data sets. The process of s pectral bandwidth broadening incorporates amplitude balancing via optical spectral whitening (OSW) and summation using single ramps to produce a composite radargram whose wavelets more closely resemble the delta function (Booth et al., 2009). The theoreti cal example of this method works very well on isolated GPR wavelets. However, the effectiveness of this method is decreased when closely spaced reflections produce multiple high frequency wavelets within the period of one lower frequency wavelet. Forward models may increase our understanding of how the degree of isolation between reflections effects the creation of spiky wavelets through summation. Accomplishing the second goal does not require an amplitude balancing scheme with the goal of spectral broade ning. When working with AGC gained data, simple time independent amplitude balancing is able to produce well balanced radargrams that illustrate key geologic phenomena. Removing the goal of spectral bandwidth broadening calls for a more conservative summat ion scheme. Not only are direct wave arrivals removed from lower frequency data, lower frequency data is not introduced into the
93 composite radargram until the next highest nominal frequency data set shows signs of attenuation. Also as lower frequency data is being introduced, the attenuating higher frequency data is being removed or suppressed by the same amount. Analysis of the amplitude envelope of ungained data allows for the summation ramp starting time to be calculated automatically. The use of autom atically calculated values for ramp length, ramp starting point, and amplitude balancing coefficients, produces a high quality composite radargram. However, the visual attributes of a composite radargram can be optimized further by subjectively choosing th e ramp length, ramp starting point, and amplitude balancing coefficients based on visual inspection of each nominal frequency data set. Even when the calculated values are not used directly, they serve as a useful starting point from which adjustments can be made. Building a composite 3D radargram presents additional challenges. The difficulty of acquiring several 2D profiles of exact length and possible dragging of the odometer wheel results in time or depth slices with undesirable artifacts. Constant GPS measurements of the GPR during acquisition instead of using an odometer wheel should eliminate this issue. The computational expense of calculating amplitude balancing coefficients and ramp starting points leads to applying the results of one 2D profile t o an entire 3D data set. Assuming minimal 3D variation in these values produced reasonable results. Clearly, composite radargrams are capable of displaying more information than nominal frequency radargrams. Lower frequency antennae are traditionally us ed in coastal studies to get a basic sense of coastal stratigraphy, e.g., horizontal beds, dipping beds, and cross bedding. Higher frequency antennae are useful in capturing the detailed
94 geometry of recent depositional events. The value of composite radar grams can be further increased by field verification to depths imaged by lower frequency antennae.
95 References Annan, A.P., and Davis, J.L. 1977. Radar range analysis for geologic materials. Report of Activities, Geological Survey of Canada Paper 77 1B, 117 124. Baker, G.S., and Jol, H.M. 2007 Stratigraphic Analyses Using GPR: Geological Society of America Special Paper 432 Baker, G.S., Jordan, T.E., and Pardy, J. 2007. An introduction to ground penetrating radar (GPR). in Baker, G.S., and Jol, H. M., eds ., Stratigraphic Analyses Using GPR: Geological Society of America Special Paper 432, 1 18. Booth, A.D., Endres, A.L., Murray, T. 2009. Spectral bandwidth enhancement of GPR profiling data using multiple frequency compositing. Journal of Applied Ge ophysics 67, 88 97. Bristow, C.S. and Jol, H.M. (eds) 2003. An introduction to ground penetrating radar (GPR) in sediments Ground Penetrating Radar in Sediments Geological Society, London, Special Publications, 211, 1 7. Daniels, J., Roberts, R., and Ve ndl, M. 1995. Ground penetrating radar for the detection of liquid contaminants. Journal of Applied Geophysics v. 33, 195 207. Davis, J.L., and Annan, A.P. 1989. Ground penetrating radar for high resolution mapping of soil and rock stratigraphy. Geophy sical Prospecting 37, 531 551.
96 Dougherty, M.E., Michaels, P., Pelton, J.R., and Liberty, L.M. 1994. Enhancement of ground penetrating radar data through signal processing. Proceedings of the Symposium on the on the Application of Geophysics to Engineer ing and Environmental Problems Horwitz, M.H., and Wang, P. 2006. Sedimentological characteristics and internal architecture of two overwash fans from Hurricanes Ivan and Jeanne. CGAGS Transactions 55, 342 352. Horwitz, M.H. 2008. Sedimentological charac teristics of 3 D internal architecture of washover deposits from Hurricanes Frances, Ivan, and Jeanne. MS Thesis, the University of South Florida Jol, H.M., Smith, D.G. and Meyers, R.A. 1996. Digital ground penetrating radar (GPR): an improved and very eff ective geophysical tool for studying modern coastal barriers (examples from the Atlantic, Gulf and Pacific coasts, U.S.A. Journal of Coastal Research 12, 960 968. Leatherman, S.P. 1987. Coastal geomorphological applications of ground penetrating radar. Jo urnal of Coastal Research, 3(3), 397 399. Miller, D.L., Thetford, M., and Yager, L. 2001. Evaluation of sand fence and vegetation for dune building following overwash by Hurricane Opal on Santa Rosa Island, Florida. Journal of Coastal Research 17(4), 936 948. Noble, T. 2001. A theoretical basis for building a composite GPR profile. Unpublished BSc Honours Thesis University of Waterloo
97 Pries, A.J., Branch, L.C., and Miller, D.L. 2009. Impact of hurricanes on habitat occupancy and spatial distribution of beach mice. Journal of Mammalogy 90(4), 841 850. Stone, G.W. and Stapor, F.W. 1996. A nearshore sediment transport model for the Northeast Gulf of Mexico Coast. Journal of Coastal Research 12(3), 786 792. Stone, G.W., Liu, B., Pepper, D.A. and Wang, P. 2 004. The importance of extratropical and tropical cyclones on the short term evolution of barrier islands along the northern Gulf of Mexico, USA. Marine Geology 210, 63.78. Thompson, K.G. 2000. Developing a procedure to build a composite of multi frequenc y GPR data. Unpublished BSc Honours Thesis, University of Waterloo Wang, P., and Horwitz, M.H. 2007. Erosional and depositional characteristics of regional scale overwash deposits caused by multiple hurricanes. Sedimentology 54:3, 545 564.