Modelling liquid water transport in snow under rain-on-snow conditions – considering preferential flow

Rain on snow (ROS) has the potential to generate severe floods. Thus, precisely predicting the effect of an approaching ROS event on runoff formation is very important. Data analyses from past ROS events have shown that a snowpack experiencing ROS can either release runoff immediately or delay it considerably. This delay is a result of refreeze of liquid water and water transport, which in turn is dependent on snow grain properties but also on the presence of structures such as ice layers or capillary barriers. During sprinkling experiments, preferential flow was found to be a process that critically impacted the timing of snowpack runoff. However, current one-dimensional operational snowpack models are not capable of addressing this phenomenon. For this study, the detailed physics-based snowpack model SNOWPACK is extended with a water transport scheme accounting for preferential flow. The implemented Richards equation solver is modified using a dual-domain approach to simulate water transport under preferential flow conditions. To validate the presented approach, we used an extensive dataset of over 100 ROS events from several locations in the European Alps, comprising meteorological and snowpack measurements as well as snow lysimeter runoff data. The model was tested under a variety of initial snowpack conditions, including cold, ripe, stratified and homogeneous snow. Results show that the model accounting for preferential flow demonstrated an improved overall performance, where in particular the onset of snowpack runoff was captured better. While the improvements were ambiguous for experiments on isothermal wet snow, they were pronounced for experiments on cold snowpacks, where field experiments found preferential flow to be especially prevalent.


Introduction
The flooding potential of rain-on-snow (ROS) events has been reported for many severe floods in the US (Kattelmann, 1997;Kroczynski, 2004;Leathers et al., 1998;Marks et al., 2001;McCabe et al., 2007), but also in Europe (Badoux et al., 2013;Freudiger et al., 2014;Rössler et al., 2014;Sui and Koehler, 2001;Wever et al., 2014b) where for example up to 55 % of peak flow events could be attributed to ROS events for some parts of Austria (Merz and Blöschl, 2003).With rising air temperature due to climate change, the frequency of ROS is likely to increase in high-elevation areas (Surfleet and Tullos, 2013) as well as in high latitudes (Ye et al., 2008).Besides spatial heterogeneity of the snowpack and uncertainties in meteorological forcing, deficits in process understanding make the consequences of extreme ROS events very difficult to forecast (Badoux et al., 2013;Rössler et al., 2014).For hydro-meteorological forecasters, it is particularly important to know a priori how much and when snowpack runoff is to be expected.Particularly, a correct temporal representation of snowpack processes is crucial to identify whether the presence of a snowpack will attenuate or amplify the generation of catchment-wide snowpack runoff.Most studies investigating ROS only consider the generation of snowpack runoff on a daily or multi-day timescale, where an exact description of water transport processes is less important than for sub-daily timescales (Wever et al., 2014a).Water transport processes are further usually described for snowmelt conditions, but not for ROS conditions, where high rain intensities may fall onto a cold snowpack below the freezing point.In this study however, we particularly focus on snowpack runoff generation at sub-daily scales with special attention to the timing of snowpack runoff which is influenced by preferential flow (PF).
Many studies have shown that flow fingering or PF is an important water transport mechanism in both laboratory experiments (Hirashima et al., 2014;Katsushima et al., 2013;Waldner et al., 2004) and under natural conditions, using dye tracer (Gerdel, 1954;Marsh and Woo, 1984;Schneebeli, 1995), temperature investigations (Conway and Benedict, 1994) or by measuring the spatial variability of snowpack runoff (Kattelmann, 1989;Marsh andPomeroy, 1993, 1999;Marsh and Woo, 1985).The variability of snowpack runoff is defined by the distribution and size of preferential flow paths (PFPs), which are dependent on the structure of the snowpack and weather conditions (Schneebeli, 1995).Beyond its importance for hydrological implications, PF may also be crucial for wet snow avalanche formation processes, where snow stability can be depending on the exact location of liquid water ponding (Wever et al., 2016a).
Most snow models describe the water flow in snow as a uniform wetting front, thereby implicitly only considering the matrix flow component.The history of quantitative modelling of water transport in snow starts with Colbeck (1972), who first described a gravity drainage water transport model for isothermal, homogeneous snow.This was done by applying the general theory of Darcian flow of two-fluid phases flowing through porous media, neglecting capillary forces.Because water transport is not just occurring in isothermal conditions and snow can therefore not be treated as a classical porous medium, Illangasekare et al. (1990) were the first to introduce a 2-D model being able to describe water transport in subfreezing and layered snow.A detailed multilayer physics-based snow model, where water transport was governed by the gravitational part of the Richards equation (RE) described in Colbeck (1972), was introduced by Jordan (1991).With the implementation of the full RE described by Wever et al. (2014a), the influence of capillary forces on the water flow was firstly represented in an operationally used snowpack model.
A model accounting for liquid water transport through multiple flow paths was developed by Marsh and Woo (1985), but was not able to explicitly account for structures like ice layers and capillary barriers.Recently, multi-dimensional water transport models have been developed, which allow for the explicit simulation of PFPs (Hirashima et al., 2014).These models are valuable for describing spatial heterogeneities and persistence of PFPs, but have not yet been shown to be suitable for hydrological or operational purposes.In general, multi-dimensional models are limited by the fact that they are computationally intensive, thus not thoroughly validated for seasonal snowpacks, and still lack the description of crucial processes such as snow metamorphism and snow settling.
In snowpack models which are used operationally, PFPs are not yet considered.The recently introduced RE solver for SNOWPACK led to a significant improvement of modelled sub-daily snowpack runoff rates.For this paper, we further modified the transport scheme for liquid water by implementing a dual-domain approach to represent PFPs.This new approach is validated against snow lysimeter measurements which were recorded during both natural and artificial ROS events.
This study aims to better describe snowpack runoff processes during ROS events within snowpack models that can be used for operational purposes such as avalanche warning and hydrological forecasting.This requires that the model results remain reliable, i.e. that improvements are not realized at the expense of a decreased model performance during periods without ROS, and that the model must not be too computationally expensive.This is the first study to test a water transport scheme accounting for PF which has been implemented in a snowpack model that meets the above requirements.
Our analysis of simulations of over 100 ROS events targets the following research questions: -Is snowpack runoff during ROS in a 1-D model better reproduced with a dual-domain approach to account for PF than with traditional methods considering matrix flow only?
-Are there certain snowpack or meteorological conditions, for which the performance specifically benefits if PF is represented in the model?
This paper is structured as follows: Sect. 2 describes the snowpack model setup, the water transport models, input data and the event definition.Results of the simulations are shown in Sect.3.This includes data of sprinkling experiments of ROS (3.1), natural ROS events (3.2) and the validation of the model on a long-term dataset from two alpine snow measurement sites (3.3).The results are discussed in Sect.4, followed by the general conclusions found in Sect. 5.

