Performance and Uncertainty Evaluation of Snow Models on Snowmelt Flow Simulations over a Nordic Catchment (Mistassibi, Canada)

Abstract

An analysis of hydrological response to a multi-model approach based on an ensemble of seven snow models (SM; degree-day and mixed degree-day/energy balance models) coupled with three hydrological models (HM) is presented for a snowmelt-dominated basin in Canada. The present study aims to compare the performance and the reliability of different types of SM-HM combinations at simulating snowmelt flows over the 1961–2000 historical period. The multi-model approach also allows evaluating the uncertainties associated with the structure of the SM-HM ensemble to better predict river flows in Nordic environments. The 20-year calibration shows a satisfactory performance of the ensemble of 21 SM-HM combinations at simulating daily discharges and snow water equivalents (SWEs), with low streamflow volume biases. The validation of the ensemble of 21 SM-HM combinations is conducted over a 20-year period. Performances are similar to the calibration in simulating the daily discharges and SWEs, again with low model biases for streamflow. The spring-snowmelt-generated peak flow is captured only in timing by the ensemble of 21 SM-HM combinations. The results of specific hydrologic indicators show that the uncertainty related to the choice of the given HM in the SM-HM combinations cannot be neglected in a more quantitative manner in simulating snowmelt flows. The selection of the SM plays a larger role than the choice of the SM approach (degree-day versus mixed degree-day/energy balance) in simulating spring flows. Overall, the snow models provide a low degree of uncertainty to the total uncertainty in hydrological modeling for snow hydrology studies.

1. Introduction

In recent decades, important developments have been made in hydrological modeling, including growth in computational power and improved understanding of the physics of the hydrological system. This has led to higher complexity of hydrological models (HM) ranging from lumped and conceptual models to distributed physics-based models. These advances have increased the need to deal with the uncertainty associated with the HM themselves. It is now broadly recognized that consideration of uncertainty in hydrological simulations is critical for both research and operational modeling [1]. In addition, the value of a hydrological simulation and associated decision-making processes are limited if estimates of uncertainty are not provided [2].

To adequately assess uncertainty in hydrological modeling, there are three aspects to be considered: understanding, quantification and reduction of uncertainty; understanding uncertainty being an essential part of both the quantification and reduction of uncertainty [3]. However, the understanding of hydrological uncertainty is still far from complete, and further efforts are required to improve this situation. As far as assessing hydrological modeling uncertainty is concerned, two main sources are worth considering: model parameters and structure—even though most studies have ruled on the greatest uncertainty being that of model structure [4,5]. Since the HM is structured in multiple components, each model component has its own source of uncertainty that contributes to total uncertainty of the HM structure. Consequently, the uncertainty associated with each HM component should be quantified in view of reducing the uncertainty associated with the HM structure. Among the various HM components, the uncertainties related to the most sensitive components to changing climatic conditions (i.e., precipitation and temperature) need to be investigated, since the reliability of streamflow simulations are strongly limited by the understanding of these uncertainties. Of the most sensitive components to changing climatic conditions, evaluating how the snow models (SM) affect the streamflow simulations might contribute to a better understanding of HM structure uncertainty. In addition, this information will help to narrow the total uncertainty of hydrological modeling in impact studies.

Many HMs contain a snow hydrology component in their basic form [6,7,8,9]. However, the physical processes involved in seasonal snowpack evolution are dealt with heterogeneously by HMs when applied over snow-covered catchments [10,11]: the SMs differ from each other based on their pathway at simulating snowmelt processes only. The amount of water held within the snowpack waiting to be melted (referred to as the snow water equivalent or SWE) can be estimated by mean of three approaches with various levels of complexity: physical energy-balance (EB) models, degree-day (DD) models, and mixed degree-day and energy-balance (DD/EB) models. The EB models explicitly simulate energy and mass exchanges between internal layers of the snowpack, as well as snowpack stratigraphy. Each of the relevant energy fluxes at the snow surface is computed from physically-based calculations using direct observations of the required meteorological variables; the SWE is then calculated as the sum of the individual fluxes [12]. Nevertheless, their applicability is often limited to studies of glacier hydrology and dynamics [13,14,15]. The DD models employ the air temperature as a proxy of the heat transfer process affecting snowmelt to compute the SWE [16]; however, additional input variables, such as incoming shortwave radiation or albedo, may be incorporated through empirical formulas based on time and location. The DD models are widely used for hydrologically-related applications due to their parsimony in data requirements compared to more sophisticated EB models [17,18,19,20,21]. The DD models’ application is therefore restricted to simulation of SWE at daily or coarser resolution, and in a lumped or semi-lumped manner for the calculation of the average SWE over a whole catchment. Recently, an increasing need for high temporal and spatial resolution simulation of the SWE for hydrological modeling purposes has prompted numerous attempts to combine the accuracy of physically-based EB models with the simplicity of DD models in order to develop mixed EB/DD models [22]. In this respect, in addition to assessing the advantages of using EB/DD models rather than DD models to obtain more robust SWE simulations over a particular Nordic catchment of interest [23], one pertinent issue remains for snow hydrological studies—being able to propose an accurate ensemble of SM-HM combinations in order to adequately depict the spring snowmelt-generated streamflow at the catchment scale. Resolving this issue will allow us to explore the uncertainty associated with the structure of snow hydrology modeling, a critical step towards the understanding, quantification and reduction of total uncertainty in hydrological simulations.

The aim of the present study is to evaluate and compare the performance and the reliability of different types of SM-HM combinations at simulating the snowmelt-dominated streamflow at the catchment scale by taking into account uncertainty from the SM-HM ensemble structure. The multi-model approach is built on the basis of 21 different SM-HM structures using seven snow models (four DD models and three DD/EB models) combined to three lumped conceptual hydrological models over the 1961–2000 period. The approach was implemented over the Mistassibi Basin situated on the Saguenay-Lac-St-Jean region of Quebec, representing typical Canadian climatic conditions and land cover types. Snowmelt processes dominate the surface hydrology of the basin.

This article is organized as follows. The experimental design, including the study area and the data, is described in Section 2. The multi-model approach is presented in Section 3. The relevant results from the multi-model approach are analyzed in Section 4. Uncertainties associated with the SMs and the HMs on the snowmelt flow simulations are discussed in Section 5. Concluding remarks appear in Section 6.

