Keywords

7.1 Introduction

Groundwater protection and contaminated site remediation remain a difficult challenge for society due to our inability to characterize subsurface conditions at a sufficiently high level of resolution and accuracy. The success of efforts to protect or remediate a site requires a good understanding of how contaminants move in the saturated subsurface. Contaminant transport is primarily controlled by two physical mechanisms, advection, the movement of dissolved solutes with flowing groundwater, and diffusion, movement of contaminants from areas of high to low concentrations as a result of solute molecular thermal motion. Except in low-permeability materials, such as clays, where flow velocity is very small, advection typically plays a much more significant role than diffusion on the physical transport of contaminants in groundwater.

Groundwater flow can be mathematically related to hydraulic conductivity K [L/T] and gradient i [dimensionless], through Darcy’s Law (Eq. 7.1), whose physical interpretation is illustrated in Fig. 7.1:

$$Q = -K \times i \times A,$$
(7.1)
Fig. 7.1
A schematic of a circular wheel labeled K has the length L. There are two tubes at the top labeled h 2 and h 1. A right arrow points from A to Q beneath the tube h 1.

Schematic of Darcy’s Law flow calculation (Liu and Butler 2019). Groundwater flow rate is proportional to K, which can vary over several orders of magnitudes in a heterogeneous system. Definition of other notations is given in text

where Q [L3/T] is the flow rate across area A [L2] of a porous medium. In groundwater hydrology, it is more common to use the Darcy flux, q [L/T], which is defined as in Eq. 7.2:

$$q = Q/A = -K \times i.$$
(7.2)

The hydraulic gradient is computed following Eq. 7.3:

$$i = \left( {h_{1} -h_{2} } \right)/L,$$
(7.3)

where h1 and h2 are the hydraulic heads at the ends of the computational domain and L is the distance between the points where h1 and h2 are measured. Note that the Darcy flux describes the average flow rate across the total cross-sectional area, which is not the same as the seepage velocity in the pore space. The seepage velocity can be related to the Darcy flux through Eq. 7.4:

$$v = q/\phi_{{\text{e}}} = -K \times i/\phi_{{\text{e}}} ,$$
(7.4)

where ϕe is the effective porosity (dimensionless), i.e., the proportion of pore space (relative to total aquifer volume) through which water can move under the hydraulic gradient. Because groundwater cannot move readily through extremely small pores (e.g., those in clays), dead-end pores, or pores that do not connect to the actual flow paths, the effective porosity is less than the total porosity (Zheng and Bennett 2002).

Groundwater flux (e.g., seepage velocity defined in Eq. 7.4) is a dynamic quantity whose value and direction will change whenever there is a change in K (e.g., borehole construction, dissolution or precipitation of solids) or i (e.g., pumping, episodic recharge events, increase or decrease of boundary heads). Compared to K, direct measurement of groundwater flux is more difficult in the field. As a result, the majority of current contaminant site characterization efforts focus on K (Boggs et al. 1992; Dagan and Neuman 1997; Fogg et al. 2000), although direct flux characterization has received increasing attention in recent years (Devlin 2020). In this chapter, we will discuss the development of high-resolution approaches for K, porosity, and flux characterization.

Conventional K characterization approaches, such as pumping tests, slug tests, or flowmeter profiling, have proven to be of limited effectiveness for obtaining information at the level of detail that is needed to accurately quantify groundwater and contaminant movement in heterogeneous media (Butler 2005; Liu and Butler 2019). There are two major drawbacks associated with the conventional approaches. The first is their dependence on existing wells, which are often sparsely distributed and not necessarily located at places of greatest interest. The second is the many theoretical and procedural limitations of well-based approaches (e.g., in-well hydraulics and short-circuiting flow through filter packs or leaky in-well packers). To address these concerns, a series of high-resolution characterization methods for unconsolidated settings based on direct push (DP) technology have been developed (Butler et al. 2002; Schulmeister et al. 2003; Dietrich and Leven 2005; McCall et al. 2005; Butler et al. 2007; Dietrich et al. 2008; Liu et al. 2009, 2012; Maliva 2016; Liu et al. 2019; McCall and Christy 2020). Compared to the conventional approaches, these DP methods can provide K measurements at a much higher level of detail, and allow the measurements to be acquired in virtually any location.

Nuclear magnetic resonance (NMR), a widely used borehole logging technique in the petroleum industry for characterizing hydrocarbon reservoirs, has been adapted for hydrological applications in recent years (Walsh et al. 2011, 2013). By measuring the responses of water molecules to perturbed magnetic fields, NMR provides information about the total amount of water (i.e., porosity in saturated zones), as well as the distribution of pore sizes in the formation. That information can then be used to estimate K (Dlubac et al. 2013; Walsh et al. 2013; Knight et al. 2016; Kendrick et al. 2021). Furthermore, based on the NMR responses of water molecules in different sizes of pores, porosity can be divided into relatively mobile and immobile portions. As discussed earlier, the effective porosity, which is similar to the mobile porosity determined by NMR, is the key parameter for calculating pore water velocity from Darcy’s flux. The immobile porosity, on the other hand, provides important information that is needed to characterize the diffusion of contaminants from low-K sediments back into high-K zones at sites where contaminant concentrations in those zones have been largely reduced by remediation efforts. Currently, the most common approach for determining porosity is to collect core samples in the field, which are then sent to a laboratory for measurements. Due to the high costs associated with this method, porosity measurements are typically kept at a minimum in most site characterization investigations. Thus, the NMR approach has great potential as a site characterization tool to obtain information on porosity as well as K estimates.

There has been a growing interest in measuring groundwater flux directly in the field as the remediation community switches from contaminant concentration-based decision-making to one that is based on contaminant mass discharge (Suthersan et al. 2010; Devlin 2020). The contaminant mass discharge, a product of groundwater flux and contaminant concentration, requires the groundwater flux to be reliably measured or estimated. Different approaches are available for obtaining groundwater flux estimates (Labaky et al. 2007; Bayless et al. 2011; Devlin 2020). These approaches can be divided into two general categories, (1) Darcy’s law calculation where both K and i need to be estimated independently, and (2) tracer-based approaches that provide a direct estimate of groundwater flux assuming a quantifiable relationship between the flux and tracer responses. For the Darcy’s law approach, the calculated flux is typically an average value over a relatively large area. As a result, the resolution is often not suitable for contaminant remediation activities. For tracer-based approaches, heat and solutes are two common tracers and responses are monitored by measuring subsurface temperature and solute concentrations (or concentration surrogates). By injecting and monitoring a small amount of tracer with an electrical conductivity different from the ambient formation fluid, centimeter scale groundwater velocity measurements can be made by a point velocity probe (PVP; Labaky et al. 2007). Initial applications of the PVP focused on deployment of the probe in unconsolidated formations, so that the probe was in direct contact with porous media and the impacts of borehole or well construction were minimized. Subsequent developments extended the PVP approach to wells and open boreholes to allow for measurements at as many depths as needed (Osorno et al. 2018), and in fractured rocks (Heyer et al. 2021).

A major limitation of tracer-based approaches is that only a few measurements can be acquired at pre-selected locations at a single time. Over the last decade, technological developments in fiber-optic distributed temperature sensing (DTS) have led to a significant improvement in using heat to measure and monitor various hydrological processes (Selker et al. 2006; Lowry et al. 2007; Moffett et al. 2008; Henderson et al. 2009; Leaf et al. 2012; Striegl and Loheide 2012; Liu et al. 2013; Read et al. 2013; Bakker et al. 2015; Maldaner et al. 2019; Munn et al. 2020; Simon et al. 2021). Due to the high-resolution temporal (sub-minute) and spatial (cm scale with wrapped DTS cable) temperature measurements by DTS, continuous information on groundwater flux distribution can potentially be obtained by tracking heat movement along the entire vertical interval of the borehole. However, a significant challenge for heat-based approaches is that density-driven buoyancy flow can affect the measurements, especially when the borehole is open and/or there is a significant amount of annular space between the measuring device and well casing. In some cases, separating the impact of thermal conduction from heat advection can also be a challenge if the rate of groundwater flux is small.

Butler (2005) provided an overview of the major methods for shallow subsurface characterization of K at that time. Many developments have occurred since then, particularly with the four approaches described in the previous paragraphs (DP, NMR, PVP, and DTS). The purpose of this chapter is to provide an overview of those developments. Specifically, we will discuss the characterization of K by DP approaches, the characterization of K and porosity by NMR, and the characterization of groundwater flux by PVP and DTS. For DP and NMR approaches, we will demonstrate their performance with data collected at a long-term research site of the Kansas Geological Survey (KGS). The DTS flux approach will be demonstrated in a laboratory sand tank while the PVP approach will be demonstrated with data from a sand aquifer at Borden, Canada. We conclude the chapter by summarizing the advantages and disadvantages of these approaches and offering some suggestions on future work to address key limitations.

