Tracking near-surface fluid saturation using compressional and
surface seismic waves

Michael West and William Menke

Introduction

The need for time-dependent methods
Seismic methods have been used successfully in water resource issues for decades. Refraction-style methods have been used to delineate aquifer boundaries such as bedrock (Haeni 1988). Reflection methods have been used to identify aquitards such as c lay layers (Dickenson 1996), and preferential flow pathways such as faults (Shtivelman et al. 1998). Under certain ideal conditions, the effects of fluid have been observed directly (Michaels and Barrash 1996).

Geophysical observations are a natural complement to precise but spatially-limited well measurements. Because the cost of well construction, development, and maintenance can be far greater than available resources, exploitation of continuous fields of soft data, such as seismic data, in conjunction with a realistically limited number of observation wells, is imperative. In recent years there has been increased interest within the scientific community regarding the incorporation of soft data into flow and transport simulations. Recent examples include, but are not limited to, the work of Rubin et al. (1992) who proposed a method for incorporating the exhaustive areal coverage offered by geophysical analyses and the coupling between hydraulic and seismi c properties and Copty and Rubin (1995) who made use of surface-based seismic surveys and well data for the characterization of lithofacies.

But directly observing fluids in the near surface environment can be difficult. The range of geophysical properties observed in the host material often exceeds the perturbations caused by the presence of a fluid. It is usually safe to assume however that fluid distribution is the only transient subsurface parameter on short time scales. Observing temporal changes in seismic response can provide insight into dynamic processes such as fluid flow, as well as information about static properties such as hydraulic conductivity and porosity. 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 - c ollect the "same" data sets at different times and analyze them for their differences. 4-D geophysics is poised to have a large impact on near-surface hydrogeologic applications. There are several reasons for this: fluid movement in the near surface ofte n occurs on short time scales of hours to weeks; data can be collected relatively inexpensively; and "ground truthing" is often able to provide critical feedback to experimental design.

Time-dependent seismic observation could be a powerful tool in several applications. Repeated observation could reveal how an aquifer responds to the competing influences of natural recharge and extraction by pumping. 4-D methods would allow planners t o monitor water extraction and develop ideal pumping schemes to maximize water production while minimizing the effect on a complicated aquifer system. Similar 4-D methods could provide near real time feedback in contaminated site remediation. Mapping cha nges in water distribution between wells could detect missed pockets of contamination and help quantify the effectiveness of remediation at highly heterogeneous sites.

However, there are several hurdles which currently limit the effectiveness of temporal geophysical data. The first is a poor understanding of how materials respond to saturation. The theories describing the seismic velocities of partially-saturated g ranular medium were developed 40 years ago (Biot 1956). These ideas have been tested extensively in controlled laboratory settings. Studies by Domenico (1976 and 1977) and Knight (1990) are representative of such efforts. The results vary considerably depending on the type of material, the pressure applied to the system, and the distribution of fluid in the pore spaces. 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. These experiments have been performed at megahertz frequencies on small clean samples or even glass beads. Bu t there have been few field studies to address the issue. The role of compaction, small amounts of clay or soil, and the strong relative pressure gradients in the very shallow earth are poorly understood.

A second factor limiting the use of time-dependent geophysics is the time and expense associated with data collection. Great advances have been made recently in the quality of shallow seismic data. 24-bit recording systems, fast personal computers a nd the simultaneous recording of hundreds of geophones have made high fold reflection data commonplace in 2-D and 3-D surveys (Bradford 1998, Buker et al. 1998). Most current methods rely on the principle that more data means more detail. But the weeks, and even months, required for such data collection preclude the possibility of repeating the surveys rapidly enough to observe changes.

The third impediment to time-dependent surveys is the lack of tools for analyzing geophysical data which varies with time. A true 3-D survey can use more sophisticated techniques than simply processing several parallel lines of 2-D data. Analogously, studies of transient phenomena must move beyond the time-consumming and imprecise method of processing multiple timeslices independently and then comparing them for differences. Techniques which make use of the similarity of repeated data sets could inc rease effeciency and precision.

Techniques for collecting and analyzing transient seismic data
Two sets of geophysical tools must be developed before time-variant wave-based methods (seismic methods and ground penetrating radar) can become practical. The first set of tools is needed to quantify and analyze basic differences in wave fields. Th e second is needed to quantitatively interpret such differences.

