# Revisiting a Simple Degree-Day Model for Integrating Satellite Data: Implementation of Swe-Sca Hystereses

Philippe Riboust, Guillaume Thirel, Nicolas Le Moine, Pierre Ribstein

► To cite this version:

Philippe Riboust, Guillaume Thirel, Nicolas Le Moine, Pierre Ribstein. Revisiting a Simple Degree-Day Model for Integrating Satellite Data: Implementation of Swe-Sca Hystereses. Journal of Hydrology and Hydromechanics, 2019, 67 (1), pp.70 - 81. 10.2478/johh-2018-0004. hal-01935399

HAL Id: hal-01935399

https://hal.sorbonne-universite.fr/hal-01935399

Submitted on 26 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Revisiting a simple degree-day model for integrating satellite data: implementation of SWE-SCA hystereses

Philippe Riboust<fnref n="1" />,<fnref n="2" />*, Guillaume Thirel<fnref n="2" />, Nicolas Le Moine<fnref n="1" />, Pierre Ribstein<fnref n="1" />

**Abstract:** Conceptual degree-day snow models are often calibrated using runoff observations. This makes the snow models dependent on the rainfall-runoff model they are coupled with. Numerous studies have shown that using Snow Cover Area (SCA) remote sensing observation from MODIS satellites helps to better constrain parameters. The objective of this study was to calibrate the CemaNeige degree-day snow model with SCA and runoff observations. In order to calibrate the snow model with SCA observations, the original CemaNeige SCA formulation was revisited to take into account the hysteresis that exists between SCA and the snow water equivalent (SWE) during the accumulation and melt phases. Several parametrizations of the hysteresis between SWE and SCA were taken from land surface model literature. We showed that they improve the performances of SCA simulation without degrading the river runoff simulation. With this improvement, a new calibration method of the snow model was developed using jointly SCA and runoff observations. Further analysis showed that the CemaNeige calibrated parameter sets are more robust for simulating independent periods than parameter sets obtained from discharge calibration only. Calibrating the snow model using only SCA data gave mixed results, with similar performances as using median parameters from all watersheds calibration.

**Keywords:** Snow model; Hysteresis parametrization; MODIS snow cover area; Rainfall-runoff model.

## INTRODUCTION

Snow accumulation in winter as well as spring snowmelt gives to mountain catchments a particular hydrological response that should be taken into account when modelling river runoff. Rainfall-runoff models are often coupled with snow models so that solid precipitation storage in snow reservoirs can be represented and released as river runoff when this accumulated snow melts. Improving the snow dynamics modelling is a necessity for improving hydrological forecasting for snow-dominated catchments. Snow-dominated areas are very sensible to temperature changes; therefore climate change can drastically impact the hydrological cycle (Beniston et al., 2017; Bernsteinová et al., 2015) and the population living in these regions (one sixth of the world, estimated by Barnett et al., 2005). For studying the impacts of climate change in mountains and other snow-dominated region, the robustness and reliability of models have to be assessed (and improved) in order to ensure that the model structure and parameters are still valid for simulating river runoff in a different climate context (Thirel et al., 2015a, 2015b).

To improve the robustness of degree-day empirical snow models, numbers of studies have demonstrated the advantages of snow cover data in addition to the runoff observations to calibrate and validate models (e.g. Franz and Karsten, 2013; Parajka and Blöschl, 2008a) or using data assimilation techniques (Andreadis and Lettenmaier, 2006; Rodell and Houser, 2004; Slater and Clark, 2006; Thirel et al., 2013). Snow models are traditionally calibrated only with runoff observations, making the obtained parameter sets dependent on the rainfall-runoff model used. Adding snow data information to the calibration procedure allows the snow model to be less dependent on the hydrological model, and therefore allows to obtain more robust snow parameter sets (Franz and Karsten, 2013). In addition to the improvement of the robustness of the snow model, this procedure also helps to better identify its structure and parameters (Grayson et al., 2002; He et al., 2014; Parajka and Blöschl,

2008b, 2008a). Most of these studies used satellite Snow Cover Area (SCA) from the MODIS sensor (Hall et al., 2006), which provides a spatialized vision of the presence of snow on the ground at a 500m resolution. The main disadvantage of MODIS satellite data is that it cannot provide SCA data when there is cloud cover. Several methods have been developed to compensate for these deficiencies, most of which are based on a spatiotemporal filter, i.e., replacing the missing pixel value with its nearest neighbors, whether it be spatial or temporal (Da Ronco and De Michele, 2014; eg. Parajka and Blöschl, 2008b; Poggio et al., 2012). Another approach, which can be additionally used for filling lacunar MODIS data, is to determine the regional snowline elevation (Gafurov and Bárdossy, 2009; Krajčí et al., 2014, 2016; Parajka et al., 2010). SCA data give only geographic information about presence or absence of snow without giving any information about the quantity of Snow Water Equivalent (SWE) present on the watershed. SWE would be more valuable than SCA for snow models, but unfortunately the accuracy of remote sensing SWE data is too low to be used effectively (Vuyovich et al., 2014) and point scale measurements are scarce and difficult to interpolate (Parajka et al., 2012).

SCA parametrization within snow models is highly dependent on the resolution selected. When the model represents an area with a heterogeneous topography, the relation between the SCA and the quantity of SWE averaged over the area forms a hysteresis (Luce and Tarboton, 2004). During an accumulation period, the SCA increases very quickly and remains stable, whereas the SWE increases slowly. On the opposite, during a melting period, the decrease in the SCA is slower and more gradual. The relation between SWE and SCA therefore differs depending on whether a snowmelt period or a snow accumulation period is considered. Snow accumulates heterogeneously in the watershed due to redistribution of snow by wind drift (Esery and Pomeroy, 2004). During the melting period, snow tends to melt uniformly with variations depending on the mountain slope and orientation (Egli and Jonas, 2009). Snow patches

<footnote n="1">¹ Sorbonne Universités, UPMC Univ., Paris 06, CNRS, EPHE, UMR 7619 Metis, 4 place Jussieu, 75005 Paris, France.</footnote>

<footnote n="2">² Hydrosystems and Bioprocesses Research Unit (HBAN), Irstea, 1 rue Pierre-Gilles de Gennes, CS 10030, 92761 Antony Cedex, France.</footnote>

<footnote>* Corresponding author. E-mail: philippe.riboust@irstea.fr</footnote>
tend to be created during the melting period, due to the combi-
nation of heterogeneous snow depths and to heterogeneous melt
rates at watershed scale (Egli and Jonas, 2009; Helbig et al.,
2015; Magand et al., 2014). This representation of the snow
cover area is an issue that has greatly interested the Land Sur-
face Modelling (LSM) community. As snow's albedo has a
strong retroaction on the earth-atmosphere coupling because of
its influence on the energy budget, precise knowledge of the
SCA is necessary (Zaitchik and Rodell, 2009). This is why
LSMs often integrate a parametrization of hysteresis, linking
SCA to SWE, whose complexity can vary (Liston, 2004; Ma-
gand et al., 2014; Nitta et al., 2014; Niu and Yang, 2007; Swen-
son and Lawrence, 2012). Parametrizations of the hysteresis
have been introduced in a few snow models (Clark et al., 2011;
Duethmann et al., 2014; Franz and Karsten, 2013; Kolberg and
Gottschalk, 2006; Magnusson et al., 2014), as the hysteresis
may be even found at small geographical scales (Luce and
Tarboton, 2004).