2. Experimental Design

2.1. Study Area

The Mistassibi River Basin is located in the Saguenay-Lac-St-Jean region of Quebec, Canada (Figure 1). The catchment covers an area of 9320 km² with a latitudinal extension from 52°N to 48°N and from 71°W to 73°W in the east-west direction. The watershed originates in an area about 80 kilometers east of Lake Mistassini and drains into the Lac-Saint-Jean. The elevations range from 102 to 758 m in the basin (Figure 1). Despite its large size, the watershed is extremely homogeneous in terms of land cover (boreal forest) and topography.

Figure 1. Location map of the Mistassibi Basin.

Figure 1. Location map of the Mistassibi Basin.

2.2. Data

Natural Resources Canada (NRCan) provided daily 10 km gridded precipitation, maximum and minimum temperatures data for the basin over the 1961–2000 period. These data have been interpolated from daily station observations using the trivariate thin-plate smoothing splines [24]. This method allows reducing spurious trends or artifacts introduced by changes in collection techniques, station relocations, or inclusions of stations with different record lengths. The interpolation and gridding technique maintains the consistency of spatial information and of time-series characteristics from the station observations. Daily precipitation, maximum and minimum temperatures for the Mistassibi Basin are shown in Figure 2.

Snow is measured by Rio Tinto Alcan (RTA) through four observing stations over the catchment (Black stars in Figure 1). RTA uses snow sampling tubes to measure snow depth and water content. The measurements are made at five points around the observing stations which are representative of the immediate area, and then averaged. The snow dataset is overwhelmingly dominated by monthly observations made on-or-near the 30^th of each month, from January to May, over the 1961–2000 period. A summary of meteorological and snow data statistics is given in Table 1.

The Centre d’Expertise Hydrique du Québec (CEHQ) provided observed daily streamflow at the Mistassibi station located at the outlet of the basin (black square in Figure 1). The hydrologic regime is highly variable with observed daily discharges ranging from 25 to 1565 m³/s with a mean of 285 m³/s over the 1961–2000 period (Table 1). Peak flows in spring are associated to snowmelt events (Figure 2). High flows observed in summer and autumn are related to precipitation events, whereas streamflow in winter becomes very low as a result of the accumulation of precipitation in the snowpack (Figure 2).

Figure 2. Mean monthly climatology at the Mistassibi Basin: (A) minimal temperature, (B) maximal temperature, (C) precipitation and (D) snow water equivalent (SWE); and (E) monthly discharge at the Mistassibi station over 1961–2000. The vertical bars represent the maximum and minimum values for each meteorological variable.

Figure 2. Mean monthly climatology at the Mistassibi Basin: (A) minimal temperature, (B) maximal temperature, (C) precipitation and (D) snow water equivalent (SWE); and (E) monthly discharge at the Mistassibi station over 1961–2000. The vertical bars represent the maximum and minimum values for each meteorological variable.

Table 1. General characteristics of the Mistassibi Basin.

	Mistassibi Basin
Province	Quebec
Drainage area (km2)	9320
Elevation range (m)	102–758
Main land cover type	Forest (84%)
Hydrology
Mean daily discharge (m3/s)	285
Minimum daily discharge (m3/s)	25
Maximum daily discharge (m3/s)	1565
Snow
Mean daily SWE (mm)	203
Maximum daily SWE (mm)	347
Climate
Annual average precipitation total (mm)	978
Annual daily temperature (°C)
Min.	−46.1
Max.	+35.1
Mean	−0.9

Table 1. General characteristics of the Mistassibi Basin.

	Mistassibi Basin
Province	Quebec
Drainage area (km2)	9320
Elevation range (m)	102–758
Main land cover type	Forest (84%)
Hydrology
Mean daily discharge (m3/s)	285
Minimum daily discharge (m3/s)	25
Maximum daily discharge (m3/s)	1565
Snow
Mean daily SWE (mm)	203
Maximum daily SWE (mm)	347
Climate
Annual average precipitation total (mm)	978
Annual daily temperature (°C)
Min.	−46.1
Max.	+35.1
Mean	−0.9

3 The Multi-Model Approach

The study’s general procedure for evaluating the performance and the uncertainties of snow models on snowmelt runoff simulations using different combinations of SM-HM is displayed in Figure 3. The multi-model approach involves the coupling of the selected range of seven snow models with three hydrological models to obtain an ensemble of 21 streamflow simulations. Each of the SMs and HMs used in this study is briefly described below.

Figure 3. Schematic description of the multi-model approach used in this study.

Figure 3. Schematic description of the multi-model approach used in this study.

3.1. Description of Snow Models

Four SMs extracted from the MOHYSE, HBV, HMETS, HYDROTEL hydrological models and three additional independent SMs (CEMANEIGE, SEB and ETI) are compared in the present study. The seven SMs can be classified into DD models (MOHYSE, HBV, CEMANEIGE and HMETS) and DD/EB models (HYDROTEL, SEB and ETI). Table 2 lists the conceptual differences between the seven SMs. In the following, only the algorithms related to the snowmelt routines in the SMs are described.

Table 2. Conceptual differences between the seven snow models. The acronyms of models used in this study are also specified.

	Model Type	Input Data a	Snow Parameters	State Variables
MOHYSE (MOH)	Degree-day	T, P	2	1
HBV (HBV)	Degree-day	T, P	3	2
CEMANEIGE (CEM)	Degree-day	T, P	2	2
HMETS (ETS)	Degree-day	T, P	10	5
HYDROTEL (HYD)	Degree-day/energy balance	T, P	5	4
SEB (SEB)	Degree-day/energy balance	T, P, ISR, α	2	1
ETI (ETI)	Degree-day/energy balance	T, P, ISR, α	2	1

^a T = temperature; P = precipitation amount (liquid and/or solid); ISR = incoming shortwave radiation and α = snow albedo.

Table 2. Conceptual differences between the seven snow models. The acronyms of models used in this study are also specified.

	Model Type	Input Data a	Snow Parameters	State Variables
