Snow density is an important measure in hydrology used to convert snow depth to the snow water equivalent (SWE). A model developed by Sturm, Tara and Liston predicts the snow density by using snow depth, the snow age and a snow class defined by the location. In this work this model is extended to include location and seasonal weather-specific variables. The model is named Weather Snow Density Model (Weather SDM). A Bayesian framework is chosen, and the model is fitted to and tested for 4,040 Norwegian snow depth and densities measurements between 1998 and 2011. The final model improved the snow density predictions for the Norwegian data compared to the model of Sturm by up to 50%. Further, the Weather SDM is extended to utilize local year-specific snow density observations (Weather&ObsDensity SDM). This reduced the prediction error an additional 16%, indicating a significant improvement when utilizing information provided by annual snow density measurements.

INTRODUCTION

Approximately 60% of the precipitation in Norway falls as snow, and the snow storage is essential for flood and avalanche forecasting and warning, climate research and hydropower production. In hydropower production snow represents energy storage, and exact estimates of the snow storage, or snow water equivalent (SWE), are important for power generation scheduling. The value of the Norwegian snow storage is approximately 4 billion US dollars in an average year, and hence small errors in snow storage estimates can represent large values. SWE is determined by snow depth and snow density. Snow depth typically has higher variability and is most strongly correlated to SWE, but is also far easier and less time-consuming to observe than snow density. SWE in Norway has traditionally been observed by measuring snow depths along snow courses selected to describe the topographical variability in the catchment. At least three density samples are recommended for each course, giving between one and two samples per 20 depth measurements. Even though this is more than reported by Sturm et al. (2010), the uncertainty related to snow storage estimates are high. The amount of snow depth data possible to collect has increased significantly over the 15 last years by use of snow radar (Bruland et al. 2001; Lundberg et al. 2006, 2010) which now are in operational use among Norwegian hydropower producers. If also LiDAR (Sturm et al. 2010) is taken into operational use among hydropower producers, the amount of snow depth observations will further increase. With better description of snow depth variability, snow density estimate errors become more predominant in the estimation of the areal SWE. There is thus an increased focus on improving snow density estimates.

The density of snow is a result of the climate during the snow accumulation period and the overlaying weight of snow. Different models have been developed to calculate the snow density. Steppuhn (1976) and Lundberg et al. (2010) indicated that there is correlation between snow depth and density, but Steppuhn found only a weak correlation at snow depths deeper than 85 cm. Elder et al. (1998) found that a linear model based on net solar radiation, elevation and slope could explain 70% of the variation in observed snow densities. Jonas et al. (2009) was able to estimate SWE with an accuracy equivalent to the variability of repeated SWE measurements at one site in the Swiss Alps by using a snow density model (SDM) based on information about season, snow depth, site altitude, and site location. CROCUS (Brun et al. 1989), Alpine 3D (Lehning et al. 2006) and SnowTran-3D (Liston & Sturm 1998; Liston & Elder 2006) are all models describing snow density as a function of snow depth, snow age and snow metamorphism. However, with the exception of SnowTran-3D, these models are all complex with an extensive need for detailed information and thus computationally intensive to run over large domains. Sturm et al. (2010) introduced a model that estimates the local snow bulk density using snow depth, snow age and snow climate classes from Sturm et al. (1995). They found that the relative errors in SWE using this density model were of the same magnitude as SWE that would be encountered at a single site due to local variability. Thus, the simplicity of the method, as well as the low requirements of input data, makes the method presented in Sturm et al. (2010) interesting for testing for operational applications for hydropower producers in Norway.

In Norway there are large differences in local climate and topography, and it is questionable whether the climate classes used in Sturm's model cover the regional and temporal variations found here. Statkraft is the largest hydropower producer in Norway and has an extensive data set of snow surveys including more than 37,000 snow density observations from 982 locations in Norway. In this study we test how well Sturm's model, with the defined climate classes, reproduces these observations. We further introduce and test new SDMs based on Sturm's model: Weather SDM. In Weather SDM, local and season-specific observed weather quantities as well as other local explanatory variables are included in the SDM. The Weather SDM, like Sturms's model, can be used for areas without snow density observations. We introduce also a version of Weather SDM where we can include year and area-specific snow density observations in the model. A Bayesian framework is chosen for modeling and inference. For a thorough introduction to Bayesian statistics relevant to this study, see e.g. Gamerman & Lopes (2006). Parameters are then considered random variables, and are given prior distributions. These distributions are updated using the data through the model for the observations given the parameters, known as the likelihood. This results in posterior probability distributions for the parameters, and it is common to report properties of posterior distribution such as the posterior mean and the posterior standard deviation. We further test whether the models reproduce the observations satisfactorily, and thus whether they can be used operationally. In this study, a selection of 4,040 of Statkraft's snow density observations from different areas in Norway between 1998 and 2011 are used to fit and test the models. This work is based on a Master thesis (Færevåg 2013), and more details can be found there.

BACKGROUND

Observed snow densities vary normally between 0.1 g/cm3 for new fallen snow to 0.5 g/cm3 and sometimes higher, for wind-packed, icy or wet snow. Sturm et al. (2010) developed a model that estimates the snow bulk density based on data from the USA, Canada and Switzerland. It takes snow depth, snow age and snow class as input variables. The snow class is found by a classification system for seasonal snow cover proposed in Sturm et al. (1995). It has six classes, where each class is defined in terms of physical characteristics of the snow and the snow layers. The classes are derived for each location by using wind, precipitation and air temperature climate variables found from a set of weather stations. The snow classes in Scandinavia are shown in Figure 1 and in Norway tundra and maritime snow are dominant.

