Upland mountain watersheds host active hydrologic systems that play important roles in the storage and release of water to the larger stream network (Hale et al., 2016; Hale and McDonnell, 2016; Liu et al., 2004; Williams et al., 2015). These groundwater systems are locally important sources of streamflow (Asano et al., 2022; Gabrielli et al., 2012; Uchida et al., 2003) and recharge regional groundwater systems, which discharge to larger streams in the lowlands (Beisner et al., 2018; Fan et al., 2013; Frisbee et al., 2011; Gardner et al., 2011; Gleeson and Manning, 2008; Smerdon et al., 2012). The amount of bedrock groundwater discharge in montane streamflow is important for conceptualizing watershed function and predicting hydrologic resource availability in the future (Condon and Maxwell, 2019; Meixner et al., 2016; Williams et al., 2015). However, groundwater flow systems are marked by extreme heterogeneity in hydraulic characteristics and the resulting flow path length and rate; as a result, groundwater discharge volumes and locations are all highly variable (Beisner et al., 2018; Gardner et al., 2011; Hale et al., 2016; Harrington et al., 2013; Smerdon and Gardner, 2021). Neither the spatial scale at which 3D groundwater flow becomes an important and predictable source of streamflow nor the seasonal variation in these processes are well-documented at present.

Subsurface water is an important source of streamflow, even in snow-dominated watersheds (Liu et al., 2004). Subsurface sources are commonly conceptualized as originating from both bedrock (weathered bedrock and saprolite) and soil reservoirs (Anderson et al., 1997a; Carroll and Williams, 2018; Hale et al., 2016; Hale and McDonnell, 2016; Salve et al., 2012; Torres et al., 1998). In this study, we refer to shallow groundwater in the soil zone as soil groundwater and deeper, fractured-bedrock groundwater as bedrock groundwater. The relative importance of soil and bedrock reservoirs depends on the soil and bedrock hydraulic characteristics (Flint et al., 2008; Gardner et al., 2020; Katsuyama et al., 2005; Pfister et al., 2017; Uchida et al., 2002, 2006). The amount of groundwater in a watershed significantly influences dynamic properties, such as hydrograph variability (Asano et al., 2022; Gardner et al., 2010; Pfister et al., 2017) and mean transit time (Asano et al., 2002; Asano and Uchida, 2012; Hale et al., 2016; Pfister et al., 2017; Uchida et al., 2006).

Soil water lateral flow and bedrock groundwater circulation are tightly coupled (Appels et al., 2015; Carroll et al., 2019; Gardner et al., 2020; Tromp-van Meerveld et al., 2007; Uchida et al., 2002). As a result of this coupling, the location and amount of bedrock recharge are controlled by subsurface permeability distribution, hillslope position, and bedrock topography, which affect the lateral soil moisture distribution and location of vertical drainage (Appels et al., 2015; Gardner et al., 2020). The relative contributions of soil and bedrock groundwater to stream water are expected to vary temporally (Spencer et al., 2021). The celerity of subsurface response depends on the soil and bedrock hydraulic properties (Anderson et al., 1997a,b).

Groundwater discharge location, rate, and flow path source can be highly localized and variable because they depend upon the local geology and surface characteristics (Beisner et al., 2018; Bergstrom et al., 2016; Gardner et al., 2011; Harrington et al., 2013; Smerdon et al., 2012). In addition, the localized bedrock and soil water contributions from hillslopes are buffered by the local alluvial groundwater system, and they can be mixed and transported within this system before reaching the stream (Jencso et al., 2010). Despite the spatial heterogeneity at local scales, it is clear that the contribution from deeper bedrock groundwater systematically increases with increasing drainage area at regional scales (Asano et al., 2022; Frisbee et al., 2011, 2012). Since there is observed heterogeneity in the amount and location of groundwater discharge occurring locally, which contrasts with predictable increases in groundwater discharge at larger, integrated scales, it appears that groundwater discharge in streams is an emergent property.

Is there, in fact, a scale at which groundwater discharge volume becomes more predictable? If so, does this scale change temporally? In this article, we investigate how the amount and location of bedrock and soil groundwater discharge vary as a function of space and time in a montane, first-order watershed in west-central Montana. The study objectives are to constrain the amount of bedrock vs. soil groundwater discharge along a stream reach as a function of seasonal changes in catchment wetness and increasing spatial scale. Leveraging a high-resolution network of gauging stations, we conducted a multi-tracer interrogation of the incremental sources of streamflow generation. We estimate the total volume of subsurface discharge and distinguish the contributions of soil and bedrock groundwater discharge along the stream. We then explore relationships of the bedrock groundwater contribution as a function of the incremental and integrated upland characteristics.

