In search of operational snow model structures for the future – comparing four snow models for 17 catchments in Norway

The seasonal snow cover has a significant impact on nature and human society. In Norway, where 30% of the annual precipitation falls as snow (Saloranta 2012), snow and snowfall affect hydropower planning and production, flood forecasting, transport and traffic, biological systems and tourism. The general public, the hydropower industry and public authorities are increasingly in need of more detailed information (in time and space) on snow states and the forecasting of these. Often, the information is needed for sites where measurements do not exist or are sparse. In addition, the need for information of snow conditions in a changing climate may test empirical snow model structures beyond the regime for which they were designed (i.e. calibrated).

Broadly, we can distinguish between two types of operational snow models: (1) empirical and usually calibrated models which are parsimonious both with respect to structural complexity, parameters and forcing (i.e. temperature) and (2) physically based models. A typical example of a calibrated empirical model is the temperature-index (‘degree-day’) model for snowmelt (Omhura 2001), where snowmelt is a linear function of the difference between air temperature and a (often calibrated) temperature threshold under which there is no snowmelt. Applied in hydrological models, the temperature-index factor is often calibrated against runoff, and may hence account for many other processes in addition to that of snowmelt. In general, if parameters in the snow models are not calibrated against snow, the snow simulations may compensate for errors in model structure, meteorological and runoff data. A consequence of using heavily calibrated hydrological models is that the snow in hydrological models and the snow in nature might be two different things. The value of the snow measurements hence decreases since the collected snow data cannot be used for updating and their representativeness is highly uncertain. Calibration is, of course, a challenge when the task at hand is to predict for ungauged basins and for a changed climate since no data are available. However, if the calibrated relationships represent simple, well known physical processes, do not suffer from equifinality issues and sufficient data is available to establish the relations, such procedures may still be useful.

The physically based models typically solve the surface energy balance in order to estimate snowmelt and other snow states. Given sufficient input data such as short- and longwave radiation, wind, humidity, precipitation and surface and air temperature, physically based models (for example Crocus; Brun et al. 1992) are in little need of calibration, but one is on the other hand faced with the difficulty of providing such data for every site of interest.

Information of snow conditions at ungauged sites and for a changing climate hence call for models that: (1) are efficient, i.e. run for areas of catchment/sub-catchment size on an operationally useful temporal scale; (2) are parsimonious with respect to forcing data, i.e. only need precipitation and temperature; and (3) have model parameters that are identifiable, are physically meaningful and constant in time and for different climates.

Through the last decades, many snow model inter-comparison projects have been conducted (Slater et al. 2001; Strasser et al. 2002; Etchevers et al. 2004; Rutter et al. 2009; Essery et al. 2013). These projects compare models for different purposes, such as land surface schemes in numerical weather prediction models, or models for hydrology or avalanche forecasting. The models compared are of different complexity and common conclusions are: (1) that there is no ‘best’ model (Rutter et al. 2009; Essery et al. 2013) and (2) that model complexity does not necessarily translate into better performance (Strasser et al. 2002; Essery et al. 2013). The comparisons have been carried out for one site (Slater et al. 2001; Strasser et al. 2002; Essery et al. 2013), or several sites (Etchevers et al. 2004; Rutter et al. 2009) where snow information is abundant.

In the literature, we find examples of snow models focused exclusively on the dynamics of snow covered area (SCA) (Mir et al. 2015; Pardo-Iguzquiza et al. 2017) or on snow water equivalent (SWE) (Sexstone & Fassnacht 2014; Collados-Lara et al. 2017). The present study is novel in that it does the comparison and validation for a number of catchments for areal values of SCA and for several measurements, located at different elevations of SWE. We believe this to be an important contribution since operational hydrology needs good estimates of snow over areas of some extent spanning different elevations and landscape types (i.e. catchments with both forests and alpine areas). In addition, such a comparison takes into account the effect of the upscaling of model parameters. Upscaling may limit the usefulness of complex, data-demanding models which may prove superior at well instrumented points, but whose performance with regionalized parameters over a catchment is less known (Kirchner 2006).

