Introduction

Seismic methods have been used successfully in water resource management for decades, usually to define significant geologic features such as bedrock (Haeni 1988), clay layer aquitards (Dickenson 1996), and faults which serve as preferential flow pathways (Shtivelman et al. 1998). In recent years there has been increased interest in the incorporation of this type of data into flow and transport simulations (Rubin et al. 1992, Copty and Rubin 1995), thus improving our ability to anticipate the nature of fluid flow.

The geologic constraints offered by seismic methods have proven useful, but the direct detection and geophysical "observation" of fluids has been elusive. Near-surface fluid detection is difficult because the wide range of geophysical properties of the host material typically exceed the perturbations caused by the presence of a fluid. In addition, the influences of material type, porosity and compaction map non-uniquely, along with fluid saturation, into geophysically detectable parameters. These two factors make it nearly impossible to discern the presence of fluids without a priori constrains.

Fluid distribution, however, can vary on short time scales. Changes in seismic response are more easily correlated with the presence of fluids than measurements from a single timeslice. Recently the exploration geophysics community has been using so-called 4-D geophysical techniques to observe these changes in petroleum reservoirs (Anderson 1998). The principle behind all 4-D techniques is similar - collect the "same" data sets at different times and analyze them for differences. There are several factors that make 4-D techniques well-suited to near-surface hydrogeologic applications: fluid movement in the near-surface often occurs on short time scales of hours to weeks; data can be collected relatively inexpensively; and "ground truthing" can often provide critical feedback to experimental design.

However, there are two hurdles which currently limit the effectiveness of temporal geophysical data. The first is a poor understanding of how materials respond to saturation. Theories describing the seismic velocities of partially-saturated granular media date back to Biot (1956). Lab experiments seeking to prove these relationships vary considerably depending on the type of material, the pressure applied to the system, and the distribution of fluid in the pore spaces (Domenico 1976 and 1977, Knight and Nolan-Hoeksema 1990). But a few trends appear consistent. Both compressional and shear wave velocities decrease as moisture is introduced into a granular matrix. As maximum saturation is approached, the compressional velocity increases dramatically while the shear velocity continues to decrease. But few field studies have been able to corroborate these experiments which were performed at megahertz frequencies on small clean samples or even glass beads. The role of compaction, small amounts of clay or soil, and the strong relative pressure gradients in the very shallow earth are still poorly understood.

The second impediment to time-dependent surveys is the lack of simple tools for analyzing geophysical data which vary with time. Studies of transient phenomena must move beyond the time-consumming and imprecise method of processing multiple timeslices independently and then comparing them for differences. In this paper we propose a signal correlation strategy which relies on the similarity of seismic waveforms from one timeslice to the next. Each timeslice is regarded as a small departure from a previously determined reference structure. Not only is this philosophy of interpretation far more precise, it also streamlines processing. For each timeslice, one need only solve for the perturbation to an already-determined Earth model.

Experiment Design

We chose a beach setting for this experiment because of the well-known temporal variations in subsurface moisture (Robinson et al. 1998) associated with tide cycles. The semi-diurnal changes allowed us to perform this experiment in one day. This eliminated errors resulting from replanting geophones and minimized the influence of changing surface conditions (Jefferson et al. 1998). While far from homogeneous, the beach provided some of the most uniform soil conditions available in a natural setting. The unconsolidated sands are typical of the environment of many unconfined aquifers and landfills where seismic surveys have been used extensively.

A section of beach was selected west of Indian Wells in East Hampton, New York. An array of twenty-four vertical component 40 Hz Mark Products geophones on 12 cm spikes were planted perpendicular to the shore at a depth of 15 cm (Figure 1). A Geometrics Strataview-R system recorded the seismograms. A short sample period of 0.0625 ms was selected because our correlation methods rely heavily on wave shape. All seismograms were 512 ms long with a 10 ms pretrigger window. The 24-bit dynamic range of the recording system allowed the simultaneous collection of high quality body and surface waves. A hammer and plate seismic source was chosen. This source offered the best combination of high frequencies, repeatability and repetition rate. Precise timing allowed us to stack 30 successive hammer blows in less than three minutes without deteriorating high frequencies. Not only did stacking average any variability from the hammer source, it also minimized contamination of the records by surf noise.

