Verification of the multi-layer SNOWPACK model with different water transport schemes

The widely used detailed SNOWPACK model has undergone constant development over the years. A notable recent extension is the introduction of a Richards equation (RE) solver as an alternative for the bucket-type approach for describing water transport in the snow and soil layers. In addition, continuous updates of snow settling and new snow density parameterizations have changed model behavior. This study presents a detailed evaluation of model performance against a comprehensive multiyear data set from Weissfluhjoch near Davos, Switzerland. The data set is collected by automatic meteorological and snowpack measurements and manual snow profiles. During the main winter season, snow height (RMSE: < 4.2 cm), snow water equivalent (SWE, RMSE: < 40 mm w.e.), snow temperature distributions (typical deviation with measurements: < 1.0 C) and snow density (typical deviation with observations: < 50 kg m) as well as their temporal evolution are well simulated in the model and the influence of the two water transport schemes is small. The RE approach reproduces internal differences over capillary barriers but fails to predict enough grain growth since the growth routines have been calibrated using the bucket scheme in the original SNOWPACK model. However, the agreement in both density and grain size is sufficient to parameterize the hydraulic properties successfully. In the melt season, a pronounced underestimation of typically 200 mm w.e. in SWE is found. The discrepancies between the simulations and the field data are generally larger than the differences between the two water transport schemes. Nevertheless, the detailed comparison of the internal snowpack structure shows that the timing of internal temperature and water dynamics is adequately and better represented with the new RE approach when compared to the conventional bucket scheme. On the contrary, the progress of the meltwater front in the snowpack as detected by radar and the temporal evolution of the vertical distribution of melt forms in manually observed snow profiles do not support this conclusion. This discrepancy suggests that the implementation of RE partly mimics preferential flow effects.


Introduction
One-dimensional multi-layer physics-based snowpack models, for example SNTHERM89 (Jordan, 1991), CRO-CUS (Brun et al., 1989;Vionnet et al., 2012) and SNOW-PACK (Lehning et al., 2002a, b), are widely used to assess various aspects of the snow cover.Recently, the SNOW-PACK model has been extended with a solver for Richards equation (RE) in the snowpack and soil, which improved the simulation of liquid water flow in snow from the perspective of snowpack runoff compared to a conventional bucket-type approach (Wever et al., 2014).In that study, a comparison of snowpack runoff measured by a snow lysimeter with modeled snowpack runoff showed a higher agreement when simulating liquid water flow with RE, especially on the sub-daily timescale.Additionally, the arrival of meltwater at the base of the snowpack in spring was found to be better predicted.However, these results were solely based on an analysis of liquid water outflow.The study raised questions to what extent the two water transport schemes differ in the simulation of the internal snowpack structure and whether the improvements in snowpack runoff estimations with RE are also consistent with simulations of the internal snowpack.
For many applications, especially in hydrological studies, the primary variables of interest are snow water equivalent (SWE) and snowpack runoff, as the first provides possible future meltwater and the latter provides the liquid water that directly participates in hydrological processes.In spite of its importance, direct measurements of SWE are relatively sparse.In contrast, snow height measurements are relatively easy to obtain either manually or automatically, and long climatological records of snow height are available.Methods have been developed to relate snow height to SWE (Jonas et al., 2009;Sturm et al., 2010).Snow density is another parameter that is variable in time and space (Bormann et al., 2013) and rather cumbersome to measure in the field.Although it is seldom of primary interest, it may serve wide applications as an intermediate parameter between a property that is observed and a property that one is interested in.For example, proper estimates of snow density will increase the accuracy of translating snow height to SWE.Snow density is also required for the conversion of measured two-way travel time (TWT) from radar applications to snow depth in dry-snow conditions (Gubler and Hiller, 1984;Lundberg and Thunehed, 2000;Marshall et al., 2007;Heilig et al., 2009Heilig et al., , 2010;;Okorn et al., 2014) or translating dielectric measurements to liquid water content, as for example with the Snow Fork (Sihvola and Tiuri, 1996) or the Denoth meter (Denoth, 1994).
Apart from bulk snowpack properties, there is also a demand for detailed snowpack models to assess the layering and microstructural properties of the snowpack, for example with the purpose of avalanche forecasting.Layer transitions within the snow cover with pronounced contrasts in for example density, grain shape or grain size can act as zones in which fractures can be initialized and slab avalanches release (Schweizer et al., 2003).The presence of liquid water can reduce the strength of a snowpack considerably (Colbeck, 1982;Conway and Raymond, 1993).Techel et al. (2011) showed that this reduction of strength depends also on the grain shape in the snow layers.When snowpack models are used to understand wet snow avalanche formation, it is important that the model can reproduce capillary barriers, at which liquid water may pond (Schneebeli, 2004;Baggi and Schweizer, 2009;Hirashima et al., 2010;Mitterer et al., 2011b).Also the arrival of meltwater at the bottom of the snowpack is considered to be a good indicator for the onset of wet snow avalanche activity.However, reliable liquid water content (LWC) measurements for the snowpack are difficult to obtain.Some attempts for continuous monitoring are promising (Schmid et al., 2014;Koch et al., 2014;Avanzi et al., 2014) but are not yet operational.Recently, Mitterer et al. (2011a) and Schmid et al. (2014) demonstrated the potential of upward-looking ground-penetrating radar (upGPR) to monitor the progress of the meltwater front and Heilig et al. (2015) present data for quasi-continuous observations of bulk liquid water content over several years and for three different test sites.Here, their results concerning the position of the meltwater front will be compared with snowpack simulations.We also consider temperature measurements taken during manual snow profiling as a reliable and precise way to determine which part of the snowpack is at the melting point (often termed isothermal) and likely contains a fraction of liquid water due to infiltration (i.e., the movement of liquid water in snow) or local snowmelt.
As with snow density, snow temperatures are rarely of primary interest in snow studies.However, a correct representation of the temperature profile of the snowpack is required, as it has a large influence on the snow metamorphism (grain shape and size) and settling rates (Lehning et al., 2002a).Temperature gradients drive moisture transport and have a strong influence on the grain growth (Colbeck, 1982;Pinzer et al., 2012;Domine et al., 2013).Furthermore, temperature profiles are an indicator of whether the combination of the surface energy balance, the ground heat flux and the internal heat conductivity of the snowpack is adequately approximated.
In this study, the SNOWPACK model is driven by measurements from an automated weather station at Weissfluhjoch (WFJ) near Davos, Switzerland.Simulations are extensively verified for several bulk properties of the snowpack and against snow profiles made at WFJ, with the aim to verify the representation of the internal snowpack structure.Time series of soil and snow temperatures, snow lysimeter measurements and upGPR data from WFJ are used to validate snowpack temperature profiles, snowpack runoff and the progress of the meltwater front within the snowpack in the simulations.This study focusses on snowpack variables that are influenced by liquid water flow with the aim of a more in-depth comparison of differences between RE and the conventional bucket scheme.The comparison is limited to snow height, SWE, liquid water runoff from the snow cover, snow density, snow temperature and grain size and shape, as for these variables validation data are available.Internally, the SNOWPACK model also uses additional state variables, like sphericity, dendricity and bond size (Lehning et al., 2002a).

Theory
The theoretical basis of the SNOWPACK model regarding the heat transport equation and snow settling has been discussed in Bartelt and Lehning (2002).The treatment of the snow microstructure and several parameterizations, as for example for snow viscosity, snow metamorphism and thermal conductivity, are presented in Lehning et al. (2002a).Some of those parameterizations have been refined in later versions of SNOWPACK.The treatment of the meteorological forcing for determining the energy balance at the snow surface is discussed in Lehning et al. (2002b).Finally, the liquid water transport schemes are presented and verified in Wever et al. (2014).Here, we will outline theoretical aspects not discussed in the aforementioned literature.

