A model for the spatial distribution of snow water equivalent parameterized from the spatial variability of precipitation

Snow is an important and complicated element in hydrological modelling. The traditional catchment hydrological model with its many free calibration parameters, also in snow sub-models, is not a well-suited tool for predicting conditions for which it has not been calibrated. Such conditions include prediction in ungauged basins and assessing hydrological effects of climate change. In this study, a new model for the spatial distribution of snow water equivalent (SWE), parameterized solely from observed spatial variability of precipitation, is compared with the current snow distribution model used in the operational flood forecasting models in Norway. The former model uses a dynamic gamma distribution and is called Snow Distribution_Gamma, (SD_G), whereas the latter model has a fixed, calibrated coefficient of variation, which parameterizes a log-normal model for snow distribution and is called Snow Distribution_Log-Normal (SD_LN). The two models are implemented in the parameter parsimonious rainfall–runoff model Distance Distribution Dynamics (DDD), and their capability for predicting runoff, SWE and snow-covered area (SCA) is tested and compared for 71 Norwegian catchments. The calibration period is 1985–2000 and validation period is 2000–2014. Results show that SD_G better simulates SCA when compared with MODIS satellite-derived snow cover. In addition, SWE is simulated more realistically in that seasonal snow is melted out and the building up of “snow towers” and giving spurious positive trends in SWE, typical for SD_LN, is prevented. The precision of runoff simulations using SD_G is slightly inferior, with a reduction in Nash–Sutcliffe and Kling–Gupta efficiency criterion of 0.01, but it is shown that the high precision in runoff prediction using SD_LN is accompanied with erroneous simulations of SWE.


Introduction
Snow is an important hydrological parameter in the Northern Hemisphere and in Norway approximately 30 % of the annual precipitation falls as snow.Snow and snow-related hydrology have a significant impact on nature and society in such regions.Seasonal snow ensures variation in outdoor activities and considerable investments in infrastructure for tourism and hydropower are dependent on stable seasonal snow.Apart from snow-related hazards such as spring melt floods and avalanches, snow may negatively affect construction safety and traffic flow at airports, roads and in urban areas.Information about snow conditions at the local, regional and national scale is therefore important for the early warning of hazards, as well as for tourism, hydropower production planning and water resources management.
Operational snow models have evolved differently for hydrology than for meteorology and avalanche warning.Whereas the model development in the latter two scientific disciplines usually include detailed, multi-layered, physically based process representations, snow models in hydrology are typically calibrated empirical relationships between snow variables and the modest model forcing at hand, i.e. snow accumulation and melt vs. precipitation and temperature.One reason for such a discrepancy in modelling approaches is that calibrated hydrological snow models have proved themselves at low temporal resolutions (i.e.24 h res-olution; Anderson, 1976) and for stationary climatic conditions.Another reason is that hydrological snow models are expected to provide simulations at the catchment scale, for which there are usually no estimates of more nonstandard hydrological model forcing such as, for example, wind and radiation.In addition, the governing equations for the physics of hydrology at the small scale have proven difficult to scale up in time and space to be relevant for catchment hydrology (Kirchner, 2006).
For predictions in ungauged basins and in a changed climate, however, calibrated empirical relations in snow models cannot be expected to give reliable and useful results.Skaugen et al. (2015) used the Distance Distribution Dynamics (DDD) model (Skaugen and Onof, 2014) for prediction in ungauged basins with model parameters estimated from catchments characteristics.When analysing the deviations in performance between the calibrated and the regionalized versions of the DDD model, the regionalized degree-day factor for snowmelt and the coefficient of variation (CV) for the spatial probability density function (PDF) of snow water equivalent (SWE) emerged as the parameters most responsible for poor regionalized results for runoff.
A realistically modelled spatial PDF of SWE is important for the temporal evolution of SWE, snowmelt and snowcovered area (SCA) (Buttle and McDonnel, 1987;Liston, 1999;Luce et al., 1999;Essery and Pomeroy, 2004;Luce and Tarboton, 2004).In the literature, many models for the PDF are proposed, especially for the period of time of maximum accumulation, such as the log-normal distribution (Donald et al., 1995, Saelthun, 1996), the gamma distribution (Kutchment and Gelfan, 1996;Skaugen, 2007;Kolberg and Gottschalk, 2010;Skaugen and Randen, 2013) and the normal distribution (Marchand andKillingtveit, 2004, 2005).Helbig et al. (2015) investigated the spatial PDF of snow depth for three large alpine areas and found that the gamma and the normal distributions were better suited than the log-normal distribution.In Alfnes et al. (2004), Skaugen (2007) and Skaugen and Randen (2013), it was demonstrated through the repeated measurements of the same snow course during the accumulation and melting seasons that the spatial PDF of SWE changed its shape continuously during the periods of accumulation and melting.During the accumulation period, the spatial distribution of SWE would become less positively skewed as accumulation progressed and increasingly more positively skewed as melting progressed.Good simulation of the evolution of SCA is especially important since it controls the runoff dynamics of the spring melt flood and is the basis for properly accounting the energy fluxes in land-surface schemes in atmospheric models (Liston, 1999;Essery and Pomeroy, 2004;Helbig et al., 2015).
In this study we will test an alternative method for parameterizing the spatial PDF of SWE.In the alternative method the spatial PDF of SWE is modelled as a dynamic gamma distribution and is hereafter denoted SD_G (Snow Distribu-tion_Gamma).The parameters of SD_G are estimated solely from observed spatial variability of precipitation; i.e. all its parameters are estimated prior to the calibration of the hydrological model against runoff.Information on the spatial variability of precipitation is available at many sites, which makes it possible to use the method for prediction in ungauged basins.Downscaled climate changes projections may also provide such information so that effects of climate change on snow conditions and hydrology may be assessed.In using such a method, the current dependency of calibration in hydrological snow models is reduced.
SD_G is described in Skaugen (2007) and has since been developed in Skaugen and Randen (2013).The method was tested at small test sites and found to model the spatial moments of SWE and SCA well (Skaugen and Randen, 2013) but has, however, not been implemented in a hydrological model and hence not been tested for larger scales and as a tool in operational hydrology.In this study, the SD_G is implemented in the DDD model and its performance is compared with the currently used snow distribution model, the Snow Distribution_Log-Normal (SD_LN) (Killingtveit and Saelthun, 1995;Saelthun, 1996).SD_LN distributes SWE log normally in space with a fixed, calibrated CV.It has been used operationally in Norwegian hydrology for many years, although it has the feature of being a calibrated model and hence not suitable for climate change studies and for predictions in ungauged basins.In addition, a fixed CV, and hence an assumption of perfect spatial correlation, is not supported by observations (Alfnes et al., 2004), and in a recent paper Frey and Holzmann (2015) show that a log-normal spatial distribution of SWE with a fixed CV of introduces so-called "snow towers".For high-elevation areas, and for the highest quantiles of the distribution, snow survived the summer and accumulated to give an overall positive trend in SWE, which was not observed.
The main objective of this paper is to evaluate if SD_G is a suitable alternative for use in rainfall-runoff models.We will compare simulated results of runoff, SWE, SCA and snow cover duration simulated with DDD using the current model, SD_LN, and with the alternative, SD_G, for 71 catchments in Norway.Time series of satellite-derived SCA from MODIS (Moderate Resolution Imaging Spectroradiometer) images are available for the catchments, so simulated runoff and SCA will also be compared against observed values.