To more accurately represent river runoff in snow-covered catchments using lumped GR rainfall-runoff models (Le Moine et al., 2007; Perrin et al., 2003; Pushpalatha et al., 2011), the CemaNeige snow-accounting routine (Valéry et al., 2014a, 2014b) was developed. The present study is based on this model, which is a degree-day model whose primary objective is to assist the hydrological model in accurately simulating mountain river runoff as parsimoniously as possible, i.e., with the minimum of free parameters to calibrate. Once its development had been completed, the number of parameters necessary for the snow module was limited to two and required only air temperature (mean temperature or minimum and maximum daily temperatures) and daily precipitation as forcing variables. The CemaNeige snow model has been tested on French, Swiss, Swedish and Canadian watersheds (Valéry et al., 2014b), and an intercomparison of multiple snow-accounting models coupled with different rainfall-runoff models has been made (Troin et al., 2016, 2015). This study showed that the CemaNeige snow-accounting model gave good performances in the Québec area when coupled with any hydrological model. In other types of climate, the CemaNeige snow model has shown some limitations. In Nepal (Pokhrel et al., 2014) the model performances were limited, certainly due to the fact that glacier melt and sublimation processes are not yet represented in CemaNeige. In the dry Andes, the CemaNeige model had to be adapted in order to take into account the sublimation process which is a major cause of snow ablation in arid mountain catchments (Hublart et al., 2016).

The objective of the study is to improve the robustness of the
CemaNeige snow-accounting model when used with a rainfall-
runoff model. For this, the use of MODIS SCA data for calibra-
tion was investigated and diverse weights in the objective func-
tion, between SCA and discharge, were tested. In order to make
full use of MODIS SCA data at the catchment scale, SWE-SCA
hysteresis were implemented in the snow model and assessed.

DATA
SAFRAN and MODIS data

The forcing data used by the GR4J model with CemaNeige
comes from SAFRAN data (Quintana-Seguí et al., 2008; Vidal
et al., 2010). SAFRAN is an atmospheric reanalysis model: it
assimilates surface observations to calculate the energy fluxes
as well as temperature, precipitation, humidity, and wind on a
regular 64-km² grid. For the needs of this study, only SAFRAN
temperature and precipitation data were used, after they were
aggregated by catchment on the daily time step.

