Data and methods

Study site and data Twelve tributary basins of the St. Lawrence River in Quebec were selected in this study (Fig. 1). The choice was based on the length of observed discharge data (> 20 years), the natural character of the hydrological regime of the ri vers and, especially, the availability of snow measurement points inside or close to the studied basin. The area of the basins varies between 367 and 4504 km2 (Table 1). They are located in four homogeneous hydrological regions, namely (i) the St. Lawrence northwest region (Batiscan, Bras du Nord, Matawin) on the north shore and characterized by a continental climate; (ii) the St. Lawrence southwest region (Nicolet, L’Acadie) characterized by a mixed maritime and continental climate; (iii) the St. Lawrence southeast region (York, Beaurivage, Bécancour, Famine, Etchemin, Ouelle) characterized by a mix of maritime and continental climate, and (iv) the St. Lawrence northeast (Godbout) characterized by a maritime climate (Assani et al., 2010a; Assani et al., 2010b; Mazouz et al., 2013). Daily historical discharge data measured at the basin outlets were extracted from the website of the Quebec Center of Water Expertise (CEHQ) (www.cehq.qc.ca). The climatic data used were extracted from the daily climate grids developed by the Atmospheric Environment Information Service (SIMAT) in collaboration with the CEHQ (Bergeron, 2015). The total daily precipitation (solid and liquid), minimum and maximum temperature at each grid point are estimated by spatial interpolation (kriging) using stations managed by the Quebec Climate Monitoring Program (Programme de surveillance du climat du Québec: PSC) and stations operated by the national hydropower company Hydro-Québec. The interpolated climate grids are considered to be of good quality in the southem part of Quebec given the high density of stations in this region. Continuous data are available for the period from 1961 to 2015 (Bergeron, 2015). The snow observations come from the permanent snow survey network maintained by the Ministry of Environment and Fight Against Climate Change (MELCC). Survey sites are located in forested areas where the depth and density of snow is measured every two weeks during the winter and spring seasons. Each SWE observation at a given survey site represents the average of ten manual measurements made with a snow tube along a 100 rn-long transect. The historical measurements of 12 survey sites located in or very close to the selected basins were used. The snow survey sites and the characteristics of the basins are shown in Table 1.

Calibration strategies In the calibration strategy, the choice of optimization algorithm and objective functions have a great influence on the identification of the parameters and their uncertainty, in addition to the structure of the model and the quality of the input data. The objective of improving the snow coyer simulation by the GR4J -Cemaneige model and the identifiability of the parameters led us to consider four calibration strategies in which three algorithms (local, global and multi-objective) were used: for the first two strategies, the calibration of the six GR4J-Cemaneige model parameters was done simultaneously using the observed discharge, with respectively a local ‘pas-à-pas’ approach (strategy 1: ‘LOCAL’) and the global SCE-UA algorithm (strategy 2: ‘SCEFLOW’). The local ‘pas-à-pas’ (step-by-step) optimization algorithm was used for the development and improvement of the two models by researchers at Cemagref, France (Perrin, 2000). The global automatic algorithm Shuffled Complex Evolution (SCE-UA) developed by Duan et al. (1992) is considered the most efficient by hydrologists. The SWE simulated by these two strategies was then compared with the SWE observed at the nivometric survey points in, or closest to, the selected basins, for the elevation band closest to the elevation of the snow survey point. In the third strategy, ‘SCE_INDEP’, the four hydrological model parameters and the two snow parameters were calibrated separately. The two snow parameters were calibrated with the observed SWE and prescribed subsequently as fixed parameters in the coupled GR4J-Cemaneige model. The four hydrological parameters were next calibrated with the observed discharge by the global SCE-UA algorithm. The fourth and final multi-objective strategy, ‘MULTI’, was finally applied in which the six model parameters were calibrated simultaneously with the observed SWE and discharge using the multiobjective optimization algorithm AMALGAM (Vrugt et al., 2007). The split-sample test calibration method (Klemes, 1986) was used to separate the observation record into two equal length periods of calibration/validation (Table 1). The Nash-Sutcliffe (1970) efficiency criterion was used as the evaluation criterion and calculated from the observed and simulated discharge (Nash-Q) and from the SWE measured at surveys point and that simulated by the model for the corresponding elevation band (Nash-SWE):

