A model for the snow water equivalent derived from stratigraphy observations in northern Sweden

The seasonal snow cover is usually described by its depth and density as well as spatial and temporal extent. The snow cover extent affects the albedo of the ground surface, the snow depth and density governs the insulation capacity, and the snow density determines the permeability of the snowpack (Gray & Male 1981). By combining snow depth and snow bulk density, snow water equivalent (SWE) can be estimated which gives information about the amount of water stored in the snowpack (Armstrong & Brun 2008). These three parameters (snow depth, snow density and SWE) show spatial variations depending on the synoptic scale weather patterns (air temperature gradients, wind and precipitation patterns), topography (Mellor 1965; Jaedicke et al. 2000; Geddes et al. 2005; Skaugen et al. 2012) and vegetation (Geddes et al. 2005). A clear example of the weather influence of snow depth and its temporal and spatial variations over 50 years is given in Dyrrdal et al. (2013). In addition, the snowpack is subjected to change due to metamorphic and sublimation processes (Armstrong & Brun 2008). The snow depth is altered by precipitation, compaction and redistribution of snow due to winds. The snow density of the layers within the snowpack changes over the season, mainly due to changes in temperature and moisture in the snowpack, which alter the size, shape and bonding of the snow particles and thereby the porosity of the snowpack. Snow density can also be altered by wind compaction and overburden pressure (Anderson & Benson 1963).

Being able to estimate the SWE is essential for the risk assessment of avalanche hazards and spring melt floods, for runoff monitoring and hydropower production planning. Observation of solid precipitation is known to be difficult (Wolff et al. 2013). Observations using common rain gauges (e.g. Geonor precipitation gauge) will miss parts of the precipitation, especially during cold and high wind speed conditions. Here, SWE can be seen as a cumulative value of precipitation and the influence from weather (e.g. redistribution due to wind). Furthermore, the structure of the snowpack (its hardness, depth, etc.) is important for grazing ability of animals. Many animals rely on scent to localize food under the snowpack. The transport of scents through the snowpack depends on the depth of the snowpack and its structure. The snowpack creates a sheltered environment from the harsh winter conditions above for plants and animals (Kausrud et al. 2008). As reported in Johansson et al. (2011), the structure of the snowpack has changed in northern Sweden since the mid-1990s, for example with increasing number of very hard snow layers. The snowpack with its depth and density insulates the ground and influences the permafrost melt rate (Zhang 2005) and phenology. Thus, a change in the snowpack will influence the permafrost beneath.

SWE gives an integrated value of the water content of the snowpack. Measurements of SWE in the field can be cumbersome and require great effort (Lundberg et al. 2009). We therefore investigate the possibility of using alternative measurements to capture variations of snow density within the snowpack. By parallel measurements of multiple parameters (snow particle size, snow layer hardness, snow density and snow depth), we can investigate their relationship by linear regression analysis.

Snow density measurements are commonly performed by weighing a known volume of snow. Detailed measurements can be obtained by using a snow cutter with a small volume (DeWalle & Rango 2008). By using this method, the stratigraphy of the snowpack can be quantified (assuming that the individual snow layers in the snowpack are thick enough for the snow cutter). Recently, a more sophisticated method of measuring snow stratigraphy has been developed that utilizes the co-variation of physical snow parameters. For example, a snow fork (LK1001, TOIKKA AS, Espoo, Finland) measures permeability which can be converted to snow density. For a more general observation, a snow corer can be used in order to determine the bulk snow density of the entire snowpack and SWE (DeWalle & Rango 2008). As a non-intrusive method, a ground penetrating radar system can be used to study snow stratigraphy (Ingvander & Brown 2013). The result can then be used for building empirical models of SWE (Sundström et al. 2012).

Geldsetzer & Jamieson (2000) investigated the possibility of measuring alternative parameters of the snowpack in order to determine the density of snow. They identified a significant correlation between snow layer hardness subdivided by snow particle size and snow density. In this study, we have further developed their work by using both snow hardness and snow particle size and their relation to snow density on a reference dataset sampled in Järämä, northern Sweden in 2009–2012 (Ingvander et al. 2012). The determined relations have then been applied to the bi-weekly snow stratigraphy observations made at Abisko Scientific Research Station (Abisko SRS) since 1961 (Johansson et al. 2011).

Jonas et al. (2009) derived relationships for snow bulk density and SWE (described below) based on snow depth by exploring an extensive dataset from the Alps. Here, we have applied the method from Jonas et al. (2009) on our snow depth measurements from Abisko SRS and compared the resulting snow layer density and SWE values. The results from Abisko SRS were then used to adjust the Jonas et al. (2009) equation for snow density and SWE in order to apply the method to a different geographical region. Abisko is considered a high latitude, low elevation and relatively dry site. The Alp site in Jonas et al. (2009) is located at lower latitude, has larger elevation and larger amount of precipitation compared with Abisko. The snow depth is up to three times as large as in Abisko.