The MODIS snow cover product (Hall et al., 2006) came
from the National Snow and Ice Data Center (NSIDC,
http://nsidc.org). The first measurements were taken in summer
2000 by the Terra satellite, followed 2 years later by measure-
ments from the Aqua satellite. The MODIS instruments meas-
ure visible radiance. A post-treatment deduces binary snow
cover data on 500-m resolution grids once or twice a day.

For this study, the MODIS data were aggregated over differ-
ent elevation zones on each catchment so that they would be
comparable to the CemaNeige simulations. This aggregation
makes it possible to obtain a snow fraction for each band. The
quality of the snow fraction provided by MODIS is considered
unsatisfactory if the proportion of pixels covered by clouds is
greater than 40% of the total number of pixels of the elevation
zone considered. If this threshold is not satisfied, these data are
not taken into account for the calibration and evaluation proce-
dures.

Selection of catchments

The catchments for this study were selected based on sever-
al criteria:

* data availability

* sufficient catchment snow cover

* low catchment’s anthropization

The availability of runoff data and MODIS data is a limiting factor. The runoff data used came from the Hydro database (http://www.hydro.eaufrance.fr/); they were available at the daily time step and extracted up to the end of 2010. In addition to this, the MODIS data begins in 2000. We decided to preserve the catchments whose runoff data contained a maximum of 1 year of missing data for the 2000–2005 period and the 2005–2010 period. Moreover, only the catchments whose runoff was nearly natural, for example not influenced by dams, were retained. This limits the number of catchments available in high-altitude mountainous areas, which are often developed for hydroelectric power production.

To ensure that the watersheds were frequently snow-
covered, we decided to retain only those whose outlet was
located above 300 m and whose annual snow cover lasted a
minimum 30 days on the catchment's highest elevation zone
(i.e. 20% of the watershed area) according to MODIS data.
Finally, 277 catchments fulfilled these conditions and were
used for this study (Figure 1).

METHODS

Description of the CemaNeige snow-accounting model

The CemaNeige snow-accounting model (Figure 2) is a degree-day snow model with two free parameters with a simple representation of the catchment’s SCA (Valéry et al., 2014b). The snowmelt generated by the snow module is considered as liquid precipitation in the GR4J rainfall-runoff model. To take into account the effects of the relief on precipitation, the CemaNeige snow-accounting model is distributed into elevation zones with equal surface area, usually five bands. Temperature and precipitation are extrapolated for each elevation zone (Valéry et al., 2014a). Both the basic version of CemaNeige and the GR4J model come from the R airGR package (Coron et al., 2017a, 2017b).

The CemaNeige model determines the fraction of solid pre-
cipitation that accumulates when temperature conditions are
adapted. To take into account the thermal inertia of the snow-
pack, melting is delayed by a formulation taking the air tempera-
<figure><img src="image_2.png" /><figcaption>Fig. 1. Map of the 277 watersheds used in this study. The color of the basins represents the mean annual solid precipitation.</figcaption></figure>

ture of the preceding days into account (Eq. 1), modulated by one parameter (cT). The snowpack’s cold content (eT<sub>G</sub>) is calculated based on weighting between the value of the internal variable eT<sub>G</sub> of the preceding time step and the air temperature of the day considered.

$$eT_{G_t} = \min(cT \cdot eT_{G_{t-1}} + (1-cT)T_{mean_t}, 0) \quad (1)$$

When the snow mantel’s temperature simulation (eT<sub>G</sub>) reaches a value of 0 and the mean air temperature (T<sub>mean</sub>) is greater than 0°C, potential melting (Melting<sub>pot</sub>) is calculated, which is dependent on the melting factor Kf (Eq. 2).

$$Melting_{pot_t} = Kf \cdot T_{mean_t} \quad (2)$$

This potential melting is modulated depending on the overall snow surface by the following formula:

$$Melting_t = (0.9 \cdot SCA_t + 0.1) \cdot Melting_{pot_t} \quad (3)$$

For a 100% SCA, melting will be identical to potential melting. The more the SCA decreases, the more the melting speed decreases, to a minimum of 10% of the potential melting. The SCA is calculated in the model as follows:

$$SCA_t = \min \left( \frac{SWE_t}{T_{h,melt}}, 1 \right) \quad (4)$$

where the SWE is the quantity of snow accumulated on the catchment in snow water equivalent (a state variable of the model, in mm) and where T<sub>h,melt</sub> is the model’s melting threshold. This T<sub>h,melt</sub> melting threshold is calculated as being equal to 90% of mean annual solid precipitation on the catchment considered (Valéry et al., 2014b).

This melting modulation stems from the fact that preferential melting is usually established based on the mountain’s slope exposure: shaded slopes tend to melt less quickly than sunny slopes. The conceptual model therefore includes a state variable close to the SCA. However, this formulation of the SCA does not necessarily give values that can be directly compared to the MODIS observations because this melting formulation only takes into account that the relation between the SWE and the SCA is not unequivocal and differs depending on whether an

<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic representation of the CemaNeige snow model (Modified from Valéry et al., 2014b).</figcaption></figure>

accumulation or a melting period is considered. To take this into account, we suggest a hysteresis that allows simulating a fast increase of the SCA at the accumulation period and a smooth decrease at the ablation period. The fast increase at the accumulation period is due to a homogeneous increase of the SCA with snowfall (although the redistribution of snow creates high disparities in snow depths). The heterogeneous SWE coupled to heterogeneous melt rates due to the topography tends to create patchy snow areas (Liston, 2004; Magand et al., 2014).

## Hysteresis

To take this melting heterogeneity into account, several hystereses were implemented and tested in the CemaNeige snow-accounting routine. The hystereses selected in the literature are presented in this section.

## Accumulation curve

The snow cover curve dependent on the SWE in case of accumulation was formulated with a linear function dependent on an accumulation threshold T<sub>h,acc</sub>. This follows the following relation:
$$ \Delta SWE_t = P_{sol_t} - Melting_t \qquad (5) $$

$$ SCA_t = \min \left( SCA_{t-1} + \frac{P_{sol_t}}{T_{h, acc}}, 1 \right); \text{ if } \Delta SWE_t \ge 0 \qquad (6) $$

where $\Delta SWE$ represents the variation of SWE at the time step considered (a positive $\Delta SWE$ corresponds to a snow accumulation, a negative $\Delta SWE$ to melting), $SCA_{t-1}$ corresponds to the SCA calculated at the preceding time step and $P_{sol_t}$ corresponds to the quantity of snowfall at current time step (in mm). This simple formulation was used in the CLSM model (Magand et al., 2014). The slope of this relation is set by the parameter $T_{h,acc}$ which defines the threshold of SWE (in mm) based from which the SCA is equal to 1. This formulation of the SCA in an accumulation period is used in this article for all the hysteresis models tested.

## Melting curve

The different formulations tested for this study are presented in Figure 3; Figure 3a presents the melting curve as it is implemented in the current CemaNeige snow-accounting model ((Valéry et al., 2014b), Eq. 4).

### Linear hysteresis (LH)

The first hysteresis tested (Fig. 3b) is a simple linear hysteresis that was implemented in the CLSM model (Magand et al., 2014).

$$ \left\{ \begin{array}{ll} SCA_t = \min \left( \frac{SWE_t}{T_{h,max}}, 1 \right) & ; \text{ if } \Delta SWE_t < 0 \\ SWE_{t,melt} = SWE_{t-1} & \text{if } SCA_{t-1} = 1 \end{array} \right. \qquad (7) $$

$$ T_{h,max} = \begin{cases} T_{h,melt} & \text{if } SWE_{t,melt} > T_{h,melt} \\ SWE_{t,melt} & \text{if } SWE_{t,melt} \le T_{h,melt} \end{cases} \qquad (8) $$

This parametrization of the hysteresis melting curve part depends on a local maximum threshold ($T_{h,max}$), which defines the inflexion point and the slope of the melting curve (in red on Fig. 3b). If the accumulated snow before snowmelt exceeds a calibrated melt threshold ($T_{h,melt}$), the local threshold ($T_{h,max}$) takes the ($T_{h,melt}$) value. On Fig. 3b, this means that when melting occurs, the SCA value decreases only if the SWE value is lower than 250 mm (this value has been taken as an example). If during the winter, the accumulated SWE never exceeds the $T_{h,melt}$ value, or in the case that there is a snow accumulation after the melting phase has started (represented by the yellow curve in Fig. 3b), the local maximal threshold ($T_{h,max}$) takes the maximal SWE value before the beginning of the melt ($SWE_{t,melt}$, Eq. 8).

### Modified linear hysteresis (LH*)

A variation of the linear hysteresis is proposed in this article, where:

$$ T_{h,melt} = \overline{P_{sol annual}} \cdot R_{sp} \qquad (9) $$

$P_{sol annual}$ represents the mean annual precipitation for the elevation zone of the catchment considered and $R_{sp}$ is a parame-

ter between 0 and 1. This formulation defines a different $T_{h,melt}$ for each of the catchment's elevation zones, which is not the case in Equation 7.

### MATSIRO hysteresis (MH)

This third hysteresis tested was used within the MATSIRO surface model (Nitta et al., 2014). The formulation is similar to the linear hysteresis formulation, except that instead of following a linear slope, the hysteresis follows a square root function for the melting curve.

$$ \left\{ \begin{array}{ll} SCA_t = \min \left( \sqrt{\frac{SWE_t}{T_{h,max}}}, 1 \right) & ; \text{if } \Delta SWE_t < 0 \\ SWE_{t,melt} = SWE_{t-1} & \text{if } SCA_{t-1} = 1 \end{array} \right. \qquad (10) $$

$$ T_{h,max} = \begin{cases} T_{h,melt} & \text{if } SWE_{t,melt} > T_{h,melt} \\ SWE_{t,melt} & \text{if } SWE_{t,melt} \le T_{h,melt} \end{cases} \qquad (11) $$

This parametrization also calls on the $T_{h,melt}$ parameter. As equation 8, the value of $T_{h,max}$ depends on the state of the accumulated SWE at the beginning of the melting period ($SWE_{t,melt}$) and of the value of the $T_{h,melt}$ parameter. This curve shape implies that the decrease in the SCA accelerates as the SWE decreases. In other words, this formulation would mean that at the beginning of melting, a decrease in the SWE affects the south-exposed snow cover only slightly. As the quantity of snow diminishes, melting affects the entire catchment and the SCA drops rapidly. The shape of this hysteresis can be consulted in Figure 3c. Like the linear hysteresis, Equation 10 updates the melting threshold with the SWE value if the SWE accumulated does not go beyond the $T_{h,melt}$ value or in the case of a snow accumulation within a snow melt period.

### Swenson hysteresis (SH)

The last parametrization tested comes from the CLM4 model (Swenson and Lawrence, 2012). It is very different from the other methods since it does not use the melting threshold to manage melting. Instead, a $N_{melt}$ parameter directly imposes the shape of the melting curve (Equations 12 and 13).

$$ SCA_t = 1 - \left( \frac{1}{\pi} \arcsin \left( 2 \frac{SWE_t}{SWE_{t,melt}} - 1 \right) \right)^{N_{melt}} ; \text{ if } \Delta SWE_t < 0 \qquad (12) $$

$$ SWE_{t,melt} = SWE_t \left[ \frac{\cos(\pi(1-SCA_{t-1})) + 1}{2} \right]^{-1/N_{melt}} \qquad (13) $$

The higher the $N_{melt}$ value, the lower the decrease in the SCA will be at the beginning of melting and it will accelerate progressively, which makes it possible to model the threshold effect as it is in the hysteresis presented above. In contrast, the lower the $N_{melt}$ value, the more the SCA will drop rapidly at the beginning of melting and will be followed by a relative stabilization of the slope until a second, sharper drop at the end of melting. Conceptually, this would mean that the SCA of the exposed slopes would drop very rapidly at the beginning of melting, whereas the SCA of the less exposed slopes would decrease less rapidly. The influence of the melting factor $N_{melt}$ on the SCA during the melting period is presented in Figure 3d, 3e, and 3f.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Schematic representation of the different types of hystereses used a) original relation in CemaNeige, b) linear hysteresis, c) MATSIRO hysteresis, d-e-f) Swenson hysteresis with three different $N_{melt}$ coefficients. The blue line represents the accumulation phase of the hysteresis while the red line represents the melting phase. The yellow curves represent an example of snow accumulation within a melting period. The $T_{h,acc}$ and $T_{h,melt}$ values have been set respectively at 20 and 250 mm to illustrate the different hystereses.</figcaption></figure>

## Calibration and optimization criteria

The four different hysteresis models implemented in the CemaNeige snow-accounting model were tested and calibrated on the 277 catchments. The calibrations were performed for each catchment for two different periods: 2000–2005 and 2005–2010. The routing and rainfall-runoff transformation were computed using the GR4J hydrological model.

Both the hydrological model's and the snow model's parameters were calibrated using the steepest-descent calibration algorithm (Edijatno et al., 1999; Perrin et al., 2001). The calibration criterion used was the KGE' (Gupta et al., 2009; Kling et al., 2012). This criterion (Eq. 14) takes into account the Pearson correlation coefficient (Eq. 15, $Cov_{so}$ being the covariance between observation and simulation, $\sigma$ being the standard deviation), the percentage bias (Eq. 16, depending on $\mu$ the average of the simulation and observations), and the ratio of the coefficients of variation between the simulated and observed temporal series (Eq. 17).

$$KGE' = 1 - \sqrt{(r-1)^2 + (\omega-1)^2 + (\gamma-1)^2} \quad (14)$$

$$r = \frac{Cov_{so}}{\sigma_s \sigma_o} \quad (15)$$

$$\omega = \frac{\mu_s}{\mu_o} \quad (16)$$

$$\gamma = \frac{\sigma_s / \mu_s}{\sigma_o / \mu_o} \quad (17)$$

To calibrate the model with the runoff and the MODIS data, the KGE' on the runoff and the KGE' on the MODIS data on five elevation zones were combined.

$$\text{Crit} = \alpha \text{KGE}'(Q) + \sum_{i=1}^{5} \beta_i \text{KGE}'(\text{SCA}_i), \quad (18)$$

