Modelling climate change impacts on the hydrology of an agricultural watershed in southern Québec

Abstract

Spring flooding and water availability during the summer are two major concerns in southern Québec. Climate change is expected to further exacerbate shortages of fresh water resources in Québec. Therefore, the impacts of climate change on the hydrology of the Pike River watershed were projected for the 2041–2070 time horizon using a version of the Soil and Water Assessment Tool (SWATqc) adapted to Québec’s agroclimatic conditions. SWATqc was driven by inputs from four climate projections created with the fourth-generation Canadian Regional Climate Model (CRCM4) and the Arpège regional climate model. CRCM4 was driven at the border of its domain by three members of the Canadian Coupled Global Climate Model (CGCM3). The Arpège model, which has a variable grid, is not driven by another model. The gas emissions scenario used to drive climate models was scenario A2. Regional climate data sets were bias-corrected with the delta method prior to being used to drive SWATqc. Projected changes in precipitation and temperatures were analysed for their effects on evapotranspiration, surface and subsurface runoff, and streamflow. Mean annual streamflow was projected to increase between 9 and 19%. All future simulations projected significant increases or decreases in evapotranspiration, surface and subsurface runoff, and streamflow at the monthly time scale. Winter streamflow, for example, increased by 2- to 3-fold compared with the historical values. In addition, future spring floods started earlier and peaked in March instead of April. The four simulations were named ADC, ACU, AFA/AFD, and ARP. All but the ACU simulation projected a decrease in peak flow, while the ADC, AFD and ARP simulations suggested an increase. Although the majority of the simulations pointed to decreases in summer streamflows, the changes were not all significantly different from historical values. Results from this study not only show that climate change has a strong potential to exacerbate seasonal water availability and flooding issues in the watershed, but also provide conservative numbers that could be used in developing improved watershed management plans adapted to future climate change impacts.

Les inondations et la disponibilité d’eau douce en été sont deux préoccupations majeures dans le sud du Québec. On craint que les changements climatiques aggravent les déficits d’eau douce au Québec. Les impacts des changements climatiques sur l’hydrologie du bassin de la rivière Aux Brochets ont donc été simulés pour l’horizon 2041–2070 avec le modèle hydrologique SWATqc (Soil and Water Assessment Tool pour le Québec) adapté pour les conditions agroclimatiques du Québec. Quatre projections climatiques produits avec le Modèle Régional de Climat Canadian (MRCC4) et le modèle de climat français Arpège ont été insérés dans SWATqc. Le MRCC4 était piloté par trois membres du Model de Climat Global Couplé (MCCG3). Arpège, dont la grille de résolution spatiale varie, n’est pas piloté par un autre modèle. Le scenario d’émission de gaz à effet de serre utilisé pour piloter les modèles climatiques est le scénario A2. Le biais des sorties des modèles climatiques a été corrigé avec la méthode du delta. Les changements en température et précipitation projetés ont été analysés pour leur effet sur l’évapotranspiration, le ruissellement de surface, l’écoulement souterrain et le débit. Les simulations obtenues avec SWATqc et les quatre projections climatiques projetèrent une augmentation significative du débit annuel moyen de 9 à 19%. Sur une base mensuelle, les débits hivernaux seront de 2 à 3 fois plus élevés que les débits actuels selon les quatre simulations. La crue printanière se fera plus hâtive avec un pic en mars plutôt qu’en avril. Les quatre simulations ont été nommées ADC, ACU, AFA/AFA et ARP. La simulation ACU suggère une augmentation du pic printanier alors que les simulations ADC, AFD et ARP suggèrent une diminution. Bien que la majorité des simulations montrent une diminution des débits en été, les changements ne sont pas significatifs pour chaque mois et chaque scénario. Les résultats de cette étude non seulement montrent que les changements climatiques peuvent fortement accentuer la saisonnalité de la disponibilité de l’eau fraîche et des inondations dans le bassin mais aussi fournit des chiffres conservateurs sur lesquels se baser pour développer un plan amélioré de la gestion de l’eau tenant compte des changements climatiques.

Introduction

Climate change is expected to negatively impact the quantity and quality of the fresh water resources in Québec (Bates et al. 2008; Ouranos 2010). Two major concerns related to water quantity in southern Québec are spring flooding and water shortages during summer. Spring rainfall and snowmelt in this region often exceeds evapotranspiration and can lead to major flood events such as the flooding of the Richelieu River valley in 2011 (Madramootoo 2011). In Québec, the flood damaged 6000 houses and 6000 ha of corn fields, and resulted in the evacuation of 1000 persons and about CAD$78 million in financial losses (Environment Canada 2013).

The potential for spring flooding notwithstanding, low summer streamflows account for the drawdowns observed in some rivers in this region (D’Auteuil 2011). Southern regions of Québec, on average, face a water deficit (evapotranspiration exceeds precipitation) ranging from 33 to 92 mm between May and October (Agroclimatic Atlas of Québec 2012). Projections for the 2050 horizon from the Canadian Regional Climate Model (CRCM) show that climate-induced changes could exacerbate this deficit by 92 to 113 mm (the lower scenario from a series of 15 simulations) or attenuate this deficit by 17 to 37 mm (higher scenario) (Agroclimatic Atlas of Québec 2012). The Atlas also projects a possible 7% increase in the intensity of annual maximum precipitation (of 12 h and 24 h duration) which could lead to increased runoff rather than soil and aquifer recharge. This might limit water availability for irrigated agriculture, especially in high-elevation areas with low water tables, as well as other parts of the Montérégie region where vegetables are grown (D’Auteuil 2011).

Mindful of the importance of water resources for the region’s human, ecosystem and economic health, the Québec government has enacted legislation aimed at the development of watershed master plans (WMPs). This initiative encourages basin stakeholders to define actions aimed at the protection and management of water resources (Auger and Baudrand 2004). Most of these WMPs, however, do not take into account the potential impacts of climate change on water availability. An understanding of the potential impact of climate change on hydrological processes is an essential first step towards the development of management practices for the protection of watershed water quantity in a future changing climate (Stager and Thill 2010). The Organisme de Bassin Versant de la Baie Missisquoi (OBVBM), a local watershed management body within our study area, seeks to initiate actions to prepare and support its population in incorporating climate change adaptation into the watershed’s WMPs (OBVBM 2011).

