A global portrait of hydrological changes at the 2050 horizon for the province of Québec

Abstract

This paper presents the methodology and results of a vast study on climate change impacts on hydrology for the province of Québec for the 2050 horizon. A climate ensemble was first built with simulations from the World Climate Research Programme (WCRP)'s Coupled Model Intercomparison Project phase 3 (CMIP3) multi-model dataset, the North American Regional Climate Change Assessment Program (NARCCAP) and the Canadian Regional Climate Model (CRCM) operational runs. Direct outputs and post-processed data from the climate simulations were used as input to the calibrated HSAMI hydrological model in order to produce a large ensemble of hydrological projections for 305 Québec watersheds. Simulations results indicate that increases in mean annual streamflow are projected for the whole province, with greater changes (up to 14%) in the north. The intra-annual distribution of streamflow also changes, with higher winter flows and lower summer flows as well as apparently earlier spring floods. The maximum height of the snow cover and the number of days with snow on the ground are likely to diminish for the 2050 horizon for the southern half of the province, while the northern half will see more snow, but a shorter snow season as well.

Cet article présente la méthodologie et les résultats d'une vaste étude portant sur l'impact des changements climatiques sur l'hydrologie au Québec à l'horizon 2050. Un ensemble climatique a d'abord été construit à partir de simulations provenant de la base de données multi-modèles de la troisième phase du World Climate Research Programme's Coupled Model Intercomparison Project (CMIP3), du North American Regional Climate Change Assessment Program (NARCCAP) et de simulations opérationnelles du Modèle Régional Canadien du Climat (MRCC). Les sorties directes ainsi que post-traitées de ces simulations climatiques ont servi à alimenter le modèle hydrologique HSAMI calibré afin de produire un large ensemble de projections hydrologiques pour 305 bassins versants du Québec. Les résultats des simulations indiquent que des augmentations du débit annuel moyen sont projetées pour l'ensemble du Québec, avec des changements plus importants (jusqu'à 14%) dans le nord. La distribution intra-annuelle du débit change également, avec des débits hivernaux plus soutenus et des débits estivaux diminués, ainsi que des crues printanières devancées. L'accumulation maximale de neige et le nombre de jours avec de la neige au sol diminueront probablement à l'horizon 2050 pour la moitié sud de la province, alors que la moitié nord verra plus de neige, mais une saison plus courte tel qu'au sud.

Introduction

The Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC 2014) confirms that climate change will significantly exacerbate the pressure on water resources. Different uses of the resource such as power generation and drinking water are particularly likely to be affected; trends in hydrometeorological variables have already been observed and further changes are projected.

In Canada, and more specifically in the province of Québec, significant trends in meteorological variables have been observed in the recent past (Yagouti et al. 2008; Desjarlais et al. 2010). In southern Québec, temperatures have increased at a rate of 0.2 to 0.4°C per decade between 1960 and 2005, with greater warming of the daily minimum temperatures than the maximum, resulting in a decrease in diurnal ranges, meaning fewer cold nights, cold days and frost days, and more warm nights and hot days (Vincent and Mekis 2006). The warming has been more important for winter months, impacting snow and ice processes. Yagouti et al. (2008) noted an increase in the number of cycles of freezing and thawing in winter, a decrease in the duration of the frost season and an increased length of the growing season for the period 1960–2005 for the province of Québec. Moreover, Vincent and Mekis (2006) estimated that there has been a decrease of about 10 to 30 frost days (days with minimum temperature below 0°C) in southern Québec over the twentieth century.

In terms of precipitation, Yagouti et al. (2008) observed an increase in annual quantities in Québec between 1960 and 2005, although several stations show a decreasing trend in summer. According to the same study, the number of rainy days and low-intensity events also increased, and the number of days with snowfall as well as snowfall amounts decreased for southern Québec, while they increased in the north.

These changes in meteorological variables directly affected hydrological variables. Whitfield and Cannon (2000) conducted a study of trends in observed streamflows for 642 hydrometric stations in Canada for the decades of 1976–1985 and 1986–1995. The authors noted that several regions, including southern Québec, showed hydrographs with earlier spring floods, more sustained winter flows and lower flows in summer. Burn et al. (2010) noted that there is no significant trend in terms of the mean annual flow in Québec between 1957 and 2006; however, they found that the spring flood is earlier in the North Shore region. Brown (2010) showed that maximum snow water equivalent and duration of snow cover decreased in the south and increased in central and northern Québec, from 1948 to 2005. Fernandes et al. (2007) showed that there was an upward trend in actual evapotranspiration for the period 1960–2000 due to higher temperatures and increased radiation across Canada. In Québec, anomalies or deviations around the mean show an upward trend ranging from 0.162 mm/year and 1.037 mm/year for the Arctic tundra and the Atlantic Canada, respectively. The trend is about 0.7 mm/year for the rest of Québec.