Models performance

For all calibration strategies, one global optimal parameter set was obtained, except for the multi-objective method MULTI where a vector of ‘Pareto-optimal’ parameter sets was obtained. The optimal point of the Pareto set that gives the best compromise between snow SWE and discharge was selected to compare with the performance of the other calibration strategies. Boxplots in Fig. 2 show the model performance for the 12 selected basins and the four calibration methods for both the calibration and validation periods. With the first method LOCAL (Fig. 2aI) the model shows a good performance in discharge simulation with a median Nash-Q value of 84% in calibration and 80% in validation but a comparatively poor performance in SWE simulation: the median Nash-SWE value is around 40% during calibration and increases to 58% in validation, but the spread of the distribution is higher (Fig. 2a2). The calibration results with the global algorithm SCE_FLOW also show a good streamflow simulation with a median Nash-Q value of 82% in calibration, similar to the first method, but with an improvement in validation with the median Nash-Q value increasing to 83% (Fig. 2bl). The SWE simulation is still weak with the SCE FLOW method but slightly improves nonetheless with the median Nash-SWE value increasing to 43% in calibration and 59% in validation (Fig. 2b2), but the spread of the distributions is still considerable and comparable to the first local method. In the third strategy SCE_INDEP, the two Cemaneige parameters were calibrated separately against the observed SWE, which improved significantly the Nash-SWE in calibration with a median value of 58%. This significant improvement in the Nash-SWE criterion is however accompanied by a degradation of the flow simulation, with the median Nash-Q value decreasing to 80% in calibration and 76% in validation, which means that the introduction of two fixed snow parameters in the GR4J-Cemaneige model led to an over-adjustment on the SWE and a degradation in the streamflow simulation (Fig. 2c).

Parameters interaction

The distribution of the optimal parameters obtained for the 12 basins is different between the calibration strategies, except for the parameter x4 (base of the unit hydrograph) and x2 (water table exchange) which are more similar (Fig. 5). The choice of the calibration method does not seem to have a clear impact on the distribution of the parameters between the basins. For the third method SCE_INDEP, where the calibration of the Cemaneige parameters is do ne separately on the SWE observations, the degree-day factor x6 which determines potential snowmelt is higher and more dispersed than for the other methods, while x5 (weighting coefficient of the thermal state), which determines the onset of melting, is minimized. Subsequently, the parameter xl (maximum capacity of the soil reservoir) approaches its lower limit in sorne basins while x3 (the routing reservoir capacity) approaches its upper limit. The optimal parameters obtained by each calibration strategy differ with each other for a given basin (Fig. 6). This shows that different optimal sets give the best simulation for each calibration method over the same period. The two parameters x2 and x4 are the only ones for which similar values are obtained from aU calibration strategies. Equifinality was studied usmg two calibration methods, namely SCE FLOW, which is the global algorithm most used in the literature and uses only observed discharge, and the method MULTI which gave the best performance (Fig. 3). The distributions of equifinal parameters obtained for the 12 basins by the procedure described in the methods section are presented in Fig. 7. The distributions for the SCE_FLOW method are narrow but with many extreme values, while the distributions for the multi-objective method are more homogeneous (Fig. 7). For the SCE FLOW method the number of sets found is very high and variable between basins. Therefore, in order to converge mathematically with a very small change in the Nash value « 1 %), the maximum capacity of the soil reservoir (xl) can, for example, vary from 200 to 450 mm for basin 1 (Batiscan). Similar variability is seen for the other basins.