Assessments of impacts of climate change on water quantity are typically performed by forcing hydrological models with time series of climate data produced by one or several regional climate models (RCMs) for a reference and future period. RCMs improve global climate model (GCM) simulations at the regional and local scales by dynamically downscaling GCM climate variables (Fowler et al. 2007). RCMs better represent the local topography, land use and interactions between the land surface and the atmosphere. The strong influence of these components on hydrologic variables such as rainfall suggests that RCM-based estimation of local precipitation should generate more accurate data. RCMs outputs, however, are known to encompass biases (Boé et al. 2009; Gagnon et al. 2009) and should preferably be corrected prior to introduction into hydrological models (Bordoy and Burlando 2012; Schoetter et al. 2012; Teutschbein and Seibert 2012). Although both simple and complex approaches to bias correction have been explored (Diaz-Nieto and Wilby 2005; Graham et al. 2007; Choi et al. 2009), most of the existing methods not only add some degree of uncertainty to the results but also have inherent limitations (Boé et al. 2009). The delta correction methods, for example, do not account for changes in precipitation frequency (Boé et al. 2009; Sulis et al. 2011). Precipitation frequency and intensity are important factors in flood and water shortage forecasting. Complex downscaling methods, on the contrary, account for changes in precipitation intensity and frequency but require high statistical expertise and computing resources that are not often available to watershed management bodies. The choice of bias correction approach in climate impact studies therefore depends on the availability of resources and the purpose of the project outcomes.

Simple bias correction methods that generate reproducible data, using limited resources, were therefore chosen for this study in order to provide small watershed management bodies with reliable and sustainable tools that could be used to incorporate climate change impacts into the development of WMPs. Robust climate change signals that stand out from modeling uncertainties are believed to be generated by this approach, which still allows the results to be used in decision making provided these uncertainties are acknowledged and assessed (Boé et al. 2009).

The aims of this study were twofold: (1) to assess SWAT’s (Soil and Water Assessment Tool) capability to reproduce past and current hydrological flows using corrected and uncorrected RCM projections and (2) to evaluate the impacts of future climate projections on key hydrological processes within the Pike River watershed, in order to inform the development of watershed management plans that incorporate climate change considerations.

Tools and methods

Description of the study site

The Pike River watershed spans the Québec/Vermont border and drains into the Missisquoi Bay located in the northern portion of Lake Champlain (Figure 1). Eighty-four percent of the watershed’s 630-km² area lies in the Eastern Townships/Montérégie region of southern Québec. The town of Bedford, which is located at the centre of the watershed, separates the geophysically distinct upstream (390 km²) and downstream (240 km²) portions of the watershed.

Figure 1. Maps of the Pike River watershed indicating the location of meteorological and hydrological stations. The land-use data are derived from a landsat image taken in July 1999 (Cattaï 2004).

The Pike River watershed is believed to be an important contributor to the recurrent contamination of Missisquoi Bay by toxic cyanobacteria, and has been the focus of several multi-scale hydrological studies. Detailed descriptions of the watershed can be found in Deslandes et al. (2007) and Michaud et al. (2007). The downstream portion of the watershed has a flat topography and fertile soils, and is mostly (75%) devoted to intensive agricultural activities. Fifty percent of the surface area under crop cultivation is devoted to annual crops, particularly grain and silage corn, and soybean production. Intensive animal operations predominate the rugged, steeper and mostly forested upstream portion of the watershed, although some annual crops are also cultivated.

The hydrological regime of the 65-km-long Pike River is characteristic of snowmelt basins, with one major peak flow during the snowmelt season (mostly in April), and a lower discharge during the summer. The mean annual discharge of 8.9 m³ s⁻¹ at the outlet can be 3 to 4 times higher during spring (OBVBM 2011). The climate is temperate, with short humid summers (mean daily maximum temperature of 20°C) and cold dry winters (mean daily maximum temperature of –10°C). Snowfall accounts for approximately 20% of the annual precipitation. Climate normals (1971–2000) for annual rainfall at the Farnham and Sutton weather stations (Figure 1) are respectively 1095 and 1272 mm. Similarly, annual normals for snowfall at the two weather stations are 204 and 281 cm. One half of the precipitation falls during the growing season, which extends from May to September (Environment Canada 2012).

Hydrological modelling

The SWAT model (Arnold et al. 1998) was chosen for use in the present study because of its sound conceptual and physical foundations as well as its ability to simulate the hydrologic response of large, heterogeneous and agriculturally intensive watersheds (Borah and Bera 2003, 2004; Deslandes et al. 2007). This model has been successfully applied to simulate different water quality processes within the Pike River watershed (Deslandes et al. 2007; Eastman 2007; Gollamudi et al. 2007). Some studies have equally used SWAT for climate change studies (Douglas-Mankin et al. 2010). The superiority of SWAT in simulating best management practices (BMP) makes it a candidate of choice for the assessment of climate change adaptation strategies (Borah and Bera 2004) in intensively cultivated watersheds like the Pike River watershed.

Although Deslandes et al. (2007) used SWAT 2000 (Arnold et al. 1998) to successfully simulate the hydrology and phosphorus (P) loss within the Pike River watershed (Deslandes et al. 2007), the model was not adapted to Québec soils and drainage conditions. In the present study, a modified version of SWAT 2005, known as SWATqc, with customized algorithms to improve subsurface tile drainage flow and surface runoff predictions in Southern Québec (Michaud et al. 2008), was used to study the impact of climate change on hydrological processes in the Pike River watershed of southern Québec.

SWAT input data

