Impact of snow distribution modelling for runoff predictions

Snow is a key component of the water cycle in cold regions. Snow, accumulated in the mountains, stores large amounts of winter precipitation that is reallocated in time to spring and summer runoff, influencing streamflow variability in the downstream regions. The knowledge of the amount of winter snow, snowmelt timing and runoff is essential for the sustainable use of this resource for human consumption, agriculture and hydropower (Viviroli et al. 2011; Sturm et al. 2017), ecological services (Wipf & Rixen 2010; Boelman et al. 2019) and flood control.

Accurate predictions of snow in this resource are challenging because processes at different scales distribute this resource unevenly (Clark et al. 2011). Orographic effects on precipitation, wind and vegetation induce snow accumulation in topographic depressions and erosion in wind-swept areas (Liston & Sturm 1998; Winstral et al. 2002). Spatially variable air temperature and exposure to solar radiation control the snow energy budget and influence snowmelt patterns and runoff dynamics (DeBeer & Pomeroy 2010; Grünewald et al. 2010). Depending on the modelling purpose and data availability, snow spatial distribution in mountain catchments has been represented in different ways. Snowpack energy balance models complemented with snow transport algorithms were used to explicitly simulate all physical processes, including the effect of wind redistribution on snow depth patterns (Lehning et al. 2006; Liston & Elder 2006). These models can reproduce snow distribution at relatively very fine (<100 m) spatial scales (Lehning et al. 2008; Mott et al. 2010; Marsh et al. 2020). However, they require spatially detailed meteorological (e.g., air temperature, wind and solar radiation) and snow observations, which are difficult to obtain and scarce in the mountains. A way to model the snow spatial distribution using only precipitation and air temperature consisted of distributing snowfall or rainfall and melt rates as a function of elevation, i.e., with the temperature-index model. Given its simplicity, this approach was found particularly suitable to be integrated into models for runoff predictions with snow processes represented at the scale of landscape units by lumping similar land uses, soils and elevations (Lindström et al. 1997; Viviroli et al. 2009; Seibert & Vis 2012). Using a simple model structure, Girons Lopez et al. (2020) evaluated the effect of modifying different temperature-index models, discovering that few modifications were most valuable.

Based on snow depth or snow water equivalent (SWE) observations, different studies have explored the link and the predictive power of topographic characteristics (e.g., elevation, slope, aspect and curvature) alone and with wind data to snow spatial distribution (Winstral et al. 2002, 2013; Grünewald et al. 2013; Helbig et al. 2015; Skaugen & Melvold 2019). Topographic characteristics are, in fact, readily extractable from digital elevation models. The combination of different topography-related factors was meant to account for the integrated effect of snow accumulation and melt processes on snow spatial distribution.

The predictive power of the topography-related factors was found to increase with the spatial aggregation of the snow observations (Grünewald et al. 2013; Helbig et al. 2015). When applied at different alpine sites, they were found to be site-specific and dependent on the spatial resolution of snow data. Studies on the relation between snow amount and topography-related factors mainly provided information on snow spatial variability at the peak of the accumulation but did not evaluate the results for runoff predictions.

Next to the spatial distribution of snow amount, snow cover variability was also explored and used to model the effect of uneven snow spatial distribution on melt dynamics (Luce et al. 1998; Luce & Tarboton 2004; Pomeroy et al. 2004). Snow cover depletion curves were found suitable to implicitly model snow processes by linking the temporal evolution of SWE to the snow-covered area on larger model units, e.g., sub-kilometre, than fine grid (Skaugen & Randen 2013; Frey & Holzmann 2015). The shapes of SWE distributions for such relations were assumed a priori and consisted, for example, of two-parameter log-normal or gamma distributions (Skaugen & Randen 2013; Frey & Holzmann 2015). In such approaches, the SWE dynamic was linked to precipitation forcing. The distribution parameters were estimated through model calibration and validation on runoff and fractional snow cover data, respectively.

In addition to the different model designs, the key point is, in fact, the availability of data to estimate or validate snow model parameters. Most studies used runoff and satellite data to estimate these parameters through calibration and to validate the evolution in time of the snow-covered area (Freudiger et al. 2017). Although spatially distributed, satellite snow cover data do not provide quantitative information about snow amount, while runoff provides integrated information on all catchment hydrological processes. Only a few studies used SWE in addition to runoff observations for parameter estimations (Whitaker et al. 2003; Girons Lopez et al. 2020; Nemri & Kinnard 2020). However, these studies mainly used SWE to calibrate one snow modelling approach, i.e., the snow routines based on temperature-index models.