MOHYSE (MOH)	Degree-day	T, P	2	1
HBV (HBV)	Degree-day	T, P	3	2
CEMANEIGE (CEM)	Degree-day	T, P	2	2
HMETS (ETS)	Degree-day	T, P	10	5
HYDROTEL (HYD)	Degree-day/energy balance	T, P	5	4
SEB (SEB)	Degree-day/energy balance	T, P, ISR, α	2	1
ETI (ETI)	Degree-day/energy balance	T, P, ISR, α	2	1

^a T = temperature; P = precipitation amount (liquid and/or solid); ISR = incoming shortwave radiation and α = snow albedo.

3.1.1. Degree-Day Models

MOHYSE

The SM derived from MOHYSE (MOdèle HYdrologique Simplifié à l’Extrême) is a simple model based on a temperature index [25]. When the mean daily air temperature (T_i) is less than a threshold temperature value (T_s), the precipitation is considered as snow and the liquid-water equivalent of the snow precipitation is added to the snowpack. The snowpack content increases with additional snowfall and decreases with snowmelt. The mass balance for the snowpack is computed as follows:

SWE_i = SWE_{i – 1} + SNOW_i + WE_SNOWMELT_i

(1)

with

WE_SNOWMELT_i = min (c_f · max (T_i – T_s, 0), SWE_{i – 1})

(2)

where SWE_i is the snow water equivalent within the snowpack on day i (mm of water); SNOW_i is the water equivalent of snow precipitation on day i (mm of water); WE_SNOWMELT_i is the water equivalent of snowmelt on day i (mm of water); and c_f is the melt factor (mm/°C/day).

HBV

