Unit fields

In Reference SkaugenSkaugen (2007), the PDF of accumulated SWE was approximated as a correlated sum of gamma-distributed unit fields, y(x), where × represents space. For the remainder of this paper the unit y(x) is denoted y The unit fields of snowfall are distributed in space according to a two-parameter gamma distribution, y = G(V ₀, α₀), with PDF

(1)

where α₀ and ν ₀ are scale and shape parameters respectively. The mean of the unit equals E(y) = ν ₀ /α ₀, and the variance equals Var(y) = ν ₀ /α ₀ ². The choice of distribution is motivated from studies reporting the gamma distribution as a suitable choice for the spatial distribution of precipitation (Reference Onof, Mackay, Oh and WeatherOnof and others, 1998; Reference Mackay, Chandler, Onof and WheaterMackay and others, 2001), SWE (Reference Kuchment and GelfanKutchment and Gelfan, 1996; Reference SkaugenSkaugen, 2007) and snow depth (Reference Egli, Griessinger and JonasEgli and others 2011). Reference SkaugenSkaugen (2007) found that the gamma distribution was in good agreement with the observed spatial distributions of SWE. It can also be noted that the presented method applies to the catchment scale, where the variability of precipitation influences the variability of snow depth (Reference ListonListon, 2004). The values for ν₀ and α₀ can be estimated from the observed spatial mean and standard deviation sampled from precipitation events. From analysis of 19year long time series of precipitation from various areas in Norway, the spatial mean, m, and standard deviation, s, of precipitation were found to follow a functional relationship of the type s = am (Reference Skaugen and AndersenSkaugen and Andersen, 2010), where a and h are determined using nonlinear regression. Once a suitable choice of the mean of the unit is chosen, the corresponding spatial standard deviation, s, can be estimated from s = am^h and ν₀ and α₀ can be determined using Eqn (2) below.

During the snow season, the snowpack may experience a series of melting and accumulation events. These two processes have different spatial frequency distributions and are differently correlated in space. Estimating the spatial variance of SWE after a series of alternating melt and accumulation events is thus a challenge and must include the covariance between the units. Furthermore, SCA varies throughout the season, which necessarily gives a non-homogeneous spatial field of SWE. In this study, SCA is set equal to 1 (full coverage) for every snowfall event, whereas a melting event implies a reduction in coverage.

The procedure for estimating the spatial PDF of SWE is to estimate the spatial conditional mean, E(z’(t)), and variance, Var(z’(t)), of accumulated SWE, z’(t), as functions of the units, y and weighted by changes in SCA. The spatial PDF of SWE is subsequently modelled as a gamma distribution with parameters

(2)

The distribution of z _’ does not contain zeros and is hereafter called conditional (i.e. conditional on snow). For the non-conditional PDF of SWE, which also includes zeros, the variable SWE is denoted z.

Moments and parameters of the gamma distribution of an individual snowfall event

We start the procedure for determining the spatial moments of SWE by investigating a simple case, namely a single snowfall event. According to Reference SkaugenSkaugen (2007), the spatial mean of a snowfall event that comprises n units, , can be written as

(3)

and the variance as

(4)

Note that we have n(n–1) covariance elements since the trace of the covariance matrix consists of the variance for each individual y We estimate the covariance between the units as the average covariance over the n(n–1) pairs of units, and equal to a fraction c(n) of the variance of the individual . This is a departure from Reference SkaugenSkaugen (2007) where c is a tuned, non-dynamic value and not a function of the number of units n. The variance of z’(t) is thus

(5)

Since c (n) is the ratio between covariance and variance, it is the average correlation for the n (n–1) pairs of units and equal to

(6)

The estimation for a single event of n units is carried out using Eqns (3) and (5). With an estimate of the spatial variance from the relationship s = am^h we can also estimate the correlation coefficient in Eqn (6).

The average correlation c(n) is a declining function of n, which corresponds with the results of Reference ZawadzkiZawadski (1973) who found that the temporal correlation of precipitation is a rapidly declining function of time. To estimate the conditional moments from individual snowfall events, we have no apparent use for the estimate of the covariance between the units. The following subsections show, however, how sequences of melting and accumulation events, together with changes in SCA, complicate the estimate of the spatial moments of SWE.

