A View Toward the Future of Subsurface Characterization: CAT Scanning Groundwater Basins

In this opinion paper we contend that high-resolution characterization, monitoring, and prediction are the key elements to advancing and reducing uncertainty in our understanding and prediction of subsurface processes at basin scales. First, we advocate that recently developed tomographic surveying is an effective and high-resolution approach for characterizing the field-scale subsurface. Fusion of different types of tomographic surveys further enhances the characterization. A basin is an appropriate scale for many water resources management purposes. We thereby propose the expansion of the tomographic surveying and data fusion concept to basin-scale characterization. In order to facilitate basin-scale tomographic surveys, different types of passive, basin-scale, CAT scan technologies are suggested that exploit recurrent natural stimuli (e.g., lightning, earthquakes, storm events, barometric variations, river-stage variations, etc.) as sources of excitations, along with implementation of sensor networks that provide long-term and spatially distributed monitoring of excitation as well as response signals on the land surface and in the subsurface. This vision for basin-scale subsurface characterization faces many significant technological challenges and requires interdisciplinary collaborations (e.g., surface and subsurface hydrology, geophysics, geology, geochemistry, information and sensor technology, applied mathematics, atmospheric science, etc.). We nevertheless contend that this should be a future direction for subsurface science research

3-D Transient Hydraulic Tomography in Unconfined Aquifers with Fast Drainage Response

We investigate, through numerical experiments, the viability of three-dimensional transient hydraulic tomography (3DTHT) for identifying the spatial distribution of groundwater flow parameters (primarily, hydraulic conductivity K) in permeable, unconfined aquifers. To invert the large amount of transient data collected from 3DTHT surveys, we utilize an iterative geostatistical inversion strategy in which outer iterations progressively increase the number of data points fitted and inner iterations solve the quasi-linear geostatistical formulas of Kitanidis. In order to base our numerical experiments around realistic scenarios, we utilize pumping rates, geometries, and test lengths similar to those attainable during 3DTHT field campaigns performed at the Boise Hydrogeophysical Research Site (BHRS). We also utilize hydrologic parameters that are similar to those observed at the BHRS and in other unconsolidated, unconfined fluvial aquifers. In addition to estimating K, we test the ability of 3DTHT to estimate both average storage values (specific storage Ss and specific yield Sy) as well as spatial variability in storage coefficients. The effects of model conceptualization errors during unconfined 3DTHT are investigated including: (1) assuming constant storage coefficients during inversion and (2) assuming stationary geostatistical parameter variability. Overall, our findings indicate that estimation of K is slightly degraded if storage parameters must be jointly estimated, but that this effect is quite small compared with the degradation of estimates due to violation of “structural” geostatistical assumptions. Practically, we find for our scenarios that assuming constant storage values during inversion does not appear to have a significant effect on K estimates or uncertainty bounds

Hierarchical Geostatistics and Multifacies Systems: Boise Hydrogeophysical Research Site, Boise, Idaho

The geostatistical structure of a heterogeneous coarse fluvial aquifer is investigated with porosity data derived from neutron logs at a research well field (Boise Hydrogeophysical Research Site, or BHRS) that was designed, in part, to support three-dimensional geostatistical analysis of hydrologic and geophysical parameters. Recognizing that the coarse fluvial deposits include subdivisions (units between bounding surfaces), we adopt a hierarchical approach and examine the porosity geostatistics of the aquifer at three scales. At the BHRS, the saturated fluvial deposits as a whole (maximum interwell spacing ~80 m, thickness ~16–18 m) are at hierarchical level 1; five subhorizontal units within these deposits (four cobble-dominated units and a channel sand) can be traced across the central area of the BHRS and are at hierarchical level 2; and subunits (patches or lenses) in one of the level 2 units (Unit 4), are at hierarchical level 3. We use variography and porosity statistics to recognize nonstationarity at hierarchical level 1 and in one of the level 2 units (Unit 4) where the means and variances of porosity differences as a function of lag are not constant between distinct units and subunits, respectively. The geostatistical structure at level 1 is modeled with different horizontal and vertical structures that have different sills (vertical sill greater than horizontal sill). The difference in sills can be explained quantitatively by the summing of weighted sills from all individual units and combined units (i.e., a given pair of different units), where the weights are the proportions of data pairs contributing to the sills at each lag from the individual and combined units. Extension of this analysis leads to a weighted, multistructure form of the variogram function whereby a global experimental variogram in a hierarchical system can be decomposed quantitatively into weighted component individual- and combined-unit (or facies) structures for any number of units or hierarchical levels. Such decomposition of the global horizontal variogram from the BHRS indicates that short-range periodicity in that structure is due to both (1) combined-unit structures associated with patches or lenses at hierarchical level 3 in Unit 4 and (2) variations in thickness of Unit 2. For hierarchical multifacies systems, structure models fit to global horizontal and vertical experimental variograms may not be useful for subsequent stochastic modeling if the system on which the structure models are based is nonstationary

