Assessing River Low-Flow Uncertainties Related to Hydrological Model Calibration and Structure under Climate Change Conditions

Low-flow is the flow of water in a river during prolonged dry weather. This paper investigated the uncertainty originating from hydrological model calibration and structure in low-flow simulations under climate change conditions. Two hydrological models of contrasting complexity, GR4J and SWAT, were applied to four sub-watersheds of the Yamaska River, Canada. The two models were calibrated using seven different objective functions including the Nash-Sutcliffe coefficient (NSEQ) and six other objective functions more related to low flows. The uncertainty in the model parameters was evaluated using a PARAmeter SOLutions procedure (PARASOL). Twelve climate projections from different combinations of General Circulation Models (GCMs) and Regional Circulation Models (RCMs) were used to simulate low-flow indices in a reference (1970–2000) and future (2040–2070) horizon. Results indicate that the NSEQ objective function does not properly represent low-flow indices for either model. The NSE objective function applied to the log of the flows shows the lowest total variance for all sub-watersheds. In addition, these hydrological models should be used with care for low-flow studies, since they both show some inconsistent results. The uncertainty is higher for SWAT than for GR4J. With GR4J, the uncertainties in the simulations for the 7Q2 index (the 7-day low-flow value with a 2-year return period) are lower for the future period than for the reference period. This can be explained by the analysis of hydrological processes. In the future horizon, a significant worsening of low-flow conditions was projected.


Introduction
Climate change is expected to have a direct impact on water resources.The increase in extreme weather events will particularly intensify both the severity and frequency of low flows (e.g., [1][2][3][4]).A wide range of applications such as water supply, irrigation, navigation, and hydropower production can be strongly affected by water shortages in rivers.Therefore, it is crucial for low-flow management that prediction about low flows and the impacts of climate change become available.However, significant uncertainties remain in climate change studies on water resources.Typically, studies evaluating the impact of climate change on low flows apply outputs (temperature and precipitation) from global or regional climate circulation models to hydrological models.Thus, uncertainties have many sources related to emission scenarios, General Circulation Models (GCMs), downscaling methods such as Regional Circulation Models (RCMs), and hydrological models (structure and parameters).Assessing different sources of uncertainty in river flow projections has become a topic of interest (e.g., [5][6][7][8][9]).Nevertheless, only a few aspects of the uncertainty related to hydrological model structure and parameters in low-flow projections have been studied.Najifi et al. [10] used eight GCMs and two emission scenarios with four hydrological models of different complexity calibrated with three objective functions to evaluate the uncertainties associated with hydrological models in a climate change impact study of a watershed in Oregon, USA.The results showed that total uncertainty and uncertainties related to hydrological models are higher for the dry season than for the wet season.Velàzquez et al. [11] analyzed the hydrological response of four hydrological models driven by five members of a RCM for a watershed in Canada and three members of a RCM for a watershed in Germany.Results showed that low-flow projections are strongly affected by the structure of the hydrological model.Vansteenkiste et al. [12] applied more than 50 runs of global/regional climate models to five hydrological models (lumped conceptual to physically based).They also compared a manual and an automatic calibration for a watershed in Belgium.They found that the structure and the calibration of the model highly influenced low-flow projections resulting in higher uncertainties related to the hydrological model for low flows than for peak flows.Dams et al. [13] studied three climate change scenarios with four physically-based hydrological models calibrated with four objective functions for a watershed in Belgium.They concluded that, although uncertainty related to climate change scenarios carry the largest source of uncertainties, the hydrological model structure contributes considerable uncertainty to low-flow projections.Parajka et al. [14] used 4 climate change scenarios with a conceptual semi-distributed model calibrated with 11 objective functions in 3 different decades to assess the uncertainty in low-flow projections for 262 watersheds in Austria.The results show that climate scenarios have the most impact on the total uncertainty; although the uncertainty related to hydrological model calibration is not negligible.They also found that relative contributions to the uncertainty vary with the hydrological regime (summer low-flow or winter low-flow).In these studies, many authors pointed out that calibrated models have difficulty in simulating low flows compared to observations.This may be related to uncertainties in the parameters of the hydrological model or in the hydrological processes and state variables simulated by the model.The objective function used to calibrate the model can also introduce uncertainty on low-flow simulation.
To our knowledge, no studies assessed parameter uncertainty related to choosing the objective function of the calibration for low-flow simulations.Parameter uncertainty stems from the equifinality of the model, where multiple hydrological model parameter ensembles produce similar or acceptable model outputs (example river flows) [15], but also from the calibration method used, which includes the objective function used [14], as well as model input and output errors [16].Different methods, each with their strengths and weaknesses [17,18], allow the estimation of the parameter uncertainties for the hydrological model during calibration, including the Bayesian approach [19], Generalized Likelihood Uncertainty Estimation (GLUE) [20,21], Markov chain Monte Carlo (MCMC) [19], sequential uncertainty fitting (SUFI-2) [22], and parameter solution (PARASOL) [23].
This study aims to answer three questions.Are hydrological model parameter uncertainties influenced by using an objective function targeted to low flows?Are the impacts of the objective function on the overall hydrological model uncertainty (model parameters, structure and calibration) the same for a global model as for a distributed model?What is the order of magnitude of this uncertainty compared with the uncertainty of GCMs in current and future climates?
To address these questions, a global hydrological model, GR4J, and a distributed model, SWAT, were used to simulate low flows for the Yamaska River watershed, in Québec, Canada.Four sub-watersheds were studied since they display differences in area and land cover.Some of these sub-watersheds currently have substantial low-flow issues.Six objective functions related to low flows were used, along with the Nash-Sutcliffe efficiency (NSE Q ) [24], which is most often used to calibrate hydrological models.A method similar to the Parameter Solution (PARASOL) [23] was used to evaluate parameter uncertainty.This method selects "behavioral" results from simulations generated through a calibration process based on an objective threshold.A group of "behavioral" parameters is then selected for each calibration process.The uncertainty in streamflows, hydrological processes and state variables resulting from this group of parameters was analyzed over a reference  and a future (2040-2070) period.Twelve climatic projections from different combinations of RCMs/GCMs were used to evaluate the uncertainty in climate projections.An analysis of variance (ANOVA) [8] was used to evaluate the relative contributions to overall uncertainty of climate model, hydrological model, objective function, hydrological model parameters, as well as their interactions.

Study Area
The Yamaska River watershed is located in the southeast of the province of Quebec, Canada [25].It covers an area of 4784 km 2 , divided between the Appalachian Mountain region upstream and the lowlands of the St. Lawrence downstream.The main rivers of the watershed are the Yamaska River, Southeast Yamaska River and Noire River.Three reservoirs are located in the Yamaska River watershed: Choiniere (4.70 km 2 ), Waterloo (1.50 km 2 ) and Brome Lake (14.53 km 2 ).These reservoirs are operated to mitigate spring flooding and to ensure a minimum streamflow for the municipal water supply.Agriculture plays a predominant role in the economy of the Yamaska River watershed.Indeed, about 52% of the territory is dedicated to agriculture and pasture land.The forested area covers about 43% of its territory, urban areas 3%, and water areas and wetlands cover 2% (Figure 1).The soil in the Yamaska River watershed consists largely of sandy silt in its upstream part and silty sand to clay in the downstream part.
Climate 2017, 5, 19 3 of 23 over a reference  and a future (2040-2070) period.Twelve climatic projections from different combinations of RCMs/GCMs were used to evaluate the uncertainty in climate projections.An analysis of variance (ANOVA) [8] was used to evaluate the relative contributions to overall uncertainty of climate model, hydrological model, objective function, hydrological model parameters, as well as their interactions.