Accumulation events on a previous snow reservoir

Let the snow reservoir, consisting of n units, be increased by a new snowfall of u units. Our task is to estimate the moments of the new spatial PDF of SWE. The snow coverage prior to the snowfall event is denoted SCA _{t _1} whereas the SCA after the event is set equal to full coverage, SCA _t = 1. We determine the mean and the variance for the previously covered part SCA _{t _1} and newly covered part (1 — SCA _{t _1}) separately. The moments for these two areas are assumed independent and the moments for the new totally covered area are estimated as

and

Hence, the task is to estimate the mean and variance for the two areas SCA_t−1 and 1 − SCA_t−1, i.e. and

The mean

The mean is simply estimated as the sum of the units times the unit mean. For the two areas, the mean is and respectively. The mean for the new totally covered area is thus

(7)

The variance

For the newly covered area, 1 –SCA _{t _1}, the variance is estimated using Eqn (5) as

Estimating the variance for the previously covered area, we need to consider the covariance matrix. The matrix is at all times symmetric, and we can view the additional snowfall as an extension of the elements of the matrix, so that after a snowfall event of u units the original n × n matrix becomes a matrix of (n + u) × (n + u) elements. We proceed to estimate the four parts of the matrix separately, Var _{n xn} , 2Cov _{n xu} and Var _{u xu} (Fig. 1), and finally estimate as the sum of these parts.

Fig. 1. Covariance matrix for n + u units. Since the matrix is symmetrical, Cov_nxu and Cov_uxn are equal.

The variance of the previous n events, Var _{n xn} , is expressed by the updated parameters (Eqn (2)) of the gamma distribution at time t—1 and equal to . This variance may differ from that obtained by Eqn (5), using n, because there may be a history of accumulation and melting events which prevents Var _{n xn} from being a straight forward function of n, as is the case in Eqn (5). The sum of the elements in the covariance matrix for Cov _{n xu} (and Cov _{u xn} ) is and the sum for Var _{u xu} is . The correlation coefficient, c _acc, is estimated as if the total variance Var[z’(t)]_SCA and Var _{n xn} were estimated using Eqn (5) with u+ n and n elements respectively. This is an approximation since we do not know if Var _{n xn} is estimated by Eqn (5) or may be the result of previous accumulation and melting events as discussed above. The equation

is then solved for c _acc. The correlation coefficient c _acc is thus estimated as

which gives

and finally

As the four parts of the covariance matrix describing the variance of the previously covered area are now estimated, we can write Var[z’(t)]_SCA as

and, finally, the spatial variance of the total snow-covered area, SCAt, as:

(8)

Melting events

Let the snow reservoir, consisting of n units, be reduced by u units after a melting event. The snow coverage before and after the melting event is SCA _t _₁ and SCA _t respectively, where SCA _t < SCA _{t _1}. We set SCA _{t _1} as 1, so that SCA _t is the relative reduction in snow coverage due to melting, and not the catchment value. The reduction in snow coverage poses a problem in that we have to separate between non-conditional (the area includes a fraction of zero values) and conditional moments. We thus have to determine the spatial moments for the area of the new coverage SCA _t (conditional moments) and for the area which includes the previously covered part, SCA _{t _1} (non-conditional moments).

The mean

The non-conditional mean after the melting event is estimated as and the conditional mean is

(9)

We note that the difference in conditional means before and after the melting event is

where u′ is the conditional number of melted units.

The variance

When assessing the conditional variance after the melting event, the melting event is seen as an extension of the elements of the covariance matrix, similarly to accumulation events. The original n × n matrix becomes, after melting u′ units, a matrix of (n + u’) × (n + u’) elements. We proceed, as in the case of accumulation, to estimate the four parts of the matrix separately (Var _{n xn} , 2Cov _{u′ ×n} and Var _{u ×u′} ; see Fig. 1 and substitute u with u) and finally estimate Var[z’(t)]_SCA as the sum

(10)

The only possible negative contribution in Eqn (10) is Cov _{u xn} , since both Var _{n × n} and Var _u ′_Xu ′ are by definition positive. Negative covariance implies that melting is negatively correlated with SWE, i.e. melting is more intense from areas with less SWE.