7.2 Characterization of Hydraulic Conductivity by Direct Push Approaches

7.2.1 Direct Push Technology

Direct push is a drilling method that advances a rod string into unconsolidated materials using hydraulic rams supplemented with vehicle weight and, depending on the application, high-frequency percussion hammers (Fig. 7.2). In contrast to conventional rotary drilling, no material is removed from the subsurface. Instead, the sediments in the immediate vicinity of the rod string are pushed aside during advancement. The major advantages of DP technology include minimal site disturbance, practical elimination of drill cuttings, reduced exposure of field personnel to hazardous materials at contaminated sites, lower time and costs (costs reduced by 40–60% compared to rotary drilling; personal communication, Wes McCall, Geoprobe Systems, Nov. 11 2021), and high mobility of the equipment (DP rigs can work under many field conditions). The main disadvantage is that DP is limited to unconsolidated sediments with penetration depths typically less than 30 m. The shoving aside and thus repacking of near-rod materials, which may lead to a change in their hydraulic properties, can sometimes have adverse impacts on field results, especially in finer-grained sediments.

Fig. 7.2
A schematic of a direct push rig has a water reservoir, data acquisition, controller, H P T pump, trunkline, and track-mounted direct push rig. Below is a line graph that plots the depth below the land surface versus the flow rate, pressure, and E C. It has a fluctuating trend.

A direct push rig and an example HPT profile (image of DP rig courtesy of Geoprobe Systems)

DP technology has been used for a wide array of subsurface investigations ranging from characterization of aquifer properties (Butler et al. 2002; Lunne et al. 2002; Butler 2005; Lessoff et al. 2010; Liu et al. 2012; Liu and Butler 2019; McCall and Christy 2020), groundwater and soil sampling (Artiola et al. 2004; Schulmeister et al. 2004; Nielsen and Nielsen 2006; U.S. EPA 2016), subsurface injection of remedial reagents (Stroo and Ward 2010), and installing wells and monitoring devices for long-term water resources management (ITRC 2006; Liu et al. 2016). In the area of hydrogeologic characterization, a series of innovative approaches have been developed for DP applications to allow information about K to be obtained at a resolution, accuracy and speed that was previously not possible (McCall et al. 2002; Butler et al. 2007; Dietrich et al. 2008; Liu et al. 2009, 2019; McCall and Christy 2020). In this section, we focus on using the DP electrical conductivity (EC) probe, DP permeameter (DPP), DP injection logger (DPIL), Hydraulic Profiling Tool (HPT), and High-Resolution K (HRK) tool to characterize vertical K profiles in unconsolidated settings. Results from the Larned field site of the KGS are used to illustrate and compare the performance of these four approaches.

7.2.2 Larned Research Site

The Larned Research Site (LRS), which was established by the KGS for studying stream-aquifer interactions and riparian zone processes in the early 2000s, is located adjacent to a US Geological Survey (USGS) stream-gaging station on the Arkansas River near the city of Larned in west-central Kansas (Butler et al. 2007; Fig. 7.3). The shallow unconsolidated subsurface at the site consists of three major units: the Arkansas River alluvial aquifer, a clay-dominated confining layer in the middle, and the High Plains Aquifer (HPA) resting on Pennsylvanian bedrock (shale and limestone) at the bottom (Butler et al. 2004). The alluvial aquifer contains many intermittent clay lenses, as indicated by the abrupt increase of EC values particularly in the lower two-thirds of the interval. The upper portions of the HPA also contain intermittent clay lenses where the EC values are elevated. The thickness of each hydrogeological unit varies significantly across the site, with the depth to water around 3–5 m. Due to extensive pumping from the HPA, groundwater levels have been declining since the late 1950s in western Kansas where the upstream reaches the Arkansas River are located. As a result, the Arkansas River near the LRS has been dry for much of the last few decades except during wet years or after intensive storms.

Fig. 7.3
A map plots the stream gaging stations, larned research site, and country line. On the top are an enlarged location of Kansas and an aerial view of the Larned research site. Below the view is a line graph of depth below the land surface versus E C.

Overview of the Larned Research Site. The image on the upper right is from Google Earth (accessed on 10/7/2013). The generalized stratigraphy on the lower right was determined from DP EC (Wenner array) profiling (Butler et al. 2004). Low EC values indicate sand and gravel, while high values indicate clay at this site. Notice the jump in electrical conductivity at the water table where the probe moved from very dry to saturated sands and gravels

7.2.3 Direct Push Electrical Conductivity

The EC of subsurface saturated zones is primarily determined by the pore-fluid chemistry, clay content (electrically conductive; not all clays are electrically conductive), and total porosity. At freshwater sites where clay mineralogy is consistent and variations in pore-fluid chemistry are negligible, such as the LRS, EC is largely a measure of the clay content (porosity variations typically have a much less significant impact than clay content) and can thus be used as a qualitative indicator of K. Because no water injection is needed during probe advancement and measurement, EC profiling can be performed with significantly less effort and more rapidly than other DP approaches.

Figure 7.4 displays three EC logs collected across the site. At LWC2, the EC values are generally the highest in the upper alluvial aquifer, indicating that the K of the alluvial zone is lower at this location than at the other two locations. The EC values are much higher for the middle confining clay layer in all profiles. The thickness of the confining layer is similar between LWC2 and LWPH9 (interval extending from depths of 10–15 m), but larger at LEC2 (interval extending from 10 to 19 m). The EC variations between the profiles indicate that the clay distribution is highly heterogeneous across the site. Below the middle confining layer, the EC values decrease to around 30 mS/m in the HPA. The thickness of the HPA is similar between LWC2 and LWPH9, but smaller at LEC2. The EC values show an abrupt increase once the probe enters the shale below the HPA.

Fig. 7.4
Three line graphs A to C of D P E C logs for L W C 2, L W P H 9, and L E C 2, respectively. The graph plot depth below the surface in meters versus electrical conductivity in milliseconds per meter. They plot a fluctuating line.

DP EC logs at three locations at the LRS: a LWC2, b LWPH9, and c LEC2. See Fig. 7.3 for the locations of LWC2, LWPH9, and LEC2. An example DP EC probe with the Wenner sensor array is shown at the upper right (image of EC probe courtesy of Geoprobe Systems)

Although EC provides a good indicator of relative K when electrically conductive clay is a dominant factor in subsurface lithology, only a small number of studies have explored the use of EC for K estimation in unconsolidated clayey formations (see the review paper by Purvance and Andricevic 2000 and cited literature). This is largely because K is very difficult to measure in clay-rich formations where K is low and hydraulic tests (typically used to obtain independent K estimates) would take a long time to complete in the field. As a result, current K characterization efforts are mainly focused on moderate to high-K zones (Bohling et al. 2012). Nonetheless, future work is needed to exploit the potential of EC as a means to estimate K in unconsolidated materials having clay as a major constituent. On the other hand, in zones that contain little clays such as the HPA (Fig. 7.4), where the pore size distribution likely has the most significant impact on aquifer permeability, the EC does not provide an effective measure of K (Schulmeister et al. 2003).

7.2.4 Direct Push Permeameter

The DP permeameter (DPP) is a tool that can be used to obtain reliable K estimates on a scale of relevance for most contaminant transport investigations (Stienstra and van Deen 1994; Lowry et al. 1999; Butler et al. 2007; Liu et al. 2008). It consists of a short cylindrical screen attached to the lower end of a DP rod and two pressure transducers inset into the rod above the screen (Fig. 7.5). During DPP advancement, water is continuously injected through the injection screen to prevent clogging. After the measurement depth is reached, advancement ceases and water injection is stopped to allow the heads in the aquifer to recover to background conditions. After the aquifer heads recover, a series of short-term injection tests are performed, and K is estimated from the spherical form of Darcy’s Law (Eq. 7.5) using the injection rate and the injection-induced pressure responses at the two transducers (Butler et al. 2007),

$$K_{{{\text{DPP}}}} = \frac{Q}{{4\pi \left( {\Delta h_{1} - \Delta h_{2} } \right)}}\left( {\frac{1}{{r_{1} }} - \frac{1}{{r_{2} }}} \right),$$
(7.5)
Fig. 7.5
A schematic of the D P P and 2 line graphs are labeled a to c. The injection line, pressure transducer, port, and injection screen are labeled in a. Graphs plot the injection rate and pressure head difference for b, and the top and bottom pressure heads for c versus time in a fluctuating trend.

a Schematic of the DPP, b injection rate and injection-induced pressure head difference between the top and bottom transducers for a DPP test at a depth of 6 m below land surface at the LRS (location LWPH9), c measured pressure heads at the top and bottom transducer for this series of tests. In a, the injection lines at the pressure transducer ports were added during the development of the HRK tool

