Main Text File
Quantifying daily versus annual contributions of snowmelt water to streamflow using graphical and geochemical hydrograph separation

Abstract

In this study, we characterize the snowmelt hydrological response of nine nested headwater watersheds in southeast Wyoming by separating streamflow into three components using a combination of tracer and graphical approaches. First, continuous records of specific conductance (SC) from 2016 to 2018 were used to separate streamflow into direct runoff and baseflow components. Then, diurnal streamflow cycles occurring during the snowmelt season were used to graphically separate direct runoff into quickflow, representing water with the shortest residence time, and throughflow, representing water with longer residence time in the soil column and/or regolith layers before becoming streamflow. On average, annual streamflow was comprised of between 22% to 46% baseflow, 7% to 14% quickflow, and 46% to 55% throughflow across the watersheds. We then quantified hysteresis at both annual and daily timescales by plotting SC versus discharge. Annually, most watersheds showed negative, concave, anti-clockwise hysteretic direction suggesting faster flow pathways dominate streamflow on the rising limb of the annual hydrograph relative to the falling limb. At the daily timescale during snowmelt-induced diurnal streamflow cycles, hysteresis was negative, but with a clockwise direction implying that quickflow peaks generated from the concurrent daily snowmelt, with shorter residence times and lower specific conductance, arrive after throughflow peaks and preferentially contribute on the falling limb of diurnal cycles. Slope aspect and surficial geology were highly correlated with the partitioning of streamflow components. South-facing watersheds were more susceptible to early season snowmelt at slower rates, resulting in less direct runoff and more baseflow contribution. Conversely, north-facing watersheds had longer snow persistence and larger proportions of direct runoff and quickflow. Watersheds with surficial and bedrock geologies dominated by glacial deposits had a lower proportion of quickflow compared to watersheds with large percentages of metasedimentary rocks and glaciated bedrock.
Keywords: snowmelt, hydrology, streamflow, watershed, tracer, specific conductance, hydrograph separation, hysteresis

INTRODUCTION

