New Insights from Legacy Seismic Data regarding Basalt Elevations and Variability on the Hanford Site

Surface-based seismic methods can provide useful information about subsurface structure as it relates to key hydrologic parameters (e.g., Hubbard and Linde, 2010). Seismic reflection methods utilize waves that reflect off subsurface boundaries to produce images of subsurface stratigraphy and structure. Seismic refraction methods utilize waves that are transmitted through the subsurface, turning back toward the surface as they encounter faster-velocity materials at depth (Steeples, 2005). Travel times of refracted waves can be used to image the seismic velocity structure within the subsurface. It is most common to measure compressional wave velocity (Vp) with refraction methods. Shear wave velocity (Vs) can be measured with dispersive Rayleigh waves (Park et al., 1999). Seismic velocities vary systematically with lithology, porosity, saturation, and pressure (e.g., Mavko et al., 2009); thus, they are useful for identifying bedrock depths, the water table, and zones that are potentially more permeable.

Several seismic reflection surveys have been conducted in the Gable Gap area of the Hanford site. During 1979 and 1980, approximately 80 km of seismic reflection data were acquired as part of the Basalt Waste Isolation Project to identify deep storage targets for spent fuel (SSC, 1979,1980). More recently, several shallow seismic investigations have been carried out to identify basalt elevations and structure for hydrologic investigations (Cummins, 2009; Hyde et al., 2011) Hyde et al. (2011) showed the utility of the seismic landstreamer for the rapid collection of seismic data on the Hanford site. Their results showed a rugged and faulted basalt surface along 12 km of seismic profiles. A subsequent landstreamer study (Sunwall et al., 2011) added 12 km of seismic profiles within the Gable Gap area and further supported a rugged and faulted basalt surface. A 2012 study that integrated all of the existing seismic reflection data in the area and compared known basalt depths to the seismic interpretations concluded that the reflection data often overestimated basalt depths, possibly due to an inaccurate seismic velocity model used to convert travel times to depth (Williams et al., 2012).

In this study, 12 km of legacy landstreamer data (Sunwall et al., 2011) are analyzed to measure Vp and Vs from the refracted p-waves and the dispersive Rayleigh wave signals. The resulting velocity models provide an estimate of basalt depth that is independent from the pre-stack depth migrated (PSDM)–derived estimates. Comparing these new results to borehole observations and integrating them with the reflection images provides an updated interpretation of depth to basalt within the Gable Gap and 200 E areas that will be used to refine the site geologic model. Zones of interest are also identified where additional characterization is needed to evaluate the possibility of vertical communication between the unconfined aquifer and deeper aquifer systems.

The Gable Gap (Figure 1) lies within the Yakima fold and thrust belt. To the north of the study area, the asymmetric east-to-west Gable Mountain anticline (GMA) folds the CRBs and the interbedded Ellensburg formation. The northern dipping limb of this anticline is steep, whereas the southern limb has a relatively gentle dip. Second-order syncline and anticline pairs have been inferred to exist along the southern limb and throughout the study area (e.g., Ault, 1981). Previous studies have hypothesized that secondary normal and thrust faults may also exist along the southern limb of the GMA (e.g., Hyde et al., 2011).

The principal CRB unit of hydrologic interest in the Gable Gap area is the Elephant Mountain member, which is considered the base of the unconfined, suprabasalt aquifer (Bjornstad et al., 2010). However, erosion related to Pleistocene flooding of Lake Missoula has locally eroded the Elephant Mountain member, exposing older basalt units and sediment interbeds to the unconfined aquifer, and may allow vertical communication between the unconfined and confined aquifers in these areas (Graham et al., 1984).

Above the CRB, the unconfined suprabasalt aquifer is hosted in the flood-deposited Hanford formation. The Hanford formation contains both gravel- (H1 and H3) and sand- (H2) dominated units, and their distribution reflects flood dynamics. The older Ringold and Cold Creek units are river deposits containing gravels, sands, and silts. However, within the Gable Gap area, they were largely eroded during the Lake Missoula floods and are unlikely targets for refraction imaging.