These observed changes in hydrometeorological variables are generally representative of the projections issued by climate models. It is projected that by 2050, the temperature of the entire province will increase, with greater changes in the north and for the cold season (Desjarlais et al. 2010). In summer, the temperature rise will be relatively similar throughout the province according to the climate models. Logan et al. (2011) indicate that the number of freeze-thaw events will increase in the province of Québec at the 2050 and 2090 horizons, with more important changes in the south for the winter season, and in the north for the spring season. In terms of precipitation, projected increases are greater for the north and the cold season (Desjarlais et al. 2010). In summer, an increase in rainfall is projected north but not south. The increase in winter precipitation would result in a greater accumulation of snow in the north, whereas in the south, a decrease of this accumulation is projected because of the rising temperatures and shortening of the cold season.

Brown and Mote (2009) found that decreases in the maximum height of snow cover would be more important in the maritime areas of Québec. They also found that the most significant increases occur in inland areas where the climate is generally dry and cold, as in central Québec. The authors also noted that the percentage of climate model scenarios indicating a decrease in the duration of snow cover increases with the time horizon. In the central and northern parts of the province of Québec, few models anticipated decreases in maximum snow cover for the 2020, 2050 and 2080 horizons.

Anticipated impacts of climate change on streamflows in Québec are periodically evaluated by Hydro-Québec (HQ) throughout the province as well as by the Centre d'expertise hydrique du Québec (CEHQ, Ministry of Sustainable Development, Environment and the Fight against Climate Change) for the southern part of Québec. The CEHQ recently published a hydroclimatic atlas (CEHQ 2013) depicting the projected changes in high and low flows at the 2050 horizon. The data used in the assessment comprise a total of 89 climate projections (direct outputs and post-processed) and the modeling domain consisted of a region called meridional Québec which includes the watersheds of the Ottawa and Saguenay rivers as well as the Saint-Lawrence Valley region. They concluded that the spring flood will occur earlier and produce smaller volumes over the modeling domain. Summer and fall high flows could increase north of the modeling domain, corresponding to the Lac St-Jean area, and low flows for the same seasons could be more severe in terms of amplitude and duration. Overall, a shift in seasonal patterns is projected, with an increase of winter flows and a decrease of summer flows.

In 2008, Desrochers et al. (2008) evaluated the hydrological impacts of climate change for the 2050 horizon on HQ's managed water systems, and for some hydrometric stations in undisturbed rivers of the north of the province. They drew their conclusions from 90 climate simulations based on a combination of global models and greenhouse gas emissions scenarios, post-processed using a perturbation method (Hay et al. 2000). They concluded that, on average, streamflow would increase across Québec, with larger increases in the north compared to the south.

Studies from Desrochers et al. (2008) and CEHQ (2013) constitute the most recent work aiming at establishing the impact of climate change on hydrological regimes for the province of Québec, and few studies have been published with regards to hydroelectricity production and climate change, as underlined in Boucher and Leconte (2013). So far, no study has been devoted to simultaneously evaluating (1) the consequences of climate change for all hydrometric stations in undisturbed rivers and managed water systems in the province, (2) the impact of climate change on key hydrologic variables such as snow cover and actual evapotranspiration over the whole province and (3) not only the median changes but also the interscenario dispersion and the direction of change. The work presented in this paper is therefore intended to depict a portrait of climate change impacts on hydrology for the province of Québec based on these aspects.

In order to achieve this, a large ensemble of climate scenarios was used as input to the HSAMI hydrological model (Fortin 2000) to produce a set of hydrological projections. The model was calibrated specifically for the needs of the study in order to obtain robust parameters in terms of water balance (Minville et al. 2014). Hydrological outputs consisting of simulated streamflow, snow accumulation and actual evapotranspiration were subsequently analyzed to draw the portrait of projected changes in these variables over Québec.

The next part of this paper briefly presents the study area and the hydrometric data used in the study. The third part describes the methodological aspects of the study, covering the building of the climate simulation ensemble as well as the hydrological simulation framework. Finally, the last section presents the resulting hydroclimatic projections at the provincial scale.

Study area

The principal objective of the study being a global assessment of climate change impacts on Québec's hydrology, a database with the greatest available number of watersheds distributed over the province had to be put together. Since meteorological gridded data are readily available at the provincial scale, the main limiting aspect for the construction of the database was the availability of good-quality hydrometric data. Overall, 305 watersheds were chosen on the basis of the length and reliability of their hydrometric time series and their strategic importance (e.g. watersheds managed by hydropower plants). The selected watersheds are depicted in Figure 1. Hydrometric data were collected from HQ, CEHQ and Rio Tinto Alcan (RTA). The database thus includes rivers managed by hydraulic structures, for which streamflow (or net basin supply) is reconstructed using a water budget equation, as well as undisturbed rivers for which streamflow is measured at hydrometric stations. Watershed areas vary between 10 and 69,000 km² and time series range from 5 to more than 100 years. However, only streamflow data concurrent with meteorological data were kept for the calibration of the hydrological model, reducing the time series to a maximum of 45 years. A more extensive description of the database can be found in Guay and Minville (2012).

Figure 1. Study area.

Methodology