The Cap Wallace Watershed (CWW) is located in west-central Montana, ~56 km northeast of Missoula on the north slope of the Garnet Range (Figure 1). The catchment drains 6 km² of forested land and is a tributary of the Blackfoot River. The CWW ranges in elevation from 1,163 m at the confluence of Cap Wallace Creek (CWC) and Elk Creek to 1,918 m at the ridge tops. Vegetation consists primarily of second-growth western larch (Larix occidentalis) on the north-slope mid and upper elevations, Douglas-fir (Pseudotsuga menziesii) on north-facing, low-elevation slopes, and Ponderosa Pine (Pinus ponderosa) on south-aspect slopes. North-aspect slopes have well-developed soils, whereas south-aspect slopes lack soil coverage and are generally covered by steep talus fields. Due to topographic variation in the radiation balance and resulting microclimatic effects (Holden et al., 2011a,b), south-aspect slopes remain significantly drier throughout the year (Hoylman et al., 2019).

At the Lubrecht Flume SNOTEL station (site #604; 1,425 m), the average annual precipitation over the record period of 10 years is 514 mm. Nearly 50% of precipitation falls as snow, classifying the Cap Wallace as a snowmelt-dominated watershed. In 2017, temperatures ranged from −27.6°C to 34.2°C with an annual average temperature of 6.1°C (Figure 2). Peak snow water equivalent typically occurs in late March, with snow melt taking place during April, May, and early June. May and June are typically cool and rainy, while July, August, and September are generally hot and dry, with precipitation beginning again in October.

The Garnet Range is part of the Northern Rocky Mountains and was formed through folding and thrusting associated with the Jurassic Sevier and Cretaceous Laramide Orogenies (140–55 Ma). Subsequent extensional normal faulting was initiated in the early Eocene as a result of Basin and Range tectonics, triggering the uplift of metamorphic core complexes in the Northern Rockies (Portner et al., 2011). The CWW is situated on the flank of a large granitic batholith called the Garnet Stock, which intruded the area during the Late Cretaceous (Lonn et al., 2010). CWC follows an east–west trending normal fault. Quaternary glacial and alluvial deposits cover the narrow valley bottom and stream floor, while the hillslopes are composed of gravelly, silt-loam soils (Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture, n.d.) that overlie several kilometers of argillites, siltites, and quartzites of the Precambrian metasedimentary Belt Supergroup. Soils are classified as Typic Haplustalfs (well-draining with a lower clay horizon), with depths ranging from 0.5 to 2 m that are thicker in hollows and thinner in upper elevation ridges (Hoylman et al., 2018). The north-aspect slope is composed of the Garnet Range Formation, a siltstone with hummocky sand lenses, which lies unconformably across the valley from the stratigraphically lower Bonner Quartzite, a massive, slightly metamorphosed, arkosic sandstone that makes up the south-aspect slopes (Lonn et al., 2010). These formations have been subject to tectonism associated with compression from major mountain-building events and subsequent active extension due to gravitational relaxation, leading to a dense fracture network within the bedrock (Brenner, 1968; Gardner et al., 2020).

We measured multiple environmental tracers in soil, groundwater, and stream water to perform hydrograph separation along CWC. We conducted synoptic sampling along CWC for discharge of the stream and dissolved ²²²Rn and major ions at four hydrologically important times: (1) early-spring baseflow (March), (2) the rising limb of spring snowmelt freshet (April), (3) the falling limb of freshet (May–early June), and (4) summer baseflow (late June–October).

Radon concentrations ranged from 0.074 to 4.25 Bq/L and showed the same general spatial pattern across varying discharges (Figure 3). During each sampling period, a significant increase in radon activity was observed approximately halfway down the observed reach (~2.5 km). This increase coincides with the location of a large contributing hillslope area from the south (Figure 1). The difference between duplicate samples ranged from 0 to 40%, with low-concentration samples exhibiting larger variability, consistent with Poisson statistics of radioactive decay. For example, at a concentration of 0.074 Bq/L (corresponding duplicate 0.111 Bq/L), there is a 40% variability, while at a concentration of 2.527 Bq/L (duplicate 2.597 Bq/L), there is a 2.7% variability.

Figure 4 shows total subsurface discharge estimates along with modeled and observed total stream discharge. The longitudinal profile of the accumulated subsurface discharge closely follows the increase in total stream discharge. For each interval, subsurface discharge accounts for 19–100% of the incremental gains in streamflow, with an average of 81% (Figure 5). In May and mid-June, the months of highest rainfall, the influence of non-subsurface discharge is strongest, with subsurface discharge accounting for, on average, 68% and 59%, respectively, of the incremental gains in streamflow. Non-subsurface discharge generation estimates predominantly occur near the center of the reach (2.5–3 km downstream of the headwaters) at the convergence of two large hillslopes.

Reprojected mean groundwater and soil water end-member concentrations are provided in Table 3. The two end members chosen, spanned the sample variance and the residual error in modeled vs. observed concentration showed no structure. The fractions of soil and bedrock groundwater in stream water samples are shown in Figure 6. The bedrock groundwater discharge fraction ranged from 0 to 0.45, with an average of ~0.2. In general, bedrock groundwater becomes a larger component of streamflow as the distance downstream increases, with the highest fraction of groundwater discharge occurring at the farthest downstream gauging location. Figure 7 shows the percentage of bedrock groundwater at the outlet as a function of time. Groundwater at the outlet fluctuated between 26 and 44% and averaged 38% of total streamflow. The fraction of bedrock groundwater to total streamflow at the outlet was highest during the late spring (2 June 2017) and lowest during mid-spring (19 May 2017).

Table 4 summarizes linear regression results for incremental and accumulated terrain characteristics with the bedrock groundwater mixing fraction. Linear regression of cumulative mean terrain indices and soil and bedrock groundwater fractions revealed strong relationships with both cumulative average elevation and cumulative upslope accumulated area (Figure 8, Table 4). R² values ranged from 0.52 to 0.98 over the 7-month study period. While the R² values ranged from 0.05 to 0.81 for the incremental terrain characteristics, they tended to switch drastically from high to low values throughout the year. No systematic trend was found relating incremental end-member discharges with incremental TWI and area (Figure 9, Table 4).

Figure 10 plots the r² values from cumulative and incremental regression with binary soil well responses as a function of time. Soil wells with high TWI values tended to stay continually saturated, whereas those with lower TWI values exhibited more transient-like behavior in response to climatic forcing. In both cases, the r² value was highest (0.98 and 0.81, respectively) in early July, ~3 weeks after the largest catchment “wet-up.”

Radon results indicate the majority of runoff in CWC is derived from subsurface sources (soil and bedrock groundwater) throughout the year, even during the snowmelt cycle, consistent with prior observations in alpine catchments (Liu et al., 2004). These results may be expected in later parts of the season (late summer and fall) when the overland flow is expected to be zero in our catchment. However, the fact that subsurface flow dominates during snowmelt (March, April, and May) indicates that direct runoff through the snowpack is not a major source of water for our stream. Radon results are consistent with soil hydrometric data. Soil wells with low TWI values (generally located at higher elevations with less accumulated area) exhibit transient behavior to catchment wet-ups (Figure 10) (Gardner et al., 2020). Some wet-ups had durations on the order of hours to days, suggesting that soil water in these landscape positions is rapidly lost through lateral inter-flow and/or infiltration to the bedrock aquifer.

While incremental discharge source plots show that subsurface sources account for the majority of streamflow gains, there are some instances where non-subsurface sources are estimated to be more important in streamflow generation (Figure 5). Non-subsurface streamflow generation occurred throughout the year but had the strongest influence on stream chemistry during the late spring to mid-summer (May through late July). These instances occurred primarily in the middle section of the reach, where multiple convergent zones from large hillslopes converge (Figure 3). Hydrometric data from soil wells indicate that the catchment was at its most saturated state during the May and June periods (Figure 10). Thus, during May and June, we hypothesize that non-subsurface sources could be due to saturation excess and overland flow in these large contributing drainage bottoms or soil and/or preferential/piping flow in the unsaturated zone. The fact that our radon modeling indicates some non-subsurface flow during the baseflow periods (late July and October) is likely due to model error. A probable source of error is the discharge of short residence time flow paths that have low Rn.

End-member mixing analysis provided an independent test of the discharge source findings, along with the ability to further divide stream-flow composition. The results indicate that the stream is generally composed of shallow saturated soil water at the headwaters and that shallow soil water is the dominant source of streamflow generation, but bedrock groundwater discharge increases with downstream distance, even during periods of snowmelt and heavy rainfall. Our bedrock groundwater discharge fractions are generally comparable with those of other mountain aquifer studies (Somers and McKenzie, 2020). The increasing influence of bedrock groundwater with downstream distance argues for a 3D watershed conceptual model consistent with those proposed by previous researchers evaluating groundwater flow in mountainous catchments in the western U.S. (Carroll et al., 2018; Frisbee et al., 2011, 2012). However, this study demonstrates that this emergent relationship exists in headwater catchments >2 km² in scale.

Temporal analysis indicates surprisingly low variability in modeled end-member proportions throughout the observed time period; however, early summer months have the highest proportion of bedrock groundwater discharge. Several large precipitation events (~60 mm/day) occurred in early June. These precipitation events could have rapidly filled the subsurface reservoirs, which could, in turn, drive increased discharge through a variety of mechanisms, including groundwater ridging, preferential recharge, and high fractured bedrock celerity. Rapid groundwater reactions to precipitation inputs are consistent with the experimental observations of Anderson et al. (1997a). However, our observed late-season bedrock groundwater ratios of streamflow are surprisingly low. Several previous hydrograph separation studies have assumed that baseflow is entirely composed of groundwater (Genereux et al., 1993; Pinder and Jones, 1969). Indeed, our results indicate that shallow soil groundwater is a substantial portion of streamflow during the late season. These results could be explained by several different mechanisms. For example, the substantial difference in the hydraulic diffusivity of soil and fractured bedrock can lead to situations where the fractured bedrock discharge exhibits a rapid response to infiltration, which is then followed by a long-term storage release from porous soil (Anderson et al., 1997b). Additionally, soil interflow and groundwater recharge and discharge have been shown to be tightly coupled at the scales of hillslopes (Gardner et al., 2020) and catchments (Carroll et al., 2019), and this interaction is capable of sustaining groundwater recharge and soil water throughflow throughout the year (Gardner et al., 2020).

Soil groundwater remains the dominant source of subsurface discharge over our entire reach throughout the year. However, bedrock groundwater discharge becomes an increasingly higher fraction of discharge as the amount of accumulated area increases. These observed trends of high soil groundwater discharge and increasing bedrock groundwater contribution with increased catchment scale imply that streamflow represents an integration of individual hillslope responses and larger-scale deeper groundwater flow paths (Figure 11). We attribute the increasing contribution of bedrock groundwater discharge with increasing accumulated area due to the integration of more heterogeneity and higher groundwater discharge zones (Asano et al., 2022; Frisbee et al., 2011). However, we demonstrate this relationship at a landscape scale that is orders of magnitude smaller than what has been observed in previous research. The relationship was generally consistent across variable streamflow states, with small temporal inconsistencies that are likely associated with dynamic catchment properties, such as climate and soil moisture, highlighting the importance of antecedent conditions.

The lack of trend in bedrock groundwater discharge with incremental sub-watershed area suggests that a minimum scale exists below which it is difficult to observe trends in bedrock groundwater discharge with terrain metrics. This lack of trend could result from two proposed processes: (1) regional groundwater flow paths could contribute to streamflow in a non-linear fashion due to geologic heterogeneities and (2) valley-bottom exchange of water. Linear trends in modeled soil and bedrock groundwater contributions against elevation and upslope accumulated area should not necessarily be expected, as a sub-catchment high in the watershed could export its regional groundwater to locations lower in the watershed. Thus, in the case of active regional groundwater flow, low-elevation sites would receive disproportionally more groundwater than high-elevation sites, potentially explaining the observed correlations with accumulated catchment area.

In this study, multiple environmental tracers were used to estimate the discharge of soil and bedrock groundwater along a 5 km reach of a snowmelt-dominated mountainous catchment at a variety of hydrologically important times. ²²²Rn modeling indicated that streamflow is generated predominantly from subsurface sources for all sampled periods. EMMA results indicated that streamflow chemistry could be explained by discharge from soil and bedrock groundwater. On average, bedrock groundwater comprised 38% of streamflow at the watershed outlet but fluctuated between a low of 26% in spring and a high of 44% in early summer. Bedrock groundwater discharge for individual stream segments was not well-correlated with landscape characteristics for incremental sub-watersheds draining to the segment (01 km²); however, it was well-correlated with the aggregated watershed landscape characteristics upstream of that segment (>2 km²). EMMA-derived mixing models showed that headwaters are dominantly composed of soil groundwater and that bedrock groundwater becomes an increasingly more important component to streamflow at larger catchment scales with lower elevations and higher UAA values, suggesting 3D topography-driven flow. For the accumulated watershed characteristics, correlations among landscape and end-member discharge were strong across variable states of catchment “wetness” conditions, and accumulated elevation or total upslope area were robust predictors of groundwater discharge during all seasons (r² = 0.52–0.98). Groundwater discharge can be expected to constitute a significant source of streamflow generation in upland and montane catchments throughout the year. Bedrock groundwater flow paths follow a 3D nested flow pattern where recharge originating in upper reaches can discharge in further downstream reaches according to the heterogeneous subsurface architecture. Groundwater discharge locations are strongly controlled by local heterogeneities and cannot be well-predicted from the incremental contributing area at scales below 1 km². As the watershed scale exceeds 2 km², the aggregated average groundwater discharge from these 3D flow paths increases with an increase in drainage area. These findings could be useful to other researchers in developing conceptual models of streamflow generation in other watersheds.