Monitoring snowpack outflow volumes and their isotopic composition to better understand streamflow generation during rain-on-snow events

Rain-on-snow (ROS) events in mountainous catchments can cause enhanced snowmelt, leading to an increased risk of destructive winter floods. However, due to differences in topography and forest cover, the generation of snowpack outflow volumes and their contribution to streamflow are spatially and temporally variable during ROS events. In order to adequately predict such flood events with hydrological models, an enhanced process understanding of the contribution of rainwater and snowmelt to stream water is needed.

In this study, we monitored and sampled snowpack outflow with fully automated snowmelt lysimeter systems installed at three different elevations in a pre-Alpine catchment in central Switzerland. We measured snowpack outflow volumes during the winters of 2017 and 2018, as well as snowpack outflow isotopic compositions in winter 2017. Snowpack outflow volumes were highly variable in time and space, reflecting differences in snow accumulation and melt. In winter 2017, around 815 mm of snowpack outflow occurred at our reference site (grassland 1220 m a.s.l. – metres above sea level), whereas snowpack outflow was 16 % less at the nearby forest site (1185 m a.s.l.), and 62 % greater at another grassland site located 200 m higher (1420 m a.s.l.). A detailed analysis of 10 ROS events showed that the differences in snowpack outflow volumes could be explained mainly by rainfall volumes and initial snow depths.

The isotope signals of snowpack outflow were more damped than those of incoming rainwater at all three sites, with the most damped signal at the highest elevation site because its snowpack was the thickest and the residence times of liquid water in its snowpack were the longest, thus enhancing isotopic mixing in the snowpack. The contribution of snowpack outflow to streamflow, estimated with an isotope-based two-component end-member mixing model, differed substantially among the three lysimeter sites (i.e. between 7±4 and 91±21 %). Because the vegetation in our study catchment is a mixture of grassland and forest, with elevations ranging from 1000 to 1500 m a.s.l., our site-specific hydrograph separation estimates can only provide a range of snowpack outflow contributions to discharge from different parts of the study area. Thus, the catchment-average contribution of snowpack outflow to stream discharge is likely to lie between the end-member mixing estimates derived from the three site-specific data sets. This information may be useful for improving hydrological models in snow-dominated catchments.

1 Introduction

Over the past 50 years, rain-on-snow (ROS) events have become more frequent in snow-dominated catchments, because rising global mean air temperatures have led to greater fractions of winter precipitation falling as rain instead of snow (Barnett et al., 2005; Beniston and Stoffel, 2016; Hartmann et al., 2013; Stewart, 2009). In Switzerland, the mean air temperature is predicted to increase by up to 1.6 ^∘C by 2050 (Swiss Academies of Sciences, 2016), and in mountain regions, including the Alps, climate warming is expected to make rain-on-snow events more prevalent in the future, and raise the elevation of the rain-to snow transition zone (McCabe et al., 2007; Surfleet and Tullos, 2013; Ye et al., 2008).

Rain on snow can either be retained in the snowpack or it can enhance snow melt; therefore, the snowpack can either reduce or amplify the volume of water reaching the ground surface compared to snow-free conditions (Kattelmann, 1987; Lee et al., 2010a). In the past, some ROS events that caused enhanced snowmelt have led to severe floods (e.g. Garvelmann et al., 2015; Kroczynski, 2004; MacDonald and Hoffman, 1995; Marks et al., 1998; Sui and Koehler, 2001; Wever et al., 2014). Although catchment models have been applied in the past to predict flood responses during ROS events, these model simulations can be highly uncertain (McCabe et al., 2007; Rössler et al., 2014) because snowpack outflow (snowmelt or a mixture of rainwater and snowmelt) is not generated homogeneously at the catchment scale (Berris and Harr, 1987; Würzer et al., 2016). Peak flows caused by ROS events result from a complex interplay of processes that mainly depend on the initial snowpack properties, rainfall characteristics and energy fluxes (Colbeck, 1977; Garvelmann et al., 2014; Würzer et al., 2016), as well as antecedent catchment wetness.

Snowpack properties such as depth, density and snow water equivalent (SWE) can vary spatially and temporally across the catchment landscape. Additionally, wind drift, landscape topography (i.e. slope, elevation and aspect) and vegetation cover (i.e. forest or grassland) affect the snowpack properties (Marks et al., 1998; Molotch et al., 2011; Stähli et al., 2000). Higher elevations are generally associated with greater snow depths due to higher precipitation rates and lower air temperatures (Beniston et al., 2003; Stewart, 2009). Compared with open grassland, forested landscapes tend to have shallower snow depths due to enhanced canopy interception of snowfall (Berris and Harr, 1987; López-Moreno and Stähli, 2007; Stähli and Gustafsson, 2006). Furthermore, water flow paths within the snowpack can be highly variable, so that calculating or measuring the snowpack outflow can be challenging (Eiriksson et al., 2013; Kattelmann, 2000; Rücker et al., 2019; Unnikrishna et al., 2002; Webb et al., 2018; Yamaguchi et al., 2018).

A detailed understanding of snowpack outflow generation is needed at both the plot scale and the catchment scale to make runoff predictions during ROS events more reliable (DeWalle and Rango, 2008; Marks et al., 1998; Šanda et al., 2014). To track the heterogeneous contribution of snowpack outflow to streamflow during ROS events, environmental tracers can be used. Stable water isotopes (¹⁸O and ²H) may be particularly useful as they allow streamflow to be separated into isotopically distinct sources (Klaus and McDonnell, 2013). Thus, if snowpack outflow is isotopically distinguishable from catchment groundwater storage, its relative contributions to streamflow during ROS events can be quantified through two-component isotope-based hydrograph separation (IHS).

In some studies, the isotopic composition of bulk snow or of individual snow layers has been used as a proxy for snowmelt isotopic composition (Cooper et al., 1993; Dinçer et al., 1970; Huth et al., 2004; Maulé et al., 1994; Sueker et al., 2000). The isotopic composition of bulk snow is known to be variable in time and space, depending on catchment characteristics such as latitude, exposure and elevation gradients (Dietermann and Weiler, 2013), as well as the structure of the forest canopy (Koeniger et al., 2008). Snowfall intercepted by the forest canopy is subject to sublimation, which is the main cause of the isotopic enrichment of winter throughfall (Claassen and Downey, 1995; Koeniger et al., 2008; Stichler, 1987). Although the spatial variations in the isotopic compositions of bulk snow and snowmelt are likely to be similar (Dietermann and Weiler, 2013), estimated meltwater contributions to streamflow can be significantly different when using bulk snow instead of snowmelt as an end-member in IHS (Moore, 1989). Numerous studies have found that IHS in snow-dominated catchments is less uncertain when snowmelt is collected by grab sampling (Obradovic and Sklash, 1986; Penna et al., 2017), with melt pans (Bales and Williams, 1993; Taylor et al., 2001) or with snowmelt lysimeters (Beaulieu et al., 2012; Buttle, 1994; Hooper and Shoemaker, 1986; Laudon et al., 2002; Schmieder et al., 2016; Shanley et al., 1995b; Unnikrishna et al., 2002; Wels et al., 1991). Some snowmelt lysimeter systems facilitate water sampling at regular time intervals (i.e. daily or sub-daily), which is recommended because the isotopic composition of snowmelt can be highly variable over time. This variability is caused by isotopic fractionation in the snowpack during phase changes (i.e. freezing/melting, sublimation/recrystallization/condensation; Judy et al., 1970; Lee et al., 2010a; Schmieder et al., 2016; Sokratov and Golubev, 2009; Stichler et al., 2001; Taylor et al., 2001, 2002a; Unnikrishna et al., 2002) and by ROS events when isotopically distinct rainwater percolates and mixes with the snowpack (Berman et al., 2009; Herrmann et al., 1981; Juras et al., 2016; Shanley et al., 1995a). Furthermore, isotopic exchange and redistribution in the snowpack can cause the snowpack outflow to be isotopically different from incoming rainfall (Judy et al., 1970; Lee et al., 2010a, b; Taylor et al., 2001).

To the best of our knowledge, only two scientific studies have estimated the contribution of snowpack outflow to streamflow during ROS events by IHS (Buttle et al., 1995; Maclean et al., 1995). However, these studies used only one or two ROS events that occurred before and during snowmelt. So far, the effects of forest cover and elevation on the generation of snowpack outflow and snowmelt contribution to streamflow have not been investigated during ROS events.

To fill this research gap, we monitored snowpack outflow and its isotopic composition in the pre-Alpine Alptal catchment in central Switzerland. Three snowmelt lysimeter sites were located between 1200 and 1400 m a.s.l., altitudes at which precipitation frequently shifts between snowfall and rainfall (Stewart, 2009). One of the lysimeter systems was installed under forest canopy, and the other two systems were installed in open grassland. We measured snowpack outflow volumes every 10 min during the winters of 2017 (1 January–7 May 2017) and 2018 (1 November 2017–6 April 2018), as well as δ¹⁸O and δ²H in snowpack outflow at daily intervals during the winter of 2017.

We hypothesise that snowpack outflow generation during ROS events is spatially and temporally variable at the catchment scale, depending on elevation and vegetation cover. Specifically, we will address the following research questions:

• What role do rainfall characteristics and initial snowpack properties play in the variability of snowpack outflow volumes?

• What is the relative contribution of snowpack outflow to streamflow during rain-on-snow events?

• What is the spatial and temporal variability of snowpack outflow contributions to streamflow?

• How does the choice of the event-water end-member (rainwater or snowpack outflow) affect the results of hydrograph separations?

2 Methodology

2.1 Field site

