4.5 Relating SC-Q hysteresis, streamflow components and watershed properties
On an annual scale, we detected a significant negative correlation between quickflow proportion and the mean annual hysteresis index at the 50th percentile of discharge values (Table 7; Pearson’s r = -0.77, p-value < 0.05). Watersheds with higher proportions of quickflow tended to have larger negative hysteresis indices implying these watersheds behaved the flashiest in response to diurnal snowmelt pulses. The other runoff generation mechanisms (baseflow and throughflow) were not correlated to any of the annual hysteresis indices, supporting evidence that the youngest fraction of runoff controls stream water solute concentration dynamics (Benettin et al. , 2017).
We found strong relationships between direct runoff partitioning (i.e. throughflow vs quickflow) and slope aspect of our study watersheds. Watersheds with more land area facing north were associated with greater annual proportions of quickflow (Table 7; Pearson’s r = 0.72, p-value < 0.05). North-facing aspects tend to receive less direct radiation, and snowpacks persist longer compared to south-facing slopes. This allows snowpacks to melt later in the year at higher temperatures and thus faster rates (Musselman et al. , 2017), which brings about conditions favorable for quickflow. Longer persistence of snowpacks on north-facing slopes leads to the development of preferential meltwater flowpaths at the snow-soil interface, resulting in as much as 170% more snow water equivalent accumulating along the length of north-facing hillslopes compared to south-facing hillslopes (Webb et al. , 2018).
The percentage of surficial geology covered with glacial deposits was significantly negatively correlated to the proportion of annual streamflow derived from quickflow (Table 7; Pearson’s r = -0.66, p-value < 0.05). Watersheds with low proportions of quickflow (NASH100, NASH200, NFLL100, NFLL200) had a larger percentage of glacial deposits covering their surficial and bedrock geologies (Tables 2, 4 & S1). In contrast, watersheds with larger proportions of quickflow tended to have bedrock geology dominated by metasedimentary and metavolcanic rocks associated with the Libby Creek Formation (Table S1) and surficial geology dominated by glaciated bedrock (Table 2). Glacial deposits were assumed to have higher infiltration capability and thus more drainage capacity thereby reducing the occurrence of quickflow relative to the bedrock associated with the Libby Creek Formation. This bedrock can be considered to have lower infiltration capacity, increasing the occurrence of quickflow.
The percentage of surficial geology covered by glacial deposits was significantly correlated to annual and daily hysteresis indices (Table 7). Watersheds with a higher proportion of glacial deposits tended to produce large, positive hysteresis indices. Waters flowing through glacial deposits, compared to glaciated bedrock, likely have longer residence times and thus more potential geochemical evolution before arriving as streamflow. These factors would lead to more variability in flowpath development and larger hysteresis indices observed at annual and daily scales.

5. DISCUSSION

5.1 Dominant controls on hydrologic response

