Assessing climate change impacts on fresh water resources of the Athabasca River Basin, Canada

Science of the Total Environment 601–602 (2017) 425–440

Contents lists available at ScienceDirect

Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv

Assessing climate change impacts on fresh water resources of the Athabasca River Basin, Canada Narayan Kumar Shrestha, Xinzhong Du, Junye Wang ⁎ Athabasca River Basin Research Institute (ARBRI), Athabasca University, 1 University Drive, Athabasca, Alberta T9S 3A3, Canada

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T

• Climate change impact analysis of the Athabasca River Basin (ARB) using the SWAT model • Explicit consideration of both the blue and green water resources • Future (mid- and late century) climate data generated by CanRCM4 for RCP 4.5 and 8.5 • Both the blue and green water resources in the ARB are likely to increase in the future. • Evidences of temporal and spatial heterogeneity of the blue and green water resources

a r t i c l e

i n f o

Article history: Received 25 February 2017 Received in revised form 8 April 2017 Accepted 2 May 2017 Available online xxxx Editor: Dr. Ouyang Wei Keywords: Athabasca River Basin (ARB) Blue and green water resources Climate change SWAT CanRCM4

a b s t r a c t Proper management of blue and green water resources is important for the sustainability of ecosystems and for the socio-economic development of river basins such as the Athabasca River Basin (ARB) in Canada. For this reason, quantifying climate change impacts on these water resources at a finer temporal and spatial scale is often necessary. In this study, we used a Soil and Water Assessment Tool (SWAT) to assess climate change impacts on fresh water resources, focusing explicitly on the impacts to both blue and green water. We used future climate data generated by the Canadian Center for Climate Modelling and Analysis Regional Climate Model (CanRCM4) with a spatial resolution of 0.22° × 0.22° (~25 km) for two emission scenarios (RCP 4.5 and 8.5). Results projected the climate of the ARB to be wetter by 21–34% and warmer by 2–5.4 °C on an annual time scale. Consequently, the annual average blue and green water flow was projected to increase by 16–54% and 11–34%, respectively, depending on the region, future period, and emission scenario. Furthermore, the annual average green water storage at the boreal region was expected to increase by 30%, while the storage was projected to remain fairly stable or decrease in other regions, especially during the summer season. On average, the fresh water resources in the ARB are likely to increase in the future. However, evidence of temporal and spatial heterogeneity could pose many future challenges to water resource planners and managers. Crown Copyright © 2017 Published by Elsevier B.V. All rights reserved.

1. Introduction

⁎ Corresponding author. E-mail address: [email protected] (J. Wang).

http://dx.doi.org/10.1016/j.scitotenv.2017.05.013 0048-9697/Crown Copyright © 2017 Published by Elsevier B.V. All rights reserved.

The Athabasca River Basin (ARB) is ecologically and economically significant to the development and sustainability of northern Alberta communities. Indeed, the multi-billion dollar oil sands industry requires a large amount of water, approximately 4.4% of the average yearly

426

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

streamflow (Sauchyn et al., 2015). Ever increasing industrialization and a growing population are putting immense pressure on the water resources of the basin and pose a notable threat to the environmental and economic sustainability of the ARB. Moreover, the negative impacts of climate change are evident all over the world. The sustained increase in the emission of greenhouse gases will alter all components of the climate system (IPCC, 2014). By the end of this century, the global surface temperature is likely to be 2 °C warmer, relative to the temperature experienced during the period of 1850–1900. Such increases are alarming, especially for the ARB with its water cycle which primarily depends on glacial melt at its headwater and spring freshet at its other sections. Some evidence of these changes have already been reported in the literature (Eum et al., 2017; Kerkhoven and Gan, 2011; Leong and Donner, 2015; Sauchyn et al., 2015). However, most studies to date have been limited to analyzing stream flow trends. This is reflective of a traditional water resource planning and management approach in which greater emphasis is put on the blue liquid water - the stream flow (Falkenmark and Rockström, 2006). Other forms of water resources, such as green water storage (soil water) and green water flow (evapotranspiration), are equally important to consider. Hence, a comprehensive climate change impact assessment should follow the “Blue and Green Water Paradigm” in which both resources are explicitly examined (Falkenmark and Rockström, 2006). Furthermore, the ARB is located in a cold-climate region with many impoundments (e.g., lakes). Thus, the hydrology of the basin is snowmelt dominated, meaning its winter flows are low and its spring flows are high. A climate change impact assessment study in such a cold region should take into account both the blue and green water resources, as both resources are affected by climate change. To our knowledge, there has not been a study which considered both water resources at different spatial and temporal scales to make climate change impact assessments of the basin. With this in mind, this study has been carried out with a primary objective of quantifying the impacts of climate change on monthly, seasonal, and annual water balances of blue and green water resources at sub-basin, regional, and basin-wide spatial scales. We addressed this objective with a two-step approach. First, we built a SWAT model of the basin using a spatial and hydrometeorological data set. We then evaluated the applicability and suitability of the model in the snow-dominated and cold climate in which the ARB basin is located. During the second step, we fed high resolution (~25 × 25 km) future climate data generated by the Canadian Center for Climate Modelling and Analysis Regional Climate Model (CanRCM4) for two emission scenarios (RCP 4.5 and 8.5) into the calibrated and validated model to assess the impacts on blue and green water resources. 2. Materials and methods

variety of processes which are primarily from rural catchments (Arnold et al., 1998). The SWAT is one of the most widely used simulators for hydrologic modelling (Arnold et al., 2012). This is partially because it is an open-source tool containing various add-ons, such as SWAT calibration and uncertainty or sensitivity analysis programs – SWAT CUP (Abbaspour et al., 2007). The simulator has also been recognized by the US Environmental Protection Agency (Abbaspour et al., 2015). Moreover, the simulator has been used to simulate hydrological processes (Leta et al., 2015), water quality processes (Santhi et al., 2001; Shrestha et al., 2017), and erosion sediment transport modelling (Shrestha et al., 2013), in addition to being utilized as a component model in the integrated modelling chain (Shrestha et al., 2014). The SWAT uses spatial datasets for elevation, land-use, and soil, along with several hydro-meteorological datasets which are typically integrated using the Geographic Information System (Winchell et al., 2010). SWAT compartmentalizes a watershed into sub-basins that are further divided into Hydrological Response Units (HRUs) which are unique combinations of soil, land-use, and slope (Arnold et al., 2011). The hydrological component of SWAT considers precipitation, infiltration, deep aquifer, channel transmission and evapotranspiration (ET) losses, surface runoff (Qsurf), and lateral and return flow (Qsub-surf) for water balance calculations (Eq. (1)). SWAT differentiates precipitation as rainfall or snowfall while comparing air temperature with a snowfall temperature parameter (SFTMP). As a result, the model keeps track of the volume and areal extent of snowpack, as well as the corresponding snowmelt as per Eqs. (2)–(4). Snow accumulation and melting are processes that can be spatially varied using elevation bands in a subbasin. A maximum of 10 elevation bands can be defined at each subbasin. The model allows for the consideration of two lapse rates, the temperature (TLAPS) and precipitation (PLAPS), in order to, respectively, vary the temperature and precipitation with elevation (Eqs. (5)–(8)). The model employs either the SCS curve number or the modified GreenAmpt method to determine the infiltration and runoff volumes for each HRU. Infiltrated water percolates through each soil layer, as estimated using a storage routing technique. SWAT offers a variable storage method or Muskingum method to route the streamflow generated as a result of runoff coming from each of the HRUs of the sub-basins (Neitsch et al., 2011). t
 SWt ¼ SW0 þ ∑ Pi −ETi −Q i;seep −Q i;surf −Q i;gw

ð1Þ

t
 SNOt ¼ SNO0 þ ∑ Pi −Ei;sub −SNOMLTi

ð2Þ

i¼1

i¼1

SNOMLTi ¼ bi;mlt  SNOCOVi þ

2.1. The study area The Athabasca River Basin (ARB) originates at the Columbia Icefields, located in the Canadian Rocky Mountains in the province of Alberta. The river flows toward the North-East and drains first into Lake Athabasca (Fig. 1) before eventually reaching the Arctic Ocean. The catchment area of the mouth of the ARB is approximately 161,000 km2. The river crosses through the municipalities of Jasper, Hinton, Whitecourt, Athabasca, and Fort McMurray (AWC, 2011). Forest is the dominant landuse type in the basin at almost 82% coverage, followed by agriculture land with about 10% coverage (Fig. S1). Major industries with activities occurring in the basin include agriculture, pulp mills, coal, and oil sand mining (AWC, 2013). 2.2. The hydrologic simulator - Soil and Water Assessment Tool (SWAT)

bi;mlt ¼

ð3Þ

    SMFMX þ SMFMN SMFMX−SMFMN 2π þ  sin ði−81Þ ð4Þ 2 2 365


 Pband ¼ Pday þ ELband −ELgauge 

! PLAPS dayspcp;yr  1000

k

ð5Þ

Pday ¼ ∑ Pband  FRband

ð6Þ

 
 TLAPS Tmx;mn;avg;band ¼ Tmx;mn;av;day þ ELband −ELgauge  1000

ð7Þ

band¼1

k

Tmx;mn;av;day ¼ ∑ Tmx;mn;avg;band  FRband band¼1