Beginning near high tide, seismograms from location A and B were recorded every half hour during one tide cycle for a total of 26 timeslices. At hour 8.5 (low tide), a complete set of five shots was recorded (Figure 1). Shots A and B were used to track variations in the complete structure obtained from the detailed reference survey at hour 8.5. One well near each end of the array monitored the actual water table. The wells were offset from the line by roughly 5 meters to avoid influencing the wave propagation. Well measurements began a few hours after the seismic surveying. As expected, the well nearest to the open water fluctuated more during the tide cycle. Fluctuations in the inland well were of smaller amplitude and were delayed by roughly one hour.

Figure 2 shows two vertical component seismograms from hour 7.5 and 12.5. The source is located at the 0 m position. Direct, refracted and Raleigh waves can be observed. Though similar, each record is systematically stretched and shifted relative to one another. To quantify these shifts, the two records of interest were split into small windows corresponding to each distance and traveltime combination. All signals outside of a small time window were tapered to zero to isolate the phase (Figure 5). After applying the identical windowing to each timeslice, the phase shift between records was calculated using a fourier-domain cross-correlation method. The time lag which yielded the highest correlation was considered the shift between the two records at each offset and traveltime.

The correlation plot (Figure 2) highlights areas of the wave field which are strongly coherent from one timeslice to the next. Smooth correlations between the compressional waves and the surface waves reflect the time-variant velocity structure. Noise dominates the lower left portion of the records creating spatially incoherent correlations. Compressional wave time shifts between the low and high tide signals are on the order of 1 millisecond. The largest timeshifts are associated with the Rayleigh surface waves. Delays of up to 4 ms can be observed. Based on this diagram, it is clear that both the compressional body waves and the surface waves contain valuable information about the time dependence of the velocity structure. Since each offers a different type of information, they are both investigated here.

Compressional velocities at the reference time

The set of five shots collected at hour 8.5 were used to determine the compressional velocity structure at the reference time. Direct and refracted compressional wave traveltimes were picked by hand and interpreted using standard refraction techniques. An appropriate model was developed using the RTMOD ray tracing program (Caress 1992). To keep the model descriptive and to avoid interpreting poorly resolved structures, we did not allow velocities to vary laterally within each layer. The sharp velocity contrast across the water table is visible in the traveltimes (Figure 3) and model (Figure 4) and follows similar observations by Michaels and Barrash (1996). The top 2 layers are used to replicate the small observed velocity gradient above the water table, presumably resulting from increased compaction with depth.

Compressional velocity changes through time

Once the reference model had been determined, the evolution of this structure during the tidal cycle was examined. Traveltime differences between each timeslice and the reference were calculated by cross-correlating the compressional waveforms. After applying identical windowing to each timeslice (Figure 5), the time lag which resulted in the best correlation with the reference signal was taken to be the traveltime difference. The evolution of refracted ray traveltimes from source A is shown in Figure 6. The traveltime variations of 1 ms (7%) are observed during the tidal cycle. Variations in the velocity structure were derived by modeling the changes in traveltime (Figure 4). Two transient features in Figure 4 are notable. The first is the change in the depth of the seismic discontinuity corresponding to the water table. The second is a decrease in velocities above the water table as high tide is approached. This same trend was observed by Bachrach and Nur (1998). They present a qualitative explanation relating the moisture content to the compressional velocity,

As fluid is introduced into a granular matrix, the density, rho, increases. The bulk modulus, kappa, and shear modulus, mu, remain unaffected because the interstitial fluid is poorly connected and has no shear strength. With this relation, the decrease in velocity above the water table can be explained by capillary forces which increase the saturation above the advancing water table.

