Modelling the spatial distribution of snow water equivalent at the catchment scale taking into account changes in snow covered area

Abstract. A successful modelling of the snow reservoir is necessary for water resources assessments and the mitigation of spring flood hazards. A good estimate of the spatial probability density function (PDF) of snow water equivalent (SWE) is important for obtaining estimates of the snow reservoir, but also for modelling the changes in snow covered area (SCA), which is crucial for the runoff dynamics in spring. In a previous paper the PDF of SWE was modelled as a sum of temporally correlated gamma distributed variables. This methodology was constrained to estimate the PDF of SWE for snow covered areas only. In order to model the PDF of SWE for a catchment, we need to take into account the change in snow coverage and provide the spatial moments of SWE for both snow covered areas and for the catchment as a whole. The spatial PDF of accumulated SWE is, also in this study, modelled as a sum of correlated gamma distributed variables. After accumulation and melting events the changes in the spatial moments are weighted by changes in SCA. The spatial variance of accumulated SWE is, after both accumulation- and melting events, evaluated by use of the covariance matrix. For accumulation events there are only positive elements in the covariance matrix, whereas for melting events, there are both positive and negative elements. The negative elements dictate that the correlation between melt and SWE is negative. The negative contributions become dominant only after some time into the melting season so at the onset of the melting season, the spatial variance thus continues to increase, for later to decrease. This behaviour is consistent with observations and called the "hysteretic" effect by some authors. The parameters for the snow distribution model can be estimated from observed historical precipitation data which reduces by one the number of parameters to be calibrated in a hydrological model. Results from the model are in good agreement with observed spatial moments of SWE and SCA and found to provide better estimates of the spatial variability than the current model for snow distribution used in the HBV model, the hydrological model used for flood forecasting in Norway. When implemented in the HBV model, simulations show that the precision in predicting runoff is maintained although there is one parameter less to calibrate.