Method
The proposed method requires an estimate of the spatial PDF of SWE at all times during the snow season.As in Skaugen (2007) and Skaugen and Randen (2013) we model the spatial PDF of Z (the accumulated positive SWE, not including zeros) as a two-parameter gamma distribution.We hence need the estimates of the mean, E(Z ), and variance, Var(Z ), in order to estimate the shape, ν, and scale, α, parameters of the gamma distribution.The following subsec-tion describes how E(Z ) and Var(Z ) are estimated for accumulation and melting events.Accumulation and melting events may change the spatial extent of SCA, which will require special consideration when estimating the E(Z ) and Var(Z ).In this study SCA is set equal to 1 (full coverage) for every snowfall event, whereas a melting event implies a reduction in coverage.With estimates of E(Z ) and Var(Z ), the parameters of the gamma distributions are calculated as In the first subsection, the model for estimating the statistical moments, E(Z ) and Var(Z ), for the accumulated sum of SWE, Z , is presented.As in Skaugen and Randen (2013), the moments are derived from the sum of correlated gamma distributed unit fields, y(x) [mm], where x represents space.For the remainder of the paper the unit field, y (x), is denoted y.
Section 2.1.1-2briefly address the estimation of E(Z ) and Var(Z ) for accumulation and melting events with a changing SCA.The derivation for accumulation events differs from that presented in Skaugen and Randen (2013) and is presented in detail.For melting events, however, only the resulting equations are presented since the full derivations can be found in Skaugen and Randen (2013).
Section 2.2 describes how change in SCA is estimated after a melting event and Sect.2.3 describes briefly the hydrological model and its current model for the spatial distribution of SWE, SD_LN.
The final subsection, Sect.2.5, describes the procedure for testing and comparing the new model for the spatial distribution of SWE, SD_G against the current, SD_LN.The data used will also be presented here.

Statistical moments of spatial SWE
The PDF of Z does not contain zeros and is hence conditional on snow.For the non-conditional distribution of SWE, which also includes zeros, the variable SWE is denoted Z.The unit fields of snowfall are distributed in space according to a two-parameter gamma distribution, y = G(ν 0 α 0 ), with PDF: where is the gamma function and α 0 and ν 0 are shape and scale parameters respectively.The mean of the unit equals E (y) = ν 0 /α 0 and the variance equals Var (y) = ν 0 /α 2 0 .When estimating the moments for the sum of n units, Z (n) = n i=1 y i , we have to take into account that the unit fields are correlated.This has no bearing on the mean, E(Z ), but affects the variance, Var(Z ), i.e.
where the function c(n) is the average correlation over n units.
From Eq. ( 4) we see that if we have perfect and constant correlation between the y's, c(n) = 1, the variance of Z equals Var(Z ) = n 2 ν 0 α 2 0 .From Eqs. ( 3) and (4)we see that the relationship between the standard deviation, σ Z , and the mean, E(Z ), is a straight line with slope equal to ν −0.5 0 , σ Z = ν −0.5 0 E Z .However, if we have no correlation between the y's, c(n) = 0, the variance equals Var(Z ) = n ν 0 α 2 0 , which gives a relationship between σ Z and E(Z ) as a curved line that departs from that of perfect correlation by n −0.5 , σ Z = (ν 0 n) −0.5 E(Z ).The variance, however, is linearly related to the mean.Correlation between the units, c(n), gives a relationship between the mean and the standard deviation that is something between the two cases described above.A typical analytical approximation to the spatial and temporal correlation function for precipitation is an exponentially decaying function with either time or space as argument.Zawadski (1973Zawadski ( , 1987) ) found exponential decorrelation for rainfall for both time and space.As n (number of summations) may be considered a variable akin to time, c(n) is approximated by an exponential correlation function: where D is the decorrelation range in which the correlation equals 1/e (Zawadski, 1973).The variance of Z can now, with Eqs. ( 4) and ( 5), be expressed as From measured, positive (i.e.not including zeros) precipitation over an area we can observe the relationship between the spatial mean and spatial variance of precipitation.Furthermore, we can estimate the two unknowns, D and α 0 , from such data by nonlinear regression.Figure 1a shows a scatter plot of spatial mean and standard deviation of positive precipitation (from the Norwegian Meteorological Institute) with a fitted function of the type Eq. ( 6).From Fig. 1b, where the spatial mean and standard deviation are plotted in log-log space, we see that the relationship is not that of a power law, as suggested in Skaugen and Randen (2013) and Skaugen and Andersen (2010), since a straight line will not represent the point cloud very well.The parameters a 0 , ν 0 and D are estimated from an analysis of the variability of precipitation as shown in Fig. 1 for the catchment of interest.A mean of the units has been chosen as E (y) = ν 0 α 0 = 0.1 mm, since 0.1 mm is the smallest precipitation value measured by the Norwegian Meteorological Institute.