Time-dependent changes in the subsurface velocity structure manifest themselves throughout a seismic record. The effect of subsurface fluid distribution on seismic wave propagation is typically small. For this reason, underground fluid has often prov ed difficult to detect. The small scale of these effects can be advantageous in time-dependent studies however. Instead of comparing seismograms produced by fundamentally different velocity structures, the time-dependence can often be regarded as a smal l modification to a reference signal. These modifications may include shifted phases, amplitude changes and altered frequency content.

Any ground motion recorded by a geophone, Y, can be represented as a convolution of a source signal S, an Earth model M, and a receiver function R.

M is usually the term of interest, but only Y can be observed. If the source and receiver functions are adequately repeatable, then observed changes in Y can be attributed to variations in the Earth structure M. Several precautions can be taken to e nsure that changes in S and R are small. To avoid irregularities in S, a source must be able to produce the same signal every time. No source is exactly repeatable. The viability of a source must be determined by the setting and objectives of the exper iment (Doll et al.).

The effect of receiver variations can be minimized by leaving geophones in place for the duration of an experiment. This is not always feasible, but care should be taken to use the same phones in the same locations without overly disturbing the surrou nding soil. Because of manufacturing differences and years of abusive field use, it would be naive to assume that geophones of the same nominal response performance identically. If geophones are to be buried multiple times, great care is advised to ensu re phones are buried consistently at the same depths. Obviously, source and receiver locations must be closely controlled. Changes in surface soil conditions are another unfortunate variable affecting both the source and receiver (Jefferson et al. 1998) . In the field, this is likely to be beyond the control of the investigator. However, detailed records of surface conditions and any precipitation can at least aid in post facto error analysis. If these precautions are taken, then changes in recorded s ignals may safely be attributed to changes in subsurface velocity structure.

If the changes in the velocity model are relatively small (as is typical of fluid movement studies) and the region is sampled frequently, the shape of the waveforms does not change significantly from one timeslice to the next. This similarity of wavef orms is a powerful analysis tool. Rather than processing and analyzing each timeslice independently, as previous studies have done (Birkelo 1987, Bachrach and Nur 1998), great attention can be paid to how each timeslice differs from a reference data set. Each timeslice can be regarded as a small departure from the 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). The water table and soil moisture are known to vary on semi-diurnal timescales associated with lunar and solar ti de cycles, as well as less predictable oceanic weather conditions. The beach location allowed us to perform this experiment in one day.

figure 1. Location of sources, receivers and water table through time in the experiment coordinate system. Full size image / color postscript version

This eliminated errors resulting from replanting geophones and minimized the influence of changing surface conditions (Jefferson et al. 1998). Our objective was the analysis of the time-dependent velocity structure, so every precaution was taken to minim ize stochastic influences. While far from homogeneous, the beach provided some of the most uniform soil conditions outside of a laboratory. The high porosity of the sand ensured that the fluctuating tide would induce rapid and significant perturbations to the moisture content. In addition to being an ideal natural lab, the granular nature of the beach resembles many settings where time-sensitive geophysical methods may offer the most benefits. Unconfined aquifers and land fills are examples of sensit ive areas often located in unconsolidated sediments. Because of their social relevance and their susceptibility to fluid and contaminant fluxes, these settings are ideal sites for techniques such as those developed here.

A section of beach was selected west of Indian Wells in East Hampton, New York. This south-facing Long Island beach was chosen for its clean sands (devoid of larger rocks and debris), deep bedrock, wide beaches and low public use. An array of twenty- four vertical component 40 Hz Mark Products geophones on 12 cm spikes were planted with one meter between phones at a depth of 15 cm. This depth protected the phones from surface noise and wind without diminishing the quality of the surface waves or inte rfering with the water table (fig 1). The line of geophones was oriented perpendicular to the water's edge. A Geometrics Strataview-R system recorded the seismograms. The sample period of 0.0625 ms is quite short but was selected because our correlatio n 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 us to record the data without having to use pre-A/D filters. This allowed the simultaneous c ollection of high quality body and surface waves.