Field work was conducted in the southern part of the 47 km² Alptal catchment, 40 km south of Zürich in the northern pre-Alps in central Switzerland (Fig. 1). The Erlenbach catchment is a 0.7 km² tributary of the Alp River with a streamflow gauging station and a meteorological station.

Figure 1 The Erlenbach catchment (red outline) in the southern part of the Alptal Valley, showing the distribution of vegetation (grassland and forest), as well as the locations of the three snowmelt lysimeter systems (MG: mid-elevation grassland site; MF: mid-elevation forest site; HG: high-elevation grassland site). A meteorological station is located near the MG site. At the Erlenbach catchment outlet, river discharge and precipitation (snowfall and rainfall) rates are measured, along with stable water isotopes in stream water and precipitation.

In the Erlenbach catchment, winter precipitation is dominated by snowfall, which accounts for up to one-third of the total annual precipitation of roughly 2300 mm yr⁻¹ (Feyen et al., 1999; van Meerveld et al., 2018). The annual average air temperature is 6 ^∘C with a distinct seasonal cycle (−2 ^∘C in February and 17 ^∘C in August; Feyen et al., 1999). The Erlenbach catchment covers an altitude range from 1080 to 1520 m a.s.l., and precipitation events frequently shift between rainfall and snowfall during the winter season (Stähli and Gustafsson, 2006). The bedrock of the Erlenbach catchment is dominated by tertiary flysch, overlain by shallow soils with low permeability (Burch et al., 1996; Fischer et al., 2015). The landscape is characterised by coniferous forests (53 %), grassland (25 %) and a mixture of both (22 %) (Burch et al., 1996; Fischer et al., 2015; Keller, 1990).

Two snowmelt lysimeter systems, hereafter referred to as the MG site (mid-elevation grassland) and MF site (mid-elevation forest), were located in the Erlenbach catchment at altitudes of 1216 and 1185 m a.s.l., respectively. The MG and MF sites were installed at field sites maintained by the Swiss Federal Institute for Forest, Snow and Landscape Research (WSL) where power supply was available. A third snowmelt lysimeter system was installed at a high-elevation grassland site at 1405 m a.s.l. (HG site), outside of the Erlenbach catchment at a location with power supply. The measurements from the HG site were assumed to be representative of the 1400 m elevation zone of the Erlenbach catchment; however, we acknowledge that the HG site was located on a flat, slightly north-facing plateau, whereas the Erlenbach catchment is characterised by a sequence of flat plateaus and west-facing slopes. The terrain at the meteorological station and the MG site was relatively flat, whereas the MF site was located on slightly sloped terrain facing west.

The forest at the MF site is dominated by Picea abies and Abies alba, whereas the MG and HG sites are located on grassland. The MG and the MF sites were located 250 m from one another with an elevation difference of only 30 m. In the following analysis, we use the MG site as a reference site because it was located in close proximity to the meteorological station in the centre of the Erlenbach catchment, which is practical for comparing the effects of vegetation cover (MF vs. MG) and elevation (HG vs. MG) on snowpack outflow generation. At the meteorological station, we measured air temperature (107 thermistor probe, Campbell Scientific, Loughborough, Great Britain) every minute, as well as snow depth (ultrasonic depth sensor, Judd Communications, Salt Lake City, Utah, USA) and precipitation (Lambrecht meteo GmbH, rain gauge 15189, Göttingen, Germany) at a 10 min temporal resolution (Stähli and Gustafsson, 2006). Additionally, air temperature was measured at the HG and MF sites every minute (107 thermistor probe, Campbell Scientific, Loughborough, Great Britain). At the Erlenbach catchment outlet, river discharge was recorded every 10 min (Burch et al., 1996) and a heated rain gauge (HOBO rain gauge, metric data logger, RG3-M, Bourne, MA 02532, USA) measured precipitation at 10 min intervals. The rain gauge at the meteorological station was not functioning during the period from 18 February to 11 March 2017; thus, we filled the data gap with measurements from the heated rain gauge installed near the catchment outlet (von Freyberg et al., 2018b).

2.2 A snowmelt lysimeter system for measuring and sampling snowpack outflow

The snowmelt lysimeter system was designed to measure the natural snowpack outflow and to collect samples for the analysis of stable water isotopes (see details in Rücker et al., 2019). The MG lysimeter system site was installed in March 2016, whereas the HG and MF lysimeter systems were installed in October and December 2016, respectively. The design of the HG and MF lysimeter systems has been slightly improved to facilitate easier access to the technical components (e.g. tipping bucket) during snow-rich periods. Each of the three snowmelt lysimeter systems (MG, MF and HG) consisted of three individual funnels (42.15 cm diameter, 1400 cm² area and 5.9 cm rim height) that were installed into the soil so that they collected the daily snowpack outflow at the snowpack–soil interface. From each individual lysimeter funnel, the snowpack outflow ran through a silicon rubber tube to a 10 L collection vessel. Every day at 05:40 UTC + 1 h, an automatic water sampler (Maxx P6L – vacuum system, Maxx GmbH, Rangendingen, Germany) pumped up to 300 mL of water sample from the collection vessel into a dry 1 L HDPE autosampler bottle. After that, pinch valves at the lower outlet of the collection vessel were opened for 10 min to drain the remaining liquid. When the snowpack outflow volumes reached the 10 L storage capacity of the water vessel, an additional bulk sample of each water vessel was collected by manually operating a pump cycle. A whole pumping and rinsing cycle took around 20 min so that the starting time for the collection of the next water sample was set to 06:00 UTC + 1 h, i.e. a 1-day sample contained the cumulative snowpack outflow that left the snowpack between 06:00 and 05:40 UTC + 1 h of the following day. The filled sampling bottles in the automatic samplers remained open until they were replaced with dry bottles once a week. One additional open sample bottle filled with a 400 mL water sample of known isotopic composition was placed in the automatic water sampler each week to test whether evaporative fractionation occurred during the 1-week sampling period. We did not find a substantial isotopic enrichment effect in these open sample bottles and thus assume that our sampling set-up (automatic water sampler inside a protection hut) provided sufficient protection against evaporative fractionation of the snowpack outflow samples.

Before the snowpack outflow reached the collection vessel, its volume was measured via a tipping bucket mechanism at 5 mL volume increments (each corresponding to 0.036 mm of outflow). The tipping bucket mechanism was installed either directly below each individual lysimeter funnel (MG) or between the end of the silicon rubber tube and the collection vessel (MF and HG). The arrangement of the tipping bucket was changed for the MF and HG systems, so that it could easily be replaced or repaired if necessary. As the tipping bucket mechanisms of the snowmelt lysimeter systems were slightly adapted to the local properties at each field site, the average measurement uncertainties of the snowmelt outflow volume measurements were determined from replicate measurements of known water volumes poured into each funnel. The average measurement uncertainties were 15 %, 7.5 % and 10 % at the HG, MG and MF sites, respectively.

In an earlier study that evaluated the performance of the snowmelt lysimeter design at the MG site, we found that the three individual funnels registered highly variable snowpack outflow volumes, thus reflecting temporal and spatial variability of the snowmelt processes at the plot scale (Rücker et al., 2019). In the analysis presented here, we calculated the site-averages of snowpack outflow volumes collected with the three individual funnels at each lysimeter site. To express the uncertainty of these measurements at the plot scale, we calculated the standard errors of these snowpack outflow site-averages considering both the relative measurement uncertainty of the tipping bucket and the spatial variability of snowpack outflow generation. Measurements at 10 min resolution were aggregated to daily resolution for the time period from 06:00 to 05:40 UTC + 1 h of the following day, which corresponds to the daily time interval that is used by the Federal Office of Meteorology and Climatology (MeteoSwiss). The 10 min data were aggregated to hourly data over the periods from HH:40 to (HH + 1):30, with HH denoting the hour of the measurement.

To prevent freezing in the tubes or in the tipping bucket mechanisms, a heating cable (Pentair, Raychem, BZV self-regulating heating band, Wisag, Fällanden, Switzerland) was attached to the silicon tubes of the MG and HG lysimeter systems in December 2017. In addition, a 12-W heating patch (110 mm × 77 mm) was attached next to the tipping bucket mechanism below the lysimeter funnel of the MG lysimeter system. The MF site was not equipped with a heating patch or a heating cable, as freezing was less problematic in the forest than at the open grassland sites. During winter 2018, the lysimeter funnels at the MG and HG sites were damaged due to the heavy snow cover. Thus, snowpack outflow volumes for the 2018 winter period were only measured at the MF site.

At all three lysimeter sites, snow depths were measured with stakes located next to the individual lysimeter funnels. A webcam at each site recorded a picture of the stakes every hour, and we used one image taken between 09:53 and 11:23 UTC + 1 h to estimate the mean daily snow depth at each lysimeter site. At the meteorological station near the MG site, the snow depth sensor provided additional information about the maximum hourly snow depth, which was used to validate the mean daily snow depths obtained from the webcam images. The snow depth sensor was not functioning during the period from 8 March 2017 16:30 UTC + 1 h to 9 March 2017 00:50 UTC + 1 h.

At the MG and HG sites, snow surveys were carried out at weekly intervals. In addition, snow surveys at monthly intervals were undertaken at the MG and MF sites, meaning that the MG site was surveyed twice in the same week roughly once a month. During each survey, snow depth and bulk snow density were measured along a ca. 30 m snow course with a snow tube (diameter 50 mm, length 1.2 m) to determine the snow water equivalent (SWE) (Stähli et al., 2000; Stähli and Gustafsson, 2006). Additionally, soil conditions were characterised (frozen or not frozen) with a 1 cm diameter aluminium stake. At the MG site, the SWE prior to a ROS event was estimated by the product of the actual snow depth recorded by the snow depth sensor and the bulk density derived from the most recent snow survey. At the MF site, only two surveys were carried out during the 2017 winter season.