Introduction
Snow is an important hydrological parameter in the northern hemisphere and quantifying the snow reservoir is necessary for water resources assessment and for mitigating 5 the potential hazard of the spring flood. In order to successfully simulate the temporal evolution of the snow reservoir, snow melt and the snow covered area (SCA), the spatial probability density function (PDF) of snow water equivalent (SWE) plays a key role (Luce and Tarboton, 2004;Essery and Pomeroy, 2004;Luce et al., 1999;Liston, 1999;Buttle and McDonnel, 1987). Furthermore, the spatial PDF of SWE is known to vary 10 throughout the snow season. This was observed by Pomeroy et al. (2004) during the melt period and through the entire snow season by Alfnes et al. (2004). The algorithms used to describe the spatial PDF of SWE in hydrological models thus have to take this feature into account.
The spatial PDF of SWE often serves as the basis for modelling SCA. The temporal 15 development of SCA is important in hydrology and in land surface schemes in atmospheric models. The dynamics of runoff is affected by changes of the area generating melt water and flux accounting must be carried out separately for snow-free and snowcovered fractions of a grid in a land-surface scheme (Essery and Pomeroy, 2004;Liston, 1999). The snow cover depletion curve (SDC) can be derived from the spatial 20 PDF of SWE and describes the relationship between SCA and spatially averaged SWE (Martinec et al., 1994;Luce et al., 1999;Luce and Tarboton, 2004). Luce et al. (1999) derived the SDC from integrating a generic PDF of SWE which shifts to the left as melting proceeds. Essery and Pomeroy (2004) assumed a log-normal distribution of SWE when they showed how the sign of the correlation between melt and SWE influences Introduction in SDC for single catchments which translates to inter-annual variability in spatial PDF of SWE. A highly relevant parameter that has to be taken into account when modelling snowmelt and the evolution of snow-free areas, is the correlation between SWE and melt. Essery and Pomeroy (2004) show that, given a PDF for SWE, the relation be- 5 tween changes in SCA and mean SWE (the SDC) varies according to the sign and magnitude of the correlation between melt and SWE. There has been some debate in the literature regarding the nature of this correlation and Faria et al. (2000) found that the spatial distribution of daily melt was negatively correlated to the distribution of SWE within a boreal forest stand. Pomeroy et al. (2001), however found no spatial co- 10 variance between melt energy and SWE in dense mature spruce forest, although this does not directly describe the correlation between melt rate and SWE. Furthermore, Pomeroy et al. (2004) found negative correlation at small scales (<100 m), and medium scales (<2000 m) and even positive correlations at the catchment scale (<200 km). In addition to these non-conclusive findings, both Pomeroy et al. (2004) and Skau-15 gen (2007) reported that the relationship between spatial mean and variance of SWE is not monotonous throughout the accumulation and melting season. At the very beginning of the melting season the spatial mean decreased whereas the variance increased slightly then later declined with the mean. This behaviour is also seen in studies in the Swiss Alps, where the spatial mean of SWE plotted against the spatial standard devia-20 tion shows that their relation is not one-to-one (Egli and Jonas, 2010;Egli et al., 2011). In the Swiss studies, this phenomenon is called hysteresis, suggesting that predicting the variance requires the history of the mean and not just the mean In Skaugen (2007) a method for estimating a temporally varying spatial PDF of SWE was introduced. This distribution has the ability to reproduce the observed variability 25 in shape of the PDF caused by accumulation and melting events. The spatial distribution of SWE presented in Skaugen (2007) could, however, be applied only to snow covered areas, and did not take into account the development of snow-free areas in a catchment. In order to take the method a step further and making it suitable for implementation in a hydrological model, changes in SCA is derived from the spatial PDF of SWE and the intensity of the melting event.
Correlation between snowfall events and also between accumulated SWE and melt plays a crucial role in the proposed method for estimating the spatial PDF of SWE. In our study we discuss how the correlation between melt and SWE and the hysteresis 5 effect may be linked.
The hydrological model used operationally in Norway for flood forecasting is the Swedish HBV model (Saelthun, 1996;Bergström, 1992;Killingtveit and Saelthun, 1995). This model has been the dominant operational hydrological model in the Nordic countries for several decades and is fitted with a snow routine in which the spatial 10 distribution of snow is modelled as the sum of uniform-and log-normally distributed snowfall events. In this study we will compare modelled results of spatial moments of SWE (mean and standard deviation) and SCA, modelled with the snow distribution routine of the HBV model and the model developed in this paper against observed values. Finally we will compare runoff predictions and simulations of SWE and SCA by the HBV 15 model using both the current snow distribution routine and the new model. The main objective of this paper is to present a method for estimating the spatial PDF of SWE at the catchment scale through estimating the temporally varying spatial moments of SWE while taking changes in SCA into account.

Unit fields
In Skaugen (2007), the PDF of accumulated SWE was approximated as a correlated sum of gamma distributed unit fields, y(x), where x 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(ν 0 , α 0 ) with PDF: Introduction where α 0 and ν 0 are parameters. The mean of the unit equals E (y) = ν 0 /α 0 and the variance equals Var(y) = ν 0 /α 2 0 . The choice of distribution is motivated partly from studies reporting the gamma distribution as a suitable choice for the spatial distribution of precipitation (Onof et al., 1998;Mackay et al., 2001), SWE (Kutchment and Gelfan, 5 1996;Skaugen, 2007) and snow-depth (Egli et al., 2011) and partly because of the practical mathematical features of the gamma distribution. The method applied in this paper 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. The spatial PDF of SWE is subsequently modelled as a gamma distribution with parameters: 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.
During the snow season, the snowpack may experience series of melting and accu-15 mulation events. It is obvious that 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 20 study, SCA is set equal to 1 (full coverage) for every snowfall event, whereas a melting event implies a reduction in coverage. The chosen approach for assessing the spatial moments of SWE is to represent melting and accumulation events separately in the covariance matrix and let the spatial moments be weighted by changes in SCA. Introduction

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 Skaugen (2007) the spatial mean of a snowfall event that comprises n units, y i , can be written as; and the variance as: 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 . This is a departure from Skaugen (2007) where c is a tuned, nondynamic value and not a function of the number of units n. The variance of z (t) is thus:

15
Since c(n) is the ratio between covariance and variance, it is thus the average correlation for the n(n − 1) pairs of units and equal to: The estimation of the moments for a single event of n units is thus quite straightforward. The values for ν 0 and α 0 can be calibrated or we can estimate the spatial conditional mean and standard deviation from historical precipitation events. From analysis of a 19 year long timeseries of precipitation, the spatial conditional mean (m) and standard deviation (s) of precipitation were found to follow a functional relationship of the type, s = am h (Skaugen and Andersen, 2010). If we choose the unit event to have a mean of m = 0.1 [mm], the corresponding standard deviation is estimated from s = am h and we can determine ν 0 and α 0 by Eq.
(2). The moments of SWE for a single event can be estimated similarly and thus the correlation coefficient in Eq. (6) can be determined. The average correlation c(n) is a declining function of n, which corresponds with the 10 results of Zawadski (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 three following subsections show, however, how sequences of melting and accumulation events, together with changes in SCA, complicate the estimate of 15 the spatial moments of SWE. The chosen way to deal with the complexity is to include estimates of the covariance between the units.

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. If the snow 20 reservoir prior to the snowfall event had a snow coverage equal to SCA t−1 , the posterior SCA 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: HESSD 8,2011 Modelling the spatial distribution of snow water equivalent Hence, the task is to estimate the mean and variance for the two areas

The mean
The mean is simply estimated as the sum of the units times the unit mean. For the two areas, the moments are E (z (t)) SCA t−1 = (n + u)ν 0 /α 0 and E (z (t)) (1−SCA t−1 ) = uν 0 /α 0 , respectively. The mean for the new totally covered area is thus:

The variance
For the newly covered area, 1 − SCA t−1 , the variance is estimated using Eq. (5) as: For the previously covered area, it is useful to consider the covariance matrix. The matrix is at all times symmetric, and we can view the additional snowfall as an extension 15 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×n , 2 Cov n×u , and Var u×u , see Fig. 1) and finally estimate Var(z (t)) SCA t as the sum of these parts. The variance of the previous n events, Var n×n is expressed by the updated parame-20 ters of the gamma distribution at time t−1 (by Eq. 2) and equal to Var n×n = ν α 2 . This variance may differ from that obtained by Eq. (5) of accumulation and melting events which prevents Var n×n from being a straightforward function of n, as is the case in Eq. (5). The sum of the elements in the covariance matrix for Cov n×u (and Cov u×n ) is u·n ν 0 α 2 0 c acc and the sum for Var u×u is u The correlation coefficient, c acc , is estimated as if the total variance Var(z (t)) SCA t and Var n×n were estimated using Eq. (5) with u + n and n elements respectively. This is an 5 approximation since we do not know if Var n×n is estimated by Eq. (5) or may be the result of previous accumulation and melting events as discussed above. The equation Var (n+u)×(n+u) = Var n×n + 2 Cov n×u + Var u×u is then solved for c acc . The correlation coefficient c acc , is thus estimated as: and finally 15 As the four parts of the covariance matrix describing the variance of the previously covered area now are estimated, we can write Var(z (t)) SCA t−1 as and finally, the spatial variance of the total snow covered area, SCA t as:

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−1 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 5 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 (nonconditional moments).