The Rohay and Brouns measurements indicate that Vp should be a reliable means of distinguishing basalt from the Hanford formation. Unsaturated Cold Creek is also easily distinguished from basalt. Vp in saturated Cold Creek and the Ringold formation overlap with the basalt measurements, indicating that it will not be possible to differentiate them using Vp measurements. Additionally, Ringold Vs overlaps with basalt Vs, suggesting that even unsaturated Ringold cannot be confidently distinguished from basalt on the basis of seismic velocity alone.

The landstreamer is a rapid seismic acquisition system that consists of a string of geophones and an accelerated weight drop (AWD) towed behind a vehicle. The geophone spacing is fixed so that after each shot, the entire system is moved forward, and the source-receiver offsets remain constant. The streamer survey geometry is suitable for reflection imaging, compressional wave refraction imaging, and Rayleigh wave dispersion analysis.

Approximately 12 km of landstreamer data were collected in April and May 2011 by the Confederated Tribes of the Umatilla Indian Reservation and Montana Technological University. The streamer comprised 96 gimballed 30-Hz geophones spaced 2 m apart. The streamer was towed by a pickup truck, which also towed an AWD 6 m in front of the streamer. The AWD was a 227-kg steel ram lifted by a hydraulic pump and accelerated by an elastic band. The AWD was used to vertically strike a 60 × 60-cm steel plate at 2-m intervals along each profile. This source produces high-amplitude seismic energy in the 5- to 150-Hz band and facilitates rapid data acquisition. The data set includes 5,564 shot gathers, and each shot record contains 2 seconds of data sampled every 0.5 ms. The minimum offset was 6 m, and the maximum offset was 196 m.

Sunwall et al. (2011) processed the reflection data and produced PSDM images. However, the raw data (Figure 3) show clear direct and refracted arrivals with apparent Vp that is consistent with suprabasalt sediments (∼1,200 ms⁻¹) and the underlying basalt (∼4,500 ms⁻¹). Dispersive surface waves are also evident and suggest suprabasalt Vs between 200 and 600 ms⁻¹. In this article, the first-arrival travel times and Rayleigh wave dispersion curves are used to image Vp down to basalt and Vs in the upper 15–30 m, respectively.

Seismic refraction methods utilize the travel time of the first arriving body wave for each source–receiver pair. Seismic velocities generally increase with depth, causing seismic rays to turn, or refract, back toward the surface as they propagate. At short source–receiver distances, the first arrival often represents the wave traveling directly from the source to the receiver, and the slope of the distance–time curve is the inverse of seismic velocity. At longer distances, it is often observed that the travel-time curve becomes flatter, indicating that the waves have traveled through a higher-velocity medium. Seismic velocities can be measured directly from the slopes of the travel-time curves. However, when there is lateral variability in the subsurface, a more robust approach is required. Tomography is an iterative method where the subsurface is represented by many small elements of constant velocity. An initial velocity structure is chosen (usually an increase in velocity with depth), and travel times for the initial model are predicted. The difference between the predicted and the observed arrival times is used to find a model update that reduces the misfit. The process of predicting travel times and updating the model is iterated until the data misfit becomes acceptable.

Finding an appropriate update to the model requires solving an inverse problem that is generally poorly constrained. The data contain errors, some regions of the model may not contribute to the travel-time prediction, and solutions are non-unique (e.g., de Wit et al., 2012). The problem can be made stable by including smoothness constraints on the model, known as regularization. Regularization constraints place a penalty on model parameters that are very different from their neighbors. There is a trade-off between data misfit and model smoothness, and it is possible to find many velocity models that explain the data.

