Analysis and modelling of snow bulk density in the Tyrolean Alps

To further reduce measurement costs, efforts are being made to develop methods for the estimation of snow density (ρ), which is required to convert HS into SWE – which is more relevant in terms of hydrology. Densification is governed by snow metamorphism, with the major physical processes being the mechanical damage of snow crystals during precipitation, the transport of water vapour caused by surface energy differences and thermal gradients within the pack, melting and refreezing as well as pressure consolidation (Sommerfeld & LaChapelle 1970; Schneebeli & Sokratov 2004). Physical snowpack models (e.g., Brun et al. 1989; Bartelt & Lehning 2002) express these processes explicitly but are demanding with regard to input data and computational resources. By contrast, several empirical methods (e.g. Rohrer & Braun 1994; Jonas et al. 2009; Sturm et al. 2010; López-Moreno et al. 2013) are far less data demanding, but still reliable. Hence, this type of model is especially suitable for practical applications.

The density of new snow is related to temperature and wind speed (Schmucki et al. 2014). A rich snow dataset from the Swiss Alps exhibited a mean new snow (HN) density of 100 kg m⁻³ (Rohrer et al. 1994). Settling curves of HN as a function of time were applied for each snowfall event to model SWE with consideration of daily records of HN and HS (Martinec & Rango1991; Rohrer & Braun 1994). Weather data were solely used for the statistical modelling of ρ and snow loads in Norway (Meløysund et al. 2007), but this application required input data comparable to physical snow models (e.g., relative humidity and wind). Bormann et al. (2013) showed that the most dominant climate metrics for annual snow densification rates and spring densities are the winter precipitation sum, followed by temperature and a parameter accounting for melt–refreezing events. In the same manner as precipitation sum SWE is highly correlated with HS. Therefore, in practical applications on regional to global scales (Jonas et al. 2009; Sturm et al. 2010) HS was the main input variable for modelling ρ and SWE. Both approaches consider the time/temperature dependency of densification. Jonas et al. (2009) pooled long-term depth–density data from the Swiss Alps in a monthly resolution and applied three different elevation zones for the characterization of ρ. Sturm et al. (2010) introduced the day of the year to account for the proceeding densification over time within snow climate regions on a global scale. Simultaneously to the present study, the Jonas model was improved regarding the temporal resolution by applying automated SWE data in daily resolution (McCreight & Small 2014) and the Sturm model was improved by the use of additional weather data (Bruland et al. 2015).

Moreover, linear depth–density models were recently applied to convert spatially distributed measurements for basin- or regional-scale studies. Such studies benefit from the finding that the spatial and temporal variability of ρ is lower than the variability of HS and SWE (Elder et al. 1998; Mizukami & Perica 2008). Correlations between HS and ρ were found to be statistically significant but are typically quite low as antecedent meteorological conditions (e.g., large amounts of new snow) can lead to varying snow conditions in the catchment (López-Moreno et al. 2013). Nevertheless, the depth–density correlation was shown to be accurate for the estimation of mean catchment snow conditions (Lundberg et al. 2006) and for end-of-winter snow conditions on glaciers (Jansson 1999). Jörg-Hess et al. (2014) applied a statistical depth–density model for the development of a gridded SWE climatology of Switzerland. Schöber et al. (2014) used a linear depth–density model (derived from snow course data) for the conversion of spatial patterns of HS (obtained from Lidar measurements) into SWE patterns of a glacierized catchment in Tyrol. These snow course data from high alpine areas (2,000–3,400 m a.s.l.) are also used in the present study to complete a rich snow dataset consisting of manual measurements from the permanently settled areas of Tyrol (500–2,000 m a.s.l.). The entire dataset consists of 18,048 HS, SWE and ρ records (mainly observed at weekly intervals) distributed over the entire investigation area allowing a systematic analysis of bulk density and SWE in Tyrol.