Using continuous records of SC and Q we showed that direct runoff (quickflow and throughflow) contributes the majority of the total annual streamflow - especially during the snowmelt period (Table 3). These results align with other investigations near our study area that used SC (Miller et al. , 2014; Rumsey et al. , 2015) or stable water isotopes (Huth, Leydecker, Sickman, & Bales, 2004; Liu et al. , 2004; Williams, Seibold, & Chowanski, 2009) for hydrograph separation. Hydrograph components were significantly correlated to aspect and surficial geology, while elevation was highly positively correlated to snow depth, mean precipitation, total runoff, and runoff ratio (Table 7). The strong link between elevation and precipitation has been well documented in other mountainous environments that span large elevation ranges (Elder, Dozier, & Michaelsen, 1991). At smaller scales, however, slope aspect impacts the ability of wind to redistribute snow (Dadic, Mott, Lehning, & Burlando, 2010) and vegetation structure significantly alters snow accumulation and ablation patterns (Varhola, Coops, Weiler, & Moore, 2010).
Streamflow and SC data from our watersheds indicated a combination of geology, topography, and precipitation inputs strongly affect the hydrologic response in snow-dominated, headwater catchments (Segura et al. , 2019). As temperatures warm, snowmelt will initiate earlier in the season at slower rates (Musselman et al. , 2017) which could bring about a slower hydrologic response and less near-surface saturation, favoring groundwater recharge and less direct runoff. We can already see evidence of this in our system where watersheds that have more south-facing area are associated with greater proportions of baseflow. Differences in slope aspect will be amplified in the future when snowmelt occurs at earlier dates and the sun angle is lower (Lundquist & Flint, 2006). In future warming scenarios, shaded hillslopes will provide habitat to native species that developed under cooler, wetter scenarios, whereas sunnier areas will be more susceptible to drought and invasion from non-native species.
For the watersheds studied here, longer snow persistence on north-facing or heavily shaded hillslopes enabled high snowmelt rates which bring soils to field capacity and facilitates rapid shallow subsurface flow and return flow (Dunne & Black, 1971). Water typically exfiltrated near the break in slope where hillslopes are adjacent to riparian areas, supporting the variable source area concept (Hewlett, 1961). Return flow exfiltrating from soil and regolith layers is an important contributor to groundwater recharge and has been shown to be resilient to drought and can buffer recharge under climate change (Carroll et al. , 2018). Conversely, south-facing slopes receive more direct radiation during snowmelt compared to north-facing slopes such that snowmelt is typically initiated earlier in the year. This likely contributes to the low percentages of quickflow produced in the watersheds with the greatest percentage of south-facing area (Table 4). This interpretation is consistent with the results from Thayer et al. (2018). Using time-lapse electrical resistivity tomography, they showed the preference for deep drainage occurring on a south-facing hillslope located within NONM100 and concluded that overland flow and lateral shallow-subsurface interflow (i.e. quickflow) were negligible on south-facing hillslopes.
Interestingly, two of our smaller study watersheds and the two lowest total streamflow producing sites (GOLD100 and NONM100; Table 4) had remarkably different runoff generation mechanisms from one another (Figures 2 & 5; Tables 3 & S5). Baseflow contributions during snowmelt were the smallest for the GOLD100 watershed while baseflow contributions were the largest during low flow conditions. This shift in baseflow dominance, from low proportions during high discharge to high proportions during low discharge, suggests that a smaller groundwater storage capacity exists for this watershed. The surficial geology of GOLD100, which is dominated by grus mixed with alluvium, residuum, slopewash and colluvium, is remarkably different than the other study watersheds and likely helps explain the hydrologic response (Table 2). Unlike the other study watersheds, the Libby Glacier only overlaid a small portion of this watershed (Atwood, 1937), resulting in few glaciated deposits which we expect have greater infiltration capacity and the ability to store more water which leads to more geochemical evolution compared to the surficial geology units present in GOLD100 which resulted in a flashy hydrograph behavior during snowmelt.
In contrast, the surficial geology in NONM100 is almost entirely glaciated deposits. The stream draining this watershed received a more constant source of baseflow throughout the year, implying groundwater is an important and stable source of stream water in this small catchment (Carroll, Deems, Niswonger, Schumer, & Williams, 2019). Reductions in snowpack, the primary water storage component in runoff-dominated watersheds like GOLD100, will be more disruptive to the natural flow regime in these types of catchments, resulting in a larger contraction of the stream network compared to streams that are more groundwater dominated (Tague, Grant, Farrell, Choate, & Jefferson, 2008). Reductions in late season flow and subsequent contraction of the stream network, for example under climatic shifts in seasonality, could potentially reduce aquatic habitat heterogeneity and stream macroinvertebrate biodiversity (Brown, Hannah, & Milner, 2007).

5.2 Specific conductance – discharge relationships reflect the storage and release of snowmelt water

While the SC-Q relationship remained consistently negative at annual and daily scales, the direction of hysteresis during snowmelt-induced diurnal cycles (clockwise) was opposite of the direction observed at the annual scale (anti-clockwise). The clockwise hysteretic behavior observed on a daily scale suggested throughflow, with its greater SC and relatively longer residence time, was first flushed out of hillslopes, followed by freshly melted snow, with lowest SC and shortest residence times, that filled in pore space occupied by the discharged throughflow. This behavior supports threshold behavior for runoff generation, where connectivity must be achieved before daily meltwater with the lowest SC contributes to streamflow (Tromp-Van Meerveld & McDonnell, 2006).
In watersheds with the largest mean daily hysteresis index (NFLL100 and NFLL200) the shape of the daily SC – Q hysteresis during snowmelt-induced diurnal cycles was relatively circular during the rising limb and became flatter during peak flow and falling limb conditions (Figure 7). This behavior was documented by Kobayashi (1986) and was attributed to an increase in the subsurface component of streamflow after a snow-free area emerged adjacent to the stream channel. Our results from this change in hysteretic shape support evidence that greater flow pathway variability was present during early snowmelt where water was moving through both relatively shallow and deep flow pathways. During the falling limb, less snow was present in the watershed and throughflow from relatively deeper flow pathways draining soil and regolith layers combined with baseflow were primarily responsible for streamflow generation. This dominance of deeper flow pathways leads to SC values measured in streamflow that vary more linearly related to discharge over daily timescales.