The conventional approach to regularizing the seismic tomography problem is to minimize the L2 norm (square root of the summed values) of the velocity gradient. This approach will not allow sharp boundaries, and because the L2 norm seeks to normally distribute velocity gradients throughout the model, it tends to produce smooth gradients in areas where the data do not constrain the model. An alternative approach is to minimize the L1 norm of the velocity gradient (sum of the absolute values), which allows sharp boundaries to develop and suppresses gradients in areas where the data do not constrain the model. In the Gable Gap area, the transition from suprabasalt sediments to basalt is likely a sharp velocity contrast for several reasons. First, borehole geophysical measurements indicate that the transition from suprabasalt sediments to basalt is sharp (Rohay and Brouns, 2007; Hyde et al., 2011; Figure 2). The raw seismic shots also contain several indicators of a sharp seismic boundary; the apparent basalt Vp measured on the raw data is ∼4,500 m s⁻¹, and the apparent sediment velocity is ∼1,200 m s⁻¹ (Figure 3). Much of the data also contain a converted shear wave phase with an apparent velocity of ∼2,200 m s⁻¹ (Figure 3). The presence of this phase indicates that Vs in the basalt is greater than the Vp of the material immediately above it (St. Clair and Liberty, 2019).

A two-dimensional tomography code (St. Clair, 2015) written in MATLAB was used to invert manually interpreted first-arrival observations. The code predicts travel times using the shortest path raytracing method (Moser, 1991), and the inverse problem is regularized with first-order derivative operators in the vertical and horizontal directions. The L1 constraint on model gradient is implemented with an iterative least squares algorithm (Ajo-Franklin et al., 2007). The high spatial density of the data results in many redundant raypaths; thus, travel times were interpreted on every second to third shot gather.

Travel-time tomography, like many geophysical methods, suffers from uncertainty and non-uniqueness. Uncertainty is due to noise in the data, the possible correlation of different model parameters, and non-uniform sensitivity of model parameters to the data. The same factors contribute to non-uniqueness of the solution. Because these issues are inherent in the method, the approach used to model the data needs to be consistent with the geologic setting, using information not contained in the travel times (e.g., the presence of converted phases) as a constraint for model selection and validation.

Non-uniqueness in the tomography solution occurs both because of errors and uncertainty in the input data and because many different velocity structures can have identical travel-time curves (e.g., Shearer, 1999). The choice of an L1 constraint on model smoothness favors models with a sharp boundary. The range of acceptable models is further limited by preferring solutions with the minimum amount of structure required to fit the data within estimated uncertainty. The picking error is estimated to be on the order of 1–3 ms; thus, the root mean square errors for the models should be ∼3 ms. Finally, regions of the models are masked if they do not contribute to the data misfit (i.e., where no rays are present) using the derivative weight sum. The derivative weight sum represents the total length of raypaths that pass through a given model cell.

The data set contains several points where profiles intersect (Figure 1b). Comparing the velocity estimates and interpreted depths to basalt at these points can give some insight into the precision of the approach. Since each profile was inverted independently, the coincident measurements can be considered as independent observations. Finally, models are validated by comparing the predicted depth to basalt to nearby well observations and, where available, to velocity profiles derived from check shot data in nearby wells.

Rayleigh waves are surface waves that have both vertical and horizontal components of motion. They travel at phase velocities, which are slightly slower than shear waves, and the lower-frequency (longer-wavelength) components are sensitive to greater depths. A typical Rayleigh wave dispersion curve will have higher phase velocities at low frequencies and lower phase velocities at high frequencies.

Like travel-time tomography, MASW requires an inverse solution to find the optimal Vs-depth profile that agrees with the measured data. The inversion is iterative and requires a method for predicting the dispersion curve for any given model. The propagator matrix approach described in Aki and Richards (2002) was used to predict dispersion curves, and first-order difference operators were used to constrain the inversions to be smooth. The Rayleigh wave phase velocity is sensitive to Vs, Vp, and density. Since the influence of Vp and density is small compared to Vs (Xia et al., 1999), the Vp/Vs ratio was fixed at 2, with density increasing linearly from 1,900 to 2,100 m kg⁻³ from the surface to a depth of 35 m. Models are parameterized as layers that increase in thickness from 1 to 3 m with increasing depth.