2.3 Identification of rain-on-snow (ROS) events

ROS events during winter 2017 were identified based on the temperature and rainfall measurements at the meteorological station and the snow properties the MG site. For winter 2018, measurements from the meteorological station and the MF site were used because the MG site was not functional. Criteria for identifying ROS events are often arbitrary and largely depend on the research purpose and data availability (Mazurkiewicz et al., 2008; McCabe et al., 2007; Würzer et al., 2016). In our study, we identified ROS events based on following criteria: rainfall rates greater than 0.1 mm h⁻¹, a total rainfall volume of at least 20 mm within 12 h, air temperatures above 0 ^∘C and an initial snowpack depth of at least 10 cm. As the beginning and the end of the snowpack outflow generation were generally delayed relative to rainfall, we defined the end of a ROS event as the point in time when the snowpack outflow volumes receded to the pre-event flow rates or when a new rainfall event started. We applied the same ROS timing for all three sites based on the measurements from the MG site during winter 2017 and from the MF site during winter 2018. To compare the volumes of snowpack outflow with those of incoming rainfall during the ROS events, we aggregated snowpack outflow volumes over these ROS event periods. We calculated the “snowpack water budget” of a ROS event as the difference between aggregated snowpack outflow volume and rainfall volume.

2.4 Sample collection and isotope analysis

We measured stable water isotopes (¹⁸O and ²H) in stream water, precipitation (snowfall and rainfall), snowpack outflow and melted bulk snow samples. All of these samples were collected once each day except for the bulk snow, which was collected roughly once each week during the snow surveys. At each of the three snowmelt lysimeter sites, daily composite samples of snowpack outflow from each of the three individual funnels were collected using autosamplers for subsequent isotope analysis (Sect. 2.2). As shown in Rücker et al. (2019), the isotopic compositions of snowpack outflows from the three individual funnels at each site were very similar; thus, we averaged the isotope values from the three individual lysimeter funnels at each site and expressed their spatial variability through the standard error of the mean.

Composite stream water samples were collected using an automatic water sampler (6712 full-size portable sampler, Teledyne Isco, Lincoln, NE, USA) that pumped 100 mL of stream water into a dry 1 L HDPE bottle four times a day (at 05:40, 11:40, 17:40 and 23:40 UTC + 1 h). Every 3 weeks, the filled autosampler bottles were replaced with empty ones. Stable water isotopes in precipitation were measured at hourly resolution at the Erlenbach catchment outlet with a wavelength-scanned cavity ring-down spectrometer (CRDS; model L2130-I; Picarro Inc., Santa Clara, CA, USA) coupled to a Picarro Inc. continuous water sampler module. The analytical uncertainty of the analyser was 0.09 ‰ for δ¹⁸O and 0.21 ‰ for δ²H (von Freyberg et al., 2018a). More details about the high-frequency sampling approach can be found in von Freyberg et al. (2017, 2018b). To obtain mean daily precipitation isotope values, the hourly precipitation isotope measurements were volume-weighted with hourly precipitation volumes for the time period from 05:40 to 05:39 UTC + 1 h of the following day. Whenever possible, this volume-weighting was based on precipitation measurements from the meteorological station at the MG site; however, during a period of instrument failure (18 February–11 March 2017), precipitation measurements from the heated rain gauge at the catchment outlet were used instead (i.e. for ROS events 4 and 5).

During the snow surveys, bulk snow was collected from the entire snow profile close to the lysimeter sites with a snow tube (diameter 50 mm, length 1.2 m) and transferred to a HDPE plastic bag ($\mathrm{300}\times \mathrm{500}\times \mathrm{0.1}$ mm, Plasti-Pac Zürich AG, Zürich, Switzerland), which was sealed immediately. At the MG site, one or two bulk snow samples were collected near the lysimeter funnels on two different days every week during the snow surveys in winter 2017. Because snow depths were generally greater at the HG site, three bulk snow samples were collected here during each weekly snow survey and the isotopic compositions of these three bulk snow samples were averaged.

All water samples that were collected in the field were stored in sealed bottles and refrigerated at 4 ^∘C in the laboratory until sample preparation. Frozen samples were melted in the laboratory at room temperature. All samples were filtered through 0.45 µm Teflon filters (DigiFilter micron Teflon, S-Prep GmbH, Überlingen, Germany) and filled into 2 mL glass vials with silicon seals. Isotopic compositions of all water samples were measured with an LGR IWA-45EP off-axis integrated cavity output spectrometer (ABB Los Gatos Research, San Jose, California, USA) at the laboratory of the Swiss Federal Institute for Forest, Snow and Landscape Research (WSL). Isotopic abundances are reported using the δ notation relative to the IAEA-VSMOW-II and SLAP-II standards. The analytical uncertainties of the analyser were 0.21 ‰ for δ¹⁸O and 0.37  ‰ for δ²H, which were estimated from replicate check-standard measurements within the same batch.

2.5 Two-component isotope-based hydrograph separation

We used two-component isotope-based hydrograph separation (IHS) to estimate the relative contribution of snowpack outflow to catchment streamflow. The calculation of the relative contribution of snowpack outflow to streamflow (F_spo) was based on the conventional mass balance equation of Pinder and Jones (1969):

$$\begin{array}{}\text{(1)}& F_{\mathrm{spo}}={\displaystyle \frac{V_{\mathrm{spo}}}{V_{\mathrm{S}}}}={\displaystyle \frac{C_{\mathrm{S}}-C_{\mathrm{pe}}}{C_{\mathrm{spo}}^{*}-C_{\mathrm{pe}}}},\end{array}$$

where V_spo denotes the volume of snowpack outflow in streamflow (V_S), and C_S, C_pe and $C_{\mathrm{spo}}^{*}$ denote the tracer signatures in stream water, pre-event water and snowpack outflow, respectively. We used daily time steps for all calculations. The pre-event tracer signature (C_pe) was represented by the isotopic composition of stream water of the day prior to the ROS event of interest, and was assumed to be constant during the event (Blume et al., 2008; Penna et al., 2016). The isotopic composition of snowpack outflow at day i was calculated as the incremental volume-weighted mean using the measured volumes of snowpack outflow or rainfall since the beginning of the event on day j (McDonnell et al., 1990):

$$\begin{array}{}\text{(2)}& C_{\mathrm{spo},i}^{*}={\displaystyle \frac{\sum _{k=j}^i\left(V_{\mathrm{spo},k}{C}_{\mathrm{spo},k}\right)}{\sum _{k=j}^i{V}_{\mathrm{spo},k}}}.\end{array}$$

Equation (3) quantifies the standard error of $C_{\mathrm{spo}}^{*}$, SE($C_{\mathrm{spo}}^{*}$), which is a combination of the measurement uncertainty of the isotope analyser (first summand on the right side) and the spatial variability of the snowpack outflow volumes collected by the three individual funnels of each lysimeter system (second summand on the right side):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\mathrm{SE}\left(C_{\mathrm{spo},i}^{*}\right)=\\ {\displaystyle }& {\displaystyle }\sqrt{{\sum }_{k=j}^i{\left({\displaystyle \frac{V_{\mathrm{spo},k}}{\sum _{k=j}^i{V}_{\mathrm{spo},k}}}\mathrm{SE}\left(C_{\mathrm{spo},k}\right)\right)}^{\mathrm{2}}+{\sum }_{k=j}^i{\left[\left({\displaystyle \frac{C_{\mathrm{spo},k}}{\sum _{k=j}^i{V}_{\mathrm{spo},k}}}-{\displaystyle \frac{\sum _{k=j}^i\left(C_{\mathrm{spo},k}{V}_{\mathrm{spo},k}\right)}{{\left(\sum _{k=j}^i{V}_{\mathrm{spo},k}\right)}^{\mathrm{2}}}}\right)\mathrm{SE}\left(V_{\mathrm{spo},k}\right)\right]}^{\mathrm{2}}}.\end{array}\end{array}$$

In Eq. (3), SE(C_spo) is the standard error of the isotope data, which is assumed to be the measurement uncertainty of the isotope analyser, and SE(V_spo) is the standard error of the mean of the three individual snowpack outflow volumes measured with the individual lysimeter funnels. Using Gaussian error propagation (Genereux, 1998), the uncertainty of F_spo was estimated as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\mathrm{SE}\left(F_{\mathrm{spo}}\right)=\\ {\displaystyle }& {\displaystyle }\sqrt{{\left({\displaystyle \frac{-\mathrm{1}}{\left(C_{\mathrm{pe}}-C_{\mathrm{spo}}^{*}\right)}}\mathrm{SE}\left(C_{\mathrm{S}}\right)\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{C_{\mathrm{S}}-C_{\mathrm{spo}}^{*}}{{\left(C_{\mathrm{pe}}-C_{\mathrm{spo}}^{*}\right)}^{\mathrm{2}}}}\mathrm{SE}\left(C_{\mathrm{pe}}\right)\right)}^{\mathrm{2}}+{\left(\mathrm{SE}\left(C_{\mathrm{spo}}^{*}\right)\right)}^{\mathrm{2}}}.\end{array}\end{array}$$

We assume that the standard errors of C_s and C_pe, SE(C_s) and SE(C_pe), are equivalent to the measurement uncertainty of the isotope analyser.