A sledge hammer and steel plate were chosen as the source. A modified shot gun source was ruled out because of the need for high repeatability. Numerous shots in the same hole throughout a day would be virtually impossible in loose sand. Though hamm er sources are often cited as having questionable repeatability (Keiswetter and Steeples 1994), on-site tests showed the signal to be quite repeatable using stacked records. We compared stacked seismograms of 1, 3, 10 and 30 traces. The signal to noise level increased significantly as more records were stacked. There was no deterioration of the high frequencies with increased stacking. This is clear evidence that the trigger timing of the recording system was sufficiently precise. We chose to use sta cks of 30 hammer blows for the experiment. Not only did this average any variability from the hammer source, it also minimized any contamination of the records by surf noise.

figure 2. Example seismogram from hour 8.5. A 70 ms agc window has been applied for viewing purposes. Full size image / color postscript version

The experiment began close to high tide. A set of seismograms was recorded every half hour for 12.5 hours, and ended just before the following high tide for a total of 26 timeslices. For one timeslice, at hour 8.5, a complete set of five shots was re corded (fig 1). This complete set was used to develop the reference profile. At all other times, shots A and B were recorded. These were used to track the evolution of the seismic wavefield.

Two wells were used to monitor the water table depth. One well was near each end of the array but they were offset by roughly 5 meters from the line to avoid influencing the seismic wave propagation. Well measurements began a couple hours after the seismic surveying. As expected, the well nearest to the open water fluctuated more during the tide cycle. Measurements in the inland well were of smaller amplitude and were delayed by roughly one hour.

Theodolite ranging provided accurate locations. The geophones, shot points, wells and beach topography were surveyed. After the field work, all measurements were converted to one coordinate system with the origin at shot point B (fig. 1). Range incr eases along the array toward the water. All plots use this coordinate system.

Figure 2 shows a representative vertical component seismogram from the reference timeslice at hour 8.5. The source is located at the 0 m position. No filters have been used but an automatic gain control with a window of 70 ms has been applied to allow all phases to be observed simultaneously. Direct, refracted and Raleigh waves can be observed. Without the gain, Rayleigh waves dominate the record. A small air wave is present but requires extreme high pass filtering to observe.

Hydrologic model of the water table

The interpretation of this experiment depends heavily on the depth of the water table through time. From the measurements, well 2 appears to track well 1 closely with decreased amplitude and a delay of 1-2 hour. From these observations one would assu me the water table depths along the array are intermediate between the levels in the two wells. It is conceivable however that the delay between the wells exceeded one tide cycle. If this were true, several oscillations of the water table might exist al ong the array.

figure 3. Circles mark observed data from wells. The gray line indicates the water table in the diffusion model. color postscript version

A model was created to help understand the dynamics of the water table. If constant diffusivity is imposed, then the system can be modeled as a diffusion of hydraulic head (water table depth) away from the tidal forcing of the open water. Water table depths from well 1 were used as a boundary condition on one end of the system. A finite-difference approach was used to estimate water table depths across the rest of the array. The model was run using a range of diffusivities until the response matched the observed water table in well 2. Using only one free parameter, two conditions had to be matched - the amplitude and time lag of the signal in well 2. Both conditions could be met because we were modeling a real system (fig 3). The existence of a fit implicitly confirms our constant diffusivity assumption.

The model demonstrated that the sand environment over our distance range, responds within a couple of hours to changes in the tide level. This indicates that the water level in well 2 is responding to the current tide and has little dependence on the previous cycle. More importantly, there should be no additional oscillations of the water table between the wells which could affect the seismic velocity structure in an unforeseen way.

Quantifying the time dependence of data

Figure 4 illustrates the time dependent nature of the waveforms in this study. The upper left is a record section sampled at a relatively low tide at hour 7.5. In the lower right is the identical source-receiver pattern at hour 12.5, high tide. Thou gh the waves are quite similar, they are systematically shifted and stretched. These differences can best be observed by cross correlating the waveforms. The bottom plot of 4 shows an example of such a correlation. To obtain the relative shift between two seismograms as a function of traveltime, a narrow time window of 10 ms was first selected. The signals outside of this window were tapered to zero. The correlation of the two waveforms was determined using a fourier-domain cross-correlation techniqu e. The time lag which yielded the highest correlation value is indicative of the time difference between the two waves at that traveltime. Figure 4 was produced by performing this process for all time windows on seismograms at all offsets. The correla tion plot is a useful way to observe the gross time-dependence of a signal. The plot reveals portions of the record