While the 30-Hz geophones used in this survey limit the amount of low-frequency content required for deeper MASW sensitivity, Ivanov (2008) demonstrated that dispersion curves can be interpreted down to about an octave lower than the natural frequency of the geophones. Frequency-phase velocity images for the data presented here show interpretable dispersion down to ∼15–20 Hz (Figure 3). To indicate where the models are sensitive to the data, the partial derivatives of data misfit with respect to model parameters are summed for each model layer, and layers that have little influence on data fit are masked. This procedure suggests that the dispersion data are most sensitive to the upper 10–20 m. Basalt depths are 50–100 m, so the data do not constrain basalt properties; however, they do provide information about shallow sediment properties.

The MASW approach produces one-dimensional models that smear geologic structure over the width of the geophone aperture of 196 m. Thus, the resulting images are unlikely to capture strong lateral changes in Vs but will highlight long-wavelength structure.

Line A North (Figure 5) is an approximately 4.1-km, north-to-south profile (see Figure 1). It is the longest profile in the data set and crosses two troughs interpreted as paleo-channels (distances of 400–800 m and 2,000–3,000 m in Figure 5a). Vp in the suprabasalt sediment section ranges from 400 and 1,500 m s⁻¹, and basalt Vp ranges between 3,000 and 4,500 m s⁻¹. The 3,000-m s⁻¹ Vp contour is intermediate between sediment and basalt velocities and typically lies at the center of the steepest vertical velocity gradient; thus, this velocity contour was selected to interpret the transition from sediment to basalt.

Lower basalt Vp is apparent along the edges of regions where the depth to basalt is not well constrained (distances of 400–800 m and 2,000–3,000 m in Figure 5a). Here, the data lack the coverage necessary to adequately constrain basalt properties. In contrast, the lower basalt velocity near 3,400–3,500 m is well constrained and nicely correlates with a discontinuity in the PSDM reflection image (yellow colors representing Vp ∼2,500 m/s; Figure 5a). This area may represent a locally permeable basalt feature. There are two regions along A North where the refraction data do not constrain depth to basalt between profile distances of 300 and 1,000 m and 2,000 and 3,000 m. Here, either the basalt may be too deep to image with the streamer offset or the basalt velocity may be lower compared to other regions. It is also noteworthy that in these areas, the PSDM image shows a less continuous reflector, suggesting that the basalt may be fractured or otherwise damaged through erosion or faulting, supporting the interpretation that basalt Vp is lower.

There are three locations along this profile that are intersected by other lines, and the depth to the 3,000-m s⁻¹ contour along those profiles is shown by red dots. Where lines B and D cross, interpreted depth to basalt on all three profiles is very close. The 3,000-m s⁻¹ contours are within 2 m at the Line D crossing and 4 m at the Line B crossing. Line C West has the depth to 3,000-m s⁻¹ at approximately the same depth as the reflections in the PSDM image.

Three well observations of basalt elevation are within 500 m of Line A North to compare to the refraction result in Figure 5a. In two cases, the well observations correspond to locations where depth to basalt is well constrained by the refractions. Well 699-55-60A did not reach basalt; thus, this well represents a minimum basalt depth. The depth to basalt in well 699-61-62 and the 3,000 m s⁻¹ are within 5 m, or ∼10 percent, of the total observed depth. The difference may be attributed to uncertainty in the seismic result or to real variation in basalt elevation between the seismic profile and the well. A third nearby well cannot be compared to the Vp result, as the refraction data do not image basalt at that location.

Two check shot velocity profiles are also available along Line A North (Figure 5a). Well 699-55-60A is ∼475 m to the east of the profile and did not reach basalt. It shows two thin, high-velocity (Vp = 2,000 m s⁻¹) layers that are not evident in the refraction result. Well 699-61-62, which did reach basalt at a depth of 54 m (Figure 5a, inset), indicates two Vp layers within the sediment section (Vp < 1,200 and Vp = 2,000 m s⁻¹). The refraction result is intermediate between these two velocities. The basalt depths and velocity estimates in well 699-61-62 closely match the refraction result.