where Q is injection rate, Δh1 and Δh2 are the injection-induced pressure head changes at pressure transducers PT1 and PT2, and r1 and r2 are the distances from PT1 and PT2 to the injection screen, respectively (Fig. 7.5). \(K_{{{\text{DPP}}}}\) is a weighted average over the interval between the screen and the farthest transducer; material outside of that interval has little influence on the estimated K (Liu et al. 2008). The major advantage of the DPP is that it only requires steady-shape conditions (constant hydraulic gradient with time) as opposed to the steady-state condition (constant head with time) required by most other pumping-based approaches (e.g., flowmeter tests). The steady-shape analysis allows for shorter test times in the field. In addition, the resulting K estimates are not sensitive to the low-K skin that can potentially develop due to material compaction during tool advancement. Background water level fluctuations due to regional stresses (e.g., well pumping, stream stage changes) also have a minimal impact on DPP results because the fluctuations are largely cancelled out using the steady-shape analysis.

Figure 7.5 shows an example of DPP test data collected at the LRS (location LWPH9) at a test depth of 6 m. Three injection tests were performed. The first test used an injection rate of ≈800 mL/min, which was started at a time of 120 s. The injection-induced pressure head difference (i.e., the pressure head difference between the two pressure transducers during injection minus the pressure head difference prior to injection) was 0.07 m. The injection-induced pressure head difference in the first test determined how the injection rate was adjusted in the next test. Because 0.07 m is a relatively small value, the injection rate was raised to ≈1600 mL/min in the second test. In the third test, the injection rate was adjusted back to the level of the first test to check if the injection-induced pressure returned to the level observed during the first test. For the three tests shown in Fig. 7.5, the calculated K values were 7.4, 7.7, and 7.4 m/d, respectively. The consistency between the K estimates from different tests is a good indication of the reliability of DPP results; a lack of consistency could be due to sensor instability (e.g., air bubbles in the injection water can cause unstable sensor readings) or formation alteration during DPP tests (e.g., vertical channeling along the probe surface).

Figure 7.6 shows the DPP K results for the shallow alluvial aquifer and the deeper HPA at LWC2, LWPH9, and LEC2. No DPP tests were performed in the middle low-K clay layer due to time constraints (with current flow instrumentation, one DPP test in this layer would take days to weeks to complete). Consistent with the EC profiles, K shows a much greater variability with depth in the alluvial aquifer than in the HPA. Except for a few thin layers, the K in the alluvial aquifer is lower than that in the HPA. The average DPP K value for the alluvial aquifer is 25 m/d, 25 m/d, and 70 m/d at LWC2, LWPH9, and LEC2, respectively, while the average DPP K value for the HPA is above 200 m/d at all three locations.

Fig. 7.6
Three line graphs plot the depth below the surface versus K for L W C 2, L W P H 9, and L E C 2. The graphs plot line and spot for D P I L-K and D P P-K in fluctuating trend.

DPP and DPIL K estimates at three locations at the LRS: LWC2 (left graph), LWPH9 (middle), and LEC2 (right). DPP values are calculated using the spherical form of Darcy’s Law. The calibrated a and b are 5.8 and −12.4 for computing DPIL K (b changes with the unit of DPIL flow/pressure head ratio; see Eq. 7.6). Neither DPP nor DPIL K estimates were obtained in most clay layers. For DPIL, the injection pressure is so large that it exceeds the pressure sensor limit in these layers. For DPP, both the time for recovery of heads to background conditions and for the injection-induced pressure difference to stabilize exceed practical time constraints in the field

The current DPP tool provides a spatial resolution of ≈0.4 m (i.e., the distance between the screen and farthest transducer; Fig. 7.5). Although a higher resolution is possible by advancing the tool a smaller distance between measurements and analyzing the set of tests together with a numerical model (Liu et al. 2008), practical time constraints (a DPP test sequence requires 10–15 min in moderate to high-K intervals) usually limit the measurement spacing to 0.4 m or larger in most cases (Butler et al. 2007; Knight et al. 2016). As such, the DPP does not provide information about conditions between test intervals. In addition, it can produce anomalous results in highly stratified formations if no information about the stratification is available (Liu et al. 2008; Zschornack et al. 2013).

7.2.5 Direct Push Injection Logger

The DP injection logger (DPIL) is a tool that can be used to rapidly obtain high-resolution information about relative variations in K (Dietrich et al. 2008). The original design of the tool consists of a short cylindrical screen attached to the lower end of a DP rod immediately behind the drive point (Fig. 7.7a). After the tool is advanced to a depth at which a K measurement is desired, advancement ceases. Water is then injected using a sequence of different rates while the injection rate and pressure are measured at the surface. Line losses between the surface and downhole injection screen can be estimated analytically (assuming laminar flow conditions) or from a regression analysis between pressures at different injection rates. The line losses are removed from the total injection pressure measured at the surface, and the ratio of injection rate over the corrected injection pressure head can then be used for K estimation. In Dietrich et al. (2008), the DPIL ratio was converted into K estimates by regression with K estimates from nearby slug tests. In the original tool design, the scale of each measurement is the vertical length of the injection screen (0.025 m). Because K measurements were made after halting probe advancement, the original DPIL profiling procedure is referred to as the discontinuous advancement mode. DPIL can also be performed while the probe is advanced into the subsurface; this profiling procedure is referred to as continuous DPIL (Liu et al. 2012).

Fig. 7.7
Two schematics of D P I L and H P T are labeled a and b. D P I L consists of a land surface, injection pressure, and water injection, and H P T consists of a D P I L injection port and an E C array.

a Schematic of the original DPIL with a screen at the lower end of the probe rod (after Dietrich et al. 2008), b artistic rendering of the HPT (continuous DPIL probe combined with an EC Wenner array; image courtesy of Geoprobe Systems)

7.2.6 Hydraulic Profiling Tool

The Hydraulic Profiling Tool (HPT) is a DP tool developed by Geoprobe Systems that combines continuous DPIL with EC profiling (Geoprobe 2007). Profiling K variations with the HPT has now been established as a standard practice by the American Society for Testing and Materials (ASTM 2016). Compared to the original DPIL, HPT consists of a single screened port inset into a DP rod at some distance above the drive point, with an EC sensor array between the drive point and the screened port (Fig. 7.7b). As the tool is advanced, water is injected continuously through the screen. The injection rate is monitored at the surface, while the pressure is measured directly behind the screen, essentially eliminating the need for a line loss correction between the pressure sensor and injection screen. The distance between the injection port and the bottom drive point (0.36 m) largely reduces the impact of the mechanical stress produced by the advancement of the drive point during continuous profiling. Because halting probe advancement is not required for injection pressure measurement (except for rod changes or pressure dissipation tests; during dissipation tests, probe advancement is halted, water injection is stopped, and the recovery of pressure to the background level is monitored; dissipation tests provide an estimate of the hydrostatic pressure at different depths and two such tests are usually sufficient per profile), HPT profiling is very rapid in the field, and up to six profiles of 12 m in length can be obtained in a day under good conditions (Bohling et al. 2012). HPT produces a K estimate every 0.015 m, thus providing unprecedentedly high-resolution information about K in the saturated zone across the entire interval traversed during probe advancement.

Figure 7.8 shows the injection rates and pressures measured by HPT at LWC2, LWPH9, and LEC2. The injection rate was set to ≈250 mL/min in the permeable zones (the injection rate decreased in lower K materials because of the high injection pressures). For the middle clay layer, the injection rate became zero when the injection pressure exceeded the upper limit of the pressure sensor (shown as gaps in the injection pressure curves). For all three profiles, the injection pressure is smallest in the HPA and in the lower portion of the alluvial aquifer, indicating K is highest at those depths. As discussed earlier, the ratio of injection rate over injection pressure head, which is also plotted in Fig. 7.8, can be used for K estimation by relating those ratios to nearby K estimates obtained from other approaches.

Fig. 7.8
Three line graphs of the depth below the surface versus injection flow and injection flow or pressure head and injection pressure head plot Q, P, and Q slash P. The graphs have three lines that are in a fluctuating trend.

Continuous DPIL logs at three locations at the LRS: LWC2 (left graph), LWPH9 (middle), and LEC2 (right). Logs are obtained using the HPT manufactured by Geoprobe Systems. The injection rate (red) and rate/pressure head (blue) reference the top axis, and the injection pressure head (green) references the bottom axis. The injection pressure head was calculated as the total injection pressure head measured by the sensor minus the hydrostatic pressure head (estimated by pressure dissipation tests at different depths). For the middle clay layer, the injection rate became zero when the injection pressure exceeded the upper end of the pressure range of the sensor (shown as gaps in the pressure curve)

The major advantage of the discontinuous DPIL mode is that the quality of K measurement is relatively high as the injection rate can be varied at each depth interval to assess if consistent K estimates are obtained across multiple rates. However, no information is available for the intervals between the discontinuous measurements. The frequent suspension of advancement for measurements also makes the discontinuous DPIL profiling less time-efficient than HPT. As a result, HPT is preferred by most practitioners.