Water retention curves
Richards equation in mixed form reads (Richards, 1931;Celia et al., 1990) where θ is the volumetric liquid water content (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 ).
To solve this equation, the water retention curve and the saturated hydraulic conductivity K sat (m s −1 ) need to be specified.For the water retention curve, the van Genuchten model is used (van Genuchten, 1980): The water retention curve is then described by several parameters: residual water content θ r (m 3 m −3 ), saturated water content θ s (m 3 m −3 ) and parameters α (m −1 ), n (−) and m (−).Sc corrects the water retention curve using the approach by Ippisch et al. (2006) to take into account the air entry pressure.As in Wever et al. (2014), an air entry pressure of 0.0058 m was used, corresponding to a largest pore size of 5 mm.Note that the residual water content in the water retention curve, which is the dry limit, is not comparable to the water holding capacity or irreducible water content in the bucket scheme, which refers to wet conditions.For the soil part, the ROSETTA class average parameters (Schaap et al., 2001) are implemented to provide these parameters for various soil types.
For snow, the parameterization for α in the van Genuchten model as proposed by Yamaguchi et al. (2010) reads where 2r g is the classical grain size (m), which is defined as the average maximum extent of the snow grains (Fierz et al., 2009).For n, the original parameterization by Yamaguchi et al. (2010) was modified by Hirashima et al. (2010) to be able to extend the parameterization beyond grain radii of 2 mm: Here, we will abbreviate this parameterization of the water retention curve as Y2010.This parameterization has been used in Wever et al. (2014).The Y2010 parameterization was determined for snow samples with similar densities.In Yamaguchi et al. (2012), an updated set of experiments was described for a wider range of snow density and grain size, leading to the following parameterization of the van Genuchten parameters: and where ρ is the dry density of the snowpack (kg m −3 ).This parameterization will be referred to as Y2012.Both parameterizations will be compared here.θ r and θ s are defined as described in Wever et al. (2014) and K sat is parameterized following Calonne et al. (2012) 1 : where ρ w and ρ ice are the density of water (1000 kg m −3 ) and ice (917 kg m −3 ), respectively, g is the gravitational acceleration (taken as 9.8 m s −2 ), µ is the dynamic viscosity (taken as 0.001792 kg (m s) −1 ), θ i is the volumetric ice content (m 3 m −3 ) and r es is the equivalent sphere radius (m), approximated by the optical radius, which in turn can be parameterized using grain size, sphericity and dendricity (Vionnet et al., 2012).
In both parameterizations and for soil layers, the van Genuchten parameter m is chosen as such that the Mualem model for the hydraulic conductivity in unsaturated conditions has an analytical solution (van Genuchten, 1980).The method to solve RE requires the calculation of the hydraulic conductivity at the interface nodes.It is common to take the arithmetic mean (denoted AM) of the hydraulic conductivity of the adjacent elements, although other calculation methods have been proposed (e.g., see Szymkiewicz and Helmig, 2011).Here, we compare the default choice of AM with the geometric mean (denoted GM), as proposed by Haverkamp and Vauclin (1979), to investigate the possible influence of the choice on averaging method on the simulations of liquid water flow.

Soil freezing and thawing
Due to the isolating effects of thick snow covers and the generally upward-directed soil heat flux, soil freezing at WFJ is mostly limited to autumn and the beginning of the winter, when the snow cover is still shallow.To solve phase changes in soil, we follow the approach proposed by Dall'Amico et al. (2011).They express the freezing point depression in soil as a function of pressure head as where T * is the melting point of the soil water (K), T melt is the melting temperature of water (273.15K), L is the latent heat associated with the phase transition from ice to water (334 kJ kg −1 ) and h is the pressure head (m).
When the soil temperature T (K) is at or below T * , the soil is in freezing or thawing state and a mixture of ice and liquid water is present.Then, the pressure head associated with the liquid water part (h w , m) can be expressed as where h is the total pressure head of the soil (m).The van Genuchten model provides the relationship between pressure head and LWC: where θ is the volumetric LWC (m 3 m −3 ).Consequently, the ice part can be expressed as In Dall'Amico et al. (2011), a splitting method is introduced to solve both the heat transport equation and RE for liquid water flow in a semi-coupled manner.We approach the problem by finding the steady-state solution for T , θ and θ i in Eqs. ( 10), ( 11) and ( 12).This steady-state solution is found numerically by using the Bisect-Secant method (Dekker, 1969), where the starting points for the method are taken as all ice melting and all liquid water freezing, respectively.In soil, liquid water flow can advect heat when a temperature gradient is present.In the soil module of SNOWPACK, heat advection associated with the liquid water flow is calculated after every time step of the RE solver, before assessing soil freezing and thawing.

Data (1): meteorological time series
The SNOWPACK model is forced with a meteorological data set from the experimental site WFJ at an altitude of 2540 m in the Swiss Alps near Davos (WSL Institute for Snow and Avalanche Research SLF, 2015b).This measurement site is located in an almost flat part of a southeasterly oriented slope.During the winter months, a continuous seasonal snow cover builds up at this altitude.The snow season is defined here as the main consecutive period with a snow cover of at least 5 cm on the ground during the winter months and is denoted by the year in which they end.The snow season at WFJ generally starts in October or November and lasts until June or July.
The data set contains air temperature, relative humidity, wind speed and direction, incoming and outgoing longwave and shortwave radiation, surface temperature, soil temperature at the interface between the snowpack and the soil, snow height and precipitation from a heated rain gauge (Marty and Meister, 2012;Schmucki et al., 2014).An undercatch correction is applied for the measured precipitation (Wever et al., 2014).Snow temperatures are measured at 50, 100 and 150 cm above the ground surface, using vertical rods placed approximately 30 cm apart.From September 2013 onwards, soil temperatures are measured at 50, 30 and 10 cm depth.The experimental site is also equipped with a snow lysimeter with a surface area of 5 m 2 , as described in Wever et al. (2014).The rain gauge and snow lysimeter measure at an interval of 10 min, whereas most other measurements are done at 30 min intervals.
In the area surrounding WFJ, field data to validate soil freezing and thawing are lacking.For modeling the snowpack, the most important influence of the soil is the heat flux that is provided at the lower boundary of the snowpack.For this purpose, we will use the temperature measured at the interface between the soil and the snowpack to validate the soil module.This temperature measurement is influenced by soil freezing and thawing.Our primary interest here is to investigate to what degree the previously described soil module of SNOWPACK is capable of providing a realistic lower boundary for the snowpack in the simulations.
SNOWPACK can be forced with either measured precipitation amounts or with measured snow height.In precipitation-driven simulations (Precip driven), measured precipitation is assumed to be snowfall when the air temperature is below 1.2 • C and rain otherwise.For these types of simulations, the study period is from 1 October 1996 to 1 July 2014 (1 week after melt-out date), consisting of 18 full snow seasons.In case of snow-height-driven simulations (HS driven), an additional threshold for relative humidity (≥ 70 %) and a maximum value for the temperature difference between the air and the snow surface (≤ 3 • C) is used to determine whether snowfall is possible.The latter condition tests for cloudy conditions, when the increase in incoming longwave radiation will warm the snowpack surface close to air temperature.Then, snowfall is assumed to occur when measured snow height exceeds the modeled snow height (Lehning et al., 1999) and, consequently, new snow layers are added to the model domain in order to match The Cryosphere, 9, 2271-2293, 2015 the measured snow height again.These layers are initialized with a new snow density dependent on meteorological conditions (Schmucki et al., 2014).In both modes, new snow layers are added for each 2 cm of new snow.An uninterrupted, consistent data set for this type of simulations is available from 1 October 1999 to 1 July 2014, consisting of 15 full snow seasons.The last snow season (2014) of the studied period has the most data available and will be used as the example snow season to explain how SNOWPACK simulates the snow cover.Results for the other snow seasons are included in the online Supplement.
Many processes in SNOWPACK are based on physical descriptions that require calibration, for example for wet and dry snow settling, thermal conductivity and new snow density.For this purpose, dedicated data sets with some additional detailed snowpack measurements from snow seasons 1993, 1996, 1999 and 2006 have been used when constructing the model.Snow metamorphism processes were mainly calibrated against laboratory experiments (Baunach et al., 2001).

