Solving Richards Equation for snow improves snowpack meltwater runo ff estimations

Introduction Conclusions References


Introduction
The presence of a snow cover has a strong impact on the hydrological cycle.The snow cover can delay the routing from precipitation to streamflow on time scales from a few hours or days to several months and many studies have shown the importance of the snow cover for accurate runoff and streamflow modelling (e.g.Seyfried et al., 2009;Koster et al., 2010;Mahanama et al., 2012).Furthermore, the effects of rain-onsnow (ROS) events on runoff from catchments are strongly dependent on the state of the snow cover (Marks et al., 2001;Mazurkiewicz et al., 2008).The use of physicallybased snowpack models in hydrological modelling is increasing, with varying degrees of model complexity.Studies have shown that by using physically-based snowpack models, the determination of (spatial distribution of) snow water equivalent and discharge at basin outlets improves (Marks et al., 1999;Zanotti et al., 2004).A correct description of water flow through a snowpack is not only important to improve meltwater runoff estimations, but also helps to improve the understanding of wet snow avalanche formation (Conway and Raymond, 1993;Mitterer et al., 2011).

N. Wever et al.: Solving Richards Equation for snow covers
An important coupling between the snow cover and surface or sub-surface flow is provided by meltwater runoff at the base of the snowpack.For an accurate assessment of the snowpack runoff, it is important to have a good understanding of water movement through a snowpack.Water transport in snow is a complex process, because in general, the snow cover consists of many layers that may vary in temperature, density as well as in grain size and type.Experimental studies have shown that microstructural properties of the snowpack strongly influence the hydraulic properties of snow (Shimizu, 1970;Colbeck, 1974;Marsh, 1991;Yamaguchi et al., 2010).Where a snow layer with small grains is on top of a layer with coarse grains, these differences in hydraulic properties will lead to the formation of capillary barriers at the interface (Jordan, 1996;Waldner et al., 2004;Peitzsch et al., 2008;Mitterer et al., 2011).Also lateral flow along these capillary barriers, ice lenses or other dense parts of the snow cover has been identified (Peitzsch et al., 2008).Field experiments also have repeatedly observed the existence of preferential flow paths, which can provide a more efficient water transport mechanisms than matrix flow alone (e.g.Kattelmann, 1985;Schneebeli, 1995;Katsushima et al., 2013).
For modelling the water flow through the snowpack, capillary effects are often neglected and many models follow a bucket-type approach (e.g.Flerchinger and Saxton, 1989;Jin et al., 1999;Marks et al., 1999;Boone and Etchevers, 2001;Vionnet et al., 2012).These water balance schemes are easy to implement and computationally very efficient, making them very suitable for distributed modelling.Richards Equation (Richards, 1931), hereafter denoted as RE, is often applied to soils to describe variably saturated matrix flow.However, its application potential for snowpacks has already been identified in literature (Colbeck, 1972(Colbeck, , 1974;;Jordan, 1983;Illangasekare et al., 1990;Hirashima et al., 2010b).Jordan (1996) and Hirashima et al. (2010b) found that RE was able to reproduce the formation of capillary barriers as observed in laboratory experiments and snow profiles.
Compared to the knowledge of hydraulic properties in various types of soil, measurements to derive hydraulic properties in natural snow are restricted to a few studies (Shimizu, 1970;Colbeck, 1974;Marsh, 1991;Yamaguchi et al., 2010).Except for the study by Colbeck (1974), these studies were focussed on wet snow.Other snow types are more difficult to investigate, because in the presence of liquid water, snow metamorphism is rapid and the transition to wet snow types almost immediate.
Several snowpack models that describe liquid water transport on the basis of RE have been developed.Some model studies were restricted to idealised snowpacks to quantify the effects of water percolation in snow (Illangasekare et al., 1990;Jordan, 1996;Daanen and Nieber, 2009).When neglecting capillary forces, the resulting gravity flow in RE, as proposed by Colbeck (1972) in the kinematic wave model, has been used in earlier versions of CROCUS (Brun et al., 1989) and the current version of SNTHERM (Jordan, 1991;Davis et al., 2001).Hirashima et al. (2010b) developed a water balance scheme where downward water flow by Darcy's law is approximated assuming stationary and equilibrium flow properties.To our knowledge, the full RE has not been used in a physically-based snow cover model to model water flow in snow in simulations of full snow seasons, forced by meteorological field measurements.In this study, the implementation of a solver for RE in the physically-based snowpack model SNOWPACK (Bartelt and Lehning, 2002;Lehning et al., 2002a, b) is described, where the snowpack and soil are treated as a continuous column.The model is applied on two alpine measurement sites with a seasonal snow cover and contrasting climatological conditions: Weissfluhjoch (WFJ) near Davos, Switzerland and Col de Porte (CDP) near Grenoble in France.In simulations of 14 and 17 snow seasons respectively, the bucket scheme, Hirashima et al. (2010b) scheme and the full RE are compared with lysimeter measurements of snowpack runoff.The focus will be primarily on the timing and magnitude of snowpack runoff.

Theory and numerical formulations
Several methods to model water flow in snow have been developed.The three that are compared here, will be discussed below.

Bucket scheme
In the so-called bucket scheme, a threshold water content (water holding capacity, θ h ) is defined, above which all the liquid water exceeding this threshold in a layer is transported downward in the snowpack or soil, regardless of the storage capacity of lower layers (Bartelt and Lehning, 2002).The downward moving water is either stored at one of the lower layers (if possible), or is drained from the model domain.In SNOWPACK, the water holding capacity varies per layer according to (Coléou and Lesaffre, 1998): where θ i is the volumetric ice content of the snow (m 3 m −3 ).

Richards equation
RE describes water movement in variably saturated porous media (Richards, 1931).For a 1-dimensional column, RE can be written in a mixed-form, which can be discretized in a finite difference approximation that ensures perfect mass balance (Celia et al., 1990): The Cryosphere, 8, 257-274, 2014 where θ is the volumetric liquid water content (LWC, m 3 m −3 ), K is the hydraulic conductivity (m s −1 ), h is the pressure head (m), z is the vertical coordinate (m, positive upwards and perpendicular to the slope), γ is the slope angle and s is a source/sink term (m 3 m −3 s −1 ).For applying Eq. ( 2), a water retention curve has to be specified that relates θ to h.For snow, it is common to take the van Genuchten (1980) model: where θ r is the residual water content (m 3 m −3 ), θ s is the saturated water content (m 3 m −3 ) and α is a fit coefficient that is related to the maximum pore size in the medium.m and n are additional fit parameters that are related to the pore size distribution.Sc is a correction factor proposed by Ippisch et al. (2006), who have shown the necessity of using an air entry pressure when n ≤ 2. Here, the correction is applied for all values of n, with an air entry pressure of 0.0058 m, corresponding to a largest pore size of 5 mm.Two parameterizations for the van Genuchten (1980) model have been published recently: Daanen and Nieber (2009) fitted the water retention curve based on several published experiments that were collected by Marsh (1991) and Yamaguchi et al. (2010) retrieved water retention curves from several snow samples in laboratory experiments.Daanen and Nieber (2009) determined the parameters α and n for the van Genuchten model to be: and where r g is the mean grain radius (mm), corresponding to the classical grain size definition (Fierz et al., 2009).In this study we (primarily) refer to this classical grain size of a snow layer as the average size of its grains, where the size of a single grain or particle is its greatest extension measured in millimetres (Fierz et al., 2009;Baunach et al., 2001).This has to be distinguished from the optical equivalent grain size that is related to the specific surface area of snow (see, for example, Carmagnola et al., 2013;Fierz et al., 2009).
The parameterization for α proposed by Yamaguchi et al. (2010) is: For n, the original parameterization by Yamaguchi et al. (2010) was modified by Hirashima et al. (2010b) to be able to extend the parameterization beyond grain radii of 2 mm: where, for numerical stability, the upper grain radius limit is set to be 5 mm in this study.
In both parameterizations, the van Genuchten parameter m is chosen as: Note that both parameterizations for the van Genuchten model for snow differ significantly.Therefore, both are taken into consideration in this study, where the Hirashima et al. (2010b) modifications of the Yamaguchi et al. (2010) parameter set will be referred to as Yamaguchi and the Daanen and Nieber (2009) parameter set as Daanen.
To apply the van Genuchten model, also θ s and θ r need to be defined.For θ s , it is common to take it equal to the pore space: Note that the correction factor at the right hand side ensures that there is enough space when liquid water freezes and thereby expands.For θ r , Colbeck (1974) experimentally found a value of 0.07.Gravity drainage experiments by Yamaguchi et al. (2010) showed values around 0.02, suggesting that additional suction could reduce θ even more.Cold fresh snow is initially completely dry and also phase changes can reduce θ below θ r , causing singularities in the van Genuchten model (Eq.3).
To circumvent these problems, a pragmatic strategy was implemented for determining θ r : first, new snow layers were initialized with a very small value for θ of θ /10, where θ is the convergence criterion (see Appendix A).This means that new snow layers are initialized in a dry state from the perspective of the convergence criterion.Then, the residual water content for the current time step (θ t r ) is scaled between 0 and 0.02 according to: where f is a tuning factor, taken as 0.75 and θ t−1 r is the residual water content from the previous time step.
For determining the hydraulic conductivity, the van Genuchten-Mualem model is used (Mualem, 1976;van Genuchten, 1980): where K sat is the saturated hydraulic conductivity (m s −1 ) and is the effective saturation: For K sat , the proposed equation by Shimizu (1970) has been widely used for a long time: where ρ w and ρ i are the density of water (1000 kg m −3 ) and ice (917 kg m −3 ), g is the gravitational acceleration (taken as 9.8 m s −2 ) and µ is the dynamic viscosity (taken as 0.001792 kg (m s) −1 ).Recently, Calonne et al. (2012) proposed a slightly different equation, based on calculations from tomography images: where r es is the equivalent sphere radius (m).
The largest difference between the parameterization of Shimizu (1970) and Calonne et al. (2012) is in the range of snow with low density and new snow or decomposing and fragmenting snow particles (Calonne et al., 2012).In the presence of liquid water, wet snow metamorphism is rapid and also densification occurs, which quickly reduces the differences between both parameterizations.We found that results marginally improve when using the equation from Calonne et al. (2012), so results presented here exclusively use Eq. ( 14).For this, we determine the optical grain radius from the simulated grain size following Vionnet et al. (2012) and use this as an approximation of r es .
In SNOWPACK simulations, the typical microstructural evolution of the snowpack results in an increase in snow grain radius when snow ages or experiences melt-freeze cycles.This will increase the saturated hydraulic conductivity.However, settling of the snow cover and subsequent melt cycles will increase the density, which will decrease the hydraulic conductivity.In the simulations, the first effect is dominating, as grain size will generally increase by a factor 3-8, thus influencing K sat by a factor 9 to 64 over a snow season.Snow density will increase by a factor 2-5, influencing K sat by a factor 5 to 25 over a snow season for typical snow densities.Field experiments have shown that different temperature regimes may result in a different microstructural evolution of the snowpack and hence of K sat (Domine et al., 2013).Hirashima et al. (2010b) developed a water balance scheme for the SNOWPACK model based on approximating the water transport by fluxes derived from Darcy's law by assuming stationary flow properties for a time step and equilibrium in water flow between two layers.RE is not explicitly solved and only downward water movement is possible.Hirashima et al. (2010b) also used the earlier discussed modified water retention curve for the van Genuchten model by Yamaguchi et al. (2010).The study focussed primarily on the internal snowpack processes and achieved the reproduction of capillary barriers on the interface between layers with different grain sizes.In this study, we will refer to this water balance scheme as NIED.

NIED scheme
3 Data and methods

Data Weissfluhjoch
The experimental site at the Weissfluhjoch (WFJ, 46.83 • N, 9.81 • E) is located at an altitude of 2540 m in the Swiss Alps near Davos.During the winter months, almost all precipitation falls as snow at this altitude.As a consequence, a continuous seasonal snow cover builds up every winter, with a maximum snow height ranging from 153-366 cm over the period 1934-2012.The measurement site is located in an almost flat part of a southeast oriented slope (γ is taken equal to 0).
The dataset contains air temperature, relative humidity, wind speed and direction, incoming and outgoing longwave and shortwave radiation, surface temperature, soil temperature at a few cm below the surface, snow height and precipitation from a heated rain gauge.The dataset has been quality checked, by replacing missing values with values from backup sensors or by applying various interpolation methods (Schmucki et al., 2013).The precipitation was corrected for undercatch during snowfall and windy conditions by using a correction proposed by Førland et al. (1996).The relationship for a Finnish H&H-90 gauge was taken, as it was empirically found to estimate precipitation at WFJ best, when comparing seasonal maximum measured snow water equivalent (SWE) in biweekly manual snow profiles to modelled SWE at the same date.
The experimental site is equipped with a lysimeter, which measures the liquid water runoff from the snowpack.The surface area of the squared lysimeter is 5 m 2 and since 2001, it is measuring at a resolution of 0.16 kg m −2 (before 0.20 kg m −2 ).Data is collected at 10 min intervals.The lysimeter is enclosed by a 60 cm high wall to reduce lateral flow effects near the soil-snow interface.However, lateral flow along capillary barriers within the snow cover or other, more impermeable layers (e.g.melt-freeze crusts) higher in the snowpack may still affect the measurements.The lysimeter is expected to collect additional meltwater advected by lateral flow mechanisms, as the device is located close to the lowest part of the site.
The studied period for WFJ is from 1 September 1996 to 30 September 2010.This period, consisting of 14 full winter seasons, is chosen based on data-availability and quality of the lysimeter measurements.

Data Col de Porte
The experimental station Col de Porte (CDP, 45.30 • N, 5.77 • E) is located in the Chartreuse range in France and is operated by MeteoFrance.The datasets of hourly measurements and the snowpack evaluation data have been published and described in Morin et al. (2012).The station is located at 1325 m altitude and, as a consequence, experiences a warmer climate than WFJ.The snow cover produces snowpack runoff more often during the snow season and ROS events are more frequent than at WFJ. Due to these frequent wetting cycles, we expect the snowpack to be less stratified at CDP than at WFJ.A snow cover builds up for at least several weeks every winter season at this site, but is interrupted by complete melt in some years.With typical maximum snow heights ranging from 68-207 cm over the period 1993-2011, snow covers at CDP are generally less deep than at WFJ.
The dataset contains a set of very similar measurements as at WFJ: air temperature, specific humidity, wind speed and direction, incoming and outgoing longwave and shortwave radiation, surface temperature, several soil temperatures, snow height and precipitation from a heated rain gauge.The precipitation data in the CDP dataset is already corrected for undercatch and the precipitation amounts are separated in rain and snow, as described in Morin et al. (2012).The SNOWPACK model accepts only a single value as precipitation input and decides based on a temperature threshold of 1.2 • C whether the precipitation should be considered rain or snow.Therefore, we adjusted the CDP dataset to an air temperature above 1.2 • C in case only rainfall was reported, and an air temperature below 1.2 • C in case snow fall was reported.Mixed precipitation was provided as a sum to SNOW-PACK and depending on the air temperature considered as rain or snowfall.
CDP is equipped with both a 5 m 2 and a 1 m 2 lysimeter.Kattelmann (2000) shows that smaller lysimeters exhibit significantly larger inaccuracies.Also in this study we found that the larger lysimeter provided highest agreement for all water balance schemes, and we will only use the data from the 5 m 2 lysimeter.The dataset from CDP is used here for the period from 1 September 1994 to 31 July 2011 (17 winter seasons), determined by the data availability of the 5 m 2 lysimeter.

Model description
The 1-dimensional mixed form of RE (Eq.2) is solved by a Picard iteration scheme (Celia et al., 1990), adapted for a variable grid spacing based on Rathfelder and Abriola (1994).The scheme has the property of an easy numerical implementation that is globally convergent.In Appendix A, several aspects of the numerical implementation are discussed.
When using RE, the snowpack and the soil are considered as a single continuous column.There are no special boundary conditions for the lowest snow layer or upper-most soil layer.The snowpack runoff, defined as all liquid water leaving the snowpack at the bottom, is calculated in a model postprocessing step by evaluating Darcy's law at the interface between the snowpack and the soil: where q is the snowpack runoff (m s −1 ), positive values denoting downward water movement.The right hand side describes the numerical implementation, where i denotes the upper-most soil element and i + 1 the lowest snow element, K i+1/2 is the hydraulic conductivity at the interface between element i and i + 1, and h i and h i+1 are the pressure heads in the top soil and bottom snow element, respectively.z is the vertical grid spacing.In the rest of the paper, snowpack runoff will be treated from a mass balance perspective, denoting downward water movement (snowpack outflow) with a negative value.
For the snowpack, the layer thickness varies per layer and with time.When there is solid precipitation, new snow layers are added on top of the domain with an initial thickness of 2 cm.Over time, settling of the snow reduces the layer thickness.When a snow layer gets smaller than a specified minimum height or the ice content decreases below a specified minimum value, it is joined with the layer above (in case it is the lowest snow layer) or below (all other cases).When two adjacent snow layers get similar properties, they are also merged (Bartelt and Lehning, 2002).
The model is run in 15 min time steps.Processes are described sequentially, assuming stationary snowpack behaviour in these 15 min.First, new snow is added to the snowpack at the beginning of a time step when solid precipitation occurs.Then, the change in temperature distribution over the time step is calculated, after which phase changes are executed based on the new temperature profile at the end of the time step.Then the water balance routine is executed.For runs with the bucket or NIED scheme, first the snowpack water flow is calculated.The sum of snowpack runoff over the time step is expressed as an average runoff from the snowpack over the time step, which is then used as a specified flux (second type) boundary condition (McCord, 1991) for RE for the soil part.For runs with RE scheme for both snowpack and soil, the soil-snow column is calculated as a continuous column.After the calculation of water transport is finished, the new snowpack state will undergo phase change again when necessary, mainly to freeze percolating meltwater and rain water.The time step is finalised by calculating the internal snowpack metamorphism, as described in Lehning et al. (2002b).Note that water transport by the NIED scheme is internally calculated with a time step of 1 min.The solver for RE uses a variable time step (see Appendix A for details).
The latent heat associated with rain, evaporation/condensation and sublimation/deposition is applied as a Neumann boundary flux for the temperature equation.In case of rain, evaporation or condensation, the liquid water flux itself is used as a specified flux boundary condition for the water balance schemes.During deposition, surface hoar www.the-cryosphere.net/8/257/2014/The Cryosphere, 8, 257-274, 2014 may be formed on the snow cover (Stössel et al., 2010).
In case of evaporation or sublimation, the implementation of the boundary condition require different treatments, depending on water balance scheme.In the bucket and NIED simulations, the latent heat flux is first used to evaporate liquid water from the snowpack.Remaining energy is then used to sublimate the ice matrix, respecting the fact that evaporation is more favoured by the lower energy required for evaporation.When RE is used for the snowpack, the specified flux at the upper boundary is actively limited.
A maximum evaporative flux is allowed, based on an imaginary limiting pressure head outside the model domain (see Appendix A).If the evaporative flux exceeds this value, the flux is limited to this value.The excess energy is used to sublimate the ice matrix.Note that the actual mass exchange differs whether evaporation/condensation or sublimation/deposition occurs, because of the difference in latent heat involved.This may cause a slightly different mass balance of the snowpack in the different simulations.

Simulation setup
The SNOWPACK model is driven by the measured precipitation amounts.Measured precipitation is assumed to be rain when the air temperature was higher than 1.2 • C and snow otherwise.The net radiation is determined from the incoming and outgoing longwave and shortwave radiation and, consequently, albedo is derived from the measured shortwave radiation components.Absorption of penetrating shortwave radiation is treated as a source term for the temperature equation in the top snow layers.The surface energy balance is evaluated using the net longwave radiation and calculated latent and sensible heat fluxes using a common form of Monin-Obukhov bulk formulation (Lehning et al., 2002a).Atmospheric stability is estimated from the difference in measured air and surface temperature.The net heat flux is used as a Neumann boundary condition for the temperature equation.
In the simulations, a roughness length (z 0 ) of 0.002 m was used for WFJ (Stössel et al., 2010) and 0.015 m for CDP (Essery et al., 2013).In all simulations, the water flow in the soil is solved using RE.Soil freezing was not considered and the soil temperature was bounded by 0 • C.This assumption is supported by snow profiles of seasonal snow covers, where ground heat creates close to melting conditions at the snow-soil interface.
For WFJ simulations, a soil of 1.5 m depth is used, divided into 30 layers of varying thickness.This setup ensures that choices for lower boundary conditions in the soil are only marginally impacting the snow cover.Typical soil properties for gravel/coarse sand were taken (θ r = 0.01, θ s = 0.35, α = 3.5 m −1 , n = 4.5 and K sat = 3.171 • 10 −6 m s −1 ).At the lower boundary, a water table is prescribed as a Dirichlet boundary condition for RE and a constant upward soil heat flux of 0.06 W m −2 is prescribed as Neumann boundary condition for the temperature equation.Although in reality, the water table is expected to be deeper at WFJ, the gravel/coarse sand material will make the influence of the water table on snowpack runoff negligible small.
At CDP, the ground heat flux is an important factor in the energy balance of the snowpack, due to the low elevation and consequently warmer soil of the site.To take this into account in the simulations, a soil of 10 cm depth was used, divided into 10 layers, and measured soil temperatures at 10 cm depth were applied as a Dirichlet boundary condition for the lower boundary.Although Morin et al. (2012) reports a substantial loam content for the soil at CDP, we used the same soil as for WFJ to prevent ponding conditions in the shallow soil.At the lower boundary, a free drainage Neumann boundary condition was applied for solving RE.

Analysis
For both sites, the model simulations and lysimeter measurements are analysed for the 24, 12, 6, 3 and 1 h time scales and, only for WFJ, also for the 0.5 h time scale.Runoff values were constructed by summing the 15 min model output resolution or the 10 min (WFJ)/1 h (CDP) lysimeter measurement resolution to the respective time scales, starting at 00:00 h (midnight) local time.We restrict the analysis in this study only to periods with a snow cover on the ground during the winter season.We will refer to these periods as snow seasons, where the snow seasons are denoted by the year in which they end (e.g.1997 denotes winter season 1996-1997).Summer snowfalls are ignored in the analysis.
To determine the beginning and end of a snow season at WFJ, it is assumed that on 1 March a snow cover is always present.Then, it is searched both forward and backward in time in the measurements and simulations until all have snow-free conditions.These mark the start and end of the snow season.The melt season at WFJ is defined here as the period from 1 March to the melt out date.In the measurements from the lysimeter, one cannot distinguish snowpack runoff with a snowpack present and rain without a snowpack.Using measured snow height to determine the end of the snow season appeared to be somewhat cumbersome as measurement inaccuracies make it difficult to determine when the surface is snow free.Moreover, the rain gauge used to derive precipitation amounts, the snow height sensor and the lysimeter are located several metres apart.Given the spatial variability of the snow cover thickness, the snow height sensor cannot be regarded as fully representative for the snow height at the lysimeter.Therefore, the snowpack runoff measured by the lysimeter was assumed to come from the snowpack as long as a snow cover was present in the simulations and consequently, the lysimeter runoff may include rainfall in snow-free conditions, if the duration of the simulated snow cover exceeds the real snow cover presence.
At CDP, the snow cover may melt completely during the winter season in some years.Therefore, we define the start of the snow season here as the first snowfall that contributed to The Cryosphere, 8, 257-274, 2014 the snow cover present at 31 December.If at May 1, no snow cover is present, it is looked backward in time to find the last day where a snow cover is present in either the measurements or one of the simulations.Else, we look forward in time, to find the last day a snowcover is present.The lysimeter data from CDP is cleaned by the providers such that it only reports values when a snow cover is present.However, if the snow cover still exists in the simulations and the lysimeter is already snow-free, rainfall will be added to snowpack runoff in the simulations, but is not reported by the lysimeter.Therefore, rainfall was added to the lysimeter series when there was no measured snow cover, but still a snow cover in one of the simulations.
To assess the added value of the water balance schemes, the performance of modelled snowpack runoff is compared to a constructed runoff series that consists of the modelled LWC production (snowmelt, refreezing and rain input).This LWC production is calculated for each time step as the average of the four simulation variants.It represents the case where the sum of snowmelt and rain input is routed to runoff immediately and it will be referred to as average LWC production (avg.LWC prod.).
To quantify the accuracy of the modelled snowpack runoff, Nash-Sutcliffe model efficiency coefficients (NSE, Nash and Sutcliffe, 1970) are calculated for all time scales.These are only determined over the winter snow season.To calculate an NSE coefficient over all snow seasons, the snow seasons are analyzed sequentially by removing the intermediate summer periods.
All simulations were run on the same desktop computer as a single-core process.The increase in computational time by using RE for the snow cover is significant but not excessive: the additional effort for water transport is of the same order as the total seasonal snow cover simulation without RE.For WFJ, the simulations with the bucket scheme for snow took on average about 1.5 min per year and the NIED scheme took about 1.7 min per year.Solving RE for snow increased the average simulation time to 3.1 min per year for RE (Yamaguchi) and 2.7 min per year for RE (Daanen).CPU times for CDP are about half as those for WFJ, as a result of the shallower snowpack.

Results and discussion
The discussion of the results will first focus on the seasonal and daily time scale and afterwards on sub-daily time scales.Table 1 shows the snow season period and the maximum measured snow height for each year.In the studied period, the snow season at WFJ is mostly starting in October, and the melt out date is mostly between mid-June and mid-July.The maximum measured snow height ranges from 182 cm in 2005 to 356 cm in 1999.At CDP, the snow season is shorter than at WFJ and is roughly between November and April/May.Also the maximum snow height reached is lower than at WFJ, and ranges from 68 cm in the years 2001 and 2007 to 207 cm in 1999.
Table 1 also shows the seasonal average difference between measured and modelled snow height.The year-toyear variability indicates that using measured precipitation to drive snowpack simulations has some inherent inaccuracies and the main problem is that the over-or underestimation of a single precipitation event will persist throughout the rest of the snow season.It is noteworthy that snow height also varies between the different schemes, although these variations are smaller than the year-to-year ones.The snow height variations between simulations can be explained by differences in mass loss due to runoff and by differences in settling, as a result of varying vertical LWC and density distributions due to the different water balance schemes.Although the differences between simulations are rather small, RE with Yamaguchi's parameterization shows the smallest deviation from measured snow height, whereas Daanen's parameterization seems to cause slower settling and higher snow heights.This is especially visible in the simulations for WFJ.
Figure 1 shows that the cumulative runoff sum from the snowpack exhibits much smaller variations between simulations than the snow height, indicating that there is only a small difference in modelled snow water equivalent.As the precipitation input is the same for all simulations, the slightly varying cumulative runoff sums are mainly caused by differences in evaporation/condensation or sublimation/deposition.The reasons for these differences are, first, that the LWC in the surface layer will influence the surface temperature, especially due to variations in time needed to refreeze at night.These differences in surface temperature in the evening hours may influence sensible and latent heat exchange and may cause differences in energy input into the snow cover between simulations.Second, the partitioning of latent heat between sublimation and evaporation is different between the water balance schemes.Because of the difference in latent heat associated with sublimation or evaporation, mass gain or loss will be smaller in case of sublimation/deposition than in evaporation/condensation.Note that in practise, the albedo of the snow cover, and thus the energy balance, is also influenced by the LWC in the surface layer.This influence is not present here, because measured albedo is used to drive simulations at both sites.
At both sites, the agreement between measured and modelled cumulative sum of snowpack runoff is limited (see Fig. 1), although most year-to-year variability is captured.The undercatch correction was verified for both sites by comparing the seasonal maximum measured SWE in snow pits to the modelled SWE at the same date.No systematic over-or underestimation was found, which suggests that the estimation of precipitation does not show a bias.The overestimation of snow height at WFJ and the underestimation of the snow height at CDP (see Table 1) are likely related to the modelling of new snow density and snow compaction over time.The discrepancies between modelled and measured runoff may also be caused by an unrepresentative snowpack state at the lysimeter or inhomogeneous flow features like flow fingering and flow over ice lenses.The accuracy of the undercatch correction for solid precipitation (Goodison et al., 1998) may also vary from year-to-year.Furthermore, the deviations are likely a general expression of the fact that snow height can vary over short distances, mainly caused by heterogeneity in wind fields (e.g.Winstral et al., 2002;Mott et al., 2010).However, the deviation between measured and modelled snowpack runoff is suspiciously large at WFJ in the years 1997 and 2000.These differences seem to be too large to be explained by lateral flow effects or inhomogeneous snow redistribution, and the agreement between mea-sured and modelled snow height for these years (see Table 1) suggest that precipitation input is also not the cause of the differences.In both cases, the consistency between maximum measured snow height and modelled snowpack runoff indicates that the most probable source of error is in the measured snowpack runoff.This likely is either a malfunctioning of the lysimeter or a strongly unrepresentative snowpack state at the lysimeter.the modelled runoff sum.Therefore, this year is used as an example.Figure 2a shows the cumulative snowpack runoff in the melt season.As can be seen, the lysimeter registers the first melt water at the base of the snowpack earlier than any of the model schemes.The simulations with RE produce runoff soon after the first measured runoff, whereas the bucket and NIED simulations show some delay.For the rest of the melt season, there are no important differences.Because the bucket and NIED simulations withhold the water too much in the snowpack compared to the lysimeter and the simulations with RE, the daily outflow near the end of the season becomes higher than in simulations with RE.General runoff dynamics and alternating phases with high and low runoff are represented well in all simulations, although in some periods, discrepancies in runoff amounts between measurements and model exist.

Daily time scale
Figure 3 shows the NSE coefficients for daily sums of snowpack runoff for both WFJ and CDP for the respectively 14 and 17 yr individually.It shows that for the daily time scale, differences between the various models are much smaller than year-to-year differences in NSE coefficients.Furthermore, the spread between simulations is larger at CDP than WFJ.For both sites, it can be seen that when the agreement between modelled and measured snowpack runoff is lower, all water balance schemes have a lower agreement with measured snowpack runoff (and vice versa).The fact that NSE coefficients for the average LWC production fol- low the same pattern is an indication that for these years, sources of error causing variation of NSE coefficients may be related to the estimation of meltwater production as determined by a positive energy balance and a possibly sometimes inaccurate partitioning of precipitation in rain and snow.The representation of internal snowpack structure and accompanying hydraulic properties in the model seem to play a less pronounced role, which may also imply that the model does not reproduce adequately enough year-to-year variations in internal snowpack structure that influence water flow.
On the other hand, lysimeter measurements are known to be notoriously difficult, as found in experiments by Kattelmann (2000).Discrepancies between measured and modelled snowpack runoff should therefore not be attributed solely to an inaccurate representation of the snow cover or energy balance in the model.For example, due to lateral flow in the snowpack, the effective area from which meltwater is collected may differ from the actual area of the lysimeter (Kattelmann, 2000).Furthermore, the walls around the lysimeter to prevent preferential flow along the base of the snowpack may influence the snow cover at the start of the snow season.The lysimeter may collect more snow due to wind effects and snowmelt can be reduced due to the shadowing effect of the wall.This may lead to a nonrepresentative snow cover inside the lysimeter.
The three years with a very low NSE coefficient for WFJ are also likely related to measurement problems.Typical causes for low NSE values are a consistent over-or underestimation (bias) and poor timing of meltwater peaks (Mc-Cuen et al., 2006).For the years 1997 and 2000, the low NSE coefficients are for an important part caused by the bias between modelled and measured seasonal snowpack runoff (see Fig. 1).For the year 2005, it appears as if the lysimeter was obstructed for quite some time, after which half the seasonal sum of snowpack runoff passed through the device in few days time (not shown).
In terms of NSE coefficient, RE does reproduce measured snowpack runoff best at both sites, although differences are marginal.At WFJ, RE achieves an NSE coefficient of 0.63 and 0.62 for Yamaguchi's and Daanen's parameterization, respectively, over all 14 snow seasons.Bucket and NIED schemes also exhibit very similar performance with NSE coefficients of 0.61 and 0.62, respectively.For CDP, NSE coefficients for the daily time scale are higher than for WFJ.RE achieves an NSE of 0.72 and 0.71 for Yamaguchi's and Daanen's parameterization, respectively.NSE coefficients for the bucket and NIED scheme are 0.70 and 0.72, respectively.The fact that the different water balance schemes have very similar performance on the daily time scale and have an NSE coefficient very similar to the one for the average LWC pro-duction is an indication that once the snowpack is isothermal, meltwater production near the surface is transported downward within the day in all models and that this is in good agreement with the measurements.This is noted already in literature (e.g.Colbeck, 1972;Brun et al., 1989;Davis et al., 2001).

Timing of seasonal runoff
To assess the performance of the water balance schemes in simulating the dynamics of the snowpack runoff over the season, it was determined at which date a certain percentage of the total cumulative snowpack runoff was reached.This was also done for the measurements.Then the difference between measured and modelled date was determined.Figure 4 shows the time delay in days between modelled and measured snowpack runoff.A positive delay means that modelled cumulative snowpack runoff arrives later in time than the measured one.For example: 20 % of total snowpack runoff is delayed by 1 day in the RE (Yamaguchi) simulations and 4 days in the bucket simulations.For CDP, the comparison is made over the complete snow season, as important melt and snowpack runoff events occur throughout the season.At WFJ, melt and subsequent snowpack runoff occurs either at the early start of the snow season (October or November), or after the beginning of March.So for WFJ, the comparison is made for the spring melt season only.
As can be seen, all models are strongly underestimating the arrival date of the first few percent of the total snowpack runoff in the melt season at WFJ, up to about 18 days for the bucket and NIED scheme.Average LWC production in the models is in fairly good agreement with measured snowpack runoff, suggesting that the first surface melt is almost directly accompanied by snowpack runoff.This fast arrival of meltwater at the base of the snowpack in the measurements likely results from the more efficient transport mechanism of preferential flow paths compared to matrix flow.Several experiments have shown that water flow in snow is not horizontally homogeneous (e.g.Conway and Benedict, 1994;Waldner et al., 2004;Katsushima et al., 2013).The 1-D approach in this study cannot resolve preferential flow paths by flow fingering, as observed in several experiments.Preferential flow paths will be able to transport water downward faster than horizontally uniform matrix flow, as simulated here.However, Fig. 4 suggest that this may involve only about 5 % of cumulative seasonal snowpack runoff.Note that this does certainly not imply that preferential flow cannot have a more pronounced effect on the internal snowpack microstructure or wet snow avalanche formation.
After the start of snowpack runoff, the simulations with RE quickly show very little delay with measured snowpack runoff.For values of about 5 % and above, the dynamical snowpack runoff during the melt season is adequately simulated with RE.The fact that the delay is fairly constant after about 25 % shows that the daily amount of meltwater leaving the snowpack is in quite good agreement with measured values.The delayed snowpack runoff in the bucket and NIED scheme persist for almost the entire season.Apparently, bucket and NIED are retaining meltwater in the snowpack too long, releasing it later in the melt season.
At CDP, the first modelled snowpack runoff in the snow season is in rather good agreement with the measurements.This is likely a result of the fact that the snowpack is still very shallow when these snowpack runoff events occur.Afterwards, there arises a marked delay, visible in all water balance schemes, but most pronounced in the bucket scheme.In this period, extending from about 8-40 % of the snowpack runoff, more snowpack runoff is measured than modelled.This may be an expression of efficient preferential flow paths.After about 40 %, all water schemes are fairly constant compared to the beginning of the snow season.It may be concluded that the variations in the RE and NIED schemes are smaller than in the bucket scheme, showing that they provide better representation of the runoff dynamics throughout the snow season.

Sub-daily time scales
Figure 2b shows the hourly flux of snowpack runoff for one week during the melt season of the example snow season 1998.A daily returning peak in melt water outflow, associ-ated with the daily cycle in melt, is visible.During 3 June a ROS event occurred with 14 mm of rain.The daily onset of runoff and timing of the peak flux is better reproduced by simulations with RE.The figure also shows that RE is able to reproduce the recession curve in the evening hours and the night.Although the NIED scheme does not reproduce this for WFJ in this case, Hirashima et al. (2010a) found that the recession curve is reproduced in the NIED scheme for warm snow regions such as the central part of Japan.
In Fig. 5, NSE coefficients for the various water balance schemes are shown for sub-daily time scales.It is clear that for smaller time scales, the NSE coefficients decrease for all water balance schemes.On the 12 and 24 h time scales, all schemes produce more or less similar results at both sites.For the 1 h time scale, RE (Yamaguchi) achieves a still reasonable NSE coefficient of 0.49 for WFJ, where the bucket and NIED scheme have NSE coefficients of 0.17 and 0.07, respectively.The decrease in NSE coefficient for the bucket and NIED schemes must be mainly caused by poor timing of the meltwater release, as the daily sums of modelled snowpack runoff do not show large differences (see Fig. 3).At CDP, the RE scheme achieves a similar NSE coefficient for the 1 h time scale as at WFJ (respectively 0.47 and 0.44 for the Yamaguchi and Daanen parameterization).On the other hand, the bucket scheme achieves a much better NSE coefficient than at WFJ (0.39).Although an exact reason for the contrasting results was not found, we think it is a confirmation of the fact that once the snowpack is isothermal, differences in water balance schemes are rather small.These snowpack conditions are more frequent at CDP than at WFJ.Moreover, the shallower and less stratified snowpack at CDP likely reduces the differences between the water transport as calculated by the different schemes.
The constructed runoff series from the average LWC production shows a strong decrease in NSE coefficient on the smaller time scales, indicating the importance of taking into account travel time and intermediate storage in the snow cover for sub-daily time scales.
Figure 6 shows that the NSE coefficients for hourly sums of snowpack runoff exhibit variation from year to year.For bucket and NIED simulations, NSE coefficients for hourly snowpack runoff are close to zero or even negative for WFJ, indicating that the model has poor performance on the hourly time scale.RE shows a much better agreement with measurements than the other two schemes for most years, with Yamaguchi's parameterization being the best.At CDP, RE has overall the best performance, although the spread between the RE and bucket scheme is much smaller than at WFJ. Year-to-year variability at CDP seems to be smaller than at WFJ, maybe caused by a smaller year-to-year variation in internal snowpack structure.

Timing of hourly runoff
The NSE coefficients for WFJ simulations on the hourly time scale were close to 0 or even negative for the bucket and NIED schemes, which was ascribed to poor timing of the meltwater release in the bucket and NIED simulations.To quantify the timing of snowpack runoff, lag correlations were calculated between each of the simulations and the measured snowpack runoff.The time span was limited to −12 and +12 h, to prevent correlations with the daily cycle.Table 2 shows the time lag belonging to the highest lag correlation, with negative values indicating the snowpack runoff is too early in the day in the simulations.For WFJ, the bucket and NIED schemes have about 1-2 h too early snowpack runoff compared to measured snowpack runoff, while both simulations with RE show fairly good agreement in timing.Yamaguchi's parameterization seems to provide the best agreement with the measurements.The time lag between the average LWC production and measured snowpack runoff is about 1-3.5 h, showing the importance of the travel time through the snow cover for the sub-daily time scale.At CDP, the bucket and NIED schemes provides the best timing of the peak in meltwater release, whereas the Daanen's parameterization is worst.As can be seen by the low time lag value for the average LWC production, the snowpack at CDP is reacting much quicker to melt than at WFJ, probably due to a shallower and less stratified snowpack.Furthermore, when the modelled density of the snowpack is too high, as suggested by the underestimation of snow height at CDP (see Table 1), hydraulic conductivity is estimated too low (Eq.14).This will result in an overestimation of travel time through the snowpack in simulations with RE.
The simulations with RE for WFJ did not show a significant time lag between lysimeter measurements and modelled snowpack runoff.The existence of preferential flow paths in snow would give the expectation that a time lag would exist when modelling matrix flow alone.Only the simulations for CDP show a tendency towards a positive time lag when using RE.Three issues may play a role here.First, the strategy for choosing θ r < θ results in a direct participation of all liquid water in water transport.This may be a simplification of reality that would compensate for the error of neglecting preferential flow paths in the model.Second, it is possible that experiments to derive parameterizations for the water retention curve and hydraulic conductivity already incorporate preferential flow effects due to the measurement setup in which average flow behaviour is observed.Finally, the simulation results can be interpreted as that the existence of preferential flow paths is not essential for correctly modelling snowpack runoff, because the amount of water involved in preferential flow is small.However, this is not supported by results from other studies (e.g.Marsh, 1999Marsh, , 2006)).