Study Area
The Yamaska River watershed is located in the southeast of the province of Quebec, Canada [25].It covers an area of 4784 km 2 , divided between the Appalachian Mountain region upstream and the lowlands of the St. Lawrence downstream.The main rivers of the watershed are the Yamaska River, Southeast Yamaska River and Noire River.Three reservoirs are located in the Yamaska River watershed: Choiniere (4.70 km 2 ), Waterloo (1.50 km 2 ) and Brome Lake (14.53 km 2 ).These reservoirs are operated to mitigate spring flooding and to ensure a minimum streamflow for the municipal water supply.Agriculture plays a predominant role in the economy of the Yamaska River watershed.Indeed, about 52% of the territory is dedicated to agriculture and pasture land.The forested area covers about 43% of its territory, urban areas 3%, and water areas and wetlands cover 2% (Figure 1).The soil in the Yamaska River watershed consists largely of sandy silt in its upstream part and silty sand to clay in the downstream part.Water quality is a serious issue for municipalities.Indeed, given the high percentage of farmland, various pollutants are released and reach the hydrographic network of the watershed.This problem becomes even more acute during periods of low flow as the dilution of contaminants is reduced.
In general, the climate of southern Quebec is humid continental with large seasonal variations.Annual precipitation averages (computed from 1980-2010) around 1100 mm.For summer months (May to October), rainfall averages around 650 mm.This region receives close to 200 cm of snowfall over an average winter, and snow usually remains on the ground from December to March.Summers are fairly warm and humid in this region with average daily temperatures near 20 °C and Water quality is a serious issue for municipalities.Indeed, given the high percentage of farmland, various pollutants are released and reach the hydrographic network of the watershed.This problem becomes even more acute during periods of low flow as the dilution of contaminants is reduced.
In general, the climate of southern Quebec is humid continental with large seasonal variations.Annual precipitation averages (computed from 1980-2010) around 1100 mm.For summer months (May to October), rainfall averages around 650 mm.This region receives close to 200 cm of snowfall over an average winter, and snow usually remains on the ground from December to March.Summers are fairly warm and humid in this region with average daily temperatures near 20 • C and average monthly precipitation around 100 mm.Spring and fall are changeable seasons, prone to extremes in temperature and precipitation.
Four hydrometric stations are located in the watershed at the outlet of each sub-watershed (Figure 1).Noire River (1414 km 2 ) and Cowansville (214 km 2 ) are upstream sub-watersheds.The Cowansville sub-watershed drains into the Farnham sub-watershed (1202 km 2 ).They all drain into the St-Hyacinthe sub-watershed (3289 km 2 ).Characteristics of each sub-watershed are presented in Table 1.

Methods
The study was carried out using two continuous hydrological models, GR4J and SWAT, which are briefly described below.

GR4J and SWAT Model Descriptions
The GR4J model (Figure 2a) is a lumped rainfall-runoff model [26] using two conceptual reservoirs (S and R) and four calibration parameters: X1, the exchange term coefficient; X2, the maximum capacity of the routing reservoir; X3, the characteristic time of the unit hydrographs; and X4, the maximum capacity of the soil reservoir.Daily time series of precipitation (P) and potential evapotranspiration (PET) are required as inputs.In this study, PET (Equation ( 1)) was computed using the Hamon model [27], and P includes the snowmelt using a simple degree-day model.In addition to the four calibration parameters of the GR4J model, three parameters were added: X5, a multiplicative adjustable factor to calculate the actual evapotranspiration (ET) from PET, and two parameters (X6, a melt temperature threshold ( • C); and X7, a melt rate factor (mm/day)) for the snowmelt model based on a degree-day approach.
where H t (day) is the average number of daylight hours per day, e s is the saturation vapor pressure (kPa), and T mean is the daily mean temperature ( • C).The Soil and Water Assessment Tool (SWAT) [28] (Figure 2b) is a semi-distributed and semi-conceptual model that simulates water, nutrient and pesticide transport at the watershed scale on a daily time step.The watershed is subdivided into sub-watersheds, river reaches and hydrological response units (HRUs).HRUs are generated by combining soil type, land use and slope classes within sub-watersheds to obtain areas which are considered hydrologically homogeneous.The SWAT simulates the hydrologic cycle based on the water balance equation (Equation ( 2)) [28]: where SW t is the soil water content (mm) at time t, SW 0 is the initial soil water content (mm) at time i, t is the time (days), R day is the amount of precipitation (mm), Q surf is the amount of surface runoff (mm), ET is the amount of actual evapotranspiration, w seep is the amount of water percolating from the unsaturated to the saturated zone (mm), and Q gw is the amount of return flow (mm).Snowmelt is calculated as a function of air temperature, snow cover area, a melt factor, a base temperature threshold, and the temperature of the snow pack.Snowmelt is added to rainfall to be infiltrated and/or to contribute to surface runoff.The Soil Conservation Service (SCS) runoff equation is used to compute runoff Q surf .The Penman-Monteith method is used to compute PET.Actual evapotranspiration varies depending on PET, interception by the canopy, and the amount of water in the unsaturated soil.
Routing is simulated with the variable storage coefficient method.SWAT can also simulate nutrient and pesticide transport, but these model components were not used in this study.Details on the hydrological processes in SWAT can be found in the SWAT manual [29].
Climate 2017, 5, 19 5 of 23 in the unsaturated soil.Routing is simulated with the variable storage coefficient method.SWAT can also simulate nutrient and pesticide transport, but these model components were not used in this study.Details on the hydrological processes in SWAT can be found in the SWAT manual [29].

Low-Flow Indices
Low-flow indices calculated from streamflow measurements were used to evaluate the impact of climate change on low flows.In the USA and Canada, the most widely used indices are the lowest average flows that occur over a consecutive n-day period for m-year return periods (nQm).7Q2, 7Q10 and 30Q5 indices are used by the Ministère du Développement Durable, de l'Environnement et de la Lutte contre les changements climatiques (MDDELCC) of Quebec [30] for the evaluation of low flows in environmental studies.These low-flow indices are related to river conditions associated with environmental issues, such as pollutant dilution or fish habitat quality.In addition, the MDDELCC limits river water withdrawals by municipalities to 15% of 7Q2 [31].The 7Q2 and 7Q10 indices used in this study were determined by fitting the 7-day lowest streamflow for each time period to a Log Pearson Type III distribution for 2-year and 10-year return periods respectively.The Log Pearson Type III distribution was found to better represent the observed distribution compared to other distributions (e.g., Gumbel, lognormal).

Calibration, Validation and Statistical Tests
In this study, the GR4J and SWAT models were calibrated using historical river flow observations for the period from 1 January 1998, to 31 December 2009.Snow water equivalent (SWE) data were available at only one station near Brome Lake (see Figure 1) at a bi-weekly frequency during the months of January-April and consequently were not used in the calibration process.The validation period spans from 1 January 1985 to 31 December 1997.A method based on the PARAmeter SOLutions procedure (PARASOL) [23] was used to perform optimization and parameter uncertainty analyses.The objective function was minimized using the shuffled complex evolution method developed at The University of Arizona (SCE-UA) [32].The simulations performed by the SCE-UA optimization were divided into "behavioral" simulations and "non-behavioral" simulations similar to the GLUE methodology [21].According to [21], behavioral models are those models that provide good fits to any observables available, while non-behavioral models do not.The SCE-UA samples over the entire parameter space with a focus on solutions near the optimum.A threshold was used to select "behavioral" simulations.This study used a bootstrap

Low-Flow Indices
Low-flow indices calculated from streamflow measurements were used to evaluate the impact of climate change on low flows.In the USA and Canada, the most widely used indices are the lowest average flows that occur over a consecutive n-day period for m-year return periods (nQm).7Q2, 7Q10 and 30Q5 indices are used by the Ministère du Développement Durable, de l'Environnement et de la Lutte contre les changements climatiques (MDDELCC) of Quebec [30] for the evaluation of low flows in environmental studies.These low-flow indices are related to river conditions associated with environmental issues, such as pollutant dilution or fish habitat quality.In addition, the MDDELCC limits river water withdrawals by municipalities to 15% of 7Q2 [31].The 7Q2 and 7Q10 indices used in this study were determined by fitting the 7-day lowest streamflow for each time period to a Log Pearson Type III distribution for 2-year and 10-year return periods respectively.The Log Pearson Type III distribution was found to better represent the observed distribution compared to other distributions (e.g., Gumbel, lognormal).