As we only measured snowpack outflow volumes and their isotopic compositions at three locations and not across the entire Erlenbach catchment, we cannot reliably estimate the catchment-wide snowpack outflow contribution during individual ROS events. Instead, we performed IHS for each ROS event and individually for each sampling site using the site-specific measurements of snowpack outflow volume and isotopic composition. Thus, we obtained the relative contributions of snowpack outflow to streamflow for three different scenarios during winter 2017, under the assumption that the catchment-wide average snowpack outflow is represented either by the measurements from the mid-elevation grassland (MG), the mid-elevation forest (MF) or the high-elevation grassland (HG) site. By comparing the IHS results of the three scenarios for each ROS event, we seek to quantify the effects of spatial variability in snowpack outflow generation due to vegetation and elevation. As no snowpack outflow could be measured at the HG and MG sites during winter 2018, a three-scenario comparison was not possible for that period.

We also quantified the relative contributions of rainwater (subscript “R”) and pre-event water to streamflow as

$$\begin{array}{}\text{(5)}& F_{\mathrm{R}}={\displaystyle \frac{V_{\mathrm{R}}}{V_{\mathrm{S}}}}={\displaystyle \frac{C_{\mathrm{S}}-C_{\mathrm{pe}}}{C_{\mathrm{R}}^{*}-C_{\mathrm{pe}}}},\end{array}$$

where V_R and $C_{\mathrm{R}}^{*}$ denote the rainfall volume and the volume-weighted isotopic composition in rainwater, respectively. The standard error of the volume-weighted rainwater isotopic composition SE($C_{\mathrm{R}}^{*}$), was estimated analogously to Eq. (5) of von Freyberg et al. (2017):

$$\begin{array}{}\text{(6)}& \mathrm{SE}\left(C_{\mathrm{R},i}^{*}\right)=\sqrt{{\displaystyle \frac{\sum _{k=j}^i{V}_{\mathrm{R},k}{\left(C_{\mathrm{R},k}-C_{\mathrm{R},k}^{*}\right)}^{\mathrm{2}}}{(j-i)\sum _{k=j}^i{V}_{\mathrm{R},k}}}}.\end{array}$$

To quantify the standard error of F_R (SE${}_{F_{\mathrm{R}}}$) we used Eq. (4), in which we replaced SE($C_{\mathrm{spo}}^{*}$) with SE($C_{\mathrm{R}}^{*}$) and $C_{\mathrm{spo}}^{*}$ with $C_{\mathrm{R}}^{*}$. We determined F_R based on the measurements collected at the meteorological station of the Erlenbach catchment. Because snowpack outflow volumes and isotopic compositions were also measured at the same location (i.e. the MG site), we can compare F_spo with F_R to study the role of snowpack storage in streamflow generation. We expect that mixing processes and storage of incoming rainfall in the snowpack result in a more damped isotope signal of snowpack outflow compared with the isotope signal of incoming rainfall. In this case, with all else being equal, the F_spo fraction will be larger than the F_R fraction.

3 Results and discussion

3.1 Variable snow conditions at the three snowmelt lysimeter sites and response of discharge

3.1.1 Spatial and temporal variability of snow properties due to elevation and vegetation

During winter 2017 (1 January–7 May 2017), most of the Erlenbach catchment was covered with a seasonal snowpack. At the beginning of January 2017, when snowfall occurred over several consecutive days during cold conditions (a snow depth of 22 cm was reached within 6 d, and the mean air temperature was −6.6 ^∘C), the seasonal snowpack established simultaneously at all three snowmelt lysimeter sites (Fig. 2a–c). Figure 2 shows that the snow depths at the three snowmelt lysimeter sites differed from one another and varied considerably over time.

Figure 2 Daily precipitation volumes (snow and rainfall) measured at the meteorological station (MG, right axis, a–c), along with hourly air temperature, snow depth, snow water equivalent (SWE) and snowpack outflow volumes measured at the high-elevation grassland site (HG, a), the mid-elevation grassland site (MG, b) and the mid-elevation forest site (MF, c) for the study period from 1 January to 22 May 2017. Panel (d) shows daily discharge at the Erlenbach catchment outlet (on log scale). Vertical grey bars indicate the six rain-on-snow (ROS) events that were analysed in this study. Asterisks (^*) indicate data gaps. Winter 2018 data are provided in the Supplement.

At the MG site, the snow depth was highest on 17 January 2017 (82.2 cm) and the SWE was greatest on 20 February 2017 (168 mm; Fig. 2b). The seasonal snow cover was established on 3 January, became discontinuous on 17 March 2017 and melted completely by 20 March 2017. Two short-term snowpacks established themselves during additional snowfall events in mid and late April 2017. The snow depth measurements at hourly and daily resolution agreed well for most of the study period except for the last 3 weeks of the seasonal snowpack (1–24 March 2017). For this period, the daily snow depth readings from three measurement stakes indicated lower mean snow depths than the readings of the snow depth sensor (hourly data). These measurement differences can be explained by small-scale spatial heterogeneities of the seasonal snowpack caused by wind drift or enhanced melt around the measurement stakes.

Compared with the MG site, the maximum snow depth at the 200 m higher HG site was reached about 7 weeks later and was 55 cm greater (137 cm; Fig. 2a). The maximum SWE (303.4 mm) occurred about 9 d after the peak in snow depth, and was almost twice the maximum SWE at the MG site. According to Stähli et al. (2000) and Stähli and Gustafsson (2006), snow depths, and thus SWE, are generally larger at higher elevations in the Alp catchment due to lower temperatures and thus a greater tendency for winter precipitation to fall as snow. Due to the greater snow depth at the HG site, the seasonal snowpack lasted around 21 d longer than at the MG site. Similar to the seasonal snowpack, the two short-term snowpacks in mid and late April 2017 reached greater snow depths and SWEs compared with the MG site.

