An empirical model to calculate snow depth from daily snow water equivalent: SWE2HS 1.0

Many methods exist to model snow densification in order to calculate the depth of a single snow layer or the depth of the total snow cover from its mass. Most of these densification models need to be tightly integrated with an accumulation and melt model and need many forcing variables at high temporal resolution. However, when trying to model snow depth (HS) on climatological timescales, which is often needed for winter tourism-related applications, these preconditions can cause barriers. Often, for these types of applications, empirical snow models are used. These can estimate snow accumulation and snowmelt based on daily precipitation and temperature data only. To convert the resultant snow water equivalent (SWE) time series into snow depth, we developed the empirical model SWE2HS. SWE2HS is a multilayer densification model which uses daily snow water equivalent as sole input. A constant new snow density is assumed and densification is calculated via exponential settling functions. The maximum snow density of a single layer changes over time due to overburden and SWE losses. SWE2HS has been calibrated on a data set derived from a network of manual snow stations in Switzerland. It has been validated against independent data derived from automatic weather stations (AWSs) in the European Alps (Austria, France, Germany, Switzerland) and against withheld data from the Swiss manual observer station data set which was not used for calibration. The model fits the calibration data with root mean squared error (RMSE) of 8.4 cm, coefficient of determination (R²) of 0.97, and bias of −0.3 cm; it is able to achieve RMSE of 20.5 cm, R² of 0.92, and bias of 2.5 cm on the validation data set from automatic weather stations and RMSE of 7.9 cm, R² of 0.97, and bias of −0.3 cm on the validation data set from manual stations. The temporal evolution of the bulk density can be reproduced reasonably well on all three data sets. Due to its simplicity, the model can be used as post-processing tool for output of any other snow model that provides daily snow water equivalent output. Owing to its empirical nature, SWE2HS should only be used in regions with a similar snow climatology as the European Alps or has to be recalibrated for other snow climatological conditions. The SWE2HS model is available as a Python package which can be easily installed and used.

1 Introduction

Seasonal snow cover is an important variable with regard to ecology, water resource management, and the tourism industry. Accordingly, a large range of models of different complexity exist to calculate various properties of the snow cover. Traditionally, snow models were emerging from the hydrological community in order to estimate water resources from snow. Therefore, the focus was set on snow water equivalent (SWE) of the snow cover for the first simple approaches such as the empirical temperature index models in which the amount of melt is estimated by the sum of positive air temperatures (see e.g. Hock, 2003). Over time, more complex models were developed which are capable of calculating snow density, snow temperature profiles (Jordan, 1991), and snow microstructure (Lehning et al., 2002; Vionnet et al., 2012). Most of these more complex, physically based models require a rich set of input parameters such as incoming shortwave and longwave radiation, wind speed, precipitation, and temperature at subdaily temporal resolution. However, when applying models on longer timescales, e.g. for climatological analyses, the required input parameters are often limited with regard to availability and temporal resolution. Accordingly, simpler empirical models are still often used for climatological analyses instead of employing more complex physically based models. Empirical models usually do not calculate snow depth (HS) and focus on SWE which is appropriate for climatological research on changes in hydrological regimes or glacier mass balance. However, when model output is addressed to stakeholders who are usually dealing with snow depth rather than snow water equivalent, such as in the traffic and winter tourism sectors, calculation of snow depth would be desirable. This applies mostly to spatially distributed model output, because due to the ease of measuring snow depth, point data are often available for specific sites.

Snow depth is the result of SWE and the bulk snow density (ρ), where $\mathrm{SWE}=\mathrm{HS}\cdot \mathit{\rho }$. Snow depth can be measured either manually by reading from snow stakes or automatically with lasers or ultrasonic devices (Kinar and Pomeroy, 2015). While modelling SWE requires the representation of snow mass accumulation and ablation, modelling snow depth needs to address different types of densification processes. These processes involve densification due to stress induced by overlying snow and metamorphic processes that change the size and shape of the snow crystals and thus affect snow density (Anderson, 1976). Metamorphic processes can be either destructive (at constant temperature), constructive (within a temperature gradient), or melt metamorphic (for melt refreeze cycles) (Sommerfeld and LaChapelle, 1970).

All densification models need to initialize the density of a snow layer or of the whole snowpack. Since there is no simple method yet for deriving new snow density from a physical snowfall model, new snow density is either parameterized or kept at a fixed value in snow models. Various parameterizations exist and are usually based on estimating new snow density as a function of wind speed, temperature, and relative humidity (see e.g. Helfricht et al., 2018; Valt et al., 2018). When applied on a daily resolution, the quality of such parameterizations is declining due to unknown timing of a snowfall event during the day and the simultaneous occurrence of settling over the course of the day (Meister, 1985).