Data (2): manual snow profiles
Every 2 weeks, around the 1st and 15th of each month (depending on weather conditions), a manual full depth snow profile is taken at WFJ (WSL Institute for Snow and Avalanche Research SLF, 2015a), following the guidelines from Fierz et al. (2009).The snow profiling is carried out in the morning hours, starting around 09:00 LT.Measurements include snow temperature at a resolution of 10 cm and snow density in steps of approximately 30 cm.Snow density and SWE are determined by taking snow cores using a 60 cm high aluminium cylinder with a cross-sectional area of 70 cm 2 inserted vertically into the snowpack.The snow core is then weighted using a calibrated spring.For comparison with the simulations, SWE values are corrected for differences in snow height at the snow pit and at the automatic weather station to eliminate the effect of spatial variability.Grain size (following the classical definition of average maximum extent of the snow grains) and grain shape are evaluated by the observer using a magnifying glass.Also snow wetness is reported in five wetness classes as well as hand hardness in six classes (Fierz et al., 2009).Because judging snow wetness has a subjective component and estimating the actual LWC is generally considered rather difficult, we consider here only three categories: dry (class 1; 0 % LWC), moist (class 2; 0-3 % LWC) and wet (class 3 or higher; ≥ 3 % LWC).

Data (3): upward-looking ground-penetrating radar
An upGPR is located within the test site at a distance of approximately 20 m from the meteorological station (Mitterer et al., 2011a;Schmid et al., 2014).The upGPR is buried in the ground with the top edge level to the ground surface and points skyward.The radar instrument and data processing is described in Schmid et al. (2014).Measurement intervals for all observed melt seasons were set to 30 min during daytime.The only difference in the processing scheme applied for this study in comparison to Schmid et al. (2014) is that for an optimized retrieval of the dry-wet transition within the snow cover, we reduced the length of the moving-window time filter to a few days (1-3) instead of 6 weeks.Since percolating water results in strong amplitude increases at the respective depth of percolation and a decrease in wave speed for electromagnetic waves traveling through wet layers, we searched for occurrences of sharp amplitude contrasts together with diurnal variations in the location of signal responses of the overlying layers.For snow layers in which liquid water is appearing during the day and refreezing during the night, or when LWC reduces through outflow, a clear diurnal cycle in TWT of the respective signal reflections can be observed.Schmid et al. (2014) describe first attempts to determine percolation depths automatically within the recorded radargrams.For this study, we manually determined all observations of the dry-wet transition in the snowpack and converted TWT in height above the radar by assuming a constant wave speed in dry snow of 0.23 m ns −1 (Mitterer et al., 2011a;Schmid et al., 2014).Data on liquid water percolation measured with upGPR have been presented in Schmid et al. (2014) for the snow seasons 2011 and 2012.Here, we present data of two more snow seasons (2013,2014) and compare all measured depths of the dry-wet transition with simulation results.In snow seasons 2011, 2013 and 2014, additional snow profiles were made in close proximity of the upGPR, with a higher frequency during the melt season than the regular snow profiles discussed in the previous section.

Methods (1): model setup
For the simulations in this study, the SNOWPACK model solves the energy balance at the snow surface.The turbulent fluxes are calculated using the stability correction functions as in Stössel et al. (2010).This is an adequate approximation for most of the snow season, when the snow surface cooling due to net outgoing longwave radiation causes a stable stratification of the atmospheric boundary layer.The surface albedo is calculated from the ratio of measured incoming and reflected shortwave radiation.The net longwave radiation budget is determined from the difference in measured incoming and calculated outgoing longwave radiation.The aerodynamic roughness length (z 0 ) of the snow is fixed to 0.002 m.The soil at WFJ consists of coarse material with some loam content, as was observed when installing the soil temperature sensors.The ROSETTA class average parameters for the loamy sand class are taken for the van Genuchten parameterization of the water retention curve for the soil (θ r = 0.049 m 3 m −3 , θ s = 0.39 m 3 m −3 , α = 3. ).For the thermodynamic properties, the specific heat for the soil constituents was set to 1.0 kJ kg −1 K −1 and the heat conductivity to 0.9 W m −1 K −1 .The total soil depth in the model is taken as 3 m, with a variable layer spacing of 1 cm in the top layers and 40 cm for the lowest layer.The dense layer spacing in the top of the soil is necessary to describe the large gradients in soil moisture and temperature occurring here.At the lower boundary, a water table is prescribed, together with a Neumann boundary condition for the heat transport equation, simulating a constant geothermal heat flow of 0.06 W m −2 .All simulations are run on the same desktop computer as a single-core process, using a model time step of 15 min.In the solver for RE, the SNOWPACK time step may be subdivided in smaller time steps when slow convergence is encountered (Wever et al., 2014).The computation time is in the order of a few minutes per year, where RE takes about twice as much time as the bucket scheme (Wever et al., 2014).Checks of the overall mass and energy balance reveal that the mass balance for all simulations is satisfied well within 1 mm w.e. and the energy balance error is generally around 0.05 W m −2 (see Table 1).We consider these errors to be well acceptable for our purpose.

Methods (2): analysis
The analysis of the simulations is done per snow season, ignoring summer snowfalls.The snow season at WFJ is characterized by an early phase at the end of autumn or beginning of winter, when the snow cover is still relatively shallow and occasionally melt or rain-on-snow events are occurring.End of November to mid-March can be defined as the accumulation period, in which snowpack runoff is virtually absent and the snowpack temperature is below freezing.This implies that in this period, all precipitation is added to the snow cover as solid mass either by rain refreezing inside the snowpack or by snowfall.Small amounts of snowmelt occurring near the surface refreeze during night or, after infiltration, inside the snowpack.Therefore, the increase in SWE between the biweekly profiles can be used to verify the undercatch correction in case the SNOWPACK model is driven with measured precipitation from the heated rain gauge or to verify the combined effect of parameterized new snow density and snow settling in case snow height is used to derive snowfall amounts.The final phase is the melting phase, starting in April in most snow seasons, when the snowpack is isothermal and wet and produces snowpack runoff.
The snow temperature sensors may be influenced by penetrating shortwave radiation in the snowpack.Therefore, snow temperature measurements are only analyzed when the measured snow height is at least 20 cm above the height of the sensor.Comparing snow temperatures between snow seasons was done by first standardizing the measurement time of the temperature series between 0 and 1 for the start and end of the snow season, respectively.Then the data were binned in steps of 0.01 and bin averages were calculated.These series were then used for calculating the average and SD of differences between snow seasons.The same procedure was followed for snow height.
To compare manual snow profiles with the model simulations, several processing steps are required (Lehning et al., 2001).The snow height at the snow pit is generally different from the simulated snow height.This is not only due to the model not depicting the snowpack development perfectly but also because the snow pit is made at some distance from the snow height sensor which is used to drive the simulations.Therefore, we scale the simulated profile to the observed profile by adjusting each layer thickness, without adjusting the density.This implies that mass may be added or removed from the modeled domain.Then, the model layers are aggregated to match the number and thickness of the layers in the observations.Model layers are assigned to observed layers based on the center height of the model layer.The typical thickness of a model layer is around 2 cm, so possible round-off errors are expected to be small.For temperature, the matching with modeled layer temperatures is achieved by linear interpolation from the measured temperature profile to the center point of the modeled layer.
The cold content of the snowpack is the amount of energy necessary to bring the snowpack to 0 • C, after which an additional energy surplus will result in net snowmelt.The total cold content Q cc (J m −2 ) of the snowpack is defined in discrete form as the sum of the cold content of each layer: where i is an index to a snow layer, n is the number of snow layers in the domain, ρ i is the density of the layer (kg m −3 ), c i is the specific heat of the layer (J kg −1 K −1 ), z i is the layer thickness (m) and T i is the temperature of the layer (K).The cold content is calculated for both the observed and modeled profiles, where the modeled profile is first aggregated onto the observed layer spacing with the procedure described above.