Statistical moments of spatial SWE after an accumulation event
From a single snowfall event of n units on a snow-free surface, the mean and the variance of the snow reservoir Z are estimated according to Eqs. (3) and (4).If there is an additional snowfall event of u units, the mean of the resulting snow reservoir is and the variance is where ν α 2 is the conditional variance prior to the accumulation event.In order to keep the notation simple we say that n is the number of units at t − 1 and u is the number of units of the event at time t.
Equations ( 7) and ( 8) are valid if SCA = 1 for both events.If SCA prior to the new event was reduced due to melting (SCA t−1 < 1), we have to scale the contributions of n and u according to the change in SCA from SCA t−1 < 1 to SCA t = 1; hence the mean is and the variance is (10)

Statistical moments of spatial SWE after a melting event
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.Reduction in snow coverage needs special attention regarding the conditional (Z ) and the non-conditional (Z) moments since we have to determine the spatial moments for the area of the new coverage, SCA t (not including zeros, i.e. conditional moments) and for the area which includes the previously covered part, SCA t−1 (including zeros, i.e. non-conditional moments).
The non-conditional mean after the melting event is estimated as and the conditional mean is We note that the difference in conditional means before and after the melting event is where u is the conditional number of melted units and describes the difference in units when the (relative) reduction in SCA is taken into account.Skaugen and Randen (2013) give a detailed derivation of the conditional spatial variance of SWE after a melting event.
Here, only the final expression is reported: where ν α 2 is the variance of Z prior to the melting event, and c mlt u is the (negative) correlation between melt and SWE and is estimated as a linearly decreasing function of u and equal to It is clear from Eq. ( 13) that estimation of the change in SCA due to melting is needed in order to assess u and consequently Var Z n−u in Eq. ( 14).The next subsection describes such a procedure.

Estimating changes in snow-covered area
After a snowfall event, the SCA for the area of interest (a catchment or a part of a catchment in the case of elevation bands) is set equal to 1.For a melting event, however, the estimate of changes in SCA is more complex.The previous subsection suggests modelling the accumulated SWE as a gamma distribution, f a , with parameters ν and α derived from the estimated mean and variance of accumulated SWE as described above.In Skaugen and Randen (2013), the spatial frequency of snowmelt, f m , was also modelled as a gamma distribution following the same principles as for accumulation, i.e. that the moments can be estimated using Eqs.( 3) and ( 4) with u replacing n.At all times u ≤ n, which implies that until the final melting event occurs, f m is more skewed to the left than f a .
Figure 2 illustrates how the reduction in SCA due to a melting event is estimated.Since the energy requirements for transforming a snowpack into snowmelt is linearly related to snow depth (Dingman, 2002), it is reasonable to assume that areas with smallest SWE are the first to become snow free.Figure 2a shows the PDFs of melt (f m , red) and accumulation (f a , blue).In and much larger than the area of SWE values less than X (the integral of f a up to X, called a, is X 0 f a = a in Fig. 2a).Consequently, the fractional area of SWE values less than X, a, becomes snow free after the melting event.In addition, there are melt values higher than X that reduce the coverage of corresponding SWE values.The sum of these bars adds up to 1 − m and equals the integral In total, the reduction of SCA after a melting event is and is seen in Fig. 2b as the sum of the cross-hatched bars.
Recall that the reduction in SCA 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 m and equal to 1.The correlation of snowmelt c(u) as a function of intensity (u) (see Eq. 14) has not yet been investigated in detail and is, in this study, modelled as that of accumulation.Skaugen and Randen (2013), however, reported empirical evidence supporting such an assumption.The observed features of f m are found to be similar to those of f a , i.e. that the spatial distribution is generally skewed to the left and becomes less skewed as the intensity of melt increases.These features for f m are confirmed by additional measurements of spatial snowmelt by Weltzien (2015).

The hydrological model
The DDD model (Skaugen and Onof, 2014;Skaugen et al., 2015;Skaugen and Mengistu, 2015) is a rainfallrunoff model written in the programming language R (www.r-project.org)and runs operationally at daily and 3-hourly time steps at the Norwegian flood forecasting service at the Norwegian Water Resources and Energy Directorate (NVE).The DDD model introduces new concepts in its description of the subsurface and of runoff dynamics and is developed with the objective of having as many as possible of its model parameters estimated directly from observed data such as maps and runoff characteristics and not through calibration against runoff.In its current version, the parameters of the modules for subsurface and runoff dynamics are all estimated prior to calibration against runoff.Input to the DDD model is precipitation and temperature.The model is semi-distributed in that the moisture-accounting (rainfall and the accumulating and melting of snow) is performed for 10 elevation bands of equal area.The catchment averages of precipitation and temperature are distributed to the elevation bands using calibrated lapse rates.The catchment averaged precipitation can be corrected by multiplying the amount with a constant in order to get the long-term water balance right.Snowmelt is estimated using a degree-day model (Ohmura, 2001;Hock, 2005), where the melted amount is a linear function of the difference between actual air temperature and a calibrated threshold temperature for melting.In the current routine in DDD for the spatial PDF of SWE (SD_LN), the PDF is mod-elled as the sum of uniform and log-normally distributed snowfall events (Killingtveit and Saelthun, 1995;Saelthun, 1996).The distribution is constant for up to a specified threshold of accumulated SWE (i.e.20 mm).Each additional snowfall event is log-normally distributed through a calibrated coefficient of variation, θ CV , and SWE is estimated for nine quantiles and added to previous quantile values.In this way, each additional snowfall event has a spatial distribution of a fixed shape (through the calibrated θ CV ) regardless of its intensity.Moreover, the method implies perfect spatial correlation in that a new snowfall is distributed such that the quantiles with highest SWE always receives most SWE so that the CV of the sum of snowfall events remains a constant.A simple example demonstrates this: if the accumulation of SWE, Z, is the sum of two snowfall events y, Z = y 1 + y 2 , where y ∼ LN(µ y , σ 2 y ) is log-normally distributed with mean µ y and variance σ 2 y , then the mean of Z is E (Z) = 2µ y and the variance is Var (Z) = σ 2 y +σ 2 y +2COV (y 1 , y 2 ).With perfect correlation the variance equals Var (Z) = σ 2 y + σ 2 y + 2σ 2 y (Haan, 1977, p. 56) and it is easily seen that the CV for Z equals that of y, i.e.
The spatial distribution of melt is constant and reduction in SCA occurs when the SWE associated with a quantile becomes 0. The fraction of snow-free areas is thus the sum of quantiles with zero SWE.The model parameters relevant for snow accumulation and melt which are estimated by calibration against runoff include θ CV , describing the spatial distribution of SWE, θ CX , which is the degree-day factor, and θ Ws , which is the maximum liquid water content in the snowpack (see Table 1

of model parameters).
Further details on the DDD model are found in Skaugen and Onof (2014) and in Skaugen and Mengistu (2015).Model parameters that can be calibrated against runoff are denoted by θ with subscripts (e.g.θ CV ) in order to clearly distinguish between estimated and calibrated parameters.From Table 1 we see that 11 model parameters have the potential to be calibrated.The next subsection shows, however, that the number of calibrated parameters for this study is reduced to five (shown in bold in Table 1).