Calibration, Validation and Statistical Tests
In this study, the GR4J and SWAT models were calibrated using historical river flow observations for the period from 1 January 1998, to 31 December 2009.Snow water equivalent (SWE) data were available at only one station near Brome Lake (see Figure 1) at a bi-weekly frequency during the months of January-April and consequently were not used in the calibration process.The validation period spans from 1 January 1985 to 31 December 1997.A method based on the PARAmeter SOLutions procedure (PARASOL) [23] was used to perform optimization and parameter uncertainty analyses.The objective function was minimized using the shuffled complex evolution method developed at The University of Arizona (SCE-UA) [32].The simulations performed by the SCE-UA optimization were divided into "behavioral" simulations and "non-behavioral" simulations similar to the GLUE methodology [21].According to [21], behavioral models are those models that provide good fits to any observables available, while non-behavioral models do not.The SCE-UA samples over the entire parameter space with a focus on solutions near the optimum.A threshold was used to select "behavioral" simulations.This study used a bootstrap method [33] applied to the best solution to define a distribution of the objective function.The threshold for the behavioral simulations was set to a confidence interval of 95%.
The calibration of the GR4J model was performed for each of the four sub-watersheds individually since the GR4J model does not simulate river routing from one sub-watershed to another as it is a lumped (i.e., global) model.For the SWAT model, the calibration was done from upstream to downstream.The Cowansville and Noire River sub-watersheds were calibrated first.The "best" parameters were selected and set to these two sub-watersheds before resuming the calibration process of the Farnham sub-watershed.The St-Hyacinthe sub-watershed was then calibrated with the "best" parameters of the Farnham, Noire River and Cowansville sub-watersheds.As the SWAT model is spatialized, streamflows coming from the Choiniere and Lac Brome reservoirs influence downstream flows.To account for reservoir management, the daily median was calculated based on streamflow observations from 1980-2010.These values were used as the reservoir output in the SWAT model.The same values were kept for reference and future climate conditions, underlying that reservoir operating rules will remain the same in a future climate.Regulation of these reservoirs depends on both the water supply and conservation of aquatic habitat and is relatively constant over the years.
A large variety of goodness-of-fit criteria have been proposed to evaluate the performance of hydrological models.These goodness-of-fit criteria can also be used as an objective function for the calibration of hydrological models.The NSE Q goodness-of-fit criterion is the widely used [34,35].Like many other goodness-of-fit criteria based on the mean square error, it puts greater weight on high flows, and therefore can adequately simulate snow accumulation and melt.A few studies have been carried out on the meaning and interpretation of goodness-of-fit criteria used to evaluate models in low-flow conditions [36] and concluded that the NSE Q is not suitable to evaluate model performance for low flows.Pushpalatha et al. [36] proposed other goodness-of-fit criteria suitable for evaluating low-flow simulations.However, they used a root mean square error calculated on root square streamflows as objective function, and the impact of the choice of the objective function on the performance of low-flow simulations was not examined.Based on these findings, six goodness-of-fit criteria were selected in this study and used as objective functions for the calibration of hydrological models.A seventh goodness-of-fit criterion based on the normalized 7Q2 (p7Q2) was also developed and tested, since this criterion is important for Yamaska River watershed managers.However, using only the p7Q2, calibration of both models tended towards a solution where there was not much accumulation of snow, resulting in low melting of the snow cover.In addition, models tended to simulate low flows during the late winter instead of summer.This can be explained by lower precipitation observed between February and March combined with low melting during this period.Thus, a combination of NSE Q and p7Q2 was tested to bring the model to better simulate the spring melt and still maintaining a good fit for low-flow.Using a threshold of 0.3 on the NSE Q and then optimizing the p7Q2, visually it was found that both calibrated models were able simulate reasonable accumulation and melting of snow.Table 2 summarizes these goodness-of-fit criteria.All goodness-of-fit criteria tend to 1 when the calibration is good.They were grouped into four categories.Group 1 is the NSE Q that focuses on high flows; group 2 focuses both on high flows and low flows (NSE sqrtQ and NSE lnQ ); group 3 applies to the summer period only (NSE Q-summer and NSE lnQ-summer ); and group 4 focuses on very low flows (NSE iQ and p7Q2).
Table 2. Description of the seven goodness-of-fit criteria used as objective functions for the model calibration and to evaluate model performance.

Goodness-of-Fit (Group)
Description Equation NSE calculated on streamflows between May and October 1 −

NSE calculated on log transformed streamflows between
May and October The 7-day low-flow value with a 2-year return period combined with a threshold on NSE calculated on streamflows Analysis of variance (ANOVA) was performed to quantify the contribution of each source of uncertainty [8,9,14] on the computed 7Q2 index as follows: where σ 2 7Q2 is the total variance on 7Q2 index, σ 2 CM , σ 2 HM , σ 2 PAR , σ 2 OF , σ 2 INT , and σ 2 ε are the variances linked to climate models, hydrological models, hydrological model parameters, objective functions, their interactions and the residual variance, respectively.Only first-order interactions (between two-factors) were considered.The ANOVA was performed for each sub-watershed in both reference and future periods.

Climate Change Projections
For the evaluation of the uncertainty related to climate models, 12 regional climate projections were used: 7 from the North America Regional Climate Change Assessment Program (NARCCAP) [37] and 5 from the Canadian Regional Climate Model version 4 (CRCM) [38] driven by the Canadian Global Climate Model, version 3 (CGCM3) [39].A bias correction on temperature and precipitation was performed to adjust the reference and future periods with the observations.This method assumes that the relative bias of climate models is the same for the future period and the reference period.For temperatures, a bias correction based on monthly averages was performed.The difference between the average monthly temperatures observed and the average monthly temperatures simulated for the reference period was applied to daily data for the reference and future periods.For precipitation, a bias correction was applied for both the monthly average precipitation occurrence (ratio of the number of days above a threshold of precipitation on the number of days in the month) and the average monthly rainfall intensity (average daily precipitation) using the Local Intensity method (LOCI).Details of this methodology can be found in of Schmidli et al. [40].
Figure 3 shows a statistically significant increase in atmospheric temperatures between the reference and future horizons.There is a mean difference of 2.5 • C for the minimal temperature and of 2.7 • C for the maximum temperature between the reference and future climate for the summer period (May to October) (Figure 3a,b).Precipitation is also significantly influenced by climate change.
Firstly, we can observe an increase in the variability of future climatic projections.Secondly, there is a decrease of the mean precipitation for the months of June and July of 5 mm (4.7%) and 9 mm (7.5%) respectively; and an increase for May, August, September and October, of 13 mm (13.7%), 8 mm (6.5%), 2 mm (1.8%, not significant) and 10 mm (10.3%) respectively (Figure 3c).Other studies support the increase in precipitation in late spring and early fall, and a decrease of precipitation during summer for southern Quebec (e.g., [41]).