We hypothesize that the snow density is mainly governed by the porosity of the snow (the ratio between snow particles size and cavities of air in the snowpack) and the snow layer hardness (the bonding strength between snow particles) (Bader et al. 1939; Pielmeier & Schneebeli 2003). By combining information on these parameters, it should be possible to estimate the snow density. Hence, snow density is determined by correlating information on snow particle size, which is affected by the level of metamorphism in the snowpack, and snow layer hardness, which mirrors the bonding strength of the snowpack.

Implementation of the method can be used where there is a lack of either observed snow density or SWE values. Our approach is based on the study by Geldsetzer & Jamieson (2000). In their study, they use a dataset of dry snow observations from two sites, Columbia and Canadian Rockies between 1993 and 2000. We will here test our hypothesis on field observations from Järämä (northern Sweden) and apply the model on the historical dataset from Abisko SRS in order to estimate snow layer density and SWE (which are not measured at Abisko SRS).

Finally, we will compare the estimated SWE values with a SWE estimation based solely on snow depth from Jonas et al. (2009).

In this paper, we use three sets of field data: a reference dataset from three field campaigns in Järämä, northern Sweden, a historical dataset from Abisko SRS and one dataset from Abisko SRS including snow layer density observations. Abisko SRS is located in the Swedish mountain range (68 °21′ N, 18 °49′ E) and runs an extensive environmental monitoring programme (Jonasson et al. 2012). The snow stratigraphy dataset from Abisko SRS is described further in Johansson et al. (2011). The data used here consist of bi-weekly data of snow depth and snowpack properties and cover the winters 1961/62–2013/14. In total, the Abisko SRS dataset consists of 524 snow stratigraphy profiles and observations from 2,552 individual snow layers. As a complement to the regular observations, snow layer density observations were made at Abisko SRS during spring 2009 (observations from 14 snow layers). The reference dataset was sampled at Järämä (67 °53′ N, 21 °57′ E) in 2009, 2011 and 2012 and consists of 20 stratigraphy profiles with a corresponding 101 individual snow layers. The dataset and site are described in Ingvander et al. (2012) where the methodology of snow particle size analysis from Abisko SRS was compared with other existing methods (Fierz 2009; Ingvander et al. 2013).

In the Abisko SRS dataset, the snow hardness and snow particle size are registered. In Järämä, these two parameters are registered together with snow layer density. This allows us to study the relationship between the morphology parameters and the snow density. The snow hardness was sampled manually by pressing the wall of the snow pit using a pen (hard), finger, fingers and fist (soft) (Fierz 2009). In the Abisko SRS sampling protocol, there are six classes of snow layer hardness (large cavities (R0), very low (R1), low (R2), medium (R3), high (R4), very high + ice hard layer (R5)). A snow layer can also be classified as an ice layer. But, analyses have shown that what one observer classified as an ice layer another observer classified as a layer with very high snow layer hardness. For this reason, these two classes are grouped together (R5 thus consists of both very high and ice layer observations).

The snow particle size was determined using visual interpretation at Abisko SRS and in Järämä, by using digital image interpretation. The two methods are compared in Ingvander et al. (2012). Snow particle size is, in the Abisko SRS protocol, sampled in six classes (nuts (E_N), peas (E_p), rice (E_R), semolina (E_S), flour (E_F) and flakes (E_Flake)). Intermediate classes sometimes occur in the Abisko SRS dataset: for example, E_FS size between flour and semolina and E_SR size between semolina and rice. Snow density was measured using a Toikka AS snow fork. An average value from three measurements in each layer has been used. The snow layer density observations at Abisko SRS during spring 2009 were obtained by weighing a fixed volume of snow. For both the Abisko SRS and Järämä datasets, the snow has been analysed by stratigraphic analysis where each layer in the snow is documented separately. The stratified analysis gives us detailed information on the snowpack and the different parameters measured.

The Abisko SRS snow stratigraphy dataset has the advantage of being made at one site and covers a time span of more than 50 years, but snow information on density is missing (except during spring 2009). The reference dataset is smaller, but contains measured snow layer density in addition to the stratigraphy observations. Modification of the method in Geldsetzer & Jamieson (2000), based on data from the reference dataset, will enable snow layer density estimations from the Abisko SRS data. From snow layer density, SWE values can be calculated given estimates of thickness of all individual layers (ΔHS) in the snowpack.

Here, i is a number from 1 to the maximum number of individual layers (n) in the snowpack. Equation (3b) specifically takes the internal variations of the snow density within the snowpack into account when estimating SWE.