Relation modelled and measured runoff
Figures 7 and 8 show scatter density plots for both the bucket scheme and RE (Yamaguchi) model for both the 24 h and 1 h time scale.In the figures for WFJ, the years 1997, 2000 and 2005 are left out, as the earlier discussed high discrepancy between modelled and measured snowpack runoff for these years appeared to be identifiable as outliers.The figures show the relative distribution of combinations of measured and modelled snowpack runoff for the 11 remaining snow seasons and for the 17 snow seasons at CDP.For the 24 h time scale, the scatter density plots for both water balance schemes look very similar, in contrast to the 1 h time scale.Combined with an almost equal r 2 value, this confirms the conclusion that all water balance schemes have almost equal performance on the daily time scale.This conclusion can be drawn for both sites.
The distribution for the 1 h time scale shows that for WFJ, modelled snowpack runoff with the RE (Yamaguchi) scheme agrees better with measured snowpack runoff than the bucket scheme.This is also expressed by the much higher r 2 value.In contrast with WFJ, the results for the 1 h time scale in terms of r 2 are almost similar for CDP for both schemes.Besides the main distribution along the diagonal where modelled snowpack runoff equals measured snowpack runoff, two distinct features are found.First, the bucket scheme seems to produce considerable amounts of snowpack runoff (−2 to −5 mm h −1 ) when there is almost no snowpack runoff measured.This effect was not present on the 24 h time scale, so it likely originates from the fact that the bucket scheme is releasing meltwater too early on the hourly time scale.Here, the neglect of travel time through the snowpack in the bucket scheme plays an important role.For WFJ, this is supported by the lag correlation (see is also visible at CDP, where no time lag was found.This apparent contradiction originates from the difference between both statistical methods.Second, another small effect is that especially in the RE (Yamaguchi) scheme, there is snowpack runoff observed on the order of −2 to −4 mm h −1 , where at the same time the modelled snowpack runoff is close to zero.Interestingly, these two features are visible in simulations for both sites, although it does not result in the same difference in r 2 value between both water balance schemes.

Conclusions
A comparison of measured snowpack runoff by a lysimeter and 3 water balance schemes for the physically-based snowpack model SNOWPACK has shown that simulating water flow through a snow cover using RE achieves the best agreement with an acceptable increase in computation effort.The water balance schemes were tested for two alpine sites with a seasonal snow cover in a different climatological regime.For WFJ, NSE coefficients for simulations with RE were higher on both the daily and sub-daily time scales when compared to bucket and NIED simulations.The strongest improvement is on the sub-daily time scales.This is also supported by r 2 values between measured and modelled runoff.At CDP, the improvement in NSE coefficient when using RE compared to the bucket scheme is modest, and absent in the r 2 value.At both sites, NSE coefficients vary from year-to-year, quite synchronously: years with lower NSE coefficients have low NSE coefficients in all water balance schemes and vice versa.This indicates that measurements of either snowpack runoff or meteorological forcing have systematic errors that also vary from year to year.On the other hand, the specific implementation of modelled processes will also introduce errors that may have varying effects from year-to-year, depending on the snowpack structure.The timing of meltwater arrival at the base of the snowpack and the runoff dynamics throughout the season also improved with RE, for both the daily and the sub-daily time scale.On the seasonal time scale, bucket and NIED simulations seem to retain the meltwater in the snowpack too long, underestimating the arrival of meltwater at the base of the snowpack in the early stages of the melt season.This was particularly found at WFJ, where there is a thick snow cover with strong layering and a clearly distinguishable dry, cool snowpack phase and a melting snowpack phase.At CDP, where melt is occurring frequently, it was found that especially in the middle of the winter, more snowpack runoff is measured than modelled.This is most pronounced in the bucket scheme.This may be caused by preferential flow paths, that provide efficient water transport to the bottom of the snowpack and which are not modelled.The results from both sites confirm that after the snow cover has become isothermal, the snowpack runoff is mainly determined by the meltwater production near the surface due to a positive energy balance and all water balance schemes route this water to snowpack runoff within the day, causing similar performance on the daily time scale.
The Cryosphere, 8, 257-274, 2014 For the sub-daily time scale, lag correlation coefficients showed that bucket and NIED simulations release meltwater too early in the day compared to measured snowpack runoff at WFJ.This is mainly because the water balance schemes in these two simulations either do not incorporate or else underestimate travel time of liquid water through the snowpack.This result is in contrast with the observations on the seasonal time scale, where it was found that these schemes retain meltwater too long.Apparently, the bucket and NIED scheme require too much time to bring the snowpack to isothermal conditions, by propagating melt water downward too slowly.But once the snowpack is isothermal and wet in these schemes, additional meltwater produced during daytime is propagated downwards too quickly.
For RE, the timing of both seasonal and daily snowpack runoff is in better agreement with measured snowpack runoff.At CDP, the bucket and NIED schemes showed a better agreement of the estimation of the peak water outflow.Of the two parameterizations of the water retention curve in the van Genuchten model, Yamaguchi's parameterization has on average the best representation of travel time, whereas Daanen's parameterization causes a little too slow travel time.
It was generally found that using RE improved modelled snowpack runoff estimations more for the high alpine site WFJ than for the lower site CDP.We propose that the shallower snowpack at CDP reduces the differences arising from using different water balance schemes.Furthermore, the frequent melt cycles occurring at CDP are homogenizing the snowpack, reducing the stratification.The absence of strong inhomogeneities in the snowpack reduces the advantage the RE would have of explicitly taking differences in capillary suction between layers into account.
This study has shown that solving RE for snow is improving several aspects of modelled snowpack runoff considerably, especially for deep, sub-freezing snow covers, which mainly form at high altitude.Yamaguchi's parameterization shows the best overall performance, especially in terms of NSE coefficients.In the simulations for WFJ, the LWC distribution in the snowpack seems to cause slower settling in SNOWPACK, overestimating snow heights when using Daanen's parameterization.This will be a drawback in some applications, especially if snow depth is used to calculate winter precipitation (Lehning et al., 2002a).Note that this study did neither consider the internal snowpack microstructure nor the LWC or density distribution, but focussed only on snowpack runoff.The use of RE may have a considerable effect on these internal snowpack properties.

Outlook
As was pointed out by Marsh (1999), there is a strong need for a better understanding and improved models to describe the complex water flow in natural snow covers.This study has shown that when solely looking at snowpack runoff, improving water balance schemes has an important consequence for the accuracy of snowpack runoff modelling.We also found that the available van Genuchten parameterizations have different effects on modelled snowpack runoff, although they are both based on experiments.Because the number of experimental studies analysing liquid water flow in snow is limited, the results in this study suggest that more experiments for different snow types are welcome as there is potential for the improvement of the model performance of water flow in snow.New measurement techniques are available (Walter et al., 2013), which allow the detailed investigation of water flows including preferential flows in snow.This will help to solve the open question on the importance of preferential versus matrix flow, which has been raised in this study.Moreover, the improvement in performance regarding snowpack runoff suggests that a deeper analysis on the effects on the internal microstructure of the snow, such as the formation of melt-freeze crusts is needed.

Appendix A
Solving RE for snow involves some numerical challenges.Generally, many layers form in deep seasonal snowpacks such as at Weissfluhjoch, which can cause strong inhomogeneities in grain size and density.For numerical models, these inhomogeneities may impact numerical performance.Here, the specific numerical implementation will be discussed with some best practice methodology.

A1 Time step control
To be able to simulate a complete snow season with optimal numerical performance, it is unavoidable to use variable time steps, as the dry snowpack in the winter months can be treated with much larger time steps than the spring snowmelt or ROS events.Infiltration fronts of meltwater in dry snow or soil require small time steps, because Picard iteration is known to have slow convergence for these situations (Paniconi and Putti, 1994).Therefore, we divide the SNOWPACK time step of 15 min in smaller time steps with variable length for solving RE, by applying the time stepping approach proposed by Paniconi and Putti (1994).Based on the number of iterations needed to achieve convergence in the Picard iteration scheme (N iter ), the new time step t new is based on the previous time step t old according to (determined by trial and error):

A2 Convergence criterion
To determine convergence of the solution, the absolute change in solution between two iterations is required to be below a certain threshold.In the proposed Picard iteration by Celia et al. (1990), convergence is checked by a threshold value for the change in pressure head ( h < h ).However, Huang et al. (1996) have found that for non near-saturated conditions, determining convergence based on a threshold for θ ( θ < θ ) will also work and will generally improve convergence considerably, especially in dry conditions as often found in the snowpack.However, close to saturation and in ponding conditions, determining convergence based on θ is not possible and the pressure head condition should be used.Therefore, it is set that where > 0.99, the pressure head criterion is used and the θ criterion elsewhere.The θ criterion at low saturation can cause very inaccurate pressure head estimations and consequently large errors in the flux determined by Eq. ( 15).To achieve an accurate estimation of the snow-soil interface flux, convergence in the upper-most soil layer and lowest snow layer is always judged by the pressure head criterion.Values for h and θ are set to 1 • 10 −3 m and 1 • 10 −5 m 3 m −3 respectively, based on Huang et al. (1996).

A3 Treatment of dry and refreezing snow layers
Fresh snow is dry below freezing and refreezing can cause the formation of dry layers.This gives a singularity in the Van Genuchten model for θ ≤ θ r , associated with an infinitely low pressure head.This problem is circumvented by initialising these new dry snow layers with a very low pressure head.First, in case θ < θ r + θ 10 , θ r was determined by setting it equal to: Then, the following algorithm was used to determine the initial pressure head for dry snow layers: for each layer, the pressure head associated with θ = θ r + θ 10 was determined, so not detectable by the convergence criterion.The smallest value found for pressure head was chosen to initialise dry snow layers with.The associated tiny amount of LWC was created by melting the ice matrix.To prevent a continuously refreezing and subsequent melting of these tiny amounts of water in the snow cover, refreezing of meltwater in snow was only allowed for LWC exceeding 0.01 %.This value is small enough not to influence other snowpack calculations (e.g.wet snow metamorphism).The initialisation value of the pressure head was also used to determine the maximum allowed evaporative flux at the top of the domain, by prescribing it for an imaginary grid point outside the model domain.

A4 Hydraulic conductivity at the interface between nodes
The finite difference approach to solve RE uses the center point of snowpack layers as nodes.The scheme therefore requires an estimation of the hydraulic conductivity at the interface between two layers.In literature, several methods to approximate hydraulic conductivity at the interface nodes between layers have been proposed (e.g.harmonic averaging, geometric averaging, see for example Szymkiewicz and Helmig, 2011).Several approaches were tested, but many yield very bad numerical performance and many tests would not complete within reasonable time.Main problems arise at interfaces where the hydraulic conductivity varies over several orders of magnitude (dry fresh snow layers on top of old snowpack).Most proposed averaging methods tend to put more weight on the lowest value for K. Arithmetic mean was found to work best, as it will effectively smooth large gradients in hydraulic conductivity.We suggest that the limited knowledge of hydraulic properties of snow is likely to influence calculation results stronger than errors arising from inaccurate approximations of the hydraulic conductivity at the interface nodes.

Fig. 2 .
Fig. 2. Cumulative snowpack runoff at WFJ for the lysimeter measurements and the various water balance schemes for the 1998 melt season (a) and hourly snowpack runoff for the lysimeter measurements and the various water balance schemes for one week during the 1998 melt season (b).

Fig. 3 .
Fig. 3. NSE coefficients (model vs. measured) for the snow seasons for the 24 h time scale for WFJ (a) and CDP (b).The NSE coefficients for snow season 2000 for WFJ are negative and not shown.

Fig. 4 .
Fig. 4. Average delay over all simulated years between modelled and observed snowpack runoff as a function of percentage of total melt season snowpack runoff (for WFJ (a), starting 1 March) or total snow season snowpack runoff (for CDP, b).

Fig. 5 .
Fig. 5. NSE coefficients (model vs. measured) over all snow seasons for the studied water balance schemes for different sub-daily time scales for WFJ (a) and CDP (b).

Fig. 6 .
Fig. 6.NSE coefficients (model vs. measured) for the snow seasons for the 1 h time scale for WFJ (a) and CDP (b).

Fig. 8 .
Fig. 8. Density scatter plots for modelled versus measured snowpack runoff (mm) for the 1 h time scale for the bucket scheme (a, c) and RE (Yamaguchi) model (b, d) for both WFJ (a, b) and CDP(c, d).In colour are shown the percentages in the specific region, cut off at 1 %.Data is divided in 15 × 15 bins.For an improved distinction between WFJ and CDP results, the colour scale for CDP is inverted.

Table 1 .
Duration a , maximum measured snow height and seasonal average difference in measured and modelled snow height b for all snow seasons.
a The snow season duration is formatted in DD-MM and denotes the period for which a snow cover existed in either the measurements or one of the model variants.b A positive value denotes an overestimation of snow height in the model compared to measured snow height (and vice versa).

Table 2 .
Lag correlation coefficients (hours) for modelled snowpack runoff compared to measured snowpack runoff.Negative values mean the measured snowpack runoff should be shifted back in time to match modelled snowpack runoff.

www.the-cryosphere.net/8/257/2014/ The Cryosphere, 8, 257-274, 2014
In colour are shown the percentages in the specific region, cut off at 1 %.Data is divided in 15 × 15 bins.For an improved distinction between WFJ and CDP results, the colour scale for CDP is inverted.

www.the-cryosphere.net/8/257/2014/ The Cryosphere, 8, 257-274, 2014 N. Wever et al.: Solving Richards Equation for snow covers and
25 t old , N iter ≤ 5 t old , 5 < N iter ≤ 10 0.5 t old , 10 < N iter ≤ 15 back − step, N iter > 15 (A1)When a back-step is performed, the time step is repeated with a smaller time step: t new = f bs t old , with f bs = (1/3) n sb n sb being the number of sequential back-steps (for 1st backstep: n sb = 1).Back steps are not only performed when N iter > 15, but also when (i) change in pressure head between iterations exceeds a prescribed value (10 3 m) or (ii) the mass balance is violated (mass balance error > 0.1 kg m −12 ).Both are early signs of growing numerical instabilities due to a too large time step.Doing a backstep immediately saves computation time by not executing all 15 iterations.