In this study, we will compare the accumulation and melting of snow from four different snow models, DDD_CX (Skaugen & Weltzien 2016), DDD_EB (Skaugen & Saloranta 2015), SeNorge (Saloranta 2012, 2014), and Crocus (Brun et al. 1992; Vionnet et al. 2011). SeNorge and Crocus are dedicated snow models and DDD_CX and DDD_EB represent snow modules implemented in a hydrological model (the Distance Distribution Dynamics (DDD) model) (Skaugen & Onof 2014; Skaugen & Mengistu 2016). These models hence represent different operational needs, but they share a common goal, namely to simulate snow and snow conditions as realistically as possible. In the suite of models we find both gridded- and catchment-based models. The models further include melt algorithms which are empirical and calibrated, and more physically based. The models hence address to various degrees the prerequisites of the future snow model listed above. Since no models are calibrated either to the SCA or SWE observations, the model comparison will inform on which model structures are best suited for operational snow modelling under a changed climate and for predictions in ungauged basins.

The four models in the comparison differ in the algorithms for snowmelt. Common features of the models are, however, that redistribution by wind and interaction with vegetation are not taken into account. These processes are clearly important in a country like Norway, which is characterized by large areas of forests and mountains where wind is important for altering the snow distribution over many scales (Melvold & Skaugen 2013; Skaugen & Weltzien 2016), but are believed to be more dominant on spatial scales smaller than that which is used here (≥1 km²).

When compared with observed data of atmospheric longwave radiation, ⁠, where is the Stefan-Boltzman constant, ⁠, Equation (3) has half the root mean square error (RMSE) to that of the previously used expression of emissivity from Unsworth & Monteith (1975).

Crocus is a detailed snowpack model originally developed for avalanche forecasting. Crocus simulates the energy exchange within the snowpack and at its boundaries, the absorption of solar radiation with snow depth, the phase changes between solid and liquid water, water flow through the snowpack, mass exchanges due to snow melt and metamorphism of snow. Snow albedo depend on the wavelength of incoming radiation and the type size and age of the snow surface (Strasser et al. 2002). Crocus uses up to 50 layers in order to simulate the internal processes and physical properties of snow and is clearly the most complex model in the ensemble. In this study, Crocus is run regionally on a 0.01° (approximately 1 km) grid within the SURFEX (Surface Externaliseé) interface (Masson et al. 2013). Once the snowpack reaches a threshold of 1 kg m⁻² SWE within a gridpoint in Crocus, the fractional snow cover fraction is assumed to be 1. In practice this means that a grid cell is either completely covered with snow or snow free.

All the models are forced with 1 hourly gridded temperature and precipitation data, provided by MET Norway for the period of 1st September 2013–31st July 2016. The dataset is based on observations from the climate database of the Norwegian Meteorological Institute (eKlima 2017). For gridded precipitation and temperature, a modified Optimal Interpolation scheme is applied (Lussana et al. 2018). The number of observations varies in time and there are regions in Norway, typically at high elevations where the density of meteorological stations is low, which leads to an underestimation for these areas. Lussana et al. (2018) provide uncertainty assessments for the 24 h product to be ±20–30%, for precipitation, depending on intensity, and +0.8–2.4 °C for temperature. In Norway the number of precipitation stations measuring at a 1 hourly interval is much lower than for 24 hours. Hence, the interpolated 1 hourly precipitation grid is corrected such that, for each grid value, the hourly values summed over 24 hours matches the 24 hourly value.

From optical MODIS (Moderate Resolution Imaging Spectroradiometer) satellite scenes, a time series of SCA has been retrieved for elevation bands for the study catchments. Each 500 × 500 meter grid cell is assigned a SCA value between 0 and 100% snow 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 (Salomonson & Appel 2004). Each catchment has 10 elevation zones of equal area. In total 691 satellite scenes over the period 2014–2015 are used which give us 7,601 data points (i.e. catchment mean and 10 elevation zones). Only scenes with less than 10% cloud cover are used.