At the forested MF site, the snowpack was generally much thinner compared with the nearby grassland MG site, and thus melted out several times during our study period (Fig. 2c). The maximum snowpack depth at the MF site was around 30 cm lower than at the nearby MG site. Based on monthly surveys, the largest SWE (72 mm) occurred on 25 January 2017, roughly 1 month earlier than the maximum at the MG site. Because of the generally shallower snow depths at the MF site, its seasonal snowpack became discontinuous on 21 February 2017, which was 24 d earlier than at the MG site. Four snowfall events (24 February, 6 March, 18 April and 27 April 2017) built up shallow snowpacks at the MF site that lasted only several days (Fig. 2c). An earlier study in the Alp catchment observed that roughly twice as much snow accumulated at grassland sites than at nearby forested sites (Stähli et al., 2000). Generally, snow accumulation under forest is significantly lower due to interception and canopy effects on radiation (i.e. lower short-wave and higher long-wave radiation (Berris and Harr, 1987; Bründl, 1997; Gustafson et al., 2010; López-Moreno and Stähli, 2007; Molotch et al., 2011; Montesi et al., 2004).

Total volumes of snowpack outflow, cumulated over the entire study period, were largest at the HG site (1319±214 mm) and smallest at the MF site (685±78 mm). At the MG site, cumulative volumes of snowpack outflow (816±128 mm) were similar to cumulative volumes of incoming rainfall (833±17 mm) and discharge at the catchment outlet (786 mm).

Weekly snow surveys at the MG site showed that the shallow soil was frozen between 28 December 2016 and 12 March 2017, likely because air temperatures were mostly below 0 ^∘C before the seasonal snowpack was established and the seasonal snow cover prevented the soil frost from thawing despite warmer conditions until mid-March (Goodrich, 1982). At the MF site, soil frost was monitored monthly and the only survey that indicated soil frost was on 28 December 2016, i.e. before the seasonal snowpack was established. On 25 January 2017, the shallow soil at the MF site was no longer frozen.

3.1.2 The spatially variable response of snowpack outflow to rain-on-snow (ROS) events

The 2017 study period was characterised by frequent ROS events which altered the snowpack properties at the three snowmelt lysimeter sites. Six ROS events are discussed in detail below to compare the snowmelt processes at the HG, MG and MF sites. Four more ROS events occurred during winter 2018, but only the MF site provided snowpack outflow volume measurements during that winter, so a site-to-site comparison was not possible.

Figure 2 and Table 1 provide overviews of six ROS events during winter 2017 (events 1–6) and four events during winter 2018 (only MF: events 7–10). During ROS event 1, the tipping bucket rain gauge at the meteorological station stopped working after 13 January 2017 03:40 UTC + 1 h, which coincided with an air temperature decrease to values below 0 ^∘C (Fig. 2). At that point, 21.6 mm of rainfall had been recorded since the beginning of the ROS event on 12 January 2017 17:40 UTC + 1 h. Our webcam images and snow depth data show, however, that the snow depth increased after 03:40 UTC + 1 h, indicating the transition from rainfall to snowfall. Also, river discharge peaked soon afterwards at 05:40 UTC + 1 h, providing further evidence of a catchment-wide transition from rain to snow. Thus, despite the malfunctioning of the tipping bucket rain gauge, we considered the rainfall measurements until 03:40 UTC + 1 h to be representative of the total volume of incoming rainfall during ROS event 1.

Table 1 Start and end times (UTC + 1 h) of 10 rain-on-snow (ROS) events that were identified in the winters of 2017 (all sites) and 2018 (MF site only, as no reliable snowpack outflow volumes were measured at the HG and MG sites in winter 2018). Rainfall characteristics (cumulative rainfall of the ROS event, maximum 4 h rainfall volume, maximum 8 h rainfall volume and rainfall duration) and snowpack outflow volumes at the three snowmelt lysimeter sites are shown. The standard errors (SE) of the snowpack outflow measurements represent the combined effects of measurement uncertainty and spatial heterogeneity of the melt process at each lysimeter site.

^* No snowpack outflow occurred.

Figure 3a compares the cumulative volumes of rainfall and snowpack outflow of the 10 ROS events during winter 2017 (MG, MF and HG sites) and 2018 (MF site only). The snowpack outflow volumes measured at the three snowmelt lysimeter sites were often associated with large uncertainties (Fig. 3a, b), mainly because snowmelt at the plot scale can be very heterogeneous (Kattelmann, 2000; Rücker et al., 2019; Unnikrishna et al., 2002). As each sampling site consisted of three individual lysimeters, we were able to estimate (at least approximately) this spatial variability of snowpack outflow.

Figure 3 A comparison of 10 rain-on-snow (ROS) events during winter 2017 (HG, MG and MF site) and 2018 (events 7–10; MF site only) indicates large spatial and temporal variability of snowpack outflow generation in response to incoming rainfall. (a) Rainfall volumes at the MG site (light blue) and snowpack outflow volumes at the HG (red, black-shaded), MG (yellow) and MF (light green: winter 2017; dark green: winter 2018) sites. Error bars indicate the standard error (SE) of the snowpack outflow measurements (combining measurement uncertainty and spatial heterogeneity of the melt process at each sampling site). (b) Comparison of snowpack outflow volumes and rainfall volumes during the 10 ROS events at the MG and MF sites (colour coding as in a). ROS events with enhanced melt plot above the 1:1 line, whereas events with rainfall retention in the snowpack plot below the 1:1 line. Error bars indicate ± SE. Please note that the MF site was already snow-free during event 6.

At the MG site, the snowpack response to ROS events was highly variable (Fig. 3b), i.e. snowpack outflow was less than incoming rainfall (events 1 and 5), similar to incoming rainfall (events 2, 3 and 4) or more than incoming rainfall (event 6). For events 1 and 5 the differences between rainfall volumes and snowpack outflow were 19 and 42 mm, respectively, which were statistically significant, i.e. larger than 2 times their pooled standard errors.

At the HG site, snowpack outflow volumes were similar to those measured at the 200 m lower MG site (within their pooled standard errors), except for events 4 and 6. During event 4, no snowpack outflow was generated at the HG site, which was probably because the local air temperature was lower and the snowpack was deeper (95 cm) and thus retained more rainwater than the snowpack at the MG site (Fig. 2a). Similarly, less snowpack outflow was recorded at the HG site than at the MG site during ROS event 6, probably because the deeper HG snowpack was not yet saturated. For event 2, however, the measurement differences between the three individual lysimeters at the MG site were particularly large, likely due to lateral flow in the snowpack (Eiriksson et al., 2013; Webb et al., 2018).

The MF and MG sites were located close to each other, so differences in snowpack outflow can be mostly attributed to the effects of forest cover. For ROS events 1 and 5, snowpack outflow at the forested MF site was larger compared with the grassland MG site, whereas it was smaller for event 6. For the remaining events 2–4, the differences between the measured snowpack outflow volumes at the MF and MG sites were only 19.3, 24.1 and 5.4 mm, and were not statistically significant. Note that for events 3 and 4 the MF site had snowpacks of only 5 and 8 cm, respectively (although these were still identified as ROS events because snowpacks at the reference site, MG, were 29 and 43 cm, respectively). Larger snowpack outflow volumes (events 1 and 5) at the MF site can be explained by the shallower snow depths below the forest canopy (Fig. 2a). Hence, the shallower snowpack saturated more rapidly during ROS events and additional meltwater was released, so that more snowpack outflow was generated compared with the MG site where the snowpack was deeper (Berg et al., 1991; Berris and Harr, 1987; Wever et al., 2014). During event 6, the MF site was already snow-free so the lysimeter funnels captured only under-canopy throughfall, which was less than the rainfall volumes measured near the MG site due to interception losses (DeWalle and Rango, 2008; Saxena, 1986). For ROS events 7–9, the snowpack outflow volumes of the MF site were larger than the incoming rainfall volumes, indicating enhanced melt. The highest snowpack outflow volumes during winter 2018 were registered during event 9, which followed 1 d after ROS event 8. The snowpack at the beginning of event 9 was thinner (event 8: 23 cm; event 9: 10 cm) and probably more saturated, with a higher snow bulk density compared with event 8.

The detailed analysis of the six ROS events at the three lysimeter sites during winter 2017 shows that incoming rainfall was attenuated differently in the snowpacks (both among sites and among events), illustrating the challenge of adequately estimating snowpack outflow volumes during ROS events at the plot and catchment scale. Previous studies used rainfall characteristics and snowpack properties to predict the effects of ROS events on catchment outflow (DeWalle and Rango, 2008; Kattelmann, 1997; Würzer et al., 2016). Thus, in the following section, we analysed the processes and properties that control the outflow response of the snowpack.

3.1.3 The effects of snowpack properties and rainfall characteristics on snowpack outflow generation during ROS events

Figure 4a–c compares the snowpack water budgets of the ROS events to the initial snow properties (i.e. bulk snow density, SWE and snow depth) and rainfall characteristics during the event period (maximum cumulative 4 h rainfall, maximum cumulative 8 h rainfall, event duration and mean air temperature) to better understand their effects on snowpack outflow generation. The snowpack water budget was calculated as the volumetric difference between snowpack outflow and incoming rainfall; therefore, positive values of the snowpack water budget indicate enhanced snowmelt, whereas negative values indicate retention of incoming rainfall in the snowpack (Table 1). Note that the MF site was already snow-free prior to ROS event 6; thus, we excluded this data point from the following analysis. At the MG site, ROS event 6 had the most positive snowpack water budget (42.0±26.1 mm), i.e. snowpack outflow was 1.6 times larger than incoming rainfall. The most positive snowpack water budgets at the MF site occurred during events 5 (winter 2017) and 9 (winter 2018).

Figure 4 Correlations between snowpack water budgets (snowpack outflow minus rainfall volume) and initial snow conditions and rainfall characteristics (measured at the MG site) of the rain-on-snow events show the strongest relationship with initial snow depth. Positive values of the snowpack water budget indicate enhanced snowmelt, whereas negative values indicate retention of incoming rainfall in the snowpack. Error bars indicate the uncertainty of the snowpack outflow, i.e. the combined effects of measurement uncertainty and spatial variability. The different scatterplots compare the snowpack water budgets at the MG site during winter 2017 with (a) initial snow density, (b) initial snow water equivalent (SWE), (c) initial snow depth, (d) rainfall intensity presented as maximum 4 and 8 h rainfall volumes, (e) rain duration and (f) mean air temperature (^∘C) during the ROS event. In (c) and (f), data from the MF site are included (light green: winter 2017; dark green: winter 2018), and the sizes of the points indicate the rainfall volume (mm) on an event basis.

Figure 4 shows that most of the relationships between the snowpack water budgets and initial snow properties or rainfall characteristics are highly scattered, indicating a large variability in snowpack outflow generation due to lateral flow and preferential flow pathways in the snowpack. We estimated a positive correlation between the snowpack water budget and initial snow depth (Fig. 4c; R²=0.50; slope = 0.04; p value = 0.02), whereas the other parameters (Fig. 4; bulk snow density, snow water equivalent, rainfall intensity and rainfall duration) were only weakly correlated with the snowpack water budget (R²<0.2; slope < 0.8; p value < 0.9). When additional data for winter 2018 from the MF site were considered, the linear relationship of the snowpack water budget and the initial snowpack became stronger (R²=0.56; slope = −1.7; p value = 0.003).

No consistent relationships emerged between the snowpack water budget and the initial SWE (Fig. 4c) or rainfall duration (Fig. 4e), suggesting that these are poor predictors for snowpack responses to ROS events. This was further confirmed by a multiple linear regression analysis (JMP 14 software, 100 SAS Campus Drive, Cary, NC 27513, USA) that estimated the effects of initial snowpack properties and rainfall conditions on snowpack outflow volumes for the 15 ROS events measured at the MG and MF sites. The best model fit (n=15, RMSE = 22.097; R²=0.8; R² adjusted = 0.76) was obtained with two predictor variables: initial snowpack depth prior to the event (p value < 0.00004) and rain volume (p value = 0.004). The other variables (rainfall intensity, rainfall duration and mean air temperature during event) were not significant predictors of snowpack outflow volumes (p value > 0.2).

A better understanding of snowmelt processes during ROS events could be obtained when individual events or event pairs were analysed in greater detail. Snowmelt at the MG site was most enhanced during event 6, when 66.9 mm of rainfall resulted in 63±39 % (42.0±26.1 mm) more snowpack outflow (total outflow 108.9±26.1 mm). This ROS event occurred during the melt-out phase of the seasonal snowpack at the MG site, when it was isothermal, ripe and already melting (high density 0.402 g cm⁻³, small snow depth 16.9 cm, low SWE 67.9 mm; Table 2). At the end of the ROS event, the snowpack was entirely melted. In addition, rainfall intensities during event 6 were the highest of all six events (maximum 4 h rainfall: 29.1 mm, maximum 8 h rainfall: 36 mm). As a result, the time lag of snowpack outflow relative to incoming rainfall was short (Fig. 5) and incoming rainfall accelerated the melt process (Berg et al., 1991; Colbeck, 1977; MacDonald and Hoffman, 1995; Marks et al., 1998; Wever et al., 2014).

Figure 5 Hourly measurements of precipitation (snowfall and rainfall) and air temperature (a), as well as snow depth and snowpack outflow measured at the HG site (b), at the MG site (c) and at the MF site (d) during the period from 1 January to 22 May 2017. Discharge at the Erlenbach catchment is also shown in (b)–(d). Vertical grey bars indicate six rain-on-snow (ROS) events that are analysed in this study. Due to a data gap in the snow depth sensor – asterisk (^*) in c, daily snow depth values of the webcam are shown for that period.

Table 2 Event characteristics and snowpack water budgets of the 10 rain-on-snow (ROS) events. Snowpack water budgets were calculated for the MG and the MF site by subtracting the snowpack outflow volume from the rainfall volume. Standard errors (SE) of the snowpack water budget were estimated from the measurement uncertainty and the spatial variability of snowpack outflow measurements at the sampling site. Manual snow surveys provided initial snow bulk density and initial snow water equivalent (SWE). Snow surveys were generally performed once or twice a week at the MG site but only monthly at the MF site during the winter of 2017, providing insufficient information at that site (–). Weekly snow surveys at the MF site in winter 2018 did not properly represent the snow properties at the MF site (–) because the snow depth under the forest canopy was much shallower and highly variable over time compared with the MG site; thus, the measured snow bulk densities at the MF site were likely not constant over several consecutive days.

The 90.2 mm ROS event 5 resulted in the most negative snowpack water budget ($-\mathrm{68.4}\pm \mathrm{13.1}$ mm) at the MG site, with 76±14 % of incoming rainfall being retained in the snowpack. By comparing event 5 to event 2, during which a similar amount of rain fell (99.2 mm) and the average snowpack water budget was less negative ($-\mathrm{31.4}\pm \mathrm{44.5}$ mm or 32±45 %), we found that the mean air temperatures during these events were similar (4.43 and 3.14 ^∘C, respectively; Table 2). However, the main difference between the events was the higher hourly mean air temperature prior to event 2 (5.0 ^∘C) compared with event 5 (1.6  ^∘C). This suggests that the cold content of the snowpack prior to event 5 could have been higher – meaning that more rainfall was retained as ice in the snowpack compared with event 2 (DeWalle and Rango, 2008; Maclean et al., 1995; Marks et al., 1998; Wever et al., 2014). It is also possible that the porosity of the snowpack prior to event 5 was higher than prior to event 2 because of a short snowfall event on the morning of 8 March 2017 (Fig. 2). A higher snowpack porosity would increase the potential snowpack storage capacity (DeWalle and Rango, 2008). However, we only have snow density measurements from 2 d prior to event 5; thus, we can only speculate whether snowpack porosity was larger during event 5 than event 2.

A similar analysis could be carried out for ROS events 1 and 3, during which 21.6 and 20.0 mm of rain fell, respectively. During event 1, snowpack outflow volumes at the MG site were 88±6 % (19.1±1.4 mm) less than rainfall, indicating the significant retention of rainwater in the snowpack. The air temperature was low before this event (0.4 ^∘C), resulting in a less dense snowpack with a high cold content, which could retain more rainfall by freezing (DeWalle and Rango, 2008). During event 3, however, the snowpack water budget was very small (12±36 % or 2.45±7.2 mm), suggesting either that most of the rainwater percolated through the snowpack, or that incoming rainwater was stored in the snowpack and replaced an equal volume of meltwater, or a combination of both. The initial conditions of event 3 differed from those of event 1, with a rainfall duration that was twice as long (22 h), a 3.7 ^∘C higher initial air temperature (e.g. mean hourly air temperature prior to the event) and a nearly 2 times denser snowpack (0.34 g cm⁻³; Fig. 4). Thus, although the initial snowpack of event 3 was slightly deeper (39.4 cm) than that of event 1 (26.2 cm), the longer rainfall duration, higher air temperatures and denser snowpack reduced the storage capacity in the snowpack and less rainfall was retained. The significant role of rainfall duration in snowpack outflow generation during ROS events has previously been shown by Kattelmann (1997) and Kroczynski (2004).

The analysis of the 10 ROS events measured at the MG and MF sites illustrates that the generation of snowpack outflow did not entirely depend on the incoming rainfall volume, but also on the initial snowpack conditions that controlled the retention of rainwater and melt processes. This was further illustrated by the hourly measurements of snowpack outflow that indicated highly variable responses and lag times (e.g. the time between the beginning of the ROS event and the first response of the snowpack outflow; the first response is defined as an increase of snowpack outflow by at least 0.05 mm relative to the previous measurement) across the lysimeter sites (Fig. 5).

For three ROS events (1, 2 and 5), the lag times of the MF site were the lowest compared with the HG and MG sites, possibly due to the generally shallower snowpack at the forested site. The snowpacks at the MG and HG sites were deeper, and for most of the events they likely had a larger buffer capacity for incoming rainwater than the shallower snowpack at the MF site. The longest time lags were observed at the MG site during event 2 (e.g. 26 h) when the ROS event occurred with relatively low rainfall volumes (21.6 mm per event) on a fresh snowpack with low density (Tables 1, 2). ROS event 6 occurred when the snowpacks at the HG and MG sites were already ripe (e.g. more than 3 mm d⁻¹ during the previous day), so that snowpack outflow increased immediately without lag times (Fig. 1; MF site was already snow-free). Our observations indicate that the magnitude and timing of snowmelt at the catchment scale strongly depended on the snow properties (such as density and depth) and the degree of ripeness of the snowpack (Berg et al., 1991; DeWalle and Rango, 2008; MacDonald and Hoffman, 1995; Maclean et al., 1995; Wever et al., 2014).

Measurements from the MG, MF and HG sites revealed that snowpack outflow generation was highly variable across space and time and as a result, the contribution of snowpack outflow to river streamflow was very heterogeneous across the catchment landscape. For instance, the streamflow response to ROS event 2 was particularly large, probably because of large snowmelt inputs from higher elevations (HG site; Fig. 5b). Daily pulses of snowmelt from the HG site in late March were also reflected in distinct diurnal variations in stream discharge, suggesting the input of snowmelt mainly from high elevations (Fig. 5b). In contrast, during events 1 and 4, no snowmelt was generated at the HG site, so that the observed discharge peak was likely caused by snowmelt from low and mid elevations (MG and MF sites; Fig. 5c,  d). However, the synchrony of responses does not allow for any conclusions to be drawn about the water sources of streamflow (McDonnell and Beven, 2014).

3.2 Contribution of snowpack outflow and rainwater to catchment outflow

3.2.1 Isotopic composition of rainwater and snowpack outflow

Figure 6 compares the isotopic composition of water samples collected with the three lysimeter systems during the 2017 winter period. Because the lysimeter funnels were permanently installed in the field, they collected snowmelt and rain-on-snow during snow-covered periods, as well as rainfall during snow-free conditions. Thus, we classified the samples as either rainwater (no snowpack), rain-on-snow or snowmelt to better quantify the effects of elevation and vegetation cover on the isotopic signatures of the different water sources. Because of the more persistent snow cover at the HG site, rainfall only occurred as rain-on-snow during the study period, so the HG lysimeter system collected predominantly snowmelt or a mixture of rain and snowmelt. Additionally, the isotopic composition of bulk snow samples at the HG and MG sites are shown (no regular bulk snow sampling was carried out at the MF site). We evaluated the isotopic differences between the water sources at each site and between the sites for the same source with an unpaired two-sample t test. At the MG site, we found that rainwater (no snowpack) was significantly different (i.e. p value < 0.01) from snowmelt, rain-on-snow and bulk snow on average. For the MF and HG sites, the average isotopic compositions of all sources were statistically similar.

Figure 6 Isotopic composition of the samples collected with the three snowmelt lysimeter systems (rainwater during snow-free conditions, snowpack outflow during rain-on-snow events and snowmelt) and bulk snow at (a, d) the high-elevation grassland (HG) site, (b, e) the mid-elevation grassland (MG) site and (c, f) the mid-elevation forest site for the study period from 1 January to 5 May 2017. Except for bulk snow, all samples were collected with the same lysimeter systems; rainwater was collected during snow-free conditions, whereas rain-on-snow and snowmelt were collected during conditions with snow cover. The isotopic composition of bulk snow is added for the HG and MG sites (a and b; grey). Panels (a)–(c) show boxplots of the isotope values δ²H (left) and δ¹⁸O (right). Panels (d)–(f) show the isotope values δ²H and δ¹⁸O plotted in dual isotope space along with the local meteoric water lines (dashed lines and equations, derived from rainwater samples collected at each field site between May and October in 2017. The slopes of the three regression lines were not statistically different, i.e. p values > 0.01).

We also evaluated the isotopic differences between the sites for the same source (unpaired two-sample t test), however, none of them were statistically significant. Nevertheless, some site-to-site isotopic differences could be observed that permit a more detailed analysis of the isotope effects due to forest cover and elevation.

Rainwater at the MG site and throughfall at the MF site had similar isotopic compositions (differences in median δ²H and δ¹⁸O were 0.4 ‰ and 0.0 ‰, respectively), however, interception and mixing of rainwater in the forest canopy resulted in a wider range of isotope values in throughfall at the MF site compared to the grassland site (Fig. 6b, c). Furthermore, our data show that rain-on-snow and snowmelt under the forest canopy (MF site) were isotopically slightly heavier than the corresponding samples from the nearby grassland (MG) site. The absolute differences in the median δ²H values between the two sites were 8.6 ‰ and 5.3 ‰ for rain-on-snow and snowmelt, respectively (the corresponding differences in median δ¹⁸O were 1.5 ‰ and 0.5 ‰ for rain-on-snow and snowmelt, respectively). This isotopic difference suggests that canopy-intercepted snow at the forest site underwent enhanced isotopic fractionation such that throughfall (and thus the snowpack) became isotopically heavier under forest cover compared with open grassland (Claassen and Downey, 1995; Koeniger et al., 2008).

At the grassland HG site, the median isotopic composition of bulk snowpack and snowpack outflow (rain-on-snow and melt) was heavier than at the lower grassland (MG) site; the median δ²H values differed by 8.4 ‰ in rain-on-snow, by 2.4 ‰ in snowmelt and by 3.2 ‰ in bulk snow (Fig. 6a, d), but the isotopic differences in median δ¹⁸O were smaller than 0.4 ‰. The isotopically heavier bulk snow and snowpack outflow at the higher site is the opposite of the expected altitude effect (Dietermann and Weiler, 2013; Moser and Stichler, 1970), but one must remember that the snowpacks at the three different sites lasted for different spans of time, during which they received different rain and snow inputs with different isotopic compositions. Additionally, the elevation range captured by our snowmelt lysimeter sites was only 220 m; therefore, one should not expect to see a conventional altitude effect (which in any case would be small) in field data like ours.

3.2.2 Temporal and spatial isotopic variation of bulk snow and snowpack outflow during rain-on-snow events and snow melt

Due to frequent melt periods and ROS events, values of δ²H in bulk snow and snowpack outflow were highly variable over time at all three lysimeter sites (Fig. 7). Similar to other studies (Gustafson et al., 2010; Taylor et al., 2002a), our data show that snowpack outflow generation can be much more variable in time than would be implied by weekly bulk snow samples alone. Nonetheless, bulk snow samples at the HG and MG sites mirrored the general isotopic pattern of the snowpack outflow samples. For instance, at the HG site, both sample types indicate a clear isotopic enrichment during melt-out of the seasonal snowpack in early April 2017 (Fig. 7a). However, bulk snow samples, which were collected only weekly or twice a week, could not capture the high temporal variability that was observed in snowpack outflow (e.g. during ROS event 6 at both HG and MG sites, and during event 5 at the MG site; Fig. 7).

Figure 7 Deuterium (δ²H) values in precipitation (snowfall and rainwater; light blue) and snowpack outflow (separated into rain, rain-on-snow and snowmelt) indicate spatial and temporal variability across (a) the high-elevation grassland (HG; red) site, (b) the mid-elevation grassland (MG; yellow) site and the mid-elevation forest (MF; green) site (d), and in stream water (grey) at the Erlenbach outlet during the study period from 1 January to 22 May 2017. Stream water isotopic composition (grey) is indicated in (a)–(c) for reference, using grey dashed lines that represent the range between −85 ‰ and −75 ‰. Error bars indicate the standard error of the isotopic composition of snowpack outflow due to spatial heterogeneity at the plot scale.

Our daily isotope measurements of rainwater and snowpack outflow across the catchment landscape allowed us to study the temporal and spatial isotopic variation of snowpack outflow during rain-on-snow events and snow melt. Our main observations are as follows:

1. The isotopic composition of precipitation affects that of snowpack outflow: during most ROS events, the isotopic composition of snowpack outflow mirrored that of incoming rainfall (events 1, 3, 5 and 6; Fig. 7; no rainfall data available for event 2). For instance, during event 6, δ²H in incoming rainwater was −143.5 ‰ and δ²H in snowpack outflow (rain-on-snow) changed from −80.0 ‰ to −137.3 ‰ and from −88.8 ‰ to −110.5 ‰ at the MG and HG sites, respectively.

2. Snow depth controls isotopic response to ROS events: similar to the snowpack outflow volumes (Sect. 3.1.3), the isotopic response of the snowpack to individual ROS events likely depended on the local initial snowpack properties and the event magnitude. Isotopic responses in snowpack outflow were more damped at the HG site compared with the signals measured at the MG and MF sites, because the snowpack was deeper at the higher elevation site. A similar effect of snow depth was also apparent at all other sites: the isotopic variability of snowpack outflow was smaller when the seasonal snowpack was relatively deep (e.g. between events 2 and 5 at the MG site), and the variability increased when the snowpack became shallower, including during the two short-term snowpacks (e.g. between 17 April and 4 May 2017 at the MG site). At the MG site, rainwater and snowpack outflow had very similar isotopic compositions during event 6 (i.e. no damping), because the ripe shallow (17 cm) snowpack enabled the vertical percolation of incoming rainwater (Fig. 7b; Kroczynski, 2004). At the HG site, however, the snowpack was deeper (91 cm) during event 6, and incoming rainwater was mostly retained in the snowpack, resulting in a damped isotopic response in snowpack outflow (Fig. 7a).

   The isotopic signal of incoming rainwater can be altered as it percolates through the snowpack, depending on snow metamorphism and isotopic exchange (Judy et al., 1970). A significant isotopic depletion or enrichment of snowpack outflow due to such rain-on-snow events has already been reported in other studies (Herrmann, 1978; Juras et al., 2016; Shanley et al., 1995a; Unnikrishna et al., 2002). The isotopic exchange in the snowpack is mainly controlled by the residence time of liquid water (snowmelt and rain-on-snow) in the snowpack, which, in turn, is determined by the depth and the density of the snowpack (Taylor et al., 2001, 2002b), the rainfall magnitude (Herrmann et al., 1981) and the flow rate of percolating liquid water. As a result, deeper snowpacks generally cause slower rainwater throughflow, which enhances isotopic redistribution in the snowpack and isotopic exchange between the liquid water and solid ice (Lee et al., 2010a, b; Taylor et al., 2001).

3. Light isotopes preferentially leave the snowpack during melt: isotopic variations in snowpack outflow can also result from freeze–melt processes in the snowpack during rain-free periods. For instance, snowpack outflow at the MG site became isotopically lighter than bulk snow after ROS event 3, despite rainwater being isotopically heavier than bulk snow during ROS event 3 (Fig. 7b).

   The isotopic contrast between the snowpack outflow of the last day of event 3 and the following day was 13.4 ‰ for δ²H and 2.2 ‰ for δ¹⁸O. This depletion signal occurred simultaneously with a decrease in air temperature to below 0 ^∘C, suggesting isotopic fractionation effects in the snowpack due to partial phase transitions of liquid water to ice (Herrmann et al., 1981; Shanley et al., 1995a; Stichler et al., 1981; Taylor et al., 2001). During partial freezing, the liquid phase becomes isotopically lighter, because the heavier isotopes preferentially transition into the solid phase, i.e. the lower free energy state (Hoefs, 2018). During this partial freezing process, lighter isotopes preferentially leave the snowpack, and over cycles of melting and refreezing, the snowpack becomes isotopically heavier and more homogeneous (Huth et al., 2004; Judy et al., 1970; Lee et al., 2010b; Schmieder et al., 2016; Taylor et al., 2002b, 2001; Unnikrishna et al., 2002). During the melt-out period of the seasonal snowpack, this fractionation effect results in snowpacks and snowpack outflows that become isotopically heavier over time. This trend could be observed at the HG site, where rising air temperatures and dry conditions between 26 March and 9 April 2017 resulted in progressive melt of the seasonal snowpack (Fig. 7a). This melt-out of the seasonal snowpack was accompanied by a gradual isotopic enrichment in snowpack outflow δ²H from −96.5 ‰ to −84.0 ‰. The isotopic composition of streamflow mirrored this isotopic trend in snowpack outflow, suggesting that snowmelt from higher elevations contributed to catchment outflow during the melt-out period.

3.2.3 The contribution of rainwater and snowpack outflow to river discharge

1. Implications of isotopic variability in snowpack outflow for end-member mixing analysis: Fig. 7 shows that the isotopic signal of Erlenbach stream water was affected by incoming rainwater and snowpack outflow during the individual ROS events. In the following section we quantified the contribution of rainwater and snowpack outflow to streamflow, using stable water isotopes as conservative tracers in two-component hydrograph separations. These analyses were carried out individually for each sampling site using the volumes and isotopic compositions of their snowpack outflows. Thus, our results reflect the relative snowpack outflow contribution to streamflow for three different scenarios that assume that the catchment-average snowpack is represented by the mid-elevation grassland (MG), the mid-elevation forest (MF) or the high-elevation grassland (HG) sites, respectively. For comparison, we also performed hydrograph separation using rainwater as the “new water” end-member. In all cases, pre-event stream water isotopic composition was used as the “old water” end-member, following conventional practice in two-component hydrograph separations. Here we present our results based on δ²H, for which the temporal variations in stream water were larger, and the measurement uncertainties were smaller, compared with δ¹⁸O.

2. IHS results are highly variable across sites and rain-on-snow (ROS) events: Fig. 8 summarises the estimated contributions of rainwater and snowpack outflow to peak streamflow during the six ROS events. Snowpack outflow contributions to peak streamflow varied among the six ROS events and the three snowmelt lysimeter sites, ranging from 34±7 % to 42±2 % at the HG site, from 13±1 % to 58±3 % at the MG site and from 7±4 % to 91±20 % at the MF site (Fig. 8, Table 3). The maximum fractions of snowpack outflow during each event were different than the contribution to peak streamflow (Table S2 in the Supplement). Maximum fractions of snowpack outflow to streamflow during events most commonly occurred the day after the peak streamflow, and were often significantly higher than during peak streamflow (for example, using δ²H from the HG site as the snowmelt source during event 6, estimated snowmelt contributions to streamflow were 34±7 % during peak streamflow but 75±18 % the day afterward). This indicates that pre-event water dominates the streamflow during peak flow, whereas snowpack outflow contributions dominate the streamflow during the recession limb (von Freyberg et al., 2018b).

   The different results among the three sampling locations reflect the highly variable isotopic compositions of the snowpack outflow across the catchment. For example, during event 6, isotope data from the HG site suggested a significant snowpack outflow contribution to discharge (34±7 %), whereas isotope measurements from the lower-elevation MG site implied a much smaller (13±1 %) contribution of snowpack outflow to streamflow. The different estimates at the two sites can be explained by the stronger retention of incoming rainwater in the higher-elevation (HG) snowpack (resulting in snowpack outflow that was isotopically closer to streamflow) and the transmission of rainwater through the snowpack at the MG site (resulting in snowpack outflow that resembled the isotopically light incoming precipitation; see Sect. 3.1.3, Fig. 8 and Table 3).

   This comparison raises an important point of interpretation. Any two-component hydrograph separation is based on the fundamental assumption that there are only two end-members (in our case, one of the snowpack outflows, and “old water” represented by pre-event streamflow). Thus one cannot interpret the results above as demonstrating that more snowpack outflow reached the stream from high-elevation sites like HG than from lower-elevation sites like MG. Instead, what these results show is that if the catchment-wide snowpack outflow resembled that of the MG site, it could only make a small contribution to streamflow (because otherwise the peak streamflow would need to be isotopically lighter than it in fact was), but if the catchment-wide snowpack outflow resembled that from the high-elevation HG site, it could plausibly make a larger contribution to streamflow.

3. Larger contribution of snowpack outflow to streamflow compared with rainwater: the contributions of rainwater to streamflow during four events (1, 3, 5 and 6) ranged between 5±2 % and 34±2 % (no estimate could be obtained for event 2 due to a gap in rainwater isotope sampling, and unrealistic hydrograph separation results were obtained for event 4 due to overlapping isotope values of stream water and rainwater). Based on snowmelt lysimeter data from the MG site (unrealistic hydrograph separation results were obtained for event 3), the contributions of snowpack outflow to streamflow were larger than those of incoming rainwater during three events (34±2 % vs. 58±3 % for event 1, 25±1 % vs. 50±5 % for event 5 and 12±1 % vs. 13±1 % for event 6; Table 3). The isotopic composition of snowpack outflow at the MG site was often more damped than that of rainwater during most events due to mixing and fractionation processes in the snowpack (Sect. 3.2.2; Fig. 8). As a consequence, the snowpack outflow was isotopically more similar to that of streamflow, which resulted in larger F_spo fractions compared with F_R.

   Although the number of ROS events in our data set is small, our results are in line with previous studies showing that the differences between hydrograph separation results obtained for rainwater and snowpack outflow can potentially be large and should be considered in snow-dominated catchments (Buttle et al., 1995). Our analysis assumes that the end-members of the different scenarios (i.e. snowpack outflow at the MG, MF or HG site) are representative for the whole Erlenbach catchment. However, the catchment is characterised by a diverse vegetation cover (22 % partially forested, 53 % forested and 25 % grassland) and surface topography (altitude 1000–1500 m a.s.l.), so that the “real” contribution of snowpack outflow to streamflow is likely to lie between the estimates derived from the three scenarios. The hydrograph separation estimates for the three lysimeter sites can only provide a probable range of snowpack outflow contributions to discharge from different landscapes of the catchment. As shown here, the estimated contributions of snowpack outflow to streamflow can vary considerably due to differences in landscape characteristics, rainfall magnitude and snowmelt processes. Future sampling strategies should take this spatial and temporal variability in snowpack outflow into account.

4. End-member mixing analyses at the snowpack scale: our data set can potentially be used for isotope-based end-member mixing analysis at the snowpack scale, to estimate the fractional contribution of rainwater to snowpack outflow. Such an analysis requires isotope measurements in the event-water end-member (rainwater), the pre-event water end-member (pre-event snowpack outflow) and the mixture of both (snowpack outflow). For one of the six ROS events (e.g. event 1), no pre-event end-member could be determined because no snowpack outflow was generated before the onset of the event at all three lysimeter sites. Due to a data gap in rainwater sampling, no analysis could be performed during event 2. For the remaining four events, we estimated the fraction of rainwater in snowpack outflow at all three sites, except for the HG site, where no pre-event snowpack outflow occurred prior to events 3 and 4. The results, summarised in Table S3, show that event water (rainfall) can comprise almost none, or almost all, of snowpack outflow; however, these results are highly uncertain and do not allow for an in-depth analysis.

Table 3 Relative contributions of rainwater or snowpack outflow to daily peak discharge based on two-component isotope hydrograph separation (IHS) using δ²H. The IHS was carried out with four different isotope data sets that were collected with snowmelt lysimeters at the HG (high-elevation grassland) site, the MG (mid-elevation grassland) site, and the MF (mid-elevation forest) site as well as with the rainwater collector at the catchment outlet.

^a No snowpack outflow occurred. ^b Data gap.

Figure 8 Relative contribution of snowpack outflow to streamflow at peak flow based on isotopic hydrograph separation for the six rain-on-snow events from winter 2017, including the incoming rainwater (blue, not filled) and the snowpack outflow of the high-elevation site (red, black-shaded), the mid-elevation grassland site (yellow) and the mid-elevation forest site (green). For some events, no data (^*) were available (no melt) or the results were unrealistic due to overlapping isotopic composition of the snowpack outflow (event 3) and rainwater (event 4) with stream water (${}^{**}$). The error bars indicate the standard error of the snowpack outflow contribution to streamflow (see Sect. 2.4).

4 Summary and conclusions

In many mountain regions, global warming is predicted to lead to more frequent rain-on-snow (ROS) events, which can enhance snowmelt and increase the risk of destructive winter floods. However, the processes leading to such enhanced melt are spatio-temporally heterogeneous, so that model-based predictions of discharge peaks during ROS events can be highly uncertain.

By using three automated snowmelt lysimeter systems, located along an elevation gradient of 1185 to 1420 m a.s.l. in a partly forested pre-Alpine catchment, we were able to capture the spatial and temporal variability of snowpack outflow generated over the winter season (Figs. 2, 5). A comparison of snowpack properties at a grassland and a nearby forested site showed that canopy interception significantly reduced incoming snowfall; thus, the maximum snow depth under forest cover was around 20 cm shallower than that of open grassland. Measurements from two grassland lysimeter sites located at different elevations (1220 and 1420 m a.s.l.) showed that the snowpack was on average 55 cm deeper and snowmelt occurred 21 d later at the higher site.

To better understand how snowpack outflow is generated during ROS events across the catchment landscape, we studied 10 ROS events in greater detail (Fig. 3). The ROS events were defined by rainfall rates greater than 0.1 mm per hour, a total rainfall volume of at least 20 mm within 12 h, air temperatures above 0 ^∘C and an initial snowpack depth of at least 10 cm. We found that the snowpack outflow volumes during ROS varied considerably across the three lysimeter sites, and that this variability was linked to rainfall characteristics and initial snowpack properties (Fig. 4). Initial snow depth and rainfall volume explained most of the event-to-event variability in snowpack outflow volumes. Overall, more rainwater was retained in the snowpack at the grassland sites (Fig. 3), which had deeper snowpacks compared to the forest site. Our data show that long and high-intensity ROS events can result in particularly high discharge peaks, even in mid-winter when the snowpack is not saturated (e.g. ROS event 2). This suggests that enhanced snowmelt during ROS events, and/or high antecedent moisture due to ongoing snowmelt, are not limited to late winter when the snowpack is mature and saturated.

We used daily stable water isotope measurements in snowpack outflow, rainwater and stream water to draw inferences about transport and mixing of rainfall within the snowpack during individual ROS events. Depending on the local rainfall characteristics and the snowpack properties, the isotopic responses in snowpack outflow could be either strongly or weakly damped, indicating large spatio-temporal variations of the snowmelt process (Fig. 7). Consequently, isotope-based two-component hydrograph separation (IHS) for estimating snowpack contributions to streamflow often yielded very different results (Fig. 8), depending on which site-specific snowpack outflow isotopic compositions were used. This range of IHS results provides reasonable estimates of relative snowpack outflow contributions to streamflow during individual ROS events, under the assumption that the three lysimeter sites are representative of the snowmelt processes at the catchment scale. Furthermore, our IHS results vary over a wide range, implying that in steep, partly forested catchments like Erlenbach, estimates of snowpack outflow contributions to streamflow derived from bulk snow samples or outflow samples collected at only one location can be highly uncertain. This is in line with the study by Fischer et al. (2017), which showed strong spatial variability in rainwater isotopic composition in the southern Alptal catchment. Using rainwater isotope data in the IHS analysis suggests that the relative contribution of rainwater to streamflow may often be much smaller than the contribution of snowpack outflow, because snowpack outflow is a mixture of both rainwater and snowmelt. Our analysis suggests that snowpack outflow can contribute substantially to streamflow during ROS events and that these contributions depend strongly on the local snowpack properties and rainfall characteristics.

In order to obtain more realistic estimates of snowpack outflow contributions to streamflow during ROS events, snowpack outflow volumes and their isotopic compositions could be interpolated across the study area using a spatially distributed snowmelt model. Recent snowmelt modelling approaches at the catchment scale do not, however, explicitly simulate snowpack outflow during rain-on-snow events (Ala-aho et al., 2017; Lyon et al., 2010; Smith et al., 2016), or use stable water isotopes to track flow pathways (Kormos et al., 2014; Marks et al., 2001; Rössler et al., 2014; Storck et al., 1998). Thus, our spatio-temporally distributed isotope measurements could be beneficial for testing and improving existing snowmelt models (Zappa et al., 2015).