Rivers draining mountainous regions around the world are vital for supplying water for both consumptive uses and ecosystem services. This is especially true in regions like the Southern Rocky Mountains where the surrounding basins are classified as semi-arid to arid and rely predominantly on mountain runoff to meet water demands. Agriculture in this area, for example, is considered at risk from changing snowmelt patterns associated with global climate change (Qin et al., 2020). In such systems, snow is the dominant form of precipitation, but warming temperatures threaten to decrease annual snowpacks and shift seasonality (IPCC, 2013). Since snowmelt is more efficient than rain at generating streamflow (Li, Wrzesien, Durand, Adam, & Lettenmaier, 2017) and recharging aquifers in many mountainous regions (Earman, Campbell, Phillips, & Newman, 2006; Jasechko et al., 2014), reductions in snowpack pose a potential disruption to the delivery of water to streams and aquifers. Increases in temperature can reduce snowmelt rates (Musselman, Clark, Liu, Ikeda, & Rasmussen, 2017) and snow water equivalent (SWE) (Hamlet, Mote, Clark, & Lettenmaier, 2005), resulting in earlier timing in the onset of streamflow in western North America (Stewart, Cayan, & Dettinger, 2005) and subsequent decreases in summer streamflow (Rood et al., 2008) – potentially leading to habitat loss for endemic species (Lytle & Poff, 2004).
Changes in snowmelt are of course only part of the story for how mountain rivers and streams respond hydrologically under climatic shifts. The time it takes for water to move through a watershed and arrive in streams, thus controlling the shape of the hydrograph, is the result of interactions between precipitation inputs and the eco-physical characteristics of the watershed (Singh, 1997). Variations in vegetation composition (Molotch et al. , 2009; Vivoni et al. , 2008), topography (Broxton, Troch, & Lyon, 2009; Webb, Fassnacht, & Gooseff, 2018), and the lithology and structure of underlying soil and geologic formations (Mueller, Weingartner, & Alewell, 2013; Nippgen, McGlynn, Marshall, & Emanuel, 2011) can influence the residence time and flow pathways of water in concert thereby impacting both the quality and quantity of water arriving in streams.
Detangling such interactions is particularly challenging in mountainous regions given the complexity and remoteness of the landscapes (e.g. Jantze, Laudon, Dahlke, & Lyon, 2015). Still, efforts to characterize older groundwater contributions to streamflow (i.e. baseflow) in snowmelt-dominated headwater mountains regions serve to increase our understanding on how groundwater interacts with the environment and has important implications for the vulnerability of water resources and ecosystems under warming conditions (Miller, Buto, Susong, & Rumsey, 2016). The effect of the underlying geologic formations becomes most apparent during dry weather and baseflow conditions, when direct runoff has drained and streamflow comes from older groundwater sources (Cross, 1949). Understanding how groundwater contributions to streamflow are affected by different geologies is needed given that increasing temperatures have shown to reduce runoff efficiency (the ratio of runoff to precipitation) in watersheds near our study area (McCabe, Wolock, Pederson, Woodhouse, & McAfee, 2017).
Hydrologists have established numerous methods for separating hydrographs for well over 150 years. For example, Hall (1968) provides a thorough review on early attempts to model baseflow using graphical or mathematical approaches. Graphical approaches continue to be popular choices in data-limited regions as they only require observation of streamflow (Mugo & Sharma, 1999; Reitz, Sanford, Senay, & Cazenas, 2017). Collection of specific source-water chemical and / or isotopic composition allows for mass-balance (tracer-based) approaches for separating hydrographs into various flow components (Pinder & Jones, 1969). Mass-balance approaches are considered more reliable, particularly in mountainous regions, compared to graphical or mathematical interpretations because they integrate watershed-specific geochemical information (Miller, Susong, Shope, Heilweil, & Stolp, 2014). However, success using mass-balance methods to separate hydrograph components hinges on the ability to accurately assign concentrations to end-members. Combining graphical and tracer techniques offers a potentially cost-effective approach that has proven beneficial in identifying dominant underlying processes that would not have been possible with one method alone (Kronholm & Capel, 2015; McNamara, Kane, & Hinzman, 1997).
One example of where both tracer and graphical hydrograph separations might be complimentary is in the snowmelt-dominated and seasonally arid Snowy Range of southeast Wyoming, USA. In such systems, specific conductance (SC), a measurement of the ionic content in water, provides a strong natural tracer due to the large differences in concentration observed between peak flow and low flow conditions (Miller et al. , 2014). Further, SC can be measured continuously at relatively low costs using automated loggers, provides the best hydrograph separation results compared to individual chemical constituents (Caissie, Pollock, & Cunjak, 1996), and has been used successfully to separate baseflow in similar watersheds in the Southern Rocky Mountains (Caine, 1989; Liu, Williams, & Caine 2004; Miller et al., 2014; Rumsey, Miller, Susong, Tillman, & Anning, 2015). In addition, given the strong diurnal responses of Rocky Mountain systems to snowmelt “pulses”, graphical techniques might have utility to isolate changes in the fast flow pathways as the snowmelt season progresses and snowpack depletes (Buttle, Webster, Hazlett, & Jefferies, 2019). Diurnal streamflow cycles have been previously analyzed to gain process understanding on the hydrologic response to snowmelt (Caine, 1992; Kobayashi, 1986; Kurylyk & Hayashi, 2017; Loheide & Lundquist, 2009; Lundquist & Cayan, 2002; Lundquist & Dettinger, 2005; Mutzner et al. , 2015; Pellerin et al. , 2012; Woelber et al. , 2018).
The goal of this study is to assess the storage and release of snowmelt by quantifying the relative proportions of quickflow, throughflow, and baseflow contributions to streamflow. We do this for nine nested watersheds of the Snowy Range in southeastern Wyoming, USA by estimating how various sources contributed to discharge under different streamflow conditions throughout the year (i.e. rising limb, peak flow, falling limb, low flow). We used a common chemical mass-balance method to separate baseflow from direct runoff (Pinder & Jones, 1969). Direct runoff was further separated into throughflow and quickflow graphically based on naturally occurring diurnal cycles in streamflow that reflect rapid additions of direct runoff to streamflow induced by daily fluctuations in temperature. Further, we compared SC-Q hysteresis at both annual and daily timescales with watershed characteristics to better elucidate runoff processes driving the streamflow partitioning. Specifically, we address the following research questions for our seasonally-arid, snow-dominated systems:
  1. How do the sources of streamflow differ between seasonal and annual timescales?
  2. How do SC-Q relationships reflect the storage and release of snowmelt water at daily and annual scales?
  3. What are the dominant watershed properties governing streamflow partitioning and SC-Q relationships?