figure 4. Record sections (top) from timeslice 7.5 and 12.5 are compared for shifted phases (bottom). Full size image of seismic records and correlation plot / color postscript version

where the signal is coherent from one timeslice to another. It also displays the effects of the time-variant velocity structure on the seismograms. Where portions of the signal are dominated by noise, this technique generates arbitrary correlations. Th ese regions are easily identified by their extreme oscillatory correlations. The first 25 ms and the lower right half of the record display this type of behavior.

Portions of the seismogram which are shifted by more than half a wavelength may be miscorrelated due to the extreme phase shift. A correlation scheme which allows user interaction could conceivably compensate for this problem. A simpler solution, use d here, is to ensure that the system is sampled adequately in time. Sufficiently high "temporal sampling" is essential to the techniques presented here since they rely on waveform similarities. The threat of temporal aliasing should be given the same at tention in project design as the more familiar time-series and spatial aliasing.

The direct and refracted compressional waves are clearly visible on this plot. Compressional wave time shifts between the low and high tide signals are on the order of 1 millisecond. The largest timeshifts are clearly associated with the surface wave s. Delays of up to 4 ms can be observed. Based on this diagram, it is clear that both the compressional body waves and the Rayleigh surface waves contain valuable information about the time dependence of the velocity structure. Since each offers a diff erent type of information, they are both investigated here.

Compressional wave methods

Variations in body wave traveltimes are one of the more readily measurable effects of subsurface changes. Any alterations to the velocities or the raypaths will affect these traveltimes. For this reason time-dependent studies are well suited to body wave analysis. Our objective is to pursue methods which take advantage of the similarity of seismograms through time. The first step is to create a velocity model which describes the subsurface at one time in detail. Then, time differences such as thos e observed in figure 3 can be used to detect how this reference model was altered through time.

Characterizing the subsurface at one reference time
The procedure used in this step is typical of many other refraction-style surveys. A set of five shots was collected at hour 8.5, when the water table was near its lowest level - one shot from each end of the

figure 5. Traveltime curves for the first arrival compressional waves from five sources. Data was recorded at hour 8.5, the reference timeslice. Full size image / color postscri pt version

array, a shot at the midpoint and a shot 10 meters beyond each end of the array. First arrival traveltimes from these records were hand picked (fig 5). A 200-1000 Hz bandpass filter was applied to enhance the arrivals. However, times picked from filter ed data were always checked against the raw records to ensure that filtering did not distort the traveltime picks. A preliminary 2-D velocity model was constructed to match the observed traveltime data using the RTMOD forward program (Caress 1992) and tr ial-and-error fitting (fig 6). As the slope of the traveltime curves indicate, the model is dominated by a two layer velocity structure. Interpreting traveltimes to obtain velocity structure is a non-unique problem. To keep our model as simple as possi ble and to avoid interpreting poorly resolved structures, we did not allow velocities to vary laterally within each layer.

The bottom layer velocity is well constrained by the slope of the refracted ray traveltimes. The top layer is not as well constrained however by first arrivals. At near offsets, the surface waves and direct

figure 6. Rays showing paths of the direct and refracted compressional waves for a source at the origin (source B). A slight gradient has been used in each of the layers to facilitate turning. Full size image / color postscript version

compressional waves arrive together. The first arrival picks are unreliable. However, the drastic contrast in velocities between the two layers presents an unusual opportunity to measure the uppermost velocity. Rays refracted along the water table inte rface must leave the source at the critical angle. To be refracted in this system, rays must depart the source between 5 and 10 degrees from vertical. The traveltime through a layer at these small angles differ from vertical traveltimes by only 0.4 to 1.5 %. In this situation, the traveltimes through the top layer can be approximated using vertical raypaths. The refracted ray traveltimes (f ig 5) display trends which mimic the topography of the site. Though each geophone lies 15 cm below the surface, in an absolute reference frame they are at different elevations. If topography is the primary cause of the variations in the refracted ray tr aveltimes, we can correct the traveltimes for the effect of topography across the array. The top layer velocity can be estimated by finding the velocity which best corrects the refracted ray traveltimes to a straight line. This procedure was applied to the data. A velocity of 145+-8 m/s was found to do the best job of bringing the data into a