Model Calibration and Validation
Both GR4J and SWAT models were calibrated using alternately each of the 7 goodness-of-fit criteria listed in Table 2 as objective functions, while the remaining six other criteria were computed to assess model performance.Figure 4 presents the relationship between the objective function values and other goodness-of-fit criteria for the best simulation (parameter set corresponding to the highest value of the objective function) for each of the four sub-watersheds for the model calibration and validation (i.e., 8 points per model for each graph in Figure 4).The closer the points are to the upper right corner of the graphs, the better the model.Values of NSE < 0 (which means that a model is a worse predictor than the average of the observations) were set to 0 for better visualization of the results.Data points located along the diagonal of a graph mean that the hydrological model performs equally well according to a given goodness-of-fit criterion as compared to the objective function used for the calibration/validation phases.Data points located above the diagonal mean that the model performs better with the objective functions in the calibration/validation phases as compared to the selected goodness-of fit criterion.
In general, the GR4J model (red) performed better than the SWAT model (blue), especially for the smallest sub-watershed, Cowansville, where low flows can be very low.The calibration using the NSEQ objective function shows values higher than 0.45 for both the calibration and validation period, but poorly represents the p7Q2 index and NSEiQ.Nevertheless, the seven objective functions for both models were less successful for both calibration and validation for this sub-watershed.The calibration using the p7Q2 converged to a value near one (perfect match between simulated and

Model Calibration and Validation
Both GR4J and SWAT models were calibrated using alternately each of the 7 goodness-of-fit criteria listed in Table 2 as objective functions, while the remaining six other criteria were computed to assess model performance.Figure 4 presents the relationship between the objective function values and other goodness-of-fit criteria for the best simulation (parameter set corresponding to the highest value of the objective function) for each of the four sub-watersheds for the model calibration and validation (i.e., 8 points per model for each graph in Figure 4).The closer the points are to the upper right corner of the graphs, the better the model.Values of NSE < 0 (which means that a model is a worse predictor than the average of the observations) were set to 0 for better visualization of the results.Data points located along the diagonal of a graph mean that the hydrological model performs equally well according to a given goodness-of-fit criterion as compared to the objective function used for the calibration/validation phases.Data points located above the diagonal mean that the model performs better with the objective functions in the calibration/validation phases as compared to the selected goodness-of fit criterion.Calibrating the GR4J model using the NSElnQ objective function gave better results for low flows than using the other objective functions for each of the four sub-watersheds simulated.This can be seen in Figure 4 (line 3) where the red points are closest to the upper right corner of the graphs.The objective functions NSElnQ-summer, NSEiQ and the one based on a combination of NSEQ and p7Q2 also gave good results for very low flows.
For the SWAT model, the calibration using NSEQ and NSEsqrtQ objective functions well represented these three goodness-of-fit criteria, but did not simulate well the very low flows, as shown by the poor goodness-of-fit criteria in the validation (group 4-NSEiQ and p7Q2).Only the NSElnQ and NSEiQ objective functions were able to successfully represent the NSEiQ goodness-of-fit.The p7Q2 value simulated using the NSEQ-summer objective function (line 4, column 7) was equal to zero for both the calibration and validation mode for the Noire River sub-watershed.This indicates that the simulated 7Q2 value was equal to 0 m 3 /s for both calibration and validation.The NSEQ-summer objective function also gave poor results for the NSE goodness-of-fit including a logarithm or inverse flow.It is likely that the model often simulated flows equal to 0 m 3 /s with the NSEQ-summer objective function for the Noire River sub-watershed.

Future Low-Flow and Uncertainty Analysis
The 7Q2 index was used to evaluate the impacts of climate change and their uncertainties for low flows.An ANOVA was performed according to Bosshard et al. [8] and Andor et al. [9] (see Equation ( 3)) to identify major sources of uncertainty on the computed 7Q2 index (Figure 5).In general, the GR4J model (red) performed better than the SWAT model (blue), especially for the smallest sub-watershed, Cowansville, where low flows can be very low.The calibration using the NSE Q objective function shows values higher than 0.45 for both the calibration and validation period, but poorly represents the p7Q2 index and NSE iQ .Nevertheless, the seven objective functions for both models were less successful for both calibration and validation for this sub-watershed.The calibration using the p7Q2 converged to a value near one (perfect match between simulated and observed 7Q2).However, the other goodness-of fit criteria were low.It is possible that the equifinality was greater for this objective function, and that the calibration did not converge to a global optimization.
Calibrating the GR4J model using the NSE lnQ objective function gave better results for low flows than using the other objective functions for each of the four sub-watersheds simulated.This can be seen in Figure 4 (line 3) where the red points are closest to the upper right corner of the graphs.The objective functions NSE lnQ-summer , NSE iQ and the one based on a combination of NSE Q and p7Q2 also gave good results for very low flows.
For the SWAT model, the calibration using NSE Q and NSE sqrtQ objective functions well represented these three goodness-of-fit criteria, but did not simulate well the very low flows, as shown by the poor goodness-of-fit criteria in the validation (group 4-NSE iQ and p7Q2).Only the NSE lnQ and NSE iQ objective functions were able to successfully represent the NSE iQ goodness-of-fit.The p7Q2 value simulated using the NSE Q-summer objective function (line 4, column 7) was equal to zero for both the calibration and validation mode for the Noire River sub-watershed.This indicates that the simulated 7Q2 value was equal to 0 m 3 /s for both calibration and validation.The NSE Q-summer objective function also gave poor results for the NSE goodness-of-fit including a logarithm or inverse flow.It is likely that the model often simulated flows equal to 0 m 3 /s with the NSE Q-summer objective function for the Noire River sub-watershed.

Future Low-Flow and Uncertainty Analysis
The 7Q2 index was used to evaluate the impacts of climate change and their uncertainties for low flows.An ANOVA was performed according to Bosshard et al. [8] and Andor et al. [9] (see Equation ( 3)) to identify major sources of uncertainty on the computed 7Q2 index (Figure 5).The major source of uncertainty in computing the 7Q2 index was the objective functions used to calibrate the hydrological model.The interaction-term between objective functions and hydrological models was also high.It shows the importance of the choice of the objective function for the calibration in low-flow analysis.The contribution of the hydrological model to the uncertainty increased in the future period for the Farnham and St-Hyacinthe sub-watersheds, while it decreased for Cowansville and Noire River sub-watersheds.It is interesting to note that both Farnham and St-Hyacinthe sub-watersheds are regulated, while Cowansville and Noire River sub-watersheds are natural.Note that SWAT can explicitly incorporate flow regulating structures such as dams, while GR4J cannot.This difference in model structure explains why hydrological model uncertainty increases in the future, as natural inflow variability is augmenting in a future climate, while hydraulic structures dampen this variability.An increase of the contribution of the hydrological models in the 2040-2070 period was also observed, along with a decrease of the contribution of the climate models in the future, and vice versa.As demonstrated by other authors [9,14], our results show that the relative importance of sources of uncertainties varies in time and with watersheds characteristics.
For each objective function and sub-watershed the ANOVA was applied to identify the contribution of the sources of uncertainty in the reference and future periods.The median, the 25%-75% interquartile range and the minimum and maximums values were also computed.Figures 6-9 present the results for the Cowansville (Figure 6), Farnham (Figure 7), Noire River (Figure 8) and St-Hyacinthe (Figure 9) sub-watersheds.The dotted line shows the observed value of the low-flow indices for the reference period.
The 7Q2 indices simulated in the reference period with climate models were compared against the 7Q2 indices observed over the same period.For the Cowansville sub-watershed with the GR4J model (Figure 6a), the 7Q2 indices were overestimated compared to the observed value of 7Q2 (dotted line) with the objective functions of groups 1, 2 and 3 for the reference period; and they were underestimated with the objective functions of group 4. The major source of uncertainty in computing the 7Q2 index was the objective functions used to calibrate the hydrological model.The interaction-term between objective functions and hydrological models was also high.It shows the importance of the choice of the objective function for the calibration in low-flow analysis.The contribution of the hydrological model to the uncertainty increased in the future period for the Farnham and St-Hyacinthe sub-watersheds, while it decreased for Cowansville and Noire River sub-watersheds.It is interesting to note that both Farnham and St-Hyacinthe sub-watersheds are regulated, while Cowansville and Noire River sub-watersheds are natural.Note that SWAT can explicitly incorporate flow regulating structures such as dams, while GR4J cannot.This difference in model structure explains why hydrological model uncertainty increases in the future, as natural inflow variability is augmenting in a future climate, while hydraulic structures dampen this variability.An increase of the contribution of the hydrological models in the 2040-2070 period was also observed, along with a decrease of the contribution of the climate models in the future, and vice versa.As demonstrated by other authors [9,14], our results show that the relative importance of sources of uncertainties varies in time and with watersheds characteristics.
For each objective function and sub-watershed the ANOVA was applied to identify the contribution of the sources of uncertainty in the reference and future periods.The median, the 25%-75% interquartile range and the minimum and maximums values were also computed.Figures 6-9 present the results for the Cowansville (Figure 6), Farnham (Figure 7), Noire River (Figure 8) and St-Hyacinthe (Figure 9) sub-watersheds.The dotted line shows the observed value of the low-flow indices for the reference period.For the simulations with the SWAT model, each objective function overestimated the 7Q2 values in the reference period when compared to the observed values (dotted line) for this period (Figure 6b).The variance (Figure 6c) was higher for the NSElnQ-summer objective function, where the contribution of the hydrological model to overall uncertainty is large.This is because SWAT produced much larger 7Q2 values compared to GR4J when calibrated with this objective function.In general, a lower uncertainty for the 7Q2 index was observed for the future period than for the reference period.For the NSEQ objective function, climate model is the main source of uncertainty, while for objective functions more related to low-flow, the main contribution to the total uncertainty is the hydrological model.
At the Farnham sub-watershed, using the GR4J model, all objective functions overestimated the 7Q2 using climatic projections of the reference period compared to the observed value (dotted line) for the same period (Figure 7a).When the SWAT model was used (Figure 7b), the NSEQ, NSElnQ, NSEiQ, NSElnQ-summer objective functions overestimated the 7Q2 value observed in the same period (while the NSEsqrtQ, NSEQ-summer, p7Q2 functions underestimated this 7Q2 value).The main source of uncertainty for the NSEQ, NSElnQ, NSElnQ-summer, and P7Q2 objective functions was the climate model, while the hydrological model was the main source of uncertainty for other objective functions.For the NSEQ, NSEsqrtQ, NSElnQ, and NSEQ-summer objective functions, the total variance in the future period was lower than the total variance in the reference period, while the opposite was observed when the NSElnQ-summer, NSEiQ, and P7Q2 objective functions were used to calibrate the models.This difference in the total variance came from the variance related to the hydrological model.The interaction term between climate models and hydrological models was not negligible, but was similar between the reference and future periods.
For the Noire River sub-watershed, the GR4J model using the climatic projections for the reference period overestimated the observed values of 7Q2 for all objective functions (Figure 8a).Similarly, the SWAT model overestimated the 7Q2 observed values using the climatic projections in the reference period for all objective functions, with the NSEQ objective function being worse than the others (Figure 7b).Model calibration with the NSEQ-summer objective function was problematic since all the simulated 7Q2 values were equal to zero.Calibrating the model using this objective function simulated very low contribution of groundwater to streamflow (see section 4.3.2),resulting of values of streamflow equal to zero.The main source of uncertainty (Figure 8c) was the hydrological model for all objective functions except for the NSElnQ objective function, where the main contribution was the climate model.This objective function also shown the lower total variance compared to the other objective functions.Furthermore, total variance was lower for the future period than for the reference period.
The SWAT model behaved similar to the GR4J model.Each objective function, except NSEQsummer, overestimated the 7Q2 values in the reference period.Note that the simulations for the St-Hyacinthe sub-watershed were influenced by the simulations of the upstream sub-watersheds.As The 7Q2 indices simulated in the reference period with climate models were compared against the 7Q2 indices observed over the same period.For the Cowansville sub-watershed with the GR4J model (Figure 6a), the 7Q2 indices were overestimated compared to the observed value of 7Q2 (dotted line) with the objective functions of groups 1, 2 and 3 for the reference period; and they were underestimated with the objective functions of group 4.
For the simulations with the SWAT model, each objective function overestimated the 7Q2 values in the reference period when compared to the observed values (dotted line) for this period (Figure 6b).The variance (Figure 6c) was higher for the NSE lnQ-summer objective function, where the contribution of the hydrological model to overall uncertainty is large.This is because SWAT produced much larger 7Q2 values compared to GR4J when calibrated with this objective function.In general, a lower uncertainty for the 7Q2 index was observed for the future period than for the reference period.For the NSE Q objective function, climate model is the main source of uncertainty, while for objective functions more related to low-flow, the main contribution to the total uncertainty is the hydrological model.
At the Farnham sub-watershed, using the GR4J model, all objective functions overestimated the 7Q2 using climatic projections of the reference period compared to the observed value (dotted line) for the same period (Figure 7a).When the SWAT model was used (Figure 7b), the NSE Q , NSE lnQ , NSE iQ , NSE lnQ-summer objective functions overestimated the 7Q2 value observed in the same period (while the NSE sqrtQ , NSE Q-summer , p7Q2 functions underestimated this 7Q2 value).The main source of uncertainty for the NSE Q , NSE lnQ , NSE lnQ-summer , and P7Q2 objective functions was the climate model, while the hydrological model was the main source of uncertainty for other objective functions.For the NSE Q , NSE sqrtQ , NSE lnQ , and NSE Q-summer objective functions, the total variance in the future period was lower than the total variance in the reference period, while the opposite was observed when the NSE lnQ-summer , NSE iQ , and P7Q2 objective functions were used to calibrate the models.This difference in the total variance came from the variance related to the hydrological model.The interaction term between climate models and hydrological models was not negligible, but was similar between the reference and future periods.
For the Noire River sub-watershed, the GR4J model using the climatic projections for the reference period overestimated the observed values of 7Q2 for all objective functions (Figure 8a).Similarly, the SWAT model overestimated the 7Q2 observed values using the climatic projections in the reference period for all objective functions, with the NSE Q objective function being worse than the others (Figure 7b).Model calibration with the NSE Q-summer objective function was problematic since all the simulated 7Q2 values were equal to zero.Calibrating the model using this objective function simulated very low contribution of groundwater to streamflow (see Section 4.3.2),resulting of values of streamflow equal to zero.The main source of uncertainty (Figure 8c) was the hydrological model for all objective functions except for the NSE lnQ objective function, where the main contribution was the climate model.This objective function also shown the lower total variance compared to the other objective functions.Furthermore, total variance was lower for the future period than for the reference period.
The SWAT model behaved similar to the GR4J model.Each objective function, except NSE Qsummer , overestimated the 7Q2 values in the reference period.Note that the simulations for the St-Hyacinthe sub-watershed were influenced by the simulations of the upstream sub-watersheds.As the NSE Qsummer objective function often simulated streamflows equal to zero for the Noire River sub-watershed, the streamflows simulated for the St-Hyacinthe sub-watershed were lower with this objective function.The NSE Q is the objective function that overestimated the 7Q2 values the most (Figure 9b).The total variance was lower for the future period than for the reference period for all objective functions except the NSE Q .The variance of NSE sqrtQ and NSE lnQ objective functions were the lowest.For all objective functions, except for these two objective functions, the main source of uncertainty was the hydrological model.
For all sub-watersheds, the NSE lnQ shown the lowest total variance and the lowest variance related to hydrological models.In general, the SWAT model shown higher 25%-75% interquartile range than the GR4J model.
Comparison of the 7Q2 index for the reference and future periods was used to evaluate whether low flows will worsen with climate change.Relative changes in the median between the reference and future periods are indicated under the x-axis for each sub-watershed (Figures 6-8).A positive change indicates that future 7Q2 is higher than reference 7Q2.Significant changes (p-value < 0.05) are shown in bold.
For the GR4J model, the changes shown in Figures 6a, 7a, 8a and 9a are all negative and statistically significant, meaning that there is a worsening, or decrease, of low flows for all objective functions and sub-watersheds assessed in the study.Current observed values of 7Q2 indices are also indicated in these figures by a dotted line.The changes vary from −11.9% to −28.8% depending on the objective functions and the sub-watersheds considered in the analysis.The Cowansville sub-watershed showed the smallest range of change depending of objective functions (from −20.6% to −24.9%), while the Noire River sub-watershed shown the largest range (from −11.9% to −28.1%).The Cowansville sub-watershed is the most forested, while the Noire River sub-watershed is the most agricultural sub-watershed in this study.The simulation using the p7Q2 objective function to calibrate the GR4J model shown the smallest change for 7Q2 between the reference and future periods, depending on the sub-watershed analyzed.
The changes in 7Q2 between current and future climates simulated with the SWAT model were not always significant (Figures 6b, 7b, 8b and 9b).Changes that were significant were almost always negative, meaning a worsening in low-flow conditions.The changes simulated for the Cowansville sub-watershed were mostly significant and were all negative.The changes in 7Q2 simulated with the SWAT model for this sub-watershed were much larger than the ones simulated with the GR4J model.The changes greatly varied according to the objective function used, from −19.8% to −59.2%.
For the Farnham sub-watershed, there were only a few significant changes in the simulated 7Q2, which were all negative.Some objective functions gave positive differences, although not statistically significant.Reservoir flow regulation for this sub-watershed (Choinière, Waterloo, and Brome Lake, see Figure 1) simulated with the SWAT model explains these results as these reservoirs control downstream flows by reducing high flows and increasing low flows.Therefore, the presence of reservoirs should help to alleviate the negative impacts of climate change on low flows.Note that the GR4J model resulted in significant negative changes for the Farnham sub-watershed.This is because GR4J does not explicitly simulate reservoir storage and routing.This illustrates the importance of a hydrological model capable of explicitly handling storage structures for low-flow impact studies in regulated watersheds.For the Noire River and St-Hyacinthe sub-watersheds, the NSE Q-summer objective function used to calibrate SWAT did not give conclusive results.The NSE Q , NSE sqrtQ , NSE lnQ and NSE lnQ-summer objective functions all resulted in a decrease of the 7Q2 value.However, the change in the 7Q2 values between the reference and future periods was not significant.This may be caused by the difficulty in simulating very low flows with SWAT for those sub-watersheds (zero flow values were obtained), which possibly related to the high number of parameters in the model, as compared to GR4J.

Hydrological Processes and State Variable Analysis
To better understand the results of the uncertainty analysis, different hydrological processes and state variables for the GR4J and SWAT models were analyzed according to the objective functions used to calibrate the models.

The GR4J Model
For the GR4J model, the hydrological processes analyzed were actual evapotranspiration (ET) and percolation (PERC).The model state variables analyzed included the soil moisture reservoir (S) and the routing reservoir (R).For the GR4J model, the dynamics of hydrological processes and temporal characteristics of state variables were found to be similar for all four sub-watersheds.
Overall, GR4J calibrated using NSE Q simulated less ET than using other objective functions (Figure 10).Potential evapotranspiration was simulated with the Hamon model and was identical for all objective functions used to calibrate the model.Actual evapotranspiration depends of the X5 parameter and the on the amount of water in the soil, represented by the soil moisture reservoir S, that depends of the X4 parameter (Section 3.1).These parameters changed depending of the objective function used to calibrate the model.The simulations obtained with GR4J calibrated with all objective functions, except p7Q2, showed a significant increase of ET during the summer-fall period (May-October).For the p7Q2, ET remained unaffected during that period between reference and future horizons.On a seasonal basis, the simulations showed an increase in ET between May and June in the future period compared to the reference period, while simulated ET remained similar or decreased slightly from July to October depending on the objective function selected for model calibration.A decrease in ET between June and October is explained by a lack of water in the soil moisture reservoir to satisfy potential evapotranspiration.The increase in ET between May and June thus can lead to a decrease in the amount of water present in the soil during the summer.The sum of actual ET from May to October significantly decreased between 8% and 11% depending on the sub-watershed and objective functions, except for NSE iQ and p7Q2, while the difference was not significant.For all objective functions, uncertainties in ET were larger in the future than in the reference period.Larger uncertainty in temperatures in the future period explains this phenomenon (see Equation ( 1)).Since the evapotranspiration uncertainty increased in the future period, it cannot explain the decreasing uncertainty for low-flow indices obtained for the future period (Figures 6a, 7a, 8a and 9a).
Total percolation simulated over the entire year varied between 40 mm to 85 mm depending on the sub-watershed and objective function.The overall evolution of percolation during the year was similar for all calibrations and sub-watersheds, showing lower percolation from May to October.The sums of percolation for the summer (May to October) are presented in Figure 11.The percolation is dependent on the X4 parameter, the maximum capacity of the soil reservoir, which relates to soil properties and land use/land cover characteristics.For all objective functions, uncertainty in percolation decreased between May and October in the future period compared to the reference period (Figure 11).The low percolation simulated is explained by a shortage of water in the soil moisture reservoir.Thus, the overall decrease in the uncertainty for the 7Q2 index in the future period can be explained by a decrease in the uncertainty for percolation during the summer period.The soil reservoir (S) was analyzed based on the percentage to which the reservoir is filled.Generally, it was found that the level of the reservoir S was high during the period from November to March, reaching 86% to 90% capacity, depending on the objective function used to calibrate the model.It gradually decreased until July when it reached a minimum and then gradually increased.The simulation with the NSE Q objective function had the smallest variation between the maximum and the minimum level of reservoir S, while simulations with NSE iQ and p7Q2 objective functions had the largest variation (Table 3).Low levels of reservoir S observed with the NSE iQ and p7Q2 objective functions mean that there was a small amount of water in the soil during the summer period.This small amount of water in the soil limited the ET and percolation between mid-June and late August.In addition, the difference between the minimum and maximum level of reservoir S in the future period was significantly lower than in the reference period.The water content in the soil was significantly lower in the 2040-2070 period.The routing reservoir (R) was also analyzed based on the percentage to which the reservoir was filled.The level of the reservoir R was typically high during the period from November to March, reaching 60% to 80% capacity depending on the sub-watershed and objective functions.It gradually decreased until July when it reached a minimum followed by a gradual rise (Table 4).The difference between the maximum and minimum level was significantly lower for the NSE Q objective function than for other objective functions.In the future period, the maximum level did not show a significant difference from the current period, whereas the minimum levels during the summer decreased significantly.In addition, in the reference period, the minimum level was reached during the winter, whereas for the future period, the minimum level was reached in the summer.

The SWAT Model
For the SWAT model, the hydrological processes analyzed were actual evapotranspiration in mm (ET), and groundwater contribution to streamflow in mm (GW).The soil water content in mm (SW), a model state variable, was also analyzed.For the SWAT model, the response of the hydrological processes differed for the four sub-watersheds and the objective function used to calibrate the model.This can be explained by the higher number of parameters in the SWAT model leading a different combination of parameters to a similar final result (streamflow).Although, intermediate results, such as hydrological processes, changed depending on the combination of parameters.As another possibility, differences in land cover/land use, drainage area, and slope between sub-watersheds may partly explain the differences in simulation results (Table 1).
The hydrological response of the Farnham sub-watershed is influenced by the regulation of the upstream reservoirs simulated by the SWAT model.In addition, the hydrological response of the St-Hyacinthe sub-watershed is highly influenced by upstream sub-watersheds.Thus, the hydrological processes and state variables of these two sub-watersheds were not analyzed here and only the Cowansville and Noire River sub-watersheds were analyzed.
Total actual evapotranspiration was computed for the summer/fall season (May to October).Potential evapotranspiration was simulated using the Penman-Monteith method and was the same for all objective functions used to calibrate SWAT.Actual evapotranspiration depends on the amount of water in the soil and a specific parameter (ESCO).The amount of water in the soil simulated by SWAT depends of parameters related to the characteristics of the soil.For the Cowansville sub-watershed, ET simulated with the NSE lnQ-summer and NSE iQ objective functions was statistically significantly higher than ET simulated with the NSE Q objective function (Figure 12a) for both current and future periods.For the Noire River sub-watershed (Figure 12b), ET simulated with the NSE lnQ objective function was significantly lower than with the other objective functions, again for current and future climates.In the future period, a statistically significant increase of median values of ET between 13% and 14% was simulated for the Cowansville sub-watershed and between 11% and 12% for the Noire River sub-watershed.The objective function used to calibrate the model therefore had little impact on the magnitude of the change of ET attributable to climate change.The total GW for the summer/fall period (May to October) was also analyzed (Figure 13).For the Cowansville sub-watershed (Figure 13a), a very low contribution from groundwater was observed for the simulation using the NSEiQ and p7Q2 objective functions.A significant difference was observed between values of GW for these two objective functions and all other objective functions.The simulation using the NSEQ objective function showed the highest contribution from groundwater to streamflow.A significant difference was observed between the value of GW for the simulation using NSEQ and all other objective functions except NSEQ_summer.For the Noire River sub-watershed a very low contribution of groundwater was simulated with the NSEQ-summer objective function.The low-flow indices simulated with the NSEQ-summer objective function were therefore equal to 0 (Figure 8b).For the future period, a significant decrease of the groundwater contribution to streamflow was observed with the NSEQ, NSEsqrtQ and NSElnQ objective function for the Cowansville sub-watershed.No significant difference was observed between the reference and future periods for the Noire River sub-watershed, as the corresponding boxplots overlap, see Figure 13b.The total GW for the summer/fall period (May to October) was also analyzed (Figure 13).For the Cowansville sub-watershed (Figure 13a), a very low contribution from groundwater was observed for the simulation using the NSE iQ and p7Q2 objective functions.A significant difference was observed between values of GW for these two objective functions and all other objective functions.The simulation using the NSE Q objective function showed the highest contribution from groundwater to streamflow.A significant difference was observed between the value of GW for the simulation using NSE Q and all other objective functions except NSE Q_summer .For the Noire River sub-watershed a very low contribution of groundwater was simulated with the NSE Q-summer objective function.The low-flow indices simulated with the NSE Q-summer objective function were therefore equal to 0 (Figure 8b).For the future period, a significant decrease of the groundwater contribution to streamflow was observed with the NSE Q , NSE sqrtQ and NSE lnQ objective function for the Cowansville sub-watershed.No significant difference was observed between the reference and future periods for the Noire River sub-watershed, as the corresponding boxplots overlap, see Figure 13b.The total GW for the summer/fall period (May to October) was also analyzed (Figure 13).For the Cowansville sub-watershed (Figure 13a), a very low contribution from groundwater was observed for the simulation using the NSEiQ and p7Q2 objective functions.A significant difference was observed between values of GW for these two objective functions and all other objective functions.The simulation using the NSEQ objective function showed the highest contribution from groundwater to streamflow.A significant difference was observed between the value of GW for the simulation using NSEQ and all other objective functions except NSEQ_summer.For the Noire River sub-watershed a very low contribution of groundwater was simulated with the NSEQ-summer objective function.The low-flow indices simulated with the NSEQ-summer objective function were therefore equal to 0 (Figure 8b).For the future period, a significant decrease of the groundwater contribution to streamflow was observed with the NSEQ, NSEsqrtQ and NSElnQ objective function for the Cowansville sub-watershed.No significant difference was observed between the reference and future periods for the Noire River sub-watershed, as the corresponding boxplots overlap, see Figure 13b.The median value of SW for the summer/fall period was analyzed (Figure 14).SW for the Cowansville sub-watershed was significantly higher than in the Noire River sub-watershed by approximately 90-100 mm for both current and future climates.These differences are attributed to differences in land cover/land use of these watersheds, the Cowansville being predominantly forested, while the Noire is mostly agricultural.For the future period, a significant decrease of SW between 1% and 2% was simulated for the Cowansville sub-watershed and between 2% and 3% for the Noire River sub-watershed for all objective functions considered.Differences attributed the choice of the objective function for model calibration were also noticed.For the Cowansville sub-watershed (Figure 14a), significantly lower values of SW were simulated with the NSE iQ , NSE Q-summer and p7Q2 objective functions.For the Noire River sub-watershed (Figure 14b), the NSE Q-summer this time showed a significantly higher value of SW.The lower ET simulated with this objective function resulted in a higher soil water content.The median value of SW for the summer/fall period was analyzed (Figure 14).SW for the Cowansville sub-watershed was significantly higher than in the Noire River sub-watershed by approximately 90-100 mm for both current and future climates.These differences are attributed to differences in land cover/land use of these watersheds, the Cowansville being predominantly forested, while the Noire is mostly agricultural.For the future period, a significant decrease of SW between 1% and 2% was simulated for the Cowansville sub-watershed and between 2% and 3% for the Noire River sub-watershed for all objective functions considered.Differences attributed the choice of the objective function for model calibration were also noticed.For the Cowansville sub-watershed (Figure 14a), significantly lower values of SW were simulated with the NSEiQ, NSEQ-summer and p7Q2 objective functions.For the Noire River sub-watershed (Figure 14b), the NSEQ-summer this time showed a significantly higher value of SW.The lower ET simulated with this objective function resulted in a higher soil water content.To summarize, analysis of the SWAT results did not clearly demonstrate a causal relationship between the objective function related to low flows used to calibrate the model and hydrological processes influencing low-flow values.Indeed, the different processes analyzed did not follow a particular trend.For example, groundwater flow in the Cowansville sub-watershed was very low for the NSEiQ and p7Q2 objective functions (see Figure 14a) and corresponding 7Q2 (see Figure 6b) were among the lowest for these objective functions.However, a similar consistency was not observed in the Noire River sub-watershed, where objective functions producing the lowest groundwater flow (NSElnQ and NSEln-Q-summer) did not produce values of 7Q2which were among the lowest (see Figure 8b).Thus, the low flows produced by the SWAT model did not result from the same hydrological processes (e.g., groundwater flow), but rather from different hydrological processes (e.g., groundwater, interflow, evapotranspiration, etc.), acting individually or in combination and which differed from one watershed to another.Only further analysis of the model results, together with observations, could help to better understand the behavior of the model as a function of different objective functions and also of land use/land cover types.Finally, the question of whether modeling of these processes was faithful to what is actually observed in the field remains unresolved due to the lack of data to make such an assessment.For example, although the overall seasonal variation of ET, SW and GW is acceptable, measurements of these processes and state variables will allow quantitatively evaluating simulation results.It is also possible that equifinality, causing several sets of parameters to produce a similar response from the model, may be the cause of the lack of consistency observed between the modeled hydrological processes and low-flow indices.To summarize, analysis of the SWAT results did not clearly demonstrate a causal relationship between the objective function related to low flows used to calibrate the model and hydrological processes influencing low-flow values.Indeed, the different processes analyzed did not follow a particular trend.For example, groundwater flow in the Cowansville sub-watershed was very low for the NSE iQ and p7Q2 objective functions (see Figure 14a) and corresponding 7Q2 (see Figure 6b) were among the lowest for these objective functions.However, a similar consistency was not observed in the Noire River sub-watershed, where objective functions producing the lowest groundwater flow (NSE lnQ and NSE ln-Q-summer ) did not produce values of 7Q2which were among the lowest (see Figure 8b).Thus, the low flows produced by the SWAT model did not result from the same hydrological processes (e.g., groundwater flow), but rather from different hydrological processes (e.g., groundwater, interflow, evapotranspiration, etc.), acting individually or in combination and which differed from one watershed to another.Only further analysis of the model results, together with observations, could help to better understand the behavior of the model as a function of different objective functions and also of land use/land cover types.Finally, the question of whether modeling of these processes was faithful to what is actually observed in the field remains unresolved due to the lack of data to make such an assessment.For example, although the overall seasonal variation of ET, SW and GW is acceptable, measurements of these processes and state variables will allow quantitatively evaluating simulation results.It is also possible that equifinality, causing several sets of parameters to produce a similar response from the model, may be the cause of the lack of consistency observed between the modeled hydrological processes and low-flow indices.

Discussion and Conclusions
This study investigated the impact of the objective function on the uncertainty of hydrological model results for simulating low flows in a climate change context.Two models were evaluated: the GR4J lumped conceptual model and the SWAT spatially distributed model.A PARASOL-based method was used to evaluate uncertainties related to hydrological model parameters.An ANOVA method was used to quantify the contribution of climate models, hydrological models, hydrological model parameters, objective functions and their interaction to the total uncertainty of an index related to low-flow (7Q2).
Our results shown that the objective function is an important source of uncertainty for low-flow simulations.This finding is partially supported by Parajka et al. [14] who concluded that uncertainty from the objective function is not negligible, but not the main contribution to total uncertainty.However, they used 11 different combinations of a weighted average between two criteria, namely NSE Q and NSE lnQ , along with a single hydrological model.Our study used a larger spectrum of objective functions related to low-flow with two hydrological models of different complexity.Our results also showed that the interaction-term between the hydrological models and the objective functions is important.The contribution of hydrological model parameters is the lowest contribution to the total uncertainty of low-flow projections.When the ANOVA was applied for each objective function, the relative contribution of climate models and hydrological models was shown to vary with watershed characteristics and time horizon (reference and future periods).This is consistent with previous studies [9,14].
The model calibrated with the NSE Q objective function, one of the seven objective functions tested in this study, and which is commonly used for model calibration in climate changes studies, was found to poorly simulate low flows.Besides, many studies assessing the impact of climate change on the hydrological response highlight that calibrated models have difficulty in simulating low flows compared to observations (e.g., [11,13]).Furthermore, total variance on the 7Q2 index was large using this objective function.Studies on low flows using NSE Q may be producing high uncertainties and therefore should be avoided.Using the logarithmic transformed NSE (NSE lnQ ) objective function, both models succeeded in adequately representing the other's goodness-of-fit related to low flows.This finding confirms results reported by other studies [14,36].In addition, using the NSE lnQ objective function, the total variance on 7Q2 index was the lowest for all sub-watersheds.
As observed by Najifi et al. [10] and Chen et al. [7], the structure of the hydrological model can influence estimates of climate change impacts on low flows.Using 12 regional climate projections from different climate models, the GR4J model calibrated with objective functions related to low flows better reproduced the observed 7Q2 index for the reference (i.e.,  period.The SWAT model hardly reproduced the observed indices for the reference period.This problem was also noticed by Dams et al. [13] using SWAT.Moreover, the 25%-75% interquartile range on 7Q2 index was found to be higher for the SWAT model than for the GR4J model.SWAT includes a large number of parameters compared to GR4J, which may amplify the effect of the model equifinality, and therefore be responsible for the larger parameter uncertainty obtained with SWAT.For the GR4J model, the 25%-75% interquartile range on 7Q2 index was lowest with the NSE iQ and p7Q2 objective functions, which are more related to low flows.
In general, the total variance was lower for the future period than the reference period. .For high flows, precipitation is an important factor in calculating streamflow.Increasing the uncertainty of precipitation caused by the climate model structure led to an increase in the uncertainty for high flows for the future period [8].In the dry period, the hydrological model structure can plays a major role and an increase in the uncertainty of precipitation did not lead necessarily to increasing uncertainty for low flows.The analysis of hydrological processes can help understanding why the uncertainty of low-flow indices for the future period was lower than for the reference period.There is no upper limit to streamflow during high flows, while streamflows are physically bounded to zero.In the future period, more streamflow values are simulated close to zero, which reduces the uncertainty on low-flow indices, such as the 7Q2.
Figures 9 and 11 show that uncertainties for actual evapotranspiration increased in the future horizon compared to the reference period for both models, even when simulated ET is different between the GR4J (from 150 to 300 mm) and the SWAT model (from 300 to 400 mm).The difference may be due to the model used to calculate potential evapotranspiration (ETP) (Hamon for GR4J, Penman-Monteith for SWAT), since comparing ETP produced by these two models shown differences.However, ETP simulated by the Hamon model for each sub-watershed was higher than ETP simulated with Penman-Monteith model, which cannot explain higher ET simulated by the SWAT model.Therefore, the difference must come from the model structure used to convert ETP to ET, which is very different between GR4J and SWAT.In addition, soil moisture is conceptualized differently in each model, which may contribute to the difference in simulated ET values.The increase of uncertainty on ET is difficult to tie to the decrease of uncertainty on the 7Q2 index, which, as mentioned above, is related to longer periods of very low flow values in a future climate.
For the GR4J model, uncertainties in percolation decrease since the available water in the soil also decreases.This can explain the lower uncertainties observed for the 7Q2 index in the future period compared to the reference period, despite the increase in the uncertainty for precipitation.
Hydrological processes in the SWAT model are more complex and the objective function used for the calibration may change the dominant hydrological processes leading to low-flow situations.Different land cover/land use between sub-watersheds, which are explicitly represented in SWAT, but not in GR4J, may also partly explain the differences in hydrological processes driving the production of low-flow.The nested method used for the calibration can also introduce uncertainty, since model calibration of upstream sub-watersheds can influence the calibration of downstream sub-watersheds.These findings highlight the need to further investigate hydrological processes involved in low-flow simulations, and hydrological models should be used with care for low-flow studies, since they show some inconsistent results.
For the Cowansville sub-watershed, there is a clear signal from the two hydrological models that low flows will be worsening in the future.It is difficult to make conclusions about the impact of climate change on lows flows in the Farnham sub-watershed since streamflows are regulated.In particular, while a significant decrease in low flows in the future is expected with the GR4J model, which does not handle reservoir storage, conflicting results are obtained depending on the objective function used to calibrate the SWAT model.Additional studies should be conducted to determine whether water released from the reservoir can be maintained in the future without affecting other usages and violating reservoir constraints.For the Noire River and St-Hyacinthe sub-watersheds, results differ depending on the hydrological model used to conduct the impact study.For example, a decrease in low flows was obtained with GR4J as shown by the decrease in both 7Q2 and 7Q10 indices for all objective functions used to calibrate the model.However, these results neglect the effect of reservoir storage in the Farnham sub-watershed located upstream, and are therefore considered suspect since GR4J cannot simulate reservoir storage.On the other hand, the SWAT model, which explicitly considers storage regulation, generated no significant changes in the 7Q2 index, while at the same time a significant worsening of low-flow conditions was obtained through a lowering of the 7Q10 index in the St-Hyacinthe watershed.Again, possible equifinality problems with SWAT may explain these contradictory results.
Great care therefore must be taken in selecting an objective function to calibrate a hydrological model which is well adapted to simulate low flows.This is particularly striking in climate change studies for watersheds which expect a worsening of low-flow conditions in the future.Moreover, this study highlights the importance of model structure in simulating low flows in climate change impact assessments.While models must explicitly consider storage retention for regulated watersheds, increasing model complexity must be addressed with care as low flows are the result of hydrological processes such as infiltration, evapotranspiration and percolation, which interact to produce subtle changes in low-flow regimes in a way that is very dependent on watershed physiographic characteristics.More than ever, climate change impact studies on hydrological regimes require conducting hydrological model uncertainty assessments.

Figure 2 .
Figure 2. Schematic diagrams of the structure of the (a) GR4J model; and (b) SWAT model.

Figure 2 .
Figure 2. Schematic diagrams of the structure of the (a) GR4J model; and (b) SWAT model.

Figure 3 .
Figure 3. Uncertainties related to the climate model for the (a) monthly minimal temperature; (b) monthly maximal temperature; and (c) precipitation for the Yamaska River watershed.Boxes represent the 25-75 percentile values.Whiskers correspond to the 5-95 percentile values.

Figure 3 .
Figure 3. Uncertainties related to the climate model for the (a) monthly minimal temperature; (b) monthly maximal temperature; and (c) precipitation for the Yamaska River watershed.Boxes represent the 25-75 percentile values.Whiskers correspond to the 5-95 percentile values.

Figure 4 .
Figure 4. Relationship between objective function values and goodness-of-fit criteria for the best simulation.

Figure 4 .
Figure 4. Relationship between objective function values and goodness-of-fit criteria for the best simulation.

Figure 5 .
Figure 5. ANOVA partitioning among the sources of uncertainty computed on the 7Q2 index for each sub-watershed in reference and future periods.

Figure 5 .
Figure 5. ANOVA partitioning among the sources of uncertainty computed on the 7Q2 index for each sub-watershed in reference and future periods.

Figure 6 .
Figure 6.Uncertainties for the 7Q2 indices for the Cowansville sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective function.

Figure 7 .
Figure 7. Uncertainties for the 7Q2 indices for the Farnham sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.

Figure 8 .
Figure 8. Uncertainties for the 7Q2 indices for the Noire River sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.For the St-Hyacinthe sub-watershed, the GR4J model overestimated the 7Q2 values for the reference period for all objective functions (Figure 9a).

Figure 6 .Figure 6 .
Figure 6.Uncertainties for the 7Q2 indices for the Cowansville sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective function.

Figure 7 .
Figure 7. Uncertainties for the 7Q2 indices for the Farnham sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.

Figure 8 .
Figure 8. Uncertainties for the 7Q2 indices for the Noire River sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.For the St-Hyacinthe sub-watershed, the GR4J model overestimated the 7Q2 values for the reference period for all objective functions (Figure 9a).

Figure 7 .Figure 6 .
Figure 7. Uncertainties for the 7Q2 indices for the Farnham sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.

Figure 7 .
Figure 7. Uncertainties for the 7Q2 indices for the Farnham sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.

Figure 8 .
Figure 8. Uncertainties for the 7Q2 indices for the Noire River sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.For the St-Hyacinthe sub-watershed, the GR4J model overestimated the 7Q2 values for the reference period for all objective functions (Figure 9a).

Figure 8 .
Figure 8. Uncertainties for the 7Q2 indices for the Noire River sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.For the St-Hyacinthe sub-watershed, the GR4J model overestimated the 7Q2 values for the reference period for all objective functions (Figure9a).

Figure 9 .
Figure 9. Uncertainties for the 7Q2 indices for the St-Hyacinthe sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.

Figure 9 .
Figure 9. Uncertainties for the 7Q2 indices for the St-Hyacinthe sub-watershed related to the climate model structure for the (a) GR4J; and (b) SWAT model; and (c) comparison of 7Q2 index variance for each source of uncertainty in reference (Ref) and future (Fut) periods for each objective functions.

Figure 10 .
Figure 10.Sum of actual evapotranspiration simulated with the GR4J model for the summer (May to October) for the (a) Cowansville; (b) Farnham; (c) Noire River; and (d) St-Hyacinthe sub-watersheds.

Figure 11 .
Figure 11.Sum of percolation simulated with the GR4J model for the summer (May to October) for the (a) Cowansville; (b) Farnham; (c) Noire River; and (d) St-Hyacinthe sub-watersheds.

Figure 10 . 23 Figure 10 .
Figure 10.Sum of actual evapotranspiration simulated with the GR4J model for the summer (May to October) for the (a) Cowansville; (b) Farnham; (c) Noire River; and (d) St-Hyacinthe sub-watersheds.

Figure 11 .
Figure 11.Sum of percolation simulated with the GR4J model for the summer (May to October) for the (a) Cowansville; (b) Farnham; (c) Noire River; and (d) St-Hyacinthe sub-watersheds.

Figure 11 .
Figure 11.Sum of percolation simulated with the GR4J model for the summer (May to October) for the (a) Cowansville; (b) Farnham; (c) Noire River; and (d) St-Hyacinthe sub-watersheds.

Figure 12 .
Figure 12.Sum of actual evapotranspiration simulated with the SWAT model for the summer (May to October) for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Figure 13 .
Figure 13.Sum of groundwater contribution to streamflow during the summer simulated with the SWAT model for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Figure 12 .
Figure 12.Sum of actual evapotranspiration simulated with the SWAT model for the summer (May to October) for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Climate 2017, 5 , 19 18 of 23 Figure 12 .
Figure 12.Sum of actual evapotranspiration simulated with the SWAT model for the summer (May to October) for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Figure 13 .
Figure 13.Sum of groundwater contribution to streamflow during the summer simulated with the SWAT model for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Figure 13 .
Figure 13.Sum of groundwater contribution to streamflow during the summer simulated with the SWAT model for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Figure 14 .
Figure 14.Median of soil water content during the summer simulated with the SWAT model for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Figure 14 .
Figure 14.Median of soil water content during the summer simulated with the SWAT model for the (a) Cowansville sub-watershed; and (b) Noire River sub-watershed.

Table 1 .
Physiographic characteristics of the study sub-watersheds.

Table 3 .
Variation of the levels of the soil reservoir S for the GR4J model.

Table 4 .
Variation of the levels of the soil reservoir R for the GR4J model.