SITE DESCRIPTION AND DATA COLLECTION

2.1 Site description

The Snowy Range is the northernmost extent of the Medicine Bow Mountain Range within the Southern Rocky Mountains. We installed stream monitoring sites at nine nested locations along Libby Creek and the North Fork of the Little Laramie and their tributaries (Figure 1; Table 1). Data considered in this study were collected during three consecutive water years (October to September) covering 2016 through 2018. Our study watersheds range in elevation from 2570 – 3619 m (8431 – 11874 ft) and areal extent from 0.53 to 63.0 km2. The watersheds drain predominantly south and east towards Centennial, WY about 50 km east of Laramie, WY.
Bedrock geology underlying the watersheds is primarily a mix of metasedimentary and metavolcanic from the Early Proterozoic and glacial deposits from the Pleistocene-Holocene (Figure S1; Table S1; Wyoming State Geological Survey, 2014). Metasedimentary rocks from the Libby Creek Group are most prevalent at higher elevations consisting of metaclastic quartzite, politic and amphibole-schist, with minor traces of marble and conglomerate. Metasedimentary and metavolcanic rocks are common in GOLD100. Lower elevations of LIBB100 contain politic-schist, gneiss, and amphibolite with minor amounts of marble and granite. Glacial deposits consisting of unconsolidated, coarse-detrital gravel, boulders, and sand are common in the region. Landslide deposits occur in portions of NFLL100 containing intermixed glacial deposits, talus, and rock-glacier deposits. The surficial geology is dominated by glacial deposits lower elevations (Figure S2; Table 2; Case, Arneson, & Hallberg, 1998). Glacial deposits consist of scattered slope wash, residuum, grus, alluvium, colluvium, and landslide deposits. Glaciated bedrock is more common at higher elevations mixed with scattered shallow eolian, grus, colluvium, and alluvium deposits. Grus and residuum mixed with alluvium dominate the surficial geology of GOLD100, an outlier compared to the other study watersheds.
No detailed soil map has been constructed for our study watersheds. Munn and Arneson (1998) performed a low-order assessment and found soils to be highly variable, both regionally and locally with changes in slope aspect, slope position, climate, and geology. Haplocryalfs and Dystrocryepts occupy most of the area. Haplocryalfs tend to occur on low-relief areas while Dystrocryepts occupy forests on glacial deposits and on steep slopes. Cryaquepts and minor histisols are present in riparian areas (Munn & Arneson, 1998). Seismic refraction surveys performed in NONM100 indicate thin regolith at the top of hillslopes (1 m) thickening toward the bottom of the slope (3-4 m; Thayer et al., 2018). A porosity transition within the regolith indicated by seismic refraction surveys suggests the top layer of upper regolith has relatively high hydraulic conductivity capable of draining quickly. A novel bulk density optimization method performed in NONM100 supports the decline in porosity with depth, which creates conditions that favor shallow lateral flow through the upper regolith as a dominant streamflow generation mechanism (Fullhart, Kelleners, Chandler, McNamara, & Seyfried, 2019).
Evergreen forest is the main land cover classification in eight of our nine study watersheds ranging from 30 % (LIBB100) to 99 % (NONM100) coverage and is more common at lower elevations (Figure S3; Table S2; Homer & Fry, 2016). Shrub/scrub and grassland/herbaceous coverages with intermittent lakes and wetlands tend to dominate at higher elevations.

2.2 Data collection

2.2.1 Snow and climatological monitoring