This rich dataset is further used to set up a new easy-to-use snow density model (1) for the whole region of Tyrol with consideration of the complete winter season, and (2) a specific model is developed for glacierized catchments with a focus on the peak snow accumulation period. Whereas Jonas et al. (2009) used a monthly clustering of data (bi-weekly resolution), the present study introduces a daily clustering of the data (weekly resolution) applying a moving window approach to overcome discontinuities and jumps at the end of months. Additional model improvements, e.g., the consideration of new snow, are tested for specific locations selected based on data availability. Furthermore, the question arose if it is possible to transfer such a model between different alpine regions. Therefore, the model developed for Tyrol was compared to the model by Jonas et al. (2009) based on Swiss data, considering an application without parameter adjustment to the regional datasets. A second spatial validation was made using the model developed for glacierized catchments in Tyrol by applying it to snow data from Italian and Swiss glaciers.

In Tyrol, manual snow measurements are available from the hydrological year 1952/53. HS and SWE are sampled weekly by staff of the Tyrolean Hydrological Service in the close proximity of meteorological stations. These ‘station data’ are the basis for a regional-scale application accounting for an area of more than 12,000 km² (Figure 1). For measuring SWE, snow pits are dug to ensure that no snow is lost when the core is extracted from the pack. Additionally, data from a snow pillow at Kühtai (1,930 m a.s.l.) are used (Kirnbauer & Blöschl 1990) to densify the data pool at higher elevations. For consistency with the other station data, only one measured pair of HS and SWE from the snow pillow per week (randomly picked) in the period 1993–2010 is added to the calibration dataset. Overall, data from 35 stations were compiled for the period of 1952 to 2010. The mean observation length is 33 years and the data cover mainly the months from November to May. The stations are located at elevations between 528 and 2,036 m a.s.l., covering approximately 60% of the area of the province of Tyrol (Figure 1). Following Jonas et al. (2009), quality checks were carried out using ρ computed from HS and SWE, with data rejected where the computed ρ was outside the range of 50–600 kg m⁻³. Table 1 provides information about the snow data, the providers and the application of the data for calibration or validation.

For six of the Tyrolean stations with measurement periods longer than 30 years, series of daily new snow (HN) are also available (Hinterriss, Jochberg, Ginzling, Boden, Obernberg and Felbertauern – see Figure 1) and used for model refinements. Additional meteorological data from Kühtai, such as wind speed (1993–1999) and air temperature (1993–2010) are the basis for the analysis of bulk density during snowfall.

Additionally, snow data from different field surveys in Tyrol were collected. The majority (93%) were measured above 2,000 m a.s.l. (above sea level) and are referred to as ‘snow course data’. This dataset was obtained in different watersheds of the Ötztal, a partly glacierized sub-region of Tyrol (Figure 1). In total, the dataset consists of 762 snow pit measurements including records of HS, ρ and SWE. Although the measurements are provided by different institutions, no differences in the variability of the data are detectable. All measurements are obtained by the gravimetric method. Roughly two-thirds of the dataset (provided by the Tyrolean hydropower company – TIWAG) were sampled between the 1970s and 1990s in basins of hydropower facilities. What is remarkable about this dataset is that between 20 and 50 snow pits were dug within short periods (1–3 days) providing a temporal snapshot of the spatial distribution of the snow cover at the watershed scale. Unfortunately, the positions of the measurements are known only approximately and the only available topographic information is the elevation of the snow pits. For the glacierized zones (Figure 1), the data are extended by measurements from the following glaciers: Vernagtferner (V) (Escher-Vetter et al. 2009), Hintereisferner (H), Gepatschferner (G) and Kesselwandferner (K) (Helfricht et al. 2014). An additional 1,080 HS soundings from that region are used for analysing spatial snow cover characteristics on the watershed scale. The larger part of this dataset (measured from January to April) is used for validations of the regional-scale models derived from the station data.

Within the snow course dataset, 256 snow pits were measured between mid-April and mid-May; the aim of these measurements was to observe the peak snow accumulation and to determine the glaciological winter mass balance. This dataset, referred to as ‘snow course (WB) data’ is the basis for the development of a statistical density model. Spatial validation is performed with additional measurements (n = 55) from glaciers in Switzerland (Huss & Bauder 2009; Huss et al. 2013) and South Tyrol/Italy (Hydrological Service Bolzano 2006–2010). An overview of the different datasets is provided in Table 1.