Snow height and snow water equivalent
Figure 1 shows the snow height for several simulation setups.Per construction, the snow-height-driven simulations provide a high degree of agreement between measured and modeled snow height.The general tendency of the precipitationdriven simulations is to follow the measured snow height, although it can be clearly seen that some precipitation events are overestimated, whereas others are underestimated.These differences are caused by inaccuracies when measuring solid precipitation with a rain gauge (Goodison et al., 1998), imperfections in the undercatch correction or the effect of ae-Table 1.Average and standard deviation (in brackets) of bulk snowpack statistics over all snow seasons for various simulation setups (bucket or Richards equation (RE) water transport scheme, snow-height-driven (HS) or precipitation-driven (Precip) simulations, Y2010 (Yamaguchi et al., 2010) or Y2012 (Yamaguchi et al., 2012) water retention curves and arithmetic or geometric mean for hydraulic conductivity) for all simulated snow seasons.Differences are calculated as modeled value minus measured value; ratios are calculated as modeled value divided by measured value.The isothermal part is only considered during the melt phase (from March to the end of the snow season).(Yamaguchi et al., 2010) or Y2012 (Yamaguchi et al., 2012) water retention curves, and arithmetic (AM) or geometric mean (GM) for hydraulic conductivity) for the example snow season 2014, from October 2013 to July 2014.Note that apart from forcing with either snow height or precipitation measurements, differences between simulation setups cause only small differences in snow height simulations, resulting in overlapping lines in the figure .olian wind transport causing either erosion or accumulation of snow at the measurement site.As drifting snow mainly occurs close to the surface, the rain gauge is rather insensitive to these effects as its installation height is higher than the typical depth of a saltation layer.However, at WFJ, drifting snow is expected to play a relatively small role.
As listed in Table 1, the RMSE of snow height for all simulated snow seasons is significantly larger for precipitationdriven simulations than for snow-height-driven ones.As snow-height-driven simulations are forced to closely follow the measured snow height, it can compensate for deviations in measured and modeled snow height due to over-or underestimated snow settling or snowmelt and occasional erosion or deposition of snow by wind.This is not possible with precipitation-driven simulations, which solely take precipitation amounts to determine snowfall.This contrast is additionally illustrated in Figs.S1 and S2 in the Supplement, where snow height for the various model setups is shown for each snow season.Typical year-to-year variability of inconsistencies in the precipitation-driven simulations are present, whereas the snow-height-driven simulations follow the measured snow height more closely.Consequently, the snowheight-driven simulations exhibit a better agreement on the melt-out date, typically within 1 day from the observed meltout date, than the precipitation-driven ones (see Table 1).
In Fig. 2, the average snow height difference is shown for all simulated snow seasons, relative to the standardized date in the snow season.Snow-height-driven simulations generally have almost no bias to measured snow height for most of the snow season.A slight positive bias in mid-winter for precipitation-driven simulations is caused by a few overesti-  2012) water retention curve and arithmetic mean for hydraulic conductivity (RE-Y2012AM).For every snow season, the first day with a snow cover is set at 0 and the last day at 1. mated snowfall events, for which the bias persists throughout the snow season (see for example snow season 2011-2012 in Fig. S2e in the Supplement).Contrastingly, in the end of the snow season (i.e., the melt season), an underestimation of the snow height occurs in precipitation-driven simulations, which is also expressed by a negative overall snow height bias in Table 1.This does not necessarily imply that the melt rates are overestimated, as snow height is the combined result of snow accumulation, settling and melt.
SWE is generally a better indicator of snow accumulation and snowmelt than snow height.A comparison between observed SWE in manual profiles and modeled SWE (Fig. 3a) shows that the agreement between both is high.The linear fits to the data points show that on average, the prediction of SWE in the model is accurate for both snow-height-driven and precipitation-driven simulations.The scatter is larger for precipitation-driven simulations and there seems to be an underestimation of low SWE values and an overestimation of high ones.
The modeled SWE is a result of several effects: (i) snowfall amounts, which rely on an accurate estimation of new snow density in case of snow-height-driven simulations or an adequate undercatch correction in the case of precipitationdriven simulations; (ii) snow settling in the case of snowheight-driven simulations; (iii) snowmelt; and (iv) liquid water flow in snow and subsequent snowpack runoff.To separate the effects of liquid water flow and snowpack runoff from the other effects, Fig. 3b shows the increase in SWE in biweekly profiles during the accumulation phase of the snow season at WFJ, when only factors (i), (ii) and (iii) play a role.The snow-height-driven simulations on average provide a high degree of agreement with the measured increase in SWE during the accumulation phase, with only a marginal difference between the bucket scheme and RE.Here, it needs to be mentioned that in snow-height-driven simulations, the snow settling formulation is able to compensate for errors in the estimation of new snow density and vice versa.For example, when new snow is initialized with a too high density, and thus too much mass is added, the snow settling will be underestimated, and consequently, the next snowfall amount is also underestimated.Because the snowfall amounts in precipitation-driven simulations are independent of the settling of the snow cover, the increases in SWE are independent of the predicted settling.From the linear least squared fit to the observed and simulated changes in SWE, it can be concluded that in the accumulation phase, the combined effect of new snow density and snow settling provides a slightly underestimated SWE increase in snow-height-driven simulations, whereas the opposite is found for precipitation-driven simulations.In the latter case, particularly a few overestimated large snowfall events can be identified to have influenced the fit.
Figure 4 shows the difference in SWE between model simulations and the snow profiles for all simulated snow seasons.The difference in snow-height-driven simulations is rather small, compared to precipitation-driven simulations.All simulations show that in the melt phase, the model underestimates SWE.This points towards either an overestimation of melt rates, a too early release of meltwater at the base of the snowpack or a combination of both.The fact that the discrepancies for the precipitation-driven simulations are larger than for the snow-height-driven ones is related to the underestimation of snow height during the melt phase.In the snow-height-driven mode, an overestimated decrease in snow height during snowmelt is compensated for by a continuous adding of fresh snow when the snowfall conditions are met.

Liquid water content and snowpack runoff
Figure 5a and b show the distribution of liquid water within the snowpack for the example snow season 2014 for the bucket scheme and RE, respectively.Here, liquid water is present during the beginning of the snow season and during the melt season, which is a typical pattern for WFJ.The simulations with RE show a quicker downward routing of meltwater from the surface, where the meltwater is produced, than the simulations with the bucket scheme.Furthermore, the latter provides a rather homogeneous LWC distribution throughout the snowpack, except for the lighter surface elements, where LWC is significantly higher.A diurnal cycle is not visible in the simulations, except for layers close to the surface.With RE, there is both a strong variation in the vertical direction as well as in time.Marked accumulations of liquid water can be seen at transitions between layers with different characteristics.These accumulations peak at around 10 % LWC and occur during the first wetting of the snowpack and above capillary barriers inside the snowpack.The appar-  ent slow downward movement of liquid water accumulations during the melt season results from snowpack settling moving the specific layers with water accumulations closer to the ground.