The mean
The non-conditional mean after the melting event, is estimated as: 15 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 20 is seen as an extension of the elements of the covariance matrix similar as for 11495 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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×n , 2 Cov u ×n and Var u ×u , see Fig. 2) and finally estimate Var(z (t)) SCA t as the sum: 5 The only possible negative contribution in Eq. (10) is Cov u ×n , since both Var n×n and Var u ×u 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×n is expressed by the updated parameters of the gamma distribution (Eq. 2) at time t − 1 and equal to Var n×n = ν α 2 . 10 The negative covariance contribution Cov n×u (or Cov u ×n between melt and SWE is estimated as: 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 15 zero, but for the latter case, the total variance becomes zero and we can estimate c mlt from Eq. (10), Var n×n + Var u ×u = 2 Cov u ×n , which gives: and: 8,2011 Modelling the spatial distribution of snow water equivalent For specific values of u (u ≤ n), we approximate c mlt (u ) as a linearly increasing function of u :

HESSD
The variance of the melting events is estimated using Eq. (5) Var u ×u = u ν 0 5 Finally, the total variance of the new conditional distribution after a melting event is computed as:

Estimating changes in snow covered area (SCA)
After a snowfall event, the SCA for a catchment or part of a catchment is set equal 10 to 1 in the HBV model. The same procedure for estimating SCA after an accumulation event is applied here. For a melting event, however, the estimate of changes in SCA is somewhat more elaborate. In Dingman (1994) the energy requirements for transforming a snowpack into meltwater is stated as: Q = Q 1 + Q 2 + Q 3 where the different energy quantities refer to warming the snowpack to a uniform temperature of 0 • C de-15 grees (Q 1 ), producing a certain fraction of meltwater contained in the snowpack (Q 2 ) and transforming the snow into meltwater (Q 3 ). 2011 Modelling the spatial distribution of snow water equivalent where h m is SWE, h wret is free water stored in the snow pack (usually a fixed percentage of SWE in hydrological models), T s and T m are snowpack -and melting point temperatures respectively and the rest of the parameters are constants. All the energy quantities are linear functions of the depth, h of SWE, so an assumption that areas with the smallest SWE are the first to become snow free, due to smaller energy re-5 quirements, 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 α as a model for the PDF of SWE. f a (z ) = 1 Γ(ν) α ν x ν−1 e −αz α, ν, z > 0. 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 (Skaugen,10 2007; Essery and Pomeroy, 2004), but a uniform distribution has also been used (Liston, 1999;Egli 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 Eqs. 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 f a . 15 In accordance with a negative correlation between snow depth and melt, 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 a = X 0 f a dx. The value of X is thus the limit where the frequencies of the melt distribution f s , are higher or equal to the frequencies of the accumulation distribution, f a . Also, since f s is not bounded 20 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 + x, will be left snow-free since p s (X +x) is smaller 25 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 HESSD Introduction Recall that the reduction in SCA red is relative, i.e. it is the reduction from the previous 5 snow-cover, which is also the probability space of both f a and f s , and thus equal to 1.

Results
In this section we first present how the results of the proposed algorithms for estimating the spatial moments of SWE and SCA (hereafter called the G model) compare with observed data. Then we present data that justifies the assumption of gamma distributed snowmelt. Finally, we show results from the implementation of the G model in a hydrological model. For five catchments we compare the model performances of the HBV model using both the G model and classic, lognormal snow distribution (LN model) currently used in the HBV model. The LN model uses a uniform spatial distribution of SWE up to a specified threshold of accumulated SWE. For additional snowfall events, 15 each snowfall event is log-normally distributed through a calibrated coefficient of variation (CV) at a specified set of quantiles, i.e. each additional snowfall event has a spatial PDF of fixed shape (through the calibrated CV) regardless of its intensity. The spatial distribution of melt is uniform and reduction in SCA occurs when the SWE associated with a quantile becomes zero. The reduction in SCA is thus the sum of quantiles with 20 zero SWE. The snow routine of the HBV model does not keep track of the spatial moments of accumulated SWE so it is not straightforward to assess the modelled spatial PDF for this model. In this study the SWE values for the different quantiles are fitted to a lognormal distribution, and the spatial moments are derived from the parameters of the fitted distribution.

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 data sets exist in Norway. The first is the data set from Norefjell in Southern Norway and was previously presented in Skaugen (2007). 5 At Norefjell snow surveys at 1000 m a.s.l. were carried out every second week along a 2 km long snow course. Snow depth was measured every 10 m and density measurements were taken twice for each snow course. 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 data 10 set is from the Norwegian Water Resources and Energy Directorate (NVE) research site for snow at Filefjell, southern Norway (Stranden and Grønsten et al., 2011). The site is situated at 1000 m a.s.l., and has a stable snow cover throughout the months of November-April. The vegetation is grass and willow thicket (<50 cm). About 45 stakes are placed at a flat stretch of 450 m, ensuring that exactly the same point is measured 15 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 four locations. Eight surveys were made at this site. Both these sites represent areas smaller than the typical catchments scale, but measurements of the spatial moments of SWE at the catchment scale is not known in Norway. 20 SCA is estimated for the sites by counting the non-zero fraction of measurements. From the observed conditional spatial mean of SWE and SCA, we can derive the nonconditional values of accumulation and melt (n and u), which are input to the snow distribution models. The output from the models is the conditional-mean and standard deviation and SCA. For calibrated values of α 0 , ν 0 and CV, the agreement between simulated and observed values (Figs. 3 and 4) is quite good for both models. The LN model has a better prediction of the conditional mean for Norefjell whereas the opposite is the case 5 for Filefjell. For the conditional standard deviation the G model has a better fit for both sites, which is also reflected in the estimated CV. SCA is reasonably well predicted for both sites by both models. When SCA reduces, the simulated conditional moments deviate from the observed for both models.
For uncalibrated values of α 0 , ν 0 (CV is still calibrated) the agreement between the results for the G model and the observed is still quite good, and do not deteriorate much from the calibrated results (Figs. 5 and 6). This result implies that the G model can be parameterized using observed (precipitation) data, and reduces hence, without much loss in precision, the number of parameters to be calibrated by one. 15 For the snow seasons 2009-2010 and 2010-2011 an attempt to measure the spatial distribution of snow melt took place at the site at Filefjell. The spatial distribution of snowmelt was estimated from differences in SWE measured at each stake after a melting event. From the two seasons of measurements seven events could be analysed. Figure 7 shows that a gamma distribution is entirely plausible and that more 20 intense melting gives a less skewed spatial distribution. These findings thus support the assumption that the spatial frequency distribution of melt can also be modelled as a sum of correlated gamma distributed variables.

Implementing the G model in the HBV model
The algorithms for estimating the spatial moments of SWE and thus the parame- is temperature and precipitation. The snow-routine is a semi-distributed degree-day model so that the moisture accounting is performed separately for ten elevation zones. The Norwegian version of the HBV uses the lognormal spatial distribution for SWE as described above (the LN model) and in Killingtveit and Saelthun (1995). Runoff dynamics are controlled by two linear reservoirs with three outlets. When testing the new 5 snow distribution model, the G model replaces the original LN model and all other routines are kept as in the original model. The HBV model with its two snow distribution models is calibrated and validated for five Norwegian catchments located at relatively high altitudes where snow is of significant hydrological importance. Two versions of the model were calibrated. The original HBV model (HBV LN) and the HBV model with the new snow distribution routine, where α 0 and ν 0 are estimated from precipitation data (HBV G). The models were automatically calibrated using a Marcov Chain Monte Carlo (MCMC) routine programmed in R (Soetaert and Petzholdt, 2010). Table 1 shows the calibration and validation results for the 5 catchments.

Complex variance-mean relationship and correlation between melt and SWE
In Figs. 3-6 we observe an increase in observed spatial standard deviation at the onset of the melting period. For both Norefjell and Filefjell the maximum standard deviation can be seen at the time of the first reduction in the spatial mean of SWE (at time 22 for 20 Norefjell and time 5 for Filefjell). The measurements at time 8 is the highest for Filefjell, but this estimate of spatial standard deviation is based on only three measurements (SCA is 0.11) and is therefore uncertain to such a degree that it can be disregarded. The standard deviation simulated by the G model captures this phenomenon to some degree in that the highest simulated standard deviation coincides in time with the ob-

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | after this time. If we consider the last two terms in Eq. 10 (Eqs. 11 and 13), we find that these terms regulate whether the changes to the variance prior to the melting event, Var n×n = ν α 2 , are negative or positive. At the start of the melting season, the negative contribution of Eq. (11) is more than compensated by Eq. (13) and the variance is increased or stable. As the melting season proceeds, the negative contribution increases 5 and the total variance decreases. A different (non-linear) function for c mlt could perhaps give a better fit to the observed data, but this was not investigated. It is interesting to note that as a mathematical consequence the sign of the correlation between SWE and melt comes out negative, since otherwise the variance would never decrease. The sign of this correlation has been debated in the literature (Skaugen, 2007;Essery and Pomeroy, 2004;Faria, 2000), but it follows from the derivation of the spatial variance of SWE in this paper, and also from physical reasoning (that less energy is required to melt smaller amounts of SWE) that correlation between melt and SWE is negative. Table 1 shows that, according to the Nash-Sutcliffe R2, both models perform well and 15 quite similarly for the five catchments. The validation results are slightly better with the G model for the catchments Atnasjø and Narsjø, and slightly inferior for the other catchments. Assessments of snow variables that are simulated (SWE and SCA), but not the objective for calibration (runoff) can be made. Figures 8 and 9 show observed and simulated runoff, simulated SWE and SCA for the catchments Narsjø and Junkerdalselv 20 respectively. In the bottom panels of the figures we also find SCA estimates derived from satellite (MODIS) optical images. Runoff and SCA (panels a and c) are plotted for spring 2007 and the snow season of 2008 and SWE (panel b) is plotted for the entire validation period . For the Narsjø catchment (Fig. 8) the G model gave the better fit to the observed, mainly due to more simulated snow. The simulations of SCA, simulations of SWE (panel b) perhaps not for the right reasons. It is clear from Fig. 9 (panel c) that the HBV LN model simulates too much snow which it fails to melt during spring. The result is an accumulation of SWE and perennial snow at the highest altitude levels (this can also be seen from the plot of SCA (panel c). It is notable that such an obvious erroneous simulation of SWE does not prevent a good precision in 5 runoff prediction. The overestimated SWE, which in summer will produce a continuous input of melt water, is compensated for by an error where output of water is reduced. Such a compensating error may, for example, be a smaller reservoir for subsurface water. A possible explanation for why HBV G gives a more realistic simulation of SWE may be that the mutual dependencies in the modelling chain of (i) spatial moments of 10 SWE conditioned on observed precipitation characteristics, (ii) the evolution of snowfree areas and (iii) runoff dynamics in spring, represent a rather strict conditioning that discourages unrealistic values of SWE. Both models underestimate SCA for this catchment, but HBV G to a higher degree than HBV LN. The underestimation of SCA (and runoff) in the 2007 season by the 15 HBV G model for the Junkerdal catchment provides a cue for a brief discussion of the sequel of this study. The current study has been devoted to estimating the spatial PDF of SWE and to developing algorithms that relates the changes in SCA to the spatial PDS of both SWE and snow-melt. In an ongoing study we have used the algorithms developed for relating the PDFs of SWE and melt to changes in SCA to update the 20 snow reservoir, i.e. the PDF of SWE, from satellite derived SCA. By manipulating the spatial PDF of melt, we can adjust the PDF of SWE so that the differences in observed and modelled SCA are matched. This can carried out both for satellite derived SCA higher and less than modelled SCA. Such an updating is thus totally dependent on the knowledge of the spatial PDF of SWE hence a very important feature of the current 25 study has been to provide its estimates.

Reducing the number of calibrated parameters in hydrological models
Ideally we would like the operational hydrological models to provide reliable predictions for more elements in the hydrological cycle than just runoff. When hydrological models are also calibrated against other variables and quite accurate predictions are obtained, the precision in runoff predictions does not necessarily increase. In fact, the opposite is seen when hydrological models are also calibrated against, for example observed SCA (Parajka et al., 2007). This is probably due to sub-optimal model structure and over-parameterisation of hydrological models. This is a challenge for the study of the interactions and interdependencies of different hydrological processes, and hinders the use of hydrological models as a tool for gaining scientific insights. In this study we have 10 developed a model that increases the number of tuning parameters by one compared with the original model if we calibrate the parameters α 0 and ν 0 . However, if α 0 and ν 0 are estimated from observed precipitation data, the number of calibrated parameters is reduced by one and equal or even better results are obtained in predicting runoff and snow variables.

Conclusions
A method for estimating the spatial frequency distribution for SWE is proposed. The distribution is dynamic in that its parameters change according to melt and accumulation events. An algorithm for estimating changes in SCA is also introduced. The model simu- 20 lates spatial standard deviation and change of snow cover well when compared with observed data. Also, when the model is implemented in a hydrological model with parameters estimated from observed precipitation data, the snow reservoir is credibly simulated. Although the number of parameters to be calibrated is reduced by one, there is no loss 25 in precision in predicted runoff. 8,2011 Modelling the spatial distribution of snow water equivalent Through the mathematical-statistical formulation of the model, the correlation between melt and SWE is necessarily negative.

HESSD
In an ongoing study we use the algorithms developed for relating the PDFs of SWE and melt to changes in SCA to update the snow reservoir from satellite derived SCA HESSD 8,2011 Modelling the spatial distribution of snow water equivalent   HESSD 8,2011 Modelling the spatial distribution of snow water equivalent   HESSD 8,2011 Modelling the spatial distribution of snow water equivalent   HESSD 8,2011 Modelling the spatial distribution of snow water equivalent   HESSD 8,2011 Modelling the spatial distribution of snow water equivalent