The major advantages of HPT are the speed and resolution at which K information can be obtained. However, the current HPT has two limitations that are important for field applications. The first is that HPT has both a lower and upper limit for estimation of K. In high-K zones, the injection-induced pressures are small and may become indistinguishable from the HPT system noise. In low-K zones, on the other hand, the injection-induced pressures are high and may exceed the upper limit of the pressure sensor. In addition, it is difficult, if not impossible, to assess if vertical channeling along the probe surface or formation alteration is occurring as a result of the elevated injection pressures. McCall and Christy (2010) estimated the K range that can be reasonably assessed by HPT to be 0.03–25 m/d under typical conditions. Note that both the upper and lower K limits can be extended by adjusting the injection flow rates; however, this requires modifications of the current flow control and measurement system (Bohling et al. 2012; Liu et al. 2012, 2018). The second limitation of HPT is that it only provides a relative indicator of K. To obtain improved K estimates, additional calibration data are needed to develop a site-specific relationship between HPT measurements and K. Dietrich et al. (2008) and McCall and Christy (2010) used K values from nearby slug tests and performed regression analyses to develop empirical relations for converting the DPIL ratios into K. Zhao and Illman (2022) combined HPT with inverse modeling and hydraulic tomography to improve the estimation of K in a highly heterogeneous site. Liu et al. (2009) used DPP tests, which were performed at the same location as the DPIL, to transform the DPIL ratios into K estimates; the resulting approach is referred to as the High-Resolution K (HRK) tool and is discussed in the following paragraphs. Borden et al. (2021) developed a physically based equation to derive K values from HPT data and an empirical hydraulic efficiency factor.

7.2.7 High-Resolution K (HRK) Tool

As mentioned above, the DPP provides reliable K estimates at a relatively coarse resolution (0.4 m), while the continuous DPIL can provide a measure of the relative variations of K at a much finer resolution (0.015 m). The HRK tool was developed to better realize the potential of these approaches by coupling the DPP and DPIL into a single probe (Liu et al. 2009). The tool has a similar appearance to the DPP (i.e., two pressure transducer ports above the bottom injection screen; Fig. 7.9). The difference is that water is injected through both pressure transducer ports during tool advancement (water also injected through the bottom screen to prevent clogging); the DPP only injected water through the bottom screen during advancement. By injecting water through the transducer ports and measuring injection responses during tool advancement, the HRK tool essentially functions in the continuous DPIL mode. Next, at selected depth intervals, tool advancement ceases and the DPP injection tests are performed. Because the DPIL and DPP measurements are collocated, the DPP data can be used to directly transform the DPIL ratios into K without the need to compare measurements at different support scales as in the regression-based approaches. In the HRK analysis, the DPIL K is calculated using a power-law relation as shown in Eq. 7.6 (Liu et al. 2009):

$$K = 10^{b} \left( {{\text{DPIL}}} \right)^{a} ,\,{\text{or}}\,\log_{10} \left( K \right) = a\,\log_{10} \left( {{\text{DPIL}}} \right) + b,$$
(7.6)
Fig. 7.9
Two photographs of high-revolution K tools are labeled a and b. A has 2 pressure transducer ports above the bottom injection screen, and b is a close-up of the pressure transducer screen.

Prototype HRK tool: a Picture showing the bottom injection screen with the two pressure transducers inset into the rod above the screen, and b expanded view of the pressure transducer screen. Water is injected through the bottom screen and both transducer ports during probe advancement

where DPIL is the ratio of injection rate (mL/min) over pressure head (m); a and b are empirical coefficients whose values are determined through calibration. Specifically, the DPIL K values calculated based on (7.6) are used as input to a numerical model that simulates the DPP test responses. The values of a and b are adjusted such that the simulated DPP responses match the observed DPP data in the field. Instead of calibrating a and b for each HRK profile individually, a single set of a and b is used to achieve the best overall match for all profiles at a site (Liu et al. 2009).

Figure 7.6 shows the DPIL and DPP K estimates from HRK profiling at LWC2, LWPH9, and LEC2; K estimates were not available for the middle clay layer as discussed earlier. Overall, the high-resolution DPIL K estimates are consistent with the DPP values calculated from the spherical form of Darcy’s Law. As a result of the high resolution, the HRK profiles provide quantitative K estimates for thin layers that would be difficult to identify with other approaches, including the DPP when used alone.

7.2.8 Summary of Direct Push Approaches

As no flow injection is involved, DP EC profiling can be used to obtain a significant amount of information about the shallow unconsolidated subsurface in a short time. However, because of the insensitivity to K variations in sediments with little clay, it is commonly applied as a site screening tool to develop a general understanding of hydrostratigraphy and identify zones for further interrogation by more quantitative approaches such as DPIL, DPP, or HRK. For zones with electrically conductive clays as a major constituent, such as the alluvial sediments at the LRS, EC could potentially be converted into semi-quantitative K estimates; further work is needed to establish the relations between EC and K in these settings.

Among the various hydraulically-based DP approaches, continuous DPIL (i.e., HPT) is one of the most powerful approaches and can be used to obtain relative K variations at many locations across a site in a time-efficient manner. However, the relations for converting DPIL data into quantitative K estimates are typically site-specific and independent data are needed for DPIL calibration. Therefore, for sites with significant heterogeneity in K, we recommend a combination of DPIL and HRK, with the DPIL profiles covering the entire site and the HRK profiles performed at a few strategic locations for detailed DPIL calibration. In addition, due to the practical limits of the instrumentation, the current HPT has a rather restrictive range for K measurement. The DPP (or the DPP component of HRK), on the other hand, can be applied in most unconsolidated settings except those with low-K (e.g., < 0.001 m/d) where the time required to complete each test may be hours to days and is too long for practical investigations.

7.3 Characterization of Hydraulic Conductivity and Porosity by Nuclear Magnetic Resonance Profiling

7.3.1 Nuclear Magnetic Resonance

Over the past decade, significant progress has been made in adapting proton nuclear magnetic resonance (NMR) profiling, a widely used borehole logging technique in the petroleum industry, to various applications in near-surface hydrology (Walsh et al. 2011, 2013; Knight et al. 2016; Krejci et al. 2018; Kendrick et al. 2021). This approach involves measuring the response of hydrogen atoms to a series of magnetic perturbations (Fig. 7.10). That response is a function of, among other things, water-filled porosity and the pore-size distribution of subsurface materials. The current NMR tools can provide K and porosity estimates that are averaged over a 0.10- or 0.25-m vertical interval, although a higher resolution can potentially be achieved by advancing the tools at a finer interval and analyzing the entire profile data with global optimization techniques.

Fig. 7.10
A schematic of N M R labeled a, highlights the sensitive shell and disturbed zone. Below are 2 schematics labeled large pores in coarse sands and small pores in sands of varying grain size. They point at graphs plotting magnetization versus time and plot a fluctuating trend.

Schematic of NMR for groundwater applications (images courtesy of Vista Clara Inc.): a The measurement domain is a thin shell around the tool that is suspended in a borehole (after Walsh et al. 2013), and b the relaxation of hydrogen atoms is affected by the pore size distribution. The yellow color in a indicates a disturbed zone from borehole drilling; the NMR measurement domain is typically outside the disturbed zone. There are two key NMR relaxation characteristics shown in b: first, the initial magnetization (M0) depends on the total amount of hydrogen atoms (total water content), and second, the relaxation rate is fast in small pores and slow in large pores

Three steps are typically involved during borehole NMR measurements (Dunn et al. 2002; Walsh et al. 2013). First, after the tool is moved to a measurement depth, the nuclear spins of hydrogen atoms in the surrounding materials are allowed to align (equilibrate) with a static magnetic field (B0) created by permanent magnets in the tool. Second, an alternating magnetic perturbation (B1) is applied at the resonant frequency, during which the hydrogen atoms absorb energy and precess toward the plane that is perpendicular to the static field. Third, the magnetic perturbation is turned off and the excited hydrogen atoms relax back to equilibrium with the static field. During the relaxation, hydrogen atoms will emit magnetic energy that can be detected by the induction coil in the tool. The amplitude of initial magnetization (M0 in Fig. 7.10) at the time of the magnetic perturbations being turned off is a function of the total water content (i.e., porosity in the saturated zone), while the decay of the NMR signal is dependent on the pore size distribution among other factors (e.g., surface geochemistry of grains). Note that the resonant frequency of the measurement volume is not a single value, but a range of values due to the magnetic gradient caused by the permanent magnets on the tool. As a result, a sequence of B1 perturbations at different frequencies (i.e., repeats of the second and third steps) is typically applied to improve the NMR data signal-to-noise ratio (Walsh et al. 2013).