The assessment of the impact of climate change on hydrology requires the use of a modeling cascade, as shown in Figure 2. Climate variables that serve as inputs to hydrological models come from different simulations conducted by global or regional climate models. The climate models themselves are forced with CO₂ concentrations from greenhouse gas emission scenarios (SRES) (Nakicenovic et al. 2000). Global models simulate the physical processes taking place in the atmosphere, oceans and land surface across the planet. Regional models, for their part, simulate the same processes, and sometimes additional phenomena, at finer resolutions and on restricted domains. They are driven by global models at their boundaries. However, temperatures and precipitations from climate models often show a bias when compared to observations from the control period. Using the climate variables without bias correction to drive the hydrological model (direct output method) has two major limitations: the potential instability of the hydrological model when using inputs different from those used for the calibration, and the interpretation of the results on a relative (future versus control period) basis only. This means that simulated future hydrological regimes possibly are biased, and cannot be used directly for most impact studies. The first limitation depends on the hydrological model structure, the robustness of its calibration parameters and the nature and amplitude of the bias between observed and simulated climate variables. In order to address both limitations, variables can be corrected for the bias, or future series must be constructed from observations and simulated changes. These two post-processing techniques are called bias correction and perturbations, respectively. Climate scenarios obtained by post-processing and the direct method are then used as inputs to hydrological models to simulate future hydrological regimes.

Figure 2. Schematics of the modeling chain.

Every step of the modeling cascade brings uncertainty to the final projected hydrological regimes, which is often described as profitable since a large set of possible outcomes allows for a comprehensive assessment of probabilities and confidence in the results (Wilby and Harris 2006; Wilby and Dessai 2010). In this study, a statistically indistinguishable interpretation of the climate ensemble uncertainty is adopted (Annan and Hargreaves 2010; Sanderson and Knutti 2012). This implies that the true outcome is considered to be drawn from a distribution of models, rather than considering each model as an approximation of the true outcome (climate).

Climatic scenarios

Simulations

A climate ensemble was produced at the Ouranos Consortium for this study and consists of simulations from the World Climate Research Programme (WCRP)'s Coupled Model Intercomparison Project phase 3 (CMIP3) (Meehl et al. 2007) multi-model dataset, the North American Regional Climate Change Assessment Program (NARCCAP) (Mearns et al. 2012) and the Canadian Regional Climate Model (CRCM) (Music and Caya 2007; deElía and Côté 2010; Paquin 2010) operational runs from Ouranos. The 87 simulations of the set are based on three greenhouse gas emission scenarios (B1, A1B, A2), 15 global models (one to five members) and four different regional models. Spatial resolutions vary from 0.4 to 5.0° (approximately 45 to 550 km). The control period is 1961–2000 (or 1970–2000 for a few NARCCAP simulations) and the future period is the 2050 horizon (2041–2070 or 2046–2065, depending on availability). The use of two different control periods is not ideal, but the two periods are considered to be representative of the end of the twentieth century. Table 1 lists the characteristics of the 87 climate simulations. For regional climate model (RCM) simulations, the global climate model (GCM) is the driver, and the countries and agencies are related to the RCM.

Table 1. Climate simulations characteristics.

No.	Ensemble	SRES	GCM	Members	RCM	Country (agency)
1–2	CMIP3	A1B, B1	GISS-AOM	1		USA (NASA)
3–4	CMIP3	A2, B1	BCM2	1		NOR (BCCR)
5–7	CMIP3	A2, A1B, B1	CSIRO-MK3	1		AUS (CSIRO)
8–10	CMIP3	A2, A1B, B1	CSIRO-MK3.5	1		AUS (CSIRO)
11–12	CMIP3	A1B, B1	CGCM3.1 (T63)	1		CAN (CCCma)
13–27	CMIP3	A2, A1B, B1	CGCM3.1 (T47)	5		CAN (CCCma)
28–30	CMIP3	A2, A1B, B1	CM2.0	1		USA (GFDL-NOAA)
31–33	CMIP3	A2, A1B, B1	CNRM-CM3	1		FRA (CNRM)
34–36	CMIP3	A2, A1B, B1	IPSL-CM4	1		FRA (IPSL)
37–38	CMIP3	A2, A1B	ECHAM4	1		DEU (MPIM)
39–41	CMIP3	A2, A1B, B1	ECHAM5	1		DEU (MPIM)
42–50	CMIP3	A2, A1B, B1	ECHO-G	3		DEU (MIUB), KOR (KMA)
51–52	CMIP3	A1B, B1	MIROC3.2 (high)	1		JPN (CCSR/NIES/FRCGC)
53–58	CMIP3	A2, A1B, B1	MIROC3.2 (medium)	2		JPN (CCSR/NIES/FRCGC)
59–73	CMIP3	A2, A1B, B1	MRI-CGCM2.3.2	5		JPN (MRI)
74	NARCCAP	A2	CGCM3.1 (CAN)	1	CRCM	CAN (Ouranos, UQAM)
75	NARCCAP	A2	CGCM3.1 (CAN)	1	RCM3	USA (UC Santa Cruz)
76	NARCCAP	A2	CM2.0 (USA)	1	RCM3	USA (UC Santa Cruz)
77	NARCCAP	A2	CCSM3 (USA)	1	CRCM	CAN (Ouranos, UQAM)
78	NARCCAP	A2	HadCM3 (GBR)	1	HadRM3	GBR (Hadley Center, MO)
79–83	Ouranos	A2	CGCM3.1 (CAN)	5	CRCM	CAN (Ouranos, UQAM)
84	Ouranos	A1B	CNRM-CM3 (FRA)	1	CRCM	CAN (Ouranos, UQAM)
85–87	Ouranos	A2	ECHAM5 (DEU)	3	CRCM	CAN (Ouranos, UQAM)