Several methods exist to model snow densification, either per layer or for the entire snowpack, which can roughly be classified into three categories. The first category is purely empirical whereby densification dynamics are described via exponential settling functions. This approach has first been proposed by Martinec (1956) and Martinec (1977) while variations of the method exist (e.g. Dawson et al., 2017; Koch et al., 2019; Essery, 2015; Aili et al., 2019; Brown et al., 2003, 2006). For example, Dawson et al. (2017) use a non-constant e-folding time of the settling rate based on air temperature with an additive overburden term, Essery (2015) use two different maximum densities for cold and melting snow where the exponential function converges, and Brown et al. (2006) use a maximum density based on snow depth. The second category of snow densification models is the semi-empirical method of Anderson (1976) which employs a two-stage compaction due to metamorphism and pressure from overlying snow. The compaction due to stress uses a parameterized viscosity coefficient based on temperature. Settling is enhanced when wet snow occurs in the snowpack. The scheme of Anderson (1976) has been adopted widely and is used in many snow and land surface models such as SNTHERM (Jordan, 1991), AMUNDSEN (Marke et al., 2015; Hanzer et al., 2016; Marke et al., 2018; Warscher et al., 2021), SNOWGRID-CL (Olefs et al., 2020), and CLM5 (van Kampenhout et al., 2017; Lawrence et al., 2019). Due to its need to determine wet snow in the snowpack, the method of Anderson (1976) has to be tightly integrated with a snowmelt model. The third and most sophisticated category of snow densification models is using a snow viscosity coefficient which is parameterized based on temperature and/or microstructure of the snow. Snow compaction is then modelled by applying stress due to weight of overlying snow. This requires a complex physical model that is able to represent the processes which affect e.g. snow microstructure and is realized by the two physical energy balance models, i.e. Crocus (Brun et al., 1992; Vionnet et al., 2012) and SNOWPACK (Bartelt and Lehning, 2002; Lehning et al., 2002).

To our knowledge, none of the densification models described above can easily be used as a standalone model independently of the snow model for transferring daily SWE to snow depth; however, many approaches exist to do the opposite, i.e. convert HS into SWE (e.g. Jonas et al., 2009; Winkler et al., 2021; McCreight and Small, 2014; Mizukami and Perica, 2008; Guyennon et al., 2019; Pistocchi, 2016). With new methods being developed to derive SWE from global navigation satellite system (GNSS) signal attenuation (Koch et al., 2019) or by cosmic ray attenuation (Gugerli et al., 2019), it would be even more desirable to be able to model snow depth from the derived SWE data (Capelli et al., 2022). For climatological applications, quantile mapping can be used to statistically correct the SWE output of simple melt and accumulation models that can be run far back in time using more complex models that require more input data and are therefore not suitable for climatological timescales up to several decades (Michel et al., 2023). To transfer the statistically corrected SWE field to HS, a standalone densification post-processing model would also be required. For these reasons, we developed a simple empirical snow densification model which uses daily SWE as sole forcing and transforms SWE to HS using exponential settling equations for individual layers inspired by Martinec (1977); Dawson et al. (2017); Koch et al. (2019); Essery (2015). We make an implementation of the model available as a Python package which can be downloaded and installed from the Python packaging index (PyPI).

Figure 1 Illustrative snowpack evolution within the SWE2HS model calculated from 35 d of synthetic SWE forcing data. The top panel shows the calculated snow depth (top solid black line) and the two generated layers (grey areas) separated by a layer boundary (middle solid black line). The dashed black line represents the exponential settling for layer 1 according to Eq. (1) without the changes in ρ_max due to the overburden from layer 2 and the decrease in SWE. The red area therefore represents the effects of changing the overburden and wetting on the densification dynamics of layer 1. The letters on the x axis indicate the following: (a) initial increase in SWE – layer 1 is created with density ρ_new at creation time. (a, b) The density of layer 1 increases according to Eq. (1), ρ_max is set to a higher value than ρ_max,init since half the weight of layer 1 is treated as its overburden. (b) SWE increases and layer 2 is created. For layer 1, ρ_max increases according to Eqs. (2) and (3). (b, c) Both layers settle after Eq. (1) with the updated ρ_max for layer 1. (c) SWE starts to decrease and ρ_max starts to increase after Eq. (4) for both layers. (c, d) Layer 2 is removed gradually in order to compensate for the decreasing SWE, both layers settle according to Eq. (1) and ρ_max is updated for each time step according to Eq. (4). (d) Layer 2 is completely removed. Subsequently, layer 1 is gradually removed in order to account for the further decrease in SWE.

The remainder of the paper is structured as follows. In Sect. 2, we describe the model and the technical implementation. In Sect. 3, we describe the data used for calibration and validation of the model alongside the used calibration methods. In Sect. 4, we show the performance of the calibrated model in alpine snow environments and discuss the scope and limitations of the model in Sect. 5.

2 Density model

2.1 Basic concept

The density model SWE2HS calculates snow depth at a daily resolution and is driven by the daily snow water equivalent of the snow cover only. In the following, we use the unit millimetre water equivalent ($\mathrm{mm}\mathrm{w}.\mathrm{e}.$) for SWE. The model creates a new layer with a fixed new snow density ρ_new for every increase in SWE such that, over time, a snowpack of individual layers builds up. The density of a layer increases with an exponential decay function towards a time-varying maximum density. The maximum density starts with an initial value at creation time of the layer and subsequently increases towards a higher value based on the overburden a layer has experienced and the occurrence of SWE losses in the snowpack. When SWE decreases, the model removes layers from the top of the snowpack. The layer number n can thus undergo changes over time based on the number of SWE increases and losses in the snowpack. The model neglects constructive metamorphism, refreezing, and is not able to capture rain-on-snow events which might lead to an increase in SWE but no increase in snow depth.

2.2 Settling mechanisms

The density of a layer at day i asymptotically converges towards the time-varying ρ_max of the layer via the following exponential function:

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_i={\mathit{\rho }}_{\text{max}}-({\mathit{\rho }}_{\text{max}}-{\mathit{\rho }}_{i-\mathrm{1}})\cdot \exp \left({\displaystyle \frac{-\mathrm{1}}{R}}\right),\end{array}$$

where ρ_i is the density of day i and ρ_i−1 is the layer density of the day before. The settling resistance (e-folding time) R is a model parameter which is optimized in model calibration.