Figure 1

Snow class distribution in Norway. Source:Liston & Hiemstra (2011).

Figure 1

Snow class distribution in Norway. Source:Liston & Hiemstra (2011).

Sturm et al. (2010) model snow density observation Y as beta distributed with an expected value that is a function of snow depth (h) in cm, the day of year (DOY) for the measurement, and the snow class parameters k1, k2, ρ0 and ρmax 
formula
1
where 0 < ρ0ρmax < 1. k1 and k2 are the densification parameters for depth and DOY, ρmax is maximum bulk density and ρ0 is the initial density of the individual snow layer. The snow class parameters are given in Table 1. The snow season begins early in October, and DOY represents the effect of snow ageing and the number of days in the winter season.
Table 1

Model parameters for each snow class (Sturm et al. 2010)

Snow class ρmax ρ0 k1 k2 
Tundra 0.3630 0.2425 0.0029 0.0049 
Maritime 0.5979 0.2578 0.0010 0.0038 
Prairie 0.5940 0.2332 0.0016 0.0031 
Alpine 0.5975 0.2237 0.0012 0.0038 
Taiga 0.2170 0.2170 0.0000 0.0000 
Snow class ρmax ρ0 k1 k2 
Tundra 0.3630 0.2425 0.0029 0.0049 
Maritime 0.5979 0.2578 0.0010 0.0038 
Prairie 0.5940 0.2332 0.0016 0.0031 
Alpine 0.5975 0.2237 0.0012 0.0038 
Taiga 0.2170 0.2170 0.0000 0.0000 
The beta distribution has two parameters, and using a parameterization with the expected value as one of the parameters, the other parameter is often called the precision and is modeled by Sturm et al. (2010) as 
formula
2
where and .

A Norwegian data set with 37,000 snow density observations from the years 1930 to 2012 is used for comparing snow density estimation methods developed in this study with traditional estimation methods. The overall characteristics of the snow density and depth in the data set are given in Figure 2. The highest observed density is 0.656 g/cm3 and the lowest is 0.052 g/cm3. The average is 0.33 g/cm3 with a standard deviation of 0.07 g/cm3. 95% of the observations are between 0.20 and 0.48 g/cm3 (Figure 2(a)). Figure 2(b) shows that both snow depth and snow density increase through the season and data analysis shows that 29% of the observed variance in snow densities in this data set can be explained by snow age and depth.

Figure 2

(a) Snow density distribution and (b) average snow densities during the winter.

Figure 2

(a) Snow density distribution and (b) average snow densities during the winter.

Calculation of SWE for a catchment is normally based on the average snow density or a density found by using the snow depth/snow density relation from the observations at the time of the measurement campaign. In the data set one campaign has on average 15 density observations in each catchment. Split sample analysis of this data set shows that the average error when estimating a snow density based on mean values is 9.9%. In a data set with less than 13 observed densities the density is better estimated when using average density than using a linear relation between density and depth found from the observations (Figure 3(d)). Figure 3 further shows that the error has increased over the years (Figure 3(a)) and decreases with snow depth (Figure 3(b)) and number of observations (Figure 3(d)) for both methods. The errors are largest for densities below 0.1 g/cm3 and above 0.5 g/cm3 (Figure 3(c)).

Figure 3

Snow density estimation errors when using traditional methods.

Figure 3

Snow density estimation errors when using traditional methods.

Apart from snow depth, the model of Sturm et al. (2010) does not include the variability in the density due to yearly variability in the climate. Using the snow classes to predict the snow density means that as long as snow depths measured on the same date in various years are equal and all other conditions are unchanged, then the estimated density will be equal too. The historical records of density measurements in Norway show that the errors in snow density estimations have increased over the years since 1930 (Figure 3(a)). As long as density observation methodology is unchanged, this indicates changing conditions and thus a need for a methodology that takes into account long-term changes in addition to year-specific conditions.

DATA AND MODELING

Snow data and meteorological data

The analysis, model development and testing in this paper are based on snow depth and snow density observations from 244 locations within 17 different areas (see Figure 4), and covers most prevailing climates in Norway. Table 2 displays the characteristics of the snow data and climate for each of the measurement areas as well as number of locations in each area and the number of snow measurements used in each area. The snow measurements are from the period February to May with the majority between mid-February and mid-April. The year-to-year variation in snow densities and concurrent snow depths for Area F can be seen in Figure 5. The measurement in this area is done in the same period each year, and the between-year variability seen here cannot be explained by snow depth and climate characteristics only. Therefore, models that also account for year-to-year variations in weather are proposed.

Table 2

Number of sub-locations, observed mean snow depth, snow density, snow water equivalent, total number of measurements (n), and information about the climate for the period November–April for each measurement area

       Temp
 