CMIP3, Coupled Model Intercomparison Project, Phase 3; NARCCAP, North American Regional Climate Change Assessment Program; USA, United States of America; NASA, National Aeronautics and Space Administration; NOR, Norway; BCCR, Bjerknes Centre for Climate Research; AUS, Australia; CSIRO, Commonwealth Scientific and Industrial Research Organisation; CAN, Canada; CCCMa, Canadian Centre for Climate Modelling and Analysis; GFDL-NOAA, Geophysical Fluid Dynamics Laboratory - National Oceanic and Atmospheric Administration; FRA, France; CNRM, Centre National de Recherches Météorologiques; ISPL, Institut Pierre Simon Laplace; DEU, Germany; MPIM, Max Planck Institute for Mathematics; KOR, Korea; KMA, Korea Meteorological Administration; JPN, Japan; CCSR, Division of Climate System Research; NIES, National Institute for Environmental Studies; FRCGC, Frontier Research Center for Global Change; MRI, Meteorological Research Institute; UQAM, Université du Québec à Montréal; UC Santa Cruz, University of California, Santa Cruz; GBR, Great Britain; MO, Met Office.

Post-processing

Two methods were chosen for the post-processing of climate variables on the basis of recommendations from Mpelasoka and Chiew (2009) and Willems and Vrac (2011) for applications in hydrology. The two methods are called daily scaling (perturbation of observations or delta method) and daily translation (bias correction of climate model output). For the purpose of post-processing, observations consist of basin averages, while simulated climate variables are computed as the average of the climate model grid points comprised in the basin, or the nearest land grid point in the case where no grid point is found inside the basin.

In the daily scaling method, the change (Δ) of the simulated climate variable (temperature, T_sim, and precipitation, P_sim) between the future horizon (h50) and the control period (ctl) is first calculated with Equation (1). A change is calculated for each month “m” of the year and each quantile “q” for the month. Quantiles are drawn from empirical distribution functions of the daily values of each month over a period. Changes in temperature for a particular month and a particular quantile are calculated as the difference between the average temperature for this month and this quantile over the future period, and the average temperature for this month and this quantile over the control period. Changes in precipitation are calculated similarly using a ratio rather than a difference. Daily temperatures for the future period (T^h50) are obtained by adding the calculated deltas to the observed temperatures (T^obs), whereas future precipitations (P^h50) are computed as observed precipitations (P^obs) multiplied by the calculated deltas, as shown in Equation (2).(1) (2)

This perturbation of observed data has the advantage of preserving the usually higher variability found in observed data compared to simulated climate data. However, the method also reproduces the temporal sequence of observed meteorological events and does not uphold the dynamics of the climate simulation.

Daily translation consists in first determining the bias (ε) of the climate variables on the control period by comparing the simulations and observations, as described in Equation (3). Again, the daily data are divided into quantiles for each month. Assuming that the bias is conserved from the control period to the future period, the climate scenario is then created by adjusting the bias of the climate simulation over the two periods, as described by Equations (4) and (5). Similarly to the daily scaling method, the correction of the temperature is additive, whereas the precipitation is multiplicative.(3) (4) (5)

When using the bias between a reference period simulation and observations to correct the future simulation bias, the different characteristics and dynamics of each climate model are conserved. Each resulting climate scenario will have its own temporal sequence of meteorological events which are different between member simulations. Both post-processing methods assume the biases to be of equal magnitude in the future and reference periods.

Hydrological simulations

HSAMI model

Hydrological modeling was performed with the HSAMI model (Bisson and Roberge 1983; Fortin 2000), which was developed by HQ and has been routinely used in forecasting natural inflows for over 30 years. HSAMI is actually used in Québec for daily forecasting of natural inflows on 87 watersheds with surface areas ranging from 160 to 69,000 km².

HSAMI is a 23-parameter, lumped and conceptual model. Two parameters account for evapotranspiration, six for snowmelt, 10 for vertical flow and five for horizontal flow. The required daily input data for the model are basin-averaged minimum and maximum temperatures, as well as liquid and solid precipitation. Sunshine hours and measured snow water equivalent are optional inputs. Vertical flows are simulated with four interconnected linear reservoirs (snow on the ground; surface water; unsaturated and saturated zones). The model takes into account snow accumulation, snowmelt, soil freezing/thawing, infiltration and runoff, and evapotranspiration. Water is transferred at the basin outlet through surface runoff, interflow and baseflow. For each time step, the model goes through the following steps: (1) estimation of potential evapotranspiration based on temperature; (2) computation of net precipitation on reservoirs (optional); (3) simulation of interception and accumulation of rainfall and snow, snowpack freezing/thawing and aging; (4) separation of available surface water between infiltration and runoff; (5) simulation of vertical flows (infiltration, interflow, evapotranspiration and interaction between the saturated and unsaturated zones); (6) simulation of horizontal flows toward the outlet; and (7) computation of natural inflows or water discharge at the outlet.

