Bedrock depth influences spatial patterns of summer baseflow, temperature and flow disconnection for mountainous headwater streams

In mountain headwater streams, the quality and resilience of summer cold-water habitat is generally regulated by stream discharge, longitudinal stream channel connectivity and groundwater exchange. These critical hydrologic processes are thought to be influenced by the stream corridor bedrock contact depth (sediment thickness), a parameter often inferred from sparse hillslope borehole information, piezometer refusal and remotely sensed data. To investigate how local bedrock depth might control summer stream temperature and channel disconnection (dewatering) patterns, we measured stream corridor bedrock depth by collecting and interpreting 191 passive seismic datasets along eight headwater streams in Shenandoah National Park (Virginia, USA). In addition, we used multi-year stream temperature and streamflow records to calculate several baseflow-related metrics along and among the study streams. Finally, comprehensive visual surveys of stream channel dewatering were conducted in 2016, 2019 and 2021 during summer low flow conditions (124 total km of stream length). We found that measured bedrock depths along the study streams were not well-characterized by soils maps or an existing global-scale geologic dataset where the latter overpredicted measured depths by 12.2 m (mean) or approximately four times the average bedrock depth of 2.9 m. Half of the eight study stream corridors had an average bedrock depth of less than 2 m. Of the eight study streams, Staunton River had the deepest average bedrock depth (3.4 m), the coldest summer temperature profiles and substantially higher summer baseflow indices compared to the other study steams. Staunton River also exhibited paired air and water annual temperature signals suggesting deeper groundwater influence, and the stream channel did not dewater in lower sections during any baseflow survey. In contrast, Paine Run and Piney River did show pronounced, patchy channel dewatering, with Paine Run having dozens of discrete dry channel sections ranging from 1 to greater than 300 m in length. Stream dewatering patterns were apparently influenced by a combination of discrete deep bedrock (20+ m) features and more subtle sediment thickness variation (1–4 m) depending on local stream valley hydrogeology. In combination, these unique datasets show the first large-scale empirical support for existing conceptual models of headwater stream disconnection based on spatially variable underflow capacity and shallow groundwater supply.

1 Introduction