The relaxation of NMR magnetization is usually described by an exponential function (Fig. 7.10). The relaxation rate is fast (short relaxation time T2) in small pores and slow in large pores (long relaxation time). The actual porous media are composed of a network of pores with different sizes, so a distribution of T2 values (typically uniformly spaced on a log scale) are predefined for the NMR data analysis. The amplitude of initial magnetization for each predefined T2 is estimated by a least squares fit between the NMR relaxation data and the exponential functions for all T2 values.

The initial magnetization provides the measurement of total porosity in saturated zones. K can be estimated from the NMR-determined porosity and relaxation time distribution. Different empirical relations between K and NMR signals have been used in the petroleum industry (Dlubac et al. 2013; Kendrick et al. 2021). One of the most commonly used of those approaches for hydrological applications is the Schlumberger Doll Research (SDR) equation (Eq. 7.7) (Knight et al. 2016),

$$K_{{{\text{NMR}}}} = b\phi^{m} (T_{{2{\text{ML}}}} )^{n} ,$$
(7.7)

where b (m/s3), m, and n are empirically determined constants; ϕ is the porosity determined from the NMR initial magnetization; T2ML is the arithmetic mean of log relaxation times weighted by the amplitudes of initial magnetization. b is often referred to as the lithologic constant and contains information about all the parameters affecting permeability except porosity and relaxation time distribution (e.g., the surface-area-to-volume ratio of the pore space). Based on comparing the NMR and DP K data at three sites, Knight et al. (2016) suggested a universal set of b, m, and n could be used for unconsolidated sand and gravel aquifers (b = 0.05 to 0.12 m/s3, m = 1, n = 2).

7.3.2 Nuclear Magnetic Resonance Application at Larned Research Site

Figure 7.11 shows the NMR relaxation data acquired at LWC2, LWPH9, and LEC2 at the LRS. Data were collected by a DP version of the Javelin tool manufactured by Vista Clara Inc, for which the measurement domain was a 0.5-m vertical shell at approximately 14.5 cm from the center of the probe (radial thickness of the shell is 2 mm). As the DP rods used for NMR tool deployment could not be advanced through the middle tight clay layer (a more powerful DP rig was not available during the field campaign), measurements were made in the upper alluvial aquifer only. At all three profiles, the lower portions of the alluvial aquifer have relatively higher amplitude of the initial magnetization with longer relaxation times (dark red colors), indicating that there are more large pores at those depths.

Fig. 7.11
Three heat maps plot the depth below the surface versus relaxation time for L W C 2, L W P H 9, and L E C 2. The high initial magnetization amplitude is between 10 and 8 meters on the y-axis.

NMR relaxation data at three locations at the LRS: LWC2 (left graph), LWPH9 (middle), and LEC2 (right). Data were acquired using a DP version of the Javelin tool (diameter 6 cm) manufactured by Vista Clara Inc. The color indicates the amplitude of magnetization fitted for each relaxation time T2 using the exponential decay function (Fig. 7.10). The red color means high initial magnetization amplitude (i.e., larger fraction of water for that T2), and the blue color means low initial magnetization amplitude (i.e., smaller fraction of water for that T2). The sum of the initial magnetization across all T2 represents the total water content at that depth; the relationship between the initial magnetization and water content is tool specific and is determined by laboratory calibration

Figure 7.12 shows the NMR porosity and K estimates at LWC2, LWPH9, and LEC2. The sharp increase in the water-filled porosity around a depth of 4 m indicates the water table. The NMR K values are calculated using Eq. 7.7 with b = 0.08 m/s3, m = 1, and n = 2. Overall, the NMR and DPP K estimates are consistent with each other. For a few low-K zones, however, the NMR K values are significantly higher than the DPP values (i.e., depth 7 m at LWC2, depths 5 and 7 m at LWPH9). This is likely because the SDR coefficients in Eq. 7.7 are primarily calibrated for relatively permeable zones. A different set of SDR coefficients may be needed to estimate K from NMR data in low-K zones where the fine-grained materials may play an important role in affecting formation permeability.

Fig. 7.12
Three line graphs of the depth below the surface versus K and N M R porosities plots N M R-K, D P P-K, and N M R porosity. The graphs depict a decreasing trend with fluctuations.

NMR porosity and K estimates at three locations at the LRS: LWC2 (left graph), LWPH9 (middle), and LEC2 (right). The DPP K estimates are also plotted for comparison with NMR results. The DPP and NMR profiles are within 1 m of each other at each location

The NMR tool used at the Larned site produced K and ϕ estimates that were averaged over a 0.5-m vertical interval; more recent NMR tools can provide measurements over a 0.25 or 0.10 m interval. Each measurement required about 5 min per interval. Higher resolution is possible by advancing the tool at a small interval (e.g., decimeter) and analyzing the entire profile data with global optimization techniques. As a subsurface characterization tool, NMR has two advantages over other approaches. First, it provides a direct measure of porosity in the saturated zone (and moisture content in the unsaturated zone), while most of the other approaches do not provide any information about porosity. Second, NMR can potentially be a very powerful tool in low-K zones, which are known to present a significant issue at many contaminated sites through the slow, persistent release of contaminants into more permeable zones. In low-K zones, hydraulic-based approaches, such as pumping and slug tests, are not as effective as in coarser materials because test durations are long (e.g., days to months or longer; Butler 2019). For the DP injection-based approaches (DPIL, DPP, HRK), a particular challenge is the significant amount of pressure generated by tool advancement when K is low (Liu et al. 2019). For example, the pressure generated by tool advancement may overwhelm DPIL injection pressure when HPT profiles are performed in silt and clay layers. The time to complete a DPP test in low-K settings can last from hours to days (typically too long for practical applications). NMR, on the other hand, can provide rapid measurements of formation properties in low-K zones because no flow injection is needed and the tool advancement pressure, if NMR is deployed by DP, has little impact on NMR responses. Future work is needed to further investigate the use of NMR for characterizing K and ϕ in silts and clays.

NMR has been increasingly used for environmental investigations in semi- and fully consolidated rock, including shale, limestone, sandstone, and fractured granite (personal communication, David Walsh, Vista Clara, Nov. 22 2021). Recent developments on combining NMR measurement with DP in unconsolidated formations have allowed NMR K and ϕ estimates to be made at a vertical resolution of under 10 cm (Vista Clara 2021).

7.4 Groundwater Velocity Characterization

7.4.1 Characterization of Velocity by Distributed Temperature Sensing

7.4.1.1 Distributed Temperature Sensing

Heat has been used extensively as a tracer to study groundwater and its interactions with other systems (Anderson 2005; Constantz 2008; Rau et al. 2014). As a temperature measurement technology, fiber-optic DTS was introduced to the hydrological community in the early 2000s (MacFarlane et al. 2002). After the mid-2000s, reductions in the cost of the instrumentation and increased data quality spurred a significant interest of using DTS to measure and monitor various hydrologic processes (Selker et al. 2006; Tyler et al. 2009). In fiber-optic DTS, a laser pulse of a certain duration (e.g., 10 ns) is sent down a fiber-optic cable. As the laser pulse travels along the fiber, it interacts with the fiber materials and produces backscattering. One group of backscattering light is known as Raman scattering. The intensity of the anti-Stokes Raman backscatter (wavelength shorter than that of the source laser) is dependent on the temperature of the cable, while the intensity of Stokes Raman backscatter (wavelength longer than that of the source laser) is largely insensitive to cable temperature. Therefore, the intensity ratio of anti-Stokes versus Stokes Raman backscatter can be used to estimate the temperature distribution along the length of the cable. Because the signal of Raman backscatter from each pulse is very weak, thousands of pulses per second are typically required to make a temperature measurement at each cable location.

One of the common DTS applications has been using heat as a tracer to investigate groundwater movement (e.g., Lowry et al. 2007; Leaf et al. 2012; Becker et al. 2013; Liu et al. 2013; Munn et al. 2020; Simon et al. 2021). For example, Lowry et al. (2007) used DTS-based temperatures to identify several discrete groundwater discharge zones along a wetland stream. Leaf et al. (2012) used DTS to monitor the vertical movement of heat in two open boreholes, which led to an improved understanding of flow processes in a hydrostratigraphic unit at the site. Becker et al. (2013) estimated the water infiltration rates from a surface recharge basin based on the propagation of the diurnal temperatures of the infiltrated water measured by DTS. Liu et al. (2013) found that the temperature responses to active heating of a probe in a well were consistent with the K profiles previously determined at the same location, and concluded that the temperature responses could be used to approximate the variations of groundwater flux at different depths in the well. The wrapping of the sensing cable around the probe was able to significantly increase the measurement resolution (to about 1.5 cm), which led to a much improved understanding of the hydrostratigraphic controls on groundwater flow processes at the site. Munn et al. (2020) used DTS with active heating to measure borehole flows in fractured rocks under different hydraulic conditions. Results indicated that flow measurements in fracture systems can be significantly affected by cross-flow between fractures along open boreholes.

7.4.1.2 Groundwater Flux Characterization Tool