Methods
All results in this study are derived from simulations with the one-dimensional physics-based snowpack model SNOW-PACK (Bartelt and Lehning, 2002;Lehning et al., 2002a, b;Wever et al., 2014a) using three different water transport schemes, described in Sect.2.2.The model was applied to four experimental sites that were set up for this study in the vicinity of Davos (Sect.2.3).These sites were maintained over two winter seasons between 2014 and 2016 where data were recorded for several natural ROS events.At the same sites, we conducted a set of six sprinkling experiments to simulate ROS events for given rain intensities (Sect.2.4).Furthermore, we conducted simulations for two extensive datasets from the European Alps: Weissfluhjoch (Switzerland, 46.83 • N, 9.81 • E, 2536 m a.s.l., WSL Institute for Snow and Avalanche Research SLF (2015), abbreviated as WFJ in the following) and Col de Porte (France, 45.30 • N, 5.77 • E, 1325 m a.s.l., Morin et al. (2012), abbreviated as CDP in the following).These datasets provide meteorological input data for running SNOWPACK as well as validation data, including snowpack runoff.Both datasets have already been used for simulations with SNOWPACK (Wever et al., 2014a) and provide data over more than 10 years each.
Below, the SNOWPACK model and the different water transport models are described first, followed by the description of the field sites for ROS observation in the vicinity of Davos.Then, we detail the setup of the artificial sprinkling experiments.After summarizing the WFJ and CDP dataset, we finally present the definition of ROS events that is used in this study.Most analyses were performed in R 3.3.0(R Development Core Team, 2016) and figures were created with base graphics or ggplot2 (Wickham, 2009).

Snowpack model setup
The setup of the SNOWPACK model is similar to the setup used for simulations in Würzer et al. (2016).For all simulations, snow depth was constrained to observed values, which means that the model interprets an increase in observed snow depth at the stations as snowfall (Lehning et al., 1999;Wever et al., 2015).Because the study focuses on the event-scale and snowpack runoff is essentially dependent on the properties of the available snow, this approach was chosen such that we have the most accurate initial snow depth at the onset of the events to achieve the best comparability between the three water transport models.The temperature used to determine whether precipitation should be considered rain (measurements from rain gauges) or snow (from the snow depth sensors) was set to achieve best results for reproducing measured snow height for precipitation driven simulations for the Davos field sites (between 0 and 1.0 • C).For WFJ and CDP, this threshold temperature was set to 1.2 • C, where mixed precipitation occurred proportionally between 0.7 and 1.7 • C. Turbulent surface heat fluxes are simulated using a Monin-Obukhov bulk formulation with stability correction functions of Stearns and Weidner (1993), as described in Michlmayr et al. (2008).At the Davos field sites (Sect.2.3) incoming longwave radiative flux is simulated using the parameterization from Unsworth and Monteith (1975), coupled with a clear-sky emissivity following Dilley and O'Brien (1998), as described in Schmucki et al. (2014).For the roughness length z 0 , a value of 0.002 m was used for all simulations at the Davos field sites and WFJ, whereas a value of 0.015 was used for CDP.The model was initialized with a soil depth of 1.4, 2.2 and 2.14 m (for WFJ, CDP and Davos field sites, respectively) divided into layers of varying thickness.For soil, typical values for coarse material were chosen to avoid ponding inside the snowpack due to soil saturation.The soil heat flux at the lower boundary is set to a constant value of 0.06 W m −2 , which is an approximation of the geothermal heat flux.