Obtaining the dispersion relation at the reference time

Rayleigh surface waves, sometimes referred to as "ground roll", dominate the seismograms in this experiment (Figure 2). Unlike compressional waves, Rayleigh wave velocities are a function of frequency. Instead of arriving in a sharp pulse like body waves, they tend to disperse with distance. Because of this phenomenon, the frequency vs. velocity relationship is often referred to as a dispersion curve. While the surface waves at different timeslices appear related, their dispersive characteristics vary considerably. Each frequency responds differently to changes in the velocity structure. Rayleigh velocity changes cannot be determined from shifted waveforms as was done with the compressional waves. This technique (Figure 2) is only sensitive to the group velocity of a specific peak. But it is the phase velocity of each frequency which can be related to sub-surface structure.

To determine the Rayleigh wave phase velocity as a function of frequency, the stationary phase method was used to create synthetic waveforms (Menke 1990, section 8.5.3). The frequency content was based on the observed spectrum of the surface wave packet. Only frequencies which had significant amplitude in the observed seismogram were used to make the synthetics. A Monte Carlo method was used to select the best fitting dispersion curve from amongst 40,000 possibilities which spanned the velocity space. Preliminarily random testing of the full parameter space was invaluable in finding the best fitting dispersion relation. Mediocre fits can be achieved with several isolated dispersion functions. Each of these corresponds to synthetic siesmograms which are in-phase with the real data but shifted by one or more wavelengths. The quality of fit was determined by cross-correlating the synthetics with the Rayleigh wave seismograms from the reference timeslice (Figure 7).

Time-variance of dispersion relation

The change in phase velocity as a function of frequency can be obtained by differencing the fourier phase spectra of Rayleigh waves from one timeslice with the corresponding spectra from the reference timeslice. The spectra must be examined carefully to ensure that any phase wrapping is corrected. If the signals are of high quality however, this is straight-forward because dispersion curves are necessarily smooth functions which vary slowly with the tide. These changes in velocity derived from phase spectra were combined with the reference dispersion relation to determine how it varied during the tide cycle (Figure 8).

The 12-hour cycle of the dispersion curve is a clear statement of the power of surface wave analysis in time dependent observations. Rayleigh wave phase velocities in a layered medium are weighted averages of the velocities in each layer. High frequency waves sample only the upper most soil while low frequency waves penetrate deeply and average over a greater depth. In Figure 8, the low frequencies show the most variation. This implies that the velocity changes associated with the tidal cycle are occurring primarily in the deeper layers sampled by the longest period waves.

Shear wave structure from Rayleigh wave phase velocities

Rayleigh waves can be conceptualized as evanescent compressional and shear waves which coexist to satisfy the zero-stress conditions imposed by a free surface (Aki and Richards 1980, section 7.1). Since Rayleigh wave velocities depend on compressional and shear velocities, any modeling attempt must include both velocities for each layer. Compressional velocities were already known from the traveltime analysis and shear velocities dominate the Rayleigh wave expressions. As a result, Rayleigh waves could be used to determine shear velocities as a function of depth. The shear velocity structure was determined using a least-squares surface wave inversion algorithm based on a computer code written by Hermann (1985). Shear velocity is an insightful parameter since it is directly related to the shear strength, mu, and the bulk density, rho:

Figure 9 (left) shows a shear velocity profile which accurately recreates the dispersion relation. To minimize the risk of non-uniqueness, discrete layers were chosen so their effects on the dispersion relation were largely independent of one another (Figure 10). The extremely low velocities observed in the top two layers reflect the dry wind blown layer of sand where contact between adjacent grains is at a minimum. Since air cannot support shear stress and this matrix of sand is poorly connected, such low shear wave velocities are understandable. The shear wave velocity increases rapidly with depth in the first 1.5 meters. As pressure increases with depth, the sands are able to settle into a compact arrangement without being continually reworked by surface processes. A laboratory test revealed that the volume of the sand when loose and dry could be decreased at least 8% simply by subjecting the sample to gentle vibrations. This compaction and pressure increase the friction between grains which increases the shear strength. A slight decrease in shear wave velocity below 3.5 m is observed in our profile. We are reluctant to interpret it as significant however, because it is well within the error bounds of this method.