HSAMI has been widely used to assess climate change impacts on local hydrology in recent years (Minville et al. 2008; Minville et al. 2009; Boyer et al. 2010; Chen et al. 2011b; Laforce et al. 2011; Beauchamp et al. 2013).

Model calibration and validation

As underlined by Brigode et al. (2013), hydrological modeling in non-stationary conditions involves additional sources of uncertainty that must be addressed. In addition to the common errors in model structure, parameter instability due to changes in catchment characteristics, weather patterns or dominant processes may occur. The calibration procedure leading to the identification of model parameters should account for such potential instabilities.

The formulation of potential evapotranspiration (PET) has been identified as a process fairly sensitive to this issue. As part of a study on the Great Lakes, Lofgren et al. (2011) showed that the use of temperature as a proxy to assess the PET in a conceptual hydrological model could lead to differences in the range of 200 mm/yr between the change in latent heat flux calculated by the hydrological model and the one simulated by the climate model. This was also emphasized by Milly and Dunne (2011) in a study of 10 US watersheds using three different climate models, where changes in PET were typically 3 times greater when simulated by the temperature-based PET included in the hydrological model compared with the change in PET issued by climate models.

In the HSAMI model, the formulation of PET is based solely on temperature. The actual evapotranspiration (AET) is calculated as a proportion of PET, according to a calibration parameter, and is constrained by the water available at the surface and in the soil reservoir. It was therefore considered necessary to develop a calibration strategy that could minimize the behavior observed in the aforementioned studies. The strategy, described in detail in Minville et al. (2014), achieves this goal by using the AET from a climate model in parallel with the traditional observed streamflow during the calibration. This is made possible by performing a constrained optimization in which the automatic optimizer works at minimizing an objective function while satisfying one or more constraints. An important assumption behind this approach is that the AET from the climate model is more reliable than the one produced by the hydrological model, since it is based on a radiative balance rather than an empirical equation depending on temperature. This assumption is consistent with the conclusions of Lofgren et al. (2011), who found that an energy budget-based approach to the calculation of PET better satisfies the conservation of energy than a formulation based on air temperature. In the present study, outputs from a regional climate model (CRCM4) driven by a reanalysis (ERA-40; Uppala et al. 2005) are used in order to reduce bias.

In this study, the objective function is the Nash–Sutcliffe Efficiency (NSE) (Nash and Sutcliffe 1970) applied to the simulated streamflow, used in conjunction with two constraints. The constraints are defined as the minimum NSE values between AET from a climate model and AET simulated by HSAMI, for two periods of the year. The daily AET being very noisy, the NSE values are calculated on the basis of annual cycles. The average annual cycle is divided into two periods (low and high AET) to compensate for the fact that the NSE assigns more weight to higher values, and therefore does not target physical processes associated with low AET such as sublimation or evaporation of snow water during melting. The AET cycle was separated by an approach adapted from Madsen (2003). According to this method, the periods of low evapotranspiration consist of days when the AET is less than half the annual average.

The splitting of the observed data into a calibration and a validation period was performed objectively using even years for the period of calibration and odd years for the period of validation. This method avoids the calibration of the model over wet or dry spells which could influence the choice of parameters, as highlighted in Minville (2010). This method is best suited for long time series, but was also applied to watersheds with shorter time series for the purpose of uniformity.

The optimization of HSAMI parameters was performed using the EAUptim tool, developed at the Institut de recherche d’Hydro-Québec (IREQ) (Leclaire 2008). EAUptim is a flexible tool that allows the connection between a hydrological model and an optimization algorithm. The Nonlinear Optimization by Mesh Adaptive Direct Search (NOMAD) algorithm (Le Digabel 2011) was used since it provides many options such as the imposition of constraints and the Variable Neighborhood Search (VNS) (Audet et al. 2008) which allows to escape from local minima. NOMAD is a C++ implementation of the Mesh Adaptive Direct Search (MADS) of Audet and Dennis (2006). A more detailed description of the use of NOMAD with the calibration strategy of the study is available in Minville et al. (2014).

Temperature, precipitation and AET from a CRCM operational run driven by the ERA-40 reanalysis (Uppala et al. 2005) were used for the computation of constraints on AET. The simulation is available at a resolution of 45 km × 45 km and was archived at daily temporal intervals.

For the computation of NSE values on streamflow, observed temperature and precipitation were extracted from interpolated gridded data produced at IREQ (Tapsoba et al. 2005), which are available at a resolution of 10 km × 10 km.

Figure 3 shows model performance with regards to streamflow and AET. Panels (a) and (b) show obtained NSE values for streamflow for the calibration and validation periods respectively, for the 305 watersheds of the study area. Values smaller than 0.5 are presented in different shades of yellow and red and are generally located in the southwest portion of the province. The subdued calibration performance for these watersheds can be due to various factors related to hydrometric data (short time series, missing data, uncertainty of calculated net basin supplies, lack of coherence between hydrometric and meteorological data, etc.) or climatic data. In the latter case, the parametric space restrained by the constraints could be inappropriate or insufficient to produce good results for streamflow. This could happen, for instance, if there were a conflict between the dynamics of evapotranspiration in the climate model and reality, caused by erroneous representation of watershed characteristics related to evapotranspiration (vegetation, land use, soil type, etc.) in the climate model.