$$\text{where } \alpha + \sum_{i=1}^{5} \beta_i = 1 \quad (19)$$

Thus a weighting $\beta_i$ can be assigned to the SCA simulation on each elevation zone i modeled by CemaNeige (here five zones were selected: $SCA_1, SCA_2, ..., SCA_5$). In addition, a single weighting $\beta$ can be assigned on the SCA criterion if the SCA is averaged over the entire watershed (i.e., all the elevation zones).

The performance presented in the Results section is calculated from the validation period complementary to each calibration period (e.g., 2000–2005 calibration, 2005–2010 validation). The performances of simulations are therefore represented twice for each catchment (one performance per validation period). To ensure that the performance criteria take into account the modifications to the snow model, the validation concerns only the months with snow. For each catchment, the months with snow correspond to the months whose upper layer mean SCA is greater than 5% according to MODIS for the whole period.

A Friedman test (Friedman, 1937) is also applied to the calibration results. This is a statistical test designed to determine whether the models give significantly different results. For each catchment, the performances of the simulations of each model are ranked, and then the ranks are summed for each model. The model with the highest ranks is considered the best according to this test.

# RESULTS

## Analysis of the different hysteresis models

The GR4J-CemaNeige coupled model contains six parameters to calibrate (four for GR4J, two for CemaNeige). The three versions of CemaNeige with hysteresis add two additional parameters: $T_{h,acc}$ and $T_{h,melt}$ or $N_{melt}$. Except for the original CemaNeige model that serves as the reference and was calibrated only on runoff ($\alpha = 1$), the other models tested were calibrated with weighting $\alpha = 0.75$ and an identical weighting $\beta$ for all elevation zones equal to 0.05. The different models used as well as their calibration weighting are presented in Table 1. The comparison of the performance of the models is available in Figure 4.
<table><caption>**Table 1.** List of the models and weightings used for analyzing the performances of the hystereses formulations.</caption><thead><tr><th>Model</th><th>Q weighting α</th><th>SCA weighting β 1 ... 5</th></tr></thead><tbody><tr><td>CemaNeige reference (C Q)</td><td>1</td><td>0</td></tr><tr><td>CemaNeige SCA (C Q SCA)</td><td>0.75</td><td>0.05</td></tr><tr><td>Linear hysteresis (LH)</td><td>0.75</td><td>0.05</td></tr><tr><td>Modified linear hysteresis (LH*)</td><td>0.75</td><td>0.05</td></tr><tr><td>MATSIRO hysteresis (MH)</td><td>0.75</td><td>0.05</td></tr><tr><td>Swenson hysteresis (SH)</td><td>0.75</td><td>0.05</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Performances of the different hystereses used with CemaNeige coupled with GR4J model. The performances presented are KGE' values calculated on validation periods. Each boxplot contains the results from both validation periods for each catchment; the whiskers indicate the 5 and 95 percentiles. The result of the Friedman analysis is represented by the letters near the boxplot. The best ranked model is indicated by letter "a" and models with identical letters are considered as not significantly different.</figcaption></figure>

It can be observed that the reference CemaNeige snow-accounting model gave the best scores in terms of runoff, according to the Friedman test, with a similar median to the models with hysteresis (0.81), however. In terms of SCA performance, the original version of CemaNeige gave much lower scores than the versions with hysteresis (median, 0.19 versus 0.79 for both linear hystereses). Calibration of CemaNeige with the MODIS data, as originally designed, showed a decrease in the runoff score of approximately 0.05 points on the median and only improved the SCA score very slightly. This shows the inability of the original CemaNeige model to take into account the calibrated SCA data without damaging the runoff score.

The implementation of hystereses in CemaNeige significantly improved the SCA score, with the best score obtained by the linear hysteresis: median, 0.81 and 0.82 for the modified version. The runoff scores modeled by CemaNeige with linear hysteresis were very slightly inferior to the reference scores. They had the same median and the same quartiles; only the mean, not represented here, differed slightly. Both of the other hystereses showed similar performance in terms of runoff, but slightly inferior in the SCA. It can be concluded that adding a hysteresis significantly improves the SCA if it is used in calibration while retaining a stable performance in runoff, which the original model did not allow.

## Analysis of runoff criterion and snow cover criterion weighting

To obtain better validation criteria on the SCA without harming the runoff criterion, several different weighting configurations were analyzed. The α weighting on runoff was tested for values of 0, 0.5, 0.75, 0.9, and 1. Weighting corresponding to the SCA criteria was 0.2, 0.1, 0.05, 0.02, and 0 for each of the layers. This analysis was based on the use of the modified linear hysteresis, which gave highly satisfactory results in the previous section. These values are summarized in Table 2 and the results are available in Figure 5.

<table><caption>**Table 2.** Weightings used for the calibrations of the snow model with the modified linear hysteresis (LH*).</caption><thead><tr><th>Model</th><th>Q weighting α</th><th>SCA weighting β 1 ... 5</th></tr></thead><tbody><tr><td>LH*</td><td>0</td><td>0.2</td></tr><tr><td>LH*</td><td>0.5</td><td>0.1</td></tr><tr><td>LH*</td><td>0.75</td><td>0.05</td></tr><tr><td>LH*</td><td>0.9</td><td>0.02</td></tr><tr><td>LH*</td><td>1</td><td>0</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Performances of the snow model with the modified linear hysteresis for different weightings of the optimization criteria. The KGE' performances presented in the boxplots have been calculated on validation periods.</figcaption></figure>

We can observe that the runoff criterion remains stable in validation as the weighting performed in calibration decreased. The SCA criterion seems, on the other hand, to improve as the weighting on runoff decreases. Nevertheless, it seems that for weighting greater than 0.25 of SCA (i.e., less than 0.75 of runoff) the gains are very low on the SCA. Calibration on runoff only (1Q) significantly deteriorates the model's simulation of the SCA. Thus, a weighting of the optimization criterion with 75% allocated to the runoff criterion and 25% allocated to the SCA criterion seems to be a satisfactory compromise between the model's performance in terms of runoff and SCA. Figure 5 also shows that runoff performance is slightly better with a 0.75 weighting rather than a calibration only on the runoff criterion. This shows that the model potentially has higher robustness using the calibrated SCA. The robustness of the model will be more fully studied in section "APPLICATION OF THE SNOW MODEL TO DIFFERENT SIMULATION CONTEXTS".

## Analysis of the type of snow cover weighting to use

In section "Analysis of runoff criterion and snow cover criterion weighting", the SCA weightings were considered to be equivalent for all layers. Based on a 0.75 runoff weighting, we studied several ways to weight the SCA simulations for each elevation zone. The first reference combination used in the preceding sections considers an identical weighting for each elevation zone (here 5% for each elevation zone). Since the highest elevation zones have a greater snow regime than the lower layers, it is possible that the low-elevation zones provide little information to calibrate the model. Therefore, a calibration mode with a weighting that increases with elevation was tested, as was a second mode with a weighting that took only the high elevation zones into account. Finally, the last type of calibration, called global calibration, was directly calculated based on the mean SCA for the entire catchment, with no distinction made on the performance for each elevation zone. This global calibration is different to the equivalent weighting for all layers since there is just one hysteresis computed for all elevation bands, this should affect the calibration of the hysteresis parameters. The different weightings in this analysis are summarized in Table 3 and the results are presented in Figure 6.
<table><caption>**Table 3.** List of the SCA weighting methods used for the snow model calibration.</caption><thead><tr><th>Model</th><th>Calibration type</th><th>α</th><th>β1</th><th>β2</th><th>β3</th><th>β4</th><th>β5</th></tr></thead><tbody><tr><td>LH*</td><td>Equally Weighted (SCA E)</td><td>0.75</td><td>0.05</td><td>0.05</td><td>0.05</td><td>0.05</td><td>0.05</td></tr><tr><td>LH*</td><td>Global (SCA G)</td><td>0.75</td><td colspan="5">0.25, the SCA from the five-layer bands are averaged and compared to the mean observed SCA</td></tr><tr><td>LH*</td><td>Increasing (SCA I)</td><td>0.75</td><td>0.025 (10%)</td><td>0.025 (10%)</td><td>0.05 (20%)</td><td>0.075 (30%)</td><td>0.075 (30%)</td></tr><tr><td>LH*</td><td>High elevation (SCA H)</td><td>0.75</td><td>0</td><td>0</td><td>0</td><td>0.125</td><td>0.125</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Performances of the snow model with the modified linear hysteresis for different types of SCA criteria weighting. The KGE' performances presented in the boxplots have been calculated on validation periods.</figcaption></figure>