figure 7. Plot showing the near-surface velocity which best relocates the geophones at each timeslice. Full size image / color postscript version

straight line. Applying this technique to each timeslice independently yielded nearly the same velocity each time, reflecting the robustness of the procedure (fig. 7). The small variations through time likely describe the error of the method; we do not interpret changes in soil properties from this signal. The depth over which this velocity is accurate cannot be ascertained from this method. The geophone elevations varied by 40 cm so this velocity can best be thought of as an average over the top 40 c m. After the velocity of 145 m/s was determined, all traveltimes were corrected by their elevation divided by this velocity to cancel the effects of variable topography.

Determining compressional velocity changes through time
Once the reference model had been determined, the time dependence of the traveltimes was examined.

figure 8. Steps used to prepare seismograms for correlation. top: raw data, middle: bandpass filtered between 200 and 1000 Hz, bottom: taper applied to isolate the signal of interest. Scale is in seconds. Full size image / color postscript version

The difference in traveltimes between each timeslice and the reference was calculated using the fourier-domain cross-correlation method. The data was first band passed filtered between 200 and 1000 Hz to highlight the compressional arrival. For each tra ce in the reference seismograms, a 5 ms window around the compressional first arrival was defined. A 3 ms taper was applied to each side of the window to isolate the signal (fig 8). After applying the identical windowing to each timeslice, the record wa s correlated with the reference seismogram to determine the time lag which resulted in the best correlation. This lag was then added to the reference to obtain the new traveltime. Figure 9 shows a plot of the time-dependence of the refracted compression al waves from source A. The mean traveltime at each distance has been removed to highlight the variance. Only those stations for which the refracted wave is the first arrival are shown. A similar plot for data from source B shares the same trend.

figure 9. Each sheet is a normalized traveltime curve for one timeslice. A mean arrival time has been removed from each station to highlight the time-dependence. color postscript version / c olor postscript version with inset

The traveltimes vary by up to 1 ms over the tidal cycle. It is notable that most of the variation is a uniform offset across the array and not a change in the slope of the refracted rays. This indicates that most of the time-dependence lies above th e water table. Surprisingly, high tide is associated with longer traveltimes. This is counter-intuitive if one only considers the water table. Bachrach et al. (1998) observed the same phenomenon in beach sands. They present a possible qualitative expl anation relying on Biot's theory. Fully saturated sand, they claims, would have a bulk modulus of the sand matrix combined with the bulk modulus of the water. The large increase in

figure 10. High and low tide compressional wave velocity models. See figure 15 for specific velocities. Full size image / color postscript version

bulk modulus outweighs the increase in density associated with higher saturation. At partial saturation however, the pore fluid exists in pockets of moisture clinging to individual grains in the interstitial spaces. Devoid of interconnectedness, the flu id contributes little to the bulk modulus of the medium. However, the density still works to decrease the velocity. The result is a decrease in velocity when the sand is partially saturated.

To loosely quantify this effect in our signals, we compared models from low and high tide seismograms. Since the time-variance was similar from either shot point, we chose to simplify the interpretation by approximating the structure as one-dimensiona l (fig 10). The upper most velocity is obtained from the geophone relocation procedure. The bottom layer velocity is calculated from the slope of the refracted ray traveltimes. The middle layer velocities and thickness are poorly constrained by our dat a however. The real velocity of this material is likely a gradient reflecting the influences of compaction, pressure and saturation. We chose to solve for the _average_ velocity of this layer by fixing the high and low tide water levels based on the wel l data and then solving for the corresponding velocities. Figure 10 shows a drop in the middle layer velocity from 325 m/s to 205 m/s when the tide increases. A higher water table would influence the velocity in this layer through two mechanisms. The s aturation increase associated with the rising capillary zone above the water table would lower the average velocity. Additionally, the shallower average depth of this layer would preferentially sample the slower less compacted sands. The decrease in hal fspace velocity reflects a small change in the slope of the refracted ray traveltimes. The decrease in velocity of the bottom layer from 1650 to 1600 m/s may be a response to the slight change in dip between the two wells. Alternatively, it may be a man ifestation of compaction. The high tide bottom layer encompasses shallower less compacted material which is subjected to less overlying pressure, potentially decreasing the bulk density. The result is a velocity decrease from 1650 m/s to 1600m/s.