The Soil and Water Assessment Tool (SWAT), developed by the United States Department of Agriculture (USDA), is a semi-distributed hydrological simulator used for continuous long term simulation of a

  Ti;snow þ Ti; max −SMTMP 2

ð8Þ

where, SWt = soil moisture at time t (mm); SW0 = initial soil moisture (mm); i = day counter; Pi = precipitation (mm); ETi =

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

427

Fig. 1. The provinces and territories of Canada (top left), the location of the Athabasca River Basin (ARB) in Alberta (top right), and the Digital Elevation Model (DEM) of the ARB showing the Athabasca Lake and the urban centers of the ARB.

evapotranspiration (mm); Qi,seep = percolation through soil profile (mm); Qi,surf =surface runoff (mm); Qi,gw = ground water return flow (mm); SNOt = snow water equivalent at time “t” (mm); SNO0 = initial snow water equivalent (mm); Ei,sub = water equivalent of snow sublimation (mm); SNOMLTi = water equivalent of snowmelt (mm); bi,mlt = melt factor (mm/°C/day); SNOCOVi = fraction of HRU covered by snow (−); Ti,snow = snowpack temperature (°C); Ti,max = maximum air temperature (°C); SMTMP = snow melt base temperature (°C); SMFMX =

maximum snow melt rate (mm/°C/day); SMFMN = minimum snow melt rate (mm/°C/day); ELband = elevation of elevation band (m); ELgauge = elevation of gauge (m); PLAPS = precipitation lapse rate (mm/km); dayspcp,yr = average number of days of precipitation (−); FRband = fraction of sub-basins that are in each elevation band (−); Tmx,mx,avg,day and Tmx,mx,avg,band = maximum, minimum, and average temperature of the day at each elevation band (°C); and TLAPS = temperature lapse rate (°C/km).

428

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

2.3. Model inputs and built-up A SWAT model requires three spatial datasets including, a Digital Elevation Model (DEM), a land-use map, and a soil map. We used a 3″ (90 m × 90 m), hole-filled SRTM DEM (Jarvis et al., 2008) Global Land Cover Characterization-based land-use map (Loveland et al., 2000) with a resolution of 1 km × 1 km and a soil map at a scale of 1:1 million, provided by the Agriculture and Agri-Food Canada (SLC, 2010). The elevation of the basin (Fig. 1) varies between 207 m and 3669 m. Similarly, 11 different land-use classes (Fig. S1) and 320 different soil types were defined and their relative databases were prepared as required by the SWAT model. We used the DEM to derive the stream network and to delineate sub-basins. To do so, we selected a threshold drainage area of 200 km2. Sub-basin outlets were also defined as required at major river confluences and at streamflow gauging stations. The entire process produced a total of 131 sub-basins. To define the HRUs, we derived a slope map from the DEM and divided it into 4 classes, introducing breaks at 5%, 10%, 15% and 20%. As suggested by Strauch et al. (2015), we used thresholds of 10%, 10% and 5% for land-use, slope, and soil, respectively. This gave us a total of 1370 HRUs. Daily meteorological data were then fed into the model. We used precipitation, and maximum and minimum temperature data from 73 stations, as recorded by GoC (2016) (Fig. 2). Depending on the network, the GoC (2016) uses different types of gauges. For

instance, rainfall is measured by using “weighing gauges” at automated stations and by using either a “standard Canadian Type B rain gauge” or a “heated Met One Tipping Bucket rain gauge” at staffed stations (GoC, 2016). Similarly, snowfall at automated stations is measured using an “acoustic snow sensor (SR-50)” while it is measured at staffed stations using a “snow ruler”, before being converted to an equivalent amount of water (GoC, 2016). As indicated in Fig. 2, there are long gaps in the precipitation and temperature data, especially during winter periods. As a first step, we attempted to fill the missing gaps using the CFSR (2016) data set. An analysis of the CFSR (2016) data and the GoC (2016) data produced a low correlation (mean r2 = 0.52, std. = 0.08, n = 52). Furthermore, the CFSR (2016) dataset tended to overestimate the precipitation, especially during the winter. Hence, we treated these values as missing and allowed the SWAT weather generator (Neitsch et al., 2011) to estimate them. For this purpose, appropriate weather generator files were prepared using the pcpSTAT program (Liersch, 2003) and supplied to the SWAT model database. A similar procedure was adopted for the temperature. Furthermore, relative humidity, solar radiation, and wind speed data which had been recorded at 230 stations by CFSR (2016) (Fig. 2) were also fed into the model. CFSR (2016) provides best estimates of several meteorological forces using a “global, high resolution, coupled atmosphere-ocean-land surface-sea ice system” which has been evaluated (Dile and Srinivasan, 2014) and found to be effective, especially in data scarce regions or in order to

Fig. 2. Location of precipitation and temperature stations operated by GoC (2016) and the wind speed, relative humidity, and solar radiation stations operated by CFSR (2016).

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

complement data with long missing series (Fuka et al., 2014). To estimate the potential evapotranspiration, we made use of the PenmanMonteith method. As several lakes of various sizes are located in the basin, we classified them as reservoirs or as lakes depending on whether or not they are regulated. For instance, the largest lake in the basin, the Lesser Slave Lake, is operated using a self-draining weir and has available a daily time series data of observed outflow though the weir. Therefore, we dedicated this lake as a reservoir. Similarly, three other lakes were implemented as reservoirs while the rest were classified as ponds. When more than one lake was located in a sub-basin, the aggregated properties of the lakes were applied to the model. As the basin receives significant amounts of precipitation in the form of snowfall, we defined an allowed maximum of 10 elevation bands within a sub-basin, a number which is typically adopted in such snow-dominated river basins (Grusson et al., 2015). This was also meant to allow for variations within snow related parameters during model calibration, as well as for temperature and precipitation differences across elevation bands. Indeed, several studies have shown an added advantage of using elevation bands to improve simulation results in similar snow-dominated catchments (Grusson et al., 2015; Pradhanang et al., 2011). 2.4. Model sensitivity, calibration-validation and uncertainty analysis We used the SWAT-CUP (Abbaspour et al., 2007) and its SUFI-2 algorithm (Abbaspour et al., 2004) for model sensitivity, calibration, and validation, as well as for uncertainty analysis. A global sensitivity analysis (Abbaspour et al., 2007) which uses Latin-Hypercube (LH) generated

429

parameters was regressed against our chosen objective function, the Nash Sutcliffe Efficiency (Nash and Sutcliffe, 1970). While SWAT offers a wide range of parameters for sensitive analysis, we selected a limited number of parameters (Table 1) that are known to influence the streamflow in similar catchments (Faramarzi et al., 2015; Grusson et al., 2015; Troin and Caya, 2014). Preceded by a two-year warm-up period, the sensitivity analysis was set to run 100 times during the period of 1990–2005. In order to identify dominant parameters at different parts of the basin, many researchers (Ahmadi et al., 2014; Cao et al., 2006; Wang et al., 2012) have suggested using such multi-site sensitivity analysis. Hence, we carried out the sensitivity analysis at Hilton, Athabasca, and just upstream of Fort McMurray (Fig. 1) to represent the upper, middle, and lower reaches of the basin, respectively. We then considered the most sensitive parameters, as suggested by the sensitivity analysis, for model calibration and validation. We dedicated the years 1980 and 1981 as a two year warm-up period, while a calibration period of 16 years (1990–2005) separated the 16 year total validation period into two smaller durations spanning the years of 1982–1989 and 2006–2013. As a result, we were able to include both the drier and wetter periods homogenously in the calibration and validation periods. The Standardized Precipitation Index (SPI) values (Fig. 3) which were calculated based on a 12-month precipitation sum helped to ensure this desired homogeneity. It should be noted that the SPI values which fell above zero indicated wetter periods, while those below the green line indicated drier periods. Both the calibration and the validation periods included at least one extreme event. These events were classified as either severely dry (SPI values b −1) or very wet (SPI values N 1.5) (WMO, 2012), ensuring a robust calibration and validation of the model.

Table 1 Streamflow related parameters considered for sensitive analysis with their maximum, minimum, and default values. Type

Name

Description

Unit

Max. value

Min. Default value value

Snow

v__SUB_SFTMP().sno v__SUB_SMTMP().sno v__SUB_SMFMX().sno v__SUB_SMFMN().sno v__SUB_TIMP().sno v__PLAPS.sub v__TLAPS.sub v__SNO_SUB.sub v__ALPHA_BF.gw v__GW_DELAY.gw v__GWQMN.gw

°C °C mm/°C-day mm/°C-day – mm/km °C/km mm days days mm

−5 −10 0 0 0 −1000 −10 0 0 0 0

5 10 10 10 1 1000 10 150 1 500 5000

1 0.5 4.5 4.5 1 0 −6 0 0.048 31 1000

v__REVAPMN.gw v__GW_REVAP.gw v__RCHRG_DP.gw r__SOL_Z().sol

Snow fall temperature Snowfall melt base temperature Maximum melt rate for snow during the year Minimum melt rate for snow during the year Snow pack temperature lag factor Precipitation lapse rate Temperature lapse rate Initial snow water content Base flow alpha factor Groundwater delay time Threshold depth of water in the shallow aquifer required for return flow to occur Threshold depth of water in the shallow aquifer for ‘revap’ to occur Groundwater revap. coefficient Deep aquifer percolation fraction Depth from the soil surface to bottom of layer