The
The formation of water accumulations on capillary barriers was also observed in natural snow covers (e.g., Techel and Pielmeier, 2011), and this process is considered to contribute to wet snow avalanche formation (Schneebeli, 2004;Baggi and Schweizer, 2009).The effect is particularly present during the first wetting, as later in the melt season, wet snow metamorphism reduces the contrast between microstructural properties, and this is at least qualitatively reproduced by the model.Furthermore, the increase in hydraulic conductivity when the snowpack below the capillary barrier gets wet, reduces its function as a barrier.RE also introduces a strong diurnal cycle in LWC in the simulations.The results for other snow seasons can be found in the Supplement Figs.S3-S5, and they illustrate that the differences occurring between both water transport schemes in the example snow season are similar for the other snow seasons as well.
Direct comparison of these model results with measurements is difficult, as continuous, non-destructive observations of the vertical distribution of LWC are not available.
However, snowpack runoff is strongly coupled to the LWC distribution.Snowpack runoff at the measurement site WFJ typically occurs in the melt season and in some snow seasons during autumn when early snowfalls may be alternated by short melt episodes or rain-on-snow events.This is illustrated by the cumulative runoff curves in the Supplement Figs.S6 and S7.Table 1 shows the ratio of modeled to measured snowpack runoff.Snowpack runoff from precipitationdriven simulations is on average 2 % less than observed, whereas snow-height-driven simulations show about 8-14 % more runoff than is observed.From the snow-height-driven simulations, simulations with RE again have higher runoff sums than the simulations with the bucket scheme.This behavior is found in most simulated snow seasons, as shown by Fig. 6.The overestimation of total runoff in snow-heightdriven simulations is caused by the previously described mechanism where the snow-height-driven simulations add snow layers in spring when the snow height decrease is overestimated.The approach is inadequate during the melt season, as these new snow layers have low densities compared to the rest of the snowpack and snow settling will quickly reduce the modeled snow height again below the measured one.As the wet snow settling is a little stronger when using RE, this effect is slightly larger for those simulations.
A common measure to quantify the agreement between measured and modeled snowpack runoff is the Nash-Sutcliffe model efficiency (NSE; Nash and Sutcliffe, 1970), which is shown in Table 1 and Figs.S8 and S9 in the Supplement for completeness.Further discussion can be found in Wever et al. (2014).NSE coefficients increase for simulations with RE, especially on the 1 h timescale, as well as the r 2 value.The decrease of performance in terms of NSE coefficient, in particular for the bucket scheme, can be mainly attributed to poor timing of meltwater release during the day.For example, the bucket scheme does not take percolation time into account, resulting in rather low NSE coefficients.The NSE coefficients and r 2 values tend to be lower for precipitation-driven simulations than for snow-height-driven ones, especially in the simulations with RE.This likely is a result of a more accurate prediction of percolation time of liquid water through the snowpack in snow-height-driven simulations.This is also indicated by the difference in time lag correlation (see Table 1) between precipitation-driven simulations and snow-height-driven ones.The best timing of snowpack runoff on the hourly timescale is achieved with snow-height-driven simulations with RE.
The NSE coefficients and r 2 values reported here were calculated over the snow-covered period from the simulations.However, this is an arbitrary choice, given the discrepancies in melt-out date from simulations and measurements, particularly for precipitation-driven simulations (see Table 1).When considering both possible definitions for snow-covered period (either determined from simulations or from measurements), differences in NSE coefficients up to 0.16 are found for individual years.This is particularly the case for precipitation-driven simulations, where the predic-tion of melt-out date is less accurate (see Table 1).However, for the average NSE coefficients, the differences are less than 0.02 for both precipitation and snow-height-driven simulations, as the year-to-year differences cancel out.The choice of calculation period has a larger influence on r 2 values, since the late melt season is associated with the highest snowpack runoff and consequently has a large effect on the r 2 values.Nevertheless, the differences between simulation setups within either snow-height-driven simulations or precipitation-driven ones are smaller than the differences between both simulation types.This implies that the same conclusions about simulation setups can be drawn regardless of the choice of calculation period.

Soil temperatures
At WFJ, soil temperatures are available at three depths but only for the last snow season in this study (see Fig. 7a).The simulated soil temperatures are satisfactorily simulated, although the soil never showed temperatures well below 0 • C.This indicates that no significant soil freezing occurred, limiting the usefulness of these data to validate the new soil module.However, it is primarily important for this study that the soil as modeled by SNOWPACK serves as an adequate lower boundary condition for the snowpack simulations.For this purpose, we examine the soil temperature in the topmost soil part at the snow-soil interface, which is available for the snow seasons 2000-2014 (see Fig. 7b).For most of the time when a snow cover is present, the interface temperature at the snow-soil interface is close to 0 • C, except in the beginning of the snow season when the snow cover is still shallow.This is common for deep alpine snowpacks due to the isolating effect of thick snow covers and the generally upward-directed soil heat flux.Figure 7b shows that the simulations capture the variability in early season soil-snow interface temperature to a high degree in most years and that the soil module in SNOWPACK is providing an accurate lower boundary for the snow cover in simulations.