The maximum density to which the density of a snow layer converges, ρ_max in Eq. (1), also evolves over time. We model the maximum density of a snow layer based on three assumptions. The first assumption is that snow which has experienced a higher load reaches a higher maximum density. The second assumption is that a snow layer is initially dry and that wet snow has a higher maximum density than dry snow. The third assumption is that the time-varying maximum density cannot decrease. Accordingly, the maximum density of a snow layer undergoes changes during its lifetime and transitions from the model parameter ρ_max,init to the model parameter ρ_max,end. At the time of deposition, the layer has a theoretical maximum snow density of ρ_max,init. Afterwards, ρ_max increases towards ρ_max,end by two mechanisms as described following.

1. If a layer experiences an overburden σ > 0 mm, its maximum density ρ_max is increased linearly with overburden. We calculate σ as a proxy for overburden stress by summing the amount of SWE above a layer and a half of the SWE of the layer itself. If the overburden weight is equal or larger than the model parameter σ_max, ρ_max is capped at ρ_max,end.

   $$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\mathit{\rho }}_{\text{max}}^{\ast }=\\ & \left\{\begin{array}{ll}\frac{({\mathit{\rho }}_{\text{max,end}}-{\mathit{\rho }}_{\text{max,init}})}{{\mathit{\sigma }}_{\text{max}}}\cdot \mathit{\sigma }+{\mathit{\rho }}_{\text{max,init}}& \text{if }\mathit{\sigma }<{\mathit{\sigma }}_{\text{max}}\\ {\mathit{\rho }}_{\text{max,end}}& \text{if }\mathit{\sigma }\ge {\mathit{\sigma }}_{\text{max}}\end{array}\right.\end{array}\end{array}$$

   If the updated ${\mathit{\rho }}_{\text{max}}^{\ast }$ is lower or equal than the value of the day before (ρ_max i−1), the value of the current day (ρ_max i) is not updated, which could otherwise cause a decrease of ρ if the density ρ_i−1 equals ρ_max i−1 (see Eq. 1).

   $$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\mathit{\rho }}_{\max ,i}=\\ & \left\{\begin{array}{ll}{\mathit{\rho }}_{\text{max}}^{\ast }& \text{if }{\mathit{\rho }}_{\text{max}}^{\ast }\in ({\mathit{\rho }}_{\max i-\mathrm{1}},{\mathit{\rho }}_{\text{max,end}})\\ {\mathit{\rho }}_{\max i-\mathrm{1}}& \text{if }{\mathit{\rho }}_{\text{max}}^{\ast }\le {\mathit{\rho }}_{\max i-\mathrm{1}}\end{array}\right.\end{array}\end{array}$$

2. SWE losses are defined by ${\mathrm{SWE}}_i-{\mathrm{SWE}}_{i-\mathrm{1}}<\mathrm{0}$. Whenever SWE in the snowpack decreases, we assume that the entire snowpack is wet since we attribute all SWE losses to runoff. In doing so, we neglect losses in SWE due to sublimation. If SWE decreases, we assume that melt metamorphism is active and the maximum snow density ρ_max of each layer is increased towards ρ_max,end by

   $$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{\rho }}_{\max i}=& {\mathit{\rho }}_{\text{max,end}}-({\mathit{\rho }}_{\text{max,end}}-{\mathit{\rho }}_{\max i-\mathrm{1}})\\ & \times \exp (-v_{\text{melt}}),\end{array}\end{array}$$

   where ρ_max i is the maximum density of day i, ρ_max i−1 is the maximum density of the day before, and v_melt is a model parameter for the speed of that transition.

At the end of every time step, the snow depth of the snowpack is calculated by summing up the thickness of all n snow layers in the snowpack as

$$\begin{array}{}\text{(5)}& \mathrm{HS}=\sum _{k=\mathrm{1}}^n{\displaystyle \frac{{\mathrm{SWE}}_k\cdot {\mathit{\rho }}_{\text{water}}}{{\mathit{\rho }}_k}},\end{array}$$

where ρ_water is the density of water with 1000 kg m⁻³, and SWE_k and ρ_k are SWE and density of layer k, respectively. All free model parameters that need to be calibrated are listed in Table 2.

2.3 Technical implementation