Mountain headwater stream habitat is affected by hydrologic connectivity along the surface channel, and connectivity between the channel and multiscale groundwater flowpaths (Wohl, 2017; Covino, 2017; Fausch et al., 2002). Discharge from shallow groundwater within the critical zone is a primary component of stream baseflow, attenuating maximum summer temperatures and creating cold water patches (Singha and Navarre-Sitchler, 2021; Sullivan et al., 2021), and shaping catchment topography (Litwin et al., 2022). In headwater stream valleys characterized by irregular bedrock topography and thin, permeable sediments, nested physical processes interact to control the connectivity of groundwater/surface water exchange (Tonina and Buffington, 2009). Between stormflow, fast soil drainage and snowmelt events, headwater streamflow (baseflow) is primarily generated by groundwater discharge due to a relative lack of soil water storage and release (Winter et al., 1998). Unlike in lower valley settings, mountain headwaters accumulate less fine soil, facilitating efficient routing of quickflow to streams through macropores and other preferential flowpaths within regolith and saprolite (Sidle et al., 2000). Recharge that does percolate vertically contributes to shallow groundwater along steep hillslopes and valley floors where groundwater flowpath depths are constrained by bedrock topography (Buttle et al., 2004). Although deeper groundwater may also represent an important contribution to summer streamflow in systems with relatively permeable bedrock (Burns et al., 1998; O'Sullivan et al., 2020), shallow, low permeability bedrock generally restricts stream–groundwater connectivity to the thin layers of unconsolidated sediments (Briggs et al., 2018b).

In addition to baseflow drainage along headwater stream networks, down-valley shallow groundwater “underflow” can be substantial when high-gradient streams lack sinuosity and flow over permeable sediment (Figs. 1a and A1). In fact, headwater stream channels may only be expected to show surface flow when the transmission ability of the underlying alluvium and colluvium is exceeded, and bedrock depth is thought to be a primary control of this underflow capacity (Ward et al., 2020). In some hydrogeologic settings, underflow can dominate groundwater export from mountain catchments compared to groundwater drainage via the surficial stream channel (Larkin and Sharp, 1992; Tiwari et al., 2017). Moreover, in addition to longitudinal transport down-valley, underflow also acts as a reservoir of exchange for hyporheic flowpaths that may mix with shallow groundwater before returning to channel flow (Payn et al., 2009), transporting buffered temperature signals back to channel waters (Wu et al., 2020). Local underflow is recharged from upgradient flowpaths and adjacent hillslopes, creating complex seasonal and interannual patterns in groundwater connectivity and discharge to surface water (Jencso et al., 2010; Johnson et al., 2017). A major challenge to understanding groundwater exchange in headwaters is that attributes of the streambed subsurface, such as the depth to the underlying bedrock contact, are often only available from limited direct measurements, coarse spatial interpolations or inferred remotely based on landscape forms. Therefore, methods that allow efficient, local measurements of the streambed subsurface are critically needed.

Figure 1 A conceptual mountain stream valley cross-section (a) and longitudinal profile (b) indicating the expected control of low permeability bedrock topography on groundwater temperature, stream–groundwater exchange, patchy stream dewatering and the underflow reservoir.

Seasonal thermal regimes of mountain headwater streams can be profoundly impacted by groundwater inflow from multiple depths (Briggs et al., 2018a). In lower valley settings, the temperature of groundwater discharge along stream networks is often assumed to be constant throughout the year and approximately equal to the average annual land surface temperature (Stonestrom and Constantz, 2003). Conversely, shallow groundwater temperature (within several meters from land surface) can show pronounced seasonality (Bundschuh, 1993; Lapham, 1989) and high spatial variability, even over small spatial extents (Snyder et al., 2015). The warming of shallow groundwater during the summer and fall seasons can limit the ability of gaining mountain streams to support cold-water fish populations during the low flow season even if baseflow (assumed to be dominated by groundwater discharge) fractions are large (Johnson et al., 2020). In systems with low permeability bedrock, thicker hillslope sediments may generate deeper, colder lateral groundwater flow to streams in summer (Fig. 1a), increasing cold water habitat resiliency (Briggs et al., 2018b). For example, a recent meta-analysis of stream and air temperature records across the contiguous United States found that a substantial fraction of shallow-groundwater-dominated streams displayed summer warming trends in recent decades, while deeper-groundwater-dominated streams were more stable (Hare et al., 2021). Steep mountain stream systems such as those found in the Blue Ridge and Cascade mountains of the USA have been found to show annual thermal regimes indicative of shallow groundwater (Johnson et al., 2020), indicating such streams may also be at risk for warming over time, contrary to assumptions based on (relatively high) elevation alone.

Beyond warm summer stream temperatures, the dewatering and disconnection of the active stream channel during summer low flows can adversely impact fish habitat by impeding fish movement (Edge et al., 2017; Labbe and Fausch, 2000; Rolls et al., 2012; Snyder et al., 2013), locally degrading water quality (Hopper et al., 2020) and increasing predation risks in isolated pools (Magoulick and Kobza, 2003). However, the physical controls on localized stream channel dewatering are not well characterized and likely involve a spectrum of nested gaining and losing flowpaths. For mountain headwater streams, previous research has documented major contractions of drainage networks during seasonal dry down (Ilja Van Meerveld et al., 2019) and general seasonal shifts in hydraulic gradients from gaining to losing, with closely coupled increased streamflow and precipitation events, indicating a dominance of shallow routing rather than deeper groundwater connectivity in maintaining streamflow (Zimmer and McGlynn, 2017). Warix et al. (2021) found that although deeper/older groundwater was found to contribute to their study streams during dry down, those sources were insufficient in preventing dewatered channel sections from occurring, also indicating the importance of shallow groundwater inflows and local geologic controls. Locally-losing sections of headwater stream channels can be associated with coarse, permeable colluvial deposits from hillslope mass wasting processes (Weekes et al., 2015; Costigan et al., 2016) as local enhancement of the total pore space under mountain streams can drive downwelling of streamwater (Fig. 1b, Tonina and Buffington, 2009). Main channel dewatering occurs when the bed sediments have a storage and transport capacity that exceeds stream discharge (Rolls et al., 2012; Ward et al., 2018), though streamwater losses can also be driven by local changes in bed morphology and slope (Costigan et al., 2016) and bedrock permeability. A downstream shallowing of the underlying bedrock contact may drive lateral underflow toward the surface causing the channel to gain water (Herzog et al., 2019, Fig. 1b), though such hypothesized dynamics are not well documented in existing literature due to a relative lack of bedrock topography data beneath and adjacent to headwater streams.

At large scales, contiguous bedrock depth layers are interpolated from a combination of relatively sparse borehole data and surface topography (Kauffman et al., 2018; Pelletier et al., 2016; Shangguan et al., 2017). However, in steep headwater systems with little borehole data, bedrock topography is difficult to predict accurately from land surface topography alone. The development of improved tools for predicting bedrock depth is an active area of research which has recently demonstrated promise when bedrock outcrop data are included (e.g., Furze et al., 2021; Odom et al., 2021). The limitations of using landform data to predict bedrock depth are compounded by inherent challenges in collecting physical data via soil pits and monitoring wells in rugged, rocky terrain, and so, direct measurement data are often limited to highly studied experimental watersheds where bedrock depth is still only inferred from piezometer installation refusal (e.g., Jencso et al., 2010; Ward et al., 2018). In more typical headwater systems, existing wells may be preferentially installed to maximize the production of water and not broadly sample the true range of bedrock depths.

Application of near-surface geophysical methods to stream corridor research has increased appreciably in recent years (McLachlan et al., 2017), and several methods are sensitive to shallow subsurface flow and geologic attributes including bedrock depth. Active seismic refraction measurements can provide high resolution (tens of centimeters) bedrock depth information along transect-based cross-sections (e.g., Flinchum et al., 2018) but are less suited for exploration throughout rugged mountain stream valleys at the many km-scale. This is due to logistical challenges in using active seismic methods to obtain a sufficient amount of data to effectively characterize important variation in bedrock depth at relatively small, ecologically relevant spatial scales.

Point-based, efficient passive seismic measurements represent a unique combination of high mobility and relative precision for measuring bedrock depth along mountain valleys. The horizontal-to-vertical spectral ratio (HVSR) method is a passive seismic technique that evaluates ambient seismic noise recorded using handheld instruments placed on the ground surface to identify seismic resonance that develops due to strong vertical changes in subsurface acoustic impedance (Yanamaka et al., 1994). While typically insensitive to variations in unconsolidated sediment permeability (i.e., clay lenses), the HVSR method is effective at identifying the depth to distinct unconsolidated sediment/bedrock interfaces at essentially the “point” spatial-scale. HVSR measurements are often not successful in settings with highly weathered bedrock surfaces such as those with pronounced epikarst and saprolite.

The control of stream to groundwater exchange (i.e., “transmission losses”) on streamflow permanence has been highlighted as an important research need by the comprehensive review of intermittent stream systems by Costigan et al. (2016). Following the conceptual model of Ward et al. (2018) for mountain stream corridors, a central hypothesis of our research was that bedrock depth along the stream corridor will act as a first-order control on stream dewatering patterns when shallow bedrock is of low permeability and streambed sediments of high permeability. Based on the concepts presented by Tonina and Buffington (2009), we postulated that relatively thick, permeable surficial sediment zones could locally accommodate the entirety of low streamflow volumes, dewatering main channel sections at varied scales when not balanced by groundwater inflow (Fig. 1b). We further hypothesized that summer stream channel thermal regimes would also be influenced by bedrock depth, as the temperature of hillslope groundwater and underflow is, in part, depth-dependent, indicated conceptually in Fig. 1a. To test our hypotheses, we extended the existing mountain headwater bedrock depth surveys from Shenandoah National Park (SNP), Virginia, USA (Briggs et al., 2018b), to seven additional subwatersheds and compared results to physical mapping of stream dewatering, multi-year stream temperature data and derived groundwater influence metrics, and baseflow separation analysis to address the following research questions.

1. Does stream corridor bedrock depth exhibit longitudinal spatial structure in mountainous streams at ecologically relevant spatial scales? Can measured bedrock depth dynamics be accurately extracted from existing large-scale datasets or inferred from high resolution soils maps?

2. Does underflow generally represent a net source or sink of summer flow for headwater streams based on observed dewatering patterns and groundwater influence metrics?

3. Does bedrock depth explain spatial variation in stream temperature and summer baseflow indices within headwater streams?

2 Study area

SNP is an 800 km² area of preserved headwater forest perched along a major ridgeline of the Blue Ridge Mountains in northern VA, USA (Fig. 2). The bedrock of the park is predominantly low permeability basaltic and granitic material in the central and northern sections, and siliciclastic along the southern section (Southworth et al., 2009), though many subwatersheds also transition in dominant bedrock type from high to lower elevation. Stream valleys of SNP are typically steep and feature a perennial channel with mainly non-perennial tributaries (Johnson et al., 2017, Fig. A1), and stream baseflow consists of less than three-year old groundwater on average (Plummer et al., 2001). In contrast, water collected from SNP hillslope wells completed in shallow fractured rock generally have higher ages of 10–20 year (Plummer et al., 2001), indicating minimal contributions from bedrock groundwater to streamflow. Previous ecohydrological research in SNP has noted that some mainstem stream channels show patchy dewatering at summer low flows (Snyder et al., 2013), though the physical controls on these patterns of stream drying were not clear.

Figure 2 This study was based in Shenandoah National Park (a) located in the Blue Ridge Mountains of northeast USA (b). Lidar hillshade cutouts of each subwatershed illustrate the rugged terrain and varied valley morphology (c). The mainstem stream channel and tributaries are traced, and passive seismic sample measurement locations noted.

In SNP, stream baseflow is thought to be predominantly generated by near-surface drainage of coarse unconsolidated alluvium and colluvium (DeKay, 1972; Nelms and Moberg, 2010). The mountain ridgeline streamflow systems are expected to drain near-surface flowpaths and accommodate substantial down-valley underflow below perennial stream channels (Fig. A1). A portion of hillslope recharge is expected to percolate downward through connected bedrock fractures into the deeper groundwater reservoir contributing to mountain block recharge along the Shenandoah River Valley. Narrow alluvium deposits mapped along the stream corridors of SNP are thought to generally range up to 6 m in thickness and be more clay rich when sourced by basaltic bedrock (Southworth et al., 2009). Data at sparse wells drilled along the SNP ridgeline indicate bedrock depth can range to over 20 m on hillslopes and be highly variable (Goodling et al., 2020; DeKay, 1972; Lynch, 1987).

Previous research has inferred summer and annual groundwater discharge patterns throughout SNP subwatersheds based on paired, local air and stream water temperature dynamics (Briggs et al., 2018a; Johnson et al., 2017; Snyder et al., 2015). Combined, these analyses indicated stream–groundwater exchange is highly variable in space along singular stream valleys and between subwatersheds, and dependent upon local- to subwatershed-scale characteristics. A combination of landform features that include stream slope and stream valley confinement operate in conjunction with seasonal precipitation to drive groundwater influence on summer stream temperatures (Johnson et al., 2017). Multi-week lags in time between streamwater and local air annual temperature signals (i.e., water signal phase shifts toward later time) were observed from dozens of the 120 total monitored stream sites, indicating a dominance of shallow groundwater discharge originating generally within approximately 3 m of land surface (Briggs et al., 2018a).

3 Methods

3.1 Passive seismic bedrock depth measurements

Periodically from the summer of 2016 to the spring of 2020, we acquired 323 HVSR measurements across SNP. The geophysical data were collected along the perennial streams of seven subwatersheds with extensive existing stream temperature and ecological datasets, and at known ridgeline and hillslope borehole locations. This effort added to previously interpreted HVSR data from 22 riparian sites collected along the White Oak Canyon subwatershed in late 2015 (Briggs et al., 2017) for a total of eight mountain streams for analysis in this study (Fig. 2). In July 2016, HVSR data were collected in the following subwatersheds: Piney River, Paine Run, Meadow Run, Jeremys Run, Hazel River and Hughes River. Some stream sections were inaccessible due to steep bedrock walls and waterfalls, resulting in poor data coverage in those areas.

Measurement locations mostly coincided with existing stream temperature monitoring stations (described by Snyder et al., 2017), and were typically made at points immediately adjacent to the stream or on larger rocks within the channel (Fig. 3). In July 2019, HVSR data were again collected along Paine Run and Piney River subwatersheds, and throughout the lower Staunton River (Fig. 3). The 2019 survey design differed in that transect measurements were made at four locations along the stream channel waterline spaced approximately 25 m apart at longitudinal locations that differed from the 2016 survey. This was done to assess potential variation in bedrock depth along short subreaches of these three streams. Finally, clustered HVSR data were collected in March 2020 in Paine Run and Piney River in zones previously observed to show channel disconnection and streamflow re-emergence. Measurement locations were chosen to test the hypothesis that the dewatering patterns were controlled locally by bedrock depth as shown conceptually in Fig. 1b.

Figure 3 Typical sections of (a) Paine Run, (b) Piney River, (c) Staunton River and (d) a section of Paine Run that was dewatered at baseflow, leaving isolated pools. The passive seismic HVSR instruments are shown deployed in panels (a), (b) and (d) (photographs by the US Geological Survey).

HVSR data were collected using multi-component Tromino seismometers (MOHO, S.R.L.) directly coupled to the land surface or placed on heavy metal plates where sediment was loose. Collection times ranged 10–20 min at either 128 or 256 Hz sampling rates. HVSR data collection locations were determined by a combination of internal Tromino GPS and external GPS units. HVSR measurements were processed to derive a resonant frequency using a commercially available program (GRILLA^® v. 8.0 (2018); further details regarding data processing are given by Goodling et al. (2020).

Resonant frequency measurements that passed a series of quality criteria were then converted to a bedrock depth estimate following Briggs et al. (2017). This conversion necessitates a shear wave velocity estimate for the unconsolidated sediments over bedrock. HVSR data collected at six spatially distributed boreholes with documented depth to varied-type bedrock along the SNP ridgeline indicated a mean shear wave velocity of 358.7±56 m s⁻¹ (Goodling et al., 2020). A similar shear wave velocity of 346 m s⁻¹ was measured at two locations along the White Oak Canyon riparian zone spaced several kilometers apart using active seismic methods (Briggs et al., 2018b). This agreement indicates a common shear wave velocity can be assumed for the unconsolidated material of SNP subwatersheds. For this study, we used the average of these spatially distributed active and passive seismic methods at 352 m s⁻¹. The average shear wave velocity calculated in this study is comparable to the mean shear wave velocity ranges in firm soils (180–360 m s⁻¹) and very dense soil and soft rock (360–760 m s⁻¹), according to National Earthquake Hazards Reduction Program (NEHRP) guidelines (Building Seismic Safety Council, 1994). As an example of measurement sensitivity to the shear wave velocity parameter for shallow bedrock contacts, a velocity change in either direction by 25 m s⁻¹ would generally shift the bedrock depth estimate by <0.2 m.

3.2 Observations of spatial dewatering patterns

Longitudinal (upstream to downstream) patterns of dewatering were determined in the summers of 2016, 2019 and 2021 during baseflow conditions over 124 total km of stream length for all surveys combined. In July–August 2016, all eight subwatersheds (Fig. 2) were surveyed. In September 2019 and August 2021, dewatering surveys were repeated in three subwatersheds (Paine Run, Piney River and Staunton River) to evaluate annual variation in dewatering patterns. Data were collected by a team of investigators walking each stream from an upstream location defined by the point along the stream draining 75 ha (assumed capture area required to generate perennial streamflow determined using watershed tools in ArcGIS) to the bottom of each watershed defined near the park boundary, and mapping transition points between three hydrologic categories: wet, dry or isolated pools, based upon investigator observation. “Wet” segments were defined as reaches where the entire channel was wet with flow between pools, “dry” segments were defined as reaches containing no water or isolated pools of insufficient depth to sustain 1+-year-old brook trout fish (Salvelinus fontinalis), and “isolated pools” were defined as reaches containing pools of sufficient depth to support brook trout but were hydrologically disconnected from other parts of the channel. An example of isolated pools is photographically depicted in Fig. 3d. Spatial coordinates of transition points were mapped using a Trimble R2 Global Navigation Satellite System (GNSS) receiver for <1 m accuracy. Surveys for each subwatershed were completed within a single day to minimize effects of temporal variation in precipitation.

In addition to local variability in bedrock depth, spatial patterns of dewatering and stream temperature are likely to be influenced by seasonal precipitation and air temperature proximate to the period of measurement (i.e., summer conditions, 2016 and 2019). We used historical weather records (1942–2020) collected from the nearby Luray Weather Station located within SNP (station no. GHCND:USC00445096) to compare weather conditions during these two study years with historical norms. Finally, 3D surface area of each subwatershed was determined from existing lidar data using the Add Surface Information tool in ArcGIS, and mean valley bottom width was evaluated from lidar data using 100 m transects measured approximately 2 m above the valley floor.

3.3 Stream channel temperature data and baseflow separation

Multi-year SNP stream temperature data were collected at hourly time intervals as described by Snyder et al. (2017) using HOBO Pro V2 thermographs (±0.2 ^∘C expected accuracy). From this larger dataset, 64 main channel locations within the eight study subwatersheds were extracted and processed for summary statistics such as the maximum and minimum of the 7 d running mean using Matlab R2019b software (Mathworks, Inc.). Only complete 7 d periods were included in the running average. Warm season data (July, August, September) were isolated and analyzed to coincide with the stream dewatering surveys and a larger body of research regarding summer cold-water brook trout habitat in SNP. We utilized stream temperature data processed to extract annual temperature signals by Briggs et al. (2018a) where dry sensor periods were identified and removed, impacting a handful of the upper stream sites. Data were visualized and downstream trends explored using Sigmaplot 14.0 software (Systat Software Inc.). Baseflow separation was conducted for the three continuously gauged streams of this study (Paine Run, Piney River, Staunton River) over summer months for the period of record (1993–2020). Following the approach of Hare et al. (2021), the daily baseflow index (BFI) was calculated using the USGS-R “DVstats” package (version 0.3.4) following methods described by Barlow et al. (2014), and dividing the calculated baseflow discharge by the corresponding stream discharge where a value of one would indicate stream discharge was entirely composed of baseflow. BFI was then averaged (mean) across each summer season, along with the mean and standard deviation of summer stream discharge.

4 Results

4.1 Stream corridor bedrock depth

Approximately 60 % of individual HVSR measurements (191 of the 323) were of high enough quality to be interpreted for bedrock depth using objective data quality metrics reported by the GRILLA software. This ratio of interpretable to total HVSR measurements was similar to the previous 2015 White Oak Canyon run study using the same instrument type (Briggs et al., 2017). For the 132 datasets that could not be interpreted, the primary reason was no identifiably resonant frequency “peak” in the multicomponent seismic data, as described in more detail in the data release of Goodling et al. (2020). The loosely consolidated, rocky surficial soils of many SNP subwatershed riparian zones likely contributed to poor instrument coupling to the land surface and, therefore, reduced measurement sensitivity/success compared to firmer soils. Additionally, bedrock depths of approximately >1 m may not be identifiable using HVSR methods based on extensive field experience with that methodology and the results of this study (Fig. 4). However, due to spatial redundancy in the measurements, the 191 locations where bedrock depth was evaluated generally covered all the intended longitudinal stream measurement locations throughout the subwatersheds.

Figure 4 Measured depth-to-rock along the stream channel and riparian zones of the eight study subwatersheds. Exposed bedrock (i.e., zero depth) observed at the intended measurement location is noted here by a value of “0.1” on the log scale. The median value is shown as a labeled horizontal line.

The median bedrock depth was smallest for Hughes River (1.52 m), and similar for Meadow Run, Jeremys Run and Hazel River (1.92, 1.94 and 1.98, respectively, Tables 1 and B1, Fig. 4). Paine Run had a median of 2.24 m, White Oak Canyon of 2.38 m and Piney River of 2.54 m. Lower Staunton River had the largest median depth-to-rock of 3.43 m (Table 1). Piney River had the largest variation in bedrock depth, including a discrete zone greater than 20 m deep, along with several zones of exposed bedrock along the channel. Visual observations of exposed channel bedrock were not incorporated into the median bedrock depth averages listed above (Fig. 4, Table 1). Simple bivariate relations were explored between the physical valley parameters presented in Table 1, and a negative relation was found between bedrock depth and mean valley bottom width (Fig. 5a) while other relations were not significant.

Table 1 The median bedrock depth along with the elevation, mean and 7 d maximum summer temperatures over the period of record collected at most downstream site location in each subwatershed.

Figure 5 Median study stream corridor bedrock depth showed a negative relation to valley bottom width (a), and mean summer stream temperature at the lower study stream boundaries was negatively related to median bedrock depth (b).

4.2 Spatial dewatering patterns and climate data

Cumulative monthly precipitation during baseflow summer (July–September) was higher than normal in 2016 and near average or lower than average (period of record 1942–2020), depending on the month, in 2019 (Fig. A2). Mean monthly air temperatures were higher than average for both study years during baseflow summer, reflecting the long-term trend of increasing air temperatures in the park (Luray weather station GHCND:USC00445096). Patches of stream dewatering were observed along five of the eight study subwatersheds between 19–27 July 2016 when over 98 km of total stream length was mapped (Fig. 6). However, for Meadow Run, Hazel River and Hughes River, stream dewatering only occurred near the upper stream origination point. In contrast, Paine Run and Jeremys Run had several discrete dewatering sections further from their origination points (examples shown in Figs. 3d and A3). During the drier period 17–19 September 2019, no dewatering was found along lower Staunton River, though Piney River had seven discrete dry patches where none were mapped in 2016, and similar patterns were observed for those two streams in 2021 (Fig. 7). Paine Run had 29 discrete zones of dewatering in 2019, distributed mainly along the central and upper sections of the stream corridor, and showed extensive dewatering in 2021 such that more than 75 % of the stream length was either dry or of isolated pools (Figs. 6–8).

Figure 6 Results from 2016, 2019 and 2021 longitudinal channel dewatering surveys conducted by physical observation where the 2019 and 2021 data are shown offset laterally from the stream channel where those surveys occurred.

Figure 7 Zoom views for the three study subwatersheds where stream dewatering observations were also collected over three summer seasons (2016, 2019, 2021).

Figure 8 The results of the 2019 stream drying survey and 2020 high spatial resolution HVSR measurements are shown over the lidar hillshade in plan view (a, b), and along a lidar-derived stream elevation profile cross-section view (c, d) for Piney River (a, c) and Paine Ru (b, d).

4.3 Stream temperature patterns

Paired air and streamwater annual temperature signals extracted exhibited a spectrum of apparent shallow groundwater influence, manifested as ubiquitous positive phase shifts (air to water) per methods described by Briggs et al. (2018a) and Hare et al.  (2021), ranging approximately 5 to 30 d with a mean of 11 d. Reduced annual temperature signal amplitude ratio (water : air) generally corresponded with increased phase shift when all SNP stream monitoring sites are plotted in aggregate, indicating greater shallow groundwater inflow was generally related decreasing stream temperature sensitivity (Fig. 9a). Staunton River stream sites cluster together and show less signal phase shift (mean of 10 d) for similarly low amplitude ratio values (mean of 0.6) observed in other subwatersheds.

Figure 9 Panel (a) shows the annual temperature signal metrics for the study subwatersheds highlighted within the larger SNP dataset with conceptual groundwater end member signature trajectories. Panel (b) displays the downstream mean summer temperature profiles and 7 d maximum and minimum temperature ranges for Staunton River and White Oak Canyon.

Although originating as similar summer temperatures at respective upstream measurement points, the longitudinal mean, 7 d maximum, and 7 d minimum stream temperature profiles differed between Staunton River and White Oak Canyon where the latter had greater temperature variation and warming with downstream distance (Fig. 9b). The mean summer stream temperature had an approximate 2 ^∘C total range over the period of record for all sites as evaluated at the most downstream measurement point along each stream (Table 1). The warmest average (19.6 ^∘C) and 7 d maximum (23.6 ^∘C) was observed for the lower Jeremys Run site which was also at the lowest elevation. However, only 23 m higher in elevation, the downstream Staunton River site had the coldest average (17.4 ^∘C) and 7 d maximum (19.9 ^∘C) summer temperature. Piney River, which has the second largest median bedrock depth (2.54 m), had the second lowest average temperature (17.4 ^∘C). No significant relation was observed between elevation and mean summer temperature at the lower stream monitoring site, but a significant negative linear relation (R²=0.52; p<0.05) was determined between median stream corridor bedrock depth and mean summer stream temperature (Fig. 5b), with strong leverage on the linear fit imparted by the Staunton River datapoint.

4.4 Baseflow separation (index)

The summer season BFI determined for Paine Run, Piney River and Staunton River over the period of flow record shows substantial variability, but the median summer BFI for Staunton River (0.62) is approximately 50 % greater than Paine Run and Piney River (0.46 and 0.41, respectively, Table 2). For the primary study years of 2016–2019, Staunton River BFI is always largest, and all sites are above their respective interquartile range in 2017 but below their interquartile range in 2018 (Fig. 10). The anomalously low 2018 BFI values can be explained by extremely high summer precipitation that year (Fig. A2), resulting in total streamflow being dominated by runoff and quickflow as determined by baseflow separation. Mean summer streamflow over the period of record was highest for Piney River and lowest for Pain Run, and overall summer streamflow was most stable for Staunton River (lowest coefficient of variation).

Table 2 The median summer baseflow index (BFI), mean summer streamflow and mean summer standard deviation (SD) streamflow for three gauged streams from 1993–2020.

Figure 10 Summer baseflow index metrics summarized from 1993–2020 for three streams with specific values from the primary study years identified.

5 Discussion

5.1 Longitudinal spatial structure in observed bedrock depth

Seminal groundwater/surface water exchange research has indicated that bedrock topography along headwater streams may be a first-order control on the arrangement of nested gaining and losing flowpaths (e.g., Tonina and Buffington, 2009), and increased depth to low permeability bedrock is recognized as a (likely) primary driver of stream disconnection during dry periods that could be exacerbated by climate change (Ward et al., 2020). However, despite their importance to a range of headwater stream physical processes and cold-water habitat, local bedrock depth data are almost universally lacking, even in heavily studied experimental watersheds. Our study provides new inferences regarding the effects of bedrock depth on groundwater exchange and consequent effects on stream dewatering and temperature patterns at ecologically relevant spatial scales in mountain streams. The combined datasets indicate stream channel bedrock depth assessments may be necessary to support stream habitat assessments and predictions of stream connectivity under drought and climate change when existing large-scale geologic datasets are not of sufficient spatial resolution to support natural resource management applications.

Bedrock depth varied substantially within and among several of the eight study SNP subwatersheds but was predominantly shallow. For half of the subwatersheds (Hughes River, Meadow Run, Jeremys Run and Hazel River), median bedrock depth along the stream channel and lateral riparian zone was less than 2 m and did not show notable variability with distance, outside of one 12.8 m depth-to-rock location at upper Hazel River (Fig. 4). This anomalous measurement at Hazel was collected lateral to the stream on a valley terrace of colluvium in the vicinity of the only cold (approximately 10 ^∘C at land surface) riparian spring that was observed during all HVSR surveys. Bedrock depths of greater than 8 m were also found along the upper White Oak Canyon riparian zone as well (Briggs et al., 2018a), also associated with surficial seepage. Two anomalous bedrock depth measurements of 22.9 and 26.8 m were collected along the Piney River channel, but instead of being associated with groundwater springs, they coincided with a discrete sections of channel dewatering at baseflow during 2019 and 2021. Therefore, it appears that discrete zones of thick surficial material are the exception along SNP streams, though they can be important to localized processes such as focused riparian discharge and streamflow disconnection (latter is discussed in Sect. 5.2).

There are several existing sources of bedrock depth data that could potentially be used to inform headwater stream modeling and habitat assessment, but the accuracy of such datasets along headwater streams (typically away from existing boreholes) has generally not been evaluated at the scales of typical stream habitats. We conducted a point-scale comparison of our relatively high-resolution bedrock depth measurements to the global bedrock depth map of Shangguan et al. (2017), and found that bedrock depths were almost universally overpredicted at the SNP by large margins (Fig. A4). Specifically, predictions from the global-scale dataset exceeded HVSR measured depths by +12.2 m (mean) or approximately four times the average bedrock depth (2.9 m). This result may not be surprising as Shangguan et al. (2017) recognizes that large-scale bedrock depth interpolations are likely to overpredict shallow bedrock contacts especially in mountainous terrain with minimal available borehole data constraints. However, given that baseflow generation is expected to be dominated by shallow groundwater sourced from unconsolidated sediment in headwater systems with low permeability bedrock, our study highlights that use of large-scale bedrock depth layers may propagate substantial uncertainty into process-based groundwater flow model predictions when used to inform model structure in the absence of local measurements.

Publicly available maps of surficial geologic materials are another potential source of bedrock depth information. High-resolution digital soils maps are now widely available, including for the catchments of SNP, and these maps do capture some of the general depth-to-rock transitions between subwatersheds observed in this study. For example, publicly available soils maps from the US Department of Agriculture Natural Resources Conservation Service (https://websoilsurvey.sc.egov.usda.gov/App/WebSoilSurvey.aspx, last access: 12 October 2020) indicates that the White Oak Canyon stream corridor is comprised of silts, loams and stony soils with a general bedrock depth of approximately 1.2 m, which is in a similar range as most HVSR measurements made along the upper stream. However, the generalized soil units may not offer needed detail regarding site-specific valley sediment thickness for hydrogeological and ecological studies where information regarding within-watershed variation is critical. Along the lower Piney River where HVSR data had depths to rock ranging from 1.4 to 3.6 m, the NRCS soils map universally indicates silt and stony material > 2 m. Along Paine Run where the stream is often scoured to bedrock, the soils map shows consistent highly permeable sandy material with >2 m thickness. This discrepancy is understandable given most of the test pits were likely substantially further downstream in better terrain for agriculture. In conclusion, analysis of large-scale patterns from existing soils maps and interpolated/predicted bedrock depth layers indicates that geophysical mapping of bedrock depth may be needed to inform stream research and management, particularly in shallow, low-permeability bedrock terrain where there are few existing wells and boreholes of known geologic detail.

5.2 Summer stream dewatering related to bedrock depth

Aligned with the conceptual model of Ward et al. (2018), our central hypothesis was bedrock depth along the stream corridor acts as a primary control on longitudinal stream dewatering and flow disconnection during summer low flows (visual example shown in Fig. A3). We postulated that permeable streambed thickness may undulate along mountain stream channels, and relatively thick sub-stream sediment zones could accommodate the entirety of low streamflow volumes, locally disconnecting channels during seasonal flow recession. We found mixed support for this simple hypothesis. Hazel River and Hughes River were two of the three subwatersheds that had dry channel zones just downstream of their respective stream origination points in 2016, and these two riparian corridors also had their deepest riparian bedrock depths in those high-elevation areas. However, as discussed above, White Oak Canyon had relatively thick, porous sediment zone near the subwatershed outlet but did not show any zones of dewatering, nor did lower Staunton River in 2016, 2019 or 2021 despite having the deepest median bedrock contact. Jeremys Run had three mapped dry zones in 2016 (not surveyed in 2019), yet depth-to-rock in those areas was only approximately 2 m, though the HVSR data collection points were not perfectly aligned with the dry patches. To address this spatial mismatch in stream dewatering and HVSR data, we used the stream dewatering maps to guide two new high-resolution HVSR surveys in March 2020 along sections of Paine Run and Piney River with dynamic patterns of channel drying, as described below.

When bedrock depth data were collected at high spatial resolution, even more variability in bedrock topography/sediment thickness was revealed than in the original larger-scale stream corridor surveys, and that finer scale of information was relevant to understanding stream dewatering patterns. For example, during summer 2019, a 291 m length section of the lower Piney River was observed to be dry and immediately preceded by 62 m of isolated stream channel pools; a nearly identical dewatering pattern was observed there in 2021 (Figs. 7 and 8a, c). The upper portion of this major feature of stream disconnection corresponded directly with a transition in bedrock depth along the channel from approximately 3 m to adjacent measurements of 27 and 23 m. This “trough” in the bedrock surface can likely act as a streamwater sink (shown conceptually in Fig. 1b), routing surface water downward to the point of draining the channel locally in the summers of 2019 and 2021 but not in 2016 when precipitation (groundwater supply) was higher than normal. Further downstream, the bedrock depth returned to approximately 3 m near the furthest downstream measurement point, and flowing channel water was again noted during the drying surveys. Such a section of stream dewatering in the lower watershed would serve to impede fish passage along Piney River during the lowest flows, likely corresponding to times of maximum thermal stress when fish mobility is critical to seeking thermal refuge (Magoulick and Kobza, 2003).

Not all variability in bedrock depth below streams associated with stream drying was as dramatic as the Piney River example but can be important in disconnecting channel habitat in summer when baseflow supply is reduced. Paine Run is a more strongly confined stream valley that had 29 discrete zones of stream channel dewatering during September 2019 and extensive dewatering in 2021 (Figs. 6, 7, and 8b and d) when numerous dead brook trout were also noted. Paine also had the greatest total exposed bedrock out of any of the SNP subwatersheds in this study, indicating a highly constrained valley underflow reservoir. High-resolution bedrock depth data were collected over a Paine Run subreach with seven discrete dry patches ranging from 17 to 185 m in channel length, with many bordered by zones of isolated pools (Fig. 8b). A comparison of these patterns with bedrock depth along the channel shows the flowing sections of stream were dominated by exposed bedrock surfaces or thin sediment. However, a notable exception is toward the upstream end of this focus reach where depth-to-rock was consistently >2 m over the run-up to a large zone of disconnected channel with some isolated pools (Fig. 8b and d). This result suggests the losses of streamwater accumulated over this approximately 80 m channel distance. In the following downstream contiguous sections of dry channel and/or isolated pools, bedrock depth averaged a larger 3.3 m, indicating the entirety of streamflow was accommodated by the subsurface, congruent with our original hypothesis. However, knowledge of bedrock depth in isolation is clearly not sufficient to predict stream channel gaining, losing and disconnection patterns as the stream with the largest average bedrock depth, lower Staunton River (median depth-to-rock 3.4 m, Fig. 4), was not observed to dewater during any of the three physical surveys (Figs. 6 and 7).

5.3 Summer stream temperature and groundwater exchange dynamics

Although headwater stream heat budgets are complex, our data indicate groundwater connectivity plays an important role when stream temperatures are already close to aquatic species thermal tolerances. The apparent dominance of shallow (<3 m depth) groundwater discharge along White Oak Canyon contributed to the Briggs et al. (2018b) prediction that the lower reaches would not provide suitable brook trout habitat by the end of the century given anticipated atmospheric warming. Jeremys Run, a long (13.4 km) stream consistently underlain by a shallow bedrock contact (median depth < 2 m), already shows a 7 d maximum summer temperature that exceeds expected brook trout tolerances (i.e., >23.3 ^∘C mean weekly average temperature, Wehrly et al., 2007) along the lowest reach.

The underflow reservoir of headwater stream valleys integrates upgradient and lateral hillslope groundwater flowpaths which accumulate with distance when bounded by low permeability bedrock. The two subwatersheds with largest median bedrock depth along their respective upstream corridors had the coldest mean summer temperatures, with Staunton River standing out as distinctly colder and having the only 7 d max temperature below 20 ^∘C (Table 1). There was a significant relation between median bedrock depth and mean summer stream temperature at the lower stream sites but not with elevation (Fig. 5b), indicating exchange with groundwater exchange had disrupted the expected elevation control on lower reach cold water habitat. Surficial hillslope contributing area is often assumed a primary control on potential groundwater discharge at the stream subreach-scale. However, Staunton River also had the second smallest drainage surface area of all study subwatersheds. Further, Staunton River had an average valley bottom width that was less than other streams that were observed to dewater and, therefore, the strong baseflow supply there cannot be explained by the lateral dimension of the underflow reservoir. This apparent conundrum indicates the importance of bedrock depth (supra-bedrock aquifer thickness) in facilitating spatially persistent baseflow generation during dry times, and we also found that the narrower headwater stream valleys of this study tended to have deeper bedrock depth (Fig. 5a).

For a more in-depth analysis of the paired bedrock depth and groundwater inflow controls on headwater summer stream dynamics, Staunton River can be contrasted with Paine Run. The latter had a similar total drainage surface area to Staunton River, with a >5 m (average) wider stream valley bottom but a 1.2 m shallower bedrock depth on average. Paine Run had dozens of dewatered stream channel sections in 2019 and 2021, and had a downstream boundary summer stream temperature that was 1.4 ^∘C warmer than Staunton River. In addition to a reduced average bedrock depth, Paine Run had numerous sections of exposed bedrock adjacent to localized pockets of stream channel alluvium and colluvium (Figs. 4 and 7), while extensive colluvial deposits along the Staunton River channel limited exposed bedrock to a few m-scale sections associated with pool steps (Fig. 4). Lower Staunton experienced major debris flows in June 1995 (Morgan and Wieczorek, 1996), events that likely created an enhanced local groundwater reservoir within coarse hillslope material compared to other SNP subwatersheds.

Based on the integrated datasets from these two SNP streams, we conclude that groundwater exchange is a critical factor determining whether headwater streams will warm and dewater in summer, which in turn is controlled in part by the thickness of supra-bedrock unconsolidated aquifer. As noted above, annual temperature metrics indicated a consistently deeper groundwater discharge influence along Staunton River while Paine Run had annual signal metrics that mainly indicated reduced and/or more shallow groundwater influence (Fig. 9a). Long-term streamflow and baseflow analysis from these streams showed Staunton River had higher but more stable summer discharge (Table 1), and substantially higher median summer BFI (0.62 vs. 0.46), indicating greater dominance of groundwater as a generator of streamflow compared to runoff and quickflow. Previous research in SNP used paired air/streamwater temperature records, precipitation and landscape characteristics to statistically model “groundwater influence” by year on a scale of 0–1 at the 100 m scale along the streams of this study, where details are described by Johnson et al. (2017). Although this previous work only extended to 2015, that year had analogous BFI scores to 2019 for Staunton River (0.88 vs. 0.84) and Paine Run (0.60 vs. 0.58). Comparing the 2019 drying survey observations to the 2015 high-spatial-resolution modeling of groundwater influence, we found that Paine Run was predicted to have groundwater influenced tributaries; but along the mainstem, where extensive dewatering was observed, there was substantially reduced modeled groundwater influence compared to the mainstem of Staunton River (Fig. 11). Johnson et al. (2017) also found a negative relation between valley bottom width and their metrics of groundwater influence on SNP streams. In the context of our finding that bedrock depth is negatively related to valley bottom width, we find further support for the hypothesis that thicker headwater stream valley sediments are influential to baseflow generation in low permeability bedrock settings.

Figure 11 The 2019 stream dewatering survey data (lines; this study) plotted offset of the mainstem, and 100 m groundwater influence predictions (points from Johnson et al., 2017) plotted along the mainstem and tributaries of Staunton River and Paine Run.

This observation and model comparison represents another line of evidence that groundwater connectivity at the subreach-scale is key in determining whether local increases in depth-to-bedrock drive channel dewatering at low flow. The impact of reduced underflow groundwater supply on stream disconnection is likely exasperated by the extensive zones of exposed bedrock along Paine Run (Figs. 4 and 7d), which locally reduce groundwater mounding in stream valley sediments as shown conceptually in Fig. 1b, such that abrupt increases in bedrock depth cause stream dewatering. Among the eight streams investigated here, Staunton River likely represents the most resilient summer cold water habitat, which could not be predicted using bedrock depth data alone but necessitated paired assessment of groundwater discharge dynamics.

6 Conclusions

In steep mountain valley stream systems underlain by low-permeability bedrock, the longitudinal underflow reservoir serves as a complex mechanism of streamflow generation, streamflow losses and stream temperature control (Figs. 1 and A1). Our study utilized complementary geophysical, temperature and hydrologic data at the scale of eight subwatersheds to highlight apparent tradeoffs in bedrock depth, shallow groundwater supply and the quality of cold-water habitat. Certain mountain stream corridor parameters may be reasonable to assume or infer from high-resolution topographic data, such as surficial sediment permeability (based on land surface roughness) and stream valley width, which can be important controls on whether underflow serves as a net source or sink of streamwater (Ward et al., 2018; Flinchum et al., 2018). However, as shown here, advances in predicting hydrologic connectivity and thermal variation along mountain stream networks may also require local evaluation of bedrock depth and stream–groundwater exchange.

When local increases in bedrock depth are not balanced by groundwater inflow, streams may be expected to dewater and disconnect under low flow conditions, and streams with reduced deep groundwater influence or shallower-sourced groundwater show warmer summer temperatures. Contrary to what might be expected, we found that mean summer stream temperature was not significantly related to elevation at the most downstream stream sites but instead was (negatively) related to average stream bedrock depth. Staunton River had the coldest summer stream temperatures and most pronounced deeper groundwater signatures. However, that subwatershed was of relatively small total surface area and average valley bottom width. The defining physical feature of Staunton River was that it had the largest average bedrock depth of all the eight SNP study streams at 3.4 m, allowing greater overall storage of recharge and baseflow generation. The other two gauged streams had substantially reduced baseflow indices, indicating streamflow generation was dominated by runoff and quickflow.

Overall, SNP streams tended to have consistently shallow bedrock depth, though a subset was more variable or had spatial trends and discrete features. Observed channel dewatering patterns during late summer baseflow periods were related to local-scale variation in bedrock depth, such as a discrete feature of greater than 20 m depth observed along Piney River that caused repeated streamflow disconnection. However, in other streams, more subtle bedrock depth variation also caused channel dewatering, indicating the importance of local hydrogeological context in determining the importance of bedrock depth on streamflow connectivity. For example, patchy 2–4 m deposits of sediment adjacent to exposed bedrock along Paine Run caused extensive summer dewatering in 2019 and 2021, and during the latter survey, several dead brook trout were noted in the disconnected sections. Paine and Piney also showed enhanced dewatering during the summers of 2019 and 2021 compared to the wetter 2016 summer, demonstrating the additional control of recent precipitation on stream disconnection in headwater systems that do not efficiently store water.

Lateral groundwater inflow through high-permeability, unconsolidated sediments is a critical component of headwater stream baseflow (Tran et al., 2020). Shallow, low permeability bedrock can constrain lateral flowpaths and underflow to the near-surface critical zone where it is highly sensitive to enhanced evapotranspiration, temperature increase and drought under climate change (Hare et al., 2021; Condon et al., 2020). As it becomes increasingly important to understand and predict the resilience of mountain cold-water stream habitat at a fine-spatial grain, continued coupled advances in geophysical characterization, stream temperature monitoring and groundwater exchange analysis are needed.

Appendix A

Figure A1 The headwater streams of Shenandoah National Park, Virginia, USA are expected to flow over coarse alluvium and colluvium, and have connectivity to shallow hillslope groundwater and underflow but reduced connectivity to deeper bedrock groundwater (modified Fig. 26 in Nelms and Moberg, 2010) US. Geol. Surv. Investigations Rep. 2010-5190.

Figure A2 Monthly precipitation and air temperature data derived from the Luray weather station (GHCND:USC00445096) located within Shenandoah National Park. Box plots show the distribution of values for the period of record (1942–2020) with the limits of the box containing 50 % of the values, whiskers containing 90 % of the values and solid line in boxes depicting the median value. The lines represent values for the four primary study years.

Figure A3 Images from the same vantage point along Paine Run during (a) high and (b) low stream flow times, the latter showing channel dewatering associated with a deposit of coarse alluvium across the channel, creating locally enhanced streambed water storage.

Figure A4 Comparison between bedrock depth modeled for the globe by Shangguan et al. (2017) at a 250 m resolution and the HVSR-calculated depths-to-bedrock in this study.

Appendix B

Table B1 Summer stream temperature metrics for each study subwatershed determined from the data set of Snyder et al. (2017), https://www.sciencebase.gov/catalog/item/get/594bdc88e4b062508e385039 (last access: 7 August 2021).