Accuracy thresholds were defined regarding the uncertainties of the automatic sensors from the Kühtai station. For this evaluation, accuracies in the order of 5 cm (HS) and 10 mm (SWE) are used. If HS and SWE are low (≤ the measurement accuracy), calculated ρ values are often unrealistic and are discarded. Additionally, ΔHS values lower than 1 cm are excluded from the analysis. From all records in the calibration period 1993–1999, 69% show SWE > 0 mm; from this ‘snow covered’ dataset, approximately 85% remain when the accuracy cut-off level of 10 mm is applied. However, only 17% of the corrected SWE dataset (i.e., 10% of all records) show ΔHS > 1 cm.

Validations of ρ_sim are performed for seven stations for which long-term ρ and HN data are available. This dataset is used to describe the influence of the additional parameter Δρ which is calculated applying either HN or ΔHS. This allows quantifying possible model errors related to the uncertainties of ΔHS. To quantify the approach's transferability within the Alps, the model of Jonas et al. (2009) is included in the modelling of these stations. For comparison with the Swiss model parameters, the daily-based model is evaluated a second time using only the 15th of each month, thus achieving a resolution equivalent to the Swiss model.

For the density model of glacierized catchments, again a linear regression approach comparable to Equation (2) is used with parameters k and d derived as ordinary least square fits. The snow course (WB) data cover a sub-area of snow region 3 (Figure 1) and focus on the peak snow accumulation. Thus, no separate OS_reg calculation applies. Furthermore, most of the data originate from altitudes above 2,000 m and are treated primarily as being located in one high alpine elevation zone. This high alpine area can only be accessed in fair weather conditions and there is a high risk of avalanches when a large amount of new snow accumulates. Thus, no model refinements regarding new snow are applied since the dataset is not influenced by considerable accumulations of new snow. The snow courses of 2011, which contain an exceptional amount of records, are used for validation. Additionally, a spatial validation is made using glacier data from Switzerland and Italy (Table 1) for different years.