The study watersheds receive most of their annual precipitation in the form of snow. Spring snowmelt, typically occurring in a relatively short period of time after the snowpack has become ‘ripe’, is a main driver of streams. Snow water equivalent (SWE) was recorded at six stations in our study watersheds (Figure 1). The Brooklyn Lake SNOTEL (sitenumber = 367) site is located within the boundaries of NASH200. According to 30-year daily median records (1981-2010), maximum SWE occurs on April 30 with a value of 59.2 cm. Mean water year total precipitation for Brooklyn Lake is 89.7 cm, which results in a mean water-year snow fraction (maximum SWE / total precipitation) of 66%. During the three-year period presented in this study, maximum SWE was remarkably consistent from 2016 to 2018 at the six measurement sites (Table S3). On average, maximum SWE ranged from about 32 cm at the lowest site (2684 m) to 72 cm at the highest site (3244 m) (Table S3).
Snow depth in 2016 was calculated at 0.5-m resolution in our study watersheds using LIDAR data based on two flights: snow-free in October 2014 and peak snow depth in April 2016 (Table S4). Monthly precipitation and temperature data from PRISM were aggregated to watershed boundaries to provide water-year (Oct. 1 – Sept. 30) mean watershed precipitation and temperature estimates (Table S4) (PRISM Climate Group, 2019). Water-year precipitation from PRISM (2016) was highly correlated to LIDAR-derived snow depth and elevation (Pearson’s r = 0.98, 0.99, respectively).

2.2.2 Stream monitoring

Stream stage (Level TROLL 500 Data Logger, In-Situ, Fort Collins, CO, USA; accuracy = ± 0.05%, resolution = ± 0.005%), electrical conductivity (ONSET HOBO U24, Onset Computer Corporation, Bourne, MA, USA; accuracy = 5 µS/cm, resolution = 1 µS/cm), and temperature (ONSET HOBO U24, Onset Computer Corporation, Bourne, MA, USA; accuracy = 0.1 °C, resolution = 0.01 °C) were measured instantaneously at 15 min intervals at nine gaging locations in the Snowy Range from 2016 to 2018 (Figure 1). Streamflow was estimated by creating rating curves at each gaging location that relate measured stream stage to repeated field measurements of discharge using a handheld electromagnetic water flow meter (OTT MF pro, OTT HydroMet, Loveland, CO, USA; accuracy = ± 2%). Electrical conductivity was measured by automated loggers and converted to specific conductance (SC) to account for differences in stream temperature following methodology from the manufacturer.
Stream stage and electrical conductivity loggers were installed each year prior to the main snowmelt period, but at slightly different dates for each gaging site. To account for variability in record length we calculated cumulative streamflow for each water year at each site and only analyzed the records after 3% of annual streamflow had occurred. This approximation allowed for consistent inter-site and inter-annual comparisons. Each year of analysis concluded at the end of the water year (September 30). Four periods were selected from the total hydrograph to assess dominant controls of runoff generation at different hydrologic conditions: the rising limb, peak flow, falling limb, and low flow conditions. The rising limb encompasses the time between the beginning of the snowmelt period of the hydrograph (defined in below in ‘Graphical quickflow separation’) and the date of maximum discharge. The peak flow time period was selected based on the highest 10% of streamflow values. The falling limb occurs between the date of maximum stream discharge and the last day of the snowmelt period of the hydrograph. It should be noted that based on these definitions peak flow partially overlaps with the rising and falling limbs. The low flow time period refers to the lowest 30% of streamflow recorded after the snowmelt portion of the hydrograph.

METHODS

3.1 Hydrograph separation

3.1.1 Tracer-based separation