mm – – mm

0 0.02 0 −10%

1000 0.2 1 10%

r__SOL_BD().sol

Soil bulk density

g/cm3

−10%

10%

r__SOL_AWC().sol

Soil available water storage capacity

−10%

10%

r__SOL_K().sol

Soil conductivity

mm H2O/mm soil mm/h

−10%

10%

v__SOL_ALB().sol

Moist soil albedo

–

0

0.25

v__CANMX.hru v__ESCO.hru v__EPCO.hru v__SURLAG.hru r__OV_N.hru r__SLSUBSN.hru r__CN2.mgt v__CH_N2.rte v__CH_K2.rte v__PND_K.pnd v__RES_K.pnd

Maximum canopy storage Soil evaporation compensation factor Plant uptake compensation factor Surface runoff lag time Manning's n value for overland flow Average slope length SCS runoff curve number for moisture condition II Manning's n value for main channel Effective hydraulic conductivity in the main channel Hydraulic conductivity through bottom of pond Hydraulic conductivity of reservoir bottom

mm – – days – m – – mm/h mm/h mm/h

0 0 0 0.05 −20% −20% −40% 0.01 0 0 0

100 1 1 24 20% 20% 40% 0.3 500 1 1

750 0.02 0.05 soil layer soil layer soil layer soil layer soil layer HRU 0.95 1 4 HRU HRU HRU 0.014 0 0 0

Elevation band

Groundwater

Soil

Soil, vegetation, slope

Stream reach Pond Reservoir

v = parameter value is replaced by given value. r = parameter value is multiplied by 1 ± a given value.

430

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

Fig. 3. Warm-up (1980–1982), calibration (1990–2005), and validation (1983–1989 and 2006–2013) periods with 12-month Standardized Precipitation Index (SPI) values.

For the calibration and validation of the model, the SWAT-CUP was run for 500 simulations using a prior uniform distribution of selected sensitive parameters which were based on the upper and lower limits (Table 1). We chose the Nash Sutcliffe Efficiency (Nash and Sutcliffe, 1970) as an objective function. Further runs were carried out when required as per the “new parameter sets”, recommended by the SUFI-2 algorithm. The new parameter set ensures fast conversion as it narrows down the ranges of parameters. We followed a sequential calibration procedure which started upstream of the basin because daily and monthly streamflow data gathered from 35 streamflow gauging stations were readily available (Environment-Canada, 2016). Environment-Canada (2016) uses “state-of-the-art monitoring equipment and methods” to collect and process the data. As a first step, water level is measured using “data loggers” and the water level measurement is converted to streamflow data using an alreadyestablished water level versus streamflow rating curve. A rigorous quality control procedure is then applied before disseminating the data (Environment-Canada, 2016). Although the streamflow data had several missing periods, we did not fill them. Rather, we treated these values as missing and did not consider them in the calibration and validation processes in SWAT-CUP. A proper calibration is assumed to be achieved when reasonable values of the selected goodness-of-fit statistics have been reached and good values for p- and r-statistics (Abbaspour et al., 2007) have been obtained. The p-statistics indicate the percentage of observations bracketed by the 95% prediction uncertainty band, while the r-factor reflects the width of the band. Ideally, we would prefer to have all observations (p-statistics = 1) bracketed in a very narrow band (r-statistics = 0). It should be noted that higher p- statistics could be obtained with an increased r-factor. 2.5. Model performance evaluation Several goodness-of-fit statistics are available to assess model performance. These statistics are generally complemented by graphical plots to assess the performance of model simulations. Thus, we chose to complement time series plots and scatter plots with three goodness-of-fit statistics: (a) the Percentage of Bias (PBIAS) – Eq. (9); (b) the Nash Sutcliffe Efficiency (NSE) – Eq. (10), and (c) the ratio of

Table 2 Range of values of goodness-of-fit statistics for a particular qualitative rating. Very Good

Good

Satisfactory

Unsatisfactory

PBIAS ≤ ±10 0.75 b NSE ≤ 1.00 0.00 b RSR ≤ 0.50

±10 ≤ PBIAS ≤ ±15 0.65 b NSE ≤ 0.75 0.50 b RSR ≤ 0.60

±15 ≤ PBIAS ≤ ±25 0.50 b NSE ≤ 0.65 0.60 b RSR ≤ 0.70

PBIAS ≥ ±25 NSE ≤ 0.50 RSR N 0.70

Root-Mean-Squared-Error to Standard Deviation (RSR) – Eq. (11). As per the ranges of values of these goodness-of-fit statistics, a qualitative rating (Table 2) was also assigned to the simulation result using four ratings (Very Good, Good, Satisfactory, and Unsatisfactory) as suggested by Moriasi et al. (2007). Finally, an aggregated rating was calculated by assigning each rating a value of 1 (Unsatisfactory) to 4 (Very Good). We believe that assigning such qualitative ratings to model simulation results can better communicate the performance of the model simulations to a wider range of stakeholders. PBIAS ¼

 n
∑i¼1 Q i;sim −Q i;obs

NSE ¼ 1−

n

∑i¼1 Q i;obs

 100

2 n
∑i¼1 Q i;sim −Q i;obs  2 n ∑i¼1 Q i;obs −Q i;obs

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 n
∑i¼1 Q i;sim −Q i;obs RSR ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 n ∑i¼1 Q i;obs −Q i;obs

ð9Þ

ð10Þ

ð11Þ

where, Qi,sim = simulated streamflow at time step i; Qi,obs = observed streamflow at time step i, Q i,obs = mean of observed streamflows; n = number of time steps. 2.6. Future climate data The calibrated and validated SWAT model in the base period (1983– 2013) was then subjected to the future climatic data and run for 500 simulations with the calibrated parameter ranges. To do so, we used daily future climate data generated by a Regional Climate Model (RCM) of the Canadian Center for Climate Modelling and Analysis (CCCMA) - the CanRCM4 (Scinocca et al., 2015). This RCM offered downscaled climate data at a spatial resolution of 0.44° × 0.44° (~ 50 km) and 0.22° × 0.22° (~ 25 km). We chose the latter option to take advantage of the spatial detail offered by the finer resolution RCM (ECCC, 2016). The CanRCM4 was developed by coordinating its parent Global Climate Model (GCM), such that all of the prognostic fields of the parent model were used to drive the CanRCM4. This, in turn, improved the quality of the model's output (Scinocca et al., 2015). Projection data from 2006 to 2100 are available for several fields. For our purposes, we downloaded and processed precipitation, maximum and minimum temperature, specific humidity (later converted to relative humidity), solar radiation, and wind speed data at 272 grids covering the ARB (Fig. 4). We chose two commonly used IPCC AR5 emission scenarios, the RCP4.5 and the RCP 8.5 (IPCC, 2014),

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

431

Fig. 4. Division of the Athabasca River Basin (ARB) into five regions. Also shown are the centroids of the CanRCM4 grids.

to define greenhouse gas concentration trajectories in the year 2100, compared to pre-industrial values of radiative forcing. Furthermore, we considered two time periods: a mid-century period with a time frame of 2021–2060, and a late-century period with a time frame of 2061–2100. In order to quantify the impact on the future climate in different spatial scales, we divided the ARB into five regions (Fig. 4). These regions were: headwaters (up to Athabasca River at Hilton), foothills (up to McLeod River near Whitecourt and the Athabasca River from Hilton to Windfall), Prairie (Pembina River at Jarvie), Lesser Slave (Lesser Slave Lake tributaries and the remaining part of the Athabasca River up to Athabasca), and boreal (all downstream sections). It should be noted that these sub-divisions are not official demarcations and that they have been used solely for our illustration purposes. 3. Results and discussion 3.1. Model validation for streamflow simulation in the base period Different parameters were found to be sensitive at different regions of the ARB. This was expected as the basin has varying topography, land-use, and soil types. In the headwater region, the lapse rates (PLAPS and TLAPS; Eqs. (3)–(8)) and snow related parameters (SMFMX and SMTMP; Eqs. (2)–(4)) were found to be the most sensitive. This sensitivity in the lapse rates follows expectations as the meteorological stations at the headwater region were sparse and variations in temperature and precipitation were generally significant with elevation. In the boreal region, however, the CANMX (maximum canopy storage) was found to be the most sensitive parameter. As this part of the basin had almost 90% of its area covered with forest (Table S1), this was also to be expected. Moreover, forest cover also played an

important role in evapotranspiration. Thus, it follows that the parameter ESCO (soil evaporation compensation factor), which controls the soil water available for evapotranspiration, was one of the sensitive parameters. Similarly, two groundwater parameters, GWQMN (threshold depth of water in the shallow aquifer required for return flow) and GW_DELAY (groundwater delay time), were also found to be sensitive in this region. In general, the snow related parameters (Table 1) were sensitive in all parts of the basin, as would be expected in a cold climate region like the ARB in which the hydrology is snow-dominated. The sensitivity ranking of the parameters were consistent with other studies of similar catchments (Ahl et al., 2008; Faramarzi et al., 2015; Grusson et al., 2015; Malagò et al., 2015; Troin and Caya, 2014). As we expected, optimized parameter values were found to vary as per HRUs and soil layers. For the purposes of clarity and due to the importance of snow-related parameters in a snow-dominated watershed within a cold climate region such as the ARB, we have presented the optimized values of snow-related parameters in this manuscript. The probability distributions of five snow-related parameters, calculated using the optimized values at all of the 131 sub-basins, are shown in Fig. S2. Comparisons of the optimized snow-related parameters to the optimized values reported by Ahl et al. (2008) and Troin and Caya (2014), indeed, show similarities (Table 3). One exception was seen for the SMTMP (Eqs. (2)–(4)), which was optimized at 5.9 °C, a slightly higher value than was observed for most of the other reported values. Streamflow simulation results during the base period (Fig. 5) using optimized parameters demonstrated the ability of the SWAT model to simulate the streamflow at various locations with varying degrees of accuracy (Fig. 6). The model results did not show any trend of systematic under- or over-estimation as indicated by the time series, scatter, and cumulative volume plots (Fig. 5). Model optimization resulted in an