Data for the spatial parameterization of the model was retrieved from Deslandes et al. (2007) and transferred from SWAT2000 to a SWAT2005 platform (Arnold et al. 1998; http://swatmodel.tamu.edu/contact/). Spatial data were entered into SWAT through a graphical interface (ArcSWAT) linked to Arcview 3.3 (ESRI, Redlands CA, http://www.esri.com/). The data consisted of a digital elevation model (DEM) with a spatial resolution of 30 ± 1.3 m (Deslandes et al. 2002), a digital hydrographic network with 99 sub-watershed outlets, a land-use layer derived from the classification of Landsat 7 ETM+ (Enhanced Thematic Mapper Plus) imagery acquired in July 1999 (Cattaï 2004), and digital soil layers and soil physical properties extracted or estimated from various soil surveys and studies conducted in this region (Thériault et al. 1943; Cann et al. 1948; Bernard 1996; United States Department of Agriculture – Natural Resources Conservation Service [USDA-NRCS] 1999; United States Department of Agriculture – Soil Conservation Service [USDA-SCS] 2006), and a soil degradation inventory (Tabi et al. 1990).

Daily precipitation and temperatures from 1970 to 2006 for the Farnham, Phillipsburg and Sutton meteorological stations closest to the watershed (Figure 1) were extracted from the National Climate Archives Database (Environment Canada 2007). The precipitation and temperature data were entered into SWAT via the ArcSWAT interface.

Information pertaining to initial soil nutrient content and field management operations (such as tillage, sowing, fertilizing, harvest and the schedule of these operations) were obtained from Deslandes et al. (2007). But for the scheduling of management operations, which was updated to correspond to the calibration and validation periods used in this study, all management parameters were maintained as reported in Deslandes et al. (2007). Probable schedules of field management were determined by overlaying information from reports on the state of the crops (La Financière agricole 2008), rainfall events and observed management practices. These reports contained general information about the evolution/state of field work, and crop growth in the fields throughout the year. The 1998–2006 schedules were applied to the period from 1971 to 1998 due to a lack of data for the years prior to 1998.

SWAT set-up and calibration

The Hydrological Response Unit (HRU; a unique combination of soil and land use type within a sub-watershed) is the calculation unit in the SWAT model. Using the automatic delineation option of SWAT, 99 sub-watersheds were delineated for the study watershed based on the input DEM and digital hydrological network data. Overlay of the land-use and soil digital maps resulted in the creation of 1872 HRUs.

Runoff and evapotranspiration were estimated by the Soil Conservation Service (SCS) curve number (USDA-SCS 2006) and Hargreaves (Hargreaves et al. 1985) methods, respectively, while flow routing through the hydrological network was performed using the Muskingum method (Overton 1966).

Model calibration and validation were performed using daily streamflow data from the four Centre d’Expertise Hydrique du Québec (CEHQ 2013) monitoring gauges: Pike River upstream (PRup), Pike River downstream (PRdw), Walbridge Creek upstream (WCup) and Walbridge Creek downstream (WCdw) (Figure 1). The available observed data were split into two data sets spanning the periods from November 2001 to May 2003 and November 2004 to May 2006, and used for SWAT calibration and validation, respectively.

Climate projections and impact assessment

The validated SWAT model was run with four climate projections provided by the Ouranos Consortium – a Quebec-based consortium of experts on climate science and modelling which provides regional climate data sets (projections) for numerous research groups (Ouranos 2013).

Climate impact assessments nowadays are usually performed using a large number of plausible climate projections in order to capture as many sources of uncertainty as possible. This approach, however, is time consuming and requires large computational resources. To circumvent these constraints, only four climate change projections were used in this study. Figure 2 shows how the four climate projections used in the present study compare to an ensemble of 140 projections from other global and regional climate models (including the well-known CSIRO (Commonwealth Scientific and Industrial Research Organisation), ECHAM (European Centre Hamburg Model) and HadCM3 (Hadley Centre Coupled Model, version 3) models driven by various emission scenarios (A2, B1 and A1B). The four projections were purposely chosen to cover a spectrum of temperature and precipitation changes (low, medium and high increases as depicted in Figure 2) and uncertainties that may arise from modelling techniques. Detailed descriptions of the climate models and attendant data sets are presented in the sections that follow.

Figure 2. Average changes in temperatures and precipitation in Missisquoi Bay watershed between 1971–2000 and 2041–2070. Projections ADC, ACU, AFD and ARP were chosen for the study. ADC, ACU, AFD and ARP are the four projections chosen for the study. The name of the projections were arbitrarily given by the consortium Ouranos provider of the climate data.

Climate data sets and impact assessment

Ouranos produced four pairs of climate data sets for the periods 1971–2000 and 2041–2070 using the French Arpège model – a global model with a variable resolution grid (Douville et al. 2002; Gibelin and Déqué 2003) – and the fourth-generation Canadian Regional Climate Model (CRCM4) (Caya and Laprise 1999; Plummer et al. 2006) driven by the third version of the Coupled Global Climate Model (CGCM3) (Scinocca and McFarlane 2004; McFarlane et al. 2005).

RCMs are subject to systematic errors (biases) transferred from the GCMs driving them (Déqué and Piedelievre 1995), or due to the RCMs themselves (Jha et al. 2004; Gagnon et al. 2009). Other uncertainties come from the connection between them (driving method) and the greenhouse gas (GHG) emission scenarios. Generally, bias correction techniques are applied to resolve these issues (Graham et al. 2007; Boé et al. 2009; Prudhomme and Davies 2009). At the time of the study, there was no widely recognized method for bias correction capable of preserving changes in the distribution of rainfall intensity and frequency. The simple delta correction described in the proceeding sections was therefore used to correct RCM biases.

For the future (2041–2070) climate simulations, the climate models were supplied with the Intergovernmental Panel on Climate Change (IPCC) Special Report on Emissions Scenarios (SRES) A2 (Nakicenovic et al. 2000), while measured GHG emission concentrations were used for the historical (1971–2000) period. The A2 scenario represents a world development scenario that produces one of the highest continual increases in GHG emissions. Other SRES scenarios present only minor differences in the global surface warming for the time horizon chosen in this study (Figure 10.4 in IPCC 2007); hence, it was assumed that the inclusion of all SRES scenarios would not result in significant improvements in the capture of uncertainties. In our view, the use of this worst-case scenario (A2) should permit the development of robust WMPs to withstand future climate-induced impacts. The characteristics of the four climate data sets, ADC, ACU, AFA/AFD and ARP (Ouranos 2010), are presented in Table 1. All simulations have distinct characteristics and equal probability of occurrence. Simulations ADC and ACU differed in their initial starting date only. AFD was produced from a different version of CRCM4 while ARP was generated with a completely different model (Arpège). Before comparing the future hydrological simulation to the historic, the combination of SWAT-RCM climate was validated by comparing the resulting historical outputs to a reference simulation driven by observed climate only. A supplemental data set known as ERA40 (Uppala et al. 2005) was also used to assess the influence of the CRCM4 on the hydrological results. The ERA40 is a reanalysis of observed data spanning the period 1971–2000. The code name for the CRCM4 climate data set derived from ERA40 is ACZ.

Table 1. Characteristics of the climate data sets driving SWAT (Soil and Water Assessment Tool). CRCM: Canadian Regional Climate Model; CGCM3: Coupled Global Climate Model version 3; SRES: Special Report on Emissions Scenarios; AMNO: Domain of North America for Climate Model.

Name of the climate simulations	Time window	Regional climate model (RCM)	Driver of the RCMs	GCM member (expressing initial conditions)	RCM domain (cells × cells)	GHG emissions scenario
ACZ	1971–2000	CRCM 4.1.1	Reanalysis ERA40	None since driven with reanalysis	Québec (112 × 88)	Historical measurements
ADC	1971–2000 and 2041–2070	CRCM 4.1.1	CGCM3 3.7.1	1	Québec (112 × 88)	SRES A2
ACU	1971–2000 and 2041–2070	CRCM 4.1.1	CGCM3 3.7.1	4	Québec (112 × 88)	SRES A2
AFA/AFD	1971–2000 and 2041–2070	CRCM 4.2.3	CGCM3 3.7.1	5	AMNO (182 × 174)	SRES A2
ARP	1971–2000 and 2041–2070	Arpège 4.4	No driver Arpège Variable grid mesh	1	T 159 (182 × 174)	SRES A2

The four data sets were chosen because of their distinct characteristics (e.g. different domain, model version or type, and different initial conditions). The use of simulations produced with existing techniques resulted in different climate outcomes and allowed for the inclusion of the main sources of uncertainty. Furthermore, the four projections (ADC, ACU, AFA/AFD and ARP) used in this project were chosen from over 140 simulations of temperature and precipitation produced with 20 recognized GCMs and the CRCM run with different parameterizations and greenhouse gas emission (GHGE) scenarios to account for various uncertainties mentioned previously.

Daily precipitation, maximum and minimum temperatures (at 2 m above ground) and the maximum half-hour rainfall for the historical and future periods were transferred from three CRCM grid points into SWAT. The grid points were chosen such that their centres coincided with the locations of the meteorological stations used to calibrate the modelled hydrology. Grid cells were therefore centred over the geographical coordinates for Farnham (–73.099, 45.3224), Sutton (–72.7086, 45.0749) and Philipsburg (–73.0591, 44.8004) (Figure 3).

Figure 3. Map of the Canadian Regional Climate Model (CRCM) grid cells providing the simulated climate data. Stars represent the meteorological stations from which observed precipitation and temperature data were obtained.

Climate data correction

The so-called perturbation or delta change method was used to correct the climate series whenever historical simulations produced with RCM data were not representative of current and past hydrology, during the validation step. To at least include changes in precipitation intensity, the maximum half-hour rainfall values produced by the RCM for historical and future simulations were used.

In the delta correction method, mean monthly changes in temperature and precipitation emanating from use of historical and future climate data sets were calculated and applied to the observed daily series (climate data used in SWAT1, long-term simulation [1971–2000] using observed climate data). Briefly, to correct for temperature, temperature changes (deltas) were calculated by averaging the daily temperature for each month for the historical and future periods and then subtracting the mean monthly temperature for the historical period from that of the future period. These deltas were subsequently added to the daily temperature for the corresponding month and used in running SWAT. Correction for precipitation was performed by multiplying the respective precipitation deltas (a ratio of the mean future monthly precipitation value to mean historical monthly value) with the corresponding observed daily values.

Summary of SWAT simulations

Figure 4 summarizes the various SWAT simulations together with their climate input data sets. SWAT1 is a long-term (1971–2000) simulation obtained with the observed climate data and it corresponds to the calibrated reference simulation. SWAT2 is the simulation obtained with the reanalysis simulation ACZ. The output for the historical period simulations with the climate model data sets was defined as SWAT3 and those for the future period (2041–2070) simulations using the delta corrected climate inputs were defined as SWAT4.

Figure 4. Schematic diagram of the simulations performed in this study (adapted from Jha et al. 2004). Members refer to an initial condition (i) defined by the starting date and the parameterization of the CGCM. X is the name ADC, ACU, AFA/AFD or ARP given by Ouranos. Members used for this project are described in Scinocca and McFarlane (2004). RCM: Regional Climate Model; GCM: Global Climate Model; SWAT: Soil and Water Assessment Tool. SWAT1 represents the hydrological results obtained with the calibrated simulation of reference. SWAT2 represents results obtained with the simulation ran with reanalysis climate data. SWAT3 represents results with simulations ran with climate model data for the reference period (1971-2000). SWAT 4 represents results obtained with simulations ran with climate model data for the future period (2041-2070).

SWAT1 was only used to determine biases arising from SWAT, by comparing SWAT1 output with observed streamflow. SWAT2 was employed to validate the use of SWAT with CRCM data and to potentially highlight sources of error in the resulting SWAT and CRCM predictions. Because reanalysis ERA40 are close to observed data, the influence (or bias) of CRCM4 on the hydrological results can be assessed by comparing SWAT2 to SWAT1. SWAT3 were used to validate historical simulations while SWAT4 were compared to SWAT 1 and used to assess the impact of climate change. Paired t-tests (or the non-parametric Wilcoxon test when data were not normally distributed) were applied to streamflow and water yield (defined as total runoff) to determine if mean changes between the future (2041–2070) and the historical (1971–2000) periods were significant. The threshold of significance was 5%. Non-significant changes were attributed to the natural variability in climate over the 30-year span.

Results and discussions

SWAT performance with observed climate inputs

The coefficients of performance for streamflow computed on a daily, monthly and annual basis for the calibration period (November 2001 to May 2003) were all above the criteria of acceptability proposed by Moriasi et al. (2007). Whereas monthly and yearly R² values exceeded 0.7, Nash Sutcliff efficiency (NSE) values were always above 0.5, with the absolute percentage deviation remaining below 25%. The coefficients of performance for streamflow for the validation period (November 2004 to May 2006) did not meet acceptability criteria, probably due to SWAT’s inability to effectively capture high-intensity, short-duration rainfall events at the daily time step. Validation of the model performance over a longer validation period (1980–2000) nonetheless resulted in relatively good coefficients (R², NSE, and % deviation coefficients of 0.6, 0.55 and 14%, respectively). This longer period was used in a second validation aimed at testing the robustness of the calibration parameters to long-term changes in climatic conditions, taking into account natural climatic variability.

The calibrated SWAT model was able to simulate components of the water budget consistent with those from other studies (den Harog 1975; Michaud et al. 2008). Twenty percent of annual precipitation, for example, occurred in the form of snow in agreement with data recorded at meteorological stations in the study area (Environment Canada 2012). Fifty percent of precipitation left the watershed through evapotranspiration while surface runoff accounted for 30% of the water yield (total runoff) (Table 2). These values were also reported on the map created by den Harog (1975). For the study, the averaged water balance error computed with Equation (1) was 3%.(1)

Table 2. Simulation of the water cycle components over the entire Pike River watershed between 1971 and 2000 for the different climate data sets. Values are 30–year annual averages in mm/yr.

	SWAT1	SWAT2	SWAT3-ADC	SWAT3-ACU	SWAT3-AFA	SWAT3-ARP
Precipitation	1147	1239	1239	1227	1267	1476
Snowfall	256	257	202	211	300	310
% of PCP	22%	21%	16%	17%	24%	21%
Snowmelt	241	245	192	202	283	295
% of PCP	21%	20%	15%	16%	22%	20%
Water yield	522	550	520	516	820	834
% of PCP	46%	44%	42%	42%	65%	57%
Evapotranspiration	597	647	662	660	417	594
% of PCP	52%	52%	53%	54%	33%	40%
Deep aquifer	68	68	79	76	112	123
% of PCP	6%	5%	6%	6%	9%	8%
Surface runoff	167	212	145	155	218	206
% of WYLD	32%	39%	28%	30%	27%	25%
Subsurface flow	354	336	375	362	502	628
% of WYLD	68%	61%	72%	70%	61%	75%

SWAT1, SWAT2, SWAT3-ADC, SWAT3-ACU, SWAT3-AFD and SWAT3-ARP are the outputs of the model SWAT run with the observed climate and climate data sets (or projections) described in section entitled Climate data sets and impact assessment. PCP: precipitation; WYLD: water yield (which is total runoff). SWAT: Soil and Water Assessment Tool.

SWAT Performance with climate model inputs for the reference period (1971–2000)

Average water cycle components

Table 2 shows annual averages for the main water cycle components for the historical period (1971–2000) simulated using observed data as well as the various climate data sets presented in Table 1.

The ACU and ARP simulations predicted the lowest (1227 mm) and highest (1476 mm) annual precipitation, respectively, in the Pike River watershed (Table 2). The historic (1971–2000) precipitation for the study watershed derived from the SWAT3 simulations driven by the four climate model data sets (ADC, ACU, AFA and ARP) were between 80 mm (7%) and 329 mm (29%) higher than the SWAT1 simulation. This precipitation overestimation by the four climate projections during the reference period is likely attributable to less-than-perfect parameterization of the climate model’s physical processes and its attendant biases. The bias might also have emanated from the spatial resolution of the RCM (45 × 45-km grid) used to model processes at the local scale (climate at a weather monitoring station).

The biases from the four simulated climates affected absolute values of the water cycle components, but not necessarily the overall hydrological processes. In fact, the proportional distribution of the watershed moisture through the water cycle was rather well simulated with SWAT2, SWAT3-ADC and SWAT3-ACU, given their close agreement with values obtained from SWAT1. Annual surface runoff was slightly overestimated by SWAT2, with 39% of the water yield (total runoff) leaving the watershed as surface runoff compared to 32% from SWAT1. The fact that other components of the water cycle were satisfactorily reproduced suggests that the overestimation was probably caused by differences in the distribution of precipitation events which affected soil moisture and water infiltration.

SWAT3-AFA and SWAT3-ARP, on the contrary, presented a relatively high water yield (60% of the total water balance) and low evapotranspiration (40% of the water balance) (Table 2). This imperfect allocation of precipitation to the water cycle components may be attributed to (1) a negatively biased temperature, especially during the summer months and/or (2) a positively biased summer precipitation coupled with insufficiently high temperatures to evaporate the resulting surplus moisture. A comparison of CRCM simulated precipitations and temperatures to measured data showed that whereas temperatures were, on average, under-predicted by 2°C throughout the year, a 60-mm average over-prediction of precipitation was observed during the May–August period. These results suggest that RCM biases were probably transferred to hydrological results. CRCM bias alone (SWAT2) did not seem to significantly impact annual hydrological processes. Nevertheless, the fact that temperature and precipitation data generated by the AFA and ARP simulations were very different from observed data suggests that the GCMs that drive the RCM may sometimes negatively affect hydrological simulations.

Annual streamflow

Figure 5 depicts annual streamflows simulated with SWAT2 (reanalysis or ACZ) and SWAT1 simulations together with the streamflow data recorded at station PRup during the period 1980–2000 (the longest available streamflow observation data). Observed flows were, on average, higher than SWAT simulated flows. The major biases in the simulation data were probably introduced by the inability of SWAT to accurately reproduce the wettest years.

Figure 5. Long-term comparison (1980–2000) of observed and simulated annual streamflow with observed climate data (SWAT1) and reanalysis data (SWAT2) at the hydrological station PRup. PRup: Pike River upstream station; SWAT: Soil and Water Assessment Tool.

Annual streamflows predicted by SWAT2 (CRCM with reanalysis) did not necessarily match observed flows or those predicted with observed climate data (SWAT1). This temporal mismatch may be explained by the fact that the spatial resolution of RCMs does not allow for the satisfactory capture of local climatic events. Summer rainfall, for example, can vary significantly between points within the same watershed; hence, there are bound to be discrepancies between a station (a point) and the 45 × 45-km tile simulated by RCMs. Furthermore, sequences of precipitation events and temperature from climate runs produced with a set of varying conditions (such as our four climate runs) may differ substantially from observed climate data (de Elia et al. 2009). Comparison of SWAT2 and SWAT3 simulated flows to observed flows produced R² and NSE values (Table 3) below the thresholds of confidence (R² and NSE < 0.5), probably because of the sensitivity of these two statistical indices to the magnitude and the order in which events occur (Moriasi et al. 2007). The NSE statistical index was thus deemed unsuitable for the assessment of the performance of these simulations.

Table 3. SWAT coefficients of performance for monthly streamflow when simulated using RCM historical runs.

		SWAT1	SWAT2	SWAT3-ADC	SWAT3-ACU	SWAT3-AFA	SWAT3-ARP
Comparison versus observed streamflow
Monthly	R^2	0.78	0.67	0.37	0.27	0.22	0.39
	NSE	0.55	0.26	–0.11	–0.28	–1.37	–1.02
Comparison versus SWAT1 streamflow
Monthly	R^2	1.00	0.84	0.62	0.49	0.40	0.56
	NSE	1.00	0.65	0.33	0.09	–1.00	–0.67

A visual inspection of the hydrographs, however, suggested that SWAT2 and SWAT1 produced similar annual flow averages and interannual variability over the past 30 years. To validate this, return periods for annual and highest monthly spring flows were analyzed for each simulation using the Gumbel probability function developed for extreme values (Bedient et al. 2008). The analysis focused on high flows because of their primary importance in determining flood occurrences. The difference in magnitude of flow predicted by SWAT1 and SWAT2 simulations for each return period fell within ± 10% (Figure 6a). SWAT3-ADC and SWAT3-ACU showed fairly close agreement with SWAT1. ADC over- or underestimated streamflows by 17% and ACU by 4%, whereas SWAT3-AFA and SWAT3-ARP systematically overestimated the magnitude of annual flows by about 50% for most return periods. These results are consistent with the water budget which showed similar trends in the overestimation of precipitation and water yield. The discrepancies recorded between simulations may thus emanate from biases generated by the parameterization of global climate models, or the use of different RCM versions (CRCM4.2.3 and Arpège 4.4 in SWAT3-AFA and SWAT3-ARP simulations, respectively).

Figure 6. Comparison of (a) annual and (b) monthly streamflow for diverse return periods for the historical simulations 1971–2000. SWAT1 is the calibrated simulation of reference. SWAT2, the simulation obtained with reanalysis climate and SWAT3-ADC, ACU, AFA or ARP are the historical simulations obtained with the 4 climate model projections.

Seasonal and monthly streamflows

SWAT1 and SWAT2 performed similarly in predicting historical average monthly streamflow at PRup and the outlet of the Pike River during most of the year, but not for the spring snowmelt months (March, April and May) (Figure 7). Furthermore, SWAT2 simulated a delayed flood period, with larger spring flood peaks (Figure 7). This phenomenon is probably explained by the fact that the mean daily temperatures for March produced by the CRCM and used to drive SWAT were too low (–3°C on average) compared to the snowmelt temperature threshold (–1°C) established during SWAT calibration. This therefore led to the underestimation of March snowmelt by SWAT2 and, consequently, an increase in snowpack and shifting of the snowmelt runoff to April and May. This probably accounted for the overestimation of simulated streamflows for April and May observed in this present study. This temperature bias may be due to the spatial resolution of the CRCM. It is indeed conceivable for temperatures simulated for one tile (45 × 45 km) to not accurately represent the exact local conditions at the weather stations or over the watershed. In order to initiate melting, the entire grid of the CRCM must receive enough energy to achieve a temperature that starts the melting (one tile equals one homogeneous climate) which takes more time than supplying energy to a smaller grid.

Figure 7. Comparison of monthly observed and simulated mean streamflow for the period 1980–2000, at the PRup hydrological station and the Pike River outlet. SWAT1 PRup is the simulation performed with observed climate data; SWAT2 PRup is the simulation performed with reanalysis climate data for the PRup gauge. PRup: Pike River upstream gauge; SWAT: Soil and Water Assessment Tool.

It is also important to note that the observed winter streamflow values shown in Figure 7 are atypically high, probably due to interferences resulting from the presence of ice on the river during measurements. This therefore suggests that the magnitude of underestimation for the simulated winter streamflows might, in reality, not be as high as it appears in Figure 7.

There were equally notable discrepancies between SWAT1 and the four simulations driven by the climate model data (SWAT3) (Figure 8). The SWAT3 simulations also predicted generally higher delayed snowmelt and summer streamflows values than those from SWAT1. This overestimation was, however, attributable to biases introduced by the CGCM3. Comparison of observed and simulated precipitation series, for example, showed that all the climate simulations used in this study over-predicted summer monthly precipitation by 25 mm (20%) to 75 mm (75%). This resulted in similar overestimation of spring and summer flows. The fact that SWAT2 did not over-predict flows during these periods suggest that systematic biases in the GCMs were solely responsible for overestimating summer flows while CRCM4 biases were responsible for the spring flood delay (Figures 7 and 8). It should be noted that, inasmuch as combination of the CRCM and CGCM can potentially reduce the impact of the bias on hydrological results (e.g. ADC and ACU peak flows), it might well amplify the bias in some instances (e.g. AFA and ARP peak flows).

Figure 8. Mean annual streamflow simulated with observed climate and the five datasets produced by the RCMs for the historical period (1971–2000). SWAT1 was produced with observed climate data, and SWAT2 with reanalyzed data of observed climate (Table 1). SWAT3 are historical simulations produced with each of the four climate projections ADC, ACU, AFA/AFD and ARP. RCM: Regional Climate Model; SWAT: Soil and Water Assessment Tool.

In order to have a rough estimate of spring flood frequency, the peak flow magnitudes for the March and April flows for various return periods were calculated using the Gumbel probability method. The differences in the magnitude of spring peak flows simulated by SWAT1 and SWAT2 for 1.1, 1.5, 2, 5, 10, 20, 50 and 100-year return periods, just like annual flows, did not exceed + 10% (Figure 6b).

It can be concluded from these results that the SWAT model driven by the CRCM was able to reproduce current interannual streamflow variability, flood frequency and average monthly and annual streamflows over a 30-year period, albeit with some biases in summer and spring which were probably introduced by the climate models. The findings from this study highlight the potential for the transfer of systematic biases from climate models into SWAT, especially in the study of a meso-scale watershed such as the Pike River (630 km²). Although the use of RCMs compared to GCMs results in improvements in the resolution of data needed for hydrological studies, results from our study showed that further downscaling or bias correction remains essential given that biases affected winter temperatures, which in turn affected winter snowpack, and later spring snowmelt and streamflow, significantly. Summer streamflow, which, unlike snowfall and snowmelt, is not governed by temperature thresholds, seems to be affected more by precipitation (precipitation overestimation similar to streamflow overestimation) than temperature, although very low temperature could potentially induce low evapotranspiration and lead to streamflow overestimation.

Impacts of climate change on the study area

Temperature and precipitation

The four climate projections chosen in this study predicted the greatest changes in precipitation and temperature in winter (December, February and March). Projection AFA predicted the greatest changes in temperature and precipitation for all seasons except spring. Projection ARP had the smallest changes in temperature and precipitation for all seasons, while ADC projected the greatest increases in temperature and precipitation for spring only. The average changes in temperatures and precipitation for the Missisquoi Bay based on ACU, ADC, AFA/AFD and ARP projections are presented in Table 4.

Table 4. Average changes in temperature and precipitation for the Missisquoi Bay based on ACU, ADC, AFA/AFD and ARP projections. ACU, ADC, AFA/AFD and ARP are the names of the climate projections.

Seasons	Average temperature changes	Average precipitation changes
Winter	+1.2°C to +4.2°C	+6% to +30%
Spring	+1.7°C to +2.9°C	12% to 21%
Summer	+1.7°C to +3.1°C	–7% to +1%
Fall	+2.2°C to +3°C	+1% to +12%

Watershed hydrology

Improvements in climate change results due to the delta correction

This section compares climate change impacts resulting from the direct use of the climate data set to those resulting from use of the bias corrected (delta method) climate data set. Figure 9 shows the range of future monthly streamflow predicted using the two methods and the four climate projections.

Figure 9. Comparison of streamflow changes simulated using corrected (delta method) versus uncorrected RCM data for (a) November to April months and (b) May to October months. SWAT: Soil and Water Assessment Tool; RCM: Regional Climate Model.

Figure 10. Mean annual streamflow simulated with RCM bias (delta) corrected climate data for the 2041–2070 period. SWAT1 was produced with observed climate data, while SWAT4 are the future simulations produced with each of the four climate projections ADC, ACU, AFA/AFD and ARP. RCM: Regional Climate Model; SWAT: Soil and Water Assessment Tool; SWAT1 is the calibrated simulation of reference; SWAT4 are the projected simulations in the future. ADC, ACU, AFD, and ARP are the name of the simulations.

The main difference between the two methods occurred in March, where streamflow predictions emanating from simulations with bias-corrected data showed a much smaller increase compared to the direct entry of uncorrected data into SWAT (Figure 9). This is due to the fact that negative biases within uncorrected data resulted in very low late-winter and early-spring temperatures which limited the melting of the snowpack compared to the simulations with the delta-corrected data. As a result, when temperatures rose above 0°C in spring, all the accumulated snowpack predicted by simulations with uncorrected data melted at once, producing higher flows in March. The predicted range of monthly minimum and maximum streamflows based on the bias-corrected data was also much narrower, indicating fewer uncertainties within the results. The fact that the use of delta-corrected climate data resulted in the prediction of slightly higher summer streamflows than the uncorrected climate data suggests a smaller climate-induced impact on water availability during this season.

Annual changes in the hydrologic cycle components

Except for snowfall and snowmelt which decreased up to 15% or 39 mm in snow water equivalent, all the components of the hydrologic cycle showed annual increases under future climatic conditions (Table 5). Annual precipitation increases of 7 to 12% led to a streamflow increase of 9 to 19%. Surface runoff increased by 11 to 21%, with the highest increase, surprisingly, predicted by the ARP simulation, which projected the smallest increase in precipitation. Interestingly, the ARP simulation equally predicted the smallest increases in total runoff (a surrogate for streamflow in this analysis), groundwater flow and recharge compared to the other three simulations. This leads us to suggest that ARP future precipitation must have been of higher intensity, leading to increased runoff and less soil infiltration. Another possible explanation for this observation could be that the lower spring temperature projected by ARP caused the soils to remain frozen for longer periods, thereby promoting runoff instead of soils and aquifer recharge. The other three simulations, by contrast, showed greater increases in groundwater flow than surface runoff. The fact that a large part of the watershed under cultivation is artificially drained, with tile drains accounting for 28 to 35% increase in water export, probably explains the greater increase in groundwater flow prediction by ADC, ACU and AFD. Evapotranspiration on an annual scale seemed to barely affect the hydrological cycle, given that all four simulations projected only small increases (2 to 7%).

Table 5. Relative changes witnessed by the water cycle components (SWAT1 vs. SWAT4). SWAT: Soil and Water Assessment Tool; SWAT1 is the calibrated simulation of reference and SWAT4 are the projected simulations in the future. ADC, ACU, AFD, and ARP are the name of the simulations.

Water cycle component	ADC	ACU	AFD	ARP
Precipitation	12%	11%	12%	7%
Snowfall	–12%	–11%	–15%	–1%
Snowmelt	–13%	–12%	–15%	–2%
Surface runoff	14%	14%	11%	21%
Lateral runoff in soil	24%	18%	24%	12%
Tile	35%	28%	31%	10%
Groundwater	17%	13%	18%	6%
Total groundwater	23%	18%	22%	8%
Deep aquifer recharge	17%	14%	19%	7%
Water yield (total runoff)	20%	16%	19%	11%
Evapotranspiration	6%	6%	7%	2%
Total recharge	17%	13%	18%	6%

Overall, the proportions of simulated components of the hydrologic cycle in the future climate remained comparable to the historical simulations, with 50% of precipitation still leaving the watershed by evapotranspiration and 45 to 50% via water yield (total runoff) while losses to the deep aquifer accounted for 6%. The projected proportions of surface to subsurface runoff were similar for all simulations but ARP.

Seasonal changes in streamflow

The majority of the four simulations typically projected significant monthly streamflow increases between September and March (Table 6). The greatest increases were generally projected in January or February (54 to 375% increase) while the smallest significant changes occurred during fall (decreases of 5 to 7%).

Table 6. Changes (%) in mean monthly streamflow for delta corrected simulations (Corr) vs. simulations driven with direct entry of RCM data (Raw). Bold values represent significant changes at p < 0.05 (paired t-test). Values in italics show changes. RCM: Regional Climate Model; ADC, ACU, AFA-AFD and ARP are the name of the simulations.

	ADC		ACU		AFA-AFD		ARP		Simulations mean
Month	Corr	Raw	Corr	Raw	Corr	Raw	Corr	Raw	Corr	Raw
J	303	158	230	344	375	202	54	128	241	208
F	375	80	220	174	336	495	47	162	244	228
M	22	59	41	67	15	148	16	49	23	81
A	–29	–3	–31	–27	–30	9	1	–17	–22	–10
M	–13	–11	–14	–20	–12	–26	3	–8	–9	–16
J	0	–31	13	0	2	–10	22	18	9	–6
J	–13	–43	–6	–13	–14	0	9	–10	–6	–17
A	–25	–13	–20	–16	–26	–37	–17	–10	–22	–19
S	4	–6	7	–21	4	0	44	11	15	–4
O	11	31	–7	–7	11	6	2	7	4	9
N	31	31	20	20	29	48	–5	–5	19	24
D	101	172	80	113	107	128	30	62	79	119
Year	19	12	14	7	17	10	9	9	15	10

March recorded the highest flows (peak) under future climatic conditions, as opposed to April peaks for historical/current conditions (Figure 10). All but the ACU simulation projected decreases in average peak flows under future conditions. As a consequence of the shift in peak flow towards March, the ADC, ACU and AFD simulations projected significant decreases in April streamflows, while ARP showed no significant changes. The highest absolute decrease in streamflow occurred in April, with a magnitude similar to the February increase. Although decreases in streamflow were also projected for the remainder of the year (May to August), the differences were not statistically significant for all simulations. The significant summer streamflow decreases projected by ADC and ACU ranged between –13 and –26%. The absolute decrease in summer streamflow never exceeded 5 mm because of the already low summer flows. This decrease, coupled with existing river drawdowns in the study watershed, however, suggests that climate change might exacerbate water scarcity; hence the need to develop WMPs that take into account this new reality.

It is interesting to note that whereas ACU projected the greatest increase in March spring flooding, it nevertheless showed temperature increases akin to ADC but slightly lower precipitation in spring. Furthermore, the increases in winter precipitation and temperature were much lower for the ACU than ADC and AFD climate projections, probably because of the higher average spring flood projected by ACU. In fact, the higher winter temperatures projected by ADC and AFD promoted rainfall and runoff, thereby diminishing water storage in the snowpack and reducing spring flood in March and April. ACU, on the contrary, suggested the storage of higher precipitation amounts on the soils during winter because the lower temperatures projected allowed for more snow accumulation, which subsequently melted and contributed to the spring flood at the onset of spring. The greatest effects of climate change on seasonal streamflow therefore did not emanate from the projections with the greatest temperature and precipitation increases.

The overall pattern of changes reported here has also been observed in studies performed in other watersheds in Québec (Minville et al. 2008; Boyer et al. 2010). Minville et al. (2008) summarized these trends as an increase in winter flow, an earlier spring flood and a decrease in summer–fall flow. Boyer et al. (2010), on the other hand, projected ± 10% changes in mean annual discharge for the entire Richelieu River valley depending on the model and GHGE scenarios used. Possible doubling or tripling of winter stream flow and a 40% increase or decrease in spring peak flow were further suggested (Boyer et al. 2010).

The changes observed in our study are mostly attributable not only to higher precipitation but also to higher winter temperatures which increase liquid precipitation and snowmelt episodes, thus reducing water storage in the snowpack. The early onset of spring flooding can be attributed to an increase in the daily mean temperatures during the month of March. These temperatures often exceeded the snowfall/snowmelt temperatures, thereby triggering earlier snowmelt runoff. The general decrease in spring peak flow is a direct result of winter changes (an increase in temperature) which shift spring snowmelt runoff from April to March. The decreases in summer streamflows can be attributed to the higher evapotranspiration, which resulted from an increase in temperature and a decrease in rainfall.

Conclusions

The present study employed the SWAT_qc model linked to regional climate models to assess the impacts of climate change on the hydrology of the Pike River watershed, an important agricultural watershed located in southern Quebec. SWAT simulations with uncorrected RCM outputs successfully reproduced the seasonality of hydrological flows and interannual variability. The transfer of RCM biases into SWAT, however, significantly impacted the model’s ability to predict not only the timing but also the magnitude of spring peak floods, suggesting the need for bias correction of RCM-generated climate data prior to use in hydrological modelling.

SWAT simulations with delta-corrected RCM data showed that the climate change projection for the 2041–2070 horizon could significantly affect the timing of spring floods, winter runoff and streamflow in the Pike River watershed. All four climate-driven simulations predicted a 9–19% increase in future annual streamflows and an 11 to 21% increase in surface runoff.

On a seasonal basis, winter runoff witnessed a 2- to 3-fold increase attributable to the temperatures, rainfall and snowmelt increase and numerous freeze–thaw events projected for this season. A noticeable shift in the peak discharge from April to March was equally envisioned by the ADC (obtained with the climate projection from CRCM4 version 1.1 driven by member 1 of the CGCM3), ACU (also obtained with the climate projections from CRCM4 version 1.1 but driven by member 4 of the CGCM3) and AFD (obtained with climate projections from the CRCM4 version 2.3 run on a larger domain than for the ADC and ACU) SWAT simulations that projected higher winter and spring temperatures and precipitations. The two climate projections (ADC and AFD) projecting the greater increase in winter temperature also projected a decrease in the magnitude of the spring peak flow attributable to more gradual snow melt processes. Projection ACU, with a moderate increase in winter temperature, on the other hand, projected a higher spring flood and consequently higher risks of flooding.

The projected decreases in late spring and summer streamflows suggests that climate change could potentially increase risks of water stress during periods when the resource is most needed for agricultural productivity and recreation in this region. An important caveat to consider when interpreting these results is the fact that changes in precipitation frequency were not incorporated into the analyses; hence, the results may represent rather conservative estimations of the potential climate change impacts on the hydrology of the study watershed.

Acknowledgements

This study was sponsored by a Natural Sciences and Engineering Research Council of Canada (NSERC) grant to C. Madramootoo. The authors also thank the Ouranos Consortium for providing data, sound advice and technical assistance as well as funding for the realization of the project. Special thanks to L’Institut de Recherche et Développement en Agroenvironnement (L’IRDA) for providing the modified SWAT model, data for its parameterization and technical assistance for hydrological modelling.