Attempts to model snow redistribution processes for a more realistic representation of snow in mountainous terrain for runoff predictions are still needed, as noted by Freudiger et al. (2017) and Frey & Holzmann (2015).

In line with the need to improve the representation of the snow distribution in mountain catchments for runoff predictions, the aim of this paper is two-fold. On one hand, we investigate the impact of different snow model structures on SWE and runoff predictions. As a novel component, a snowfall distribution (SF) function based on wind shelter factors and wind direction data was implemented in the semi-distributed hydrological model HYPE (hydrological predictions for the environment, Lindström et al. 2010). On the other hand, assess the value of using distributed snow observations in addition to runoff and non-distributed snow observations in hydrological model calibration using snow survey observations and runoff observations from the Lake Överuman Catchment, Northern Sweden.

Meteorological forcing and runoff data covered the periods 2016–2020. The meteorological forcing data consist of precipitation, air temperature and fields of wind direction. Precipitation data were provided by a gridded product based on gauge data, Precipitation Temperature Hydrologiska Byråns Vattenmodell (PTHBV) (Johansson & Chen 2003). The PTHBV product is available at a spatial resolution of 4 × 4 km² and a daily time resolution. Hourly air temperature and wind direction data were provided by the MESAN meteorological analysis of the AROME meteorological model and observations (Häggmark et al. 2000). The spatial resolution of air temperature and wind direction data is 2.5 × 2.5 km². They were aggregated to a daily time scale, the same temporal resolution of precipitation data. Precipitation data were resampled to the same air temperature and wind direction data spatial resolution. Both products are currently used by the Swedish Meteorological and Hydrological Institute (SMHI 2022).

The daily average runoff at the catchment outlet is measured as the residual between reservoir storage change and reservoir outflow. These data were provided by the water regulation company, Vattenregleringsföretagen AB.

Snow surveys were performed in the Överuman Catchment close to the snow accumulation peak in the years 2017–2020 along survey lines distributed in the catchment area (Figure 1). Snow depth measurements were acquired along the survey lines using a ground penetrating radar (GPR) system (Sensor & Software, Noggin Plus with a 500 MHz antenna) connected to a global positioning system (Trimble) pulled by a snowmobile. SWE was obtained by combining the two-way travel time (TWT) of the radar signal, between the snow surface and the ground, with point density measurements. The latter were taken approximately every kilometre along the survey lines (Marchand 2003). The final spatial resolution between SWE measurements was about 1 m distance. SWE estimates were aggregated to the spatial resolution of 25 m distance and used for model calibration and validation.

A digital elevation model (DEM) with a resolution of 25 × 25 m was obtained by bilinear resampling of a 32 × 32 m DEM provided by ArcticDEM (Porter et al. 2018). The digital elevation model was used to extract topographic characteristics of the catchment (i.e., elevation, aspect) and to calculate the wind shelter factor (Winstral et al. 2002). Elevation and aspect were used to discretize the catchment landscape in classes of the model domain. The wind shelter factor was used to model the snow spatial variability induced by wind transport processes in the catchment (Section 4.1). Land-use characteristics were derived from the European Space Agency (ESA) land cover climate change initiative data aggregated into surface water, forest, non-forest, bare soil and glacier classes.

This study used the semi-distributed hydrological model HYPE (Lindström et al. 2010) to simulate SWE and runoff in four hydrological years (2016–2020). The model consists of several subroutines where the different hydrological processes are described. The landscape can be divided into different units. Each is further divided into classes representing different soil and land use properties. Model parameters link landscape characteristics to the physical processes through the model parameters. Two different model spatial configurations were considered, lumped and gridded (Figure 1). In the lumped configuration, the Överuman Catchment consisted of one model unit and meteorological data were aggregated at the catchment scale. In the gridded configuration, the catchment was divided into 2.5 × 2.5 km² units preserving the spatial variability of meteorological forcing at that spatial scale. Each unit was further divided into three elevation zones, four aspect zones (north, east, south and west facing slopes), five land use classes (water, forest, non-forest, bare soil and glacier) and one soil class. In general, in the HYPE model, hydrological processes in the classes are simulated independently within each model unit, except for the SF model introduced in this study, which will be described in more detail in the next section. Runoff from the land classes is aggregated into water classes representing river and lake networks in each unit. It is further routed into the river classes in the next downstream model unit.

In this model version, snow is assumed to cover the entire land class, if present.

The model, including a snow DC to calculate the snow-covered area separately for each land class, is referred to as the DC model (Table 1). The snow DC is used to scale the estimated snow melt and evaporation.

The effective melt is then calculated as in the base model.