For the runoff validation, few significant differences were observed for the different methods taking the SCA selected into account. However, in terms of SCA validation, the results clearly show that use of the highest elevation zones only to calibrate the model is not sufficiently informative and gives unsatisfactory results compared to the other calibration methods. The three other types of calibration provide similar results in SCA validation. Identically weighting each elevation zone used in the preceding sections therefore seems to be a viable choice to calibrate the snow model with hysteresis and will be retained for further analysis.

**Overall results: hydrographs and analysis of the model's parameters**

The analyses above demonstrate that use of modified linear hysteresis (Eqs. 5, 6, 7, 8 and 9) with an optimization criterion with a 0.75 weighting on runoff and a 0.05 weighting on each of the elevation zones gives one of the best compromises in runoff and SCA performance over most of the catchments studied. The difference in performance between the reference CemaNeige routine calibrated on runoff and the model developed in this study, for each catchment, is presented in Figure 7.

The map in Figure 7 indicates that the performance in the runoff simulations of the reference CemaNeige model with hysteresis is similar for most catchments. The majority of the catchments gaining in performance seem to be the small catchments in the Massif Central. However, this figure does not show spatial coherence in the deterioration of the runoff performance criterion.

Concluding this analysis on various catchments, hydrographs as well as the SCA and the SWE over time are presented in Figure 8. The three catchments presented are the Durance at Embrun (X0310010, 2283 km², the Alps) and the Arac at Southern (O0384010, 170 km², the Pyrenees). The Durance alpine catchment has its performance improved by the new model (+0.06 points in KGE'), while the performance of the simulations of the Arac catchment remained stable (identical KGE').

These hydrographs show relatively few differences between the reference CemaNeige simulations and the simulation with CemaNeige improved by hysteresis. It can nonetheless be observed that the winter flow peaks were better reproduced on the Durance (X0310010). A clear improvement in the SCA simulation by CemaNeige with hysteresis, compared to the reference CemaNeige, can be observed for the two catchments. Figure 8 also shows that SWE was modified over time by the new SCA formulation in the CemaNeige model. Overall, it seems that the improvement tends to reduce the SWE accumulated. This shows that this new SCA parametrization also influences melting in the model. The impact of hysteresis on the model's parameters can be observed in Figure 9.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Improvement (or degradation) of KGE' performances on river discharge for each studied watershed. The value indicated is the difference between the performances of the modified CemaNeige with the hysteresis formulation minus the performance of the original CemaNeige model for the 2005–2010 validation period. 16 watersheds have significant degraded runoff performance (red), 104 have slightly degraded runoff (orange) and 51 watersheds have similar performances (yellow), while 75 have slightly improved performances (green) and 31 have seen significant improvement (dark green).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Hydrograph, evolution of SCA and SWE for three snow-dominated watersheds. These two watersheds were calibrated on the 2000–2005 period and the simulations presented in this figure are coming from the 2005–2010 validation period.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Analysis of the melt coefficient (Kf, in mm.°C⁻¹ d⁻¹), the cold content parameter (cT) and the hysteresis parameters (T<sub>h,acc</sub> and R<sub>sp</sub>) for the reference model (C Q 1) and the snow model with the modified linear hysteresis calibrated with a weighting of the runoff and SCA criteria (LH* 0.75), with a calibration on the runoff criterion only (LH* 1) and with a calibration on the SCA criterion only (LH* 0).</figcaption></figure>

Figure 9 shows that the melting parameter changes considerably between the reference model and the model with hysteresis calibrated by weighting the two objectives. Using the MODIS data and hysteresis constrains the melting parameter to much lower values than for calibration with runoff only. With quartiles between 2 and 4.5 mm °C⁻¹ d⁻¹, these values seem more plausible than what calibration of the reference model proposes (quartiles 4.5 and 38 mm °C⁻¹ d⁻¹), as values for degree day coefficient varying between 0.3 and 7.6 mm °C⁻¹ d⁻¹ can be found in the literature (Martinec and Rango, 1986). As snow melt rate decreases with SCA in the model (due to Eq. 3), the implemented hysteresis directly impacts the melt rate and the

melting parameter. Globally, the model with hysteresis has much higher SCA values compared to the original formulation (Fig. 8). This could explain the decreased value of the degree-day melt coefficient as higher SCA values tend to increase this melt rate (Eq. 3).

The cold content parameter cT seems to be less restricted by the model calibrated with MODIS than with the reference. The reference model's cT values cluster around a very low value (median, 0.04), indicating that for a consequential number of catchments, thermal inertia is considered to be very low. Figure 9 also shows that the parameters with T<sub>h,acc</sub> and R<sub>sp</sub> hysteresis require MODIS data to be constrained. The values are indeed
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Performances for runoff and SCA simulations, calculated on the whole validation period and calculated on the winter month of the validation period only.</figcaption></figure>

highly dispersed with calibration based only on runoff. To limit the addition of a parameter to the CemaNeige SCA, it seems cautious to set the accumulation threshold ($T_{h,acc}$). Its value varies little from one catchment to another (quartiles at 5.8 and 15 mm) and since the increase of SCA is very fast in accumulation periods, its influence on the model (whether for the SCA or runoff) should be less important than the melting threshold ratio ($R_{sp}$).

Finally, to check that the improvement of the winter runoff performance does not occur at the expanse of runoff performance for the rest of the year, Figure 10 allows one to verify the performance over the validation periods for criteria calculated only on winter runoff and criteria calculated on the whole year runoff.

The performance over the whole period is close to the winter period performance, considering that overall, taking the whole period into account tends to improve both the runoff and SCA scores. This remains true for the snow module with hysteresis and the reference CemaNeige SCA, meaning that the introduction of the hysteresis has a negligible influence on non-snow dependent runoff modelled by GR4J.

## APPLICATION OF THE SNOW MODEL TO DIFFERENT SIMULATION CONTEXTS

### Using the calibrated model on a period from the past

Following calibration of the snow model with the MODIS observation data, several questions come to mind. For example, is the model with hysteresis more robust in terms of climatic variability than the reference CemaNeige model? For this analysis, both models were run for three 10-year periods from the past using the parameters calibrated on the 2000–2005 and 2005–2010 periods. Only the catchments with runoff measurements for these three periods were used (representing approximately 50% of the sample used in this article). The three periods used were 1960–1970, 1970–1980, and 1980–1990. The boxplots in Figure 11 present the performances of the simulations using the model with hysteresis and the reference model for each of the past periods. The SCA performance criteria could not be calculated as the MODIS data were not available for these periods.

These results are encouraging because they show that the model calibrated on the SCA tends to keep slightly more stable performance for simulations far from the calibration period, as observed in the boxplots, which are slightly decentered toward a positive difference in performance. These preliminary results on the use of a set of parameters on a distant period should make the model with hysteresis more adapted than the reference model for analyses such as impact studies on climate change.

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Performances of the snow model with the modified linear hysteresis and the original CemaNeige model for three past periods: 1960–1970 (in red), 1970–1980 (in green), 1980–1990 (in blue). The parameters used here come from the calibrations of the 2000–2005 and 2005–2010 periods.</figcaption></figure>

### Sequential calibration of the snow model and the hydrological model