Surface wave analysis

Rayleigh waves, sometimes referred to as "ground roll", dominate the seismograms in this experiment (fig. 2 and 4). 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. While the surface waves at different timeslices appear related, the waveforms are not simply time shifted. Each frequency responds differently to changes in the velocity structure. It is tempting to take the time delays observed in figure 4 and apply them directly to determine how the surface wave velocities change through time. A clear distinction must be made however between the velocity at which each frequency propagates (phase velocity) and the velocity at wh ich the wave packet propagates (group velocity). The cross-correlation procedure in figure 4 measures the shift of a given peak in the seismogram. Therefore it is sensitive to variations in group velocity. It is the phase velocities however which can b e related to subsurface structure.

Rayleigh wave modeling to determine dispersion relation
To determine the Rayleigh wave phase velocity as a function of frequency, a 1-D modeling effort was undertaken. The stationary phase method was used to approximate the waveforms (Menke 1990). While

figure 11. Top: Observed and synthetic Rayleigh waves. Observed data has been tapered to zero before the Rayleigh arrival for clarity. Bottom: Dispersion curve which best fits the observed seismograms. Seismograms: Full size image / color postscript version Disp. curve: Full size image / color postscript version

no setting can truly be represented by a 1-D model, the P-wave imaging and a priori well data support the assumption that the system is dominated by vertical gradients. The stationary phase method assumes that the seismogram is dominated at any instant b y one frequency. The amplitude of the seismogram at this point depends on the distance from the source, the strength and phase of this stationary frequency and how rapidly the group velocity is changing. The frequency content was based on the spectrum o f the surface wave packet. Only frequencies which had a measurable presence in the real seismogram were used in the creation of synthetics. The relationship between velocity and frequency is often referred to as a dispersion curve because it governs how the wave train disperses with time. A forward approach was used to test dispersion curves. A second order function of frequency was chosen to represent the dispersion relation. The parameters of this expression were then varied over all conceivable va lues. 40,000 second order functions, spanning all possible velocities, were tested. The synthetic seismograms generated from each dispersion relation were cross-correlated with the Rayleigh wave seismograms from the reference timeslice to determine the quality of the fit (fig 11). In an approach such as this, it is essential to test numerous unrelated dispersion relations. Mediocre fits can be achieved with several vary different Rayleigh wave velocity functions. Each of these correspond to synthetic siesmograms which are roughly in-phase with the real data but are shifted by one or more wavelengths. Randomly testing the full parameter space ensures that the best fitting dispersion relation is found.

Time-variance of dispersion relation
The change in phase velocity as a function of frequency can be obtained in a straight forward manner. By differencing the fourier phase spectra of two similar records, the change in phase velocity as a function of frequency can be obtained directly. The spectra must be examined carefully to ensure that any phase wrapping is corrected. If the signals are of high quality however, this should be straight forward. The lowest frequencies are less likely to be wrapped because of their longer wavelength . Since the dispersion relations are smooth, the phase shift function is continuous and any wrapping should be evident. The changes in velocity derived from phase spectra were combined with the previously determined dispersion relation from the referenc e time. This allowed us to determine the dispersion relation as a function of time.

figure 12. Each sheet displays the change in velocity as a function of frequency for one timeslice Full size image / color postscript version

The 12-hour cycle of the dispersion curve (figure 12) is a clear statement of the power of surface wave analysis in time dependent observations. Rayleigh wave 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 12, the low frequencies show the most variation. This implies that the velocity changes associated with the tidal cycl e are occurring primarily in the deeper layers sampled by the long period waves.

Shear wave structure from Rayleigh wave phase velocities
Rayleigh waves can be conceptualized as evanescent compressional and shear waves which must coexist to satisfy the zero-stress conditions imposed by a free surface (Aki and Richards 1980). Since Rayleigh wave velocities depend on compressional and sh ear velocities any modeling attempt must include both velocities for each layer. Shear velocities dominate these expressions. As a result, Rayleigh waves are an excellent indicator of shear velocities, and hence, shear modulus of a material.