432

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

Table 3 Average value of snow related parameters (calculated based on the optimized values of all 131 sub-basins of the ARB), compared with the optimized values reported in the literature. Parameters

SFTMP

SMTMP

SMFMX

SMFMN

TIMP

Default values Optimized values

1 −0.1 1.0 1.3 (1.5) 2.0 0.5(−0.5)

0.5 5.9 2.5 2.0 (−0.5) −0.1 2.0 (0.5)

4.5 1.9 3.0 5.0 (3.1) 5.8 4.5 (4.0)

4.5 2.6 2.9 3.2 (5.8) 1.2 6.0 (5.0)

1 0.19 0.06 0.14 (0.64) 0.05 1.0 (0.5)

This study Ahl et al. (2008) Grusson et al. (2015)a Troin and Caya (2014) Malagò et al. (2015)b

SFTMP: snowfall temperature (°C). SMTMP: snow melt base temperature (°C). SMFMX: maximum snow melt rate during the year (mm/°C·day). SMFMN: minimum snow melt rate during the year (mm/°C·day). TIMP: snow pack temperature lag factor (−). a Optimized values outside the parentheses are for snow parameter project and those inside the parentheses are for elevation band project. b Estimate mean values outside the parentheses are for Iberian Peninsula and those inside the parentheses are for Scandinavian Peninsula.

acceptable total predictive uncertainty band (p-statistics N 0.5 and rstatistics ~1) at the majority of the stations (Table S2) in both the calibration and the validation periods. In general, streamflow at the

headwater region was simulated with a “Very Good” accuracy rating. At the foothill region, the accuracy remained “Very Good” at the gauging stations located on the main river. This could have been due to the

Fig. 5. Time series plots (first column), scatter plots (second column), and cumulative plots (third column) of observed and simulated streamflow in the calibration (1990–2005) and validation (1983–1989 and 2006–2013) periods at (a) Headwater (Athabasca River at Hilton), (b) Foothill (McLeod River near Whitecourt), (c) Prairie (Pembina River at Jarvie), (d) Lesser Slave Lake (Athabasca River at Athabasca), and (e) Boreal (Athabasca River below Fort McMurray) regions. Also shown in the time series plots are the 95% total predictive uncertainty bands.

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

433

Fig. 6. Qualitative ratings of streamflow simulation results at 35 gauging stations across the Athabasca River Basin (ARB).

contribution of the “Very Good” quality of the simulated streamflow from the headwater region. However, the accuracy of streamflow results at tributaries varied from “Satisfactory” to “Good”, which was consistent with the findings of Eum et al. (2017). At the prairie region, the quality of streamflow results at upstream reaches were “Good” and were “Very Good” at downstream reaches. In the Lesser Slave region, while the quality of some tributaries' streamflow simulation results varied from “Satisfactory” to “Good”, a “Very Good” quality of streamflow was obtained from the majority of the headwater, foothill, and prairie regions. The model, however, underperformed at the lower parts of the boreal region, especially the North-Eastern part. The quality of streamflow simulation results in that region were “Satisfactory” at best. A streamflow simulation rating of “Very Good” was obtained

only at stations where a major contribution to the streamflow came from the upstream regions (e.g., the section of the Athabasca River below Fort McMurray). This underperformance can be attributed to the lack of a sufficient number of meteorological stations and to a large missing series of meteorological input data, such as precipitation and temperature (Fig. 2). Overall, the quality of streamflow simulation results were comparable to the reported values present in the literature (Betrie et al., 2015; Eum et al., 2017; Eum et al., 2014a; Eum et al., 2014b; Faramarzi et al., 2017; Leong and Donner, 2015) (Table 4). Availability of a calibrated and validated SWAT model during the base period formed the basis of our assessment of the impacts of future climate on different freshwater resources of the basin. This assessment will be presented in the following sections.

Table 4 Comparison of the Nash-Sutcliffe Efficiency (NSE) value obtained in this study, as an indicator of the accuracy of streamflow simulation results, with reported values from the literature. Station

Athabasca river at Hilton Athabasca river at Winfall Pembina river at Jarvie Athabasca river at Athabasca Athabasca river at Fort McMurray Muskeg river near Fort MacKay Firebag river at mouth All stations average

Region

Headwaters Foothills Prairie Lesser Slave Boreal Boreal Boreal –

This study

Literature values (Eum et al., 2017)1

(Eum et al., 2014a)2

(Eum et al., 2014b)3

(Betrie et al., 2015)4

(Faramarzi et al., 2017)5

(Leong and Donner, 2015)6

0.86 0.82 0.83 0.91 0.91 0.27 0.29 0.57

0.90 0.81 0.55 0.78 0.79 0.63 0.70 –

0.71 0.86 0.67 0.87 0.84 – 0.61 –

0.57 0.65 0.66 0.82 0.84 – 0.70 –

– – – – 0.79 – – –

– – – – – – – 0.21

– – – 0.50 0.35 – – –

1,2,3 Model used: process based Variable Infiltration Capacity (VIC). 4,5 Model used: process based Soil and Water Assessment Tool (SWAT). 6 Model used: land-surface model - Integrated BIosphere Simulator (IBIS) coupled with Terrestrial Hydrology Model with Biogeochemistry (THMB).

5.2 4.3 7.5 3.4 5.5 6.1 5.4 8.2 4.0 6.9 6.6 4.5 7.3

8.9

3.3. Climate change impacts PCP: precipitation in mm, change in %; TMP: average temperature in degree Celsius, change in degree Celsius.

5.5 4.9 7.4 3.0 6.8 1.0 1.2 3.2 −2.3 2.1

TMP

Relative to the base period, the ARB was projected to receive greater amounts of precipitation and experience warmer future conditions (Table 5). Especially in the high emission scenario (RCP 8.5), the latecentury period (2061–2100) was expected to be the wettest and warmest. However, there was evidence of marked temporal and spatial variability (Fig. 7; Table 5). On an annual time scale, relative to the base period, the mid-century (2021–2060) mean temperature was projected to increase by 2–3 °C, while the late-century period (2061–2100) was expected to increase by 2.5–5.4 °C. These findings were consistent with the observations made by Leong and Donner (2015) who reported an increase, of 0.9–3.1 °C for the mid-century period and of 0.5–7 °C for the late-century period, depending on which GCMs were used. Similarly, Golder-Associates (2009) reported an increase in annual mean temperature of 1.8 to 4.4 °C for the mid-century period, while Toth et al. (2006) reported an increase of 1.6 to 6.7 °C for the same period. Kerkhoven and Gan (2011) reported an increase of up to 9 °C for the warmest GCM output, while Eum et al. (2017) reported an increase of up to 3.3 and 5.6 °C for the mid- and late-century periods, respectively. The average annual precipitation was projected to increase by 21%–29% in the mid-century period while the late-century average annual precipitation was projected to increase by 24%–34%, relative to the base period. However, these increments were slightly higher than the values that have been reported by other studies. For example, Leong and Donner (2015) reported an increase of 12% and 6% for the respective mid-century and late-century periods, while Kerkhoven and Gan (2011) found an increase of 20% for their wettest GCM. Using the ensemble mean of several GCMs, Eum et al. (2017) reported an increase of 7%–10% and 13%–14% for the mid-century and late-century periods, respectively. In another study, Golder-Associates (2009) reported an increase of up to 12% in mean annual precipitation during the midcentury.

5.8

2.2 0.9 4.2 1.1 2.5 3.2 1.8 4.5 1.9 4.8 3.7 −0.8 −4.6 0.2

−2.5

−1.9

5.1

0.9

3.6

1.4

2.8

5.0 2.5 5.4

2.1

2.8 1.1 5.1 2.0 3.1 3.8 2.0 5.5 2.5 5.0 4.2 0.2 −4.1 0.6

−2.3

−1.4

5.4

1.4

4.5

1.6

3.2

5.9 3.1 5.4

2.5

1.8 0.8 3.9 1.0 1.3 2.8 1.8 4.2 1.6 3.5 3.2 4.6 2.2 −1.3 −4.9 −0.8

−2.7

−2.4

3.9

0.6

3.2

1.3

2.3

4.1

2.1

2.5 3.3 15.3 −11.3 2.8 2.3 2.8 14.2 −10.3 2.5 2.7 3.3 14.9 −10.4 3.0 1.9 2.5 12.7 1.9 −9.5 0.4 0.7 11.2 −10.6 0.4

41% 70% 21% 89% 7% 20% 36% 0% 76% −5% 25% 85% N100% 44% 75%

70%

−16% 64%

8%

32%

21%

5% 82% −2%

34%