Here we wish to answer the following question: is it possible to calibrate the snow model completely independently of the rainfall-runoff model? This question seeks to determine which part of the information comes from MODIS SCA and whether it is sufficient to calibrate the snow model. Independently calibrating CemaNeige would facilitate the elaboration of a regionalization strategy for these parameters. Therefore, the snow model was first calibrated only on the SCA and then the rainfall-runoff model was calibrated on runoff observations, with the parameters from the snow model set to the values obtained in the first calibration (sequential calibration). This result was compared to the calibration of the rainfall-runoff model coupled with the snow model (reference or with hysteresis) whose parameters are set to the median value of the 277 catchments for each period (LH* med and C Q med) as well as a joint calibration of the snow model with hysteresis and the rainfall-runoff model (as was done above, LH* 0.75 here LH* ref).

The sequential calibration strategy applied is summarized in Table 4; the results are available in Figure 12.

<table><caption>**Table 4.** List of the calibrations realized for the analysis of the independency of the snow model from the rainfall-runoff model.</caption><thead><tr><th colspan="2">Snow model</th><th>Rainfall-runoff model</th></tr></thead><tbody><tr><td>LH* ref</td><td colspan="2">Coupled model optimization on Q and SCA, with α = 0.75 and SCA E</td></tr><tr><td>C Q ref</td><td colspan="2">Coupled model optimization on Q</td></tr><tr><td>Sequential calibration (LH* seq)</td><td>Optimization on SCA only</td><td>Optimization on Q only</td></tr><tr><td>LH* med</td><td>Fixed parameters at median value</td><td>Optimization on Q only</td></tr><tr><td>C Q med</td><td>Fixed parameters at median value</td><td>Optimization on Q only</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Performances of the model simulations used for sequential calibration analysis and the two references calibrated on runoff data presented in Table 4.</figcaption></figure>

The performances between the model calibrated sequentially and the original CemaNeige reference model are very similar for runoff (median performance of 0.81). Using median param-
eters for the snow models gives results that are also very similar to the sequential and reference calibrations. It therefore seems that the hydrological model can compensate for a large part of the modifications resulting from change in parameters from the snow model to the sequential calibration. Sequential calibration therefore contributes little to the rainfall-runoff model compared to the use of median parameters for the snow model.

## Application of the model to a ungauged context

Many mountain catchments are not gauged, and runoff measurements are therefore not available to calibrate the snow model or the hydrological model. On the other hand, the MODIS SCA data are available everywhere, meaning that, as above, they provide a calibration of the snow model for the SCA even if the catchment is not instrumented. The difference compared to the preceding section is that the rainfall-runoff model cannot be calibrated; the parameters are then set to their median value (the median of the parameters between all the catchments for each calibration period). This calibration of the snow model on SCA only is compared to the performance of the simulations with the reference CemaNeige model with hysteresis, with their parameters set at a median value. A summary of this analysis is available in Table 5 and the results are presented in Figure 13. The performance data presented are the KGE' deltas, the difference in performance of the catchments considered to be ungauged, and the same catchments calibrated on runoff with the reference model.

<table><caption>**Table 5.** List of the simulations made for analyzing the performance of the snow model for ungauged watersheds.</caption><thead><tr><th></th><th>Snow model</th><th>Runoff model</th></tr></thead><tbody><tr><td>LH* ref</td><td colspan="2">Coupled model optimization on Q and SCA with α = 0.75 and SCA E</td></tr><tr><td>C Q ref</td><td colspan="2">Coupled model optimization on Q</td></tr><tr><td>Sequential calibration (LH* SCA)</td><td>Optimization on SCA only</td><td>Fixed parameters at median value</td></tr><tr><td>LH* med</td><td>Fixed parameters at median value</td><td>Fixed parameters at median value</td></tr><tr><td>C Q med</td><td>Fixed parameters at median value</td><td>Fixed parameters at median value</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Performances of the models presented in Table 5 used in an ungauged watershed situation and the two reference simulations.</figcaption></figure>

Using the median parameters for the GR4J rainfall-runoff model implies a much greater loss in performance compared to the results shown in Figure 12 in section “Sequential calibration of the snow model and the hydrological model”. Setting parameters of the hydrological model at a median value is relatively harsh. However, we observe that calibrating the snow model on SCA with hysteresis gives slightly better results than using the median parameters of the reference CemaNeige model but the LH* model with median parameters still performs a little better

for flow simulation. It remains difficult to calibrate the snow model as well as the hydrological model without flow data.

## CONCLUSION

Many studies have shown the utility of calibrating the snow model with MODIS snow cover data (SCA). The objectives of the present study were to improve the SCA formulation of the CemaNeige snow-accounting model, to test several model calibration strategies using SCA MODIS data, and to assess their performance in a variety of application conditions.

The results of this study demonstrated that the original CemaNeige model gives highly satisfactory runoff simulations when it is combined with a rainfall-runoff model, but also that its internal snow cover surface variable cannot capture the MODIS SCA data. Calibrating CemaNeige in its original version with an optimization criterion taking runoff and SCA into account shows a clear reduction in performance in terms of runoff simulation, with only a slight improvement in performance with SCA.

Introducing a hysteresis in CemaNeige made it possible to go beyond the restrictions of the original model by improving the performance of the SCA simulation, without deteriorating the performance for runoff. Different SWE-SCA hystereses were tested: the modified simple linear formulation with a melting threshold dependent on annual solid precipitation presented the best results and was retained. Several analyses also determined that the most advantageous compromise for the model's evaluation criterion was a 75% weighting of the flow criterion, and an identical weighting of each of the elevation zones with 5% weighting on the SCA criterion.

These different analyses improved the CemaNeige model and allowed us to design an adapted calibration method taking the MODIS snow cover data into account. Introducing this hysteresis adds two parameters to the CemaNeige snow-accounting model, which brings the total number of parameters to four. Even though this improvement of the model requires additional parameters, it is constrained by both runoff and SCA instead of only runoff. Further development of the model could also use SCA data for correcting solid precipitation inputs as done by Shrestha et al. (2014), as snowfall measurements are generally underestimated due to precipitation gages under-catchment.

This CemaNeige snow-accounting model improved by hysteresis was compared to the original CemaNeige model for simulating time periods far in the past. The model with hysteresis gives slightly better performances than the original model and seems more robust. This is encouraging in the study of climate change impacts where the model's uncertainty is great in the context of climate non-stationarity. The snow model was also submitted to different calibration strategies, namely comparing a sequential calibration of its two components (ie. by calibrating the snow model only with SCA and the rainfall-runoff model only with runoff) or a calibration in ungauged conditions (ie. by calibrating the snow model only with SCA and using uncalibrated parameters for the rainfall-runoff model) to the combined calibration of the snow model and the rainfall-runoff model used in previous sections. The independent calibration of the snow module, i.e., performed only on the SCA, within a sequential calibration or the estimation of an ungauged catchment, gives fairly mixed results. These analyses were carried out using the median values of the parameters. This was the simplest parameter regionalization method that allowed testing the snow module in the most unfavorable condition possible. For the hydrological model, more advanced parameter
regionalization techniques could be tested to improve the model's performance on ungauged catchments.

The MODIS data used underwent a very simple treatment: no spatiotemporal data filtering system was applied. This could improve the model's performance by improving the calibration data (Parajka and Blöschl, 2008b). One of the future perspectives could also be the use of single SWE measurements and integrating them into the model calibration. It is highly probable that the conceptual model should be modified to provide a response that is comparable to the observation. The analysis of the sensitivity of the hysteresis parameters would be useful given that it would inform even more on their influence on runoff simulation. It is likely that the snow accumulation threshold could be set easily: this sensitivity analysis would make it possible demonstrate this and find the best adapted value.

*Acknowledgements.* The authors wish to thank Franck Masclet, intern, who conducted the preliminary analyses on incorporating hysteresis into the CemaNeige model. They would also like to thank the different entities that provided data, without which this analysis would not have been possible: Météo-France for the temperature and precipitation data, SAFRAN; the Banque Hydro for the runoff observation data; and NSIDC for making the MODIS SCA data available. Finally, we would like to thank the two reviewers, Jan Magnusson and Ryan Teuling, for their comments on an early version of the manuscript, which helped improving its quality.