We followed a chemical mass-balance methodology from Pinder & Jones (1969) for a tracer-based approach to separate baseflow and direct runoff from the total hydrograph:
\(C_{T}Q_{T}=\ C_{\text{BF}}Q_{\text{BF}}+C_{\text{DR}}Q_{\text{DR}}\)(1)
Where C and Q represent concentration and discharge, respectively, and subscripts T , BF , and DR refer to total, baseflow, and direct runoff, respectively. For this study, we use SC as the “concentration” of interest in Eq. (1) (e.g. Caine, 1989; Kobayashi, 1986; Miller et al. , 2014; Pilgrim, Huff, & Steele, 1979). The quantity of baseflow can be solved by rearranging Eq. (1):
\(Q_{\text{BF}}=\ Q_{T}\ \left[\frac{C_{T}-C_{\text{DR}}}{C_{\text{BF}}-C_{\text{DR}}}\right]\)(2)
Accurate interpretation of Eq. (2) requires that the following assumptions are made (Sklash & Farvolden, 1979; Buttle, 1994): (1) direct runoff and baseflow have different chemical composition, (2) direct runoff and baseflow chemical compositions are constant in time and space or the variability can be accounted for, (3) soil water contribution to streamflow is small or has similar composition of direct runoff, and (4) surface storage has a minor contribution to streamflow. Total stream concentration (CT ) and discharge (QT ) are obtained using in-situ measurements. The baseflow concentration (CBF ) is typically estimated from low flow stream conditions when total streamflow is assumed to be entirely derived from groundwater. We followed suggestions from Miller et al. (2014) and calculated CBF at each watershed by selecting the 99th percentile of the SC data to account for potential outliers. We assignedCDR a value of 21.6 µS/cm based on the mean value measured from 13 snow and 10 snowmelt samples. The SC of these samples ranged from 9.8 µS/cm to 42.1 µS/cm and the standard deviation across all 23 samples was 10.5 µS/cm.
Uncertainty was quantified following methods presented by Genereux (1998):
\(\text{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }W_{\text{BF}}=\sqrt{\left(\frac{f_{\text{BF}}}{C_{\text{DR}}-C_{\text{BF}}}W_{C_{\text{BF}}}\right)^{2}+\left(\frac{f_{\text{DR}}}{C_{\text{DR}}-C_{\text{BF}}}W_{C_{\text{DR}}}\right)^{2}+\left(\frac{-1}{C_{\text{DR}}-C_{\text{BF}}}W_{C_{T}}\right)^{2}}\)(3)
where WBF represents the uncertainty in baseflow discharge to streamflow at the 95% confidence interval,fBF and fDR are the fractions of total streamflow from baseflow and direct runoff, respectively, \(W_{C_{\text{BF}}}\) is the standard deviation of the highest 1% of measured SC concentrations multiplied by the t-value from the Student’s t distribution, \(W_{C_{\text{DR}}}\) is the standard deviation associated with observations of CDR (21.8 µS/cm) multiplied by the t-value, and \(W_{C_{T}}\) is the analytical error in the SC measurement multiplied by the t-value.

3.1.2 Graphical separation

After separating total streamflow into baseflow and direct runoff using the tracer-based approach, we applied a graphical hydrograph separation technique based on the presence of natural snowmelt-induced diurnal cycles to differentiate direct runoff into two components: throughflow and quickflow. Quickflow runoff represents water that moves most rapidly to the stream channel resulting in sharp rises of the daily hydrograph (Hewlett & Hibbert, 1967; Buttle et al. , 2019). Quickflow is generated during the snowmelt period as the result of daily fluctuation in temperature creating diurnal snowmelt pulses (Caine, 1992). For this reason, we constricted quickflow to only occur at times of peak snowmelt based on when the moving-average seven-day runoff was greater than 0.5 times mean water year daily runoff (i.e. the snowmelt period). This procedure approximated the time when snowmelt-induced diurnal cycles were present. The residual of the direct runoff water not contributing to quickflow was assumed to be throughflow, which represents delayed flow interacting more with the subsurface, resulting in a longer travel time.
A computer code was written to identify diurnal cycles by selecting times when instantaneous streamflow curvature was at a daily maximum. Quickflow volumes were calculated by subtracting the area created above the diurnal snowmelt cycles from the total hydrograph. The diurnal cycle (%) was then calculated by dividing the difference between daily maximum and minimum streamflow by the daily maximum streamflow for each day in the snowmelt period. Visual inspection of each hydrograph was performed to ensure accurate snowmelt period selection and quickflow volume calculation.

3.2 Hysteresis index