28% 41% 18% 59% 5% 13% 5% 11% 39% −14% −3% 19% 47% 15% 69% 45% 76% 78%

62%

−32% 27%

6%

11%

−26%

27%

33% 65% 21% 56% 3% 17% 27% 9% 43% −9% 26% 20% 47% 26% 53%

85% 71% 77%

68%

−23% 24%

26%

20%

−14%

30%

26% 44% 21% 44% −3% 9% 5% 9% 28% −17% −2% 20% 39% 16% 42%

69% 69% 83%

60%

−27% 25%

−2% 9%

−22%

30%

235.6 79.6 467.4 91.0

Base period (1983–2013) Future: RCP 4.5 2040's (2021–2060) Future: RCP 4.5 2080's (2061–2100) Future: RCP 8.5 2040's (2021–2060) Future: RCP 8.5 2080's (2061–2100) Base period (1983–2013) Future: RCP 4.5 2040's (2021–2060) Future: RCP 4.5 2080's (2061–2100) Future: RCP 8.5 2040's (2021–2060) Future: RCP 8.5 2080's (2061–2100) PCP

Boreal

61.2 272.1 87.9 518.9 99.4 59.5 482.9 97.5 57.2 156.4 115.1

314.5 164.0 750.0

66.5

115.3 290.8 98.8

571.4

239.7 88.5

MAM JJA Annual DJF MAM JJA

Prairie

SON MAM JJA

Foothills

Annual DJF SON JJA MAM DJF

Headwaters Variables Periods

3.2. Climate change projection

Annual DJF

SON

Lesser Slave

SON Annual DJF

MAM JJA

SON Annual

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

Table 5 Regional averaged seasonal and annual changes in precipitation and temperature in mid- (2021–2060) and late-century (2061–2100) periods for two emission scenarios (RCP 4.5 and RCP 8.5), relative to the base period (1983–2013).

434

3.3.1. Blue water yield and blue water flow (streamflow) Due to a projected increase in precipitation, both the surface and sub-surface flow are expected to increase significantly on an annual basis (Table 6). It is clear that the increase in total surface and subsurface flow was more significant than the increase in precipitation, as was indicated by a non-linear relationship between the two factors. For every 1% rise in annual precipitation, an increase of 1.4% (r2 = 0.97, n = 4) in total flow was estimated. Similarly, evidence such as the increase in the ratio of streamflow to precipitation and the decreased contribution of sub-surface flow to the total flow, implied that this projected precipitation increase would be primarily due to intense rainfall events which would generate greater amounts of surface runoff. Hence, more frequent incidences of flash-floods could be expected in the future, as suggested by results of the streamflow frequency analysis plots (Fig. S3). On a finer temporal (seasonal) and spatial (sub-basin and regionaggregated) scale, the variability of future changes in precipitation and temperature (Fig. 7) has been reflected in both blue water yield (Fig. 8) and blue water flow (Fig. 9). Despite the expected decrease in future winter precipitation everywhere except in the headwater region, the blue water yield was expected to increase (Table 7) due to increased snowmelt driven by the projected increase of winter temperatures. This would decrease the snowpack volume for the spring freshet. Hence, we would expect the blue water yield during the spring season to decrease, despite a projected increase in temperature and precipitation during this time. In contrast, we would expect the lower regions (Lesser Slave and Boreal regions) of the basin to experience a substantial increase (N 100%) in spring season blue water yield. We speculate this would be due to an increased sub-surface contribution to the total blue water yield in these regions. Runoff from the snowmelt would

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

435

Fig. 7. Seasonal changes in precipitation and (mean) temperature for two emission scenarios (RCP 4.5 and RCP 8.5) during mid- (2021–2060) and late- (2061–2100) century periods, relative to the base period (1983–2013).

have been expected to infiltrate because the topography of the basin in the lower regions was flat and forest was the main land-use type. Thus, a significant delay, controlled by the parameter GW_DELAY

(groundwater delay), was expected for the infiltrated water to reappear in the stream. This particular parameter was one of the most sensitive parameters in this region of the basin. During the summer

Table 6 Basin-wide averaged annual water balance components in the base (1983–2013) and future periods (mid-century: 2021–2060 and late-century: 2061–2100) for two emission scenarios (RCP 4.5 and 8.5). Water balance components

Precipitation (mm) Surface runoff (mm) Sub-surface flow (mm) Total flow/precipitation Sub-surface/total flow ET (mm) ET/precipitation PET (mm)

Base

RCP 4.5

% Change (RCP 4.5)

% Change (RCP 8.5)

(1983–2013)

2040's (2021–2060)

2080's (2061–2100)

2040's (2021–2060)

RCP 8.5 2080's (2061–2100)

2040's (2021–2060)

2080's (2061–2100)

2040's (2021–2060)

2080's (2061–2100)

510 40 60 20% 60% 378 74% 538

636 61 84 23% 58% 413 65% 559

676 67 88 23% 57% 437 65% 590

648 63 87 23% 58% 417 64% 564

697 73 90 23% 55% 462 66% 658

25% 52% 40% – – 9% – 4%

33% 67% 46% – – 16% – 10%

27% 58% 44% – – 10% – 5%

37% 83% 49% – – 22% – 22%

436

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

Fig. 8. Seasonal changes in the blue water yield (WYLD), green water storage (soil water - SW), and green water flow (evapotranspiration – ET) in two emission scenarios (RCP 4.5 and RCP 8.5) during mid- (2021–2060) and late- (2061–2100) century periods, relative to the base period (1983–2013).

season, both precipitation levels and temperature were expected to increase, potentially causing an offsetting effect. In the headwater region, which had the highest proportion of barren land (Table S1), the effect of

increased precipitation seemed to outweigh the effect of the increased temperature, demonstrated by an increased summer blue water yield of up to 63%. Because the foothill region experienced up to 95% forest

Fig. 9. Comparison of monthly averaged streamflow (first and third columns) and percentage change in streamflow (second and fourth columns) during the base period (1983–2013), and during the mid- (2021–2060) and late- (2061–2100) century periods for two emission scenarios (RCP 4.5 and 8.5) at (a) Headwater (Athabasca River at Hilton), (b) Foothill (McLeod River near Whitecourt), (c) Prairie (Pembina River at Jarvie), (d) Lesser Slave Lake (Athabasca River at Athabasca), and (e) Boreal (Athabasca River below Fort McMurray) regions. Error bars in the first and third columns show the maximum and minimum range of monthly averaged streamflow.

Table 7 Seasonal changes in blue water yield and blue water flow (streamflow) during the mid- (2021–2060) and late-century (2061–2100) periods for two emission scenarios (RCP 4.5 and RCP 8.5), relative to the base period (1983–2013). Also presented is the comparison of annual changes in streamflow with values found in the literature. Variables

Periods

Headwaters DJF

WYLDa

Eum et al. (2017) Leong and Donner (2015) Kerkhoven and Gan (2011) Toth et al. (2006)

Foothills JJA

SON Annual DJF

Base period 39.7 112.4 318.6 88.9 559.6 (1983–2013) RCP 4.5 40's 15% −47% 45% 21% 20% (2021–2060) RCP 4.5 80's 26% −37% 59% 33% 33% (2061–2100) RCP 8.5 40's 18% −43% 50% 20% 24% (2021–2060) RCP 8.5 80's 34% −2% 63% 70% 49% (2061–2100) Base period 27.2 111.3 423.9 95.2 164.4 (1983–2013) RCP 4.5 40's −4% −61% 44% 25% 21% (2021–2060) RCP 4.5 80's 3% −55% 63% 38% 37% (2061–2100) RCP 8.5 40's −2% −58% 50% 23% 26% (2021–2060) RCP 8.5 80's 9% −17% 76% 99% 61% (2061–2100) Mid-century – – – – 7% End-century – – – – 12% −10% to +53% depending on the GCMs used

Prairie JJA

MAM

JJA

65.9

106.9 39.5

231.8

3.2

32.2

30.9 9.5

−1% 0%

39%

−14% 77%

24%

11%

85%

−1%

−15% −12% 14%

1%

2%

36%

5%

20%

1%

16%

9%

N100% 9%

17%

14.5

68.0

92.1

34.2

52.2

6.4

46.8

18%

−10% 23%

20%

11%

48%

57%

2%

30%

48%

26%

18%

−7%

27%

19%

45%

15%

6%

– –

– –

4%

−4%

MAM

JJA

SON

Annual DJF

MAM

JJA

SON

Annual

75.8

10.5

33.9

73.4

27.3

145.1

10.1

33.0

47.7

19.2

109.9

53%

33%

9%

19%

17%

9%

15%

N100% 83%

96%

129%

99%

56%

77%

35%

23%

52%

16%

25%

27%

N100% N100% 83%

145%

107%

−11% 63%

40%

28%

17%

38%

21%

14%

23%

N100% N100% 94%

119%

103%

66%

24%

38%

67%

4%

20%

24%

N100% N100% 87%

134%

115%

48.5 15.7

117.5

137.2 377.8

48%

48%

48%

44%

0%

−18% 37%

19%

97%

9%

55%

N100% 45%

18%

0%

50%

14%

43%

−3%

67%

56%

36%

4%

−8%

27%

15%

114%

17%

14%

58%

27%

18%

– –

– –

6% 11%

– –

– –

– –

– –

8% 23%

– – – –

– – – –

19% – – –