In this study, an explicit SF function was introduced to represent the influence of local topography on the spatial variability of snow.

The SF function based on the wind shelter factor (Winstral et al. 2002) was introduced. The wind shelter factor was calculated on a 25 × 25 m² resolution DEM for eight wind directions and five searching distances (from 150 up to 2,500 m). A preliminary investigation of the predictive power of the wind shelter factor and the GPR-derived SWE data (Section 3.2) at 25 m spatial resolution indicated a correlation between the two variables (Figure S1). To further investigate the type of relationship between SWE data and wind shelter factor that could have been used to model the snow distribution, we derived the ratio between the SWE in each GPR point and the mean SWE along the survey line. A log-linear function best represented the relation between the snow spatial distribution and the wind shelter factor (Figure S2). This optimal searching distance of 300 m was chosen as a trade-off among the different analysed searching distances based on the correlation between GPR-derived SWE and wind shelter factor for the eight wind directions (Figure S1). The wind shelter factors were then aggregated on the land classes of the model.

The model calibration was based on the generalized likelihood uncertainty estimation (GLUE) technique (Beven & Binley 1992) and a split-sample approach (Klemes 1986). The GLUE method was implemented as follows: a Monte Carlo sampling from a parameter space was used to generate a prior distribution of model realizations (Beven & Freer 2001). Pseudo-likelihood transformations of performance criteria (goodness-of-fit between observations and predictions) were introduced and used to obtain the posterior distributions, rescaling the prior distribution with the pseudo-likelihood values of each model realization as weight. This rescaling was performed by resampling the prior distribution (with replacement) with the estimated pseudo-likelihood values as a probability. The pseudo-likelihood transformations consisted of scaling performance criteria linearly between zero and one. For this purpose, a minimum and a maximum threshold were defined for each performance criterion. In this way, a strict subjective decision of a behavioural cut-off value that excludes the model realizations below that threshold, as often done in the implementations of the GLUE methodology, is avoided. In addition, the pseudo-likelihood transformations help to combine different performance criteria to rescale the prior distribution into the posterior distribution of model realizations. Due to the limited snow data available, model calibration and validation can be sensitive to a specific year. This limitation is overcome with the split-sample approach that divides the study period into an equal number of years both for calibration and validation (cross-calibration/validation).

Around 10,000 simulations were performed with a Monte Carlo sampling from a feasible parameter space with uniform distribution for both the lumped and gridded HYPE models. Model parameters and their ranges were defined based on previous HYPE model applications in cold environments (Strömqvist et al. 2012; Gelfan et al. 2017), as reported in Table 2.

The minimum and maximum thresholds were 0.5–1 for q₁, ± 20% for q₂ and s₂ and above 0.5 for s₂ performance criteria.

Posterior distributions were then generated by resampling model realizations (with replacement) with a probability given by the likelihood estimated as described above.

The 4-year study period was divided into two equally long 2-year periods. Information from the two different periods was combined to generate two separate final posterior distributions (one posterior distribution for calibration and one for validation).

The correlation of the SWE distributions, s₂ and the runoff Nash–Sutcliffe efficiency, q₁, were then used to evaluate the different models about catchment runoff and snow distribution for the different performance criteria.

The posterior distributions of the model's performances in snow distribution (s₂) and catchment runoff (q₁) were tested for significant differences. For each of the four considered snow models (base, DC, SF and SF + DC, Section 4.2), the non-parametric Kruskal–Wallis test was used to detect whether significant differences occurred among all the posterior distributions. The post-hoc Dunn's test (Dunn 1964) was then used to identify the distribution pairs that were significantly different at the 95% significance level.

This statistical analysis was performed for both lumped and gridded HYPE model configurations in cross-calibration and cross-validation periods.

The results with the gridded model configuration were similar to those with the lumped model configuration (Figure 2(b)). However, the performance of the SF and the SF + DC models showed a higher spread than in the lumped configuration (Figures 2(a) and (b)). The best model performances were found when models were calibrated based on the snow distribution (s₂). The combination with the runoff criteria (q), in particular with Nash–Sutcliffe efficiency (q₁), did not decrease the model's performance (Figure 2(b)).

However, this difference is less marked than the model performances in snow distribution (Figure 2). In contrast with the results of snow distribution, similar model performances were found with both the lumped (Figure 5(a)) and gridded (Figure 5(b)) configurations for a given snow model.