This sensitivity of Rayleigh waves to shear strength means changes in the dispersion relation can be used to calculate changes in shear velocity as a function of depth. The shear velocity structure was determined with the help of a standard surface wa ve forward modeling program (Hermann 1985). This program allows the user to predict the dispersion relation for a given layered velocity structure. Compressional velocities were constrained by the traveltime analysis so the shear velocity was the domina nt free variable. Figure 13 shows such a shear velocity profile which accurately recreates the dispersion relation. To minimize the risk of non-uniqueness, layers were chosen so that their effects on the dispersion relation were largely independent of o ne another (fig. 14). The drastic velocity jump across the water table which was observed for the compressional waves, is not present here. Instead, shear wave velocity increases smoothly with depth. Note the extremely low velocities observed in the to p two layers. This is the dry wind blown layer of sand which is very loosely packed. In this top layer, contact between adjacent grains of sand is at a minimum. Since air cannot support shear stress and this

figure 13. Left: Shear wave velocity structure. Right: Changes in this structure over the tide cycle. Full size image / color postscript version

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 depth increases, the sand matrix is likely to be more compacted. These sands are abl e to settle into a dense arrangement without being constantly reworked by surface processes. A laboratory test revealed that the volume of our sand when loose and dry could be decreased at least 8% simply by subjecting the sample to gentle vibrations. T his compaction would result in a large increase in grain contact, and thus shear strength. The pressure induced by overlying sand would further increase the shear strength of the sand matrix. Not only does pressure tend to enhance compaction, it also in creases the friction between adjacent grains. The shear strength increases due to compaction and pressure overwhelm any small changes in the density related to the grain packing.

The concept that compaction and pressure dominate the determination of shear velocity is further supported by the absence of a jump in velocity at depths associated with the water table. In fact, though the water table at the reference time lies at ro ughly 1.9 m, there are only small changes in shear velocity below 1.5 meters. This is reasonable since fluids do not, in theory, contribute to shear strength. Michaels and Barrash (1996) claim to observe an abrupt velocity increase of over 250% across the water table in unconsolidated sands. Their "dry" sand velocity of 161 m/s is essentially the same as ours. However we find no evidence for the 422 m/s "saturated" velocity they observe, though our method can resolve velocities several meters under t he water table. Biot predicts decreased shear wave velocities as saturation increases. Fluid in the pore spaces increases the bulk density of the material while having little effect on the shear strength. It is important to note however, that this theo ry overlooks the possibility of changing grain contacts which may accompany increased saturation.

Unlike the compressional wave data, the surface wave approach has considerable resolution below the water table interface. A slight decrease in shear wave velocity below 3.5 m is observed in our profile. We are reluctant to interpret it as significan t however, because it is well within the error bounds of this method.

Determining shear velocity changes through time
Our goal is to examine changes through time without processing each timeslice individually. A thorough understanding of one snapshot in time is a powerful tool for understanding the time-dependent behavior of the system. In this case, we have alread y developed a shear wave velocity structure which describes the subsurface at the reference time. The variations of the surface wave phase velocities are small. They can be thought of as perturbations to the reference timeslice. A finite difference lin ear inversion was constructed to determine the changes in the shear velocity model that are needed to match the dispersion relation through time (Menke 1989). The kernel functions for this inversion are obtained by perturbing the model one layer at a tim e (Fig 14). These functions relate small changes in subsurface velocity structure to changes in the predicted dispersion relation. They are valid for any model which differs only slightly from the reference velocity profile. Stated formally, the forwa rd and inverse problems are

dVi=Gij*dBj

dB=[GTG]-1 * G * dV

where dV is the difference between the Rayleigh wave phase velocities at the time of interest and the

figure 14. Kernel functions for inversion. Each function shows the degree to which frequencies are controlled by a given layer. Full size image / color postscript version

reference time; dB is the difference between the perturbed shear velocity structure and the reference velocity structure; and G is the matrix of dV/dB obtained by perturbing the reference velocity structure one layer at a time. If the layers are chosen s o that each is largely responsible for a different frequency range in the dispersion relation, the inversion can be performed without any damping. In the example above, the dispersion curves predicted by this method were nearly indistinguishable from the data.