The coefficient of variation (CV; standard deviation (σ) divided by the data's mean) is used as a measure for the spatial and temporal variability of the data. The correlation between different data are evaluated by means of Pearson's correlation coefficient r, for which a significance test (p < 0.05) is applied. Mean absolute errors (MAE) for ρ and SWE are calculated in order to check the quality of the model results. Additionally, relative mean absolute errors (MAE_rel) are calculated according to the standardized mean absolute error presented by López-Moreno et al. (2013), which is computed as the MAE divided by the mean of all respective measurements.

The regression models (fitted to the Tyrolean station data) are validated using the remaining SWE data from the Kühtai snow pillow (1,930 m a.s.l.), which is not used for calibration. The snow course data are used as a completely independent validation set (Table 1). For the SWE of the Kühtai snow pillow, a MAE of 27.6 kg m⁻² and a MAE_rel of 13.25% is calculated. For the snow course data of the Ötztal Alps, a MAE of 66.5 kg m⁻² and a MAE_rel of 13.9% is calculated. Coefficients of determination (SWE_obs vs. SWE_sim) are considerable for both validation datasets: R² is 0.95 for the snow pillow and 0.88 for the snow course data (data not shown). However, all of these records were measured in region 3 above 1,400 m a.s.l. Thus, validation of the regression models with more than 4,000 snow datasets, which are not included in the model fitting yields comparable or even better results than for the ‘calibration’ data of the upper elevation zone of the largest snow climate region 3.

Overall improved or equal predictive skills are observed when the model is calibrated with regional data. Comparing the models ‘Figure 11(a),Jonas et al. (2009)’ and ‘Figure 11(b) Tyrol monthly’ shows improvements at three stations (Jochberg, Ginzling and Küthai) and an increase of MAE_rel at the stations Boden and Hinterriß. When changing to a daily resolution, the predictions at the stations Jochberg und Ginzling are again improved. An exception is the station Küthai, which exhibits an increase of the MAE_rel. Further adaptations (OS_regand Δρ) show slightly improved or equal skills for the majority of the stations. Results of the transferred model ‘Figure 11(a) Jonas et al. (2009)’ are equally good or slightly better for the stations Obernberg, Boden and Hinterriß than the Tyrolean models, supporting the transferability of the approach.

The winter mass balance subset of the snow course data are used to characterize peak snow accumulation conditions (mid-April to mid-May) in glacierized regions of the Ötztal, Tyrol. In this area, no regularly observed SWE and ρ data are available. The mean HS of the snow course (WB) data is 2 m (σ = 0.75 m), the mean SWE value is 620 mm (σ = 320 mm) with a mean bulk density of ρ = 370 kg m⁻³ (σ = 62 kg m⁻³). The HS records of this dataset range between 0.23 and 3.7 m with SWE ranging from 101 to 1,485 mm. Bulk densities vary between 232 and 538 kg m⁻³.

In addition, the WB model is spatially validated using measurements obtained at the Langenferner and Weißbrunnferner glaciers in South Tyrol, Italy (Hydrological Service Bolzano 2006–2010), and the following Swiss glaciers: Findelen (years 2012–2013), Plainemorte (2011–2013), Tsanfleuron (2013), Claridenfirn (nine measurements between 1993 and 2008) and Silvretta (6 years between 1960 and 1980) (Huss & Bauder 2009; Huss et al. 2013). Missing years in the long-term snow stake series on the Claridenfirn and Silvretta glaciers were measured outside the defined temporal range of the WB model (mid-April to mid-May). Based on the entire validation data from Italy and Switzerland, a MAE = 86.7 mm and a MAE_rel = 7% are calculated for SWE. As an example for the Swiss validation data, the snow courses on Findelen glacier for the years 2012 and 2013 are plotted in Figures 14(b) and 14(d). The satisfying results of the Tyrolean WB model regarding SWE (MAE_rel = 8.3% in 2012 and 7.3% in 2013) show that the model is suitable to be transferred for estimations of ρ and SWE at the end of the accumulation.

The presented mean density–time curves of the Tyrolean region cover data from the period 1952–2010. With regard to the temporal distribution of the number of records (n), 80% of the data have been collected since the late 1970s, with an almost constant number of records per year collected until 2010. With a mean station observation length of 33 years, the results may be interpreted in terms of climatological mean values of SWE and ρ. The observed higher values of SWE and HS in the mid-1960s might correlate to the higher mean station elevation at that time. Mean daily densification rates (0.6–1.9 kg m⁻³ d⁻¹) extracted from the long-term mean density curves of the Tyrolean station data are in the range of the alpine snow class of a density dataset with an almost global spatial coverage (Bormann et al. 2013). However, a shift to an increased densification rate, as found by Bormann et al. (2013) or Mizukami & Perica (2008), is visible in Tyrol only in the elevation zone between 1,400 and 2,000 m a.s.l.

In contrast to the Swiss model by Jonas et al. (2009), the regression gradient k is occasionally negative in the Tyrolean model. Below 1,400 m a.s.l., almost the complete season faces negative gradients k. Negative correlations may be caused by a large amount of new snow in the Tyrolean data (Table 4) and differences in the sign of k between Tyrol and Switzerland may be attributable to the different sampling scheme (Tyrol = weekly, Switzerland = bi-weekly). At timescales of a few days, HS and ρ are negatively correlated due to new snow which decreases the bulk density. Afterwards, when HS decreases due to settling, ρ increases, which leads to negative correlations as well. This is illustrated in Figure 10. Melt implies a negative correlation between HS and ρ at several days to several weeks. A positive correlation between HS and ρ is the case at a monthly scale during the accumulation period (McCreight & Small 2014). In the long-term mean values of Tyrol, mean HS typically increases from November to March along with a slight increase of ρ (i.e., positive correlation). However, the gradient k is mostly negative during that time. As shown in Figure 8, the negative correlations are significant from November to March. The plateau of the mean HS (0.9 m) from March to May (above 1,400 m a.s.l.) indicates a weak correlation with an ongoing increase of ρ (Figure 6). Still, with only a few exceptions, the correlation is significant and during that time the approach involves the lowest uncertainties (Figure 9). The high resolution data used by McCreight & Small (2014) enabled the predictor to be derived separately for individual years in a temporal window of several days. Differences may occur due to the lower resolution of the Tyrolean data. In contrast, in the same temporal window the predictors represent mean values, aggregating data from all available years. According to McCreight & Small (2014), a 30 day timescale faces mainly negative correlations between HS and ρ. This is only partly the case in the present study which may be due to the different data resolution and the way of deriving the predictors as multi-year mean values.

The daily representation of the depth/density–time relation in the present study is an improvement compared to the original model of Jonas et al. (2009), overcoming data jumps when changing parameters between months. However, the comparison between the Tyrolean models in daily versus monthly resolution reveals only a slight improvement on average. For two stations, it improves the results distinctly (Figure 11). At the Kühtai snow pillow station, however, the shift to the daily resolution decreases the performance measures. The analysis of this behaviour revealed that the data of this station (randomly picked data in weekly resolution) include an especially large number of November and May records where the confidence intervals of the parameters k and d indicate uncertain daily values. Hence, the daily model resolution may include a performance reduction in periods where the confidence interval of parameters becomes too large.

The present study also applies a regional offset parameter (OS_reg; comparable to Jonas et al. 2009) to account for regional characteristics such as a horizontal precipitation gradient. However, the overall model performance is only partly improved by that. The change in the predictive performance appears site specific (Figure 11).

Regarding the possible negative correlation on a daily timescale (Figure 10), a simple parameterization of the decrease of ρ during snowfall is introduced. The bulk density model is extended by an additional offset Δ ρ as a function of the new snow height. The applied regression is based on data from the Kühtai snow pillow. The advanced method shows improvements at most stations with HN records available (Figure 11). However, a limited transferability of this parameterization might be the case as the largest improvement was found for the Kühtai station itself. Δρ was not sensitive to the different new snow data inputs and led to almost similar results when using either manually observed HN or ΔHS (extracted from automated HS). Thus, the approach could be improved in the future by using other automatic sensor sites to derive a regional mean relation for Δρ.

At elevations above 2,000 m a.s.l., mean values of SWE and ρ were evaluated for the peak snow accumulation period at the end of winter, with 15 years of data available during the period 1977–2011. That dataset alone does not make it clear whether the data cover the full range of possible snow accumulation (minimum and maximum SWE) during the observation period. An additional data source is available in the form of three cumulative precipitation gauges in the area which were analysed by Helfricht et al. (2014). A comparison of the winter precipitation statistics with the peak snow accumulation data of the Ötztal Alps collected here reveals that the extremes of winter precipitation within the period 1983–2012 are included in the SWE measurement series (minimum in 1996 and maximum in 2009; Figure 12). Moreover, the mean SWE of 620 mm corresponds well to the long-term mean winter precipitation (1983–2012) of 632 mm.

In the glacierized areas of Tyrol, site elevation (Z_m) as an additional predictor in addition to HS did not improve the SWE estimates considerably. Correlations between ρ and elevation are typically not significant in Tyrol. This finding is similar to results of other studies, e.g., Braun (1984) or López-Moreno et al. (2013). However, in years with large accumulations of SWE (1999, 2005, 2006, 2008 and 2009), the use of Z_m seems to increase the model performance (compare Figure 12 and Figure 15). In contrast, in the winters of 2007 and 2011, which were characterized by large positive deviations from the long-term mean temperature (ZAMG 2014), the model without Z_m is significantly better. Thus, in the current model set-up it is not possible to recommend the use of Z_m for estimating density in the glacierized catchments. In terms of a simpler model the glacierized catchments are seen as a third elevation zone (>2,000 m a.s.l.) in the Tyrolean region. However, site-specific Z_m may be useful in combination with additional weather data as predictors of ρ (Bruland et al. 2015).

The accuracy of SWE measurement generally depends on snowpack characteristics, the ground surface, measurement tools and techniques and the human factor. Manual SWE measurements are expected to have an accuracy in the order of ±10% (Stuefer et al. 2013). The within-site variability of ρ was found to be in the order of 10–20% in Switzerland (Jonas et al. 2009).

With respect to these findings, the results of the density model for glacierized regions are encouraging; a mean relative error of 12.7% with more than half of the annual data having relative errors below 10% (Figure 15) is close to the measurement accuracy of SWE and clearly at the lower bound of the within-site variability. The model benefits from the especially low spatial variability of snow accumulations on glaciers (Schöber et al. 2014) and the low year-to-year variability of ρ at the end of the accumulation season (Figure 12(f)). López-Moreno et al. (2013) found a decreased uncertainty of ρ and SWE estimates (simple function of HS) at the basin scale for ripened snow conditions in April with relative errors in the order of 15%. A more complex model (including additional weather data) which was used to simulate ρ of numerous snow courses resulted in an average error of 9.4% (Bruland et al. 2015).

At the regional scale of Tyrol, model uncertainties are larger. Model errors are about 10% higher for the lower elevation zone (Figure 9) where mean HS is generally lower. Especially in November, a high MAE_rel is found, attributable to the large variability of ρ for shallow snow covers (Jonas et al. 2009; Sturm et al. 2010). This is emphasized by the increased confidence intervals of k and d in November and April/May (Figures 7 and 8). Still, the model performance of the regional scale model is encouraging in the context of the literature. For instance, Rohrer & Braun (1994) applied a model which is forced by daily time series of HN at seven SWE stations in Switzerland. Mean absolute errors of SWE ranged from 25 to 91.2 mm. For the full model of the present study (Figure 11(e), Tyrol daily + OS_REG + Δρ), MAEs range from 16.3 to 35.7 mm at different stations. Egli et al. (2009) used the Jonas et al. (2009) approach to model bi-weekly SWE measurements of a Swiss snow research site with a root mean squared error (RMSE) of 48 mm. Results of the model in Figure 11(e) are equivalent to a RMSE range of 22.1 mm at the Ginzling station and up to 48.9 mm at the Felbertauern station. The model of McCreight & Small (2014) – involving time series in daily resolution – resulted in RMSE for SWE being approximately 63 mm.

The transfer of regression models calibrated for Swiss snow data to simulate SWE at Tyrolean stations led to almost similar MAE values as achieved by the regionally fitted Tyrolean model in monthly resolution. As discussed above, the taking into account of locally fitted model extensions (OS_reg, Δρ) leads to only minor improvements. The transferability of the WB model from glacierized sites in Austria to Switzerland and Italy is found to be reliable without any model extension. This emphasizes the applicability of the presented bulk density models for applications where SWE data are more advantageous than HS, such as the calculation of the glaciological winter mass balance or the calibration of a snow-hydrological model (Schöber et al. 2014).

Long-term mean density–time curves are calculated for Tyrol based on data of the period 1952–2010. Using this comprehensive dataset, an easy-to-use method for estimating snow bulk density and SWE is developed. Daily regression parameters for elevations below and above 1,400 m a.s.l. are compiled (see Figures 7 and 8) and can be used to simulate snow bulk density based on observed snow depth. The presented approach denotes an improvement regarding the temporal resolution especially in the middle of the winter season when confidence intervals of the daily regression parameters are narrow. Model extensions taking into account regional snow characteristics or the impact of new snow improve the model performance locally. Considering new snow is especially important for a model operating on short timescales. However, most reliable model results occur for peak snow accumulation characteristics in spring in glacierized catchments with results close to the measurement accuracy (10%). The model for glacierized catchments is not significantly improved by the application of elevation as additional predictor. It is demonstrated that the simple depth–density models are transferable to other, comparable regions in the Alps. Overall, the developed models are beneficial for (1) decreasing manpower for field work and (2) providing SWE for hydrological and glaciological studies where only HS data are available.

This work has been undertaken in the context of the alpS Comet Project ‘H05 HoPI II – Enhancing Preparedness through Improved Flood Forecasting’ and the Austrian Academy of Science Project ‘HyPCC embedded in the International Strategy for Disaster Reduction (ISDR) programme. The authors would like to thank the Federal State of Tyrol, the Austrian Research Promotion Agency (FFG) and TIWAG – Tiroler Wasserkraft AG who support the project. Furthermore, the authors would like to thank the Hydrographic Service of Tyrol (Wolfgang Gattermayr, Georg Raffeiner), TIWAG (Helmut Schönlaub), the Commission for Glaciology of the Bavarian Academy of Sciences and Humanities (Heidi Escher-Vetter, Ludwig Braun), the Institute for Interdisciplinary Mountain Research (Andrea Fischer), the Institute of Meteorology and Geophysics of the University of Innsbruck (Stephan Galos), and the Unit of Geography of the University of Fribourg (Matthias Huss) for supporting the study with the large number of snow measurements. The authors would like to express their thanks to the two anonymous reviewers for their helpful comments.