Water transport models
The two previously existing methods for simulating vertical liquid water movement within SNOWPACK are either a simple so-called bucket approach (BA) (Bartelt and Lehning, 2002) or solving the RE, a recently introduced method for SNOWPACK (Wever et al., 2014a, b).
The BA represents liquid water dynamics by an empirically determined irreducible water content θ r which defines whether water stays in the corresponding layer or will be transferred to the layer below.This irreducible water content varies for each layer according to Coléou and Lesaffre (1998).The RE represents the movement of water in unsaturated porous media.Its implementation in SNOW-PACK and a detailed description can be found in Wever et al. (2014a).
The PF model presented in this study is based on the RE model, but follows a dual-domain approach, dividing the pore space of the snowpack into a part representing matrix flow and a part representing PF.For both domains the RE is solved subsequently.The PF model is described by (i) a function for determining the size of the matrix and preferential flow domain, (ii) the initiation of PF (i.e., water movement from matrix flow to PF) and (iii) a return flow condition from PF to matrix flow.
The area of the preferential domain (F ) is as a function of grain size (Eq.1), which has been determined by results of laboratory experiments presented by Katsushima et al. (2013): where r g is grain radius (mm).F is limited between 1 and 90 % for reasons of numerical stability.The matrix domain is then accordingly defined as (1 − F ). Water is transferred from the matrix domain to the preferential domain if the water pressure head for a layer in the matrix domain is higher than the water entry pressure of the layer below, which can, according to Katsushima et al. (2013), also be expressed as a function of grain size.This condition is expected to be met if water is ponding on a microstructural transition (i.e.capillary barriers, ice lenses) inside the snowpack.Additionally, saturation was equalized between the matrix and the preferential domain, in case the saturation of the matrix domain exceeded the one in the preferential domain.To move water back into the matrix part, we apply a threshold in saturation of the PF domain and water will flow back to the matrix domain once this threshold is exceeded.This threshold is used as a tuning parameter in the model.Refreezing of liquid water in the snowpack is crucial for modelling water transport in subfreezing snow and may also be important for modelling PF.The presented PF model has also been used to simulate ice layer formation under the presence of PF by Wever et al. (2016b).Thereby, a sensitivity study on the role of refreeze in the PF domain and the return flow condition from PF to matrix flow was conducted.It was found that neglecting refreeze led to the best results for reproducing ice layer formation, but did not significantly affect the performance in reproducing measured hourly snowpack runoff.Therefore, refreeze in the preferential domain is neglected in the presented study.The threshold in saturation for PF (return flow condition) was also determined by the sensitivity study described in Wever et al. (2016b).While they determined a threshold in saturation of 0.1 to reproduce ice-layer observations at WFJ best, a value of 0.06 was determined to reproduce observed seasonal runoff best.We therefore used the value of 0.06.In contrast to Wever et al. (2016b), we did not set the hydraulic conductivity in soil to 0, because this can lead to an inaccurate representation of observed lysimeter runoff due to modelled ponding on soil, which is not expected to happen on a snow lysimeter.Further details on the implementation of the PF model and its performance can be found in Wever et al. (2016b).
In summary, the PF model accelerates liquid water transport in the preferential domain by concentrating water mass in a smaller area, representing the area fraction of flow fingers in the snowpack.The saturation in the preferential domain is hence higher and unsaturated conductivity is larger.Further acceleration is achieved by disabling refreeze in the preferential domain.

Davos field sites
Four field sites have been installed within an elevational range of 950 to 1850 m a.s.l. in the vicinity of Davos, Switzerland, with one meteorological station and 3-4 snow lysimeters each (15 in total, 0.45 m diameter).The meteorological stations provided most data necessary for running the SNOWPACK model and missing parameters were estimated as described in Sect.2.1.Lysimeters were installed at ground level with an approximate spacing of 10 m horizontal distance.The lysimeters consisted of a funnel attached to a precipitation gauge buried in the ground, which monitored snowpack runoff with a tipping bucket.To block lateral inflow at the snow-soil interface, each lysimeter was equipped with a rim of 5 cm height around the inlet.The multiple snow lysimeter setups allowed analysing the spatial heterogeneity of snowpack runoff.Snowpack properties (SWE, LWC, HS, TS) were manually measured directly before each natural ROS event so that the initial conditions of the snowpack are known in detail.LWC was measured with the "Denoth meter", a device introduced by Denoth (1994).The onset of runoff was defined as the time when cumulative snowpack runoff (measured and simulated, respectively) has reached 1 mm.

Sprinkling experiment description
During winter 2014/15, a total of six artificial sprinkling experiments were performed on all four Davos field sites described above to be able to investigate snowpack runoff generation for different snowpack properties (Table 1).For each experiment, a sprinkling device was placed above a snow lysimeter, covered by an undisturbed natural snowpack, i.e. each lysimeter was only used for one experiment.The device used for sprinkling was a refined version of the portable sprinkling device described in Juras et al. (2013Juras et al. ( , 2016a)).The water used for sprinkling was mixed with the dye tracer Brilliant Blue FCF (concentration 0.4 g L −1 ) to be able to observe PFPs within the snowpack.Sprinkling was performed in four bursts of 30 min each, interrupted by 30 min breaks.Sprinkling was conducted over a 2 × 2 m plot centred above the lysimeters, and with an intensity of 24.7 mm h −1 , leading to a total of 49.4 mm artificial rain in each of the experiments.The intensities were determined by calibration experiments on lysimeters not covered by snow and are valid for a certain distance between the nozzle and the sprinkled surface and water pressure at the nozzle.Despite the fact that this value still represents a very intense ROS event, it is within range of natural ROS events and similar or much lower compared to previous studies (19 mm h −1 ; Eiriksson et al., 2013; 48-100 mm h −1 ; Singh et al., 1997).For the sprinkling experiments, the exact timing of rain and intensities are known and the snowpack runoff measured at 1 min intervals allowed precise analysis of the performance of model simulations.Figure 1 shows a vertical cut of a snowpack after the sprinkling experiment and a top view of the lysimeter after the snowpack was removed for cold and wet conditions, respectively.The blue colour indicates where water transport took place and where sprinkled water was held by capillary forces or refrozen.