The Anomalous Behavior of SH‐Waves Across the Water Table

Most theoretical studies of seismic wave propagation in a porous medium do not predict a significant increase in SH-wave velocity with increasing water saturation. Although that type of behavior is commonly predicted for P-waves (and confirmed by countless observations), the expectation for SH-waves is a slight decrease in propagation velocity with increasing water saturations. While published measurements of SH-wave velocity in laboratory studies have been supportive of such a slight decrease in velocity, the data have been biased towards high pressures (typical of oil reservoirs at large depths of burial). On the other hand, the few published low pressure laboratory measurements have revealed significantly different results. The authors’ in-situ measurements of seismic wave velocities in a shallow, coarse grained, unconfined alluvial aquifer document a significant SH-wave velocity increase in the transition from the vadose zone to the water table. In one vertical seismic profile (VSP), the P-wave velocity increases by a factor of 4.2 and the SHwave velocity increases by a factor of 2.6. What is not clear at this point is the true nature of the increase. Is the velocity increase an expression of the presence of water in the pores, or does water alter the rigidity of the soil matrix? In addition to the broad-band velocity increase, we have also observed changes in the attenuation of SHwaves across the water table. After correcting for geometric spreading, the amplitude decay observed in the vadose zone has been found to be larger than that observed below the water table. However, the variation in amplitude decay as a function of frequency and the measurements of body wave dispersion were found to be larger below the water table than above. That is, the water saturated soil behavior is consistent with a Voigt solid, but the dry material appears to follow a different model. The authors will discuss these observations in the context of the current debate on Poisson’s ratio and the Vp/Vs ratio

Crosshole Radar Tomography in a Fluvial Aquifer Near Boise, Idaho

To determine the distribution of heterogeneities in the saturated zone of an unconfined aquifer in Boise, ID, we compute tomograms for three adjacent well pairs. The fluvial deposits consist of unconsolidated cobbles and sands. We used a curved-ray, finite-difference approximation to the eikonal equation to generate the forward model. The inversion uses a linearized, iterative scheme to determine the slowness distribution from the first arrival traveltimes. The tomograms consist of a layered zone representing the saturated aquifer. The velocities in this saturated zone range between 0.06 to 0.10 m/ns. We use a variety of methods to assess the reliability of our velocity models. Finally, we compare our results to neutron-derived porosity logs in the wells used for the tomograms. The comparison shows that the trends in porosity derived from the tomograms match the trends in porosity measured with the neutron probe

The Effects of Sparging on P- and SH- Vertical Seismic Profiles

While the introduction of pressurized air into an unconsolidated, coarse-grained fluvial aquifer might well be expected to affect the P-wave velocity profile below the water table, we have found that S-waves are also sensitive to changes induced by air sparging. In a study spanning over a year of sparging, observations of both P- and S-waves were conducted by Vertical Seismic Profiling (VSP). While the primary objective was to characterize the aquifer, we have found that air sparging has significantly affected both P- and S-wave propagation. Below the water table we have observed as much as a 54% decrease in P-wave velocity, and as much as a 31% increase in S-wave velocity after continued sparging. Above the water table, we observe only small changes in both P- and S-wave velocities. This pattern of velocity change (decreasing P, increasing S) may be due to an increase in the amount of trapped air below the water table. Published laboratory studies in the small strain regime have shown P-wave velocities to be sensitive to void ratio, fluid content, and confining stress. On the other hand, most similar studies of S-waves have only been conducted on either dry or saturated samples. However, one recent laboratory study suggests that shear modulus and shear velocity may increase significantly at partial water saturations (due to capillary forces). Data from our in-situ survey supports this more recent lab work. We have observed that S-wave propagation may be significantly altered by fluid content when soils are partially saturated with water (where trapped air may exist, producing a 3-phase fluid-frame system). In addition, we have observed changes in the propagating wavelet. This may be an indication that viscous damping is also affected by partial water saturation. We conclude by observing that S-waves may prove to be an attractive alternative for mapping the effects of air sparging

Hydraulic Conductivity Imaging from 3-D Transient Hydraulic Tomography at Several Pumping/Observation Densities