SWE data are collected routinely at the peak of the accumulation season and sometimes also at mid-winter by hydropower companies. The usual procedure is to measure snow depth every 10–20 meters at a snow course of length 1–2 km, at a certain elevation, and snow density at 2 or 3 representative points (e.g. where snow depth is close to the average). For this study we use the average SWE value measured in the actual, or in a nearby catchment. 193 such averaged SWE values, collected during the period 2014–2016 are used in this study.

The hydrological models are calibrated against hourly discharge data provided by the HYDRAII database of NVE. The measured water stage is transformed to runoff using estimated stage-runoff relationships. Many of the catchments have problems of ice clogging the cross-section where the water stage is measured during winter. The observed daily runoff is corrected manually if suspicious values due to ice clogging appear. The hourly values are then corrected automatically so that hourly values summed over 24 hours match the corrected 24-hour value. This production chain produces uncertainty at every step, so that winter runoff at mountainous catchments in Norway must be treated with caution. In addition to precipitation and temperature, the Crocus model needs additional forcing data (relative humidity, incoming short wave and long wave radiation, wind speed and wind direction). Forcing for these variables were derived from the numerical weather prediction model AROME MetCoOp (AROME) used operationally by MET Norway (Luijting et al. 2017).

The model comparison is carried out for catchments located in central- and mid-Norway. 17 catchments are evaluated for the models ability to simulate SCA, whereas a subset of 10 catchments are evaluated for SWE. The gridded models, SeNorge and Crocus, simulate regionally or nationally for areas which include these catchments, whereas DDD_EB and DDD_CX simulate for the catchments and their elevation zones. Figure 1 shows the location of the catchments and Table 1 provides geographical information in addition to average annual maximum SWE. The average annual maximum SWE is the simulated catchment value over 58 years as simulated by the HBV model (Hydrologiska Byråns Vattenbalans model; Bergström 1992) which is, together with DDD, used at the Norwegian flood forecasting service at NVE. The HBV model is calibrated against daily runoff and precipitation is adjusted too so that the volume error is minimized (no bias). The values of average annual maximum SWE are hence not necessarily the truth, but indicate the snow regime associated with each catchment.

All models use a threshold temperature of to separate between solid and liquid precipitation. When applicable (for SeNorge and DDD_CX), the threshold temperature for snowmelt is set equal to ⁠. All models use the hourly gridded observations of precipitation and temperature as input data. SeNorge aggregated the 1 hourly data to 3 hourly data. This, because SeNorge is an operational tool providing gridded snow information for all of Norway and it is not set up to run on a 1 hourly time interval since such a fine temporal resolution does not, at the moment, meet any operational needs.

During the calibration of the hydrological model DDD_CX, there was no possibility to adjust precipitation and temperature in order to improve calibration results or to better close the water balance. In using such a set-up for comparison, we believe that we can focus the comparison on the snow algorithms of the different models and better diagnose the results. In addition, such a set-up mimics the situation faced when we need to estimate the snow conditions at ungauged sites and for a changed climate where only scenarios of meteorological forcing is at hand.

With imperfect input, such as underestimated precipitation, we expect the models to underestimate SWE and possibly also SCA. For DDD_CX this is hazardous since through calibration, all model parameters will be optimized in order to try to make up for the deficiency provided by the underestimation of precipitation. The number of parameters to be calibrated against runoff for DDD_CX is six. Table 2 show the calibration parameters of DDD_CX in addition to parameters common to all the models and we note that two of the calibration parameters are associated with snow, and ⁠, but through the calibration process these parameters may have attained values not necessarily associated with snow processes.