The significance values (p-values) of the determined regression coefficients are listed in Table 1. For Equation (1), the regression coefficient k_E for classes of larger snow particle sizes (E_N – E_SR) are not as well determined (p < 50%) as for the classes including smaller snow particle sizes (E_S – E_F) (p ≥ 95%). This result is most likely linked to the uneven distribution of the number of observations within each snow particle size class. Within the snow particle size class E_Flake, no regression line could be determined since observations were only available for one snow layer hardness class (R1). Instead, an average value m_E was determined. For Equation (2), the significance for the regression coefficient k_R are all above 50% (R1, R2 and R4 p ≥ 99%) except snow layer hardness class R0, i.e. the change in sign, slope of the regression lines from softer to harder snow layers are well determined. The slightly lower significance level for k_R for snow layer hardness class R3 is not unexpected since the regression line changes sign from negative slope in R0–R3 and positive slopes for R4 and R5. However, the p-values for m_R’ (determined at the centre of each distribution) are high p ≥ 99% (especially important for R3 since the slope of the regression line is weak) for all snow layer hardness classes. The low p-value for snow layer hardness class R0 is most likely due to the low number of observations.

Another difference between the two equations, which can be hard to see from Figure 3, is the number of missing data points (i.e. conditions when snow density data cannot be calculated). Both equations need simultaneously observations of snow particle size and snow layer hardness. If either is missing, this automatically leads to missing data. Nevertheless, applying Equation (1) leads to missing data in 7% of the observations, whereas Equation (2) gives missing data in less than 1% of the observations. This can be attributed to the limited ranges of snow particle sizes in the reference dataset, which has larger impact on Equation (1). In the reference dataset there are no observations of the snow particle size nuts (E_N); hence, no regression estimation for this size class was possible. When using Equation (2), the regression lines have been extrapolated and the equation has been used for E_N for all classes except snow layer hardness R4 (Table 1).

A comparison of observed snow layer density made at Abisko SRS (spring 2009) with modelled snow density values with Equations (1) and (2) show that Equation (1) tends to give higher snow density values than Equation (2). The root mean square deviation from observed snow layer density using Equation (1) is 80.8 kg/m³ and with Equation (2) 59.4 kg/m³. The mean values are significantly different from each other on a 95% confidence level.

Our conclusion is that Equation (2) gives a better estimation of the snow density both concerning values and reliability. In the following calculations only results using Equation (2) will be shown.

As a consequence of the lower snow bulk density values obtained by Equation (4a), the SWE values from Equation (4b) are also lower. The root mean square error for Equation (4b) compared with SWE _Equation2 is 31.7 mm and by using Equation (6b) reduces to 21.0 mm.

By modifying and combining the methods developed by Geldsetzer & Jamieson (2000) and Jonas et al. (2009), an expression for SWE at Abisko SRS has been derived (Equation (6)), which in turn is based on Equation (2) for modelling snow layer density.

From a statistical point of view, Equation (2) and its regression coefficients are determined with a higher accuracy compared with Equation (1). Starting with the reference dataset (used to calibrate the regression models), the results on snow layer density are very similar comparing Equations (1) and (2) (see Figure 2). Both equations result in fewer numbers of observations than observed in densities below 200 kg/m³. Hence from this point it is difficult to say which equation gives the best results. A formal validation of the model is currently only possible for a limited dataset. The root mean square deviation is reduced by 20 kg/m³ by using Equation (2) compared with using Equation (1). Even though the dataset is small it gives an indication that Equation (2) gives a better result.

A comparison of results from Equations (1) and (2) applied to the large Abisko SRS dataset shows a large difference in snow density distribution (see Figure 3). Equation (1) results in a large number of values in the range 300–400 kg/m³ (Figure 3(a)) giving a distribution with a distinct peak, whereas Equation (2) gives a distribution (Figure 3(b)) that is normally distributed up to densities 600 kg/m³, comparable to the shape of snow density distribution presented in Jonas et al. (2009).