## REFERENCES

Andreadis, K.M., Lettenmaier, D.P., 2006. Assimilating remotely sensed snow observations into a macroscale hydrology model. Adv. Water Resour., 29, 872–886. DOI: 10.1016/j.advwatres.2005.08.004.

Barnett, T.P., Adam, J.C., Lettenmaier, D.P., 2005. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 438, 303–309. DOI: 10.1038/nature04141.

Beniston, M., Farinotti, D., Stoffel, M., Andreassen, L.M., Coppola, E., Eckert, N., Fantini, A., Giacona, F., Hauck, C., Huss, M., Huwald, H., Lehning, M., López-Moreno, J.-I., Magnusson, J., Marty, C., Moran-Tejéda, E., Morin, S., Naaim, M., Provenzale, A., Rabatel, A., Six, D., Stötter, J., Strasser, U., Terzago, S., Vincent, C., 2017. The European mountain cryosphere: A review of past, current and future issues. Cryosphere Discuss, 2017, 1–60. DOI: 10.5194/tc-2016-290.

Bernsteinová, J., Bässler, C., Zimmermann, L., Langhammer, J., Beudert, B., 2015. Changes in runoff in two neighbouring catchments in the Bohemian Forest related to climate and land cover changes. J. Hydrol. Hydromech., 63, 342–352. DOI: 10.1515/joh.h-2015-0037.

Clark, M.P., Hendrikx, J., Slater, A.G., Kavetski, D., Anderson, B., Cullen, N.J., Kerr, T., Hreinsson, E. Ö., Woods, R.A., 2011. Representing spatial variability of snow water equivalent in hydrologic and land-surface models: A review. Water Resour. Res., 47, 7. DOI:10.1029/2011WR010745.

Coron, L., Perrin, C., Michel, C., Andréassian, V., Brigode, P., Delaigue, O., Le Moine, N., Mathevet, T., Mouelhi, S., Oudin, L., Pushpalatha, R., Thirel, G., Valéry, A., 2017a. airGR: Suite of GR Hydrological Models for Precipitation-Runoff Modeling. R package version 1.0.5.12. IRSTEA, Antony, France. https://irsteadoc.irstea.fr/cemoa/PUB00052697

Coron, L., Thirel, G., Delaigue, O., Perrin, C., Andréassian, V., 2017b. The suite of lumped GR hydrological models in an R package. Environ. Model. Softw., 94, 166–171. DOI: 10.1016/j.envsoft.2017.05.002.

Da Ronco, P., De Michele, C., 2014. Cloud obstruction and snow cover in Alpine areas from MODIS products. Hydrol. Earth Syst. Sci., 18, 4579–4600. DOI: 10.5194/hess-18-4579-2014.

Duethmann, D., Peters, J., Blume, T., Vorogushyn, S., Güntner, A., 2014. The value of satellite-derived snow cover images for calibrating a hydrological model in snow-dominated catchments in Central Asia. Water Resour. Res., 50, 2002–2021. DOI: 10.1002/2013WR014382.

Edijatno, Nascimento, N.D.O., Yang, X., Makhlouf, Z., Michel, C., 1999. GR3J: a daily watershed model with three free parameters. Hydrol. Sci. J., 44, 263–277. DOI: 10.1080/02626669909492221.

Egli, L., Jonas, T., 2009. Hysteretic dynamics of seasonal snow depth distribution in the Swiss Alps. Geophys. Res. Lett., 36, L02501. DOI: 10.1029/2008GL035545.

Essery, R., Pomeroy, J., 2004. Implications of spatial distributions of snow mass and melt rate for snow-cover depletion: theoretical considerations. Ann. Glaciol., 38, 261–265. DOI: 10.3189/172756404781815275.

Franz, K.J., Karsten, L.R., 2013. Calibration of a distributed snow model using MODIS snow covered area data. J. Hydrol., 494, 160–175. DOI: 10.1016/j.jhydrol.2013.04.026.

Friedman, M., 1937. The use of ranks to avoid the assumption of normality implicit in the analysis of variance. J. Am. Stat. Assoc., 32, 675–701. DOI: 10.1080/01621459.1937.10503522.

Gafurov, A., Bárdossy, A., 2009. Cloud removal methodology from MODIS snow cover product. Hydrol. Earth Syst. Sci., 13, 1361–1373. DOI: 10.5194/hess-13-1361-2009.

Grayson, R.B., Blöschl, G., Western, A.W., McMahon, T.A., 2002. Advances in the use of observed spatial patterns of catchment hydrological response. Adv. Water Resour., 25, 1313–1334. DOI: 10.1016/S0309-1708(02)00060-X.

Gupta, H.V., Kling, H., Yilmaz, K.K., Martinez, G.F., 2009. Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. J. Hydrol., 377, 80–91. DOI: 10.1016/j.jhydrol.2009.08.003.

Hall, D.K., Salomonson, V.V., Riggs, G.A., 2006. MODIS/Terra, MODIS/Aqua Snow Cover Daily L3 Global 500m Grid, Version 5. NASA National Snow and Ice Data Center, Boulder, Colorado, USA. DOI: 10.5067/63NQASRDPDB0.

He, Z.H., Parajka, J., Tian, F.Q., Blöschl, G., 2014. Estimating degree-day factors from MODIS for snowmelt runoff modeling. Hydrol. Earth Syst. Sci., 18, 4773–4789. DOI: 10.5194/hess-18-4773-2014.

Helbig, N., van Herwijnen, A., Magnusson, J., Jonas, T., 2015. Fractional snow-covered area parameterization over complex topography. Hydrol. Earth Syst. Sci., 19, 1339–1351. DOI: 10.5194/hess-19-1339-2015.

Hublart, P., Ruelland, D., García de Cortázar-Atauri, I., Gascoin, S., Lhermitte, S., Ibacache, A., 2016. Reliability of lumped hydrological modeling in a semi-arid mountainous catchment facing water-use changes. Hydrol. Earth Syst. Sci., 20, 3691–3717. DOI: 10.5194/hess-20-3691-2016.

Kling, H., Fuchs, M., Paulin, M., 2012. Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios. J. Hydrol., 424–425, 264–277. DOI: 10.1016/j.jhydrol.2012.01.011.

Kolberg, S.A., Gottschalk, L., 2006. Updating of snow depletion curve with remote sensing data. Hydrol. Process., 20, 2363–2380. DOI: 10.1002/hyp.6060.

Krajčí, P., Holko, L., Perdigão, R.A.P.P., Parajka, J., 2014. Estimation of regional snowline elevation (RSLE) from MODIS images for seasonally snow covered mountain basins. J. Hydrol., 519, 1769–1778. DOI: 10.1016/j.jhydrol.2014.08.064.

Krajčí, P., Holko, L., Parajka, J., 2016. Variability of snow line elevation, snow cover area and depletion in the main Slovak basins in winters 2001–2014. J. Hydrol. Hydromech., 64, 12–22. DOI: 10.1515/johh-2016-0011.
Le Moine, N., Andréassian, V., Perrin, C., Michel, C., 2007. How can rainfall-runoff models handle intercatchment groundwater flows? Theoretical study based on 1040 French catchments. Water Resour. Res., 43, W06428. DOI: 10.1029/2006WR005608.

Liston, G.E., 2004. Representing subgrid snow cover heterogeneities in regional and global models. J. Clim., 17, 1381–1397. DOI: 10.1175/1520-0442(2004)017<1381:RSSCHI>2.0.CO;2.

Luce, C.H., Tarboton, D.G., 2004. The application of depletion curves for parameterization of subgrid variability of snow. Hydrol. Process., 18, 1409–1422. DOI: 10.1002/hyp.1420.