To illustrate the potential of DTS as an approach for characterizing groundwater flux, the Groundwater Flux Characterization (GFC) tool developed by Liu et al. (2013) is assessed here in a controlled laboratory setting. The GFC tool was constructed by wrapping a DTS fiber-optic cable and resistance heating cable around a sealed hollow PVC pipe (Fig. 7.13). For groundwater flux profiling, the GFC tool is first deployed to a measurement interval in an existing well. DTS temperatures are then monitored for a period of time (30 min or more) until the thermal disturbance from tool deployment has dissipated. Once temperatures return to background conditions, heating begins by flowing an electric current through the resistance cable. The rate of groundwater flux in the surrounding aquifer is proportional to the average temperature increase during heating as expressed in Eq. 7.8:

$$\Delta T_{{{\text{ave}}}} = \frac{1}{{t_{1} - t_{0} }}\mathop \int \limits_{{t_{0} }}^{{t_{1} }} \left[ {T\left( t \right) - T_{0} } \right]{\text{d}}t,$$
(7.8)
Fig. 7.13
Two views of the G F C tools and a line graph. A highlights the wrapped heating and F O cable and 8.9 centimeters P V C. B highlights the well screen, water, heating cable, F O cable, support P V C, air, and aquifer. The line graph with a decreasing trend plots temperature versus time.

a Schematic of the GFC tool, b planar schematic view of the GFC tool in a well, and c heat-induced temperature increase at different flux rates (assuming a constant rate of thermal conduction) (after Liu et al. 2013)

where T0 is the temperature [K or °C] before heating starts at time t0; t1 is the time when heating ceases; and T(t) is the temperature at time t during heating (Fig. 7.13). The heating duration (t1 − t0) is typically between 5 and 10 h and kept constant for comparisons between measurement intervals. The larger the rate of groundwater flux, the faster the movement of heat away from the probe by groundwater advection, and the smaller the temperature increase computed by Eq. 7.8 during the heating tests.

Two key assumptions are invoked using Eq. 7.8 to estimate horizontal groundwater flux. The first is that vertical flow is negligible, as the vertical flow will cause heat to move vertically along the probe and the resulting temperature responses will be difficult to separate from the temperature responses produced by horizontal flux at different depths. This limitation can be potentially addressed by zoned heating (i.e., discrete sections of the tool are heated while the temperatures of the entire probe is monitored), which requires a significant modification of the heating component of the tool. The second assumption is that the vertical variations in the thermal conductivity of the materials near the test well are negligible, so that the temperature differences between depths are mainly a result of groundwater flux instead of thermal conduction. This assumption appears to be valid when the well is backfilled with an artificial filter pack at the time of well installation. For wells without artificial filter packs, caution is needed when using Eq. 7.8 to predict groundwater flux.

For the laboratory tests discussed here, a smaller version of the GFC probe was constructed (Knobbe et al. 2015). The wrapped PVC pipe was reduced to a length of 0.98 m, with the total lengths of fiber-optic and heating cables at 43 and 40 m, respectively. As a result, each DTS measurement (1 m of fiber-optic cable) is equivalent to a 0.02-m vertical interval on the probe. There are a couple of differences between the laboratory probe and the prototype tool developed by Liu et al. (2013). The wall thickness of the PVC pipe is smaller (about 50% less), thus reducing the thermal mass and improving the sensitivity of the probe to heating. The heating cable is placed inside the fiber-optic cable (the original tool has the heating cable outside of the fiber-optic cable), diminishing the potential of heating-induced buoyancy flow in the annular space between the GFC tool and the well screen.

7.4.1.3 GFC Laboratory Test Results