Figure 3. Model performance. (a) Nash–Sutcliffe Efficiency (NSE) values for streamflow per basin for the calibration period. (b) NSE values for streamflow per basin for the validation period. (c) Boxplot of NSE values for streamflow (Q) over the calibration and validation periods (all basins). (d) Boxplot of NSE values for actual evapotranspiration (AET, low and high) for the calibration and validation periods (all basins).

Panel (c) presents the distribution of the NSE values for streamflow (Q) for all basins using a boxplot for each period. The median performances are 0.75 and 0.72 for the calibration and validation periods, respectively. The dispersion around the median for both periods is 0.15.

Panel (d) also shows distributions of NSE values using boxplots, but for AET. NSE values are calculated with the AET from the CRCM simulation as observed values, and the AET simulated by HSAMI forced with CRCM meteorological data as simulated values. As described earlier, the constraints were applied separately on the low and high portions of the annual AET cycle. NSE values are thus presented for both portions. Performance for the low-AET portion of the year is generally less than for the high-AET portion, and there isn't much difference between calibration and validation periods for both portions of the cycle. Median NSE values for the low AET are 0.58 and 0.61, whereas the high AET shows values of 0.73 and 0.74 for the calibration and validation periods, respectively.

Hydrological simulations

Hydrological simulations for each of the 305 watersheds were conducted with the calibrated HSAMI model and climatic inputs from the previously described 87 climate simulations (direct outputs), 87 climate scenarios produced with the daily scaling method and 87 scenarios produced with the daily translation method, leading to a total of 261 scenarios for each watershed.

Results and discussion

The results presented below are summarized using three indicators allowing a global portrait of the 261 scenarios for each watershed. The first indicator, shown in Figure 4a, is the median of the empirical distribution of the 261 change values simulated for a watershed. It provides a central estimation for change and should not be used as an accurate value without consideration of the uncertainty it carries. The second indicator is called the interscenario agreement on the direction of change and is shown in Figure 4b. It is the percentage of scenarios projecting an increase. Based on recommendations from the IPCC'S fourth assessment report (IPCC 2007), the following likelihood ranges, or probabilities of occurrence, are used: very probable increase, > 90%; probable increase, > 66%; no signal, between 33% and 66%; probable decrease, < 33%; very probable decrease, < 10%. Finally, the third indicator is the dispersion of the empirical distribution (Figure 4c) calculated as the difference between the 75th percentile and the 25th percentile. The color scales of the figures indicate more water or cold in blue shades, while the red and yellow shades are associated with increases in heat or aridity.

Figure 4. Illustration of the three indicators for change. (a) Median change; (b) interscenario agreement on the direction of change; (c) interscenario dispersion.

Climate change

Figure 5 illustrates changes in minimum and maximum temperatures and precipitation over the study area. In every case, the direction of change indicates a very probable increase over the entire territory.

Figure 5. Median change (left), interscenario agreement on the direction of change (center) and interscenario dispersion (right) for (a) minimum daily temperature, (b) maximum daily temperature, (c) annual precipitation and (d) annual streamflow.

The minimum temperature increases more than the maximum temperature with values around 4.5°C, in comparison to 2.5°C. The increase in maximum temperature follows a north–south gradient whereas the minimum temperature does not show any particular spatial pattern. The interscenario dispersion for the maximum temperature is greater in the north, where dispersion values even exceed values of median change. The projected increase in precipitation in northern Québec reaches 15%, and the southern portion of the province shows an increase slightly below 10%. Dispersion values for the change in precipitation are slightly smaller than the median values of change.

Hydrologic changes

Streamflow

Change in mean annual streamflow is also presented in Figure 5. Changes in streamflow follow closely changes in precipitation in terms of spatial distribution, and are spreading from a 2% increase in the south to a 14% increase in the north. In this region, changes in precipitation are thus directly reflected in changes in streamflow, whereas in the south, the increase in streamflow appears attenuated by a probable increase in evapotranspiration. The dispersion does not follow the same spatial pattern, though, being probably affected by multiple factors other than the increase in precipitation like the seasonality and timing in the changes in temperature and precipitation. Finally, the direction of change indicates an overall very likely increase in the north, and a likely increase for the south and the North Shore region. For a few watersheds near the US border, results indicate no signal, with, however, a relatively high percentage of scenarios projecting an increase (60 to 64%). Ratios of annual streamflow to precipitation (not shown in the figure) indicate an overall decrease throughout the province ranging from –5% (south) to –1% (north), meaning that the proportion of annual precipitation actually ending in the rivers will be slightly less at the 2050 horizon.

Figure 6 shows the median change in monthly contribution to annual streamflow. For each basin “b”, each climate scenario “s” and each month “m”, the ratio of monthly flow Q_m to annual flow Q_y is computed over the control and future periods, and the change (Δ) is calculated with Equation (6) in which n is the number of years in the period (typically 20 to 40): (6)