Wind
 
 
ID Area No. of locations Depth (cm) Density (g/cm3SWE (mm) n Mean Max Min Mean Max Precip total 
Adamselv 74.1 0.345 26.4 42 −5.7 15.0 −26.0 4.5 21.1 49 
Alta 51.7 0.266 13.4 18 −8.4 11.3 −38.8 2.2 16.7 101 
Aura/Grytten 24/3 119.3 0.349 43.2 439 −0.9 20.0 −30.0 2.5 19.6 437 
Tysso-Folgefonn 12 219.5 0.390 87.7 168 −2.1 15.0 −23.7 2.8 27.0 912 
Innset 32 90.5 0.282 27.1 556 −5.8 15.0 −40.0 2.2 13.9 366 
Jostedalen 190.3 0.327 64.1 39 −2.6 12.2 −25.3 1.6 22.0 631 
Kobbelv 162.5 0.414 68.8 188 −3.8 12.8 −30.9 5.1 27.7 514 
Nea-Nidelv 10 108.8 0.338 36.6 144 −4.1 16.9 −34.5 2.7 19.9 273 
Nore 19 81.8 0.295 25.2 306 −6.4 16.2 −39.2 2.8 26.1 276 
10 Rana 107.4 0.315 35.3 254 −4.4 14.5 −46.7 1.2 17.7 457 
11 Røssåga 15 103.3 0.353 37.7 273 −3.0 15.0 −26.9 1.6 9.8 1145 
12 Skjomen 24 123.6 0.34 43.4 324 −3.4 18.9 −25.7 3.0 50.0 273 
13 Svorka/Trollheim 3/9 126.5 0.376 50.3 243 −1.8 18.1 −34.4 1.6 50.0 599 
14 Tokke 40 99.9 0.306 31.9 535 −3.4 17.6 −34.9 1.8 24.6 604 
15 Ulla-Førre 10 138.6 0.398 57.6 243 −0.7 21.7 −26.5 3.6 29.9 1143 
16 Vik/Høyanger 5/2 184.8 0.384 75.4 191 −3.4 14.4 −45.0 4.4 23.2 668 
17 Sira-Kvina 13 82.3 0.339 28.8 77 −3.0 16.2 −29.0 2.2 11.5 658 
Total  244 118.2 0.337 42.5 4040 −3.7 21.7 −46.7 2.7 50.0 536 
       Temp
 
Wind
 
 
ID Area No. of locations Depth (cm) Density (g/cm3SWE (mm) n Mean Max Min Mean Max Precip total 
Adamselv 74.1 0.345 26.4 42 −5.7 15.0 −26.0 4.5 21.1 49 
Alta 51.7 0.266 13.4 18 −8.4 11.3 −38.8 2.2 16.7 101 
Aura/Grytten 24/3 119.3 0.349 43.2 439 −0.9 20.0 −30.0 2.5 19.6 437 
Tysso-Folgefonn 12 219.5 0.390 87.7 168 −2.1 15.0 −23.7 2.8 27.0 912 
Innset 32 90.5 0.282 27.1 556 −5.8 15.0 −40.0 2.2 13.9 366 
Jostedalen 190.3 0.327 64.1 39 −2.6 12.2 −25.3 1.6 22.0 631 
Kobbelv 162.5 0.414 68.8 188 −3.8 12.8 −30.9 5.1 27.7 514 
Nea-Nidelv 10 108.8 0.338 36.6 144 −4.1 16.9 −34.5 2.7 19.9 273 
Nore 19 81.8 0.295 25.2 306 −6.4 16.2 −39.2 2.8 26.1 276 
10 Rana 107.4 0.315 35.3 254 −4.4 14.5 −46.7 1.2 17.7 457 
11 Røssåga 15 103.3 0.353 37.7 273 −3.0 15.0 −26.9 1.6 9.8 1145 
12 Skjomen 24 123.6 0.34 43.4 324 −3.4 18.9 −25.7 3.0 50.0 273 
13 Svorka/Trollheim 3/9 126.5 0.376 50.3 243 −1.8 18.1 −34.4 1.6 50.0 599 
14 Tokke 40 99.9 0.306 31.9 535 −3.4 17.6 −34.9 1.8 24.6 604 
15 Ulla-Førre 10 138.6 0.398 57.6 243 −0.7 21.7 −26.5 3.6 29.9 1143 
16 Vik/Høyanger 5/2 184.8 0.384 75.4 191 −3.4 14.4 −45.0 4.4 23.2 668 
17 Sira-Kvina 13 82.3 0.339 28.8 77 −3.0 16.2 −29.0 2.2 11.5 658 
Total  244 118.2 0.337 42.5 4040 −3.7 21.7 −46.7 2.7 50.0 536 
Figure 4

Locations of the study 17 areas grouped in six main groups (A–F). Circles mark positions of meteorological stations.

Figure 4

Locations of the study 17 areas grouped in six main groups (A–F). Circles mark positions of meteorological stations.

Figure 5

Observed snow density in the Ulla-Førre area in the period 1983–2012. Corresponding snow depths are shown as bars at the bottom of the plot.

Figure 5

Observed snow density in the Ulla-Førre area in the period 1983–2012. Corresponding snow depths are shown as bars at the bottom of the plot.

In the proposed models, weather variables that are based on observed temperature, precipitation and wind are used. This meteorological information is obtained from meteorological stations within the areas and covers both data from the Norwegian Meteorological office (Met.no) and Statkraft's meteorological station network (Figure 4).

Temperature normally decreases with elevation, and location-specific temperature is calculated from the closest weather station observations using Equation (3): 
formula
3
where T is the air temperature to be found for the location of interest at elevation z and T0 and are the air temperature and elevation (masl) of the representative weather station and a is the temperature gradient from a nearby calibrated precipitation–runoff model.

The precipitation measurements are corrected for catch loss using local wind measurement and corrections suggested by Førland et al. (1996).

Construction of weather variables

In addition to the two explanatory variables used in Sturm et al. (2010), snow depth and DOY, seven new explanatory variables are introduced. These variables are denoted x1, x2, … , x9. Most of the proposed explanatory variables are based on weather observations and state variables from nearby calibrated precipitation–runoff models. The rationales behind choosing these explanatory variables, as well as details of how they are derived, are given below.

Several accumulated weather variables are introduced, and these are calculated for each snow density observation. All accumulated weather variables are calculated based on hourly observations from a snow accumulation start time t =A0 to the day of the snow density observations (DOY). Snow accumulation calculations in HBV models established for catchments in the vicinity are used to set the accumulation start time A0. If there was snow falling prior to A0, the model would have found that this has melted away. The snow accumulation start time A0 is thus the most recent date snow accumulation started, i.e. when the observed snow layer started to accumulate. The diurnal snow melt and accumulation have to be calculated on an hourly basis, but can be summarized and represented as daily values for numerical convenience in the further calculations.

Snow depth, snow age and elevation

The snow density increases with age and higher compaction (snow depth). x1 represents the observed snow depth of the measurement, which is the accumulated sum of new and old snow. x2 is the DOY the measurement is carried out, representing the snow age. DOY starts at −92, the 1 October and runs to +181, 30 June. If DOY is 1, it means that the snow depth and density is measured on 1 January. x3 is the elevation, or the height, in meters above sea level (masl) of the location of the measurement.

Positive degrees

In order to include the influence of high temperatures on the snow density, x4 is the accumulated sum of positive degrees.

Wind

Wind speed less than 2 m/s is here assumed to have no effect on the snow density. The variable x5 is a measure of the accumulated sum of wind velocities above 2 m/s while temperature is below freezing point.

Snowfall and wind

x6 is the amount of accumulated precipitation falling as snow when there is wind 
formula
4

Precipitation type

Under which conditions the precipitation is falling can influence the density of the snow. This is quantified by classifying precipitation in three categories and using the ratio of each class as explanatory variable: x7 is ratio of the total precipitation that comes when there is snow, x8 mixed snow and rain and x9 when it is raining 
formula
5
is the total precipitation at the location of each snow density observation. Ts and Tr are threshold temperatures for precipitation as snow and rain, respectively, found from the nearest located HBV models.

Weather snow density model (Weather SDM)

The model of Sturm et al. (2010) is extended by allowing more explanatory variables in the model of the expected snow density. A general form of the model is given by 
formula
6
where xv is the explanatory variable number v introduced in the section Construction of weather variables for each observation and kv the associated densification parameter to be estimated. Model performance will define which variables are to be included in the final model. To complete a Bayesian model specification, prior distributions for all parameters both in Equations (6) and (7) and Equation (2) need to be specified. Following Sturm et al. (2010), uniform prior distributions are used, but the range needs to be wider due to other covariates: , , .

Observed density extended weather SDM (Weather&ObsDensity SDM)

The Weather SDM presented above can be used to estimate snow density only using weather variables and snow depth. It is likely that there exist effects caused by weather and local conditions that Weather SDM does not capture. These effects can probably only be captured by local snow density observations. A reasonable hypothesis is thus that if there are snow density observations available for the specific year and area of interest, a better estimate can be obtained by using these in the model. Weather&ObsDensity SDM, a version of Weather SDM that allows us to utilize snow density observations for a given area and year, is therefore introduced. This is done by extending the Weather SDM in Equation (6) to also include random effects that are unique for each year (m) and area (j). 
formula
7
The j different areas are likely to have different overall responses for each year m. The model accounts for this by including a term that can be interpreted as a correction factor for area j in year m. The parameters are given the same prior distributions as for the Weather SDM, and the random effects are given uniform prior distribution .

Inference and evaluation

It is only possible to analytically derive posterior probability distributions in a few special cases. In most cases, including the models presented above, numerical approximations have to be used. A powerful method for doing Bayesian inference is the Markov Chain Monte Carlo (MCMC) methods (Gilks 2005). A Markov chain is then constructed such that it, after a burn in, provides us with samples from the posterior distribution. These MCMC samples can then be used to estimate, for example, posterior mean for each parameter of interest. All inference is performed by MCMC simulations using the WinBUGS (Lunn et al. 2000) and OpenBUGS (Spiegelhalter et al. 2007) software. Each Markov chain ran for 15,000 iterations, and the first 5,000 were discarded as a burn-in period.

Models' performance is based on their predictive ability of snow density. The mean error, the mean absolute error (MAE), the root mean squared error (RMSE) of the posterior mean estimates and mean continuous ranked probability score (CRPS) for the posterior density distributions are used to evaluate the models. Since a predictive distribution is considered, and CRPS measures both sharpness of the distribution and reliability of the predictions, CRPS will be the preferred score. To make the results easier to compare, MAE and RMSE are also calculated with weights compensating for different numbers of observations between the areas.

If models are fitted and tested on the same data set, a more complex model will always have a better fit. But the better fit might be due to over-fitting; some of the variables should not be in the model, and the better fit is due to fitting random errors. Over-fitting will result in poorer predictive performance when used to new observations. To make the evaluation realistic (cross-)validation schemes are used. Data from 1998 to 2005 are used to train the models (i.e. fit the parameters), and data from 2006 to 2011 are used to test the models (with some adjustments for areas where the complete data set is not available). This means that all parameters for the Weather SDM and all parameters except the random effects for the Weather&ObsDensity SDM are estimated using the observations in the training set (1998–2005). These parameters are applied for testing the models using snow density observations in the test set (year 2006 to 2011). To estimate the random effect for a specific year and area in the test set, snow density observations for this year and area are needed. The parameters for the Weather&ObsDensity SDM are estimated in two different ways, denoted Method-2a and Method-2b. For Method-2a the corresponding Weather SDM is first fitted using the training set. For Method-2b the Weather&ObsDensity model is fitted using the training data set, i.e. both parameters and random effects are estimated. For both methods the parameters obtained are hold fixed when the models are tested against snow density observations for 2006–2011. To estimate random effects for the years and areas in the test set, a cross-validation scheme is used. The observations done in the same year and area are randomly partitioned k = 5 times into a training set and a test set. The random effect is fitted k = 5 times using the training sets, and the model is tested k = 5 times using the snow densities in the corresponding test set. The mean of the scores for the k = 5 training sets are reported. To test how many snow density observations that are useful, the number of observations in the training sets is changed. For each training set size the procedure above is followed.

RESULTS

First, different configurations of the explanatory variables of the Weather SDM are compared to choose which explanatory variables to include. Next, the corresponding model including observed densities are fitted and tested. For comparison, Sturm's model described in the section Background is also used, and corresponding predictive scores given.

Model choice Weather SDM

The explanatory variables and their correlation to snow density and internal correlation are summarized in Table 3. Snow depth has a high correlation to snow density, and therefore snow depth is included in all models. Ten different models with different combinations of the explanatory variables were tested. These are named Model A to Model J, and they are specified in Table 4. These models were chosen based on the estimated k parameters for Model J (the full model), the estimated parameters and predictive performance for Model A (including snow depth and precipitation type ratios) and Model B (including snow depth and accumulated weather variables).

Table 3

Correlation between the nine explanatory variables constructed for the model and between snow density (y) and the covariates

Covariate Description y x1 x2 x3 i4 x5 x6 x7 x8 x9 
x1 Snow depth 0.48 1.00         
x2 Snow age (DOY) 0.39 0.25 1.00        
x3 Elevation (masl) 0.20 0.39 0.00 1.00       
x4 Degree days when Tt >0 °C 0.16 −0.05 0.40 −0.30 1.00      
x5 Wind days when Tt < 0 °C, W > 2 m/s 0.24 0.21 0.42 0.32 −0.06 1.00     
x6 Snowfall when W > 0 m/s 0.29 0.20 0.23 0.34 −0.02 0.53 1.00    
x7 Ratio snow precipitation −0.05 0.08 0.09 0.34 −0.53 0.15 0.09 1.00   
x8 Ratio mixed precipitation 0.07 −0.04 −0.09 −0.17 0.34 −0.02 0.08 −0.75 1.00  
x9 Ratio rain precipitation 0.03 −0.08 0.03 −0.35 0.60 −0.19 −0.19 −0.69 0.33 1.00 
Covariate Description y x1 x2 x3 i4 x5 x6 x7 x8 x9 
x1 Snow depth 0.48 1.00         
x2 Snow age (DOY) 0.39 0.25 1.00        
x3 Elevation (masl) 0.20 0.39 0.00 1.00       
x4 Degree days when Tt >0 °C 0.16 −0.05 0.40 −0.30 1.00      
x5 Wind days when Tt < 0 °C, W > 2 m/s 0.24 0.21 0.42 0.32 −0.06 1.00     
x6 Snowfall when W > 0 m/s 0.29 0.20 0.23 0.34 −0.02 0.53 1.00    
x7 Ratio snow precipitation −0.05 0.08 0.09 0.34 −0.53 0.15 0.09 1.00   
x8 Ratio mixed precipitation 0.07 −0.04 −0.09 −0.17 0.34 −0.02 0.08 −0.75 1.00  
x9 Ratio rain precipitation 0.03 −0.08 0.03 −0.35 0.60 −0.19 −0.19 −0.69 0.33 1.00 
Table 4

Different test models

Model x1 x2 x3 x4 x5 x6 x7 x8 x9 
✓      ✓ ✓ ✓ 
✓   ✓ ✓ ✓    
✓   ✓ ✓  ✓ ✓ ✓ 
✓   ✓ ✓     
✓  ✓ ✓ ✓     
✓  ✓ ✓  ✓    
✓ ✓  ✓ ✓     
✓  ✓ ✓      
✓  ✓ ✓ ✓ ✓    
✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ 
Sturm ✓ ✓        
Model x1 x2 x3 x4 x5 x6 x7 x8 x9 
✓      ✓ ✓ ✓ 
✓   ✓ ✓ ✓    
✓   ✓ ✓  ✓ ✓ ✓ 
✓   ✓ ✓     
✓  ✓ ✓ ✓     
✓  ✓ ✓  ✓    
✓ ✓  ✓ ✓     
✓  ✓ ✓      
✓  ✓ ✓ ✓ ✓    
✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ 
Sturm ✓ ✓        

x1: snow depth, x2: day of year, x3: elevation, x4: accumulated wind, x5: accumulated positive degrees, x6: accumulated snow when there is wind, x7: ratio snow, x8: ratio mixed snow and rain, x9: ratio rain.

In Figure 6, posterior mean snow density estimates from Sturm's model and Model A, Model E and Model G are plotted together with the observed snow density for four locations to give an impression of how the predictions look. The performances of the different selected combinations are presented in Table 5. Model J (the full model) is the only Weather SDM that performs more poorly than Sturm's model. Further, Model C, Model D and Model E have very similar scores, but Model E has a slightly better CRPS score. The only difference between Model D and Model E is that the covariate elevation (x3) is included in Model E. As the CRPS implies a slightly better estimate for Model E, and as the cost of including elevation is also low, this model is recommended. From this point Model E is referred to as the Weather SDM. This model has expected value: 
formula
8
or written in another form 
formula
9
Table 5

Test result of the different models

Model MAE Weighted MAE RMSE Weighted RMSE CRPS 
0.0508 0.0537 0.0607 0.0648 0.03776 
0.0465 0.0491 0.0567 0.0602 0.07472 
0.0460 0.0487 0.0556 0.0590 0.03440 
0.0460 0.0486 0.0555 0.0589 0.03440 
0.0462 0.0488 0.0558 0.0590 0.03418 
0.0493 0.0520 0.0596 0.0632 0.07675 
0.0467 0.0490 0.0561 0.0590 0.03451 
0.0477 0.0504 0.0572 0.0608 0.03538 
0.0465 0.0490 0.0567 0.0600 0.06931 
0.1487 0.1497 0.1577 0.1593 0.08074 
Sturm 0.0617 0.0636 0.0719 0.0743 0.06170 
Model MAE Weighted MAE RMSE Weighted RMSE CRPS 
0.0508 0.0537 0.0607 0.0648 0.03776 
0.0465 0.0491 0.0567 0.0602 0.07472 
0.0460 0.0487 0.0556 0.0590 0.03440 
0.0460 0.0486 0.0555 0.0589 0.03440 
0.0462 0.0488 0.0558 0.0590 0.03418 
0.0493 0.0520 0.0596 0.0632 0.07675 
0.0467 0.0490 0.0561 0.0590 0.03451 
0.0477 0.0504 0.0572 0.0608 0.03538 
0.0465 0.0490 0.0567 0.0600 0.06931 
0.1487 0.1497 0.1577 0.1593 0.08074 
Sturm 0.0617 0.0636 0.0719 0.0743 0.06170 

Framed values indicate the best score for the different evaluation criteria.

MAE: mean absolute error; RMSE: root mean square error; CRPS: continuous ranked probability score.

Both MAE and RMSE scores are also weighed by the number of observations.

Figure 6

Posterior mean estimates and observed density in four different locations for four different models without random effects. (a) Vinjerui: area 14, maritime; (b) Illekleiv: area 4, tundra; (c) Stordalen: area 3, tundra; (d) Stearuvuggi: area 7, tundra.

Figure 6

Posterior mean estimates and observed density in four different locations for four different models without random effects. (a) Vinjerui: area 14, maritime; (b) Illekleiv: area 4, tundra; (c) Stordalen: area 3, tundra; (d) Stearuvuggi: area 7, tundra.

where h is snow depth, z is elevation, T* are accumulated degree days when Tt >0 °C and W* are accumulated wind days when Tt < 0 °C and W > 2 m/s. Compared to the model from Sturm et al. (2010), the DOY is excluded from the model. This variable is relatively highly correlated with the snow density, but also with the accumulated weather variables T* and W*, and including DOY gave worse predictions. An explanation for this is that the effect of ageing is already embedded in the measurement of snow depth and the accumulation of weather variables. The parameter estimates from the MCMC simulations are given in Table 6. The average absolute error of the prediction based on Weather SDM is 0.046 g/cm3 or 9.4% for the test period from 2006 to 2011. Table 7 shows the errors from the different areas. There are three areas which have a remarkably high error: Innset (area 5), Jostedal (area 6) and Nore (area 9). For Nore, the errors are about average for all years except 2006 which has large errors. For Innset and Jostedal the errors are generally large.

Table 6

Posterior mean and standard deviation for model parameters and number of Markov Chain Monte Carlo samples

Variable Mean Std dev Sample size 
ρ0 0.1481 0.0148 10,000 
ρmax 0.4720 0.0128 10,000 
k1 0.00503 0.000548 10,000 
k3 0.00018 0.00463 10,000 
k4 0.00477 0.000682 10,000 
k5 0.00042 0.0000741 10,000 
Variable Mean Std dev Sample size 
ρ0 0.1481 0.0148 10,000 
ρmax 0.4720 0.0128 10,000 
k1 0.00503 0.000548 10,000 
k3 0.00018 0.00463 10,000 
k4 0.00477 0.000682 10,000 
k5 0.00042 0.0000741 10,000 
Table 7

Prediction error in g/cm3 and average absolute relative error for each area in the test period

Area/year 2006 2007 2008 2009 2010 2011 Abs rel error (%) 
– – – – 0.015 0.03 
– 0.046 −0.033 0.022 0.024 0.062 12 
−0.045 0.007 0.034 −0.019 −0.005 
– −0.016 −0.005 −0.027 −0.018 0.009 
−0.065 −0.02 −0.052 −0.049 −0.076 −0.02 18 
−0.026 −0.05 −0.08 −0.087 −0.064 −0.047 19 
0.034 0.064 0.065 0.052 0.025 −0.003 
−0.056 −0.037 −0.025 −0.014 0.004 −0.014 
−0.134 0.018 −0.009 −0.039 −0.041 −0.028 18 
10 −0.017 −0.015 −0.022 0.036 −0.013 −0.027 
11 0.02 0.019 0.002 0.054 0.004 −0.001 
12 – 0.021 −0.029 −0.017 −0.037 −0.026 
13 0.01 0.038 −0.03 0.072 −0.033 −0.009 
14 – – – −0.021 −0.024 −0.023 
15 0.062 0.041 0.027 0.035 0.026 0.052 11 
16 −0.009 0.007 0.003 0.016 0.004 −0.02 
17 – – – – – 0.046 12 
Area/year 2006 2007 2008 2009 2010 2011 Abs rel error (%) 
– – – – 0.015 0.03 
– 0.046 −0.033 0.022 0.024 0.062 12 
−0.045 0.007 0.034 −0.019 −0.005 
– −0.016 −0.005 −0.027 −0.018 0.009 
−0.065 −0.02 −0.052 −0.049 −0.076 −0.02 18 
−0.026 −0.05 −0.08 −0.087 −0.064 −0.047 19 
0.034 0.064 0.065 0.052 0.025 −0.003 
−0.056 −0.037 −0.025 −0.014 0.004 −0.014 
−0.134 0.018 −0.009 −0.039 −0.041 −0.028 18 
10 −0.017 −0.015 −0.022 0.036 −0.013 −0.027 
11 0.02 0.019 0.002 0.054 0.004 −0.001 
12 – 0.021 −0.029 −0.017 −0.037 −0.026 
13 0.01 0.038 −0.03 0.072 −0.033 −0.009 
14 – – – −0.021 −0.024 −0.023 
15 0.062 0.041 0.027 0.035 0.026 0.052 11 
16 −0.009 0.007 0.003 0.016 0.004 −0.02 
17 – – – – – 0.046 12 

Weather&ObsDensity SDM

The next question is whether the performance achieved using the Weather SDM can be improved by utilizing snow density measurements for the current year and area of interest. Further, how a potentially improved predictive performance relates to the number of snow density measurements, and how many snow density measurements are beneficial to do. These questions are answered by fitting the Weather SDMs including density observations (Weather&ObsDensity SDM) with the selected explanatory variables for different number (j) of area and year-specific (training) data.

The model has expected snow density: 
formula
10
where ρ is the snow density for each snow measurement in area j for year m. j ∈ (1, … ,17) and m ∈ (1998, … ,2012), and is the corresponding random effect.

The two different methods for fitting the Weather&ObsDensity SDM are tested against reduction in prediction error when using more snow density measurements from the current year and area. Figure 7 displays the reduction in the MAE as a function of number of annual measurements n = {0,5,10,15,25} in year 1998 (n = 0 corresponds to Weather SDM). Method-2a gave on average 6.5% reduction with 20 observations within an area compared to the Weather SDM (non-year and area-specific snow density observations), while Method-2b already with five density samples gave a reduction of 9.3% increasing to 16% with 20 observations. These results show that: (1) if year and area-specific measurements are to be used, the parameters should be estimated taking that into account, i.e. with Method-2b; and (2) area and year-specific snow density measurements can contribute to better snow density predictions, but only the first 10–15 measurements. This result supports the hypothesis that including area and year-specific snow density measurements improves the model compared to the Weather SDM which do not use (or need) year and area-specific snow density observations. But the improvement stagnates after 10–15 snow density observations.

Figure 7

Error reduction in absolute error by collection of n = {0,5,10,15,20,25} number of year-specific measurements. Data from year 1998.

Figure 7

Error reduction in absolute error by collection of n = {0,5,10,15,20,25} number of year-specific measurements. Data from year 1998.

DISCUSSION AND CONCLUSION

Snow densities predicted from a model presented by Sturm et al. (2010) were compared to a Norwegian data set. By using the snow classification data set in Sturm et al. (1995) most of the 244 locations in the Norwegian data set were classified as maritime (58) and tundra (175) snow class, and a few were classified as taiga snow. In most cases, the model underestimated the Norwegian snow densities.

Based on the previous work of Sturm et al. (2010), a new model that directly includes the variability in weather and climate is developed. This model is also an extension of the model of Sturm et al. (2010), using other explanatory variables and new model parameters.

Several models consisting of different explanatory variables were tested against Norwegian snow data. Models A–J are the best combinations of variables of a wider selection of combinations where the effect of each variable was explicitly tested. Therefore, the internal correlation between variables (Table 3) was ignored in this study, but for further refinement of the models this could be considered. The final model included the following predictor variables: (1) measured snow depth; (2) elevation (masl) at the location; (3) temperature – accumulated positive degrees; and (4) wind – accumulated wind speeds above 2 m/s when the temperature is below freezing point. These variables need to be known for each location where snow density is to be predicted. The explanatory variables are derived based on weather data from weather stations nearby the location of interest/measurement locations.

Weather SDM includes the seasonal variability by the variables wind and temperature. The elevation of the location did not influence the model much, but was chosen based on the CRPS score. An explanation of the small effect is that the effect of location altitude is already explained in the model through the correction of temperature due to altitude.

Weather SDM gives prediction errors, MAE, that are under 10%, and thus in the same range as what was found in historical records that can be expected when using the relatively sparse information about densities from one snow storage observation campaign. Weather SDM gives about 30% lower errors than Sturm's model, and the CRPS score is improved by about 50%. This illustrates the effect of including year-to-year variability of the weather in the model.

Even if the variability among areas and years is considered through the use of weather variables, it is conceivable that there is some other kind of variation in the properties of the snow densities that does not emerge through the explanatory variables. Thus, a random year-area effect is added in the Weather&ObsDensity SDM. This random effect is estimated by year-specific measurements in two different ways. The best predictions came from a model with random effects. Weather&ObsDensity SDM used information from year and area-specific measurements already at the start while estimating model parameters from the training data set. The use of year and area-specific measurements assisted in bias removal and improved the model by up to 16%. The improvements are the largest for the five first measurements and stagnate after 10–15 observations. This illustrates the effect of including information about the specific year. Nevertheless, the improvements are low compared to the effort of collecting these samples. Thus, it can be questioned how many samples should be recommended, at least after the five first samples have provided information about year-specific conditions beyond what the weather data can tell. Using the same approach with Sturm's model directly would also make this model able to catch interannual variations in density due to different weather conditions. This approach could be used in areas with sparse meteorological information.

The snow density prediction model, Weather SDM, developed in this study, can be useful for estimating the snow density instead of manual measurements. As it takes wind speeds, temperature, elevation and snow depth as input variables, it is easy to use and applicable wherever weather stations are available and representative. This model is compared to the original model on which it is based (Sturm et al. 2010) and gives more reliable predictions for Norwegian snowpacks. Analyses show that information gained by collecting five to 10 annual snow densities in each area can improve the predictions.

REFERENCES

REFERENCES
Bruland
O.
Sand
K.
Killingtveit
Å.
2001
Snow distribution at a high Arctic site at Svalbard
.
Nord. Hydrol.
32
,
1
12
.
Brun
E.
Martin
E.
Simon
V.
Gendre
C.
Coleou
C.
1989
An energy and mass model of snow cover sidtable for operational avalanche forecasting
.
J. Glaciol.
35
,
333
342
.
Færevåg
Å.
2013
Predicting Snow Density. Master Thesis, NTNU
,
Trondheim
,
Norway
.
Førland
E. J.
Allerup
P.
Dahlström
B.
Elomaa
E.
Jónsson
T.
Madsen
H.
Vejen
F.
1996
Manual for Operational Correction of Nordic Precipitation Data
.
DNMI Report 24/96, Met.no, Oslo
.
Gamerman
D.
Lopes
H. F.
2006
Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference
.
CRC Press
,
Boca Raton, FL
.
Gilks
W. R.
2005
Markov Chain Monte Carlo
. In:
Encyclopedia of Biostatistics
(P. Armitage & T. Colton, eds)
,
4th edn. John
Wiley & Sons, Chichester
.
Jonas
T.
Marty
C.
Magnusson
J.
2009
Estimating the snow
.
J. Hydrol.
378
,
161
167
.
Lehning
M.
Völksch
I.
Gustafsson
D.
Nguyen
T. A.
Stähli
M.
Zappa
M.
2006
ALPINE3D: A detailed model of mountain surface processes and its application to snow hydrology
.
Hydrol. Process.
20
(
10
),
2111
2128
.
Liston
G.
Sturm
M.
1998
A snow-transport model for complex terrain
.
J. Glaciol.
44
,
498
516
.
Liston
G. E.
Elder
K.
2006
A distributed snow-evolution modeling system (SnowModel)
.
J. Hydrometeorol.
7
,
1259
1276
.
Lundberg
A.
Richardson-Näslund
C.
Andersson
C.
2006
Snow density variations: consequences for ground-penetrating radar
.
Hydrol. Process.
20
(
7
),
1483
1495
.
Lunn
D. J.
Thomas
A.
Best
N.
Spiegelhalter
D.
2000
WinBUGS – a Bayesian modelling framework: concepts, structure, and extensibility
.
Stat. Comput.
10
(
4
),
325
337
.
Spiegelhalter
D.
Thomas
A.
Best
N.
Lunn
D.
2007
OpenBUGS User Manual
,
version 3.0. 2
.
MRC Biostatistics Unit
,
Cambridge
.
Steppuhn
H.
1976
Areal water equivalents for prairie snowcovers by centralized sampling
. In:
44th Annual Western Snow Conference
,
Calgary, Alberta
, 1976
.
Sturm
M.
Holmgren
J.
Liston
G. E.
1995
A seasonal snow cover classification system for local to global applications
.
J. Climate
8
(
5
),
1261
1283
.
Sturm
M.
Tara
B.
Liston
A. G. E.
2010
Estimating snow water equivalent using snow depth data and climate classes
.
J. Hydrometeorol.
11
,
1380
1394
.