– – – –

– – – –

– – – –

– – – –

4% – – –

17%

SON

Boreal

Annual DJF

−20% −13% 11%

SON

Lesser Slave

Annual DJF

19.4

MAM

878.4 300.7 423.5

218.4

528.4

1124.2 432.3 575.9

19%

34%

19%

56%

48%

44%

37%

34%

51%

43%

64%

64%

58%

43%

19%

24%

39%

35%

64%

51%

52%

22%

45%

47%

38%

52%

66%

58%

67%

61%

– –

– –

– –

– –

5% 13%

– –

– –

– –

– –

4% 14%

– – – –

– – – –

– – – –

– – – –

9% – −12% –

– – – –

– – – –

– – – –

– – – –

7% – −8% –

−5% to +37% depending on the GCMs used and SRES adopted

Mid-century End-century Golder-Associates Mid-century (2009) End-century

– – – –

– – – –

– – – –

– – – –

28% – – –

– – – –

– – – –

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

Streamflowb

MAM

1

Eum et al. (2017): model used: process based Variable Infiltration Capacity (VIC), 6 CMIP5 GCM projection for RCP 4.5 and 8.5 scenarios. Leong and Donner (2015): model used: land-surface model - Integrated BIosphere Simulator (IBIS) coupled with Terrestrial Hydrology Model with Biogeochemistry (THMB), 3 CMIP5 GCM projection for RCP 2.5 to 8.5 scenarios. At a stations downstream of Athabasca river below Fort McMurray. 3 Kerkhoven and Gan (2011): modified Soil-Biosphere-Atmosphere (MIBSA), seven GCMs were used and four SRES (A2, B1, B2, A1F1). 4 Toth et al. (2006): distributed hydrological model - WATFLOOD, seven GCMs used. 5 Golder-Associates (2009): process based Hydrologic Simulation Program - Fortran (HSPF) model, five GCMs and five SRES (A1F1, A1, A2, B1, B2). a WYLD = water yield (mm). b Streamflow in cumecs; Streamflows measured at: Headwater - Athabasca river at Hilton, Foothill - McLeod river near Whitecourt, Prairie - Pembina river near Jarvie, Lesser Slave - Athabasca river at Athabasca, Boreal - Athabasca river below Fort McMurray. 2

437

438

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

coverage (Table S1), an 8% increase in summer precipitation during the end-century period in the high emission scenario (Table 5) led to only a 1% increase in blue water yield. This suggests that incremental temperature changes might have affected blue water yield more in this region during the summer season as a result of increased green water flow (evapotranspiration) from the mostly forested landscape (Table S1) of the foothill region. During the fall season, changes in precipitation and temperature have been duly reflected in the blue water yield. The exception to this was observed in the middle and lower regions of the basin (e.g., prairie regions) where even projected decreases in autumn precipitation resulted in increased blue water yields. The effect of precipitation and temperature changes on the blue water yield could be synergetic or offsetting, depending on the seasons and regions. During the spring season, increases in precipitation and temperature seemed to have a synergetic effect on the blue water yield across all regions, as indicated by positive slopes of the trendlines (Fig. S4) drawn between these two variables and the blue water yield. During the summer season, however, contrasting slopes of trend-lines indicated that an increase in precipitation and temperature had an offsetting effect in all regions except for the headwater and boreal regions. In the headwater region, for the same increase in autumnal precipitation (85%) in the late century period for the two emission scenarios, a difference in temperature changes (−2.3 °C in the lower emission scenario and + 1.2 °C in the higher emission scenario) (Table 5) resulted significant differences in the hydrological responses. The difference in blue water yield (+33% versus +70% in the lower and higher emission scenarios, respectively) (Table 7) and blue water flow (+ 38% versus + 99%) (Table 7) indicated that this region was more sensitive to the projected temperature changes than it was to the projected precipitation changes. A similar observation was made by Kerkhoven and Gan (2011). This can be attributed to several factors such as the change in precipitation type corresponding to changes in the temperature. The average autumnal temperature for the base period was 0.7 °C (Table 5). A projected decrease in temperature by 2.3 °C in the lower emission scenario would mean that most of the precipitation would occur in the form of snow. The changes in the blue water flow (streamflow) dynamics have, as expected, followed the trend of the blue water yield at the representative stations of each region. In general, the winter streamflow tended to increase throughout the basin (Table 7). Due to earlier spring freshet, the spring streamflow tended to decrease at the upper reaches while it increased at the lower reaches. A significant increase in the summer and autumn streamflow was expected for both future periods and scenarios in all reaches of the basin. The average annual streamflow was expected to increase by 36%, 16%, 38%, 29%, and 54% at the outlet of the headwater, foothill, prairie, Lesser Slave, and at a representative station of boreal region, respectively. Hence, it is evident that the boreal region will likely be the main source of blue water flow in the future, a result that was congruent with the findings of Rood et al. (2015). Comparing the changes in streamflow from our study with those found in the literature revealed contrasting scenarios (Table 7). It should, however, be noted that the results of our study are not directly comparable to other studies due to key differences in the model we used, future climatic forcing, and the baseline period we chose. Our finding of a basin-wide average increase of 35% in annual streamflow with a 54% increase occurring at the farthest downstream location, seemed to be consistent with the findings of Leong and Donner (2015) with their “wettest” GCM forcing, for which the study reported an increase of 53%. Similarly, Toth et al. (2006) obtained results at headwater and foothill regions which were also comparable to our findings. Nonetheless, our results were consistently higher than those present in other studies (Eum et al., 2017; Golder-Associates, 2009; Kerkhoven and Gan, 2011). These researchers reported mean ensemble results which tended to have an averaging effect, giving projections of wetter and drier, and warmer and colder future climates. Furthermore, there exists a great deal of uncertainty related to climate change projections (Latif, 2011). Moreover, different

sources of uncertainties such as input data and model parameters (Leta et al., 2015) are inherent in hydrological models. Although using a multi-model ensemble would reduce the overall uncertainty of climate change projections (Scinocca et al., 2015), we believe that our decision to use the CanRCM4 data set, the result of coordinated global and regional climate modelling which “offers a number of practical and conceptual simplifications that are generally not available to independent RCM modelling centers” (Scinocca et al., 2015), was a reliable one in the case of the ARB. Indeed, a recent study by Diaconescu et al. (2016) showed that the CanRCM4 dataset was one of the best of all North American models currently available for use. 3.3.2. Green water storage and green water flow As has previously been stated, an increased annual blue water yield to precipitation ratio and lower contributions of the sub-surface flow to the total future flow (Table 6) indicated a lower residence time of precipitated water in the basin. This has the potential to affect the green water storage and green water flow dynamics of the ARB. It was found that the proportion of annual green water flow (evapotranspiration) to precipitation was projected to decrease from an average of 74% to 65% in the future, despite a projected increase in the annual green water flow of up to 22%. Hence, a major portion of an annual precipitation increase in the future would also be observed in the increase of the yield of blue water. While analyzing the impacts of climate change on green water storage on a regional scale across seasons (Fig. 8; Table 8), it was clear that the changes observed in green water storage in the upper regions (headwater and foothill) of the basin were less sensitive to the projected precipitation changes than were those in the downstream regions. This can be attributed to the region's soil profile. At the headwater region, the soil profile was shallow and 55.8% of the region's area was exposed bed rock. Hence, a lower amount of water could infiltrate and be stored in the soil stratum. Furthermore, stored soil water would be depleted quickly due to a projected increase in green water flow. In the lower Prairie and Lesser Slave Lake regions, agricultural activities were significant due to the fact that 45% and 20%, respectively, of the total areas of the regions were made up of agricultural land. The associated tillage activities would, thus, allow for a significant portion of precipitated water to infiltrate these regions. Similarly, the boreal plain had a flat topography, 95% covered by forest, allowing for high levels of infiltration. Furthermore, increased amounts of plant biomass, with the exception being the headwater region, (Table S3) could have decreased the runoff coefficient (Kerkhoven and Gan, 2011). Therefore, both the projected increase in precipitation and in temperature might have favoured an increased amount of soil water storage in these downstream regions. Additionally, a deep soil stratum would also have been expected to dampen the short-term variations of climatic forcing. With the availability of greater green water storage and a projected increase in temperature, the green water flow increased substantially in all regions across seasons. While analyzing the water resource components on a monthly time scale (Table S4) in the prairie region, it was clear that that the green water flow would decrease by an average of 5% in May during the mid-century period, despite a projected increase in precipitation and in temperature for that month, due to the fact that a high proportion of the region was made up of agricultural land. During the month of May, evaporation from top soil on the surface contributed mote to the green water flow as most of the sown crops, such as Spring Wheat with sowing dates falling between May 1 and May 30 (AGRI-FACTS, 2013), would have just germinated. Due to earlier spring freshet, infiltrated water would have percolated to deeper soil stratums, leaving the top soil layer with lower amounts of green water storage. As a result, the evaporation from the top layer would decrease, along with the overall amount of green water flow. Similarly, a serious depletion of green water storage (23%) was observed for the month of August. During this period, sown crops would be in in their mature stage and the leaf area index (LAI) would be at its highest. Soil evaporation would,

20% 19% 15% N100% 28% 20% 11% N100% 48% 23% 24% 14% N100% 34% 29% SW: Green water storage in mm, change in %, ET: Green water flow in mm, change in %