Figure 6. Change in monthly contribution to annual streamflow.

Figure 6 thus illustrates the projected shift in the intra-annual distribution of streamflow, independently of the change in total annual flow. As a complement, Figure 7 presents the direction of change for each month. The months are grouped in threes, starting with December, so that each line represents a season (winter, spring, summer and fall).

Figure 7. Interscenario agreement on the direction of change for monthly contribution to annual streamflow.

For the entire province, winter months are characterized by increases in their contribution to annual flow ranging from slightly more than 0% to 5% with the greatest changes in the south of the province. Although the median increases are fairly small, the direction of change globally indicates probable and very probable increases for these months. This result is consequent with projected milder temperatures and more frequent freeze/thaw events for these months as well as with the projected increase in precipitation.

From March to June, results indicate that the spring flood might happen earlier at the 2050 horizon. For the southernmost basins, for instance, the volume usually occurring in April is rather seen in March, resulting in a decrease in the monthly contribution of April, and a concurrent increase in the contribution of March. The same phenomenon gradually applies northward throughout the season, with the northernmost basins showing a shift of volume from May to June. The direction of change shows a progressive shift from very probable increases to very probable decreases during the season. Further analysis of the simulated hydrographs indicates that the mean spring flood volume (not shown in the figure) decreases by about 5% at the 2050 horizon for the southernmost basins, North Shore and Gaspésie, whereas it increases by about 5% for the rest of the province, yet with a weak consensus for the direction of change for the whole province.

Summer months show an overall very probable decrease in their contribution to the annual volume, ranging from –0.5 to –3% (less probable for the southernmost basins). The usual summer contribution to annual flow already being fairly small, this indicates an increased aridity for these months. This phenomenon extends to September, and even October for the south of the province, but with a decrease in the signal for many basins. The contribution of late fall months to the annual flow for the basins of the north increases, which can be explained by more frequent temperatures above the freezing points and thus more liquid precipitation.

Overall, although median values of changes in monthly contributions are rather small, the consensus with respect to the direction of change is strong. Annual flow distribution shows a subtle shift from summer months to winter months that could be explained by more aridity and dry spells projected for the summer, and more precipitation and temperate days during winter. The spring flood also appears to happen earlier at the 2050 horizon for the whole province.

Snow cover

Changes in snow cover are presented from two different perspectives: the number of days with snow on the ground, and the maximum accumulation of snow (Figure 8). The upper part of the figure provides information on the length of the cold season (although days are not necessarily successive), and the second part integrates the amount of precipitation concurrent with extended spells of cold temperatures, which fosters snow accumulation.

Figure 8. Median change (left), interscenario agreement on the direction of change (center) and interscenario dispersion (right) for (a) number of days with snow on the ground and (b) maximum snow accumulation.

Figure 8 shows decreases in the number of days with snow on the ground varying between 18 and 35 days, the most important changes occurring in the south of the province. The decrease is very probable for the entire study area, as shown in Figure 8, and interscenario dispersion is on the same order of magnitude as the median change itself. The results thus indicate that at the 2050 horizon, the season with snow on the ground will very likely be 2 weeks to a month shorter.

The change in maximum snow accumulation, though, is strongly polarized from north to south, with a projected increase of roughly 5% for the north, evolving southward to a decrease of as much as 30%. The increase in the north could be explained by an increase in precipitation combined with temperatures that remain below the freezing point. The direction of change for this area indicates a probable increase, whereas there is no signal in central Québec, and a strong probability of decrease in the south. The rather high dispersion values and the absence of signal for central Québec could be induced by the high, and presumably increasing, interannual variability of snow processes, and the variability between projections on the rate of migration of the snow-line from south to north.

Evapotranspiration

The change in evapotranspiration is presented along with the hydrological variables since it is simulated by HSAMI as a part of the water budget. Figure 9 depicts the change in the annual ratio of AET to precipitation, which is a convenient way of expressing the proportion of incoming water withdrawn from the surface, and therefore not available for infiltration, runoff and/or streamflow.

Figure 9. Change in annual actual evapotranspiration (AET) to precipitation ratio. (a) Median change; (b) interscenario agreement on the direction of change; (c) interscenario dispersion.

The changes in this ratio are rather small for the whole province, ranging from 0.5 to 2.5%. The increase is likely for most of the basins, apart for some located in the south, in Gaspésie and North Shore where it is highly probable. This is in accordance with the decrease in summer contribution to annual flow presented in Figure 6. Dispersion values are the same order of magnitude as the median change values, but generally higher than the change itself. The small changes in the AET/P ratios are consequent with the climate of Québec where AET is limited by the availability of surface water and energy. If more water is available from precipitation, more will be evaporated, given the right amount of energy. Absolute changes in AET are in the range of 10 to 22% for the province, with greater changes in the north, which is the same order of magnitude as the changes in precipitation and streamflow.

Conclusion