The variance of the n events prior to the melting, Var _{n xn} , is expressed by the updated parameters of the gamma distribution (Eqn (2)) at time t — 1 and equal to .

The negative covariance contribution Cov _{n xu} ′ (or Covu’xn) between melt and SWE is estimated as

(11)

The nature of the correlation c _mlt is unknown, but we can estimate the limiting values for no melt (u′ = 0) and complete melt (u′ = n). When (u′ = 0), the covariance is obviously zero, but for the latter case the total variance becomes zero and we can estimate c _mlt from Eqn (10), Var _{n ×n} + Var _{u ×u} ′ = 2Cov _u′_xn, which gives

and

For specific values of u′ (u’ ≤ n), we approximate c _mlt(u’) as a linearly increasing function of u′:

(12)

The variance of the melting events is estimated using Eqn (5):

(13)

Finally, the total variance of the new conditional distribution after a melting event is computed as

(14)

Estimating changes in snow-covered area (SCA)

Recall that after a snowfall event, SCA for the area of interest is set equal to 1 (the same procedure as applied in the HBV model). After a melt event, however, the estimation of changes in SCA is somewhat more elaborate. In Reference DingmanDingman (1994) the energy requirements for transforming a snowpack into meltwater are stated as Q = Q₁ + Q₂ + Q ₃ where the different energy quantities refer to warming the snowpack to a uniform temperature of 0°C (Q1), producing a certain fraction of meltwater contained in the snowpack (Q ₂) and transforming the snow into meltwater (Q3):

where h_m is SWE, h _wret is free water stored in the snowpack (usually a fixed percentage of SWE in hydrological models), c_i is the heat capacity of ice, p_w is the density of water, T _s and T_m are snowpack and melting point temperatures, respectively, and λf is the latent heat of fusion. All the energy quantities are linear functions of the depth, h, of SWE, so an assumption that areas with the least SWE are the first to become snow-free due to smaller energy requirements appears reasonable. This assumption is used to estimate the reduction in SCA after a melting event. Previous sections propose a gamma distribution, f_a, with parameters ν and a as a model for the PDF of SWE,

We also assume that the spatial frequency of melt has a gamma distribution, f_s. Various studies suggest gamma- or log-normal distribution for melting (Reference Essery and PomeroyEssery and Pomeroy 2004; Reference SkaugenSkaugen, 2007), but a uniform distribution has also been used (Reference ListonListon, 1999; Reference Egli and JonasEgli and Jonas, 2009). It is further assumed that the parameters of f_s follow the same principles as for accumulation, i.e. that the moments can be estimated using Eqns (3) and (5) with u′ replacing n. At all times u′ ≤ n, which implies that until the final melting event occurs, f_s is more skewed to the left than fa.

Based on the above, and also on the theoretically established negative correlation between melt and SWE, we state that all points with SWE values less than some value X will be left snow-free. This gives us a reduction in the spatial extent of SCA equal to . We furthermore argue that the value of X is the value where the frequencies of the melt distribution, f_s, are equal to the frequencies of the accumulation distribution, f_a. This is because the areal coverage of the melt values less than is greater than the areal coverage of the accumulated values less than . In addition, since f _s is not bounded to the right, some areas with higher values of SWE than X will be left snow-free after a melting event. For example, if we consider discrete PDFs of the accumulation and melt distribution (p _a and p _s), then a fraction of the total snow-covered area will contain SWE values in the interval defined by, say, SWE=X+x. A smaller fraction of the area with values in the interval SWE = X+× will be left snow-free since p _s(X+x) is smaller than p _a(X+x). When we consider the total snow-covered area, a fraction of all the frequencies of the accumulation distribution f _a for SWE values higher than X will be left snow-free. If these frequencies are summed they will thus represent the area for SWE values higher than X that are left snow-free, The reduction in SCA after a melting event is thus

(15)

Recall that the reduction in SCA_red is relative, i.e. it is the reduction from the previous snow cover, which is also the probability space of both f _a and f _s, and thus equal to 1.

Comparing estimated spatial mean, standard deviation and SCA with observed data