With Equation (1), the resulting regression model (see lines in Figure 1(a)) is similar to Geldsetzer & Jamieson (2000), but our snow density values are higher (as also seen compared with Jonas et al. 2009). Whether this is due to the location of the observations (Geldzetser and Jamieson's observations are made in Purcell, Selkirk, Monashee, and Rocky Mountain Ranges in western Canada) resulting in different snow density values compared with those at Abisko SRS, or the fact that they use average values in their model, is difficult to say.

For Equation (1) our reference dataset has a limitation. When using Equation (1), where snow density is correlated to snow hardness, subdivided by snow particle size, there is no separation in the size class E_Flake. All observations have the same hardness, very low (R1). The result of this is more missing data from Equation (1) compared with Equation (2).

In general, Equation (1) gives increasing snow layer density with increasing snow layer hardness index, with some slight modifications depending on snow particle size (Figure 1(a)). Equation (2) on the other hand shows that snow density of a soft snow layer is decreasing with increasing particle size index (corresponding to decreasing snow particle size up to E_F) and that snow density of a hard snow is increasing with increasing particle size index (see Figure 1(b)).

Equation (2) agrees well with the system established by Bader et al. (1939), linking 12 types of snow with different metamorphic and transformation states, to particle size and snow layer hardness. According to Colbeck (1982), the numbers of bonds between snow particles influence the hardness of the snow layer. The snow particle size on the other hand, for a given hardness, influences the pore space in the snow layer, and thus the snow density. That is, we argue that Equation (2) takes the physical processes within the snowpack into account in a better way than Equation (1). A general weakness in regression models can be that temporal development is not taken into account. By developing a model based on six different classes, temporal changes in the snowpack will be captured. For example, if rain falls on snow and refreezing occurs, hard layers may form inside the snowpack. As long as the change in snow layer property is observed (for example that a layer has become harder) modelled snow layer density will also change. Since the models are based on the behaviour of a given class of either snow particle size (Equation (1)) or snow layer hardness (Equation (2)), temporal changes within the snowpack properties and hence the change in snow layer density will be captured by the models.

Our best estimate of SWE from Abisko SRS is given by first modelling the individual snow layer density by Equation (2), and then calculating SWE by Equation (3b). The shape of the distribution of SWE_Equation2 looks similar to the SWE distributions in Jonas et al. (2009) (their Figure 2(a)).

A less detailed way to estimate SWE is to use a bulk snow density value. When applying Equation (2) to the Abisko SRS dataset we obtain consequently larger bulk densities than when using Equation (4a) from Jonas et al. (2009) (see Figures 5 and 6(a)). A comparison of average snow layer values from Equations (1) and (2) does not indicate that Equation (2) on average gives larger values than Equation (1) (seen from both the reference and Abisko SRS dataset). That means that the higher values are not likely to be caused by the choice of regression model. Rather these systematically larger values are probably due to the different geographic settings between Abisko and the Alps. Abisko is located at low elevation (388 m above sea level (a.s.l.)) and high latitude (67 degrees north). The snow depth is generally lower in Abisko because it is situated in a rain-shadowed region behind the Scandinavian mountain range (Kohler et al. 2006). Our best estimation of SWE is based on individual layers, and both the reference and historical Abisko SRS datasets are from the same region (100 km between the field sites) and we thereby have a solid dataset for locally adjusting the Jonas et al. (2009) equation. The Jonas et al. (2009) equation, Equation (4b), has thereby been adjusted into an Abisko-corrected SWE equation (SWE_{Abisko SRS}), Equation (6b).

A regression model for snow layer density based on snow layer hardness and snow particle size has been developed. The snow layer density values have been used to develop a SWE expression from snow depth. Our method shows the ability to use alternative parameters for estimating snow density and SWE. This information can be used to homogenize datasets or analyse datasets that lack information on either snow density or SWE; for example, in estimating SWE from operational snow depth measurement or field campaigns. The results can also be compared with existing SWE models such as the data provided by the Globe Snow project or other products that model SWE on a spatial or temporal scale. A recommended improvement of the SWE_{Abisko SRS} model is, in addition to collecting more data for the snow density regression model, to apply a temporal (seasonal) adjustment to it such as shown by McCreight & Small (2014), in order to capture seasonal changes of the snowpack caused by metamorphism. Furthermore, Equation (6b) is not yet applicable to really deep snow depths such as glacier accumulation areas or similar. This is because the equation is near to linear.

We have found that snow density in northern Sweden is higher than previous results from the Alps (Jonas et al. 2009) and western Canada (Geldsetzer & Jamieson 2000), and the models from these areas are thus not directly applicable to northern Sweden. We therefore recommend using Equation (6b) instead in which the coefficients have been scaled to northern Scandinavian conditions.

To conclude, by applying Equation (6b) to daily snow depth observations, daily values of SWE can be obtained for northern Sweden.

The authors would like to thank the crew, especially Jan Åge Riseth, Hans Tømmervik, Eirik Malnes, Rickard Pettersson and Lars-Thomas Labba, participating in the field campaigns, the Swedish Polar Research Secretariat, and the staff at Abisko SRS. The project was funded by Swedish National Space Board (D70/07:1-2, D116/09:1-2, 139/10:1) and the Research Council of Norway (216434; headed by Jarle W. Bjerke). We would also like to thank Hans Bergström for valuable discussions.