5.3 Potential limitations and future work

Based on field observations in our study watersheds, it is obvious that the vast majority snowmelt infiltrates the ground surface for a period of time. We make no distinction between surface (i.e. overland flow) and subsurface runoff in the nomenclature of our selected hydrograph components based on extensive field observations. During peak snowmelt conditions, water stored in hillslopes was rapidly transported via lateral flowpaths arising from macropores and / or interfaces between depositional units with large density contrasts (Fullhart et al., 2019; McNamara, Chandler, Seyfried, & Achet, 2005; Roberge & Plamondon, 1987; Uchida, Tromp-Van Meerveld, & McDonnell, 2005; ) indicative of threshold runoff response (Penna, Tromp-Van Meerveld, Gobbi, Borga, & Dalla Fontana, 2011; Spence, 2007; Tromp-Van Meerveld & McDonnell, 2006). This quick moving lateral runoff was observed in NONM100 with low SC, exfiltrating from the base of hillslopes where fully saturated riparian areas could not accommodate additional subsurface flow (Figure S4).
As mentioned before, interpretation of our SC-based baseflow hydrograph separation hinges on the ability to accurately assign concentrations to end-members. Due to the fact that SC concentrations are effectively constant throughout the low flow period, we believe the selection of the baseflow end-member is fairly accurate. However, relatively few snow and snowmelt samples were used to assign the direct runoff component. Extensive field mapping and SC measurements of the stream network expansion and contraction in NONM100 during snowmelt revealed that exfiltrated shallow subsurface flow (i.e. return flow) often had an SC concentration similar to or lower than the value assigned to the direct runoff component (21.6 µS/cm; Figure S4). We were surprised by the lack of geochemical evolution of this fast-moving return flow and hypothesize that water flowing along these flowpaths had the shortest residence times. The value we assigned to the direct runoff component is in between two values chosen assigned by Miller et al. (2014) (10 and 33 µS/cm). Assigning a lower value to our direct runoff end-member would have had little implications to our results since the difference between direct runoff and baseflow SC is so large in our study watersheds.
Persistent streamflow occurs in all watersheds throughout the winter, except perhaps GOLD100 and NONM100 (not shown). In most of the watersheds stream stage and SC data are constant, indicating stable sources of baseflow contribution during the winter. Unfortunately, like many other similar studies, we collected data only during the relatively warm period of the year when ice is less likely to affect data quality and freezing is less likely to damage equipment. Our results are, for example, underrepresenting the amounts and proportions of baseflow by only using data from the relatively warm period of the year where data collection is more feasible.
While SC is a robust measurement of concentration and automated loggers can record instantaneous values at relatively low costs, information is lost when solute concentrations are aggregated together. Many studies have shown that concentration-discharge relationships vary for different individual solutes (e.g. Lewis & Grant, 1979; Godsey, Kirchner, & Clow, 2009). Exploring individual solute concentration-discharge relationships may help refine our knowledge on the number of different flow pathways water may take before arriving as streamflow. Further, stable water isotope composition sampled from our study watersheds did not permit the use of hydrograph separation due to event and pre-event waters having essentially the same concentration (analysis not shown). This is consistent with previous work demonstrating that stable water isotopes have proven less useful in snowmelt-dominated, seasonally arid environments where recharge is dominated by snowmelt (Earman et al., 2006), and thus indistinguishable from event water (Jin, Siegel, Lautz, & Lu, 2012). Regardless, other process-relevant information may be gained by exploiting the isotopic signatures and variability in the translation of snowmelt to streamflow during different discharge conditions to elucidate controls on runoff generation.

6. CONCLUSIONS

Specific Conductance (SC) recorded in our study watersheds showed large amounts of dilution in peak streamflow periods indicating that seasonal snowmelt dominated annual streamflow contributions. Streamflow and SC also showed distinct hysteretic patterns at annual and daily scales. Annually, most hysteresis indices that could be identified were negative, suggesting faster flow pathways dominated streamflow on the rising limb of the annual hydrograph compared to slower flow pathways occurring during the falling limb. During snowmelt-induced diurnal cycles, SC-Q hysteresis was consistently negative, indicating the daily melt water with the lowest SC preferentially contributed to streamflow on the falling limb of diurnal streamflow cycles after throughflow water was displaced. Hysteresis indices derived from SC-Q relationships were significantly correlated to proportions of quickflow and surficial geology, supporting the idea that the youngest fraction of runoff controlled stream water solute concentration dynamics in these systems.