In order to assess the performance of the G_model, we need to compare the results against observed data. Two such datasets exist in Norway. The first is the dataset from Norefjell in southern Norway and was previously presented in Reference SkaugenSkaugen (2007). At Norefjell, snow surveys at 1000 m a.s.l. were carried out every second week along a 2 km long snow course during the winters 2002/03 and 2003/04. Snow depth was measured every 10 m, and density measurements were taken twice for each snow course. Average snow density was measured from two snow pits at locations with average snow depth. This provided a time series of snow course data covering an entire snow season from the start of accumulation to the end of the melting period. Twenty-five surveys were made at this site. The second dataset is from the Norwegian Water Resources and Energy Directorate (NVE) research site for snow at Filefjell, southern Norway (Reference Stranden and GrønstenStranden and Grønsten, 2011). The site is located at 1000ma.s.l., and has a stable snow cover throughout November-April. The vegetation is grass and willow thicket (<50cm). About 45 stakes are placed at a flat stretch of 450 m, ensuring that exactly the same point is measured for each sample. At the stakes snow depths and densities were usually measured once a week throughout the melting season. For each survey of snow depths, snow density was measured at every tenth stake using a snow tube (Reference DingmanDingman, 2002, p. 174). Seven surveys were made in 2011 and eight in 2012. Both sites represent areas smaller than the typical catchment scale, but assessments of SWE at catchment level are usually carried out using snow surveys of a similar spatial scale. SCA is estimated for the sites by counting the zero fraction of measurements. From the observed conditional spatial mean of SWE and SCA, we can derive the non-conditional values of accumulation and melt (n and u), which is input to the snow distribution models. The output from the models is the conditional mean, standard deviation and SCA. Figures 2–4 show observed and estimated (by Eqns (7–9) and (14) for the G_model) conditional mean and standard deviation, SCA and CV for Norefjell and Filefjell. For the G_model the parameters α₀ and ν ₀ are estimated from precipitation data according to the procedure described previously and in Reference Skaugen and AndersenSkaugen and Andersen (2010). The spatial mean and standard deviation of precipitation were sampled for an area of 40 × 40 km² for Norefjell and 30 × 30 km² for Filefjell, and the relationships s = amh were established. The mean of a unit was chosen to be E(y)=m=0.1 (mm).

Fig. 2. Observed and simulated (by G_model and LN_model) conditional mean (a), standard deviation (b), SCA (c) and CV (d) for Norefjell.

Fig. 3. Same as Figure 2 but for Filefjell 2011.

Fig. 4. Same as Figure 2 but for Filefjell 2012.

The parameter CV for the LN_model was calibrated to optimize the estimation of the spatial standard deviation of SWE. The optimal value calibrated for the 2011 season at Filefjell was used to estimate the 2012 season.

Table 1 shows the root-mean-square error (RMSE) for the simulated spatial mean, standard deviation, SCA and CV. The G_model has better RMSE scores for CV than the LN_model for all series. For spatial standard deviation and SCA the G_model gives better RMSE scores for two of the three series. For the spatial mean, the LN_model is better for two of the three series. Note especially that for the Filefjell 2012 series, for which the LN_model is not calibrated, the G_model is only slightly inferior in estimating the spatial mean and considerably better in estimating the other parameters (see also Figs 2–4).

Table 1. RMSE for the G_model and the LN_ model. Bold values indicate the better estimate of G_model and LN_model. The LN_model for 2012 at Filefjell uses a CV calibrated using the data from 2011

Spatial distribution of snowmelt

For the 2009–12 snow seasons, attempts were made to measure the spatial distribution of snowmelt at Filefjell. The spatial distribution of snowmelt was estimated from differences in SWE at the stakes if consecutive measurements showed a decline in SWE. Figure 5 presents three melt events where we have plotted the empirical cumulative distribution functions together with estimated gamma distributions. The moments needed for estimating the parameters are estimated using Eqns (3), (5) and (6), and Eqn (2) is used to estimate the parameters. Note that the input to the estimation of the parameters is just the average melted amount, u. When applying a Kolmogorov–Smirnov test, we find that for 8 of the 12 observed melting events we could not reject the hypothesis that the empirical spatial distribution was gamma-distributed. Choosing the gamma distribution as the melt distribution, f _s, is hence justified.

Fig. 5. Spatial empirical and gamma cumulative distribution function for three measured melting events at Filefjell.