Figure 13 displays the changes in shear velocity with time. Shear velocities in the top two layers vary little during the cycle. The same feature was found in the P wave velocities. This region sat well above the water table, even at high tide. The third layer also sits above the fluctuating water table. However, capillary forces draw moisture into the region. The moisture content decreases away from the water table. We expect layer 3 to have a higher moisture content at high tide when it is clos er to the water table. As any sandcastle enthusiast knows, small amounts of moisture greatly increase the cohesion of sand. This cohesion increases the shear strength of the sand matrix and could easily outweigh the tiny effect of increased density. Th e result is an increase in shear wave velocity. This interpretation is speculative. However, we consider the absence of cohesion in prior theoretical calculations to be a serious shortcoming.

Layer four should be interpreted cautiously. At high tides, the water table is in this layer, while at low tides the region is entirely above the water table. The layer is a composite of two very different regimes - one is fully saturated, the other is partially saturated. Since the water table boundary is not stationary, it is impossible to separate the effects of the two regimes. It is worth noting that despite the changing water level within this layer the velocity does not display great variat ions. This is further evidence that shear wave velocities are not directly sensitive to saturation.

The bottom two layers display the strongest time dependence. When the tides are high (hour 1 and 12.5), the velocities decrease. 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 changes since the sand is saturated at all times. The only aspect of the subsurface which does change at these depths is the hydrostatic pressure. As the water level rises, the hydrostatic, o r pore fluid, pressure increases. Though the pore fluid pressure increases, the pressure exerted on the sand matrix by the overlying material remains the same. This reduces the effective pressure on the matrix. The result is a decrease in the shear str ength of the system.

From these observations, we propose that shear wave velocities in unconsolidated granular media are governed primarily by changes in grain contacts. Compaction increases the shear strength by increasing the number of grain contacts; pressure increases the shear strength by increasing the frictional force across each grain contact; at low saturations, moisture increases the shear strength by increasing the cohesive attraction between grains; in fully saturated sands, pore fluid pressure lowers the effe ctive pressure and results in decreased shear strength. These effects can change the shear wave velocity by orders of magnitude. On the field scale, these considerations are far more significant than the small changes in density associated with moisture content.

Comparison of compressional and shear wave velocity changes
Based on the compressional and shear velocity structures presented above, several material parameters can be estimated. It should be noted that the layered model approach averages velocities over a range of depths. The precise value of the velocitie s may differ considerably from the averages stated in figure 15. However, the time-variance of these parameters through the tidal cycle is well constrained using the techniques presented above.

Full size image / color postscript version

Full size image / color postscript version

. depth, m 0 0.15 1.2 > 2.0

LOW TIDE compressional velocity,m/s shear velocity, m/s Vp/Vs ratio Poisson's ratio 145 11.0 13.1 0.497 145 40.5 3.6 0.458 325 115 2.8 0.429 1646 174 9.5 0.494

HIGH TIDE compressional velocity,m/s shear velocity, m/s Vp/Vs ratio Poisson's ratio 145 11.2 12.9 0.497 145 38.3 3.8 0.463 205 114 1.8 0.278 1600 161 9.9 0.495

figure 15. Top: Compressional and shear velocities at high and low tide. Bottom: Material parameters calculated at selected depths.

Conclusions

Fluid distribution in the subsurface manifests itself through measureable, though complex, velocity changes. Such velocity changes are evident throughout a seismic record. Any seismic method of monitoring fluids in the subsurface must rely on disting uishing these subtle waveform changes. Often, these changes are too small to be accurately quantified by eye. Our approaches capitalize on the similarity of shifting waveforms. Signals which are nearly the same from one timeslice to the next can be com pared far more accurately than signals with wildly varying shapes. We discuss three techniques which make use of this similarity:

1. Waveform correlation provides a precise and repeatable method for discerning small variations in body wave traveltimes. When applied to surface waves, this technique can be used to calculate changes in group velocity.

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

3. Velocity structures which change with time can be accurately tracked with a finite difference inverse approach, provided the variations from a reference model are sufficiently small.

The ability to interprete small waveform changes has several implications for tracking transient features such fluid in the shallow subsurface. The most important is the need to consider the "temporal sampling" of a signal. To ensure signals c an be accurately compared, seismic surveys of a site must be repeated often enough to avoid aliasing between timeslices. To track complicated systems, significant attention will have to be given to making surveys repeatable on short timescales. This con cern will have to be balanced carefully with the ubiquitous necessity for dense spatial sampling.

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 geopho nes 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.