[1] 3-D Hydraulic tomography (3-D HT) is a method for aquifer characterization whereby the 3-D spatial distribution of aquifer flow parameters (primarily hydraulic conductivity, K) is estimated by joint inversion of head change data from multiple partially penetrating pumping tests. While performance of 3-D HT has been studied extensively in numerical experiments, few field studies have demonstrated the real-world performance of 3-D HT. Here we report on a 3-D transient hydraulic tomography (3-D THT) field experiment at the Boise Hydrogeophysical Research Site which is different from prior approaches in that it represents a “baseline” analysis of 3-D THT performance using only a single arrangement of a central pumping well and five observation wells with nearly complete pumping and observation coverage at 1 m intervals. We jointly analyze all pumping tests using a geostatistical approach based on the quasi-linear estimator of Kitanidis (1995). We reanalyze the system after progressively removing pumping and/or observation intervals; significant progressive loss of information about heterogeneity is quantified as reduced variance of the K field overall, reduced correlation with slug test K estimates at wells, and reduced ability to accurately predict independent pumping tests. We verify that imaging accuracy is strongly improved by pumping and observational densities comparable to the aquifer heterogeneity geostatistical correlation lengths. Discrepancies between K profiles at wells, as obtained from HT and slug tests, are greatest at the tops and bottoms of wells where HT observation coverage was lacking

Recognizing and modeling variable drawdown due to evapotranspiration in a semiarid riparian zone considering local differences in vegetation and distance from a river source

Riparian zones in semiarid regions often exhibit high rates of evapotranspiration (ET) in spite of low-soil moisture content due to the presence of phreatophytic vegetation that is able to withdraw water from shallow aquifers. This work seeks to better define the relationship between ET, the saturated zone and the river boundary by comparing observed water table drawdown records to analytically modeled drawdown in fully penetrating wells of an unconfined aquifer in response to daily ET flux. ET at the Boise Hydrogeophysical Research Site (BHRS), a riparian zone in a temperate, semiarid environment, is calculated using a radiation-based method to provide ET values at four different wells with different vegetation densities. Analytically modeled drawdown response to ET forcing shows that drawdown magnitude increases with increasing distance from the river edge even as the surficial ET forcing remains constant. This behavior is also observed in well hydrographs and shows the buffering effect that flow from the river has on drawdown in fully penetrating riparian wells in highly permeable, unconfined aquifers. Relative contributions of river water to aquifer storage are calculated for ET-induced diurnal fluctuations of the water table at increasing distances from the river boundary. Failure to account for these spatial differences in drawdown related to the river source may explain some errors associated with estimating ET from well hydrographs alone

Hydraulic Tomography: 3D Hydraulic Conductivity, Fracture Network, and Connectivity in Mudstone

We present the first demonstration of hydraulic tomography (HT) to estimate the three‐dimensional (3D) hydraulic conductivity (K) distribution of a fractured aquifer at high‐resolution field scale (HRFS), including the fracture network and connectivity through it. We invert drawdown data collected from packer‐isolated borehole intervals during 42 pumping tests in a wellfield at the former Naval Air Warfare Center, West Trenton, New Jersey, in the Newark Basin. Five additional tests were reserved for a quality check of HT results. We used an equivalent porous medium forward model and geostatistical inversion to estimate 3D K at high resolution (K blocks m3), using no strict assumptions about K variability or fracture statistics. The resulting 3D K estimate ranges from approximately 0.1 (highest‐K fractures) to approximately 10−13 m/s (unfractured mudstone). Important estimated features include: (1) a highly fractured zone (HFZ) consisting of a sequence of high‐K bedding‐plane fractures; (2) a low‐K zone that disrupts the HFZ; (3) several secondary fractures of limited extent; and (4) regions of very low‐K rock matrix. The 3D K estimate explains complex drawdown behavior observed in the field. Drawdown tracing and particle tracking simulations reveal a 3D fracture network within the estimated K distribution, and connectivity routes through the network. Model fit is best in the shallower part of the wellfield, with high density of observations and tests. The capabilities of HT demonstrated for 3D fractured aquifer characterization at HRFS may support improved in situ remediation for contaminant source zones, and applications in mining, repository assessment, or geotechnical engineering

Crosshole Radar Tomography in an Alluvial Aquifer Nearboise, Idaho

To determine the distribution of heterogeneities in an unconfined aquifer in Boise, ID, we compute radar tomograms for three adjacent well pairs. The fluvial deposits consist of unconsolidated cobbles and sands. We used a curved‐ray, finite‐difference approximation to the eikonal equation to generate the forward model. The inversion uses a linearized, iterative scheme to determine the slowness distribution from the first arrival traveltimes. The tomograms consist of a sequence of layers representing the saturated aquifer. The velocities in this saturated zone range between 0.06 to 0.10 m/ns. We use a variety of methods to assess the reliability of our velocity models. Finally, we compare our results to neutron‐derived porosity logs in the wells used for the tomograms. The comparison shows that the trends in porosity derived from the tomograms match well with the trends in porosity measured with the neutron probe