Snow temperatures
Figure 8a and b show the simulated temperature distribution within the snowpack for the example snow season 2014 for the bucket scheme and RE, respectively.The other snow seasons are shown in the Supplement Figs.S10-S12.For each snow season, the snowpack temperature at WFJ is below freezing for an extended period of time and for these periods no noticeable differences are found between simulations with the bucket scheme or RE.As a result of the differences in liquid water flow depicted in Fig. 5a and b, the parts of the snowpack that are isothermal differ significantly.Table 1 shows that the r 2 value between the relative part of the snowpack that is isothermal, as determined from measurements in the observed snow profiles and from the simulated ones, increases from 0.74 to 0.87 when solving liquid water flow with RE.
The temperature distribution of the snowpack is strongly related to the combination of the net energy balance of the snowpack and snow density.The latter influences the snow temperature through the thermal inertia of dense snow layers and through the strong density dependence of thermal conductivity (e.g., Calonne et al., 2011).Errors in either the energy balance or snow density may result in errors in snow temperatures.The cold content of the snowpack may be considered a more robust method to verify the simulated energy balance of the snow cover.Table 1 shows that the RMSE in cold content in the snow-height-driven simulations is larger for the bucket scheme than RE, with a RMSE of around 630 kJ m −2 , which is equivalent to 2 mm w.e.snowmelt.This shows that the estimation of cold content in the simulations is adequate when, for example, estimating the onset of snowmelt and refreezing capacity inside the snowpack.Larger RMSE for precipitation-driven simulations can be associated with the larger discrepancy between measured and modeled snow height.The bias in the cold content is small compared to the RMSE, denoting that the average simulated energy input in the snowpack is accurate compared to its temporal variation.This conclusion is only valid for the period when the snowpack temperature is below freezing, as in the melt season the cold content is by definition 0 kJ m −2 .
Figure 9 shows the measured and modeled snow temperature time series at three heights for the example snow season.The change of snow temperature over the snow season is adequately captured.There is almost no difference between simulations with the bucket scheme and RE except for the timing when the snowpack gets isothermal, associated with the meltwater front moving through the snowpack.For this example snow season, simulations with RE seem to better capture when the snowpack becomes 0 • C, suggesting a better prediction of the movement of the meltwater front through the snowpack.In the Supplement Figs.S13 and S14, results for each snow season are shown.In most snow seasons, simulations with the RE provide a better agreement with measured temperatures in spring than the bucket scheme.However, in some snow seasons (e.g., 2001 and 2011), simulations with RE show an increase in snow temperature before the measured temperature increases, which suggests a simulated progress of the meltwater front that is too fast.
In Fig. 10a and b, the average and SD, respectively, of the difference between modeled and measured temperatures are shown, including snow surface and snow-soil interface temperatures, determined over all 15 snow seasons of the snowheight-driven simulations and plotted as a function of the relative date in the snow season.During the main winter season, the temperatures at 50 and 100 cm height are on average up to 0.5 • C lower in the model than in the measurements, whereas the temperature at 150 cm is on average up to 1. perature at the highest snow temperature sensor is overestimated in the simulations.The contrasting result suggests that the snow layers near the top of the snowpack have a too low density in the simulations, impacting both thermal conductivity and heat capacity of those layers, or the thermal conductivity is underestimated for typical snow densities found close to the surface.These effects provide a stronger isolation of the snowpack, causing heat from inside to escape at a slower rate and allowing the surface to cool more.This offers an explanation why the underestimation of the snow surface temperature particularly occurs at night (not shown).In contrast, errors in diagnosing the snowpack energy balance (i.e., in net shortwave or longwave radiation or in turbulent fluxes) would be expected to influence all temperature sensors in the same direction.
The SD of the difference between modeled and measured temperatures shows an increase with height above the ground.This can be attributed to higher temporal variations in temperature in the upper snowpack due to highly variable surface energy fluxes.The SD for the snow and snow-soil interface temperature typically is less than 1.0 • C and decreases towards the melt season.For the surface temperature, the SD is typically high in the beginning and the end of the snow season.In the beginning of the snow season, lower snow densities, low air temperatures and reduced incoming shortwave radiation allow for a strong radiative cooling of the snow sur-face, which is delicate to simulate correctly and may result in errors in simulated snow temperatures up to 10 • C. In the melt season, discrepancies in the duration the snow surface needs to refreeze at night may contribute to the increase in SD between modeled and measured surface temperatures.
Figure 10a also shows that in the beginning of the melt season, the difference between snow temperatures simulated with RE and measurements is on average smaller than with the bucket scheme at 0, 50 and 100 cm depth.Although this suggests a better timing of the movement of the meltwater front through the snowpack and the associated temperature increase to 0 • C, heat advection through the ice matrix and preferential flow and subsequent refreezing inside the snowpack may also increase the local snowpack temperature to 0 • C. The reason why the results from the temperature series at 150 cm contrast those at 0, 50 and 100 cm depth remains unclear.denser parts, whereas other layers remain less dense because less meltwater is retained.This type of stratification is known to happen, although verification is difficult, because density is sampled at a low spatial resolution in the manual snow profiles.In Fig. 12a, average snow density as observed in the manual profiles is compared with the modeled snow densities for the snow-height-driven simulations for the period 1999-2013.Generally, the seasonal trend in snow density is captured well in the model.Discrepancies between modeled and observed profiles are larger than the differences aris-ing from the different water transport schemes.In general, SNOWPACK overestimates the density near the base of the snow cover, while it underestimates the density of the upper part of the snowpack.The bucket scheme, which was used to calibrate the wet snow settling, keeps higher densities near the surface than RE, which is in closer agreement with observed snow density.These observations are consistent for all simulated snow seasons, as illustrated in Supplement Fig. S18.It supports the argument in the previous section.These over-and underestimations are larger than the differences between water transport schemes.In Fig. 12b For every snow season, the first day with a snow cover is set at 0 and the last day at 1.The statistics are determined over the 15 snow seasons of the snow-height-driven simulations using the bucket scheme or Richards equation using the Yamaguchi et al. (2012) water retention curve and arithmetic mean for hydraulic conductivity (RE-Y2012AM).

Snow density
the average and SD of the difference between simulated and observed density is shown, determined over the 15 snow seasons of the snow-height-driven simulations.Average discrepancies in snow densities are less than 25 kg m −3 , increasing to 50-100 kg m −3 shortly before melt out.The SD of the discrepancies is less than 50 kg m −3 , increasing to 100-150 kg m −3 near the end of the melt season.This illustrates that the new snow density parameterization and the snow settling formulation are able to provide accurate predictions of snow density.During the snowmelt season, the deviations between observed and simulated snow density increase as a result of new snowfall events that are simulated to compensate for the overestimated SWE depletion.
The depletion rate is the result of many interacting processes.First of all, it is strongly coupled to snowmelt, and thus dependent on the surface energy fluxes.Given the high agreement in cold content in the main winter season, errors in diagnosing the surface energy balance due to uncertainties in atmospheric stability and measurement errors in radiation, wind speed or air temperature seem to be small on average.However, a consistent or incidental overestimation of the energy input in the snow cover during the snowmelt period may result in overestimated snowmelt.Once the meltwater leaves the snowpack, the mass associated with it is definitely lost.Additionally, we would argue that an insufficient simulation of the densification during spring, under the influence liquid water flow, may also be important here.A too low snow density will result in a deeper penetration of shortwave radiation, effectively providing heat transport into the snowpack.Furthermore, heat conductivity will be underestimated, with the consequence that the simulated snowpack in spring is too isolated to be able to release heat during night.

Grain size
Grain size plays an important role in liquid water flow, as it has a strong influence on the water retention curves (Eqs.3-6).Figure 13a and b show modeled grain size profiles for the example snow season 2014 for the bucket scheme and RE, respectively.Differences between the schemes are mainly found in the melt season where the bucket scheme produces slightly larger grains.This is associated with the typically higher liquid water content using that scheme (Fig. 5a) compared to RE (Fig. 5b).This results in a stronger wet snow grain growth rate.Figure 13b also illustrates the cause of the liquid water accumulation found near a height of 120 cm in the beginning of April in Fig. 5b.The layer below the ponding water consisted of significantly larger grains and was creating a capillary barrier for the liquid water.In the Supplement Figs.S19-S21, results are shown for each snow season Figure 14a and b show the average and SD of the grain sizes from the manual profiles and the simulations for the snow seasons 2000-2014.Most distinguishable is the steady increase in grain size towards and during the melt season.Both simulations show an increase in grain size towards the end of the snow season, although the average observed grain size is often underestimated.The underestimation of grain size in simulations with RE compared to the bucket scheme is consistent for most snow seasons.It results from generally lower LWC values in the snowpack in simulations with RE and, consequently, lower wet snow grain growth rates.This contributes to a reduced r 2 value for grain size (see Table 1).Most of the variation in grain size that exists before the initial wetting of the snow remains present throughout the snow season in the simulations.However, the vertical variation of grain size typically decreases during the melt season, shown in Fig. 14b.However, opposite trends can be found, mainly caused by snowfalls during the melt season.The simulations tend to provide a decrease of the SD in the melt season and the agreement with the observations varies from year to year.Especially large variations in grain size in the profiles are not captured in the simulations.