An index was calculated to quantify specific conductance-discharge (SC-Q) hysteresis at annual and daily scales during snowmelt-induced diurnal cycles. We adopted methodology from Lloyd, Freer, Johnes, and Collins (2016) who suggest rescaling the SC and Q data through minimum-maximum normalization to allow for an index that represents changing dynamics of an event and permits comparison between events:
\(\text{Normalized\ }Q_{i}=\ \frac{Q_{i}-\ Q_{\min}}{Q_{\max}-\ Q_{\min}}\)(4)\(\text{Normalized\ }\text{SC}_{i}=\ \frac{\text{SC}_{i}-\ \text{SC}_{\min}}{\text{SC}_{\max}-\ \text{SC}_{\min}}\)(5)
where Qi and SCi represent discharge and specific conductance at time i ,Qmin and SCmin are the minimum discharge and specific conductance, andQmaz and SCmax are the maximum discharge and specific conductance. The hysteresis index (HI) is calculated by subtracting the normalized SC on the falling limb from the normalized SC on the rising limb for a particular flow percentile:
\(\text{HI}_{Q_{i}}=\ \text{Normalized\ SC}_{\text{i\ }_{\text{Rising\ limb}}}-\ \text{Normalized\ SC}_{\text{i\ }_{\text{Falling\ limb}}}\)(6)
At the annual scale, we calculated a hysteresis index for the 10th, 25th, and 50th percentile of discharge values. For example, HI10 was calculated by subtracting the tenth percentile of discharge occurring during the falling limb from the tenth percentile of discharge occurring during the rising limb.
The HI varies between -1 and 1, where the larger the absolute value indicates a more circular hysteresis and the sign represents the direction of the loop (positive for clockwise, negative for counter-clockwise) (Lloyd et al. , 2016). To avoid spikes in the data affecting the hysteresis index, we averaged the normalized SC for values within 1% of each annual hysteresis index percentile (i.e. the 10th percentile of normalized discharge was expanded to include the 9th through 11thpercentile values).
Similar methodology was applied to create a hysteresis index at the daily scale. For this analysis, each date was rescaled from noon on the particular date to noon on the subsequent date because snowmelt-induced diurnal cycles begin rising in the late afternoon following a delay from the time of peak daily solar maximum. The daily hysteresis index was calculated for the 25th, 50th, and 75th percentile of normalized discharge values. Due to many fewer data points compared to the annual datasets, we expanded the buffer by averaging the normalized SC for values within 10% of each daily hysteresis index percentile (i.e. the 25thpercentile of normalized discharge was expanded to include the 15th through 35th percentiles).

4. RESULTS

4.1 Tracer-based hydrograph separation

Baseflow contributions to streamflow varied spatially and temporally but were greatest at low flow conditions for all watersheds (Figure 2; Table 3). On average, total baseflow contributions to streamflow in the study watersheds ranged from 22.1% (GOLD100) to 45.5% (NONM100). This range of baseflow contribution is similar to other studies in the Upper Colorado River Basin (Miller et al., 2014; Rumsey et al., 2015). During the rising limb portion of the annual hydrograph, baseflow contributions ranged from 14.2% (GOLD100) to 33.9% (NFLL200). Baseflow contributions were lowest during peak flow conditions, ranging from 10.8% (GOLD100) to 29.0% (NASH100). During the falling limb, proportions of baseflow increased slightly and ranged from 13.0% (GOLD100) to 32.8% (NASH100). Baseflow contributions to streamflow were greatest at low flow conditions and ranged from 65.6% (NASH100) to 93.4% (GOLD100) (Table 3). The close proximity of NASH100 to Brooklyn Lake, which is one of the largest lakes in the region and is located upstream of NASH100, likely explains the smaller contribution of baseflow to total streamflow at low flow conditions. Uncertainty was greater when proportions of direct runoff were greater due to a smaller, and more variable, number of snowmelt/snowpack samples used to define CDR in Eq. 3 (Table 3).
It is important to emphasize that most of the total streamflow in our study watersheds was generated during a relatively brief period associated with peak snowmelt. The highest 10% of recorded streamflow values (i.e. peak flow conditions) were responsible for 37% (NASH200) to 47% (GOLD100, LIBB200, LIBB400) of total streamflow in our study watersheds (Table 4). During this time, streamflow was overwhelmingly derived from direct runoff (Table 3). In contrast, the lowest 30% of streamflow values recorded after peak flow (i.e. low flow conditions) accounted for only 1.4% (LIBB100, LIBB200) to 3.0% (NFLL100) of total streamflow (Table 4) and were primarily sourced from baseflow (Table 3). However, baseflow contributions were responsible for a majority of total streamflow generation on any particular day outside of the relatively short snowmelt period.