Vs along this profile ranges between 150 and 650 m s⁻¹ (Figure 5b). The 1,000-m s⁻¹ Vp contour closely resembles the transition between Vs less than 500 m s⁻¹ and Vs greater than 550 m s⁻¹, suggesting that this transition differentiates the H2 from the H3 formation.

Line B is an approximately 2.7-km-long northwest-to-southeast trending profile in the northern part of the study area (Figure 1). Here, a similar distribution of velocities is observed in the suprabasalt sediment section as Line A North. The top of the basalt shows a prominent dip at profile distance ∼1,350 m, mirroring the structure depicted in the PSDM image (Figure 6a). Recovered velocities in that area are lower compared to adjacent sections. Combined with the discontinuous structure indicated by the PSDM image, it suggests that this is an area where the basalt surface is irregular and possibly fractured and permeable. On the eastern edge of the profile, either the basalt is too deep for imaging or the velocity of the upper basalt is lower compared to elsewhere along the line.

Four well observations of basalt elevation are within 500 m of Line B (Figure 6a). The well projected onto profile distance 2,400 m is 446 m away and suggests that the refraction result underpredicts the basalt depth and agrees more closely with the PSDM image. The well is far away and may not represent the basalt depth at the seismic profile. Elsewhere along the line, the nearby wells are in better agreement with the refraction result.

Well 699-61-62 is 149 m to the south of Line B. It is evident that the refraction result represents a smoothed version of the velocity profile obtained from the check shot. The refraction result here shows a smooth transition into basalt velocities, which may explain why it overpredicts basalt depth and demonstrates the uncertainty inherent in the refraction method. It is also possible that basalt elevation is different at the well site than it is along Line B.

Vs along this profile shows a similar correlation to Vp structure, indicating that the H1 and H2 formations are not uniformly thick. Near the center of the profile, a high-velocity anomaly in both Vp and Vs suggests a thinning of the H2 layer.

Williams et al. (2011) compiled all the available seismic reflection data in the Gable Gap area and compared them to well observations and found that the reflection data consistently overpredicted the estimated depth to basalt. Much of the data that Williams et al. (2011) reviewed were time-stacked images converted to depth using a smooth-velocity model derived from check shots. Williams et al. (2011) suggested that the depth prediction errors were likely due to a poor time-to-depth-conversion velocity model. That argument does not apply to the PSDM images, as the PSDM approach simultaneously and iteratively produces a high-resolution velocity model directly from the data (e.g., Sheriff and Geldart, 1995). The final output is a reflection image in depth that can be directly compared to the refraction-derived Vp images.

Compared to the interpreted refraction Vp images, basalt depths predicted using the PSDM images do tend to be greater. However, this is not universally true. Along the northern-most end of Line A North (Figure 5a) and along Line C East (Figure S3), the refraction Vp results are remarkably similar to the PSDM images. Ignoring the depth discrepancies, a close similarity between topography along the basalt surface interpreted from the PSDM images and from the refraction Vp images was observed. This suggests that the structural complexity indicated by both approaches can be reliably interpreted even if the absolute elevations are uncertain.

Coincident measurements of basalt elevation interpreted from the Vp results along Lines A North, B, and D suggest that the precision of the refraction method is on the order of ±2–4 m. This likely applies to regions where the basalt surface is laterally uniform. In areas where the basalt has rough topography, the refraction data will tend to smooth out the small-scale variations associated with sediment velocity transitions. A good example of this is the comparison of the Vp image and the check shot from well 699-61-62 (Figure 6a). The Vp tomogram shows a smooth gradient from sediment to basalt velocities, whereas 150 m away, the check shot shows a hard boundary. The reflection image shows rough topography on the basalt at this location. Since the refraction data can resolve only a smoothed version of the structure, the depth to basalt is poorly constrained.