25% 54%

63%

53% 31% 60%

N100% 36%

16%

18% 20% 14% N100% 13% 18% 19% 26%

32%

25% 7% 62%

N100% 9%

8%

22% 30% 18% N100% 21% 27% 18% 28% N100% 18% 36%

44%

35% 12% 51%

ET

thus, decrease as less sunlight would be able to reach the soil surface. Furthermore, crops would respond to the increased temperatures of up to 9.9 °C for August during the late century in a higher emission scenario by closing their stomata (Schulze et al., 1973). This would be especially true in water stressed conditions, as would be typically observed during the summer season in this region (Table S3). All of these processes would have had the potential to contribute to the observed 8% decrease in green water flow during the month of August. The projected increase in temperature seemed to have positive effect on plant biomass, as the temperature stress was lowered (Table S3) across all seasons for all regions. Higher green water flow, due to increased plant biomass and a projected increase in temperature resulted in a decrease in green water storage, especially during the summer months of the late-century period for a higher emission scenario. Thus, an increase in water stress (Table S3) was evident, especially in the prairie and Lesser Slave regions where agriculture activity levels were the highest. Such an increase in water stress during the growing and flowering stages of crops would result in decreased crop yield (Shrestha and Shrestha, 2017). Hence, either seedling dates need to be shifted or additional water via irrigation must be supplied in the summer months to stabilize crop yields in the future.

14%

7% 6% 7% N100% −1% 10% 7% N100% 7%

6%

12% 13% 10% N100% 10% 14% 10% N100% 18%

13%

6% 5% 7% N100% −2% 8% 5% 7% N100% 3% 17% 20% 15% N100% 11% 18% 8% 20% 24%

25%

22% 8% 41%

N100% 4%

386.0 255.9 57.0 70.4 113.1 39.9

175.6

2.6 20.7 2.0

67.5

220.9

61.0

352.0

1.6

92.8

254.1 53.7

402.3

264.9 58.8

391.0 65.3 2.0

2.7

26% 27% 15% 29% 1% 1%

−2%

−2%

−3% −4% −5%

−3%

−17% −9% −7%

4%

6%

−9%

−7% −1%

−7%

6% 15% 12%

34%

31% 36% 25% 29% 19% 5%

4%

1% −3% −1%

−10%

−6%

1%

1%

−4%

5%

4%

10%

13%

8%

11%

16% 16% 16%

35%

32% 40% 21% 30% 5%

3%

0%

3% −4% −2%

−1%

2%

11%

3%

16%

12%

10%

21%

14%

7%

22%

17% 20% 21%

38%

30% 37% 25% 28% 33% 11%

19%

14% 15% 15% 9% 15% 14% 5% 5% −3% 3% 4% −6% −12% 6%

4%

1% −4% −2%

Annual SON MAM JJA

134.8 124.2 113.5 123.1 119.9 126.9 119.1 105.9 116.2 112.7 110.6 129.5 106.6 98.4 107.9 67.8 57.3 62.5 69.6 135.7 141.5 134.3 128.1 134.9

81.6

439

4. Conclusions

Base Period (1983-2013) Future: RCP 4.5 2040's (2021–2060) Future: RCP 4.5 2080's (2061–2100) Future: RCP 8.5 2040's (2021–2060) Future: RCP 8.5 2080's (2061–2100) Base Period (1983–2013) Future: RCP 4.5 2040's (2021–2060) Future: RCP 4.5 2080's (2061–2100) Future: RCP 8.5 2040's (2021–2060) Future: RCP 8.5 2080's (2061–2100) SW

Boreal

Annual DJF SON MAM JJA

Lesser Slave

Annual DJF SON MAM JJA

Prairie

Annual DJF SON MAM JJA

Foothills

Annual DJF SON MAM JJA DJF

Headwaters Variables Periods

Table 8 Regional averaged seasonal and annual changes in green water storage (SW) and green water flow (ET) during the mid- (2021–2060) and late-century (2061–2100) periods for two emission scenarios (RCP 4.5 and RCP 8.5), relative to the base period (1983–2013).

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

In this study, we built, calibrated, and validated a Soil and Water Assessment Tool (SWAT) model of the Athabasca River Basin (ARB) in order to quantify the impacts of climate change on the blue and green water resources of the basin at different temporal and spatial scales. The model was able to reproduce the historical streamflow dynamics at 35 gauging stations across the basin with varying degrees of accuracy. Future climate data projected wetter and warmer climatic conditions for the basin. These factors seem to have both synergetic and offsetting effects on the water resources of basin, depending on the region and season. On average, the water resources in the ARB are likely to increase in the future with the highest contribution expected to come from the lower regions of the basin, especially the boreal region. Significant increases (16–54%) in annual streamflow held the potential to pose flooding problems across the basin. Increment in plant biomass was observed across the basin as a result of decreased temperature stress occurring during all seasons. The biomass increment, as well as a warmer and wetter future climatic condition, led to higher green water flow (9–22%) from the basin. Consequently, the green water storage was projected to decrease, especially during the summer and autumn seasons for the late-century period in the middle regions of the basin where agricultural activity levels were significant. Plants would be expected to experience increased water stress which might require a solution in the form of the artificial supply of water. There was ample evidence of temporal and spatial heterogeneity of the blue and green resources of the basin for the future. This could pose many challenges to water resource planners and managers working to ensure smooth future development and sustainability of northern Alberta communities, as well as to the future environmental and economic sustainability of the river basin. Acknowledgements The authors would like to thank the Alberta Economic Development and Trade for the Campus Alberta Innovates Program Research Chair (No. RCP-12-001-BCAIP). We would also like to thank Mr. Jim Sellers for the proofreading. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.scitotenv.2017.05.013.

440

N.K. Shrestha et al. / Science of the Total Environment 601–602 (2017) 425–440

References Abbaspour, K.C., Johnson, C.A., van Genuchten, M.T., 2004. Estimating uncertain flow and transport parameters using a sequential uncertainty fitting procedure. Vadose Zone J. 3 (4), 1340–1352. Abbaspour, K.C., et al., 2007. Modelling hydrology and water quality in the pre-alpine/ alpine Thur watershed using SWAT. J. Hydrol. 333 (2–4):413–430. http://dx.doi. org/10.1016/j.jhydrol.2006.09.014. Abbaspour, K.C., et al., 2015. A continental-scale hydrology and water quality model for Europe: calibration and uncertainty of a high-resolution large-scale SWAT model. J. Hydrol. 524:733–752. http://dx.doi.org/10.1016/j.jhydrol.2015.03.027. AGRI-FACTS, 2013. Dryland agronomic management of soft white spring wheat and Canada prairie spring wheat in Alberta. http://www1.agric.gov.ab.ca/$Department/ deptdocs.nsf/all/agdex14523 (In: Available, o. (Eds.)). Ahl, R.S., Woods, S.W., Zuuring, H.R., 2008. Hydrologic calibration and validation of SWAT in a snow-dominated Rocky Mountain Watershed, Montana, U.S.A.1. J. Am. Water Resour. Assoc. 44 (6):1411–1430. http://dx.doi.org/10.1111/j.1752-1688.2008. 00233.x. Ahmadi, M., Ascough Ii, J.C., DeJonge, K.C., Arabi, M., 2014. Multisite-multivariable sensitivity analysis of distributed watershed models: enhancing the perceptions from computationally frugal methods. Ecol. Model. 279:54–67. http://dx.doi.org/10.1016/ j.ecolmodel.2014.02.013. Arnold, J.G., Srinivasan, R., Muttiah, R.S., Williams, J.R., 1998. Large area hydrologic modeling and assessment part I: model development. J. Am. Water Resour. Assoc. 34 (1), 73–89. Arnold, J.G., et al., 2011. Soil and Water Assessment Tool Input/output File Documentation, Version 2009. 2011. Agrilife Blackland Research Center, Temple, Texas, p. 76502. Arnold, J.G., et al., 2012. SWAT: model use, calibration, and validation. Am. Soc. Agric. Biol. Eng. 55 (4), 18. AWC, 2011. Athabasca Watershed Council - State of the Watershed Report: Phase 1. Athabasca Watershed Council, Hilton, Alberta. AWC, 2013. State of the Watershed Report: Phase 3 - Water Quantity and Basic Water Quality in the Athabasca Watershed. Athabasca Watershed Council. Betrie, G.D., Deng, B., Wang, J., 2015. Integrated Modelling of the Athabasca River Basin Using SWAT, 2015 Proceedings of Science and Technology Innovations. Athabasca University Press, Alberta, Canada, pp. 27–38. Cao, W., Bowden, W.B., Davie, T., Fenemor, A., 2006. Multi-variable and multi-site calibration and validation of SWAT in a large mountainous catchment with high spatial variability. Hydrol. Process. 20 (5):1057–1073. http://dx.doi.org/10.1002/hyp.5933. CFSR, 2016. In: G.W.D.f.S. (Ed.), Climate Forecast System Reanalysis http://globalweather. tamu.edu/. Diaconescu, E.P., Gachon, P., Laprise, R., Scinocca, J.F., 2016. Evaluation of precipitation indices over North America from various configurations of regional climate models. Atmosphere-Ocean 54 (4):418–439. http://dx.doi.org/10.1080/07055900.2016. 1185005. Dile, Y.T., Srinivasan, R., 2014. Evaluation of CFSR climate data for hydrologic prediction in data-scarce watersheds: an application in the Blue Nile River Basin. J. Am. Water Resour. Assoc. 50 (5):1226–1241. http://dx.doi.org/10.1111/jawr.12182. ECCC, 2016. Climate Data and Scenarios for Canada: Synthesis of Recent Observation and Modelling Results. Environment and Climate Change Canada, Quebec, Canada. Environment-Canada, 2016. In: A.a. (Ed.), Historical Hydrometric Data http://www.ec.gc. ca/rhc-wsc/. Eum, H.-I., Dibike, Y., Prowse, T., Bonsal, B., 2014a. Inter-comparison of high-resolution gridded climate data sets and their implication on hydrological model simulation over the Athabasca Watershed, Canada. Hydrol. Process. 28 (14):4250–4271. http:// dx.doi.org/10.1002/hyp.10236. Eum, H.-I., Yonas, D., Prowse, T., 2014b. Uncertainty in modelling the hydrologic responses of a large watershed: a case study of the Athabasca River basin, Canada. Hydrol. Process. 28 (14):4272–4293. http://dx.doi.org/10.1002/hyp.10230. Eum, H.-I., Dibike, Y., Prowse, T., 2017. Climate-induced alteration of hydrologic indicators in the Athabasca River Basin, Alberta, Canada. J. Hydrol. 544:327–342. http://dx.doi. org/10.1016/j.jhydrol.2016.11.034. Falkenmark, M., Rockström, J., 2006. The New Blue and Green Water Paradigm: Breaking New Ground for Water Resources Planning and Management (DOI:10.1061/ (ASCE)0733-9496(2006)132:3(129)). Faramarzi, M., et al., 2015. Setting up a hydrological model of Alberta: data discrimination analyses prior to calibration. Environ. Model. Softw. 74:48–65. http://dx.doi.org/10. 1016/j.envsoft.2015.09.006. Faramarzi, M., et al., 2017. Uncertainty based assessment of dynamic freshwater scarcity in semi-arid watersheds of Alberta, Canada. J. Hydrol. 9:48–68. http://dx.doi.org/10. 1016/j.ejrh.2016.11.003. Fuka, D.R., et al., 2014. Using the Climate Forecast System Reanalysis as weather input data for watershed models. Hydrol. Process. 28 (22):5613–5623. http://dx.doi.org/ 10.1002/hyp.10073. GoC, 2016. In: G.o.C. (Ed.), Environment and Natural Resources; Weather, Climate and Hazard http://climate.weather.gc.ca/. Golder-Associates, 2009. Hydro-Climate Model Selection and Application on the Athabasca and Beaver River Basins. Oil Sands Environmental Management Division, Alberta Environment, Calgary, Alberta, Canada. Grusson, Y., et al., 2015. Assessing the capability of the SWAT model to simulate snow, snow melt and streamflow dynamics over an alpine watershed. J. Hydrol. 531 (Part 3):574–588. http://dx.doi.org/10.1016/j.jhydrol.2015.10.070.