We provide an implementation of the model as a Python package under GNU General Public License v3.0 (GPLv3). One-dimensional station data and two-dimensional model grids of daily SWE time series can be transformed to snow depth with the snowpack evolution described above. Additionally, a step-by-step processing mode with caching of the model state variables for two-dimensional SWE grids of consecutive days is available for operational applications. Python, being a high-level, interpreted general-purpose programming language has been chosen due to its easy-to-read syntax, growing user base, and community support for scientific computing and data analysis. Our implementation is using the just-in-time Python compiler Numba (Lam et al., 2015) for increasing runtime efficiency. Additionally, it depends on the libraries NumPy (Harris et al., 2020) for numerical computations, Pandas (Reback et al., 2022) for one-dimensional input series, and xarray (Hoyer and Hamman, 2017) for multidimensional input grids. The multidimensional, distributed versions of the model can make use of Dask (Dask Development Team, 2016), which makes it possible to execute the model in parallel on standalone computers or high-performance computing environments. Processing 23 years of daily SWE data from the Swiss 1 km × 1 km domain (8401 pixel × 365 pixel × 272 pixel), generated with the method of Michel et al. (2023), took ∼ 10 min on a desktop PC (8 cores, Intel Core i7-4790 CPU @ 3.60 GHz, 24 GB RAM) including file IO. The model implementation can be installed from the official third-party software repository for Python, the Python Package Index (PyPI: https://pypi.org/project/swe2hs, last access: 20 October 2022); the source code of SWE2HS is hosted on a Gitlab instance of the Swiss Federal Institute for Forest, Snow and Landscape Research WSL (https://code.wsl.ch/aschauer/swe2hs, last access: 20 October 2022); and the software version which was used for this publication is archived at https://doi.org/10.5281/zenodo.7228066 (Aschauer, 2022).

3 Model calibration and validation

As for every empirical model, parameters in our density model need to be calibrated. Calibrated parameters may differ depending on the station, snow type, and snow climatological setting. Here, we try to find one single generic optimal parameter set which suits most snow climatological conditions in Switzerland and the European Alps in general. We do so by calibration of a data set which covers a large range of different altitudes and climatologic settings in Switzerland (see Sect. 3.3.1) and test the found parameters on two other independent data sets compiled from automatic weather stations (AWSs) in the European Alps (see Sect. 3.3.2) and from withheld data of the calibration data set which was not used for calibration.

3.1 Calibration and validation methods

Our model has 6 model parameters which need to be calibrated. Before calibration, we define upper and lower bounds of possible values for each model parameter (see Table 2) and apply the constraints that ρ_max,init needs to be smaller than ρ_max,end and ρ_new needs to be smaller than ρ_max,init. The chosen upper and lower parameter bounds are based on literature (e.g. for new snow density ρ_new) or based on our previously gathered experience with the model for the model specific parameters such as the settling resistance R.

For parameter calibration, we use the differential evolution algorithm which is a stochastic population-based method for minimizing nonlinear and non-differentiable continuous space functions as implemented in SciPy (Storn and Price, 1997; Virtanen et al., 2020). We chose differential evolution due to its gradient-free nature and ability to overcome local minima (Storn and Price, 1997). After an initial Sobol' sequence sampling (Sobol', 1967), the algorithm draws parameter candidate samples from the parameter space by mutating the current best member of the sample population with the difference of two other randomly chosen members. After the global optimization with the differential evolution algorithm, the result is refined by the L-BFGS-B method of Byrd et al. (1995) which is a quasi-Newton method that estimates the Hessian of the objective function based on the recent parameter sample history and can handle bound constraints.

We optimize the model by minimizing the root mean squared error (RMSE) which is a measure of the distance between the predicted values from the model $\widehat{\mathit{y}}$ to the reference y. It is defined as

$$\begin{array}{}\text{(6)}& \mathrm{RMSE}(\mathit{y},\widehat{\mathit{y}})=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\left(\widehat{y_i}-y_i\right)}^{\mathrm{2}}},\end{array}$$

where y_i and $\widehat{y_i}$ are the ith element of the $i=\mathrm{1},\mathrm{\dots },n$ elements in y and $\widehat{\mathit{y}}$, respectively. Additionally, we use the two statistical error measures, coefficient of determination (R²) and bias, to evaluate the model. The R² score is representing the proportion of variation in the data y that can be predicted from the model and is defined as

$$\begin{array}{}\text{(7)}& R^{\mathrm{2}}(\mathit{y},\widehat{\mathit{y}})=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(y_i-{\widehat{y}}_i\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}},\end{array}$$

where $\widehat{y_i}$ is the predicted value of the ith sample, y_i being the associated reference value for total n samples and $\stackrel{\mathrm{‾}}{y}$ being the mean of y. The bias is a measure for the systematic tendency of a model to over- or underrepresent the reference data. Therefore, it has large implications in climatological contexts. We calculate the bias for a sample of size n as follows:

$$\begin{array}{}\text{(8)}& \mathrm{bias}(\mathit{y},\widehat{\mathit{y}})={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{y}_i-{\widehat{y}}_i,\end{array}$$

where $\widehat{y_i}$ is the predicted value of the ith sample and y_i is the associated reference value. All presented score values RMSE, R² and bias are calculated only for the subset of $\widehat{\mathit{y}}$ and y where any of the two vectors is not zero.

3.2 Sensitivity analysis methods

In order to assess the importance of individual model parameters on the result, we perform a sensitivity analysis on the validation data set by calculating R² and bias on 114 688 parameter sets sampled after the method of Saltelli (2002). This method expands on the low-discrepancy quasi-random Sobol' sequence and generates uniformly distributed samples of the parameter hypercube space. Afterwards, we run the model with the sampled parameter sets on the validation data set from the automatic stations (see Sect. 3.3.2), calculate the error measures for R² and bias for each parameter set, and compute global sensitivity indices (S_Ti) after Sobol' (2001). These indices give an estimate about the proportion of variance in R² and bias that can be attributed to a model parameter and all its interactions with other model parameters. We perform the sensitivity analysis within the Python framework SALib of Herman and Usher (2017).

3.3 Calibration and validation data

3.3.1 Data from Swiss manual observer stations

To calibrate the SWE2HS model, we use data from 58 Swiss manual observer stations between 1080 and 2620 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ operated by the WSL Institute for Snow and Avalanche Research SLF (Marty et al., 2017). Snow depth is measured daily with a snow stake and SWE is measured every 2 weeks in a snow pit on the same site. In order to get daily SWE data, HS data are transformed to SWE with the Δsnow model of Winkler et al. (2021); Δsnow is a semi-empirical multilayer snow compaction model that uses a continuous time series of snow depth to calculate SWE and bulk snow density. The model calculates snow compaction by treating snow as a Newtonian viscous fluid, is able to include transient settlement due to new snow loading, and accounts for measurement errors in the input data. Once a defined maximum density is reached in a layer, the snow is melted and distributed from top to bottom in the snowpack. Runoff is triggered when all layers reach the maximum snow density. The model has 7 free parameters which were calibrated in Winkler et al. (2021) with data from 14 stations at different elevations and snow climatological regions in the European Alps. The Δsnow model performs very well in modelling the temporal evolution of SWE on the daily scale (Fontrodona-Bach et al., 2023). Additionally, the accuracy of the Δsnow model is already high for the original unified Alpine-wide parameters as shown by Winkler et al. (2021) and confirmed by Fontrodona-Bach et al. (2023). In order to improve this accuracy of the daily SWE time series, the Δsnow model parameters were optimized for each station individually using the biweekly SWE measurements from the manual observer profiles. Due to its destructive nature, the snow pit is not at the exact same location as the snow stake and consequently the profile cut height can deviate from the measured height at the snow stake. Therefore, the biweekly SWE data were corrected by calculating the bulk density from the profile and applying it to the measured height from the snow stake. Δsnow model parameter optimization was done by minimizing the RMSE between modelled SWE and corrected SWE from the profiles while we allowed the Δsnow parameters ρ_max (maximum density) to vary between 300 and 600 kg m⁻³, ρ₀ (new snow density) to vary between 65 and 135 kg m⁻³, and the remaining parameters to vary by ± 25 % from the optimized value found in Winkler et al. (2021). For optimization, we again used differential evolution as described in Sect. 3. In order to further increase the reliability of the calibration data set, we only kept station winters with more than 2 SWE profiles and RMSE below 7.5 mm in the resulting daily SWE data set from the Δsnow model. Since we did not want certain stations with long SWE and HS records to bias the calibration, we shortened the length of station records longer than 15 years by randomly selecting 15 water years from the full station record. The resulting set consists of 741 station years from 60 stations. Compared with the biweekly manual SWE measurements, the modelled daily SWE calibration set has an RMSE of 30.0 mm and bias of −1.09 mm in the calibration data set. However, we cannot assess the uncertainty for the dates between the biweekly SWE measurements.

While we are aware that it might be preferable to calibrate a model on measured data instead of output from another model, we still chose the approach described above in order to have an exhaustive calibration data set which (a) covers a wider range of altitudes, expositions, and snow climatic settings in our target region; (b) does not have problems of potential over- and under-measurement from automatic SWE measurement devices (Johnson and Marks, 2004); and (c) does not contain measurement noise which is not properly tackled in the SWE2HS model.

As a first validation data set, we use the remaining years of the long station records which were shortened to 15 water years for the calibration data set (see above). This calibration set of manual observer station data contains 1279 station years from 42 different stations.

Figure 2 Locations of the snow measurement sites in the European Alps from which snow water equivalent and snow depth data were used to compile the calibration and two validation data sets. Background colouring resembles elevation. GMTED2010 elevation data, courtesy of the U.S. Geological Survey, was used to create the map (Danielson and Gesch, 2011).

Lejeune et al. (2019)Krajci et al. (2017)Hagen et al. (2023)

Table 1 Automatic weather stations from which we used snow water equivalent and snow depth data for validation of the model. The number of years refers to complete hydrological years (September–August) included after data cleaning, the average snow depth ($\stackrel{\mathrm{‾}}{\text{HS}}$) is calculated in the winter months from November to April.

¹ WSL Institute for Snow and Avalanche research SLF. ² Bavarian avalanche warning centre.

3.3.2 Data from automatic weather stations in the European Alps

As a second validation data set, we gathered data of 10 different automatic weather stations (AWSs) in Austria, France, Germany, and Switzerland that automatically measure SWE with either a snow pillow or a snow scale and measure snow depth with an ultrasonic measurement device at subdaily resolution (see Table 1). The raw SWE and HS data with a temporal resolution ranging between 15 min and 1 h were resampled to daily resolution by taking the median of all measurements between 06:00 and 08:00 LT. Any systematic offset errors from raw sensors were corrected by subtracting the mode of the summer months (MJJASON) from the SWE or HS time series of each hydrological year. Missing data gaps shorter than 5 consecutive days in SWE have been filled by linear interpolation. In total, 77 short SWE gaps were filled, including 39 1 d gaps, 17 2 d gaps, 8 3 d, and 4 d gaps, and 5 5 d gaps. For longer gaps, the time period before the gap in the respective hydrological year is included and data after the gap are discarded. Short gaps in snow depth are accepted since it is not required to drive the density model but for evaluating the quality of the model. Missing HS data points will thus not be included when calculating any score metrics. All hydrological years that were included in the final validation data set have been quality checked by visual inspection and by ensuring the bulk density stays below 700 kg m⁻³.

Figure 3 Schematic modelled snowpack evolution for 6 different station years from the validation data set from automatic stations with different altitudes and snowpack thicknesses (see Table 1). In the top panels, the dotted red line is the measured snow depth (HS), the solid black line bounding the coloured area is the modelled snow depth, the thin black lines represent the layer boundaries within the modelled snowpack, and the colouring refers to the modelled layer densities. The lower panels show the daily snow water equivalent time series used to force the model in each station year.

Table 2 Parameters of the model, lower and upper bounds during calibration, and optimized value.

4 Results

The model calibration on the data set described in Sect. 3.3.1 yielded the optimized parameters listed in Table 2. Figure 3 shows the temporal evolution during 6 example winters for stations in the validation data set from AWSs calculated with the optimized parameter set. Looking at the temporal development of the density and layering in the modelled snowpack, the density rapidly increases in the first few days after layer creation leading to enhanced settlement in this period. Additionally, the density of a layer reacts to changes in SWE in the overlying layers (see e.g. bottom three layers, end of January 2021 for station Kühroint). For the deep snowpacks in the top two panels in Fig. 3, the lower layers are no longer able to respond to the overburden and no further compaction occurs once a certain amount of snow is on top (see Eq. 2). The same settling dynamics in the snowpack can be observed for the calibration data set and validation data set from manual stations (see Figs. A2 and A1).

Figure 4 Scatterplots of modelled versus measured snow depth values for (a) the calibration data set from Swiss manual observer stations, (b) the validation data set of the remaining years from the Swiss manual observer stations which were not used for calibration, and (c) the validation data set from automatic snow stations. The dashed red line is a linear fit to the data, the solid black line represents perfect predictions. In the insets, the same data are shown as bivariate histograms indicating the density of the scatter points. Score values of the shown data are listed in Table 3.

Table 3 Score values of RMSE, R², and bias for the calibration data set and the two validation data sets after parameter calibration. The accompanying data are visualized in Figs. 4 and 5.

With the optimized parameter set, the model is able to fit the calibration data with RMSE of 8.4 cm, R² of 0.97, and negligible bias of −0.3 cm (see Table 3 and Fig. 4, left panel). The seasonal evolution of the bulk density can be reproduced well on the calibration data set with monthly R² scores larger than 0.93 and June being the only month with considerable underestimation of the median snow depth (Fig. 5, left panels). On the validation data set from manual stations, the model is able to achieve the same performance as on the calibration data set with RMSE of 7.9 cm, R² of 0.97, and bias of −0.3 cm. Monthly R² values for the validation data set from manual stations are above 0.94 for the winter months of November to May. The R² values are lower for the months of October and June with 0.88 and 0.89, respectively. On the validation data set from AWSs, the performance is weaker than for the calibration data set and validation data set from manual stations with RMSE of 20.5 cm, R² of 0.92, and bias of 2.5 cm. The model slightly underestimates the median snow depth in February and March on the validation data set and overestimates the median snow depth in the ablation months April, May, and June. Monthly R² score values are above 0.79 for all months except October, where R² is below 0.4. On the calibration data set, R² is above 0.95 for 75 % of the stations and below 0.8 for only two stations (see Fig. 6). The R² scores per station in the validation data set from manual stations are similarly distributed as for the calibration data set; however, for none of the stations the R² score is below 0.8. On the validation data set, R² per station varies between 0.15 and 0.95 and is larger than 0.75 for 75 % of the stations. On the calibration data set, the bias per station is uniformly distributed around 0 and for all but two stations, it is smaller than ± 10 cm. On the validation data set, the bias per station ranges from −7 to 22.7 cm with four stations having a positive bias larger than 10 cm.

According to the sensitivity analysis, the settling resistance factor R is the most important model parameter with a global sensitivity index of 0.44 and 0.43 for R² and bias, respectively (Fig. 7). This means that within the 114 688 samples drawn during the sensitivity analysis, 44 % and 43 % of the proportion of variance in R² and bias can be attributed to the settling resistance factor R, respectively (see Sect. 3). For R², the second most influential model parameter is new snow density ρ_new followed by the final maximum snow density ρ_max,end. For bias, the model is less sensitive to ρ_new than for R². The model is relatively insensitive to changes in the model parameters ρ_max,init, v_melt, and σ_max with total sensitivity indices below 0.1 for both R² and RMSE.

Figure 5 Boxplots comparing the distributions of measured and modelled data in the months from October to June for (a) the calibration data set from Swiss manual observer stations, (b) the validation data set of the remaining years from the Swiss manual observer stations which were not used for calibration, and (c) the validation data set from automatic snow stations. A box spans the lower and upper quartile of the data with a line at the median. The whiskers extend to the last datum within 1.5 times the interquartile range while the points represent outliers past the range of the whiskers. The lower three panels show monthly R² scores for the modelled data to the reference, calculated for each data set.

Figure 6 Boxplots of the scores R² and bias calculated individually on the data from each station in the calibration and two validation data sets. The black dots show the underlying data from which the boxplots were calculated.

Figure 7 Global Sobol' (2001) sensitivity indices (S_Ti) calculated for R² and bias of the model predictions on the validation data set from automatic stations for 114 688 samples drawn with the method of Saltelli (2002).

5 Discussion

5.1 Model design and complexity selection

On our way towards the model presented here, we tried models of different complexity and we included and removed processes while iterating back and forward. Some prototype model versions additionally included daily temperature as input forcing, which we tried to use for parametrization of new snow density and onset of the wetting from the top of the snowpack by using the cold content parameterizations used in Scheppler (2000) and Szentimrey et al. (2012). Other versions parameterized the settling resistance R based on the overburden or density of a layer or a combination of the two factors. In order to do an objective model selection for a final model version, in addition to the scores defined in Sect. 3, we calculated the Akaike information criterion (AIC) by using the RMSE as an estimator for the maximum value of the likelihood function. The AIC is a statistical error measure that penalizes larger numbers of free model parameters (Akaike, 1998; Cavanaugh and Neath, 2019). We then ranked the four different scores for the optimized model of each version and averaged the score ranks over the calibration and validation data set. This allowed us to make an informed decision regarding which model version to use. In order to avoid overfitting on the calibration data set, we gave more focus to the validation data set.

We assumed that it would be beneficial to use temperature in the beginning of model conceptualization and one could argue that when using SWE from an accumulation and ablation model, there is always at least daily mean temperature available used to drive a melt model. However, when quantitatively assessing the model versions with parameterized new snow density or cold content parameterizations, we did not see model improvement from the daily mean temperature inclusion and thus decided to only use SWE. This additionally comes with the asset that the model can be plugged in as a post-processing tool to any snow model which outputs daily SWE. Besides the best performance, another important factor to keep the model simple was to reduce the risk of equifinality, meaning that an optimal solution can be achieved through different states, i.e. parameter combinations of the model (Beven and Freer, 2001). In a very early version of the model, we had implemented 3 modes of layer removal: from the top, proportional, and from the bottom. In a deep snowpack, we assume that most melting occurs at the top due to energy input from shortwave and longwave radiation. There is also secondary melting due to geothermal heat fluxes, but we assume that this is negligible compared to the energy input from the top. During melting, SWE is not actually “removed” from the top of the snowpack, unless it is sublimating directly into the air. In fact, melted liquid water is assumed to percolate down through the snowpack where it can potentially refreeze if the snowpack has not yet become isothermal. However, in a simple empirical model such as SWE2HS, it is not possible to implement such processes in a meaningful way and for shallow snowpacks, other processes and predominant melt mechanisms might be relevant. Nevertheless, predominant melting at the top of the snowpack could imply that the top layer should be removed first for deep snowpacks. These considerations, together with the fact that the layer boundaries looked most realistic in comparison to measured settlement curves for both shallow and thick snowpacks, led us to the choice of layer removal from the top presented here. This has the effect that new snow events late in the season followed by an immediate decrease in SWE have a short-lived effect on the bulk snow density. We would also assume that the top-down removal of layers reduces the sensitivity of the model to new snow density in general.

5.2 General remarks and limitations

As shown in Sect. 4, the model is able to fit the calibration data very well. The calibration data have been compiled from manually measured snow depth data and modelled SWE data with the Δsnow model of Winkler et al. (2021). Therefore, e.g. the occurrence of rain-on-snow events cannot degrade the model skill since the Δsnow model is also not able to represent rain-on-snow events. We still consider it a valid approach to use modelled data for calibration as we hold back a set of measured data for independent validation of our optimized model parameter set (see Sect. 3.3.2). Additionally, the main scope of application of SWE2HS is post-processing output from simple accumulation and melt models. The model performs better on the Δsnow data set compared to the data set from AWSs. This is due to several reasons. First and foremost, the model has been calibrated on the data in the Δsnow data set and not on the data in the data set from AWSs. Accordingly, the fitted parameters do not necessarily suit the data in the data set from AWSs. Another reason is that the calibration data set and validation data set from AWSs cover different spatial domains. The scores on the validation data set from manual stations, which are comparable to those on the calibration data set, show that the model is in principle able to perform on unseen data not used for calibration. Another potential reason why the scores on the validation data set from AWSs might differ from those on the calibration data set is the fact that the calibration of the SWE2HS model is made with modelled data from the Δsnow of Winkler et al. (2021) (see Sect. 3.3.1). These modelled input data are also attributed with a considerable uncertainty compared to the biweekly manual SWE data (RMSE of 30.0 mm) which can influence the calibrated parameters in a way that degrades the model skill on the data set from AWSs. Other reasons for weaker model performance in the data set from AWSs potentially arise from noise and measurement uncertainties. One source of these uncertainties are problems of over- and under-measurement from the automatic SWE measurement devices (Johnson and Marks, 2004). This uncertainty increases with time during the winter and could explain the overestimation in the ablation season. Additionally, the measurement uncertainty of the automatic SWE and HS data can cause small changes in SWE and HS, which are not physically based (Capelli et al., 2022). In this regard the SWE and HS in the calibration data set are much more consistent. We could not include a mechanism to deal with measurement uncertainties analogous to Winkler et al. (2021) since a SWE time series does not contain any information on settlement which could be used to correctly distinguish a signal from noise. A last source of uncertainty in the data set from AWSs is that the automatic SWE measurements are not necessarily located at the exact same place as the snow depth measurements and we did not have a way to correct this error in the same way as we did for the manually observed SWE and HS measurements. The low R² score in October for the data set from AWSs (see Fig. 5) could be explained by the increased importance of new snow density in the beginning of the snow season. In this data set, October snow cover is often characterized by ephemeral snowfall events, where new snow density is more important due to the lack of settled older snow layers. New snow density can have a high variability (see e.g. Helfricht et al., 2018), which is not predictable if new snow density is kept fixed or independent of temperature or similar, as is done in the SWE2HS model.

Other sources of uncertainty are due to inherent limitations of our empirical modelling approach. As mentioned above, the model is not able to represent rain-on-snow events. In the exemplary snowpack evolution of the winter of 2020–2021 at station Kühroint, an increase in SWE causes modelled snow depth to increase although the measured snow depth is constantly decreasing during this time (middle of March 2021, Fig. 3, middle left panels). This could be an example of either erroneously measured SWE or a rain-on-snow event which caused an increase in SWE but not in HS. The latter seems more likely since, along with the increase in SWE, settling is enhanced for the measured HS and (as additional not-shown data demonstrate) the temperature is rising above 0 ^∘C in combination with precipitation. Morán-Tejeda et al. (2016) show that such events are rare and contribute to maximum 100 mm for elevations above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ With a changing climate, rain-on-snow events might become more likely above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ but might decrease for low altitudes as a decrease in rainfall and shorter snow cover duration are thought to counteract increased temperatures. In the SWE2HS model, we assume that the snowpack is completely wet when SWE decreases. However, prior to complete wetting, a wetting front propagates through the snowpack from top to bottom over time (Marsh and Woo, 1984). The choice to use only SWE as a forcing leads to the limitation that we can only detect the wetting front that reaches the bottom, which will be observable as a negative change in SWE. Therefore, the model likely misses the onset of melt metamorphism in the upper snowpack and overestimates HS during this time. In a post-processing setting, this could be addressed by looking forward in time and assuming wetting prior to the SWE decrease, but this design choice would make it impossible to use the model in an operational setting when no information about future snowpack evolution is available. We try to partly compensate for this flaw by increasing the maximum density with increasing overburden. As R² is not decreasing in spring (Fig. 5) for the calibration and validation data set, the model nevertheless seems able to predict snow depth during the ablation period reasonably well. By assuming that the entire snowpack is wet when SWE decreases, we neglect processes other than runoff that lead to decreases in SWE. Sublimation, the direct phase transition from solid water to water vapour, is one of these processes that can have a remarkable impact on the alpine water budget, especially at wind-exposed high elevations and in forested areas Strasser et al. (2008). The inability to represent sublimation in the SWE2HS model can lead to an underestimation of snow depth because SWE losses lead to an increase of the maximum snow density (ρ_max) in Eq. (1). This limitation should be kept in mind when applying the model to conditions where sublimation processes play an important role.

Since the model is of empirical nature, the parameter set which is presented here for the European Alps might not be suited for other regions on earth with different climatologic conditions. If applied to other regions, the model parameters need to be calibrated again. However, as we never tested the model in e.g. Arctic regions, we cannot make any statements about whether the model is able to represent settling dynamics in these snow climatologic conditions even if it would be calibrated on data from there. Although the calibrated parameter set presented in this paper is thought to be representative for the European Alps, it may not be suitable for some stations in the validation data set with bias of up to 22.7 cm (see Fig. 6). This high positive bias is occurring at station Davos (only 1 year of data) and station Laret (2 years of data). We suspect that in all 3 years, a short period of potential positive measurement errors from the SWE sensor caused the snow cover evolution to become defective. In order to achieve the best model skill for a single location, recalibrating the model to data from that location is necessary.

Simple empirical models that try to conceptualize processes in a non-physical way are often subject to the risk of potential model equifinality (Beven and Freer, 2001). While we chose to present a single calibrated parameter set in Sect. 4, there is still the risk that other parameter sets might lead to equally good predictions. However, an in-depth analysis of potential equifinality is out of scope of this model description paper and might be the subject of future work.

The model is less sensitive to changes in the model parameter ρ_new (new snow density) for bias than for R². This is likely due to the nature of the two statistical error metrics used. The R² measures the proportion of variance the model is able to explain in the data and the bias measures whether the model is on average under- or overestimating the data to be predicted. The new snow density is mainly affecting the model result for the times shortly after increases of SWE and thus is not as important for model bias as for R². Due to the sensitivity to new snow density and the fixed new snow density approach in the model, it is not reasonable to derive climatologic indices related to the amount of new snow such as the maximum increase in HS during 3 d from model output.

The SWE2HS model is tailored for use with daily resolution SWE data. When attempting to use the model with higher temporal resolutions such as hourly, additional processes to those considered in the model become increasingly important and additional parameters such as radiation and temperature are likely to be required to satisfactorily represent densification. For example, the new snow density will be much more variable on shorter timescales, and it is likely that the fixed new snow density approach used in the SWE2HS model will not be sufficient at hourly resolutions. In addition, the empirical transition rate from dry to wet snow (v_melt) would have to be changed when adapting the model to higher temporal resolutions.

6 Conclusions

We present a simple snow density model which can be used to transfer continuous daily snow water equivalent data to snow depth. The empirical multilayer model uses exponential settling equations, a fixed new snow density, and assumes a changing maximum snow density over time based on overburden and SWE losses. The model was calibrated with a gradient-free evolutionary algorithm on a data set from the Swiss Alps that was generated from biweekly SWE and daily HS records. Prior to calibration, the biweekly SWE records were converted to daily values with the Δsnow of Winkler et al. (2021). On the calibration data, the model is able to reproduce the measured snow depth with RMSE of 8.4 cm and bias of −0.3 cm. The SWE2HS model is validated on multiyear data from 10 automatic snow stations between 1100 and 2500 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ in the European Alps where it can reproduce the measured data with RMSE of 20.5 cm and bias of 2.5 cm. In addition, the model is validated against withheld data from the Swiss manual observer station data set not used for calibration, on which it achieves results as good as on the calibration data set, with RMSE of 7.9 cm, R² of 0.97, and bias of −0.3 cm. Due to its simplicity, SWE2HS can be used for climatological use cases where input data for more sophisticated densification models are sparse. Since the only input needed to drive the model is daily SWE, it can also be used to post-process model output from any other snow model or to transfer SWE data obtained from automatic SWE sensors.

Appendix A: Additional figures

Figure A1 Schematic modelled snowpack evolution for 6 different station years from the manual station validation data set. Winters from stations with different elevations and with differing snowpack thicknesses are shown. For an explanation of the figure, the reader is referred to the caption of Fig. 3.

Figure A2 Schematic modelled snowpack evolution for 6 different station years from the manual stations calibration data set. Winters from stations with different elevations and with differing snowpack thicknesses are shown. For an explanation of the figure, the reader is referred to the caption of Fig. 3.