The Vp tomograms show variations in basalt velocity that may indicate areas where the basalt is fractured and likely to be permeable. The highest basalt Vp (∼5,000 m s⁻¹) on Line A North (Figure 5a) occurs where the PSDM image shows a strong coherent reflection, suggesting a smooth basalt surface (between x = 100–300 m and 1,200–2,000 m). On the edges of these sections, Vp is not as well constrained due to a lack of data coverage. However, basalt Vp is lower at transitions in reflection coherence. These areas may represent a rubbly basalt surface or a region of increased fracture density where communication between the unconfined and confined aquifers could potentially occur. However, well observations are needed to confirm this.

At profile distances greater than 3,000 m for Line A North, Vp in the basalt is ∼3,500–4,500 m s⁻¹, and the reflection image indicates several discontinuities. Here, the low Vp along the top of the basalt is well constrained and may indicate fractured and permeable basalt. The combination of low basalt Vp and discontinuous structure in the PSDM image supports an interpretation of folded and/or fractured basalt.

The 3,000-m s⁻¹ Vp contour marks the transition from suprabasalt sediments to basalt and shows good agreement with well observations (Figures 5 and 6). There are a few sections of Line A North where the top of the basalt is too deep to be imaged with the 196-m streamer offset, and these are areas where paleo-channels carved through the basalt have been previously interpreted (Sunwall et al., 2011). Here, the PSDM images were used to interpret basalt elevation with the understanding that they may be biased toward overly deep estimates.

Most of the variation in basalt elevation is supported by the seismic results; however, well observations also show large variations over short distances. For example, two wells north of Line C East show elevation differences of 24 m over a distance of 190 m.

Both the Vp- and the Rayleigh wave–derived Vs images indicate lateral variations in shallow (<20 m deep) sediment properties. These variations likely represent changes in grain size distribution or sandy versus gravelly deposits of the Hanford formation. For example, the dashed lines in Figures 5b and 6b indicate an interpreted boundary between the sandy H2 and the gravel-dominated H3 unit. Given the different grain size distribution in each unit, a difference in porosity and hydraulic conductivity is also expected, and variations in the distribution and thickness of these two units may impact vertical fluid infiltration and unsaturated flow within the vadose zone. Future studies with lower-frequency geophones could potentially image deeper Vs structure and provide the possibility of using Vp/Vs ratios to discriminate between dry and saturated conditions.

The seismic landstreamer is an effective tool for characterizing near-surface hydrostratigraphy at the Hanford Site. While previous work (Hyde et al., 2011; Sunwall et al., 2011) focused on the reflected wave field, this analysis demonstrates that the refraction and Rayleigh wave data can provide additional information on basalt elevations and properties.

The integration of reflection, refraction, Rayleigh wave data, and borehole information has revealed previously undetected features within the study area. The refraction data provided an updated estimate of basalt elevation and allowed the identification of low-velocity zones that coincide with discontinuities in the reflection images. These low-velocity zones may be indicative of fractures supporting vertical communication between the confined and unconfined aquifers. Integrating well observations with the seismic interpretations suggests that the refraction approach provided basalt elevation estimates that agree more closely with well observations compared to PSDM. Both methods support the interpretation of a highly variable basalt surface. Vs structure derived from Rayleigh wave analysis was broadly consistent with Vp structure in the upper 10–20 m, indicating variations in the distribution of Hanford formation units. This information, along with continuing geophysical investigations, will be used to refine the site geologic model.

This document was prepared by the Deep Vadose Zone–Applied Field Research Initiative at Pacific Northwest National Laboratory. Funding for this work was provided by the U.S. Department of Energy (DOE) Richland Operations Office. The Pacific Northwest National Laboratory is operated by the Battelle Memorial Institute for the DOE under contract DE-AC05-76RL01830.