DDD_CX was calibrated using the Particle Swarm Optimization and Dimensioned Search (PPSO) method (https://www.rforge.net/ppso/) for the period 1st September 2014–31st December 2015. As can be seen from Table 2, DDD_EB uses model parameters calibrated with DDD_CX, except, for the parameter CX which obviously does not apply when using the energy-balance approach for snowmelt.

Table 3 shows the result for the calibration of DDD_CX and comparison with DDD_EB. Crocus and SeNorge are dedicated snow models and are therefore not validated against runoff. We can note from the mean values of the statistics that the two models give very similar results with the Kling-Gupta efficiency criterion (KGE; Kling et al. 2012) and bias. The estimate of bias, the ratio between simulated and observed runoff shows that, with a few exceptions, the meteorological grids seem to underestimate precipitation with, on average, nearly 20%. This finding corresponds with those of Lussana et al. (2018) as reported above. The skill scores are, in general, not very high, although 5(6) out of the 17 catchments have KGE values over 0.8 for DDD_CX (DDD_EB). Several reasons for this can be put forward and lack of ability to correct the precipitation input to close the water balance together with high uncertainty in the runoff data are the two most prominent candidates.

From the 691 scenes of satellite derived SCA, we could abstract between 25 (133.7) up to 53 (48.5) values of SCA for each catchment and their elevation zones for the two snow seasons 2013/14 and 2014/15. Figure 5 shows the MAE for SCA organized by elevation and by catchment. We can see that the errors in estimating SCA are not sensitive to elevation. What is striking is the clear difference in performance between Crocus and DDD_CX on the one hand and SeNorge and DDD_EB on the other hand. The two latter models have very similar performance with MAE quite close to zero and low RMSE (figure not shown). Figures 6 and 7 show simulated and observed catchment values of SCA for the two catchments 48.5 and 122.9. From Figure 6 (catchment 48.5) we see that both DDD_EB and SeNorge simulates the catchment SCA very well. DDD_CX underestimates the SCA in 2014 but has a similar performance to DDD_EB in 2015. For both years Crocus has a too steep decline in SCA and starts the melting too early. This latter feature is common for all the models except DDD_EB. From Figure 7 (catchment 122.9) we see that the observed SCA has a more varied behaviour during the snow season which is not captured very well by the models. DDD_EB is overestimating SCA, whereas SeNorge, Crocus and DDD_CX underestimate SCA in the melting period of 2014. All models capture the timing of melt-out reasonably well.

Table 4 shows the overall results of validating the models against SWE and SCA. For both validations the SeNorge and DDD_EB emerge as the best models. We see for example that on average, the error of SeNorge for SWE is 31% less than that of DDD_CX.

The aim of this study is to learn what principles in snow modelling model are best suited for operational snow modelling under a changed climate and for predictions in ungauged basins. When validating the models against observed SWE and SCA, SeNorge and DDD_EB, quite convincingly stand out as best suited, and especially so for the validation against observed SCA. These two models represent the most commonly used snow melt algorithms; the temperature index (although extended with a radiation term) and the energy balance approach (i.e. Equations (6) and (2)). The inferior models, DDD_CX and Crocus also use these snow melt algorithms, so, clearly, the choice of empirical/physically based approach has not been the key issue for differences in results, but more on how these methods are applied. We have set rather strict rules for the model comparison set-up, with the intention of isolating the workings of the models' snow algorithms. In addition, such a set-up is intended to be as close to situations of predicting for ungauged basins and for effects of climate change as possible.

Crocus is quite demanding in respect of forcing variables. In this study, these forcing variables are abstracted from AROME. One could argue that using the interpolated meteorological grids of observed precipitation and temperature instead of the physically consistent precipitation and temperature from AROME would provide better results. However, in a recent study by Luijting et al. (2017), simulated snow depths were compared against observed and it was quite clearly shown that using Crocus with the interpolated meteorological grids performed much better than with only using forcing from AROME. In this study, Crocus was originally intended to be the benchmark model, the role Crocus has had in many model comparison studies (Strasser et al. 2002; Essery et al. 2013). However, this study differs from previous ones in that we have run Crocus regionally and conducted the validation on spatially averaged values from elevation zones of catchments and entire catchments, and not at points as has been customary. The spatial distribution of snow over an area has been the focus of several investigations and is seen as crucial for the dynamics of snowmelt, the evolution of SCA and for properly accounting for the energy fluxes in land-surface schemes (Essery & Pomeroy 2004; Helbig et al. 2015). Crocus, being a point model, does not aim especially at representing a realistic spatial distribution of snow, which is in contrast to the other models.

It may be instructive to look at the behaviour of the models in more detail in order to better understand why a well-established model such as Crocus shows deficiencies in simulating at a regional scale. Figure 8 shows accumulated melt rate for the catchment 122.9 for the winter 2013/14 (Figure 8(a)) and the winter 2014/15 (Figure 8(b)).

We see from Figure 8 that Crocus tends to start melting early for both years, and has finished melting earlier than all other models. In order to further investigate the models' discrepancies we have, in Figures 9 and 10, plotted net radiation and albedo as simulated by Crocus and DDD_EB (DDD_CX and SeNorge do not simulate such elements being temperature-index based models). Figure 9 shows that the general pattern of net radiation for DDD_EB and Crocus is quite similar except that the net radiation of Crocus increases much earlier than DDD_EB. This can be explained by the differences in albedo shown in Figure 10 and suggests that the difference in net radiation is caused by increased SW radiation of Crocus. The loss in albedo of Crocus is consistent with the higher fraction of snow-free areas found early in the season for the Crocus simulations. We also note that, in general, the albedo of Crocus is lower than that of DDD_EB which suggests increased melting. At this point, we do not find any obvious reasons as to why the melt starts so early for Crocus and it is beyond the scope of this study to investigate and perhaps trace errors to the AROME input. This study represents a first attempt at regional simulations by Crocus in Norway, and further investigations will have to address questions such as the need for a more detailed land-surface description and the scaling up of model parameters. Figure 10 also shows that the albedo of DDD_EB displays some strange temporal fluctuations. The albedo seems to fluctuate mainly between 0.8–0.87. DDD_EB uses the albedo algorithms of the Utah Energy balance Snow Accumulation and Melt Model (Tarboton & Luce 1996), and clearly the implementation needs adjustment. The most probable cause is associated with the illumination angle of radiation. The effects of the errors, however, are not large enough to seriously affect the results of SWE and SCA simulations.

The results of DDD_CX present a very convincing case of the problems encountered in calibrating a parameter (CX) as a member of a set of parameters that simultaneously are to be calibrated against a variable (i.e. runoff). The relationship between parameter and variable may vary over time and, in addition, non-physical combinations of parameters may work together to produce good results for the variable. The result can be that we may have sensible, reliable results as well as results where, for example, the simulation of snow variables have been sacrificed in order to provide good results for runoff. In Figures 2 and 5 we see that DDD_CX quite often show the poorest performance for simulating SWE and SCA, and the summary of results in Table 4 points at DDD_CX as the inferior model for simulating both SWE and SCA. Figure 4 shows how SWE is consistently underestimated by DDD_CX compared to the other models and to observations, whereas in Figure 3, the DDD_CX simulations do not differ much from the other models. The process parameterization itself (temperature-index) appears to be quite capable of providing good results for the simulation of snowmelt (see also the results of the SeNorge model), but the procedure of calibrating the temperature-index (as part of a parameter set) against runoff instead of against snow observations is hazardous and cannot be recommended. It is a paradox that such a calibration set-up has been the usual practice in operational hydrology in the Nordic countries for decades, although a wealth of snow measurements exists which could be used for calibration of the temperature-index directly.

DDD_EB shows surprisingly robust results across the 17 (10) catchments for SCA and SWE, when we take into consideration that only precipitation and temperature is used to force the energy balance calculations. The catchments represent a variety of climatic regimes of Norway from relatively dry inland catchments (2.268 and 12. 171), moist lowland catchments (122.9, 133.7) and to alpine and moist (located on the west coast) catchments. It is promising, with respect to the overall aim of providing model structures useful for predicting in ungauged basins and for a changed climate, that the proxy models used to calculate the energy balance seem to work well for such different areas. We can also note, from Table 3 that there is no significant difference in skill for simulating runoff between DDD_CX and DDD_EB, indicating that the introduction of the proxy models for energy balance calculations, and the associated loss of one calibration parameter (CX), did not deteriorate the runoff predictions of the DDD model.

The SeNorge model shows the best performance of all the models. The fact that it is run at a 3-hourly resolution instead of 1-hourly as the other models, may account for the results being less sensitive to high variability in the meteorological input (precipitation and temperature) but the validation against data for which it has not been calibrated shows that the model is able to capture the essential snow characteristics at different regions and in different climates in Norway. The results of SeNorge also demonstrate that simple empirical relations such as the extended temperature index method work well and can be used regionally if calibrated against relevant data (snow melt measurements).

The Crocus model is primarily an alpine model and needs additional forcing data for simulating snow. In this study, Crocus simulates regionally over many different landscape types and the additional forcing, provided by the numerical weather prediction model AROME, may not be of sufficient quality to give us the full potential of the Crocus model. The temperature-index approach in the DDD_CX model represents a very common approach of simulating snowmelt in operational hydrological models. It has been clearly demonstrated in this study that calibrating the temperature index factor for snowmelt against runoff is not optimal for simulating the coverage and dynamics of SWE. The proxy models for estimating the energy balance in DDD_EB have shown quite acceptable results although the approach is coarse, i.e. neglecting spatial variability in topography and landscape types and temporal variability in wind speed and air pressure. There might be dependencies between the choice of the constants in the in the model and the chosen temporal resolution (1 h) that have to be investigated. The SeNorge model applies the extended temperature-index method which includes a radiation term. The snow melt algorithm is calibrated against observed snow melt rates which is both reasonable and gives good results. In addition, the models discriminate between areas below and above the tree-line. The temperature-index factor is dependent on the chosen temporal resolutions and needs, in principle, recalibration for other resolutions. There is a potential for improving the model by applying a more realistic spatial distribution of snow.

The overarching aim of this study is to provide suggestions for snow model structures suited for predicting snow conditions at ungauged sites and under a changed climate. Common to the results of several other such model comparisons we find that well-established empirical parameterizations give as good (and poor) results as more physically based relationships (Essery et al. 2013). Some quite clear conclusions can, however, be drawn from this exercise.

Empirical relationships, such as the temperature index method, work well, also regionally if the temperature index is calibrated against the variable it is supposed to represent, namely the melting of snow as a function of temperature. The SeNorge model has shown itself to be the best model when validated against SWE and SCA. The temperature-index based method calibrated against runoff (DDD_CX), had the poorest results when validated against SWE and SCA.

Proxy modelling of the energy balance to simulate snowmelt (DDD_EB) works surprisingly well, and almost as good as SeNorge when validated against SWE and SCA. It is especially promising with respect to predictions in ungauged basins that no extra tuning was performed for the different catchments. The only difference in parameters from one catchment to the next (besides those calibrated with DDD_CX) was that the location of the centre of the catchment was provided so that the SW calculations could be performed.

Both SeNorge and DDD_EB show promise as candidates for model structures suited for simulating snow for ungauged basins and for a changed climate. The model parameters in the snow module of DDD_EB need no calibration, and the consistently good results of SeNorge over all the catchment indicate that the calibration of the extended temperature index melt model is physically meaningful and thereby possibly independent of climate. Running Crocus regionally, probably need some scheme for realistically representing the spatial distribution of SWE. The temperature-index method for snowmelt calibrated against runoff, cannot be recommended.

The support of the Norwegian Research Council, Glommen and Laagen Brukseierforening, E-CO, Trøndelag Energiverk and HYDRO Energi through the project ‘Better SNOW models for predictions of natural Hazards and HydropOWer applications’ (SNOWHOW) (244153/E10), is gratefully acknowledged.