Weather SDM: estimating snow density with high precision using snow depth and local climate

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.

Observed snow densities vary normally between 0.1 g/cm³ for new fallen snow to 0.5 g/cm³ 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.

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/cm³ and the lowest is 0.052 g/cm³. The average is 0.33 g/cm³ with a standard deviation of 0.07 g/cm³. 95% of the observations are between 0.20 and 0.48 g/cm³ (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.

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/cm³ and above 0.5 g/cm³ (Figure 3(c)).

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.

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.

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).

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

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 x₁, x₂, … , x₉. 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.

The snow density increases with age and higher compaction (snow depth). x₁ represents the observed snow depth of the measurement, which is the accumulated sum of new and old snow. x₂ 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. x₃ is the elevation, or the height, in meters above sea level (masl) of the location of the measurement.

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

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

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.

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.

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).

where h is snow depth, z is elevation, T* are accumulated degree days when T_t >0 °C and W* are accumulated wind days when T_t < 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/cm³ 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.

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 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.

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.