A global assessment of climate change impacts on Québec's hydrology was carried out over 305 watersheds. Eighty-seven climate simulations from a set of global and regional models served as input to the HSAMI hydrological model. Direct and post-processed outputs from the climate simulations were used, providing a large number of scenarios. The hydrological simulations were carried out with the HSAMI model, which was calibrated with special attention to the reproduction of the AET annual cycle.

Overall, the projected increases in temperatures and precipitation result in a rather proportional increase in annual flow volumes and annual evapotranspiration at the 2050 horizon, with differences in the way water is distributed throughout the year. The results underline the typical north–south gradient of the province, and illustrate the migration of the processes generally experienced in the south towards the northern part of the province. The increase in winter flows and the probable decrease in summer flows will constitute a challenge for water management, especially in the south of the province. Substantial decreases in snow amounts and the number of days with snow on the ground are also projected for the south. The northern part of Québec will also see a decrease in the number of days with snow on the ground, but total snow accumulation will be greater. The changes in the hydrograph pattern (earlier spring flood, increased annual streamflow and higher flow in winter) are consistent with the findings of other studies on snowmelt-dominated Québec watersheds (Dibike and Coulibaly 2004; Chen et al. 2011b). These conclusions are also overall consistent with Ouranos's most recently published synthesis report (Ouranos 2014) and CEHQ's hydroclimatic atlas (CEHQ 2013) in terms of changes in streamflow, from south to north. However, the projected changes in maximum amounts of snow for the north obtained in the present study are not in accordance with Ouranos's report, in which a small decrease in maximum snow water equivalent is projected for the 2050 horizon. The most probable reason for this difference is that the results in the Ouranos report are based on 19 climate scenarios using solely the most severe representative concentration pathway (RCP8.5) (Vuuren et al. 2011) issued by the IPCC. In the present study, the most severe SRES scenario used is A2, and it is presented along with the other, and less severe, scenarios. This highlights the potential influence of the methodological choices made during the simulation process shown in Figure 2 on the results obtained.

One of the methodological choices made in this study was to attribute equal weight to hydrological variables issued from regional and global climate models during the result analysis. Spatial resolution is likely to have an impact on the simulated variables, especially when using a distributed hydrological model (which was not the case here) or when modeling small watersheds. Further analysis of the hydrological change signal with regards to the spatial resolution of the climate model used and the size of the watershed would be of interest to validate the methodology used in this study, and to assess whether this parameter has an impact at the temporal and spatial scales involved. The over-representation of some global models with multiple members in the ensemble should also be studied.

Results from direct and post-processed outputs from climate models were also shown without differentiation in this study to simplify the presentation and to benefit from a large ensemble of scenarios. However, bias correction (or not) of the climate variables can also influence the simulated hydrological variables. Figure 10 shows side by side the median change in mean annual streamflow projected by the three methods. Although small differences appear, especially with the daily scaling method in the south, overall amplitude and pattern of change are fairly similar for the three methods. Greater differences between the post-processing methods are observed for the change in maximum snow accumulation, as shown in Figure 10. For this variable, the direct output method produces rather subdued changes compared with the other methods, and even projects probable decreases where the others project probable increases. The sensitivity of snow accumulation to bias correction can be explained by the strong link between this process and temperature. Nonetheless, these figures illustrate the need to further assess the impact of the post-processing method on the resulting change in various hydrological indicators.

Figure 10. Median change by post-processing method (direct outputs on the left, daily scaling in the center and daily translation on the right) for (a) annual streamflow and (b) maximum snow accumulation.

Another factor impacting the results is the choice of the hydrological model. In this study, only one model with one set of parameters for each watershed was used, making uncertainty analysis related to the model choice or its parameterization impossible. The structure of the model and the state variables it simulates are obviously related to its ability to project realistic hydrological scenarios for the future. For instance, although the HSAMI model simulates soil thawing/freezing, it does not account for complex permafrost dynamics, which impact the hydrological cycle in the northern part of the province. The choice of a variety of models, accounting for a variety of processes, could be of interest to assess the sensitivity of the change signal to the hydrological model and its structure over the study area.

Finally, the projections presented in this paper were for the 2050 horizon only, based on the availability of the climate simulations at the time the study was conducted. Recent availability of the Coupled Model Intercomparison Project phase 5 (CMIP5) (Taylor et al. 2011) simulations will trigger the production of updated hydrological projections for the 2050 horizon, but also for 2080. Furthermore, the CMIP5 ensemble provides continuous simulations up to 2100 that will allow the follow-up of hydrological variables through time.

Acknowledgements

The authors would like to acknowledge the Ouranos Consortium for providing climate scenarios from the CRCM and building the climate ensemble for this project. CEHQ and Rio Tinto Alcan are also acknowledged for providing streamflow data. Special thanks are addressed to Frédéric Guay (IREQ) for the support provided during the building of the watersheds database and to James Merleau (IREQ) for his expertise in statistical analysis of the results.

We also acknowledge the modeling groups, the Program for Climate Model Diagnosis and Intercomparison (PCMDI) and the WCRP's Working Group on Coupled Modelling (WGCM) for their roles in making the WCRP CMIP3 multi-model dataset available. The authors would also like to thank the North American Regional Climate Change Assessment Program (NARCCAP) for providing the data used in this paper.