Comparison of simulated dry-wet transition with upGPR
Detailed comparisons of radar-determined dry-wet transitions with simulations of the water transport schemes for the snow seasons 2011 through 2014 are presented in Fig. 15.Measured snowpack runoff (by the snow lysimeter) is included in this presentation together with grain shapes observed in snow pits, both of which are indicative of water flow processes in snow.The dry-wet transition is only plotted when the upGPR signal indicated that parts of the snowpack were wet (see Sect. 3.3) or, for the simulations, when the modeled snowpack was partly wet.Due to beam divergence, a preferential flow path that forms in the vicinity above the upGPR could potentially be detected, although generally the upGPR would be particularly sensitive to matrix flow.However, liquid water accumulations above ponding layers are clearly visible in radargrams independent of matrix or preferential flow that formed such accumulations.It is impossible to discriminate from the radar data which flow regime caused the respective liquid water accumulations.In addition, layer transitions within the resolution limit of the radar (≈ 0.07 m for dry-snow conditions; Schmid et al., 2014)  ble to discriminate as well and as a consequence, percolation depths of the wetting front close to the ground surface (< 10 cm above the ground) cannot be accurately allocated anymore.Interferences with the reflection signal from the cover box of the radar prevent an accurate location of such signals.
From the four snow seasons presented in Fig. 15, the following observations can be made: (i) snowpack runoff measured by the snow lysimeter consistently starts earliest in the snow season; (ii) the progress of the meltwater front is always faster in the simulations with RE compared to the bucket scheme; (iii) the radar-derived meltwater front progresses generally slower through the snowpack than in both water transport schemes in the model; (iv) the manual snow profiles mostly show melt forms in parts of the snowpack that have been wet according to the radar data, whereas the simulations often show larger parts of the snowpack becoming wet earlier than indicated by the profiles.These observations will now be discussed in more detail.
(i) Since preferential flow can route liquid water efficiently through the snowpack (Kattelmann, 1985;Waldner et al., 2004;Techel and Pielmeier, 2011), upGPR-determined depths of dry-wet transitions are not necessarily linked to the onset of measured snowpack runoff (Heilig et al., 2015).Studies by Katsushima et al. (2013) and Hirashima et al. (2014) found that ponding plays a crucial role in forming preferential flow in both laboratory experiments as well as model simulations.The ponding of liquid water in the simulations for WFJ (see Fig. 5b) suggests that preferential flow may have developed.The amount of snowpack runoff measured before the arrival of the meltwater front is highly variable.From 1 to 8 April in snow season 2011, large amounts of snowpack runoff were observed, most likely due to lateral flow processes, whereas in snow season 2014 only marginal amounts were observed.In the latter snow season, there is a strong increase in observed snowpack runoff close to the time of the arrival of the radar-derived meltwater front at the snowpack base.This variability between years is not necessarily caused by different preferential flow path structures but may also result from the limited capturing area of the snow lysimeter (Kattelmann, 2000).(iii, iv) The vertical distribution of the melt forms in the observed snow profiles may be considered particularly representative for matrix flow and for the 4 presented years it generally corresponds well with the parts of the snowpack that may be considered wet from the upGPR signal.(ii) As the bucket scheme shows a higher correspondence with the upGPR data than RE, the convenient improvement in the accuracy of simulated snowpack runoff with RE, as found in Wever et al. (2014)  or by the specifics of the implementation of RE in SNOW-PACK, remains unclear.(ii, iii) Although the bucket scheme may seem to better coincide with the meltwater front in the upGPR data, it may as well be argued that the differences between both water transport schemes are smaller than the discrepancies with the upGPR data.It is likely that the limits of one-dimensional models with a single water transport mechanism will prevent a correct simulation of both snowpack runoff as well as the internal snowpack structure at the same time.
In the beginning of the melt season, observations contrasting to the main melt phase discussed above can be made.The initial melt phase is characterized by a regularly disappearing meltwater front at night.During this period, the depth to which the liquid water infiltrates the snowpack is underestimated in the simulations.Here, the RE scheme shows larger infiltration depths, which are in better agreement with the upGPR data, although again differences between both simulations are smaller than the discrepancies with the upGPR data.This result is contradictory with the main melt phase, where the speed with which the meltwater front progresses through the snowpack is largely overestimated in the simulations.Furthermore, the distribution of melt forms in the snow profiles does not always coincide with the deeper infiltration depths detected by the upGPR.
An exception to the discussion above is snow season 2012, for which the results are consistent to a high degree.The progress of the meltwater front through the snowpack is accurately modeled by RE and only slightly less accurately by the bucket scheme for this snow season when comparing with the upGPR signal.The snow lysimeter measurements show runoff almost directly at the time the meltwa-ter front as detected from the upGPR reaches the soil.In the first snow profile made afterwards, melt forms were found for most parts of the snow cover.However, it is important to note that the progress of the meltwater front is much quicker than in the other snow seasons.Firstly, due to large snowfalls in that snow season, the snow stratification was rather homogeneous, limiting the amount of possible capillary barriers or impermeable layers in the snowpack that could hinder the liquid water flow.The relatively homogeneous stratification can be found in snow density (Supplement Fig. S17e, f) as well as in grain size (Supplement Fig. S21e, f) and grain type (Supplement Fig. S25e, f).Secondly, the onset of the snowmelt was initiated by a very warm period, leading to sufficient snowmelt to infiltrate the complete snowpack in a short amount of time.These factors all provide fewer challenges for the model.
Figure 15 also illustrates the effects of the choice of averaging method for the hydraulic conductivity at the interface nodes.The progress of the meltwater front follows a stepwise pattern.The arithmetic mean reduces the contrast in hydraulic conductivity, causing a smearing of liquid water between layers as well as over microstructural transitions inside the snowpack.The geometric mean puts more weight on the lowest hydraulic conductivity, which is found in dry snow.This results in a strengthened capillary barrier, indicated by the flatter temporal position of the meltwater front compared to the arithmetic mean.

Outlook
The extensive validation of the SNOWPACK model presented here has indicated several areas for future research and development.When focussing on processes directly impacted by liquid water flow, we can identify grain size and snow density as important properties, since they also influence hydraulic properties.It was found that during the spring melt season, both water transport schemes underestimated grain growth.Furthermore, indications were found that snow density in the melt season, which depends on the wet snow settling, is underestimated.This could be a result of either not fully representative parameterizations of these process in SNOWPACK or an underestimation of LWC in the snowpack.The latter hypothesis is supported by the comparison of bulk LWC from simulations and upGPR measurements, which has revealed an underestimation of bulk LWC in both water transport schemes on the flat site WFJ (Heilig et al., 2015).Interestingly, this underestimation was not found on slopes, which leads to the proposed hypothesis that on a flat field, capillary barriers and ice lenses may introduce stronger ponding of liquid water inside the snowpack than on slopes, where water can flow laterally.
It was also identified here that SWE depletion rates in the SNOWPACK model for the measurement site WFJ are overestimated.The SWE depletion in spring is dependent on many factors, such as snow density and wet snow settling, influencing the heat capacity, internal heat fluxes and penetration of shortwave radiation, as well as the surface energy balance and liquid water flow.These processes are difficult to investigate separately and errors could also be introduced by errors in the meteorological measurements that are used to force the model.For verifying the surface energy balance, ideally, repeated cold content measurements could be performed using the calorimetric method.However, this type of measurement is rather cumbersome to perform in the field.Accurate turbulent flux measurements would allow us to verify the parameterizations for latent and sensible heat.Snow compaction (settling) could be assessed with in situ snow harps or snow profiles at a higher temporal resolution than only biweekly.In addition, recent advances in snow micropenetrometry (SMP) are also highly promising, allowing us to achieve density measurements at high temporal and spatial resolution with relatively little effort (Schneebeli and Johnson, 1998;Proksch et al., 2015).A drawback of that method is that SMP measurements are not suitable for wet snow conditions.
The simulation of liquid water flow in snow currently only considers a one-dimensional component, assuming homogeneity in the horizontal dimension.However, this is a very strong simplification.In reality, liquid water flow exhibits strong variation in three dimensions due to preferential flow paths or flow fingering (Waldner et al., 2004;Techel and Pielmeier, 2011).The comparison of modeled liquid water flow with upGPR data and snowpack runoff measurements has identified that this simplification is indeed introducing representation errors.Numerical experiments (Hirashima et al., 2014) and laboratory observations (Katsushima et al., 2013) have provided promising indications that these processes could be described using Richards equation in three dimensions.At the same time, several processes that do appear in one-dimensional simulations, as for example the ponding of liquid water on capillary barriers, seem to be essential in forming preferential flow paths.This possibly allows for a parameterization of preferential flow in the SNOWPACK model that is closely linked to physical processes.Validation could be achieved by more detailed snow lysimeter studies, for example from measurement sites with multiple neighboring lysimeter, improved laboratory experiments or further exploiting the upGPR data.