Figure 7.14 shows the setup of the laboratory sand tank for testing the GFC probe (Knobbe et al. 2015). The inner dimensions of the rectangular box (steel container) are 1.83 m by 1.14 m by 1.14 m. Rigid foam insulation boards were installed inside the container to minimize the thermal interactions between the sand tank aquifer and the ambient surroundings. Commercial medium grade sand (Quikrete, Medium No. 1962, #20 - #50, 0.8–0.3 mm) was used for creating the synthetic aquifer. A 1.3 m long 10.2 cm inner diameter test well (schedule 40 PVC with screen slot width 0.025 cm) was installed at the center of the box. Two reservoirs, which were constructed using perforated PVC pipes as space retainers to provide additional support to the screened reservoir walls, were used to establish the flow field in the sands. Flow was from right to left (parallel to the longest side of the box) in the sands by pumping water from the left reservoir into the right one with a peristaltic pump. The sand aquifer has an average K of 218 m/d and an effective porosity of 33%.

Fig. 7.14
A schematic of a sand tank of length, 1.83, and width, 1.14 meters. It is divided into 2 parts and has reservoirs on both sides. A left arrow in the left reservoir is labeled flow.

Lab sand tank setup for testing the groundwater flux characterization probe. The diagram shows the cross-section of the tank along the flow direction (flow in the sands created by pumping water directly from the left reservoir to the right reservoir). The width of the tank (perpendicular to the cross-section) is 1.14 m. Perforated PVC pipes are used as space retainers in the reservoirs to provide additional support to the screened reservoir walls against the sands. Black color represents the impermeable, thermally-insulating foam boards that minimize the thermal interaction between the sand aquifer and room surroundings

Figure 7.15 shows the results of GFC tests in the sand tank at different flow rates. When the flow rate was zero, the DTS-measured temperature increase during heating was the largest at different depths. The average temperature increase over the heated section was 1.68 °C. As the flow rate increased, the temperature increase became smaller. The average temperature increase reduced to 1 °C when the Darcy’s flow rate increased to 0.78 m/d. A repeat heating test was performed at a flow rate of 0.78 m/d for quality assurance, and the temperature responses were nearly identical to those from the original test. For the tested range of flow rates (0–0.78 m/d), the relation between the average temperature increase and flow rate appears to be approximately linear (Fig. 7.15b).

Fig. 7.15
Two line graphs labeled a and b, plot the depth below the top of the probe versus T increase and T increase versus Darcys velocity. A plots for no flow, 0.26, 0.52, 0.78, and 0.78 m over d and has an increasing trend with fluctuations, and b has a decreasing trend.

The GFC temperatures from lab sand tank tests: a Time-averaged temperature increase during GFC heating at different flow velocities, and b temperature increase averaged along the length of the probe. The temperature data between depths of 0.18 and 0.63 m in (a) are used to compute the depth-averaged temperature in (b); small portions of the upper and lower heated sections are not included in the average temperature calculation due to vertical boundary impacts. In a, the temperature increase was larger at shallow depths; this might be an indication of impacts from heat-induced buoyancy or vertical variation in formation thermal conductivity (the more compacted sands at deeper locations likely have higher thermal conductivity)

The GFC laboratory test results show the promise of the tool for characterizing vertical variations in groundwater velocity. Because of the significant impact of thermal conduction, the relationship between temperature response and groundwater velocity is expected to vary in different formations or wells with different construction specifications. Further tests of the approach under different field settings are needed before it can be widely applied as a velocity characterization tool. Heating-induced buoyancy effects should always be taken into account, especially when the annular space between the tool and well screen is large (e.g., larger than several millimeters).

7.4.2 Characterization of Groundwater Velocity by Point Velocity Probe

7.4.2.1 Point Velocity Probe

The point velocity probe (PVP) uses localized tracer tests to measure the magnitude and direction of groundwater velocity at the centimeter scale (Labaky et al. 2007). A PVP operates by injecting a small amount of tracer solution into the formation (typically less than 1 mL) and measuring the tracer movement around the probe surface and across two or more detectors (Fig. 7.16). As the detectors on a PVP are a pair of parallel electrical wires that measure EC, the injected tracer needs to have a significantly different EC from the formation fluid. The most common tracer used in PVP tests is a dilute saline solution (e.g., NaCL with concentration < 1 g/L), which provides an electrical conductivity signal well above background values in freshwater systems (Devlin 2020). In some cases where the formation fluid has high salinity, deionized water can be used as the tracer, and the decrease in the measured EC is used to quantify tracer breakthrough at the detectors.

Fig. 7.16
A photograph and a schematic of the point velocity probe are labeled a and b, respectively. A highlights the injection ports and detectors, and B has a right arrow that points to the probe.

Overview of the PVP: a Photo of a multilevel arrangement of probes showing injection ports and detectors, and b illustration of tracer movement around the probe surface during a test (Ozark Underground Laboratory 2021). Injection ports can be mounted in up to three locations on the probe surface 120° apart, to ensure detection of tracer solution no matter how the probe is oriented with respect to the ambient groundwater flow field (Gibson and Devlin 2018). Each detector consists of two parallel electric wires. The measured EC between the wires provides information on tracer solution breakthrough. In b, the blue-shaded arrow indicates the general groundwater flow before the probe is inserted into the formation; the dashed lines indicate the altered groundwater flow lines produced by probe installation. The red color indicates the movement of injected tracer around the probe surface. The diameter of the PVP can vary from 2.5 to 15 cm, depending on the application

Both the direction and magnitude of the ambient groundwater velocity at the PVP measurement interval can be determined from the tracer breakthrough data at the two detectors, d1 and d2, using Eqs. 7.9 and 7.10 (Fig. 7.17),

$$\alpha = \tan^{ - 1} \left( {\frac{{v_{1} \gamma_{1} \left( { \cos \gamma_{2} - 1} \right) + v_{2} \gamma_{2} \left (1- {\cos \gamma_{1} } \right) }}{{v_{1} \gamma_{1}\sin \gamma_{2} - v_{2} \gamma_{2} \sin \gamma_{1} }}} \right),$$
(7.9)
$$v_{{\text{g}}} = \frac{{v_{1} \gamma_{1} }}{{2\left( {\cos \alpha - \cos \left( {\alpha + \gamma_{1} } \right)} \right)}}$$
(7.10)
Fig. 7.17
A schematic of the tracker breakthrough highlights the position of the injection port and detectors. A right arrow is labeled flow direction. On the right is a line graph that plots E C versus time since injection for detectors 1 and 2. The graph has a decreasing trend.

Typical tracer breakthrough during a PVP test. The diagram on the left shows the position of injection port (i), detector 1 (d1), and detector 2 (d2) relative to the ambient groundwater flow: α is the angle between ambient flow and i, γ1 is the angle between i and d1, and γ2 is the angle between i and d2

where vw is the magnitude of the ambient groundwater seepage velocity prior to the installation of the PVP, α is the angle between ambient flow and the injection port i, γ1 is the angle between i and d1, γ2 is the angle between i and d2, and v1 and v2 are the apparent velocities determined from tracer breakthrough at d1 and d2. Note that after the installation of PVP, the flow field will be altered in the immediate vicinity of the probe as groundwater has to move around the probe surface. As a result, the tracer breakthrough responses from PVP detectors measure apparent velocities on the probe surface, not the ambient velocity itself. The apparent velocities at the detectors are computed from the breakthrough data as in Eq. 7.11 (Fig. 7.17),

$$v_{1} = \frac{{r \gamma_{1} }}{{\Delta t_{1} }}, v_{2} = \frac{{r \gamma_{2} }}{{\Delta t_{2} }},$$
(7.11)

where r is the radius of the probe (m), and Δt1 and Δt2 are the durations between the start of injection and arrival of tracer concentration peaks at detectors 1 and 2 (s).

The ambient velocity calculated using Eqs. 7.97.11 is the horizontal component of groundwater flow. If there is a vertical flow component, it can, in principle, be monitored by the breakthrough data at the vertical detectors. Limited experimental data indicate the vertical detectors function as designed. However, further work is needed to fully assess the accuracy of vertical velocity measurements.

Since the PVP was first introduced into the groundwater community, early applications of the approach have focused on direct measurement of groundwater velocity by installing the probe in dedicated boreholes completed in unconsolidated, non-cohesive sediments, where the sediments can collapse and be in direct contact with the probes. This practice eliminates biases from well construction issues on the measured flow. The dedicated borehole PVP has been used for characterizing the transience of a flow field at a site undergoing bioremediation of hydrocarbons (Schillig et al. 2011, 2016), estimating contaminant mass discharge across streambanks and streambeds (Rønde et al. 2017; Cremeans et al. 2018), and measuring groundwater flow in horizontal wells for passive in situ remediation (Cormican et al. 2021).

During recent years, the PVP has also been modified for measuring groundwater flow in wells (Osorno et al. 2018, 2022) and at the groundwater-surface water interface (GWSWI) (Cremeans and Devlin 2017). In the former case, due to the significant impacts of well construction, as well as the impacts of the PVP on annular flow, empirical relationships, in the form of calibration curves, are relied upon to convert the in-well velocity measurements to ambient groundwater velocity under different well and PVP construction parameters. The in-well PVP has been used for collecting hundreds of groundwater velocity measurements across an alluvial aquifer to identify contaminant preferential flow paths (sand and gravel lenses distributed discontinuously in tight clays), and characterizing groundwater flows in fractured rocks (Ozark Underground Laboratory 2021; Heyer et al. 2021).

In the case of GWSWI studies, the PVP was adapted by miniaturizing the probe and equipping it with a hyporheic shield to isolate vertical flow. Using this technology, Cremeans et al. (2018) mapped a stream bed to identify localized zones of high discharge, facilitating the characterization of contaminant discharge zones associated with a plume of chlorinated solvents. The adapted probe was also used to measure groundwater discharge at the base of a small lake near Bemidji, MN (French et al. 2021). Preliminary results suggested that there was rapid, upward flow into the base of a thick muck layer; the flow then followed a path along the lake bottom without much mixing with the bulk of the lake water body due to the muck layer acting as a flow barrier.

7.4.2.2 PVP Case Study

Many in situ remediation technologies rely on biological processes to break down contaminants, during which aquifer bioclogging can occur as a result of biomass accumulation in the pore space. Despite the wide recognition of its significance, few studies have investigated bioclogging and its impact on groundwater flow under field conditions. This is mostly due to the difficulty of making repeated K or groundwater velocity measurements at a scale that is sufficiently small for assessing aquifer property changes caused by biological processes. The development of the PVP provides an effective means of addressing this concern.

Schillig et al. (2011) presented a PVP case study for investigating the transience of groundwater flow in a sand aquifer undergoing benzene, toluene, ethylbenzene, and xylene bioremediation. Figure 7.18 shows the site setting. The unconfined sand aquifer under study was isolated by a clay aquitard on the bottom and two sheet pilings on the sides. Hydrocarbon was released into the aquifer between depths 2.5 and 4.0 m in the source zone. Five dedicated multilevel PVP stands, each equipped with four probes, were installed along a transect 13 m downstream of the hydrocarbon source. Five fully screened remediation wells were used to administer dissolved oxygen, using Oxygen Release Compound (ORC®) 1 m upstream of the PVP transect. Oxygen was released to the aquifer for stimulating aerobic biodegradation of the dissolved hydrocarbons. The ORC® releases oxygen with diminishing strength over time, so the wells were recharged with additional oxygen three times during the experiment in September 2005, February 2006, and September 2006.

Fig. 7.18
A schematic of an aquifer consists of a clay aquifer 7 meters long. In the aquitard are a petroleum source and an up-arrow labeled average ambient flow direction. In the front are P V P multilevel designations with the probe numbers.

Study site setting at Borden, Canada by Schillig et al. (2011). The aquifer is comprised of 7 m of well-sorted fine- to medium-grained sands bounded by sheet piling on the east and west and a clay aquitard underlying the sand. The depth to the water table fluctuates with precipitation and is generally within 1 m of ground surface. Groundwater generally flowed from the south to the north, with a hydrocarbon source 13 m upstream of a fence of PVP probes. Each PVP borehole had 4 measurement intervals with a vertical spacing of 0.8 m. Five fully screened remediation wells (not shown here) with dissolved oxygen release compound were installed 1 m upstream of the PVP fence to stimulate aerobic biodegradation of the released hydrocarbon. The objective of the study was to document changes to the groundwater velocity in response to microbial growth and activity during the aerobic biodegradation of the hydrocarbons

Figure 7.19 shows the groundwater velocity changes at the PVP transect at the different times; August 2005 represents the pre-oxygen addition background condition for the experiment, so the differences plotted are zero (Fig. 7.19a). The plotted values are interpolated from the velocity changes at the 20 measurement locations (Fig. 7.18). After the first addition of oxygen in September 2005, groundwater velocity showed a clear decrease across the transect, with the largest reduction occurring near the lower right corner of the plot (Fig. 7.19b). Flow directions (not shown) remained largely similar to the background field, except in the lower right area where the direction of flow changes from northeast to northwest. Figure 7.19c shows the measured velocity three months after the second addition of oxygen in February 2006. Compared to October 2005, groundwater velocity in May 2006 increased significantly—possibly the result of seasonal changes of flow in the aquifer—and flow directions changed more significantly across a large portion of the transect. Most notably, the largest changes from background were found to occur in the locations coinciding with the highest concentrations of hydrocarbon, where bioactivity might be expected to be maximized. Following the May measurements, the ORC® was permitted to deplete to exhaustion, allowing the site to return to ambient conditions. The resilience of the system was tested with a final addition of oxygen in September 2006. One month later, in October 2006 (Fig. 7.19d), the measured groundwater appeared to undergo a large decrease in velocity in a fashion resembling the response a year before. Samples of aquifer material recovered in core and examined for microbial biomass near the PVP transect, and at a separate location removed from the biostimulation, showed significant increases in microbial numbers in the biostimulated zone, lending support to the notion that biomass accumulation might have been responsible for at least some of the observed velocity transience.

Fig. 7.19
Four heat maps labeled a to d, plot elevation, above the aquifer bottom versus distance from left sheet piling for 5 August and 5 October of 2005, and 6 May and 6 October of 2006, respectively.

Changes in velocity at the PVP transect at different times before and after adding oxygen to the hydrocarbon contaminated aquifer: a August 2005, difference from background is zero, b October 2005, following O2 addition in September 2005, c May 2006, after established biodegradation, and d October 2006 re-addition of O2 following a 2–3 month waning of the O2 source (after Schillig et al. 2011). For comparison, the hydrocarbon concentrations measured in October 2006 (based on water samples from 36 multilevel sampling wells across the site) are plotted as contours on the maps (the two centers of the contours at elevation 3.6 m represent peak concentrations). The changes in groundwater velocity were believed to be at least partially caused by biological growth following the additions of oxygen at different time. Note the coincidence of the greatest changes in velocity with the most concentrated portion of the hydrocarbon plume where biodegradation was expected to be the most significant

It should be pointed out that the plotted values in Fig. 7.19 provide point assessments of velocity responses, and the interpolations presented should not be considered true representations of the detailed velocity distribution. To obtain such a picture of the aquifer structure and velocity distribution, the density of sampling would have to be increased considerably—a testimonial to the challenges in aquifer characterization. Further effort is needed to determine the density of sampling needed to achieve practical uses for such data, for example the estimation of water flow across the transect, or contaminant mass discharges and the associated uncertainty. Nonetheless, the sampling density adopted for this project was sufficient to identify and quantify possible hydraulic responses to biostimulation.

7.5 Summary and Conclusions

In this chapter, we discussed high-resolution characterization of the physical properties in shallow, saturated, unconsolidated subsurface conditions using DP technology, NMR logging, DTS, and PVPs. These developments have led to a significantly improved ability to obtain information about subsurface properties (K and porosity) and groundwater velocity at a speed and resolution that has not previously been possible. These advances have allowed the discipline to gain important insights into the fundamental controls on groundwater flow and transport processes under field conditions (Schillig et al. 2011; Bohling et al. 2012; Fiori et al. 2013; Liu et al. 2013; Dogan et al. 2014; Knight et al. 2016; Munn et al. 2020).

The performance of these approaches was demonstrated using examples from field and laboratory settings. DP EC profiling could potentially provide a good indicator of relative K when electrically conductive clay is a dominant control on permeability, although more work is needed to explore this use of EC. The DPP is a tool that can provide reliable K estimates in moderately to highly permeable zones. It has several major advantages over conventional hydraulic testing approaches (such as slug and flowmeter tests) due to the steady-shape requirement for the analysis. In addition, the DPP is not impacted by low-K skins, which are known to be a major challenge for slug tests. The practical time constraints, however, limit the DPP vertical resolution to 0.4 m or larger under most field conditions. The continuous DPIL is a high speed and resolution (0.015 m) approach, but only characterizes relative variations in K. A general empirical model may be used to estimate K values within a narrow, but useful, K range (McCall and Christy 2010). Alternately, regressions with additional data, such as slug test K values from nearby locations, are used to convert the DPIL data into site-specific K estimates. The HRK tool was designed to exploit the advantages of the DPP and DPIL by coupling them into a single probe. Because DPP tests are collocated with DPIL logs, they can be used to directly transform the DPIL data into K without the need to compare measurements at different support scales as in the regression-based approaches.

In field investigations, the DPIL and HRK can be applied in a complementary fashion (Bohling et al. 2012). Because the continuous DPIL is rapid, many DPIL profiles can be performed across a site in a short period. The HRK tool, on the other hand, takes more time (due to halting tool advancement for DPP tests) and can be applied at far fewer locations. By coupling the HRK tool with DPIL, the DPIL K power-law relation determined from the HRK profiles can be applied to all the DPIL logs, resulting in a high-resolution characterization of K for the entire site at only a fraction of time and costs of other approaches (Bohling et al. 2012).

The current DP tools (including the DPP, DPIL, and HRK tool) do not perform well in low-K settings. The DPP tests take a long time to complete when K is low (e.g., < 0.001 m/d). For DPIL, it remains a challenge to inject water at the small rate (e.g., 1 mL/min) that is needed in low-K formations to avoid high injection pressure and formation alteration. The pressure generated from tool advancement can also significantly impact the measured injection pressure signal. In addition to the low-K limit, the DPIL also has an upper limit as the injection-induced pressure responses become very small in highly permeable zones. Future work is needed to improve the range of reliable K estimates by the DPP, DPIL, and HRK tool.

Logging with NMR technology, a widely used borehole technique in the petroleum industry, has been adapted to various applications in near-surface hydrology (Walsh et al. 2011, 2013). Compared to the DP approaches, a major advantage of NMR is that it provides an estimate of porosity in addition to K. Because the NMR signal is directly a function of the total number of water molecules in the measurement zone, the accuracy of estimated porosity is much higher than that of K. Given the challenges of DP approaches in low-K settings, NMR logging holds great potential for characterizing K and porosity in low-K formations.

There has been a growing interest in measuring groundwater flux directly in the field as the remediation community switches from contaminant concentration-based decision-making to one based on contaminant mass discharge (Suthersan et al. 2010; Devlin 2020). In recent years, significant progress has been made in using heat and other tracers to measure groundwater flux. Liu et al. (2013) developed a GFC tool by wrapping a DTS fiber-optic cable and resistance heating cable around a sealed hollow PVC pipe. The GFC tool was assessed in a series of laboratory tests in a sand tank, and results showed a linear relationship between the heating-induced temperature increase and the ambient groundwater flux over the range of 0–0.78 m/d (Darcy’s flux). One key assumption in using the GFC is that the vertical variations in the thermal conductivity of the materials near the test well are negligible, so that the impacts of thermal conduction can be ignored. However, when thermal conduction is important, the detailed characteristics of temperature change with time may be used to quantify formation thermal conductivity in addition to groundwater flux (Simon et al. 2021).

The PVP uses localized tracer tests to measure the magnitude and direction of groundwater velocity at the centimeter scale (Labaky et al. 2007). It operates by injecting a small amount of tracer solution into the formation and measuring the movement of tracer around the probe surface via multiple detectors. Both the direction and magnitude of the ambient groundwater velocity can be determined from the tracer breakthrough data at the two detectors. Early applications of the approach focused on direct measurement of groundwater velocity by installing the probe in dedicated boreholes, while recent adaption of the tool allows it to be used for flow measurement in wells, at the groundwater-surface water interface, and in horizontal wells (Ozark Underground Laboratory 2021). The PVP has been used in a variety of site characterization projects, and its high-resolution measurement can provide improved understanding of groundwater flow under transient conditions.

The approaches discussed in this chapter were primarily developed for use in intervals saturated with groundwater. When non-aqueous phase liquids (NAPLs) are present, additional research is needed to investigate the performance of these approaches as NAPL-water movement will not only depend on the properties of the NAPL and groundwater but also on their interactions. One of the most common approaches for characterizing the transmissivity of LNAPL (NAPL lighter than water) intervals is slug testing (Charbeneau et al. 2016; Butler 2019). However, this only provides an estimate on the mobility of LNAPLs as a free phase in the immediate vicinity of the test well. To characterize LNAPL movement throughout the groundwater system, detailed information about relative permeabilities, NAPL saturations and their spatial distribution is needed. Furthermore, a practical challenge for characterizing sites with LNAPL is the potential risk of vertically spreading the contaminants across different geologic units by well construction (Newell et al. 1995). This challenge can be largely overcome in shallow unconsolidated sediments with DP-based approaches such as the Membrane Interface Probe (MIP; ASTM 2018) and laser-induced fluorescence (LIF) probes presented in Chap. 8. As discussed in Chap. 9, NMR logging could be used to map the distribution of NAPL because of the abundance of hydrogen in petroleum NAPL, although the presence of NAPL can affect the utility of NMR for characterizing porosity and K. For PVP-based studies, preliminary testing has demonstrated that velocity measurements are sensitive to the presence of gases in the pore space (Cormican et al. 2021), which is analogous to the NAPL case since tracer responses are mainly controlled by pore water flow that is restricted when NAPLs are present. Additional research is needed to extend the preliminary work and investigate how the presence of NAPL affects the measured velocity under different subsurface conditions.

There has been an increasing recognition of the need for “focused” site remediation with remedial measures targeted at the processes and locations that are crucial to remediation success (ITRC 2004). Identifying these critical processes and their locations, however, remains a difficult challenge at many sites as the conventional approaches (e.g., the well-based pumping and slug tests) often fail to provide information at sufficient details to make those identifications. The high-resolution approaches discussed in this chapter hold great promise for addressing this challenge. Each of these approaches has its own advantages and limitations, and future studies are needed to refine their performance as discussed above. Finally, when site conditions are complex, a combination of different approaches may provide the best solution for site investigators to obtain the high-resolution data that are needed to develop and implement a successful site remediation program.