Model performances also differed with calibration criteria, although these differences are less marked than for snow distribution (Figures 2 and 5). The performance of the base model was high when the model was calibrated with the runoff (q) and snow volume (s₁) criteria, not with the snow distribution criteria (s₂). The results of the DC model were similar to the base model. The performance of the SF model varied less with the calibration criteria compared to the performances of the base and the DC models. In particular, the model performance did not decrease with snow distribution observations as for the other two models. The best performances of the SF were found when calibrated with snow and runoff (q + s). Similar results were found for the SF + DC model.

The performance of the base and DC models with the lumped configuration slightly decreased in cross-validation compared to cross-calibration (Figure 7(a)). This decrease in model performance occurred with calibration criteria based on snow distribution observations (s₂). Low performances were already found with the same criteria in cross-calibration (Figure 5(a)).

The use of distributed meteorological forcing (gridded configuration) did not increase the performances of models based only on precipitation partition with elevation, e.g., base and DC models (median among criteria equal to 0.60 in Figure 8(a)). The performances of models with the SF function increased (median among criteria equal to 0.84 in Figure 8(a)). At the same time, the larger spread of the model performances found with the gridded (Figure 2(b)) compared to the lumped configuration (Figure 2(a)) indicated higher model uncertainty. Since the model calibration with lumped and gridded configurations was performed with the same approach, this uncertainty might result from the uncertainty in the gridded meteorological data. This uncertainty is mainly related to the precipitation data that originate from a topography-driven spatial interpolation of point gauge data relatively sparse in the model area (Johansson & Chen 2003). The distributed meteorological data provide some additional information. However, the uncertainty in model forcing for single years or grid cells can lead to increased model uncertainty, i.e., a large spread of the model performances. On the other hand, the assessment of model performance change from cross-calibration to cross-validation indicated a more stable model performance with the gridded models (Figure 4(b)) than in the lumped models (Figure 4(a)) – almost irrespective of snow model types or calibration strategies. The models without SF function, i.e., base and DC, are the most unstable from calibration to validation among the different criteria, especially for the lumped configuration (Figure 4(a)).

Regarding catchment runoff, high model performances (median among the criteria equal to 0.71 in Figure 8(b)) were found for models with the SF function (SF and SF + DC models). The base and DC models showed similar performances (median among the criteria equal to 0.71 in Figure 8(b)), except for the decreased performance for the snow distribution criteria (s₂). The models with the SF function (SF and SF + DC models) better captured both the runoff dynamics and the snow distribution in the catchment compared to the others. In contrast with the model's performances in snow distribution, the gridded configuration did not improve the runoff predictions. However, the model performance is still more stable from calibration to validation than in the lumped (Figure 7).

The analysis of the model performances thus suggests that, beyond spatial model configurations, models with the snowfall function are more suitable to predict both SWE and runoff models than those without this function. This finding is in line with the previous modelling attempts to improve snow volume estimations and thus snowmelt magnitude and timing by including snow redistribution for runoff simulations in mountainous catchments (Frey & Holzmann 2015; Freudiger et al. 2017).

Similar to the previous studies (Girons Lopez et al. 2020; Nemri & Kinnard 2020) that used snow observations in model calibration and showed an overall positive impact, in this study, model performances did not decrease in terms of runoff. In addition, models with SF function proved to be more stable across calibration criteria for snow and runoff predictions and able to include the information from distributed snow observations.

For operational applications where the objective is spring flood forecasting (e.g., for hydropower), an implicit way to represent snow distribution, such as the DC, could also be adopted.

Investigating how to use the information from the distributed snow observations to improve SWE simulations with the DC, understanding whether these results are site-specific, consistent across more snow seasons and how to transfer to other mountainous areas are relevant to the research objectives of future work.

In this study, we used different approaches to represent the spatial distribution of snow in a hydrological model and investigated their impact on SWE and runoff predictions. At the same time, we also investigated the impact of criteria such as snow volume and spatial distribution, in addition to runoff only in model calibration. Models including an SF function based on wind characteristics better predicted SWE distribution at the accumulation peak than models without.

For runoff predictions, models using the SF function and the DC showed good performances. In addition, the gridded model configuration showed similar model performances to the lumped model. To predict SWE and runoff, models with the SF function are more suitable than those without. Despite differences among the calibration criteria for the different snow process representations, the inclusion of snow observations in model calibration provided added value for snow and runoff predictions. As a general conclusion, the use of explicit representation of snow distribution, either by a spatially gridded model configuration or by the explicit SF function, provided stable improvements of the runoff simulation compared to the base model with a lumped configuration.

The authors would like to thank the editor, Bjørn Kløve, Pertti Ala-aho and an anonymous reviewer for comments and suggestions that helped to improve the manuscript significantly. The authors would like also to thank Energimyndigheten, Sweden, for supporting this study (Grant No. 46424-1).

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.