Conclusions
The one-dimensional physics-based multi-layer SNOW-PACK model has been evaluated against measured time series and manual snow profiles for the measurement site WFJ in the Swiss Alps near Davos.Two water transport schemes, the bucket scheme and RE, were taken into consideration as well as two modes to provide the precipitation forcing for the simulations: snow height driven (15 snow seasons) and precipitation driven (18 snow seasons).Along with the im-plementation of the solver for RE, the soil module of SNOW-PACK has also been updated.Comparing simulated and measured temperatures at the snow-soil interface confirmed that the updated soil module can provide a correct lower boundary for the snowpack in the model.
The snow-height-driven simulations provide good agreement with measured snow height (RMSE around 4 cm) and, during the accumulation phase of the snow cover, with SWE.This indicates that the model adequately simulates the combination of snow settling and new snow density.In precipitation-driven simulations, the SWE in the accumulation phase exhibits a slightly larger error than in snow-heightdriven simulations, which is mainly caused by deficiencies in the precipitation undercatch correction and possibly snow drift effects.This results in a lower RMSE for snow height (20-23 cm).For the simulations at WFJ, SNOWPACK consistently overestimates the depletion rate of SWE during the spring melt season, resulting in an underestimation of SWE of typically 200 mm w.e.near the end of the snow season, accompanied by an underestimation of snow height up to 30-40 cm.In snow-height-driven simulations, this is compensated for by simulating regular snowfalls in order to match measured snow height.This procedure has as a drawback that too much mass is added to the snowpack in spring, resulting in an about 8-14 % overestimation of cumulative runoff over the snow season, whereas precipitation-driven simulations provide on average 2 % less snowpack runoff than measured.
The comparison of simulated snow density with snow density measurements made in snow profiles has shown that both the average snow density and the seasonal trend is well simulated in SNOWPACK during the main winter season.Average bias is around 25 kg m −3 and the density of deep snow layers is slightly overestimated, whereas the density of upper layers is slightly underestimated.In snow-height-driven simulations, the discrepancies grow in the melt season, when SNOWPACK underestimates snow density on average by up to 100 kg m −3 as a result of new snowfall events that are simulated to compensate for overestimated SWE depletion.The model provides simulations of grain size which are consistent with observations in manual snow profiles.Although RE causes a slight underestimation of grain size compared to the bucket scheme, snow density and grain size are adequately simulated for the parameterization of the water retention curves.
Modeled and measured snow temperatures showed a satisfying agreement with average discrepancies of around 0.5 • C. The discrepancies in the surface temperature were found to be larger, likely associated with the above mentioned underestimation of snow density in the upper layers and consequently the effect on thermal conductivity.The discrepancy in the cold content of the snow cover from simulations and field measurements was found to be small, suggesting that the surface energy balance and the soil heat flux are on average satisfactorily estimated.However, this conclusion only holds for the main winter period, as the defined cold content can only be used to assess energy budgets of snow that is below freezing.
The temporal evolution and the vertical distribution of the LWC in the snowpack differ significantly between the bucket scheme and RE.The latter provides a faster downward propagation of the meltwater front.This is accompanied by a higher r 2 value and NSE coefficient between simulated and measured snowpack runoff for the simulations with RE compared to the bucket scheme.RE also provides a higher r 2 value for the isothermal part of the snowpack compared to the manual snow profiles as well as a closer agreement with snow temperatures during the melt season.These results suggest a more accurate simulation of the progress of the meltwater front through the snowpack with RE.Although the data from the upGPR support the deeper meltwater infiltration in the snowpack in the early melt phase as simulated with RE, the opposite is found for the main wetting phase.Additionally, the distribution of melt forms in the observed snow profiles shows a higher agreement with the upGPR signal than with the simulations.Both type of observations may be considered particularly representative of matrix flow processes.The high agreement between simulations with RE and snowpack runoff therefore suggests that the use or implementation of RE is unintentionally mimicking preferential flow effects.However, the differences between both water transport schemes are relatively small compared to the differences between simulations and the observed meltwater front in the upGPR data.The results suggest that the ability of a one-dimensional approach to correctly estimate both snowpack runoff as well as the internal snowpack structure in wet snow conditions is rather limited.As the simulation of ponding of liquid water on capillary barriers and crusts is only captured with RE and not with the bucket scheme, RE seems promising however for the ability of SNOWPACK to assess wet snow avalanche risks.Future studies may also focus on the possibilities to assimilate radar-derived vertical snowpack structure (e.g., density, ice layers, liquid water) into the SNOWPACK model.This would allow us to better understand to what extent discrepancies between simulations and radar data are caused by deviations in the simulated snowpack state at the onset of snowmelt or by an insufficient process representation in the model.
The validation has shown that SNOWPACK has sufficient agreement with measurements for snow temperatures, snow density and grain size in the main winter season for a wide range of applications.When using RE, we found that the Y2012 water retention curve provides better results than the Y2010 parameterization, whereas different averaging methods to determine the hydraulic conductivity at the nodes between layers seem to have little influence.In general, several aspects of the simulations related to liquid water flow improve with RE, although often the differences between simulations tend to be smaller than differences between the simulations and the observations and the improvements are often inconsistent with the representation of the internal snowpack structure as indicated by the upGPR data.
The Supplement related to this article is available online at doi:10.5194/tc-9-2271-2015-supplement.

Figure 1 .
Figure 1.Measured and modeled snow height for different model setups (bucket or Richards equation (RE) water transport scheme, snow-height-driven (HS) or precipitation-driven (Precip) simulations, Y2010(Yamaguchi et al., 2010) or Y2012(Yamaguchi et al., 2012) water retention curves, and arithmetic (AM) or geometric mean (GM) for hydraulic conductivity) for the example snow season 2014, from October 2013 to July 2014.Note that apart from forcing with either snow height or precipitation measurements, differences between simulation setups cause only small differences in snow height simulations, resulting in overlapping lines in the figure.

Figure 2 .
Figure2.Difference in modeled and measured snow height relative to the snow season for both snow-height-driven (HS) and precipitation-driven (Precip) simulations determined over 15 and 18 years, respectively, using the bucket scheme or Richards equation withYamaguchi et al. (2012) water retention curve and arithmetic mean for hydraulic conductivity (RE-Y2012AM).For every snow season, the first day with a snow cover is set at 0 and the last day at 1.

Figure 4 .
Figure 4. Difference in modeled and observed SWE in the biweekly profiles for both snow-height-driven (HS) and precipitation-driven (Precip) simulations, using the bucket scheme or Richards equation with Yamaguchi et al. (2012) water retention curve and arithmetic mean for hydraulic conductivity (RE-Y2012AM).

Figure
Figure 11a and b show simulated snow density profiles for the bucket scheme and RE, respectively, for the example snow season 2014.In Supplement Figs.S15-S17, the other snow seasons are shown.Differences in density mainly arise when liquid water is involved.The accumulation and subsequent partial refreeze of meltwater at some layers form

Figure 14 .Figure 15 .
Figure 14.Average (a) and SD (b) of observed and modeled grain size (mm) from snow-height-driven (HS) simulations using both the Bucket scheme and Richards equation using the Yamaguchi et al. (2012) water retention curve and arithmetic mean for hydraulic conductivity (RE-Y2012AM).