Test of SD_G in DDD
We will evaluate the performance of SD_G, parameterized from observed spatial variability of precipitation, by implementing it in DDD (DDD_G) and compare performance with DDD_LN, in which SD_LN, with its calibration parameter θ CV , is implemented.The parameters D and α 0 for SD_G are estimated for each catchment by analysing the spatial mean and spatial standard deviation of positive precipitation (excluding zero values).The precipitation data, provided by Table 1.Parameters of the DDD model with description and method of estimation.Some parameters (denoted with a * ) have values obtained through experience in calibrating DDD for gauged catchments in Norway.These values are within the recommended range for the HBV model (Saelthun, 1996).Other parameter values are assigned standard values as suggested in the literature.The GIS analyses are carried out using the national 25 × 25 m digital elevation model (http://www.statkart.no).Parameters in bold have been calibrated in this study, by either dataset V1 or V2.

R
Parameter defining field capacity (Skaugen and Onof, 2014) the Norwegian Meteorological Institute, are daily precipitation values from precipitation gauges (a minimum of two stations) located close to the catchment in question and are from the period 1990-2011.
DDD_G and DDD_LN are run for 71 catchments distributed across Norway (see Fig. 3).The catchments vary in latitude, size, elevation and surface cover (see histograms of selected catchment characteristics (CC) in Fig. 4) and constitute thus a varied, representative sample of Norwegian catchments.The runoff data are provided by NVE and we use the period of 1 September 1985 to 31 August 2000 for calibra-tion of DDD_G and DDD_LN and the period of 1 September 2000 to 31 December 2014 for validation.
The following procedure was conducted: the models were initially calibrated using long time series of precipitation and temperature to simulate runoff using a Monte Carlo Markov chain method (Soetart and Petzhold, 2010)   grid was developed, denoted V2 (Lussana et al. 2014a, b), which reduced much of the positive bias in precipitation characteristic of V1 (see Saloranta, 2012).The new meteorological grid (V2) in DDD gives reasonable simulated values of runoff without the need for a calibrated correction of the amount of precipitation (θ Pc ; see Table 1 for parameters of the DDD model).Areal averages of precipitation and temperature values are extracted for 10 elevation zones, which makes it possible to eliminate calibrated precipitation and temperature gradients (θ Plr and θ Tlr ).Three parameters associated with snow accumulation and melt -the correction factor for solid precipitation (θ Sc = 1.0), the threshold temperature for snowmelt (θ Ts = 0 • C) and the threshold temperature for solid and liquid precipitation (θ TX = 0.5 • C) -were fixed, thereby reducing the number of calibration parameters from 11 to 5. For the remaining five parameters, the calibrated values (from using V1 as input) are retained for three parameters (θ Ws , θ v r and θ cea ), whereas for the DDD_LN model, θ CX and the parameter of interest for this study, θ CV , are recalibrated using V2 as input data.In using such a procedure we assume that the three parameters which are calibrated using the V1 data (and, most likely, not optimal using the V2 data as input) will not favour either of the two compared model structures (DDD_LN and DDD_G).When recalibrating the θ CV with V2 data, we attempt to make it as difficult as possible to accept the new spatial frequency distribution of SWE (SD_G).
If we calibrated all three parameters (θ Ws , θ v r and θ cea ) using V2, we could risk that errors associated with the structures of SD_G and SD_LN were compensated by the other three parameters, such that we could not isolate and evaluate the effect of implementing SD_G.In addition, for the DDD_G model, the degree-day factor θ CX was calibrated since correlation between this parameter and θ CV was revealed.It would hence be probable that a θ CX optimized using SD_LN with V1 would not be optimal for testing SD_G.
From almost 1500 optical satellite scenes from MODIS during the period 2001-2015, SCA for each elevation band have been estimated for 69 of the 71 catchments (for two of the catchments SCA observations were not retrieved).Many scenes are discarded due to insufficient light caused by the low solar angle during the Norwegian winter, but for each catchment, about 150 estimates of SCA during the 15 years can be used for validation of the snow distribution models' ability to simulate a realistic evolution of snow-free areas during the snowmelt period.For each MODIS satellite scene, each pixel (500 × 500 m) is assigned a SCA value between 0 and 100 % coverage using a method based on the Norwegian linear reflectance to snow cover (NLR) algorithm (Solberg et al., 2006).The input to the NLR algorithm is the normalized difference snow index (NDSI) signal (Salomonsen and Apple, 2004).

Results
With the procedures and data described in the previous section, we can compare the performances of the DDD model with calibrated PDF of SWE (DDD_LN) and the DDD model with estimated PDF of SWE (DDD_G) with respect to runoff, SWE, SCA and duration of the snow cover for the validation period (1 September 2000-31 December 2014).In Table 3 we present the significant Spearman correlations (with p value < 0.01) between simulation results for these variables and catchment characteristics such as catchment size, areal percentages of lakes, bogs, bare rock and forest and mean elevation of catchment in order to investigate if the results are stratified with respect to the physiography of the catchments.

Runoff
Figure 5 shows different skill scores obtained for runoff simulations for the 71 catchments with DDD_LN (red crosses) and DDD_G (blue circles).The figure is organized such that the catchments are sorted geographically starting from the south-east (S-E), then follows the south-west (S-W) to mid-Norway (M-N) and finally Northern Norway (N-N).Figure 5a shows the Nash-Sutcliffe efficiency criterion (NSE; Nash and Sutcliffe, 1970) and Fig. 5b the Kling-Gupta efficiency criterion (KGE; Gupta et al., 2009;Kling et al., 2012) and Fig. 5c-e the three components of the KGE, correlation, bias and variability error respectively.The variability error is given by the ratio of the coefficients of variation of simulated and observed runoff as suggested in Kling et al. (2012).The mean values of the skill scores for DDD_LN and DDD_G are shown in Table 2 and as straight lines in the plots.We have also added a moving average of the results for enhanced readability.We see from the Fig. 5 and Table 2 that little precision in predicting runoff is lost when using DDD_G.The mean values for NSE, KGE and the different elements of KGE are practically identical.Differences between runoff simulations for DDD_G and DDD_LN are mostly pronounced in the south-east, where, especially for NSE, DDD_LN appears to be consistently better.
Table 3 shows that significant correlation between NSE and CC was only found for catchment area.Such a correlation was not found for KGE; rather, significant negative correlation were found for both models between KGE and the areal fraction of bare rock.The Cryosphere, 10, 1947Cryosphere, 10, -1963Cryosphere, 10, , 2016 www.the-cryosphere.net/10/1947/2016/deviates increasingly as time proceeds.Figure 7a shows a scatter plot of the mean simulated SWE (averaged over the time series) for the 71 catchments by the two models and it is clearly seen that SWE simulated by DDD_LN is higher than simulated by DDD_G although both the precipitation and temperature input are identical for the two models.From linear regression between SWE, precipitation and temperature with time we can estimate simple annual trends.The regression slopes of SWE for both models were correlated with CC and for DDD_LN no significant correlations were found.Significant correlation was, however, found between the slopes of SWE for DDD_LN and the parameter values of θ CV , r S_SWE , θ CV = 0.45, which in turn is significantly correlated with skill score KGE, r KGE_LN , θ CV = 0.40.For DDD_G significant correlations were found between the slopes and lakes, bare rock, bogs and forest.

Snow-covered area and snow cover duration
Figure 8a shows the root mean square error (RMSE) between observed and simulated catchment values of SCA for 69 catchments.Although the mean RMSE does not differ much between the two models (mean(RMSE) = 0.14 for DDD_G and mean(RMSE) = 0.15 for DDD_LN) we can note that SCA is better estimated using DDD_G for 46 out of 69 catchments (67 %).DDD_LN appears to be better in the southwestern part of Norway whereas DDD_G performs better in the other regions.The mean elevation of catchments was found to be significantly correlated to RMSE for simulated SCA using DDD_LN and nearly significantly correlated using DDD_G.The correlation implies that the errors in estimating SCA are, for both models, reduced as the mean elevation of the catchments increase.Figure 8b shows the mean absolute error (MAE) and we see that DDD_G is the superior method with respect to MAE for all regions except for the south-west.The errors are mostly positive indicating a general overestimation of SCA, although underestimation is also found in south-western Norway.The mean value over all the catchments is mean(MAE) = 0.03 for DDD_G and mean(MAE) = 0.06 for DDD_LN.For both models, MAE was significantly correlated to the areal percentage of lakes and the size of the catchment, but not the mean elevation.
The mean annual snow cover duration was calculated as the mean number of days with snow present in the catchment and is shown in Fig. 9.There is a striking difference in this results between DDD_LN and DDD_G.The mean duration of the snow cover of DDD_LN shows almost no variability, is very high and suggests an almost perennial snow cover.This result is consistent with the positive trends for SWE associated with DDD_LN.From Table 3 we see that the snow cover duration are, for both models, significantly correlated with catchment size, fraction of forest and bare rock and the mean elevation of the catchment.

Discussion
Table 2 and Fig. 5 show that, according to the NSE and KGEs, the models are almost identical with respect to the simulation of runoff.This implies that little performance is lost in simulating runoff by introducing the new procedure for modelling the spatial frequency distribution of SWE although there are one parameter less to calibrate against runoff.A reduction in the number of parameters to calibrate reduces the dimensions of the parameter space and thus the parameter uncertainty.In addition, it makes the model less flexible and more dependent on its structure so that possible structural deficiencies can be more easily identified (Kirchner, 2006).These are very important points when the demands on hydrological models move from just predicting runoff to reliable predictions for more elements in the hydrological cycle such as for example SWE and SCA.In addition,  to properly assess the hydrological effects of climate change and to provide useful predictions for ungauged basins, we have to move towards the use of hydrological models with a minimum of calibration parameters.
The major objective of this study is to investigate whether DDD_G gives a more realistic simulation of snow properties, such as a realistic temporal evolution of SWE and SCA during the snow season.Figures 6 and 7 show that DDD_LN gives a pronounced positive trend for simulated SWE, whereas DDD_G gives a small positive trend in SWE that corresponds roughly to that of precipitation (recall that SWE is the accumulated solid precipitation during a period of time).It is notable that such an obvious erroneous simulation of SWE using SD_LN has so little impact on the precision of runoff predictions.A possible reason is that the surplus of snow, located at the highest elevations and for small areal fractions, will not have temperatures high enough, even during summer, to generate intense snowmelt affecting the precision of runoff simulations.In over-parameterized rainfall-runoff models, the optimal runoff simulation is often a system of compensating errors in states (i.e.soil moisture and SWE) and updating one of the states from observations may, in fact, deteriorate the simulation result because the balance of errors is disturbed (Parajka et al., 2007).It is, however, of concern that the method itself introduces trends that could easily be interpreted as a trend in SWE in a climatic study.This problem of "snow towers" in models using a log-normal distribution for SWE with a fixed, calibrated CV has recently been addressed in the literature (Frey and Holzmann, 2015).In Norway, using such a snow distribution model with the well-known Swedish rainfall-runoff model, HBV (Hydrologiska Byråns Vattenbalansmodell; Bergström, 1992), has led to the operational procedure of deleting the remaining snow reservoir at the end of summer.Such a procedure clearly constitutes an example of a model working well (with respect to runoff) but not for the right reasons.This point is further illustrated when we focus on one of the catchments that gives better NSE values using DDD_LN than DDD_G.
The Masi catchment (5543 km 2 ) is located in Northern Norway and is relatively flat (90 % of its area is located below 600 m a.s.l. and its minimum and maximum elevation is 370 and 1085 m a.s.l.respectively) so that the snow melt season is quite short and intense.Figure 10a shows the simulation of SWE using SD_LN with optimized CV (θ CV = 0.88), which gave a NSE value for runoff of NSE = 0.75, and using SD_G, which gave a NSE value for runoff equal to NSE = 0.72.In Fig. 10b we have adjusted the CV value from θ CV =0.88 to θ CV =0.1, and the simulation of SWE using SD_LN no longer exhibits the very strong positive trend seen in Fig. 10a; it looks more realistic and very similar to that of SD_G.The precision of runoff simulation was, however, affected and the NSE value dropped from NSE = 0.75 to NSE = 0.60.A reasonable conclusion may thus be that the slightly higher values for NSE and KGE using SD_LN is at the expense of unrealistic values of SWE.The correlation analysis supports this conclusion (see Table 3).The increase in SWE with time of DDD_LN is not correlated to any CC but to the parameter values of the method for spatial distribution of SWE, θ CV .The parameter θ CV is found to be significantly correlated to the skill score for predicting runoff, KGE; i.e. high values of θ CV gives high values of KGE.The high skill scores for DDD_LN are clearly not due to a realistic process description of snow but rather to an inadequate model structure that gets it right for the wrong reasons.
Figure 8 shows that, in general, SCA is better simulated using DDD_G than DDD_LN.In the south-western region, however, DDD_LN performs better than DDD_G, which underestimates SCA.This region is characterized by a very rugged topography, which may influence both the estimates of SCA derived from the MODIS satellite and the quality of the meteorological grid, as very few meteorological observations are located at high altitudes (Dyrrdal et al., 2012;Saloranta, 2012).Further investigations relating the errors in estimating SCA to other CCs, such as catchment gradients, may explain the difference in results between DDD_LN and DDD_G.Figure 11 shows a typical example where DDD_G has estimates of SCA close to the observed especially during late spring.Naturally, the problem of "snow towers" of DDD_LN influences its ability to simulate a realistic decrease in SCA since small fractions of the catchments remains snow covered at all times.The heavy tails of the optimized accumulation distribution produced by DDD_LN make a complete melt-out of the snow reservoir very dif- ficult.DDD_G, in contrast, provides an accumulation distribution without the heavy tail, which appears as a better choice with respect to the simulation of both SWE and SCA.The difference between the two methods with respect to the modelling of SCA became very clear when we compared the average annual duration of the snow cover.DDD_ LN, due to the positive trends in SWE, ended up with an almost perennial snow cover for most of the catchments (see Fig. 9), whereas DDD_G showed a variability in snow cover durations that is more consistent with the varying climate in Norway.For both models the correlation analysis between snow cover duration and CC showed that the duration of snow cover was positively correlated to catchment size, mean elevation and areal fraction of bare rock (area above the treeline) and negatively correlated to the areal fraction of forest.Since the areal fraction of forest and bare rock is highly correlated, these are expected relations illustrating that both models have a realistic snow distribution with respect to elevation.
A more realistically simulated SCA is important for many applications.Updating of snow and hydrological models using observed SCA is dependent on realistic simulations of SCA.A realistic simulation of SCA is also necessary for the properly accounting of energy fluxes over an area partly covered by snow (Liston, 1999;Essery and Pomeroy, 2004) and is hence important for the assessment of hydrological impacts of climate change.Without realistically simulated SCA, we cannot expect credible simulations for climate projections for neither runoff dynamics nor energy fluxes.
SWE is represented here as the sum of correlated (in time) spatial variables (solid precipitation).Precipitation events such as snow are assumed to be gamma distributed in space with parameters varying with intensity.The parameter scale, α 0 , and decorrelation length, D, are estimated from observed spatial moments of precipitation.Recall that the shape parameter ν 0 , is just set as one-tenth of α 0 through the relation E (y) = ν 0 α 0 = 0.1 mm.From Fig. 1 we see that the variance levels off, and even decreases, for increased spatial mean intensity.The presented model captures this observed feature since the variance will cease to increase as the correlation decreases with intensity (the number of summations).As correlation approaches 0, we will have a sum of independent events.According to the central limit theorem, such a sum will have a normal distribution.The shape parameter of y, ν 0 and the correlation determines the rate of the convergence to a normal distribution.For example, if the decorrelation range is long, then more summations are needed for the spatial frequency distribution of SWE to approach a normal distribution.The literature shows that empirical spatial distribution of SWE has a tendency to be positively skewed.This is especially the case for many observations of SWE in Norway in high alpine areas (Alfnes et al., 2004;Marchand and Killingtveit, 2004;Marchand and Killingtveit, 2005).For more lowland and forested areas, the distribution tends to be more normal (Alfnes et al, 2004;Marchand and Killingtveit, 2004;Marchand and Killingtveit, 2005).In our modelling framework, this would imply that we would expect small shape parameters and long decorrelation lengths for mountain areas and larger shape parameters together with short decorrelation lengths for lower-lying forested areas.Table 4 show correlations and their significance (p values) between the parameters α 0 and D and the CC fraction of bare rock, fraction of forest, mean elevation and catchment area.We see that α 0 is significantly correlated to the mountain/forest and highland/lowland indices as expected.The decorrelation length D is weakly correlated to the mean elevation in a way implying shorter correlation lengths at high altitudes, i.e. contrary to what is expected from reported shapes of the PDF of SWE, and uncorrelated to the other indices.It is promising, and somewhat unexpected, that correlation between α 0 (ν 0 ) and catchment characteristics supports our theory so clearly since the location of Norwegian precipitation gauges, which is has a very poor representation at high elevations (Dyrrdal et al., 2012;Saloranta, 2012), was not expected to discriminate this behaviour very well.The somewhat confusing results of the decorrelation length suggest a dedicated study using a more dense network of precipitation gauges.
As mentioned in the introduction, many models for the spatial PDF of SWE have been proposed in the literature (i.e.normal, gamma, log normal, mixed log normal).The diversity in distributions is often addressed to the physical processes responsible for the shape of the spatial distribution of SWE, which include wind, during and after the snowfall, spatial variability of precipitation and topographic features.This topic is continuously debated in the literature (Liston, 2004;Skaugen, 2007;Lehning et al., 2008;Clark et al., (Blöschl, 1999;Alfnes et al. 2004;Liston, 2004;Marchand and Killingtveit, 2004;Marchand and Killingtveit, 2005).A major problem is that the spatial distribution of snow and SWE is very hard to measure at the appropriate scale, i.e. the catchment scale, which often covers different elevations and both forested and open (alpine) areas.Various airborne observation techniques such as laser scan (Melvold and Skaugen, 2013) and passive microwave (Vuyovich, 2014) are promising but restricted by landscape features such as vegetation and topography and the state of the snow (wet/dry).Consequently, investigations on the spatial distribution of SWE has to rely on in situ measurements, which seldom covers entire catchments.In this study, in situ information (the spatial variability of solid and liquid precipitation), is obtained from the station network of precipitation gauges of the Norwegian Meteorological Institute, which measures precipitation at 2 m above ground.It is highly probable that the observed spatial variability, measured near to the surface, captures information of the influence of the wind on precipitation in general and on snowfall in particular.This assumption is justified by the significant and relatively high correlations seen in Table 4 between the scale parameter, α 0 (and hence, in our case, the shape parameter, ν 0 ), and landscape features such as elevation and vegetation and suggests a sensitivity to the exposure of wind.Johansson and Chen (2003) demonstrate the influence of wind speed on the spatial distribution of precipitation and Mott et al. ( 2011) and Lehning et al. (2008) show that near-surface wind fields highly influence snow distribution patterns through preferential deposition.
The method presented in this study does not include redistribution of SWE due to wind as a driving force for shaping the spatial frequency distribution of SWE at the catchment scale.Some authors suggest that this process occurs on a spatial scale much smaller than the catchment scale (Liston, 2004;Melvold and Skaugen, 2013).In Figure 11 we see that DDD_LN shows a better simulation of SCA for the start of the melting period than DDD_G for at least two of the years (2011 and 2014).The reason to why DDD_LN simulates the initial development of snow-free areas better than DDD_G is probably that SD_LN produces a generally more positively skewed distribution of SWE than SD_G, and has, hence, a higher frequency of small values of SWE that melts quickly.Whereas the distribution of SD_G, which in general seems to be more appropriate, should perhaps have a fraction of the catchment populated with small values of SWE in order to simulate this observed initial development of snowfree areas.By including redistribution due to wind, we might produce areas of shallow SWE, such as over wind-exposed ridges which are known to become free of snow rather early in spring.
Finally, it is important to keep in mind that this study aims at determining the spatial frequency distribution of SWE for elevation bands for a catchment.These areas may comprise several square kilometres.The spatial distribution of SWE for distributed hydrological modelling, i.e. simulating the amount of SWE at specific locations, is another, much more challenging, task which involves taking into account very small-scale (< 25 m according to Lehning et al., 2008) landscape features and their complex relation to accumulation, melting and redistribution of SWE.

Conclusions
In this paper a method for estimating the spatial frequency distribution of SWE is implemented in the parameter parsimonious rainfall-runoff model DDD.The new method, first described by Skaugen (2007) and further developed by Skaugen and Randen (2013) and here, has its parameters estimated from observed spatial variability of precipitation measured from precipitation gauges.The new method (SD_G) has hence no parameters to be optimized from calibration against runoff unlike the current operational snow distribution routine (SD_LN), which has one calibration parameter.The new method gives a dynamic presentation of the distribution of SWE, which, at the start of the accumulation season, may be positively skewed but converges to a more symmetrical distribution as the accumulation season progresses.The parameters of the method show significant correlations with catchment characteristics discriminating between sheltered and wind-exposed areas.
DDD_G is tested for 71 catchments in Norway and little loss in precision of predicted runoff is seen when skill is measured with the NSE and KGE criteria.SWE is simulated more realistically in that the seasonal snow is melted out every year and no trend in SWE is observed, which is consistent with the absence of trends in precipitation and temperature.The current operational routine for snow distribution (SD_LN), however, displays a tendency to produce ever increasing "snow towers" (Frey and Holzmann, 2015), which in turn gives the erroneous impression of an increasing trend in SWE and unrealistic annual durations of snow cover, which for most catchments approach a full year.Such behaviour can be remedied by adjusting the optimized parameter value for the spatial snow distribution, θ CV , but at the expense of the precision of simulated runoff.The simulated SCA for both SD_G and SD_LN is compared to MODIS-

Figure 1 .
Figure 1.Scatter plot of the spatial mean and spatial standard deviation of observed precipitation over a catchment.Equation (6) is fitted to the data by nonlinear regression (red line).Bottom panel shows the scatter plot in log-log.

Figure 2 .
Figure 2. Schematic of how changes in SCA are estimated.(a) f m and f a are the spatial frequency distributions (PDF) of snowmelt and accumulation respectively.m, 1 − m, a and 1 − a are partially integrated values of the PDFs.(b) The integral of the PDFs for successive intervals of SWE and melt and their spatial coverage.The cross-hatched bars constitute the reduction in SCA.
Fig. 2b we have plotted the integral of the PDFs for successive intervals of SWE, so each horizontal bar represents a fractional area (see the x axis) of SWE or melt values.The horizontal bars for each integrated PDF sum up to unity, i.e. the entire area covered by snow.The figure illustrates that melt values less than X cover a large area (the www.the-cryosphere.net/10/1947/2016/The Cryosphere, 10, 1947-1963, 2016 T. Skaugen and I. H. Weltzien: Spatial distribution of SWE integral of f m up to X, called m, X 0 f m = m in the Fig. 2a

Figure 3 .
Figure 3. Location of the 71 catchments used to evaluate the new subsurface routine.The results presented in the next section are organized with respect to the regions south-east (S-E), south-west (S-W), mid-Norway (M-N) and Northern Norway (N-N).

Figure 4 .
Figure 4. Histograms of catchment characteristics for the 71 catchments: (a) mean of the hillslope distance distribution, d; (b) areal percentage of lakes; (c) areal percentage of bogs; (d) catchment area; (e) mean elevation; (f) areal percentage of glaciers; (g) areal percentage of forests; and (h) areal percentage of bare rock.

Figure 6 .
Figure 6.Time series of simulated SWE using DDD_G (blue line) and DDD_LN (red line) for catchment Tansvatn in Southern Norway.
c, d  show plots of the slopes of the regression lines.Whereas both precipitation and temperature show very modest annual rates of change, both models simulate increasing SWE with time, but DDD_LN, on average, 5 times as much as DDD_G.If 100 days a year may serve as an estimate of days with solid precipitation, the increase in SWE due to the positive trend in precipitation comes very close to the trend in SWE found for DDD_G.Positive trends of SWE greater than 5 mm year −1 was found for 26 out of 71 (37 %) catchments for DDD_LN model and 7 out of 71 catchments (10 %) for the DDD_G model.

Figure 7 .
Figure 7. Scatter plots of mean SWE simulated with DDD_G and DDD_LN for 71 catchments (a), annual slope of SWE (b), annual slope of precipitation (c) and temperature (d).

Figure 8 .
Figure 8.(a) Root mean square error (RMSE) for simulated SCA for DDD_G (blue) and DDD_LN (red).(b) Mean absolute error (MAE) for simulated SCA for DDD_G (blue) and DDD_LN (red).Moving averages and mean values of RMSE and MAE are shown with the same colour code.

Figure 9 .
Figure 9. Mean annual days of snowcover in the catchments for DDD_LN (red) and DDD_G (blue).Moving averages are shown using the same colour code.

Figure 10 .
Figure 10.Time series of simulated SWE for the Masi catchment in Northern Norway with DDD_G (blue) and DDD_LN (red).In (a) SWE is simulated with optimized CV = 0.77, which gives a NSE = 0.75.In (b) SWE is simulated with adjusted CV = 0.1, which gives a NSE = 0.60.Using DDD_G gives a NSE = 0.72.

Figure 11 .
Figure 11.Time series of simulated SCA with DDD_G (blue) and DDD_LN (red) together with MODIS derived SCA (green circles) for catchment Tansvatn in Southern Norway.
100.Derived from GIS. θ Ws [%] Maximum liquid water content in snow.Calibrated (V1).Hfelt Mean elevation of catchment.Derived from GIS. θ Tlr [ • C/100 m] Temperature lapse rate (pr 100 m).Not used in this study.θ Plr [mm/100 m] Precipitation lapse rate (pr 100 m).Not used in this study.θ Pc Correction factor for precipitation.Fixed at value 1.0 (see text).θ Sc Correction factor for precipitation as snow.Fixed at value 1.0 (see text).θ TX [ • C] Threshold temperature rain/snow.Fixed at value 0.5 (see text).θ TS [ • C] Threshold temperature melting/freezing.Fixed at value 0.0 (see text).θ CX Mean of distance distribution for bogs.Derived from GIS. Bogfrac Fraction of bogs in catchment.Derived from GIS. Zsoil Areal fraction of zero distance to the river network for soils.Derived from GIS. Zbog Areal fraction of zero distance to the river network for bogs.Derived from GIS.
NOL Number of storage levels.Fixed at value 5 written in R. The time series for precipitation and temperature are mean areal catchment values extracted from the current, operational meteorological grid (1 × 1 km 2 ) which provides daily values of precipitation and temperature for Norway from 1957 to the present day (see www.senorge.no).This meteorological grid is denoted V1.Recently, a new, improved meteorological www.the-cryosphere.net/10/1947/2016/TheCryosphere,10, 1947Cryosphere, 10,  -1963Cryosphere, 10,  , 2016

Table 2 .
Mean values of skill scores for the validation period 2000-2014 simulated with DDD_G and DDD_LN for 71 catchments.KGE_r measures correlation, KGE_b, the bias error and KGE_g the variability error.All skill scores have an ideal value of 1.
using DDD_G (blue) andDDD_LN (red).This example illustrates that SWE simulated with DDD_LN tends to survive the summers at the highest elevations, which results in a positive trend for SWE.Seasonal SWE simulated by DDD_G and DDD_ LN is similar at the start of the time series but www.the-cryosphere.net/10/1947/2016/TheCryosphere,10, 1947Cryosphere, 10,  -1963Cryosphere, 10,  , 2016

Table 3 .
Spearman correlations between simulated model results and catchment characteristics for the validation period 2000-2014.Only significant correlations are shown (p value < 0.01) except for the correlation marked * , which has a p value slightly larger than 0.01 (p value = 0.013).

Table 4 .
Spearman correlations between model parameters and catchment characteristics indicating alpine and lowland areas where the spatial distribution of SWE is expected to vary.The bracketed numbers indicate significance (p value).
www.the-cryosphere.net/10/1947/2016/The Cryosphere, 10, 1947-1963, 2016 derived SCA and SD_G has the lower RMSE.The difference in simulated SCA between the two models is especially seen for median to low values of SCA.SD_LN can be seen to simulate better SCA at the beginning of the melt season, suggesting that SD_G has a too-low frequency of low SWE values.NVE supports an open data policy; real-time and near-realtime data are available at http://www.nve.no/en/15Water/Data-databaser/Real-time-hydrological-data/ and historical data are freely available at request to hydrology@nve.no.