Magand, C., Ducharne, A., Le Moine, N., Gascoin, S., 2014. Introducing hysteresis in snow depletion curves to improve the water budget of a land surface model in an Alpine catchment. J. Hydrometeorol., 15, 631–649. DOI: 10.1175/JHM-D-13-091.1.

Magnusson, J., Gustafsson, D., Hüsler, F., Jonas, T., 2014. Assimilation of point SWE data into a distributed snow cover model comparing two contrasting methods. Water Resour. Res., 50, 7816–7835. DOI: 10.1002/2014WR015302.

Martinec, J., Rango, A., 1986. Parameter values for snowmelt runoff modelling. J. Hydrol., 84, 197–219. DOI: 10.1016/0022-1694(86)09123-X.

Nitta, T., Yoshimura, K., Takata, K., O'ishi, R., Sueyoshi, T., Kanae, S., Oki, T., Abe-Ouchi, A., Liston, G.E., 2014. Representing variability in subgrid snow cover and snow depth in a global land model: offline validation. J. Clim., 27, 3318–3330. DOI: 10.1175/JCLI-D-13-00310.1.

Niu, G.-Y., Yang, Z.-L., 2007. An observation-based formulation of snow cover fraction and its evaluation over large North American river basins. J. Geophys. Res. Atmospheres, 112, D21101. DOI: 10.1029/2007JD008674.

Parajka, J., Blöschl, G., 2008a. The value of MODIS snow cover data in validating and calibrating conceptual hydrologic models. J. Hydrol., 358, 240–258. DOI: 10.1016/j.jhydrol.2008.06.006.

Parajka, J., Blöschl, G., 2008b. Spatio-temporal combination of MODIS images - potential for snow cover mapping. Water Resour. Res., 44, 3. DOI: 10.1029/2007WR006204.

Parajka, J., Pepe, M., Rampini, A., Rossi, S., Blöschl, G., 2010. A regional snow-line method for estimating snow cover from MODIS during cloud cover. J. Hydrol., 381, 203–212. DOI: 10.1016/j.jhydrol.2009.11.042.

Parajka, J., Haas, P., Kirnbauer, R., Jansa, J., Blöschl, G., 2012. Potential of time-lapse photography of snow for hydrological purposes at the small catchment scale. Hydrol. Process., 26, 3327–3337. DOI: 10.1002/hyp.8389.

Perrin, C., Michel, C., Andréassian, V., 2001. Does a large number of parameters enhance model performance? Comparative assessment of common catchment model structures on 429 catchments. J. Hydrol., 242, 275–301. DOI: 10.1016/S0022-1694(00)00393-0.

Perrin, C., Michel, C., Andréassian, V., 2003. Improvement of a parsimonious model for streamflow simulation. J. Hydrol., 279, 275–289. DOI: 10.1016/S0022-1694(03)00225-7.

Poggio, L., Gimona, A., Brown, I., 2012. Spatio-temporal MODIS EVI gap filling under cloud cover: An example in Scotland. ISPRS J. Photogramm. Remote Sens., 72, 56–72. DOI: 10.1016/j.isprsjprs.2012.06.003.

Pokhrel, B.K., Chevallier, P., Andréassian, V., Tahir, A.A., Arnaud, Y., Neppel, L., Bajracharya, O.R., Budhathoki, K.P., 2014. Comparison of two snowmelt modelling approaches in the Dudh Koshi basin (eastern Himalayas, Nepal). Hydrol. Sci. J., 59, 1507–1518. DOI: 10.1080/02626667.2013.842282.

Pushpalatha, R., Perrin, C., Le Moine, N., Mathevet, T., Andréassian, V., 2011. A downward structural sensitivity analysis of hydrological models to improve low-flow simulation. J. Hydrol., 411, 66–76. DOI: 10.1016/j.jhydrol.2011.09.034.

Quintana-Seguí, P., Le Moigne, P., Durand, Y., Martin, E., Habets, F., Baillon, M., Canellas, C., Franchisteguy, L., Morel, S., 2008. Analysis of near-surface atmospheric variables: validation of the SAFRAN analysis over France. J. Appl. Meteorol. Climatol., 47, 92–107. DOI: 10.1175/2007JAMC1636.1.

Rodell, M., Houser, P.R., 2004. Updating a land surface model with MODIS-derived snow cover. J. Hydrometeorol., 5, 1064–1075. DOI: 10.1175/JHM-395.1.

Shrestha, M., Wang, L., Koike, T., Tsutsui, H., Xue, Y., Hirabayashi, Y., 2014. Correcting basin-scale snowfall in a mountainous basin using a distributed snowmelt model and remote-sensing data. Hydrol. Earth Syst. Sci., 18, 747–761. DOI: 10.5194/hess-18-747-2014.

Slater, A.G., Clark, M.P., 2006. Snow data assimilation via an Ensemble Kalman Filter. J. Hydrometeorol., 7, 478–493. DOI: 10.1175/JHM505.1.

Swenson, S.C., Lawrence, D.M., 2012. A new fractional snow-covered area parameterization for the Community Land Model and its effect on the surface energy balance. J. Geophys. Res. Atmospheres, 117. DOI:10.1029/2012JD018178.

Thirel, G., Salamon, P., Burek, P., Kalas, M., 2013. Assimilation of MODIS snow cover area data in a distributed hydrological model using the particle filter. Remote Sens., 5, 5825–5850. DOI: 10.3390/rs5115825.

Thirel, G., Andréassian, V., Perrin, C., 2015a. On the need to test hydrological models under changing conditions. Hydrol. Sci. J., 60, 1165–1173. DOI: 10.1080/02626667.2015.1050027.

Thirel, G., Andréassian, V., Perrin, C., Audouy, J.-N., Berthet, L., Edwards, P., Folton, N., Furusho, C., Kuentz, A., Lerat, J., Lindström, G., Martin, E., Mathevet, T., Merz, R., Parajka, J., Ruelland, D., Vaze, J., 2015b. Hydrology under change: an evaluation protocol to investigate how hydrological models deal with changing catchments. Hydrol. Sci. J., 60, 1184–1199. DOI: 10.1080/02626667.2014.967248.

Troin, M., Arsenault, R., Brissette, F., 2015. Performance and uncertainty evaluation of snow models on snowmelt flow simulations over a nordic catchment (Mistassibi, Canada). Hydrology, 2, 289–317. DOI: 10.3390/hydrology2040289.

Troin, M., Poulin, A., Baraer, M., Brissette, F., 2016. Comparing snow models under current and future climates: Uncertainties and implications for hydrological impact studies. J. Hydrol., 540, 588–602. DOI: 10.1016/j.jhydrol.2016.06.055.

Valéry, A., Andréassian, V., Perrin, C., 2014a. "As simple as possible but not simpler": What is useful in a temperature-based snow-accounting routine? Part 1 – Comparison of six snow accounting routines on 380 catchments. J. Hydrol. 517, 1166–1175. DOI: 10.1016/j.jhydrol.2014.04.059.

Valéry, A., Andréassian, V., Perrin, C., 2014b. "As simple as possible but not simpler": What is useful in a temperature-based snow-accounting routine? Part 2 – Sensitivity analysis of the Cemaneige snow accounting routine on 380 catchments. J. Hydrol., 517, 1176–1187. DOI: 10.1016/j.jhydrol.2014.04.058.

Vidal, J.-P., Martin, E., Franchistéguy, L., Baillon, M., Soubeyroux, J.-M., 2010. A 50-year high-resolution atmospheric reanalysis over France with the Safran system. Int. J. Climatol., 30, 1627–1644. DOI: 10.1002/joc.2003.

Vuyovich, C.M., Jacobs, J.M., Daly, S.F., 2014. Comparison of passive microwave and modeled estimates of total watershed SWE in the continental United States. Water Resour. Res., 50, 11, 9088–9102. DOI: 10.1002/2013WR014734.

Zaitchik, B.F., Rodell, M., 2009. Forward-looking assimilation of MODIS-derived snow-covered area into a land surface model. J. Hydrometeorol., 10, 130–148. DOI: 10.1175/2008JHM1042.1.