IPCC, 2014. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, IPCC. IPCC. Jarvis, A., Reuter, H.I., Nelson, A., Guevara, E., 2008. Hole-filled SRTM for the Globe Version 4, Available From the CGIAR-CSI SRTM 90m Database. CGIAR http://srtm.csi.cgiar.org. Kerkhoven, E., Gan, T.Y., 2011. Differences and sensitivities in potential hydrologic impact of climate change to regional-scale Athabasca and Fraser River basins of the leeward and windward sides of the Canadian Rocky Mountains respectively. Clim. Chang. 106 (4):583–607. http://dx.doi.org/10.1007/s10584-010-9958-7. Latif, M., 2011. Uncertainty in climate change projections. J. Geochem. Explor. 110 (1): 1–7. http://dx.doi.org/10.1016/j.gexplo.2010.09.011. Leong, D.N.S., Donner, S.D., 2015. Climate change impacts on streamflow availability for the Athabasca Oil Sands. Clim. Chang. 133 (4):651–663. http://dx.doi.org/10.1007/ s10584-015-1479-y. Leta, O.T., et al., 2015. Assessment of the different sources of uncertainty in a SWAT model of the River Senne (Belgium). Environ. Model. Softw. 68:129–146. http://dx.doi.org/ 10.1016/j.envsoft.2015.02.010. Liersch, S., 2003. The Program pcpSTAT. http://swat.tamu.edu/media/83150/manual_ pcpSTAT.pdf (Berlin). Loveland, T.R., et al., 2000. Development of a global land cover characteristics database and IGBP DISCover from 1 km AVHRR data. Int. J. Remote Sens. 21 (6–7): 1303–1330. http://dx.doi.org/10.1080/014311600210191. Malagò, A., Pagliero, L., Bouraoui, F., Franchini, M., 2015. Comparing calibrated parameter sets of the SWAT model for the Scandinavian and Iberian peninsulas. Hydrol. Sci. J. 60 (5):949–967. http://dx.doi.org/10.1080/02626667.2014.978332. Moriasi, D.N., et al., 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50 (3), 885–900. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I — a discussion of principles. J. Hydrol. 10 (3):282–290. http://dx.doi.org/10.1016/00221694(70)90255-6. Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., 2011. Soil & Water Assessment Tool Theoretical Documentation, Version 2009. Grassland, Soil and Water Research Laboratory-Agricultural Research Service, Blackland Research Center-Texas AgriLife Research. Pradhanang, S.M., et al., 2011. Application of SWAT model to assess snowpack development and streamflow in the Cannonsville watershed, New York, USA. Hydrol. Process. 25 (21):3268–3277. http://dx.doi.org/10.1002/hyp.8171. Rood, S.B., Stupple, G.W., Gill, K.M., 2015. Century-long records reveal slight, ecoregionlocalized changes in Athabasca River flows. Hydrol. Process. 29 (5):805–816. http:// dx.doi.org/10.1002/hyp.10194. Santhi, C., et al., 2001. Validation of the SWAT model on a large river basin with point and nonpoint sources. J. Am. Water Resour. Assoc. 37 (5), 1169–1188. Sauchyn, D.J., St-Jacques, J.-M., Luckman, B.H., 2015. Long-term reliability of the Athabasca River (Alberta, Canada) as the water source for oil sands mining. Proc. Natl. Acad. Sci. 112 (41):12621–12626. http://dx.doi.org/10.1073/pnas.1509726112. Schulze, E.D., Lange, O.L., Kappen, L., Buschbom, U., Evenari, M., 1973. Stomatal responses to changes in temperature at increasing water stress. Planta 110 (1):29–42. http://dx. doi.org/10.1007/BF00386920. Scinocca, J.F., et al., 2015. Coordinated global and regional climate modeling. J. Clim. 29 (1):17–35. http://dx.doi.org/10.1175/JCLI-D-15-0161.1. Shrestha, L., Shrestha, N.K., 2017. Assessment of climate change impact on crop yield and irrigation water requirement of two major cereal crops (rice and wheat) in Bhaktapur district. Nepal. J. Water Clim. Chang. 8 (2):320–335. http://dx.doi.org/10. 2166/wcc.2016.153. Shrestha, N.K., Leta, O.T., De Fraine, B., van Griensven, A., Bauwens, W., 2013. OpenMIbased integrated sediment transport modelling of the river Zenne, Belgium. Environ. Model Softw. 47:193–206. http://dx.doi.org/10.1016/j.envsoft.2013.05.004. Shrestha, N.K., et al., 2014. Modelling Escherichia coli dynamics in the river Zenne (Belgium) using an OpenMI based integrated model. J. Hydroinf. 16 (2), 354–374. Shrestha, N.K., Leta, O.T., Bauwens, W., 2017. Development of RWQM1-based integrated water quality model in OpenMI with application to the River Zenne. Belgium. Hydrol. Sci. J. 62 (5):774–799. http://dx.doi.org/10.1080/02626667.2016.1261143. SLC, 2010. In: A.a.A.-F.C. (Ed.), Soil Landscapes of Canada Version 3.2 http://sis.agr.gc.ca/ cansis/nsdb/slc/v3.2/index.html. Strauch, M., Otto, R., Volk, M., 2015. HRU Aggregation and Its Effects on Model Outputs, 2015 International SWAT Conference, Sardinia, June 24–25. Toth, B., Pietroniro, A., Conly, F.M., Kouwen, N., 2006. Modelling climate change impacts in the Peace and Athabasca catchment and delta: I—hydrological model application. Hydrol. Process. 20 (19):4197–4214. http://dx.doi.org/10.1002/hyp.6426. Troin, M., Caya, D., 2014. Evaluating the SWAT's snow hydrology over a Northern Quebec watershed. Hydrol. Process. 28 (4):1858–1873. http://dx.doi.org/10.1002/hyp.9730. Wang, S., et al., 2012. Multi-site calibration, validation, and sensitivity analysis of the MIKE SHE model for a large watershed in northern China. Hydrol. Earth Syst. Sci. 16 (12):4621–4632. http://dx.doi.org/10.5194/hess-16-4621-2012. Winchell, M., Srinivasan, R., Luzio, M.D., 2010. ArcSWAT Interface for SWAT2009 User's Guide. Soil and Water Research Laboratory-Agricultural Research Service, Blackland, Texas. WMO, 2012. Standardized Precipitation Index User Guide. World Meteorological Organization, Geneva.