The HBV SM (Hydrologiska Byråns Vattenbalansavdelning; [26] is based on a DD relation. The precipitation is assumed to accumulate as snow when the mean daily air temperature drops below a threshold temperature value (T_s). Snowmelt begins if the mean daily air temperature (T_i) is above the threshold value, calculated with an empirical DD equation as follows:

SNOWMELT_i = c_f · (T_i – T_s)

(3)

where SNOWMELT_i is the snowmelt on day i (mm of water).

The snowpack is assumed to retain meltwater until the amount exceeds a certain fraction of SWE. When the mean daily air temperature decreases below the threshold temperature value, liquid water within the snowpack refreezes based on the following algorithm:

REFSWE_i = c_fr · c_{f max} · (T_s – T_i)

(4)

where REFSWE_i is the amount of water that refreezes on day i (mm of water) and c_fr is the refreezing melt factor (mm/°C/day).

CEMANEIGE

The DD-based CEMANEIGE SM can be viewed as a complexified version of the MOHYSE SM to improve its snowpack dynamics simulation. Parsimony, performance and robustness were the main objectives of the CEMANEIGE SM development [27].

CEMANEIGE accounts for the daily evolution of the thermal state of the snowpack (eT_G) which is used to determine potential snowmelt (PSNOWMELT_i) on day i as:

$$e{T}_{Gi}=c_{TG} · e{T}_{Gi} + (1 - c_{TG}) · T_i$$

(5)

PSNOWMELT_i = c_f · (T_i – T_s) if eT_G = 0 and T_i > T_s

(6)

PSNOWMELT_i = SWE_i if PSNOWMELT_i > SWE_i

(7)

where c_TG is the coefficient of the thermal state of the snowpack.

The mass balance for the snowpack is then computed as follows:

SWE_i = SWE_{i − 1} + SNOW_i − ASNOWMELT_i

(8)

with

ASNOWMELT = (0.9.SNOWCOVER_i +0.1) · PSNOWMELT_i

(9)

where ASNOWMELT_i is the actual snowmelt (mm of water) and SNOWCOVER_i is the percentage of the catchment covered by snow on day i.

HMETS

The SM derived from HMETS (Hydrological Model of Ecole de Technologie Supérieure) is based on the work of Vehviläinen [28] that allows accounting for melting and refreezing processes within the snowpack. The SM includes three successive steps: overnight refreezing process, snowmelt, and snowpack water retention capacity.

Freezing of liquid water in the snowpack is effective if the mean daily diurnal temperature (T_dt) is below the freezing temperature threshold (T_bf). The potential amount of overnight refreezing (POR_i) on day i is estimated as:

POR_i = c_fr · (T_bf − T_dt)^fe

(10)

where fe is the empirical exponent. The actual amount of refreezing cannot exceed the amount of liquid water in the snowpack.

The potential snowmelt (PSNOWMELT_i) on day i is then calculated as:

PSNOWMELT_i = max (0, ddf_i (T_i – T_bm)) if T_i > T_bm

(11)

with

ddf_i = ddf_min · (1 + k_cum · CSM_i)

(12)

where ddf_i is a DD factor that can vary between a minimum (ddf_min) to a maximal value as a function of CSM on day i (mm/°C); T_bm is the base melting temperature (°C); k_cum is the empirical parameter (mm⁻¹ of water); and CSM_i is the cumulative amount of snowmelt on day i (mm of water). The snowpack aging is accounted for in the model through CSM. This variable is reset to zero if the snowpack disappears. The actual snowmelt cannot exceed the amount of snow on ground.

The snowpack water retention capacity (WRF_i) varies as a function of the snowpack aging represented by CSM_i on day i as follows:

WRF_i = max (fc_min, fc_max · (1 – c_cum · CSM_i))

(13)

where fc_min and fc_max are the minimum and maximal fractions for the snowpack water retention capacity, respectively, and c_cum is an empirical parameter (mm⁻¹ of water). The water in the snowpack is supplied by the actual snowmelt and liquid precipitation.

3.1.2. Mixed Degree-Day and Energy Balance Models

HYDROTEL

The SM derived from HYDROTEL is based on a DD/EB approach [29]. The model uses partially empirical relationships to simulate air-snow and ground-snow interface melt, compaction, albedo evolution and the liquid water retained by the snow cover. The SM defines the precipitation phase from air temperature data by using a fixed threshold set by the user. In case of rain-on-snow events, the water is added to SWE, and the heat deficit is adjusted. The model simulates convection losses if the mean daily air temperature is below the melt temperature threshold. The heat deficit is influenced by melt at the soil-snow interface and the air-snow interface. The latter is calculated using a modified DD equation that requires an empirical melt factor. The albedo of fresh snow is assumed to be 0.8 and, in the absence of new snowfall, decreases exponentially over time until it reaches a minimum value, typically set to 0.5. Snow compaction and maximum density are constant and calibrated by the user. The HYDROTEL SM is based on the evolution of two variables, SWE and heat deficit. The SWE is calculated using a mass balance as follows:

SWE_i = SWE_{i – 1} + P_liq,i + SNOW_i − M_i + WATER_i

(14)

where WATER_i is the water equivalent of snow on the ground (mm of water) and M_i is water from melting, calculated after obtaining the heat deficit using an energy balance model:

U_i = U_{i – 1} + U_n,i + U_p,i − U_c,i − U_a-s,i − U_g-s,i

(15)

where U_i is the calorific deficit on day i; U_n,i is the calorific deficit from solid precipitation on day i; U_p,i is the calorific deficit from liquid precipitation on day i; U_c,i is the calorific deficit from heat losses by convection on day i; U_a-s,i is the calorific deficit from snowmelt at the air-snow interface on day i; and U_g-s,i is the calorific deficit from snowmelt at the snow-ground interface on day i. All terms are expressed in J/m². Melting is estimated by dividing the heat surplus by the product of the melting heat (C_f), the water density (ρ_w):

$$M_i=\frac{-U_i}{C_f · {\rho }_w }$$

(16)

Additional details about the equations of the HYDROTEL SM can be found in Turcotte et al. [29].

SEB

The SEB (Simplified Energy Balance) model was developed to provide a mixed model between the physically based energy balance and the empirical degree-day models [30,31]. The mass and energy balances for the snowpack are provided by:

SWE_i = SWE_i + SNOW_i − SNOWMELT_i

(17)

with

$${SNOWMELT}_i=\frac{-Q_{S,i}}{\rho · L_f }·1000$$

(18)

and

Q_S,i = ISR_i · (1 − a_i) + C₀ + C₁ · T_i

(19)

where Q_S,i is the energy available for snowmelt on day i (W/m²); ρ is the snow density (kg/m³); L_f is the latent heat of fusion (J/kg); ISR_i the incoming shortwave radiation on day i (W/m²); and C₀ and C₁ are two empirical factors that account for the temperature-dependent energy fluxes (net longwave radiation and turbulent heat fluxes) (W/m² and W/m²/°C, respectively). No threshold temperature is assumed as discriminant to the snowmelt occurrence.

The incoming shortwave radiation data was not available in this study because no pyrheliometer or pyranometer has been installed in the basin or its surroundings due to its relative isolated location. Instead, ISR was estimated using the “Air–Sea Matlab Toolbox” from the United-States Geological Survey’s Woods Hole Science Center (http://woodshole.er.usgs.gov/operations/sea-mat/index.html). ISR values were then compared to reanalysis data taken from the ERA-Interim product for the study area to ensure they were representative of expected ISR values [32].

Snow albedo is calculated as a logarithmic function of accumulated daily maximum positive temperature since snowfall:

a_i = p₁ − p₂ log₁₀ Ta_i

(20)

where Ta_i is accumulated daily maximum temperatures above 0 °C since snowfall on day i (°C); and p₁ and p₂ are empirical coefficients, where p₁ is the albedo of fresh snow (for Ta = 1 °C).

ETI

The ETI (Enhanced Temperature-Index) model is described in detail by Pellicciotti et al. [12], so only the mean features related to snowmelt will be recalled here. The mass and energy balances for the snowpack are estimated as follows:

SWE_i = SWE_{i − 1} + SNOW_i − SNOWMELT_i

(21)

with

SNOWMELT_i = TF · T_i + SRF · (1 − a_i) · ISR_i if T_i > T_s

(22)

where TF is the temperature factor (mm/°C/day) and SRF is the shortwave radiation factor (m²/mm/W/day). Snowmelt is reset to 0 if the mean daily air temperature is below to the threshold temperature value.

3.2. Description of Hydrological Models

3.2.1. MOHYSE

MOHYSE (MOdèle HYdrologique Simplifié à l’Extrême) is a lumped and conceptual model that runs at the daily time step [25]. It was initially designed for academic purposes but soon found applications in research as well. It has been used in multi-model ensemble studies [33] as well as in streamflow prediction in ungauged sites [34] MOHYSE is composed of a production store and a routing store and has 10 calibration parameters (Figure 4A). The climate variables required to run MOHYSE are daily precipitation and daily mean temperature.

In MOHYSE, precipitation falls on the snowpack and the melt water is added to precipitation before evaporating, infiltrating or flowing towards the stream. The infiltrated water passes through the reservoir of the vadose zone and can either return to the atmosphere (by plant transpiration), flow to the stream or to the aquifer. Finally, the routing store allows transforming the process of the catchment’s production into discharge towards the outlet.

3.2.2. HMETS

The Hydrological Model of Ecole de Technologie Supérieure (HMETS) is a lumped and conceptual model using two connected reservoirs for the vadose and saturated zones. It has been designed to perform best in cold climates by taking advantage of its more complex snow module which requires 10 calibration parameters by itself. It has been used in multi-model averaging projects [35] and in climate change impact studies [36]. The model simulates the basic hydrological processes: evapotranspiration, infiltration, snow accumulation and melting as well as flow routing to the catchment outlet (Figure 4B).

Figure 4. Structure of the (A) MOHYSE; (B) HMETS; and (C) GR4J models.

Figure 4. Structure of the (A) MOHYSE; (B) HMETS; and (C) GR4J models.

HMETS only requires precipitation, minimum and maximum temperatures data at the daily time step. HMETS has up 21 parameters that can be optimized during the calibration. Several of these parameters can be fixed if a more parsimonious, less sensitive to equifinality model is needed.

3.2.3. GR4J

The GR4J model (modèle du Génie Rural à 4 paramètres Jounaliers - daily rural engineering model with four parameters) is a four-parameter lumped rainfall-runoff model that operates at the daily time step [37]. It has been used extensively in a wide array of research applications [38,39]. GR4J is divided into two stores: a production store and a routing store (Figure 4C). Inputs for the model for a given day are precipitation and potential evapotranspiration.

The production store has three functions. The first one is a neutralizing function of precipitation by potential evapotranspiration. The second one is an output function controlled by the routing store. This is fed by a portion of liquid precipitation after interception and is drained either by evaporation or percolation. The last one is a function of exchange with the outside of the catchment (groundwater flows). The routing store is divided into two components: the first hydrograph manages the transfer of 90% of effective rainfall towards the routing store and the second hydrograph manages the transfer of the 10% remaining towards the catchment outlet.

3.3. Models Calibration Strategy and Evaluation Method

The parameters for each SM-HM combination are calibrated with observed daily discharge data on even years, and models’ combination performance is then validated with the daily observations on odd years over the 1961–2000 period. Ten automatic calibrations are performed with different combinations of model parameters and the optimal combination of parameter values are selected for each SM-HM combination based on the Nash Sutcliffe efficiency (NS) criterion. The Shuffled Complex Evolution optimization algorithm developed at the University of Arizona (SCE-UA) is used to obtain optimal parameter values for the models [40]. SCE-UA was selected based on a comparative study of 10 algorithms for the type of model in this paper [41].

Models performance is evaluated using three standard statistical measures of agreement between observed and simulated daily streamflow. The first statistic is the Pearson’s coefficient of determination (R²) given by:

$${{\displaystyle R}}^2={\left[\frac{{\displaystyle \sum _{i=1}^N\left(O_i-\overline{O}\right)\left(S_i-\overline{S}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^N{\left(O_i-\overline{O}\right)}^2}}\sqrt{{\displaystyle \sum _{i=1}^N{\left(S_i-\overline{S}\right)}^2}}}\right]}^2$$

(23)

with

$$\overline{O}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{O}_i}$$

(24)

where O_i refers to the observed streamflow for day i; $\overline{O}$ is the mean daily observed streamflow; S_i is simulated streamflow for day i; $\overline{S}$ is mean daily simulated streamflow; and N is the number of observed or simulated streamflow values. The second statistic, the Nash-Sutcliffe model efficiency (NS; [42]), follows

$$NS=1-\left[\frac{{\displaystyle \sum _{i=1}^N{(O_i-S_i)}^2}}{{\displaystyle \sum _{i=1}^N{\left(O_i-\overline{O}\right)}^2}}\right]$$

(25)

and is a normalized statistic determining the relative magnitude of the residual variance compared to the observed data variance. The last statistic is the percent deviation of streamflow volume (D_V) given by

$$D_V=\left[\frac{{\displaystyle \sum _{i=1}^N(S_i-O_i)}}{{\displaystyle \sum _{i=1}^N(O_i)}}\right]\times 100$$

(26)

which measures the accumulation of differences in streamflow volume between simulated and observed data and is commonly used to quantify water balance errors. The optimal value is 0 with low-magnitude values indicating accurate simulation of streamflow volume. Positive values indicate model overestimation bias and negative values indicate model underestimation bias [43].

4 Results

4.1. Performance of the SM-HM Combinations in Snowmelt Flows

4.1.1. Calibration

The observed and simulated mean daily discharges and SWEs during the calibration period are graphically displayed in Figure 5 for the Mistassibi Basin. Overall, the simulated mean daily discharges by the 21 SM-HM combinations follow the mean daily trend of observed discharges in the basin. However, only the HMs with the CEMANEIGE SM are able to reproduce the magnitude of the spring peak flow. Both the SM-MOH and SM-ETS combinations underestimate the winter low flows, while the SM-GR4J combinations overestimate the winter low flows; the flows during the onset of spring snowmelt are underestimated in all the SM-HM simulations. According to the annual hydrographs of some SM-HM combinations, the timing of all snowmelt generated peak flows is successfully captured during the calibration (Figure 6). However, even though the lowest peak flow in year 1988 (i.e., Julian day 4878 in Figure 6) was captured by the SM-HM combinations, a slight underestimation (35%) of the highest peak flow in year 1976 (i.e., Julian day 2698 in Figure 6) is observed. The simulated SWEs by the SM-HM combinations generally follow the observed trend in the catchment: the maximum snow accumulation is reached at the beginning of April (i.e., Julian day around 100 in Figure 5) and the snowmelt period ends at the end of May (i.e., Julian day around 150 in Figure 5). Overall, all the SM-HM combinations are able to reproduce the shape of the observed snowgraphs over the basin. Although a generally good accordance with the observations can be observed over the study basin, some apparent weaknesses of both the DD and DD/EB models in simulating the seasonal snowpack processes can be underlined: the SMs underestimate the SWE during the snow accumulation period, while they overestimate the SWE over the snowmelt period, with higher differences between the models with the SM-GR4J combination. It is noted that the overestimation of the peak snow accumulation directly causes the overestimation of the magnitude of the spring peak flow in the basin; this is the case of the event for the year 1996 (i.e., Julian day 6356 in Figure 6).

Figure 5. Mean annual hydrographs of observed (OBS), simulated discharges and SWEs using the 21 SM-HM combinations over the calibration period (even years over 1961–2000).

Figure 5. Mean annual hydrographs of observed (OBS), simulated discharges and SWEs using the 21 SM-HM combinations over the calibration period (even years over 1961–2000).

Figure 6. Annual hydrographs of observed (OBS), simulated discharges and SWEs using the 12 SM-HM combinations over the calibration period (even years over 1961–2000).

Figure 6. Annual hydrographs of observed (OBS), simulated discharges and SWEs using the 12 SM-HM combinations over the calibration period (even years over 1961–2000).

To summarize the results of the 21 calibrations, Figure 7 presents the scatterplots between the observed and simulated daily discharges and SWEs along with 1:1 line (line of perfect fit) for the Mistassibi Basin. It can be observed that the simulated daily discharge values by the 21 SM-HM combinations are distributed uniformly about the 1:1 line for low values of observed discharge. For high values of observed discharge, the majority of the simulated daily values are below the line of perfect fit, indicating that all the SM-HM combinations underestimate the high values of daily discharge. A close observation of Figure 7 shows that the clusters of data scatter for high values of daily discharge for SM-GR4J differ from those of SM-MOH and SM-ETS combinations. From the scatterplots of snow, an over-prediction of low SWE values and an under-prediction of high SWE values by all the SM-HM combinations can be deduced.

Figure 7. Scatterplots of observed, simulated daily discharges and SWEs with the 21 SM-HM combinations during the calibration period (even years over 1961–2000). The coefficient of determination (R²) obtained for each SM-HM is also specified.

Figure 7. Scatterplots of observed, simulated daily discharges and SWEs with the 21 SM-HM combinations during the calibration period (even years over 1961–2000). The coefficient of determination (R²) obtained for each SM-HM is also specified.

The statistics of the 21 SM-HM combinations over the calibration period are presented in Table 3. The ensemble of SM-HM combinations performs well in simulating daily mean discharges with NS (R²) values above 0.78, 0.77 and 0.66 (0.89, 0.88 and 0.82) for SM-MOH, SM-ETS, and SM-GR4J, respectively, and a low value in volume biases (Dv < 2%). The best performance of the 21 SM-HM combinations at simulating daily discharges over the calibration period is obtained for CEM-MOH and CEM-ETS. It is noted that the HM combined with the degree-day CEMANEIGE snow model stands out among all SM-HMs for its ability at simulating snowmelt flows.

Figure 8. Mean annual hydrographs of observed (OBS), simulated discharges and SWEs using the 21 SM-HM combinations over the validation period (odd years over 1961–2000).

Figure 8. Mean annual hydrographs of observed (OBS), simulated discharges and SWEs using the 21 SM-HM combinations over the validation period (odd years over 1961–2000).

Figure 9. Annual hydrographs of observed (OBS), simulated discharges and SWEs using the 12 SM-HM combinations over the validation period (odd years over 1961–2000).

Figure 9. Annual hydrographs of observed (OBS), simulated discharges and SWEs using the 12 SM-HM combinations over the validation period (odd years over 1961–2000).

4.1.2. Validation

A visual comparison of the observed and simulated mean daily discharges and SWEs over the validation period for the Mistassibi Basin is presented in Figure 8. The results show a good agreement between observed and simulated mean daily discharges by the 21 SM-HM combinations, except for the events around Julian day 100 which are captured by none of the SM-HM combinations. However, the snowmelt generated peak flow is underestimated in all the SM-HM simulations. Further, the winter flows are alternatively under- and over-estimated depending of the HM considered and the events at the onset of spring are captured by none of the SM-HM combinations. Interpretation of the annual hydrographs of some SM-HM combinations during the validation period is similar to that of the results in the calibration period: the magnitude of the highest peak flows is generally underestimated by the SM-HM combinations but the timing of all spring peak flows is quite well reproduced (Figure 9). The shape of the observed snowgraphs is still quite well captured by the SM-HM combinations over the basin, with the peak SWE observed at the beginning of April (i.e., Julian day around 100 in Figure 8) and the end of the snowmelt period occurring at the end of May (i.e., Julian day around 150 in Figure 8). The snow accumulation period is underestimated and the snowmelt period is overestimated by both the DD and DD/EB models. Again, the overestimation of the peak snow accumulation in year 1995 (i.e., Julian day 5988 in Figure 9) explains the abnormally high peak flow obtained by the SM-HM combinations.

Figure 10. Scatterplots of observed, simulated daily discharges and SWEs with the 21 SM-HM combinations during the validation period (odd years over 1961–2000). The coefficient of determination (R²) obtained for each SM-HM is also specified.

Figure 10. Scatterplots of observed, simulated daily discharges and SWEs with the 21 SM-HM combinations during the validation period (odd years over 1961–2000). The coefficient of determination (R²) obtained for each SM-HM is also specified.

The scatterplots of observed and simulated daily discharges and SWEs are shown in Figure 10. Although the data scatters lying below the line of perfect fit show an overall discharge underestimation by the ensemble of SM-HM, the closeness of the data scatters to the line of perfect fit indicates a quite good performance of the 21 SM-HM combinations in estimating daily discharges for the study catchment. Similar to the calibration period, the low values of daily discharges and SWEs are overestimated, while the high values of daily discharge and SWEs are underestimated by all the SM-HM combinations. The clusters of data scatter for high values of daily discharge are comparable for both SM-MOH and SM-ETS combinations.

Table 3. Performance of the 21 SM-HM combinations over the calibration (even years) and validation (odd years) during the 1961–2000 period. The performance is expressed with the Nash-Sutcliffe model efficiency criterion (NS), the Pearson’s coefficient of determination (R²) and the percent deviation of streamflow volume.

	Calibration			Validation
	NS	R2	DV (%)	NS	R2	DV (%)
MOH-MOH	0.81	0.90	−0.57	0.82	0.91	−4.27
HBV-MOH	0.78	0.89	0.07	0.80	0.90	−4.13
CEM-MOH	0.83	0.91	−0.90	0.87	0.93	−3.54
ETS-MOH	0.81	0.90	−0.48	0.80	0.90	−4.54
HYD-MOH	0.80	0.90	−0.96	0.80	0.90	−4.97
SEB-MOH	0.81	0.90	−0.50	0.82	0.91	−4.27
ETI-MOH	0.80	0.90	−1.47	0.82	0.91	−4.64
MOH-ETS	0.82	0.91	−0.52	0.85	0.92	−2.74
HBV-ETS	0.77	0.88	0.43	0.79	0.89	−3.63
CEM-ETS	0.83	0.91	−0.26	0.87	0.93	−3.47
ETS-ETS	0.80	0.89	−1.60	0.81	0.90	−3.99
HYD-ETS	0.79	0.89	−0.44	0.80	0.90	−3.84
SEB-ETS	0.82	0.90	−0.80	0.84	0.92	−3.58
ETI-ETS	0.80	0.89	0.05	0.83	0.91	−2.95
MOH-GR4J	0.74	0.86	−0.54	0.77	0.89	−5.34
HBV-GR4J	0.66	0.82	−1.40	0.73	0.86	−5.33
CEM-GR4J	0.79	0.89	−1.06	0.85	0.92	−4.79
ETS-GR4J	0.67	0.83	1.84	0.74	0.86	−2.73
HYD-GR4J	0.71	0.85	−0.43	0.75	0.87	−5.22
SEB-GR4J	0.74	0.86	−0.59	0.79	0.89	−5.44
ETI-GR4J	0.71	0.85	−1.82	0.76	0.88	−5.00

Table 3. Performance of the 21 SM-HM combinations over the calibration (even years) and validation (odd years) during the 1961–2000 period. The performance is expressed with the Nash-Sutcliffe model efficiency criterion (NS), the Pearson’s coefficient of determination (R²) and the percent deviation of streamflow volume.

	Calibration			Validation
	NS	R2	DV (%)	NS	R2	DV (%)
MOH-MOH	0.81	0.90	−0.57	0.82	0.91	−4.27
HBV-MOH	0.78	0.89	0.07	0.80	0.90	−4.13
CEM-MOH	0.83	0.91	−0.90	0.87	0.93	−3.54
ETS-MOH	0.81	0.90	−0.48	0.80	0.90	−4.54
HYD-MOH	0.80	0.90	−0.96	0.80	0.90	−4.97
SEB-MOH	0.81	0.90	−0.50	0.82	0.91	−4.27
ETI-MOH	0.80	0.90	−1.47	0.82	0.91	−4.64
MOH-ETS	0.82	0.91	−0.52	0.85	0.92	−2.74
HBV-ETS	0.77	0.88	0.43	0.79	0.89	−3.63
CEM-ETS	0.83	0.91	−0.26	0.87	0.93	−3.47
ETS-ETS	0.80	0.89	−1.60	0.81	0.90	−3.99
HYD-ETS	0.79	0.89	−0.44	0.80	0.90	−3.84
SEB-ETS	0.82	0.90	−0.80	0.84	0.92	−3.58
ETI-ETS	0.80	0.89	0.05	0.83	0.91	−2.95
MOH-GR4J	0.74	0.86	−0.54	0.77	0.89	−5.34
HBV-GR4J	0.66	0.82	−1.40	0.73	0.86	−5.33
CEM-GR4J	0.79	0.89	−1.06	0.85	0.92	−4.79
ETS-GR4J	0.67	0.83	1.84	0.74	0.86	−2.73
HYD-GR4J	0.71	0.85	−0.43	0.75	0.87	−5.22
SEB-GR4J	0.74	0.86	−0.59	0.79	0.89	−5.44
ETI-GR4J	0.71	0.85	−1.82	0.76	0.88	−5.00

From Table 3, it can be seen that the daily discharges are well captured over the validation years, with NS (R²) values above 0.80, 0.79 and 0.73 (0.90, 0.89 and 0.86) for SM-MOH, SM-ETS, and SM-GR4J. All the SM-HM combinations underestimate the daily mean discharges (Dv < −6%). Again, the CEM-HM combination provides the best performance in simulating daily discharges over the validation period. The ensemble of SM-HM combinations over the validation period shows comparable or slightly better performances to the calibration period. As indicated by the values of NS, R², and the visual interpretation of the simulated hydrographs in Figure 9, the SM-HM combinations have a better overall performance in the validation years than in the calibration years. This might result from the observations of discharge, with a high proportion of missing data and inconsistencies for the winter months over the entire 1961–2000 period.

Figure 11. Relative difference (%) of the investigated indicators (the mean high flow (MF), the annual extreme flood (AF) and the timing of the annual extreme flood (TAF)) for each snow model over the 1961–2000 period. The boxes indicate the ensemble simulations and the color dots specify the three hydrological models; (A) MF; (B) AF; (C) TAF;

Figure 11. Relative difference (%) of the investigated indicators (the mean high flow (MF), the annual extreme flood (AF) and the timing of the annual extreme flood (TAF)) for each snow model over the 1961–2000 period. The boxes indicate the ensemble simulations and the color dots specify the three hydrological models; (A) MF; (B) AF; (C) TAF;

4.2. Uncertainty Analysis

This section provides a closer look at the model structure uncertainties on the snowmelt flow simulations by comparing the effects of the individual SM component with the uncertainty range introduced by the entire HMs. For this purpose, three hydrologic indicators are selected: the mean high flow (MF) estimated during the snowmelt period from March to June, the annual extreme flood (AF) defined as one-day maximum flow that occurs in spring due to snowmelt and the timing of the annual extreme flood (TAF).

4.2.1. Uncertainty of Snow Models on Snowmelt Flows

Figure 11 presents the relative differences of the simulated indictors compared to the observed MF, HF and TAF over the Mistassibi Basin. Each SM can be distinguished within each SM approach (DD and DD/EB). The results show that the SMs’ ensemble leads to quite similar relative difference values of TAF when compared to the observations over the basin, with a low variability between the SMs. However, for MF and AF, even though the SMs’ ensemble does almost reach the observed values of relative difference, the SMs are somewhat less consistent between them over the basin. This is particularly noticeable for AF, where the large variability of the indicator is caused by individual models in the SM approach. Additionally shown in Figure 11 is that the SMs generally provide comparable spread of the relative difference in TAF over the 1961–2000 period in the basin. The similar representation of the timing of the spring peak flow by the SM-HM combinations over the study period confirms this aspect (Figure 5 and Figure 8). Nevertheless, the same conclusion cannot be transposed to the MF and AF relative difference values, where the SMs display higher variability over the 1961–2000 period in the basin.

The choice of the SM approaches (DD versus DD/EB models) has an influence on the simulated discharges depending of the hydrologic indicator. For the study basin, the highest magnitude of variability in the simulated discharges for both the DD and DD/EB models is observed for MF and AF, whereas the lowest range of variability in the simulated discharges for both the DD and DD/EB models is found for TAF. The uncertainty associated with the two SM approaches on the simulated discharges is therefore more notable in the statistics of quantification than in the statistics of duration over the study basin. Overall, we can state that the DD/EB models generally provide the lowest range of variability in the simulated relative difference of discharge over the 1961–2000 period.

The uncertainty associated with the choice of the model within each SM approach is therefore higher than that associated with the choice of the SM approach in simulating snowmelt flows. This difference is illustrated by the MF and AF values of relative difference over the basin, where individual models in the SM approach give different results against the other models in the approach considered. In general, the simulations derived from the individual SMs within the DD/EB approach are more consistent than those within the DD approach. This study clearly emphasizes the importance of using several models within the SM approach to adequately depict the spring-snowmelt-generated peak flow. Finally, it is important to specify that the SMs contribute to the overall uncertainty in the hydrological simulations, but with a lower intensity than expected over such typical Nordic environment.

Figure 12. Relative difference (%) of the investigated indicators (the mean high flow (MF), the annual extreme flood (AF) and the timing of the annual extreme flood (TAF)) for each hydrological model over the 1961–2000 period. The boxes indicate the ensemble simulations and the color dots specify the seven snow models, (A) MF; (B) AF; (C) TAF.

Figure 12. Relative difference (%) of the investigated indicators (the mean high flow (MF), the annual extreme flood (AF) and the timing of the annual extreme flood (TAF)) for each hydrological model over the 1961–2000 period. The boxes indicate the ensemble simulations and the color dots specify the seven snow models, (A) MF; (B) AF; (C) TAF.

4.2.2. Uncertainty of Hydrological Models on Snowmelt Flows

Figure 12 presents the relative differences of simulated indicators compared to observed MF, AF and TAF over the Mistassibi Basin. Here, the distinction is made for each HM. Regarding the results from individual models, GR4J encompasses the observed values of relative difference for the three indicators in the basin. MOHYSE and HMETS give different relative values of MF and AF against the observations. The three models show comparable relative difference value of TAF. As mentioned previously, the similar results obtained for TAF by the HMs are supported by the fact that the models capture adequately the timing of the spring peak flow in the basin (Figure 5 and Figure 8).

The analysis of the hydrologic indicators also shows a comparable behaviour between MOHYSE and HMETS, both in the magnitude of variability and in the direction of simulated discharges. This finding cannot be extended to GR4J, where the HM displays higher variability of MF and AF relative difference values by almost 5% over the 1961–2000 period in the basin. This means that the hydrological simulations performed with GR4J are made with a greater uncertainty than those performed with MOHYSE and HMETS. From this analysis, we can conclude that the choice of the HM has an influence on the snowmelt flow simulations, with a higher uncertainty related to the statistics of quantification (MF and AF) than to the statistics of duration (TAF). However, the present study is only based on three HMs, which naturally limits the HM estimates’ statistical robustness. In this aspect, this study emphasizes the need to consider additional HMs in order to provide a more comprehensive multi-model approach. This new experiment design should allow reinforcing our results about the exploration of the model structure uncertainty on snowmelt flows—uncertainty that has to be taken into account in hydrological studies.

5. Discussion

As mentioned in the introduction, the understanding, quantification and reduction of uncertainty are three critical aspects to be considered in order to adequately address uncertainty in hydrologic modeling uncertainty [3]. Considering these elements, our work demonstrates that the SM components are not an important source of uncertainty when compared to the entire HMs for snowmelt flow simulations over the study Nordic catchment. This is in agreement with the companion work of Troin et al. [22], where a comparison of different responses in SWE simulations by an ensemble of eight SMs is conducted over three Canadian catchments. The authors found a low degree of uncertainty of the eight SMs in simulating SWE under present and future climates. However, additional studies will need to investigate if the conclusions still hold for the simulated SWE over high-elevation snow-covered catchments and, by extension, for the streamflow generation over such typical Nordic environments.

As the other hydrological model component that is directly associated with changing climatic conditions (i.e., precipitation and temperature) is potential evapotranspiration [44], a complementary step of this study will be to compare a set of potential evapotranspiration (PET) computation methods. Evaluating how the SMs and PET computation methods affect the streamflow simulations might contribute to a better understanding of HM structure uncertainty. In addition, this information will help limit the total uncertainty in hydrological modeling, and probably confirm the greatest uncertainty of model structure compared to model parameters [4,5]. A future work conducted with the present experimental design of the multi-approach will be dedicated to this aspect.

6. Conclusions

The multi-model approach is becoming a common procedure for considering uncertainty in hydrological simulations. Of the many sources of uncertainty that may affect a hydrological simulation, uncertainty from the structure of the hydrological model is usually pinpointed as the most crucial uncertainty. This has led to one popular method to implement hydrological ensembles: pooling the outputs of a group of hydrological models for accounting for the uncertainties related to the structure of the models. This study uses seven snow models (four DD models and three DD/EB models) combined to three hydrological models to evaluate snowmelt-dominated streamflow simulations over a Canadian catchment. The source of uncertainty related to the model structure is considered, either for individual component (SM) or entirely (HM). The analysis is carried out on the 1961–2000 historical period based on a comprehensive hydrometeorological dataset. The main results obtained from the present multi-model approach are summarized as follows.

(1)

For the study period and the basin considered, the 21 SM-HM combinations have similar satisfactory discharge simulation performances during the calibration and validation periods, with R² values above 0.8 and NS values above 0.7.

(2)

The ensemble of 21 SM-HM combinations simulates daily discharges and SWEs that follow the trends of observed values in the Mistassibi Basin. All the SM-HM combinations underestimate the spring-snowmelt-generated peak flow due to the underestimation of the peak snow accumulation by the SMs.

(3)

Both the DD SM-HM and DD/EB SM-HM approaches agree fairly well regarding the simulation of snowmelt flows over the Mistassibi Basin. This level of agreement implies that the simple DD SM models can be considered to be a tool as effective as the more sophisticated DD/EB SM models for simulating snowmelt flows at the catchment scale.

(4)

The present investigation has permitted an estimation of the source of uncertainty associated to model structure in discharge simulations. The analysis of specific hydrologic indicators shows that the main uncertainties result from the hydrological models followed by the snow models. The uncertainty related to the choice of the given HM cannot be neglected in snow hydrological studies. The selection of the SM plays a more important role than the choice of the SM approach (DD versus DD/EB) in simulating snowmelt flows.

The results obtained in this study provide key directions for evaluating uncertainty associated with model structure in hydrological simulations. In the future, it will be necessary to perform additional tests to support the results and conclusions drawn from the present multi-model approach. The comparison of uncertainties related to the most sensitive components of hydrological models to changing climatic conditions (SMs and PET computation methods) will be particularly useful to increase the reliability of the model structure, and thus to maximally reduce uncertainty in hydrological impact studies.