Extensive dataset for in situ validation
Two long-term datasets from two study sites in the European Alps providing snow lysimeter data and high-quality meteorological forcing data for running the energy balance model SNOWPACK were chosen to validate the different water transport models systematically.Datasets of both study sites used for the extensive in situ validation are publicly available.The CDP site, located in the Chartreuse range in southeastern France, has been described in Morin et al. (2012) and the Weissfluhjoch site (WFJ) in the Swiss Alps has been described in Wever et al. (2015).WFJ (46.83  snowpack produces snowpack runoff more often throughout the entire snow season and ROS events are more frequent than at WFJ.A multi-week snowpack builds up every winter season at CDP, but is, in contrast to WFJ, interrupted by complete melt in some years.The WFJ site is equipped with a 5 m 2 snow lysimeter, which measures the liquid water runoff from the snowpack.It has a 60 cm high rim to reduce lateral flow effects near the soil-snow interface (Wever et al., 2014a).CDP is equipped with both a 5 and a 1 m 2 lysimeter.Here we use data from the 5 m 2 lysimeter, but include data from the 1 m 2 lysimeter to discuss the uncertainty associated with measurements of the snowpack runoff.The studied period for WFJ is from 1 October 1999 to 30 September 2013 (14 hydrological years).Because of possible errors in the lysimeter data in the winter seasons of 1999/00 and 2004/05 as described in Wever et al. (2014a), these data were excluded from the study.For CDP the studied period is from 1 October 1994 to 31 July 2011 (17 winter seasons) according to the data availability from the 5 m 2 lysimeter.The temporal resolution of lysimeter data is 1 h for CDP and 10 min for WFJ.Simulation results for CDP and WFJ as well as lysimeter data for WFJ were aggregated to an hourly timescale.

CDP+WFJ event definition
As the number and characteristics of ROS events are strongly dependent on the event definition, special care needs to be taken to determine beginning and end of a ROS event.Being interested in the temporal characteristics of snowpack runoff during ROS, we need to include the entire period from the onset of rain to the end of ROS-induced snowpack runoff.Here we use an event definition according to Würzer et al. (2016) with slightly decreased thresholds to identify ROS events.According to this definition, a ROS event requires a minimum amount of 10 mm rainfall to fall within 24 h on a snowpack with a height of at least 25 cm at the onset of rainfall.While the event is defined to begin once the first 1 mm of rain has fallen, the event ends once there is less than 3 mm of cumulative snowpack runoff recorded within the following 5 h.This definition resulted in a selection of 61 events at CDP and 40 events at WFJ.The model simulations were subsequently evaluated over a time window that extends the event length by 5 and 10 h at the beginning and end, respectively (Fig. 2).These extended evaluation periods allowed us to also investigate a possible temporal mismatch between modelled and observed snowpack runoff.

Experimental sprinkling experiments
During the winter period 2014/15, six sprinkling experiments (Ex1-Ex6) were conducted on four different sites to be able to investigate snowpack runoff generation for different snowpack properties.With distinct differences in snowpack properties but controlled rain intensities, these experiments were expected to reveal the influence of snow cover properties and differences between the water transport models best.For all experiments, initial snow height (HS), snowpack temperature (TS) and LWC profiles were measured (Table 1 and Fig. 3).
According to these measurements, the snowpack conditions on which the sprinkling experiments were conducted can be  For all winter experiments (Figs. 4 and 5a, b, c), both modelled and observed total event runoff remained below the amount of sprinkling water.Energy input estimated by the SNOWPACK simulations suggests that snowmelt was insignificant for the winter experiments, but refreeze led to significant retention of liquid water.Additionally some sprin-kled rain was retained as LWC at the end of the experiments.During Ex3 no snowpack runoff was observed, visual inspection afterwards revealed an impermeable ice layer covering both the lysimeter and the adjacent ground.During spring conditions, on the other hand, snowmelt (5.1, 8.4 and 27.4 mm respectively) led to snowpack runoff exceeding total sprinkling input, except for measured snowpack runoff in Ex6 (Figs. 4 and 5d, e, f).
Additionally, Fig. 5 shows that only the PF model was able to reproduce all four peaks of observed snowpack runoff for winter conditions (Ex1 + 2), and even the magnitude of the first peak of Ex1 was captured well.For spring conditions however, all three models managed to represent four peaks corresponding to the four sprinkling bursts, but the PF model showed best correspondence with observed snowpack runoff (Figs. 4 and 5d, e, f; Table 1).Regarding the onset of snowpack runoff, the PF model especially led to faster snowpack runoff for the first two winter experiments, where the RE and BA models showed delayed snowpack runoff onset.For spring conditions the faster snowpack runoff response of the PF model led to a slightly early snowpack runoff.Maximal snowpack runoff rates for dry and cold conditions were generally overestimated by all models, whereas wetter conditions led to a minor underestimation (except for Ex3, where no snowpack was measured).
Regarding the overall correlation between measured and simulated snowpack runoff, PF outperformed the other models (Table 1), in particular during winter conditions.Summarizing, this initial assessment suggests that the PF approach has potential advantages in particular (a) as to the timing of snowpack runoff and (b) for cold snowpacks which are not yet entirely ripened.

Natural occurring ROS events
In January 2015, two ROS events occurred in the vicinity of Davos.They were observed over an elevational range of 950 to 1560 m a.s.l. on the same sites on which also the sprinkling experiments were conducted.Figure 6 shows the course of cumulative rainfall and snowpack runoff for both dates and all sites.Pre-event conditions (HS, LWC, TS) were measured shortly before the onset of rain for both events and are shown together with coefficients of determination (R 2 ) for hourly snowpack runoff of the different models Table 2.
For the event of 3 January 2015 (Fig. 6, upper row) the lower sites Serneus and Klosters (950 and 1200 m a.s.l.) showed a similar snowpack runoff dynamics regarding the delayed onset and the total amount (cumulative sum averaged over the three corresponding lysimeters: 20.3 and 21.1 mm, respectively).Also, the heterogeneity between data from the individual lysimeters was relatively low (Range of 3.1 and 3.9 mm, respectively).For the highest located site (Davos), however, the snowpack runoff measured by all four lysimeters showed a greater variability (Fig. 6c) in the delayed onset of snowpack runoff (0 to 7 h) and the total amount of snowpack runoff (mean 24.7 mm; range of 57.9 mm).The snow cover mostly built up within 1 week before the event.
Cold temperatures led to a light melt refreeze crust at the top, but no distinct ice layers were observed.For the lower sites (Serneus and Klosters), the PF and RE models generated snowpack runoff too early (PF: approx.3 h; RE: 0.2 to 1.4 h).The BA model generated snowpack runoff rather too late (1.3 to 2 h), but still within range of the variability of observed snowpack runoff for Serneus.However, the cumulative lysimeter snowpack runoff showed good accordance with modelled PF and RE snowpack runoff at Serneus, whereas PF led to an overestimation at Klosters and BA to an underestimation of cumulative snowpack runoff at all sites.At the higher-elevation site Davos, the RE model led to a better representation of mean observed snowpack runoff amount, when compared with BA and PF.The mean observed snowpack runoff onset however was represented best by the PF model (0.3 h early) when compared to the BA (3.4 h delay) and RE (1.1 h delay).
For the event of 9 January 2015 (Fig. 6, bottom row) the lower sites showed again little temporal and spatial heterogeneity in lysimeter runoff (range of 1 and 2.2 mm, respectively), whereas this was more the case for Davos again (range of 13.3 mm) probably owing to ice layers that were formed after the event on 3 January.Observed mean event snowpack runoff was more diverse for all elevations, where Klosters had the highest cumulative snowpack runoff (Serneus 13.3 mm; Klosters 17.7 mm; Davos 7.8 mm).If compared to observed total snowpack runoff, the PF model overestimated snowpack runoff for Serneus and Klosters, whereas the RE and especially the BA model underestimate event snowpack runoff for both sites.For Davos, all models were overestimating event snowpack runoff and led to early snowpack runoff.Apart from the RE model, which represented onset of snowpack runoff correctly for Serneus, none of the models were able to model snowpack runoff onset correctly for any of the sites.

Modelled and observed snowpack runoff for the whole dataset
Given the partly contradictory findings on the performance of the three model variants based on the above assessment for artificial ROS simulations under controlled conditions (Sect.3.1), as well as natural ROS events (Sect.3.2), further more systematic model tests were needed.Therefore we validate the different models based on extensive datasets from the two sites WFJ and CDP, as described in Sect.2.4.Before we focus on the specific performance of the PF model for a large number of individual ROS events, we first analysed the overall model performance throughout the whole study period, i.e. over entire winter seasons.For this, we analysed observed and modelled hourly snowpack runoff  provided snow heights exceeded 10 cm to ensure that lysimeter runoff was caused by snowpack runoff and not rainfall.For both sites, R 2 values for PF were slightly higher than for RE (Table 3), which both clearly outperformed the BA.The root mean squared errors (RMSEs) of the PF model were also lower compared to RE and BA.We can therefore con-clude that the implementation of the PF approach slightly improves water transport over entire winter seasons.

ROS event characteristics of the extensive dataset
Median characteristics of the individual ROS events at CDP and WFJ are summarized in Fig. 7.The temporal course of median rain and snowpack runoff rates of all events at WFJ (40 individual events) and CDP (61 individual events) are shown in Fig. 7a, b.ROS events at WFJ showed generally higher maximum rain intensities than at CDP, leading to higher median snowpack runoff intensities at the beginning of the events.Whereas at WFJ, ROS events tended to be short and intense, at CDP the event rainfall extended over a longer period of time.Interestingly, we observed relatively high initial snowpack runoff rates before the actual beginning of the ROS event, especially for WFJ, which suggests that many ROS events at this site occurred during the snowmelt period.Median snowpack runoff reached a peak after 1 and 3 h after the onset of rain for WFJ and CDP, respectively.At WFJ snowpack runoff and rain rates at the beginning of the events were generally higher than at CDP.The course of the median air temperature during ROS events at both sites is shown in Fig. 7c.Especially for WFJ, median air temperature (TA) dropped with the onset of rain and median TA was higher than at CDP.The mean initial ROS event snow height (HS) for WFJ was 95 cm, which is approximately the average snow height during mid-June (for 70 years of measurements).The mean initial HS for CDP is 67 cm.With a SD of 42 cm, the variability of initial HS for WFJ was higher than for CDP (29 cm).

Modelled and observed snowpack runoff at the event scale
Below we investigate the performance of the three water transport schemes at the event scale.Modelled snowpack runoff was assessed against observations by the coefficient of determination (R 2 ) and the RMSEs.To further analyse the representation of snowpack runoff timing, we defined an absolute time lag error (TLE) as the difference between the onsets of modelled and observed snowpack runoff in hours.The onset of snowpack runoff is defined as the time when cumulative snowpack runoff has reached 10 % of total eventsnowpack runoff.
Figure 8 shows box plots of R 2 (a, d), RMSE (b, e) and absolute TLE (c, f) for all 40 ROS events at WFJ (a, b, c) and 61 events at CDP (d, e, f), respectively.For both sites, R 2 values show that the BA model performance was inferior to the RE model which was in turn slightly outperformed by the PF model.The interquartile range of R 2 values for CDP was generally higher than for WFJ and increased from BA to RE, whereas it was decreasing for PF.The PF also led to a re-duction in RMSE by approximately 50 % if compared to the BA, but less (9 % for WFJ and 25 % for CDP) if compared to the RE model.Whereas the median of TLEs for all models at WFJ was 0 and therefore all models reproduced the onset of snowpack runoff very well, the interquartile range decreased from BA to the RE and PF models.The same behaviour in interquartile range decrease could be observed for CDP, where the magnitude of TLE was higher than for WFJ and mostly negative.The median TLE was again 0 for the PF and −1 h in the case of BA and RE, indicating that for these models, snowpack runoff was on average a bit delayed compared to the observations.For WFJ, TLE for BA was more often positive (early modelled snowpack runoff), which led to a very good median for BA, but also a larger interquartile range.Hence, the PF model showed the most consistent results, especially if regarding the interquartile range.For CDP we added the comparison between the 1 and 5 m 2 lysimeters installed at CDP (Sect.2.5) as a reference to Fig. 8, referred to as RL.This comparison can be seen as a benchmark performance, as it represents the measurement uncertainty of the validation dataset.As expected, RL shows the highest overall performance measures, but while the results for both PF and RE were reasonably close to those of RL, the BA model performed considerably worse.
The results shown in Fig. 8 may be influenced by both a time lag as well as the degree of reproduction of temporal dynamics.To separate both effects, we conducted a crosscorrelation analysis, allowing a shift of up to 3 h to find the best R 2 value.Figure 9 shows both the time lag, as well as the best R 2 value achieved.Interestingly, the BA model showed best correlations if the modelled snowpack runoff was shifted by 1 or 2 h (consistently too early compared to observations).The RE model, on the other hand, showed best correlations for a shift in the other direction (consistently too late compared to observations).Neither was the case for PF with lags centred around 0.
The R 2 of the cross-correlation analysis gives some indication of how well the temporal dynamics of the observed snowpack runoff can be reproduced, neglecting a possible time lag.The results in Fig. 9 show an improvement in R 2 values for both sites and all models if a time lag is applied.Greatest improvements were observed for the BA model for both sites.The good timing with the PF model is confirmed by almost no lag for WFJ and only a small lag for CDP needed to maximize R 2 .For CDP, both RE and PF had maximized R 2 values in range of the lysimeter comparison (RL).

Discussion
Even though PF of liquid water through snow is a phenomenon that has been known and investigated for a long time, it has not yet been accounted for in 1-D snow models that are in use for operational applications.The results of this study show that including this process into the water trans-Hydrol.Earth Syst.Sci., 21,[1741][1742][1743][1744][1745][1746][1747][1748][1749][1750][1751][1752][1753][1754][1755][1756]2017 www.hydrol-earth-syst-sci.net/21/1741/2017/  port scheme can improve the prediction of snowpack runoff dynamics for individual ROS events as well as for the snowpack runoff of entire snow seasons.Moreover, the represen-tation of the onset of snowpack runoff is improved.This is particularly important at the catchment scale, where a delay of snowpack runoff relative to the start of rain may affect the catchment runoff generation, especially if the time lag varies across a given catchment.
During the sprinkling experiments, sprinkling intensities were higher than average rain intensities during ROS but still within range of peak rain intensities during naturally occurring ROS events in the Swiss Alps (Rössler et al., 2014;Würzer et al., 2016) and the Sierra Nevada, California (Osterhuber, 1999).The use of the PF model clearly led to a better representation of the runoff dynamics for all experiments, including shallow and ripe snowpacks during spring conditions as well as cold and dry snowpacks representing winter conditions.The improvements were strongest for winter conditions, suggesting that under these conditions accounting for PF is most relevant.This is supported by observations of PFPs during winter conditions (Fig. 1a), which were not visible after the spring experiments.During winter conditions just a fraction of the lysimeter area was coloured with tracer, indicating PF of the sprinkled water (Fig. 1b), whereas spring conditions left the whole cross-section of the lysimeter coloured (Fig. 1c).While a fast runoff response can be expected for wet and shallow snowpack and may be easier to handle for all models tested, it is the cold snowpacks that both RE and BA models did not manage to represent well: runoff from these models was more than 1 h delayed (Ex1  1), the total event runoff of both RE and PF models is very similar for most conditions.Notably, the total event runoff for dry snowpacks is mostly overestimated by all models, suggesting an underestimation of water held in the capillarities.In cold snowpacks, dendricity of snow grains may still be high, such that water retention curves developed for rounded grains underestimate the suction.Additionally, high lateral flow was observed during the experiment for those conditions (Fig. 1a).This leads to an effective loss of sprinkling water per surface area of the lysimeter, which of course cannot be reproduced by the models.Therefore, observed snowpack runoff likely underestimates the snowpack runoff that would have resulted from an equivalent natural ROS event and we assume that the performance of the PF and RE models to capture the event runoff is probably better than reported in Table 1.Note that neglecting refreeze in the PF model should not be accountable for differences in the total event runoff between the RE and PF model, if we assume that the cold content is depleted by the end of the event.
Interestingly, despite having the coldest snowpack, time lag for the first natural ROS event at Davos was shorter than for the other two sites.This relationship where a cold and non-ripe snowpack with low bulk density led to smaller lag times was also found during sprinkling experiments conducted by Juras et al. (2016b).We assume that this is an indication for the presence of pronounced PFPs under those conditions, which is also supported by the high spatial variability of snowpack runoff.Glass et al. (1989) state that the fraction of PF per area is decreasing with increasing permeability, which itself was found to be increasing with porosity (Calonne et al., 2012).Therefore, with a decreasing PF area due to lower densities, the cold content of a snowpack loses importance, but saturated hydraulic conductivity is reached faster within the PFPs.The combination of those effects then is suspected to lead to earlier runoff.This behaviour should be ideally reproduced by the PF model and indeed the onset of runoff is caught well for this event.Here, our multilysimeter setup raises the awareness that the observed processes can show considerably spatial heterogeneity as documented, for example, in Fig. 6.The formation of ice layers also underlies spatial heterogeneity.Moreover, the creation of PFPs is strongly dependent on structural features like grain size transitions leading to capillary barriers.Unfortunately, no detailed information about grain size is available in the observations to verify this.
The PF model led to improvements in reproducing hourly runoff rates at CDP and WFJ for a dataset comprising several years of runoff measurements.This is an important finding, demonstrating that the new water transport scheme aimed at a better representation of PF during ROS events, did not negatively impact on the overall robustness of the model.To the contrary, the overall performance over entire seasons could even be improved.All three models represent the overall seasonal runoff better for WFJ than for CDP (Table 3), which was also found on the event scale (Fig. 8).Moreover, the CDP simulations exhibit a larger interquartile range in R 2 values and are therefore generally less reliable.The observed differences in model performance between both sites may either be caused by differences in snowpack or meteorological conditions or by issues with the observational data.Moreover, SNOWPACK developments have in the past often been tested with WFJ data, which could lead to an unintended calibration favouring model applications at this site.Despite an obvious contrast in the elevation of both sites, the average conditions during ROS events seem to vary. Figure 7 suggests that at WFJ short and rather intense rain events dominate.The higher maximum rain intensities at WFJ, compared to CDP, are probably due to the later occurrence of ROS at this site (May-June), where air temperatures and therefore rain intensities are usually higher than earlier in the season (Molnar et al., 2015).Regarding mean intensities over the event scale, data shown in Fig. 7 further imply that short and intense ROS events typically attenuate the rain input (ratio runoff to rain < 1), whereas long ROS events rather lead to additional runoff from snowmelt, which is in line with results presented in Würzer et al. (2016).
Snow height is generally higher at WFJ where the average initial snow height for the ROS events analysed was approximately 30 cm higher than at CDP. Ideally, the performance of the water transport scheme in the snowpack should not be affected by the snow depth.At both sites, the snowpack undergoing a ROS event is mostly isothermal with a mean initial LWC of 1.8 % vol (CDP) and 3.0 % vol (WFJ).The initial snowpack densities at both sites were quite different.At WFJ, densities for all ROS events are around Hydrol.Earth Syst.Sci., 21, 1741Sci., 21, -1756Sci., 21, , 2017 www.hydrol-earth-syst-sci.net/21/1741/2017/ 450-500 kg m −3 , whereas for CDP densities are spread from around 200 kg m −3 up to 500 kg m −3 .This suggests that the variable performance of all models at CDP (Fig. 8d) may be associated with early season ROS events.At CDP, a linear regression fit suggests a positive, albeit weak correlation between snowpack bulk densities and event-R 2 for the RE (R 2 of 0.2), but no correlation for both the PF and the BA model.It seems that the RE model had some difficulties with low-density snow, which was not the case for the PF model (Fig. 10).This may explain why PF outperformed RE at CDP, but not for WFJ.
Remaining inaccuracies in the representation of runoff for low densities for both models applying the RE may be explained by the fact that the water retention curve have been derived by laboratory measurements with high-density snow samples (Yamaguchi et al., 2012).The parameters defining the PF area (F ) have also been developed from snow samples with a density mostly above 380 kg m −3 (Katsushima et al., 2013).
We further analysed snowpack stratigraphy derived from the SNOWPACK simulations, such as marked grain size changes (bigger than 0.5 mm) and density changes (bigger than 100 kg m −3 ) in two adjacent simulated layers as well as the wet layer ratio (percentage of layers exceeding 1 % vol over layers below 1 % vol) and the percentage of melt forms (Table 4).These stratigraphy measures represent possible capillary barriers having implications on the single event-R 2 and might help understanding the advantages and disadvantages of the different models.Any considerable correlation between the abundance of stratigraphy features and event-R 2 would be indicative of potential errors in the respective model.Negative albeit small correlations could be found between the number of grain size changes and the event-R 2 for WFJ.Similar correlations were noted with regards to the number of changes in density between layers for the RE and PF model.In both cases correlations were less negative for the PF model indicating a more balanced and ultimately less degraded performance with increasing number of potential capillary barriers.While at WFJ most events occurred with ripe snow this was not the case for CDP.There, positive correlations were found between the ratio of melt forms and the wet layer ratio with event-R 2 for the RE model (Pearson's R of 0.57 and 0.66) and for the PF model (Pearson's R of 0.37 and 0.39).In this case the PF model also showed more balanced results that were less influenced by the initial LWC, which is in line with our findings of the sprinkling experiments.
System input rates (sum of melt rates and rain rates) are known to significantly affect water transport processes.For example, the area of PF (Eq. 1) is likely to depend on the water supply rate.Data using sandy soils from Glass et al. (1989), shown in DiCarlo ( 2013), suggest that with increasing system input rates the finger width of PF is increasing.Even though we have used the lowest influx rates from Katsushima et al. (2013), these rates still exceeded what seems representative of natural ROS events.We therefore analysed the effect of system input rates on the performance of our water transport models.Positive, albeit weak correlations (R 2 of 0.07 to 0.21) could be observed between event-R 2 and system input rates for all models, suggesting that they generally performed (slightly) better for higher influx rates.For the PF model this could probably be explained by the PF parameters depending on laboratory measurements with high influx rates.
In combination with the hydraulic properties for lowerdensity snow samples, additional laboratory experiments might be able to determine the number and size of PFPs for lower input intensities and snow densities.Especially the calibrated parameters threshold for saturation ( th ) and the number of PFPs for refreeze (N) could benefit from such experimental studies.Even though CDP and WFJ provide long-term measurements on an adequate temporal resolution, these data give little information about spatial variability of snowpack runoff limiting further validation opportunities.Large area multi-compartment lysimeter setups might help to improve estimating size, amount and spatial heterogeneity of flow fingers.Sprinkling experiments with preferably low sprinkling intensities on such a device could fill a knowledge gap about water transport in snow under naturally occurring conditions.

Conclusions
A new water transport model is presented that accounts for PF of liquid water within a snowpack.The model deploys a dual-domain approach based on solving the Richards equation for each domain separately (matrix and preferential flow).It has been implemented as part of the physics-based snowpack model SNOWPACK which enables us for the first time to account for PFPs within a model framework that is used operationally for avalanche warning purposes and snow melt forecasting.
The new model was tested for sprinkling experiments over a natural snowpack, dedicated measurements during natural ROS events, and an extensive evaluation over 101 historic ROS events recorded at two different alpine long-term research sites.This assessment led to the following main conclusions.
Compared to alternative approaches, the model accounting for preferential flow (PF) demonstrated an improved overall performance, particularly for lower densities and initially dry snow conditions.This led to smallest interquartile ranges for R 2 values and considerably decreased RMSEs for a set of more than 100 ROS events.When evaluated over entire winter seasons, the performance statistics were superior to those of a single domain approach (RE), even if the differences were small.Both PF and RE models, however, outperformed the model using a bucket approach (BA) by a large margin (increasing median R 2 by 0.49 and 0.48 for WFJ and 0.53 and 0.48 for CDP).In sprinkling experiments with 30 min bursts of rain at high intensity, the PF model showed a substantially improved temporal correspondence to the observed snowpack runoff, in direct comparison to the RE and BA models.While the improvements were small for experiments on isothermal wet snow, they were pronounced for experiments on cold snowpacks.
Model assessments for over 100 ROS events recorded at two long-term research sites in the European Alps revealed rather variable performance measures on an event-by-event basis between the three models tested.The BA model tended to predict too early onset of snowpack runoff for wet snowpacks and a delayed onset of runoff for cold snowpacks, whereas RE was generally too late, especially for CDP.Combined with results from a separate cross-correlation analysis, results suggested the PF model to provide the best performance concerning the timing of the predicted runoff.
While there is certainly room for improvements of our approach to account for PF of liquid water through a snowpack, this study provides a first implementation within a model framework that is used for operational applications.Adding complexity to the water transport module did not negatively impact on the overall performance and could be done without compromising the robustness of the model results.
Improving the capabilities of a snowmelt model to accurately predict the onset of snowpack runoff during a ROS event is particularly relevant in the context of flood forecasting.In mountainous watersheds with variable snowpack conditions, it may be decisive if snowpack runoff occurs synchronously across the entire catchment, or if the delay between onset of rain and snowpack runoff is spatially variable, e.g. with elevation.In this regard, accounting for PF is a necessary step to improve snowmelt models, as shown in this study.

Figure 1 .
Figure 1.(a) Vertical cut of a snowpack after the sprinkling experiment Sertig Ex3 (28 February 2015).Lateral flow and the presence of PFP were observed.PFP were generated at regions with rain water ponding at ice layers and layer boundaries with a change in grain size (creating capillary barriers).(b) Lysimeter area after sprinkling during winter conditions (Serneus Ex1, 26 February 2015): coloured areas indicate the area where water percolated due to PF. (c) Lysimeter area after sprinkling during spring conditions (Klosters Ex4, 26 March 2015): coloured area shows that water percolated uniformly, indicating dominating matrix flow.

Figure 2 .
Figure 2. (a) Example of a ROS event occurring at WFJ.The entire extent of the x axis refers to the evaluation period; the bar above the x axis refers to the event length.(b) Cumulative version of the plot.

Figure 3 .
Figure 3. Snow temperature and LWC profiles measured directly before the sprinkling experiment started.The lines represent observed ice layers (blue) and crusts (orange).

Figure 5 .
Figure 5. Rain and snowpack runoff displayed as hydrographs for the six sprinkling events.Ex1 (a)-Ex3 (c) were conducted during winter conditions, Ex4 (d)-Ex6 (f) were conducted during spring conditions.

Figure 7 .
Figure 7. Temporal course of median rain (a), measured snowpack runoff (b) and air temperature (c) for WFJ (dotted) and CDP (solid) aggregated over all 40 and 61 events respectively.The thinner lines represent the lower and upper quartiles, respectively.The displayed period is extended by 5 h prior to event commencement according to the event definition (0 h).

Figure 8 .
Figure 8. RMSE, R 2 and TLE for simulations of 61 ROS events at the CDP site and of 40 ROS events at the WFJ site for all models (BA, RE, PF) and the reference lysimeter (RL) available only for CDP.

Figure 9 .
Figure 9. Best R 2 values and corresponding lags using a crosscorrelation function allowing a time shift (lag) of max ±3 h.

Figure 10 .
Figure 10.Distribution of event-R 2 for CDP events for the PF (a, c) and RE (b, d) model.The sample is split into initial bulk snow densities above 350 kg m −3 (a, b) and below 350 kg m −3 (c, d).

Table 1 .
Snowpack pre-conditions and execution dates for the sprinkling experiments as well as R 2 values for the different model simulations.Measured values are snow height (HS), bulk liquid water content (LWC), bulk snow temperature (TS).No snowpack runoff measurements were available for Sertig (Ex3).

Table 2 .
Snowpack pre-conditions and R 2 for hourly snowpack runoff for natural events on 3 and 9 January.

Table 3 .
R 2 and mean absolute errors for hourly snowpack runoff for 17 and 14 years, for CDP and WFJ, respectively., and missed approx.10mm of snowpack runoff within the first hour of observed runoff.This can partly be explained by the fact that BA and RE need to heat up the subfreezing snowpack before they can generate snowpack runoff, whereas refreezing is neglected in the preferential domain of the PF model and runoff can occur even in a not yet isothermal snowpack.Adjusting parameters like the irreducible water content θ r for the BA model could probably lead to earlier runoff under these conditions, but thereby lead to earlier runoff, for example for WFJ events, where TLE already is positive for several events.Despite the improved representation of the temporal runoff dynamics of the PF model (Table

Table 4 .
Pearson correlation coefficients between event-R 2 and stratigraphic features at WFJ and CDP.Stratigraphic features are marked grain size changes (bigger than 0.5 mm) and density changes (bigger than 100 kg m −3 ) in two adjacent simulated layers as well as the wet layer ratio (percentage of layers exceeding 1 % vol over layers below 1 % vol) and the percentage of melt forms.