The drastic velocity jump across the water table, observed for the compressional waves, is not present here. Instead, shear wave velocity increases smoothly with depth. Theory (Biot 1956) and lab work (Domenico 1977) predict a decrease in shear wave velocity as saturation increases. Fluid in the pore space increases the density of the material while having little effect on the shear strength. It is important to note however, that these explanations overlook the changes in grain contacts and hydrostatic pressure which may accompany increased saturation in the field. Using in-situ VSP methods, Michaels and Barrash (1996) claim to observe an abrupt shear velocity increase of over 250% across the water table in unconsolidated sands. Though our method resolves velocities several meters under the water table, we find no evidence for the much faster saturated velocities they observe.

Determining shear velocity changes through time

The temporal variation in shear wave velocity in each of the layers is determined by considering them as small perturbations to the reference structure calculated at timeslice 8.5. This approach is justified by the small (<15%) temporal variations in the observed dispersion curves. The kernel functions that relate perturbations in the dispersion curve to perturbations in the layer velocities (Figure 10) need, therefore, be computed only once, making the inversion process very efficient.

Figure 9 (right) displays the changes in shear velocity with time. Shear velocities in the top two layers vary little during the cycle. This stability, observed for compressional waves as well, is expected because the region sits well above the water table, even at high tide. The third layer also sits above the fluctuating water table. At high tide, however, we expect layer 3 to have an increased moisture content. This layer progressively includes more of the capillary zone which precedes the advancing water table. As any sandcastle enthusiast knows, small amounts of moisture greatly increase the cohesion of sand. This cohesion increases the shear velocity of the sand matrix and could easily outweigh the density change. Further velocity modeling is needed to address the role of cohesion in partially saturated granular media.

We don't interpret layer 4, which contains the transient water table, because the effects of saturation, water table level and hydrostatic pressure cannot be separated. It is worth noting that despite the changing water level in this layer, the net velocity does not vary considerably. This is further evidence that shear wave velocities are not directly sensitive to saturation.

The bottom two layers display the greatest velocity changes. These velocities are at a minimum near high tide (hour 1 and 12.5). Both of these layers are below the water table at all times, yet their velocities are closely linked to the tide cycle. These changes cannot be explained by density or even cohesion since the sand is saturated at all times. Only the fluid pressure changes at these depths. As the water level rises, the pore pressure increases due to the weight of the additional ~1 meter of overlying water. The Darcian pressure gradient induced by the flow leads to a secondary, and less predictable head, on the order of 0.1 meter. Both of these work to reduce the effective pressure on the matrix. The resulting loss of shear strength is clearly visible in the decreased shear wave velocities.

Conclusions

We present several techniques which allow the effects of fluid saturation to be determined from small waveform variations. Our method treats each timeslice as a small deviation from a reference model. This approach is both faster and more precise than processing individual timeslices from scratch. Three procedures were presented which use this approach.

1. Waveform correlation provides a precise and repeatable method for discerning small variations in body wave traveltimes.

2. Phase spectra of surface waves can be differenced to determine small variations in phase velocity as a function of frequency.

3. Velocity structures that change with time can be accurately tracked with a linearized inverse approach.

These strategies were used to determine how partial and full fluid saturation affect seismic velocities in shallow unconsolidated sands. Several effects on shear and compressional waves were observed.

1. On the field scale, shear velocity changes are controlled predominantly by variations in the shear strength of the matrix, not the small density changes which accompany fluid saturation.

2. Shear velocities are sensitive to moisture in the unsaturated zone. At low saturation levels, moisture increases the cohesive attraction between grains, increasing the shear velocity by at least 10%.