4.2 Graphical hydrograph separation

Direct runoff produced during diurnal snowmelt cycles (Figure 3) resulted in quickflow contributions to total streamflow that ranged from 7% (NASH100) to 14% (LIBB100, LIBB200) annually (Table 4). Quickflow was responsible for 15.4% (NASH100) to 27.7% (LIBB100) of direct runoff during the snowmelt period and 11.4% (NASH100) to 22.1% (LIBB200) of direct runoff annually for our study watersheds (Table 5). When present, the mean snowmelt-induced diurnal streamflow cycle ranged from 19.7% (NASH100) to 33.5% (LIBB100) of total streamflow during the snowmelt period (Table 5). The mean magnitude of the diurnal cycle was significantly (p < 0.05) greater for Libby Creek (LIBB100, LIBB200, LIBB400) compared to Nash Fork and the North Fork of the Little Laramie (NASH100, NASH200, NFLL100, NFLL200). Throughflow was the main contributor to streamflow at all times except during low flow conditions (Figure 4; Tables 4 & 5) when baseflow was dominant. The close proximity of NASH100 to Brooklyn Lake, which drains directly to the Nash Fork main channel, likely dampened the amplitude of the diurnal fluctuations measured in the stream and resulted in the lowest proportion of quickflow of all the watersheds considered.

4.3 Annual SC-Q hysteresis

Discharge (Q) and SC generally mirror each other in our study watersheds with the highest Q values corresponding with the lowest SC (and vice versa) implying dilution from snowmelt water with low SC at high Q (Figure 5). Time lags between SC and Q at an annual scale result in hysteresis which varies systematically across the different runoff components. In NONM100, LIBB100, LIBB200, and LIBB400, SC-Q relationships have a negative, concave, and anticlockwise hysteretic behavior at annual timescales, resulting in negative annual hysteresis indices (Figure 5d; Table 6), implying that SC concentration is greatest in baseflow, followed by throughflow, and then quickflow (Evans & Davies, 1998). The lack of annual hysteresis observed at GOLD100, NASH100, and NFLL100 implies little to no time lag between annual peak SC and annual peak Q suggesting relatively inert subsurface material properties that result in little chemical evolution as water travels through the soil and /or regolith layers. Thus throughflow and quickflow result in similar geochemical alteration in these watersheds (Figure 5c).
All watersheds had a negative relationship between SC and Q on an annual scale, implying recently melted snow (i.e. from the current snowmelt season) is responsible for the majority of streamflow generation in a given year. Seasonal snowmelt has been shown to be the dominant streamflow generation mechanism in the western United States (e.g. Li et al. , 2017) but differs from what has been shown using isotope tracers in snowmelt-dominated watersheds located in lower-elevation, humid areas (e.g. Buttle, 1994; Laudon, Sjoblom, Buffam, Seibert, & Morth, 2007) where snowmelt mainly recharges groundwater and displaces previously stored older water.

4.4 Daily SC–Q hysteresis

When distinct snowmelt-induced diurnal streamflow cycles are present, Q and SC often demonstrated daily hysteresis. The shape of the diurnal streamflow cycle was typically front-loaded, where streamflow rises rapidly followed by more gradual decline (Figure 6a). The opposite was true for SC in streamflow which exhibited a sharp drop followed by a gradual increase (Figure 6a). This phenomenon was also demonstrated by Kobayashi (1986) during snowmelt in the headwaters of the Uryu River, Japan. In most observed cases for our study watersheds, SC was lower on the falling limb of diurnal cycles than on the rising limb, resulting in clockwise hysteresis and positive daily hysteresis indices (Figure 6b-c; Table 6). The daily hysteresis index was largest on the rising limb of the annual hydrograph for most watersheds (Figure 6b; GOLD100, NASH100, NASH200, NFLL200, NONM100) and became smaller throughout the snowmelt period during the falling limb of the annual hydrograph (Figure 6c). The opposite was true for Libby Creek watersheds (LIBB100, LIBB200, LIBB400), where the hysteresis index became larger as the snowmelt period progresses.