3. In fully saturated sands, pore fluid pressure lowers the matrix shear strength, decreasing the shear velocity by at least 15%. This suggests that shear velocities could be used to monitor pressure changes in fluid reservoirs.

Acknowledgments. We would like to thank the trustees of the town of East Hampton, NY for allowing us to conduct this experiment on their quiet beaches. Josh Radoff and Kyle Anders deserve special thanks for their tireless field work and 2500+ hammer blows. The geophones used in this experiment were loaned by the USGS Water Resources Division, Branch of Geophysical Applications and Support. We also thank Ross Bagtzoglou, Roelof Versteeg, Marc Spiegelman and Robert Stoll for valuable discussions.

References

Aki, K., and P. Richards. 1980. Quantitative Seismology: Theory and Methods, Vol I, W.H. Freeman and Co., San Fransisco.

Anderson, R. A., 1998. Oil production in the 21st century: Sci. Am., v. 278, n. 3, p. 86-91.

Bachrach, R. and A. Nur. 1998. High-resolution shallow-seismic experiments in sand, part I & II, Geophys. v. 63, n. 4, p. 1225-1240.

Biot, M. A., 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range and II. High-frequency range: Acoust. Soc. Am., v.28, p.168-178.

Caress, D., Burnett, M. S., & Orcutt, J. A., 1992. Tomographic image of the axial low velocity zone at 12:50N of the East Pacific Rise, J. Geophys. Res. v. 97, 9243-9263.

Copty, N., and Y. Rubin. 1995. A stochastic approach to the characterization of lithofacies from surface seismic and well data, Water Resources Research, v.31, n.7, p. 1,673-1,686.

Dickenson, O., G. Steensma, T. Byod. 1996. Mapping of impediments to contamination flow using multicomponent reflection seismology at the Savannah River site, SAGEEP conf. proc., p. 121-127.

Doll W., R. P. Miller, and J. Xia. 1998. A noninvasive shallow seismic source comparison on the Oak Ridge Reservoir, Tennessee, Geophysics, v. 41, n. 5, p. 1318-1331.

Domenico, S. N., 1976. Effect of brine-gas mixture on velocity in an unconsolidated sand reservoir, Geophysics, v. 41, n. 5, p. 882-894.

Domenico, S. N., 1977. Elastic properties of unconsolidated porous sand reservoirs, Geophysics, v. 42, n. 7, p. 1339-1368.

Jefferson R., D. Steeples, R. Black, and T. Carr. 1998. Effects of soil-moisture content on shallow-seimic data, Geophys., v. 63, n.4, p.1357-1362.

Haeni, F.P., 1988, Application of seismic-refraction techniques to hydrologic studies: U.S. Geological Survey Techniques of Water-Resources Investigations, book 2, chap. D2, 86 p.

Hermann, R.B. 1985. Computer programs in seismology: vol. III, St. Louis University, St. Louis, MO.

Knight, R., and R. Nolan-Hoeksema. 1990. A laboratory study of the dependence of elastic wave velocities on pore scale fluid distribution. GRL, v.17, n.10, p.1529-1532.

Menke, W., and D. Abbott. 1990. Geophysical Theory, Columbia University Press, NY.

Michaels, P. and W. Barrios. 1996. The anomalous behavior of SH-waves across the water table. SAGEEP conf. proc., p. 137-145.

Robinson, M., D. Gallagher, and W. Reay. 1998. Field observations of tidal and seasonal variations in ground water discharge to tidal esturine surface water, GWMR, winter 1998, p. 83-92.

Rubin, Y., G. Mavko, and J. Harris. 1992. Mapping permeability in heterogeneous aquifers using hydrological and seismic data, Water Resources Research, v.28, n.7, p. 1,809-1,816.

Shtivelman, V., U. Frieslander, E. Zilberman and R. Amit. 1998. Mapping shallow faults at the Evrona playa site using high-resolution reflection method, Geophys., v. 63, n.4, p.1257-1264.