Spatial and Temporal Change in the Hydro-Climatology of the Canadian Portion of the Columbia River Basin under Multiple Emissions Scenarios

Abstract

The Canadian portion of the Columbia River basin, or Upper Columbia River basin (UCRB), which is made up of the sub-basins of the Columbia and Kootenay rivers, has a hydrologic regime dominated by snow and glacier melt that makes it particularly susceptible to climate change. To investigate the influence of climate change on the water balance and streamflow of the UCRB, simulations are made for a 20-member ensemble made up of seven global climate models (GCMs) run under three emissions scenarios (A1B, A2, and B1 (with one exception)) for three future periods (2020s, 2050s, and 2080s) downscaled using Bias Corrected Spatial Disaggregation and simulated via the spatially distributed Variable Infiltration Capacity (VIC) model. To maximize the utility of the VIC model, a process-based distributed model, it is calibrated to 24 sub-basins. Snow water equivalent (SWE) decreases at lower elevations primarily in the southern portion of the UCRB and increases progressively at high elevations in the north from the 2020s to the 2050s to the 2080s based on the mean of seven GCMs run under the A2 scenario. Evaporation increases where adequate moisture is available in the summer months out to the 2080s. Runoff increases strongly in winter and spring and increases moderately in fall in the UCRB. Mean annual streamflow is projected to increase twice as fast in the northern, higher elevation, Columbia River basin than in the southern Kootenay River basin. Projected increases in mean annual flow are greatest under the A2 emissions scenario, ranging from a 9% increase in the 2020s to a 27% increase in the 2080s in the Columbia River at Keenlyside Dam. Changes in the hydrograph for four sub-basins of the UCRB are most pronounced in the 2080s under the A2 emissions scenario with streamflow increases projected in November to May, decreases in July to September, and mixed responses in June and October. The magnitude of peak flows and low flows is projected to increase for most sub-basins, time periods, and scenarios. Such projected changes could be expected to have implications for water resources management in the basin.

RÉSUMÉ [Traduit par la rédaction] La partie canadienne du bassin du fleuve Columbia, ou bassin du Haut–Columbia, qui comprend les sous-bassins du Columbia et de la rivière Kootenay, possède un régime hydrologique dominé par la fonte de la neige et des glaciers, ce qui le rend particulièrement sensible au changement climatique. Pour étudier l'influence du changement climatique sur le bilan hydrologique et l’écoulement fluvial dans le bassin du Haut–Columbia, nous avons fait des simulations avec un ensemble de 20 membres issus de sept modèles climatiques mondiaux exécutés sous trois scénarios d’émissions (A1B, A2 et B1 (avec une exception)) pour trois périodes futures (les décennies 2020, 2050 et 2080) avec une échelle réduite au moyen de la désagrégation spatiale corrigée du biais et simulées par le modèle VIC (Variable Infiltration Capacity) spatialement distribué. Pour maximiser l'utilité du modèle VIC, un modèle distribué fondé sur des processus, celui-ci est étalonné par rapport à 24 sous-bassins. L’équivalent en eau de la neige diminue aux élévations plus faibles, principalement dans la partie sud du bassin du Haut–Columbia, et augmente progressivement aux grandes élévations dans le nord, de la décennie 2020 à la décennie 2050 et à la décennie 2080, selon la moyenne de sept modèles climatiques mondiaux exécutés sous le scénario A2. L’évaporation augmente là où l'apport en humidité est adéquat durant les mois d’été jusque dans les années 2080. Le ruissellement augmente fortement en hiver et au printemps et augmente modérément à l'automne dans le bassin du Haut–Columbia. Il ressort que l’écoulement fluvial annuel moyen devrait augmenter deux fois plus vite dans le bassin du Columbia au nord, où les élévations sont grandes, que dans le bassin de la Kootenay au sud. Les accroissements prévus dans l’écoulement annuel moyen sont les plus grands sous le scénario d’émissions A2, allant d'un accroissement de 9% dans les années 2020 à un accroissement de 27% dans les années 2080 au barrage de Keenlyside sur le Columbia. Les changements dans l'hydrographe pour quatre sous-bassins du bassin du Haut–Columbia sont les plus prononcés dans les années 2080 sous le scénario d’émissions A2 avec des accroissements d’écoulement fluvial prévus de novembre à mai, des diminutions de juillet à septembre et des réponses mixtes en juin et en octobre. L'intensité des débits de pointe et d’étiage est prévue s'accroître dans la plupart des sous-bassins, des périodes et des scénarios. On peut s'attendre à ce que ces changements prévus aient des répercussions sur la gestion des ressources en eau dans le bassin.

Keywords:

• Columbia River
• Kootenay River
• streamflow projections
• climate change
• extremes
• SWE

1 Introduction

The Columbia River is the largest hydroelectric producing river in North America, and roughly 50% of the electricity for the province of British Columbia is produced by 18 hydroelectric generating facilities in the Canadian portion of the basin (CBT, 2012), the Upper Columbia River Basin (UCRB). The UCRB also provides the domestic water supply for its 145,000 inhabitants (CBT, 2012). Additionally, its water resources sustain forestry, agriculture, fisheries, and recreation in the basin. The operation of reservoirs on the Columbia River is coordinated between Canada and the United States to achieve flood control and power generation based on the precepts of the international Columbia River Treaty ratified by Canada in 1964 (Treaty Relating to Cooperative Development of the Water Resources of the Columbia River Basin, United States–Canada, 1961–1964); since the 1960s much has been invested to develop this resource. Changes to the hydrologic regime driven by climate change can be expected to impact UCRB water resources and demand.

Historical hydroclimatic change in western North America (WNA) has been attributed to anthropogenic climate change, predominantly from widespread regional warming (Barnett et al., 2008; Bennett, Werner, & Schnorbus, 2012; Pierce et al., 2008). Throughout most of British Columbia, including the UCRB, seasonal runoff is either dominated by snow or a hybrid of snow and glacial melt (Eaton & Moore, 2010). Within these types of regimes recent hydrologic trends include a widespread decline in snowpack, earlier onset of snowmelt, and a decrease in summer flows (Chang, Jung, Steele, & Gannett, 2012; Dery, 2009; Mote, 2006; Regonda, Rajagopalan, Clark, & Pitlick, 2005; Stewart, Cayan, & Dettinger, 2004). Recent climate projections indicate further warming is expected within WNA throughout the remainder of the twenty-first century at a rate that is also likely to exceed the rate of projected global-average warming (Christensen et al., 2007). Climate projections also indicate a consistent trend of increasing precipitation in three of four seasons (projections are ambiguous for the summer) in the region (Christensen et al., 2007). Further warming suggests that hydrologic trends observed over recent decades are likely to continue through the coming century (Elsner et al., 2010; Vano et al., 2010). Nevertheless, snow and glacier changes are sensitive to changes in both temperature and precipitation (Hamlet, Mote, Clark, & Lettenmaier, 2005; Serquet, Marty, Dulex, & Rebetez, 2011). Consequently, the possibility of increased precipitation, as well as increased temperatures, complicates the extrapolation of recent trends into the future. In addition, the magnitude of future hydroclimatic trends is affected by uncertainty between different global climate models (GCMs) and differences between scenarios of the evolution of future greenhouse gas emissions, especially later in the next century (Bennett et al., 2012).

The impact of climate change on the hydrologic regime of the Columbia River has been a subject of ongoing research. Hamlet and Lettenmaier (1999) projected increases in winter runoff, on average, across the entire Columbia River basin as a result of increased winter precipitation, warmer winter temperatures, and resulting reductions in snowpack using the Variable Infiltration Capacity (VIC) model at 1/8° resolution in combination with four GCMs from the first Coupled Model Intercomparison Project (CMIP). More recently, streamflow projections have been made using GCMs from the third CMIP (CMIP3). Schnorbus, Werner, and Bennett (2012) projected increased runoff in fall, winter, and spring and decreased discharge in summer, on average, over the UCRB based on the median change for eight CMIP3 GCMs for the A1B and A2 emissions scenarios and seven in the case of the B1 scenario. The regime of the Mica sub-basin of the UCRB was projected to remain dominated by a hybrid of snow and glacier melt with a shift to increased discharge in the winter and spring and decreased discharge in the summer and early fall for the 2050s under the A1B scenario. Bürger, Schulla, and Werner (2011) looked at streamflow changes, including extremes, in the 2050s for the Columbia River above Donald from four CMIP3 GCMs run under the A2 scenario. They projected that the hydrograph will shift towards a more rain-fed regime with peak flows occurring in June instead of July. Annual peak flows were projected to decrease moderately in three of the four models while August low flows were projected to decrease in all models. Hamlet (2011) provided an overview of the hydrologic impacts of climate change in the entire Columbia River basin for the 2020s, 2050s, and 2080s for the median of multiple CMIP3 GCMs run under B1 and A1B scenarios. He projected that, overall, the UCRB would remain dominated by snowmelt, while changes to hydrologic extremes would be spatially heterogeneous.

Other studies in the Pacific Northwest have found high spatial variability in hydrologic response to climate change (e.g., the Willamette basin in Oregon (Chang & Jung, 2010), Washington State watersheds (Elsner et al., 2010), and the Fraser River in British Columbia (Shrestha, Schnorbus, Werner, & Berland, 2012)). For example, changes in snow storage resulting from climate change are strongly affected by elevation; therefore, large spatial variations occur over regions of high relief (Hamlet et al., 2005; Kim, 2001; Knowles & Cayan, 2004; Mote, 2006; Mote, Hamlet, Clark, & Lettenmaier, 2005). The future changes could also have high temporal variability, as shown by previous studies using multiple future periods (Chang & Jung, 2010). Such spatial and temporal variability of the future hydrologic changes will have implications on general water availability in the UCRB such as for water supply and hydropower generation, and extremes such as flood and drought (Cunderlik & Simonovic, 2005; Hamlet, 2011; Hamlet & Lettenmaier, 2007; Lehner, Döll, Alcamo, Henrichs, & Kaspar, 2006).

Although much work has been done to explore projected changes to streamflow in the UCRB a few key pieces remain. Aside from Schnorbus et al. (2012), who studied changes in runoff from a single outflow point, no study has considered projections of streamflow in the UCRB on a detailed sub-basin scale. This work expands upon the overview provided by Schnorbus et al. (2012) by providing results for more sites, an analysis of extremes, and a more spatially explicit look at annual and monthly runoff changes for three future time periods. The objective of this work is to illustrate the spatial and temporal variability of future projections of hydrologic changes in the UCRB. Building on the previous study by Schnorbus et al. (2012), an ensemble of streamflow projections is provided for a range of representative sub-basins of the UCRB for the 2020s, 2050s, and 2080s including monthly and annual means, 7-day low flows, and annual daily maxima. Multiple GCMs from CMIP3, run under three emissions scenarios, are used to characterize the range of uncertainty.

2 Study area

The UCRB includes the 104,000 km² drainage area upstream of the confluence of the Kootenay and Columbia rivers. It is located predominantly in southwestern British Columbia but also includes northern Idaho and northwestern Montana (Fig. 1). The study area is centred on the Columbia Mountains and is bordered by the Rocky Mountains to the east and the Shuswap–Okanagan Highlands to the southwest. Elevation ranges from 420 m to 3550 m, with a median elevation of 1600 m. The UCRB has a continental climate (Demarchi, 1996) with an annual average temperature of 1.9°C and monthly average values that range from −9.4°C in January to 13.4°C in July based on gridded observations from 1961 to 1990 (see Section 3). Precipitation rates display some slight seasonal variation and are highest in the winter. Approximately 65% of the annual precipitation falls as snow, with snowfall possible throughout the year in some parts of the basin. Consequently, rivers in the UCRB are snow dominated with a prominent spring freshet followed by a relatively steep recession transitioning into low-water availability in the late summer and early fall. In addition, approximately 4% of the area is covered by glaciers such that natural streamflow at many locations in the UCRB exhibits a glacial–nival regime, where the hydrologic cycle is dominated by the spring freshet, with a gradual recession in flow during the late summer and fall and with the lowest flow occurring during the winter.

Fig. 1 Study area of the UCRB. The numbers 1–5 are the sub-basin outlets defined in Table 1.

The study area has been divided into the Columbia River and Kootenay River sub-basins to demonstrate the diversity across the UCRB (Fig. 1). The outlets of these sub-basins correspond with two BC Hydro project sites: Keenlyside Dam (at the outlet of Arrow Lakes; BCHAR), and Kootenay Canal (BCHKL), for the Columbia and Kootenay sub-basins, respectively (Table 1). Two additional sites have been selected for investigation in each sub-basin, each corresponding to Water Survey of Canada (WSC) gauge sites. These smaller headwater stations are the Columbia River at Donald (COLAD) and Elk River at Fernie (ELKAF) (Table 1). Response at the two BC Hydro project sites demonstrates the cumulative effect of temperature and precipitation changes on discharge and other water balance components over a relatively large area, whereas response at the small nested headwater stations illustrates the impact of climate change on natural flows at higher elevations and more localized settings. The higher elevation sub-basins have colder annual temperatures and less precipitation on average than the larger basins. Note that although BCHAR (36,659 km²) has a smaller drainage area than BCHKL (46,398 km²), BCHAR has more runoff (890 mm) on average over its basin than BCHKL (515 mm) and a higher runoff ratio (0.72 versus 0.61), defined as runoff divided by precipitation (R/P) indicating that evaporation rates in BCHKL are higher. The BCHAR drainage area contains more area at high elevation and has the majority of the glacier cover (63% of total glacier area) in the UCRB, whereas the BCHKL sub-basin has a smaller range in elevation and less glacier cover (37% of total glacier area).

Table 1. Study area physical description and 1961–90 annual hydroclimatology.

	Basin
Meta Data	UCRB	BCHAR	COLAD	BCHKL	ELKAF
Site Number	1	2	3	4	5
Water Survey of Canada #	NA	NA	08NB005	NA	08NK002
Project Site or Station Name		Columbia River at Keenlyside Dam	Columbia River at Donald	Kootenay River at Kootenay Canal	Elk River at Fernie
Flow Regime	Regulated	Regulated	Natural	Regulated	Natural
Drainage Area (km^2)	104,000	36,659	9710	46,398	3110
Minimum Elevation (m)	400	423	771	497	995
Median Elevation (m)	1606	1716	1845	1549	1900
Maximum Elevation (m)	3552	3552	3433	3483	3186
Annual Precipitation (mm)	1016	1223	803	840	671
Annual Temperature (°C)	1.6	0.6	0.1	2.3	0.6
Annual Runoff (mm)	686	890	550	515	379
Annual Evaporation (mm)	305	278	258	322	286

3 Methods

Our study applies GCM output contributed by the World Climate Research Programme through the CMIP3 multi-model dataset (Meehl et al., 2007). A subset of seven of the 25 participating models was selected with the goals of 1) using the most robust models, based on model performance over the globe and WNA and 2) sampling the range of plausible projections of climate change represented in CMIP3 in a reasonably comprehensive manner (Werner, 2011). Simulations were forced by the B1, A1B, and A2 emissions scenarios described in the Intergovernmental Panel on Climate Change (IPCC) Special Report on Emissions Scenarios (SRES) (Nakićenović & Swart, 2000), which were selected to look at the low, medium, and high ranges of possible future emissions, respectively. A simulation forced by the B1 emissions scenario was not available for HadGEM1; therefore, seven of the selected GCMs were run under A1B and A2 and six were run under B1 for a total of 20 climate projections, as summarized in Table 2. By using an ensemble of simulations from multiple models under three emissions scenarios we explicitly address both emissions and GCM structural uncertainty to some extent. Recent research suggests that differences in global climate response to greenhouse forcings between different GCMs is the largest source of uncertainty, and uncertainties attributed to downscaling, emissions scenarios, and hydrologic modelling are of lesser magnitude (Bennett et al., 2012; Bloschl & Montanari, 2010; Prudhomme & Davies, 2009).

Table 2. Global climate model and SRES scenario selection summary.

Model Abbreviation	Modelling Centre, Full Model Name, and Model ID	SRES Scenarios	Primary Reference
CGCM3	Canadian Centre for Climate Modelling and Analysis (Canada), Third Generation Atmospheric General Circulation Model, CGCM3.1 (T47)	A2, A1B, B1	Scinocca, McFarlane, Lazare, Li, and Plummer (2008)
CCSM3	National Center for Atmospheric Research (USA), Community Climate System Model, version 3.0, CCSM3	A2, A1B, B1	Collins et al. (2006)
ECHAM	Max Planck Institute for Meteorology (Germany), Fifth Generation European Centre Hamburg Model (ECHAM5) and Max Planck Institute Ocean Model (MPI-OM), ECHAM5/MPI-OM	A2, A1B, B1	Roeckner et al. (2006)
GFDL2.1	US Dept. of Commerce/NOAA/Geophysical Fluid Dynamics Laboratory (USA), Geophysical Fluid Dynamics Laboratory-Coupled Climate Model 2.1,GFDL-CM2.1	A2, A1B, B1	Delworth et al. (2006)
HadCM	UKMO Hadley Centre for Climate Prediction and Research (UK), Hadley Centre Coupled Model version 3, HadCM3	A2, A1B, B1	Collins, Tett, and Cooper (2001)
HadGEM1	UKMO Hadley Centre for Climate Prediction and Research (UK), Hadley Centre Global Environmental Model version 1, HadGEM1	A2, A1B	Martin et al. (2006)
MIROC3.2	Center for Climate Systems Research (The University of Tokyo), National Institute for Environmental Studies, and Frontier Research Center for Global Change (JAMSTEC), Model for Interdisciplinary Research on Climate version 3.2 medium-resolution, MIROC3.2 (medres)	A2, A1B, B1	K-1 model developers, Hasumi and Emori (2004)

Daily time series of gridded temperature and precipitation were generated for the period 1950 to 2098 for each GCM emissions scenario combination run at the 1/16° resolution of the hydrologic model using the Bias Corrected Spatial Disaggregation (BCSD) statistical downscaling technique (Werner, 2011). The BCSD technique was used because of its computational efficiency (Salathé, 2005; Wood, Maurer, Kumar, & Lettenmaier, 2002; Wood, Leung, Sridhar, & Lettenmaier, 2004) and effectiveness in projecting changes to temperature and precipitation that is on a par with several more complex techniques (Bürger, Murdock, Werner, Sobie, & Cannon, 2012; Maurer, Hidalgo, Das, Dettinger, & Cayan, 2010). Monthly GCM climate projections of temperature and precipitation were spatially downscaled via bias correction of monthly GCM precipitation and temperature onto gridded observed data (aggregated to the scale of the GCM) using quantile mapping over the 1950 to 1999 period. Temporal disaggregation from monthly to daily data is accomplished via random sampling of historical months, where each day in the selected month is rescaled (multiplicative for precipitation and additive for temperature). Note that this procedure assumes that the daily statistics (frequency, duration of wet spells, duration of hot spells, etc.) of precipitation and temperature do not change in projected future climates.

Gridded historical climate data are the basis for calibrating the hydrologic model and the statistical downscaling approach. Daily gridded surfaces of maximum and minimum temperature, daily precipitation accumulation, and daily average wind speed were produced following the techniques of Maurer, Wood, Adam, Lettenmaier, and Nijssen (2002) and Hamlet and Lettenmaier (2005) at a spatial resolution of 1/16° (see Schnorbus, Bennett, Werner, & Berland, 2011). The technique involves spatial interpolation of daily temperature and precipitation. Station data were gathered from multiple networks including those of Environment Canada, the British Columbia Ministry of Forests, Lands and Natural Resource Operations, BC Hydro, and the US National Weather Service Co-operative Observer Program. The raw gridded fields were temporally homogenized to remove interpolation artifacts introduced by using a temporally varying mix of stations and corrected for topographic effects using Climate WNA, a 1961 to 1990 high-resolution climatology for western Canada based on the Parameter-elevation Regressions on Independent Slopes Model (PRISM) (Daly, Neilson, & Phillips, 1994; Wang, Hamann, Spittlehouse, & Aitken, 2006).

Hydrologic projections were derived using the VIC model (Liang, Lettenmaier, & Wood, 1994; Liang, Wood, & Lettenmaier, 1996). The VIC model is a spatially distributed macroscale hydrologic model that was originally developed as a soil-vegetation atmosphere transfer scheme (SVATS) for general circulation models. It has been used to evaluate climate change impacts on global river systems (Nijssen, O'Donnell, Hamlet, & Lettenmaier, 2001) and in the mountainous western United States and British Columbia (Elsner et al., 2010; Hamlet & Lettenmaier, 1999; Hamlet et al., 2005; Schnorbus et al., 2012; Shrestha et al., 2012). Its spatially distributed (gridded) nature makes it suitable for capturing regional variation in the hydrologic cycle because of variation in topographic, physiographic, and climatic controls. The VIC model is a process-based model allowing for a more plausible extrapolation of hydrologic processes into future climate regimes (Leavesley, 1994; Ludwig et al., 2009).

The VIC model was applied at a resolution of 1/16° (approximately 27–31 km², depending on latitude) and run at a daily time step (one-hour time step for the snow model). Subgrid-scale topographic variability was parameterized using up to five elevation bands per grid cell, with elevation data derived from a post-processed version (version 3) of the Shuttle Radar Topography Mission 90 m resolution digital elevation model (Jarvis, Reuter, Nelson, & Guevara, 2006), obtained from the Consultative Group on International Agricultural Research–Consortium for Spatial Information (http://srtm.csi.cgiar.org/). Input values of temperature and precipitation for each grid cell were interpolated to each elevation band. Precipitation was adjusted for elevation using a precipitation gradient derived from a 4 km resolution 1961–90 PRISM climatology (Daly et al., 1994) downscaled to 400 m resolution using Climate WNA (Wang et al., 2006). Temperature was lapsed to each elevation band using a rate of 6.5°C km⁻¹, applied over the difference between the mean band elevation and the mean grid cell elevation. Soil classification and parameterization for each VIC cell was interpolated from the 1/12° Soils Program in the Global Soil Data Products CD-ROM (Global Soil Data Task, 2000; Robock et al., 2000). Subgrid-scale variability in land cover was assigned by partitioning each grid cell into one or more land cover types based on 25 m resolution land cover from the Canadian Forest Service's Earth Observation for Sustainable Development of Forests dataset (Wulder et al., 2003). The forest cover parameterization represents circa 2000 conditions throughout the simulations to 2099. It is recognized that changes in vegetation related to climate change and forest management practices are likely to occur over the relatively long time horizon of these projections (e.g., Carroll et al., 2006). However, it is beyond the scope of the current project to implement dynamic vegetation. Greater details on the VIC model set-up for the UCRB can be obtained from Schnorbus et al. (2012). Surface routing between grid cells was done using the linearized Saint-Venant equations through a network specified by a flow direction and distance for each 1/16° grid cell (Lohmann, Nolte-Holube, & Raschke, 1996).

Terrain variability plays a major role in the hydrology of mountainous regions, particularly in the way that it affects the spatial variability of snow accumulation and ablation. Variation in relief is a major source of spatial variability (because of its effects on temperature and precipitation) and is specifically addressed using subgrid elevation bands. The VIC model, however, does not explicitly incorporate the effects of other terrain variables, for example, slope, aspect, and shading, that govern the local energy balance and the redistribution of snow by wind and gravity (Elder, Dozier, & Michaelsen, 1991; Lehning, Löwe, Ryser, & Raderschall, 2008). By using the VIC model one implicitly assumes that these ignored local-scale terrain effects cancel out when integrated over the grid-cell scale. Schnorbus et al. (2012) have shown that this approach is reasonably effective at replicating peak seasonal snow water equivalent (SWE) measurements and capturing regional variation in SWE with elevation in the UCRB. Nevertheless, ignoring several components of spatial variation driven by subgrid terrain remains a source of uncertainty and simulated results are not expected to reflect precisely the local reality of snow accumulation, ablation, evaporation, soil moisture, and runoff at the grid-cell scale (27–31 km²). However, this is not considered a serious limitation because the spatial results for snow, evaporation, and soil moisture are not intended to be interpreted on a cell-by-cell basis but rather from a regional perspective.

Because glaciers are present in the UCRB, in this implementation the VIC model is initialized in a modified way to simulate glaciers. A simple conceptual representation of glacier mass balance is used that models glacier mass using perennial snow in conjunction with the VIC model's built-in snow modelling routines (Schnorbus et al., 2012). A glacier state was estimated for circa 1995 using information from the 1:250,000 Landsat 5-based Baseline Thematic Mapping (version 1) land cover dataset (Government of British Columbia, 2010). Because of the lack of a historical glacier inventory, this glacier state was used to initialize the VIC transient runs starting in 1950 and then to reinitialize glacier state information in 1995. It is recognized that this process does not explicitly capture the historical trajectory of glacier area and volume within the UCRB. However, the work of Debeer and Sharp (2007), which investigated changes in glacier area between the years 1951/1952 and 2001 in southeastern British Columbia, indicates a trend in glacier area of approximately −220 km² in the UCRB region. Tennant, Menounos, Wheate, and Clague (2012) estimated a trend in glacier area of approximately −306 km² between 1950 and 1995 in the Canadian Rocky Mountains (estimated from a trend of −6.8 km² y⁻¹ from 1919 to 2006). Consequently, the use of a 1995 glacier mask underestimates 1950 glacier area in the UCRB by roughly only eight to eleven grid cells, assuming 28 km² per cell. It should be noted that glacier dynamics are not explicitly represented; regionally, simulated changes in glacier area occur only through the local loss (or emergence) of perennial snow within individual glacier grid cells (and elevation bands within grid cells). Further discussion on the implementation and limitations of this approach can be found in Schnorbus et al. (2012).

Calibration of the UCRB application of the VIC model is described in detail by Schnorbus et al. (2012) and is only summarized herein. Initial calibration involved manual adjustment of parameters controlling snow accumulation and melt, and based on an analysis of 1 April snow course data, the VIC model simulates the regional variation of snow accumulation in the UCRB with reasonable accuracy (Schnorbus et al., 2012). To take advantage of its distributed nature, the VIC model was subsequently calibrated for 24 streamflow locations corresponding to WSC hydrometric stations and BC Hydro project sites (Fig. 1). The VIC model was calibrated using the automated Multi-Objective Complex Evolution (MOCOM) method (Yapo, Gupta, & Sorooshian, 1998). With some exceptions, because of data limitations (see the notes for Table 3), calibration and validation at all sites used daily discharge data for the periods 1990 to 1995 and 1985 to 1989, respectively (these two periods represent the period of highest concurrent streamflow data density in the UCRB). Observed data were used at the unregulated WSC stations, whereas naturalized discharge was used for regulated sites (BC Hydro, unpublished data, 2009). Nash–Sutcliffe (NS) values for the calibration period range from 0.65 to 0.99, with a mean of 0.83 over the 24 calibration basins (Table 3). Performance across calibration sites is similar overall, for the Nash–Sutcliffe coefficient of log transformed flow (LNS), with a domain average of 0.80 and values ranging from 0.27 to 0.99. The percentage mean volume bias (%VB) errors during the calibration period range from −15% to 10%, with a mean of −4%, and the majority of absolute %VB values (23 of 24 sites) are less than or equal to 10%. As might be expected, statistics are slightly lower overall during the validation period than during the calibration period because of some overfitting to specific conditions experienced during the calibration period. Nevertheless, the VIC model parameterizations seem robust because basin by basin, the NS, LNS, and %VB performance measures do not degrade substantially during the validation period (Table 3). No observed data are available to compare evaporation and soil moisture directly to modelled results; hence, these components of the hydrologic cycle were not validated explicitly. Nevertheless, our evaporation results are broadly consistent with those of other global (Vinukollu, Meynadier, Sheffield, & Wood, 2011) and regional (Liu, Chen, & Chihlar, 2003) modelling efforts.

Table 3. Summary of calibration and validation results for Columbia River sub-basins for three performance measures. NS is the Nash–Sutcliffe values, LNS is the Nash–Sutcliffe coefficient of log-transformed discharge, and %VB is the percentage mean volume bias. BC Hydro project sites are indicated with bold text.

Basin	Calibration 1990–1994			Validation 1985–1989
	NS	LNS	%VB	NS	LNS	%VB
BCHAR ^a ^, ^b	0.78	0.84	−2	0.80	0.65	−8
BCHDN ^a	0.65	0.54	−7	0.73	0.59	−8
BCHKL ^a ^, ^c	0.67	0.60	−7	0.72	0.74	4
BCHLB ^a ^, ^d	0.97	0.80	−4	—	—	—
BCHMI ^a	0.89	0.83	−9	0.88	0.79	−7
BCHRE ^a ^, ^b	0.97	0.97	−4	0.92	0.81	−10
BCWAT ^a	0.76	0.75	−10	0.74	0.67	−12
BLAAW	0.72	0.86	10	0.75	0.87	−3
BULNW	0.83	0.77	−4	0.81	0.72	−17
COLAD	0.94	0.94	−2	0.91	0.88	−1
ELKAF	0.99	0.97	−3	0.81	0.69	−13
ELKNN	0.89	0.92	−3	0.75	0.77	−8
FORAM	0.66	0.70	−6	0.72	0.74	−1
KICAG	0.77	0.87	−8	0.77	0.86	−2
KOOAF	0.99	0.99	3	0.85	0.80	−4
KOOAK	0.78	0.77	−4	0.75	0.68	17
KOOCF	0.91	0.91	−4	0.84	0.86	−4
KOONS	0.98	0.99	0	0.85	0.83	−6
MICBN	0.75	0.80	−5	0.76	0.79	−15
PALIL	0.80	0.77	−4	0.80	0.80	−10
SALNS	0.74	0.27	−15	0.73	0.38	−12
SLONC	0.78	0.66	−2	0.78	0.72	−2
SPINS ^e	—	—	—	0.67	0.82	−1
STMAW	0.99	0.99	−2	0.84	0.63	−16
STMNM	0.76	0.46	−10	0.82	0.47	−20
Average	0.83	0.80	−4	0.79	0.73	−7
^aCalibrated and validated to naturalized discharge.
^bValidation based on monthly streamflow.
^cValidation based on monthly streamflow from 2003 to 2006.
^dCalibrated based on 2003 to 2006 period.
^eUncalibrated; validation based on 1980 to 1984 period.

Daily observed (naturalized) and modelled streamflow from 1985 to 1994, which includes the validation and calibration periods, are presented for the example basins Elk River at Fernie (ELKAF), Columbia River at Donald (COLAD), Kootenay River at Kootenay Canal (BCHKL) and Columbia River at Keenlyside Dam (BCHAR). These four sites were selected to represent the more glacier-influenced northwestern domain of the UCRB (COLAD and BCHAR) and the snow-dominated southeastern portion of the UCRB (ELKAF and BCHKL). The BCHAR and BCHKL basins drain the entire Upper Columbia and Kootenay rivers, respectively, near the international border between British Columbia and Washington State. The longer recession in the more glacier-dominated BCHAR and COLAD is apparent in the calibration and validation periods (Fig. 2). Mid-winter events are displayed in the observed flow and captured in the simulated flow for BCHKL. These events are over-simulated for ELKAF. Calibration values for NS and LNS are highest for COLAD and ELKAF.

Fig. 2 Modelled daily streamflow compared with observed (COLAD and ELKAF) or naturalized (BCHAR and BCHKL) streamflow over the calibration and validation periods.

Extreme flow statistics were calculated from daily streamflow for the baseline (1970s), 2020s, 2050s, and 2080s. The Generalized Extreme Value (GEV) distribution was fit to annual daily maximum and 7-day minimum streamflow time series for each time period and basin by a maximum likelihood parameter estimation (Hosking & Wallis, 1993; Katz, Parlange, & Naveau, 2002) using the “ismev” package for R (Heffernan, 2012). The 7-day consecutive low flow with a return frequency of 10 years (7Q10) and the 25-year return period peak-flow value were estimated for each basin and future period from the fitted GEV distribution. Additionally, the peak-flow dates for each 30-year period and basin were extracted.

4 Results and discussion

a Temperature and Precipitation Projections

Results in this section report the ensemble mean of the mean for each scenario as the consensus estimate of projected temperature and precipitation changes for each time period and season. Annual and seasonal temperature increases are projected for all time periods, under all emissions scenarios (Table 4). Annual temperature in the UCRB in the 1970s was 1.6°C based on gridded observations. Under the A2 emissions scenario, annual temperatures are projected to increase from 1.2°C to 1.6°C per time period for the 2020s, 2050s, and 2080s. By the 2080s, mean annual temperature is projected to be 6.0°C averaged over the UCRB. Increases are similar under A1B in the 2020s and 2050s. Temperature increases are projected for all seasons, for all time periods and emissions scenarios at a rate of roughly 1°C per 30-year period, except for summer when increases are closer to 2°C (Table 4). Although winter temperatures are projected to increase, the ensemble means remain below freezing over the UCRB in all scenarios. The largest difference in projected change is between A2 and B1 in the 2080s in summer when increases were 5.7°C under A2, but only 3.2°C under B1 (Table 4). Precipitation is projected to increase annually under all time slices and emissions scenarios, for each 30-year period to a maximum of a 15% increase in the 2080s versus 1970s values for the A2 emissions scenario (Table 5). Projected increases in annual precipitation in the 2050s are largest for the B1 scenario and largest for the 2080s under the A2 scenario. Precipitation is projected to increase monotonically for each 30-year time period in winter (DJF), spring (MAM), and fall (SON) to a maximum of 25% in the 2080s under A2 and progressively decrease in summer (JJA) by up to 17% in the 2080s under A2 (Table 5). Thus, all seasons are projected to become warmer and wetter, except summer which is projected to become warmer and drier.

Table 4. Mean annual (ANN) and seasonal (DJF, MAM, JJA, and SON) temperatures by time period (1970s, 2020s, 2050s, and 2080s) and emissions scenarios (A2, A1B and B1).

	Units	Historical	A2			A1B			B1
		1970s	2020s	2050s	2080s	2020s	2050s	2080s	2020s	2050s	2080s
ANN	(°C)	1.6	3.0	4.2	6.0	3.0	4.5	5.5	2.9	3.6	4.3
DJF	(°C)	−8.6	−7.2	−6.3	−4.3	−7.1	−6.0	−5.0	−7.2	−6.5	−5.9
MAM	(°C)	1.2	2.3	3.3	4.7	2.5	3.6	4.4	2.3	3.0	3.7
JJA	(°C)	12.0	13.7	15.3	17.7	13.9	15.7	16.9	13.4	14.4	15.2
SON	(°C)	1.9	3.2	4.4	6.0	3.0	4.7	5.6	3.0	3.7	4.4

Table 5. Mean annual and seasonal precipitation by time period (1970s, 2020s, 2050s, and 2080s) and emissions scenarios (A2, A1B, and B1).

	Units	Historical	A2			A1B			B1
		1970s	2020s	2050s	2080s	2020s	2050s	2080s	2020s	2050s	2080s
ANN	(mm)	1033	1074	1099	1187	1079	1121	1161	1077	1131	1141
DJF	(mm)	351	373	391	437	374	403	413	365	393	403
MAM	(mm)	204	221	232	253	220	229	245	217	234	233
JJA	(mm)	218	203	192	180	200	187	187	212	209	201
SON	(mm)	261	278	286	318	285	303	317	284	297	305

On an annual basis, basin-averaged temperature increases projected for the mid- and late twenty-first century over the Columbia River basin above the Dalles, Oregon, by Hamlet and Lettenmaier (1999) are similar to those projected for the UCRB in this study. Projected precipitation increases in winter are similar between the two studies, but projected decreases in summer precipitation are greater in this study. Discrepancies likely occur because Hamlet and Lettenmaier (1999) analyzed changes in precipitation and temperature over a larger domain, using an older subset of the GCMs investigated here, even though they used the same statistical downscaling approach applied in this study. Annual precipitation and temperature increases projected in this study were greater than those found in neighbouring Washington state (Elsner et al., 2010). Differences likely result from the application of a different downscaling method, the delta method, and the use of slightly different GCMs from CMIP3.

b Spatial Response under the A2 Emissions Scenario

To limit repetition and maintain brevity, results and discussion in this section focus only on the A2 scenario. Because recent emissions have already exceeded the most pessimistic SRES projection (Raupach et al., 2007), we feel that the A2 scenario currently reflects the most likely projection of emissions at the end of this century. Hence, results for this scenario are used to exemplify the spatial response of the UCRB to climate change. Nevertheless, results for the A1B and B1 scenarios are qualitatively similar. Again, discussion throughout this section uses the A2 ensemble mean as the consensus estimate of future changes.

The effect of temperature changes on snowpack is assessed using the ratio of peak (i.e., 1 April) SWE to total cold-season precipitation (October through March), P_w , (Barnett et al., 2008). Using this ratio, hydrologic regimes are classified as nival (SWE _p /P_w  > 0.5), hybrid pluvial-nival (0.1 ≤ SWE _p /P_w ≤0.5) and pluvial (SWE _p /P_w  < 0.1) (Elsner et al., 2010). In the historical period, precipitation falls as snow over most of the UCRB, and the majority of the area is classified as nival (Fig. 3). High-elevation areas of perennial snow and ice, where SWE _p /P_w is in excess of 1.0, are consistent with documented glacier areas (Bolch, Menounos, & Wheate, 2010). Through the 2020s, 2050s, and 2080s, there is a progressive decline in SWE _p/P_w over the study area, signalling a gradual transition from nival to hybrid to pluvial regimes in large parts of the UCRB. A transition to pluvial regimes is especially notable in valley bottom areas, particularly in the BCHKL basin. Nevertheless, a substantial proportion of the UCRB is at a high enough elevation (where winter temperatures remain below freezing), so that the majority of winter precipitation is projected to continue to fall as snow, even into the 2080s (Fig. 3). The area of SWE _p/P_w  > 1 decreases as we move through the 2020s, 2050s, and 2080s, but large areas remain, particularly in BCHAR. However, the extent of glacier cover remaining in the future should be interpreted with caution because Schnorbus et al. (2012) point out that the persistence of glaciers (SWE _p /P_w  > 1) throughout the 2020s, 2050s, and 2080s is likely overestimated by the VIC model. The lack of glacier dynamics in the current VIC model application to the UCRB results in unrealistic accumulations of snow and ice at high elevations because the snow and ice are not subject to gravitational redistribution to lower elevations, where they would be subject to higher melt rates. Nevertheless, the effect of the subtle interplay between precipitation and temperature changes on snow and ice accumulation at high elevations is aptly illustrated by the current results.

Fig. 3 Mean (A2) 1 April SWE versus October to March precipitation and 1 April SWE for 1961–90 and changes for 2011–40, 2041–70, and 2071–2100 for the mean of the A2 scenario.

Change in 1 April SWE is a function of both precipitation and temperature changes. Through the 2020s, 2050s, and 2080s, there is a pattern of decreased 1 April SWE at low elevation despite higher precipitation (because of higher temperatures) and increased SWE at higher elevation (because of increased precipitation) (Fig. 3). Nevertheless, with progressive warming from the 2020s through to the 2080s under the A2 emissions scenario, the area of decreased SWE expands (with a rising snowline), the magnitude of the SWE decrease becomes larger, and the region that experiences increased SWE becomes smaller (i.e., is confined to higher and higher elevations). By the 2080s, SWE decreases are projected to be as much as 100% in some areas, particularly in the BCHKL sub-basin.

Evaporation, as described here, includes evaporation from the canopy, transpiration, evaporation from bare surfaces, sublimation from the canopy, and sublimation from snow. Negative evaporation represents condensation, which occurs when warm air flows over cold snow or ice surfaces. Summer evaporation ranges from −10 mm to 400 mm over the basin over the historical 1970s period (Fig. 4). In the 2020s, evaporation is projected to decrease in the southern portion of the BCHKL basin and increase at high elevations and in the more northern BCHAR basin, based on the mean of seven GCMs run under the A2 emissions scenario. This general trend continues into the 2050s, when evaporation decreases at low elevations and in valley bottoms and increases at higher elevations. Projected changes in summer evaporation in the 2080s are similar in spatial pattern to those for the 2050s, but higher in magnitude with projected changes ranging from −80 mm to +80 mm (Fig. 4). In the southern part of the UCRB, predominantly in the BCHKL sub-basin, evaporation tends to decrease quite uniformly over the region. However, as one moves further north (most of the BCHAR and northern parts of the BCHKL) evaporation changes tend to correlate with the magnitude of historical summer evaporation, in that evaporation tends to increase where historical summer evaporation is low (less than 150 mm; typically at higher elevations) and tends to decrease where historical evaporation is high (greater than 150 mm; typically at lower elevations).

Fig. 4 Summer evapotranspiration and monthly soil moisture flux (June, July, and August) for 1961–90 and changes for 2011–40, 2041–70, and 2071–2100 for the mean of the A2 scenario.

Historically, soil moisture ranges from 0 to 700 mm in June, July, and August (Fig. 4). Low soil moisture, which is generally a function of soil depth, which decreases with increasing elevation, is generally found at high elevations and high soil moisture is prevalent at low elevations in June and July (Schnorbus et al., 2011). In August, the valley bottoms become dry, along with most of the Kootenay River basin. In early summer (June), conditions are projected to become wetter in the north at high elevations and drier in the south at low elevations. Each consecutive 30-year period has a similar pattern of change but with greater amplitude. Projected increases are greater than 95 mm in the north by the 2080s. Drying becomes much more prevalent in July, with the exception of some high elevation areas in the northern Columbia River basin, and is almost ubiquitous in August when some areas are projected to dry by more than 95 mm.

In future, potential evapotranspiration (PET) is projected to increase with increasing temperature (not shown), with actual evapotranspiration (AET) constrained by changes in moisture availability. Thus, those areas where summer (JJA) evaporation is projected to increase match those areas where early summer (June) soil moisture is also projected to increase (Fig. 4). Similarly, those areas that are projected to experience soil moisture declines of larger absolute magnitudes in August tend to also correlate with areas experiencing increased evaporation over the summer. This indicates that the northern portion of the basin will experience higher AET with increasing temperatures (i.e., energy limited), whereas the southern portion of the UCRB will experience less evaporation with drier initial summer soil conditions (i.e., water limited).

The spatial variation in projected runoff changes is assessed by comparing seasonal and annual runoff changes through the drainage network for the 2080s using the A2 ensemble mean, with results given in Fig. 5. Results are given only for channel networks of second order or higher (which represents a drainage area threshold of approximately 300 km²). This is considered a sufficient scale over which to integrate runoff accurately over individual grid cells. In winter, runoff is projected to increase throughout the UCRB. All increases are greater than 100%, except in the headwaters of the Kootenay River where increases range from 50% to 100%. These increases are attributed to increased precipitation and warmer temperatures causing increased snowmelt and more precipitation to fall as rain (see Section 4c). In spring, runoff increases of more than 100% are projected for the majority of the UCRB, including the main branch of the Columbia River. However, some tributaries in the far south of the UCRB show decreased runoff (between 0 and −10%). Overall, runoff in the Kootenay River increases on the order of 50 to 100%. Warmer temperatures result in a decreased SWE _p/P_w ratio (Section 4b) and earlier snowmelt, causing runoff to increase over spring (see Section 4c). The summer season experiences greater spatial variability in projected runoff changes. With the exception of one tributary, runoff is projected to decrease in all stream reaches upstream of BCHKL, with changes ranging from 0 to −50%, including a decrease on the order of 25 to 50% on the Kootenay River at the BCHKL site. Low soil moisture, decreased evaporation, and decreased snowmelt in the summer months (see Section 4c) result in an overall reduction in water availability, especially in the Kootenay River basin. Although the upper reaches of the Columbia River (predominantly above COLAD) also experience decreased summer runoff on the order of 10 to 50%, the majority of the main stem experiences increased summer runoff. However, once the river reaches BCHAR, runoff is again switched to a regime of decreased runoff on the order of 0 to 10%. As the Columbia River already receives greater summer runoff than the Kootenay River during the historical period, the larger projected decrease in summer runoff in the Kootenay River relative to the Columbia River suggests that this disparity will increase in the future.

Fig. 5 Multivariate drainage plot of seasonal and annual change in runoff (%) for the 2080s versus the 1970s (line colour) and 1970s seasonal and annual runoff relative to the study area total (line thickness) for the mean of seven GCMs run under the A2 scenario. Note that the modelled drainage network only shows 2nd order and higher channel segments. Red dots denote outlets for the COLAD, BCHAR, ELKAF, and BCHKL basins. The total 1970s study area streamflow for each season is given in the bottom left of each panel as area-average runoff (mm) and discharge equivalent (m³ s⁻¹).

For the UCRB as a whole (i.e., downstream of the Columbia–Kootenay confluence), runoff is projected to decrease on the order of 10 to 25%. The flush of increased summer runoff experienced in the middle reaches of the Columbia River (downstream of COLAD to about 100 km upstream of BCHAR) is attributed to increased glacier runoff from the northern region of the UCRB. However, these results should be interpreted with caution, because Schnorbus et al. (2012) point out that the simple glacier scheme in the VIC model likely underestimates the rate of glacier retreat and overestimates the volume of glacier runoff in the 2080s. For instance, Stahl and Moore (2006) suggest that glaciers in the UCRB are already in the declining phase of runoff generation in response to regional warming. Also, a recent study in the nearby leeward side of the Rockies in Alberta (Marshall et al., 2011) suggests that much of the ice mass has already melted over the past number of decades and that the decline will continue for perhaps two more decades. Hence, future runoff in the Columbia River may be overestimated and runoff changes may be less positive (or more negative) than shown. Nevertheless, results still serve to illustrate the sensitivity of runoff within the Columbia River to changes in glacier melt, precipitation, and temperature.

Projected changes in fall runoff also exhibit some spatial variation. Runoff in the upper reaches of the Kootenay River is expected to decline (e.g., between 0 and 10% at ELKAF), likely because of residual moisture deficits carried over from the summer months. Nevertheless, the middle reaches and tributaries of the Kootenay River show increased runoff, such that runoff near the outlet of the Kootenay River at BCHKL is expected to increase by 10 to 25%. Fall runoff changes along the Columbia River are generally positive, ranging from 0 to 10% at COLAD, increasing to 10 to 25% farther downstream at BCHAR. For the study region as a whole, total fall runoff is projected to increase on the order of 10 to 25%. Although early fall runoff is affected by glacier melt (more so in September, less so in October), the majority of runoff changes upstream of BCHAR are the result of increased rainfall throughout September, October, and November (see Section 4c); hence, fall runoff projections are expected to be more robust to uncertainty in glacier runoff than those for summer. On an annual basis, runoff in the 2080s is projected to increase from 0 to 10% in the Kootenay River and between 10 and 25% in the Columbia River.

Note that the majority of runoff in the UCRB (and, hence, the largest absolute runoff changes) occurs during the spring–summer half of the year, such that annual runoff changes generally reflect changes that occur during the spring and summer. Some of the increased annual runoff in the 2080s is attributed to projections of increased glacier runoff that may be overestimated by as much as 6% for the upper Columbia River at Mica Dam (see BCHMI in Fig. 1; Jost, Moore, Menounos, & Wheate, 2012). Errors would be less for the Kootenay River sub-basin where glacier area is considerably smaller.

c Temporal Response under the A2 Emissions Scenario

To demonstrate the influence of temperature, precipitation, evaporation, and snowmelt on runoff on a monthly time step for three future periods, combined results are shown for BCHAR for the 1970s, 2020s, 2050s, and 2080s (Figs 6 and 7). Individual plots show the historical and future distribution of each variable for the pooled seven-run A2 ensemble (i.e., 7 × 30 values for each month for each period). The historical period (1970s) distributions for each variable predominantly represent interannual variability because most aspects of GCM bias have been corrected. The distributions for the future periods represent potential changes in both interannual variability and variability in GCM response to the A2 radiative forcing. Any historical and future monthly distributions that exhibit statistically significant differences (two-sided Wilcoxon rank-sum test, α = 0.05, Helsel & Hirsch, 2002) are denoted in Figs 6 and 7.

Fig. 6 Temperature and precipitation for the historical period (1970s) and the 2020s, 2050s, and 2080s for seven GCMs for A2 emissions scenarios for the Columbia River at Keenlyside Dam. The solid lines indicate the median of the historical (black), 2020s (blue), 2050s (green), and 2080s (red) and the dashed lines indicate the 5th and 95th percentiles of the seven GCMs, 30-year ensemble for each month. Light grey shows the range of the 25th to 75th percentiles, and dark grey shows the 5th to 95th percentiles for the historical period. The blue diamond, green plus-sign, and red dot symbols denote significant differences (α = 0.05) between the 1970s and 2020s, 2050s and 2080s periods, respectively.

Fig. 7 As in Fig. 6 but for rainfall, snowfall, evaporation, snowmelt, and runoff.

Warmer temperatures are projected in all months and significant increases begin as early as the 2020s (Fig. 6). Nevertheless, median basin-wide temperature is projected to remain below freezing throughout winter, even in the 2080s. Precipitation is projected to increase between September and May for most months and time periods, except January and February when significant increases do not occur until the 2080s (Fig. 6). Precipitation decreases in July and August in all time periods. June is the only month when precipitation is not projected to change. As a result of increased precipitation and temperatures, rainfall is projected to increase throughout fall, winter, and spring (Fig. 7). Snowfall is also projected to increase in November and December because of increases in precipitation, whereas snowfall is projected to decrease in May through September for all time periods because of increasing temperatures and decreasing precipitation (Fig. 7).

Evaporation increases in all months for all time periods, except August when it starts to decrease in the 2080s (Fig. 7) because of limited soil moisture (Fig. 4). Snowmelt is projected to increase throughout the winter and spring, primarily at lower elevations where temperatures are warmer (not shown). Runoff increases in October through June, coincident with rainfall and snowmelt increases. During July through September, decreases in runoff generally coincide with decreases in snowmelt and rainfall and increases in evaporation.

d Temporal Response by Emissions Scenario

1 Annual streamflow

Projected annual discharge from 1950 through 2098 is shown for BCHAR and BCHKL in Fig. 8. The figure plots the mean and range of annual discharge for the full 20-member ensemble, as well as the individual ensemble means for the 7-member A1B and A2 and the 6-member B1 ensembles. Also shown are the 30-year scenario means for the 2020s, 2050s, and 2080s. The full range of annual discharge as shown reflects interannual variability, the variability in individual GCM response to prescribed emissions and, particularly for the 2080s, divergence in response among the three emissions scenarios. When integrated over a sufficiently long period, annual discharge changes reflect changes in the balance of precipitation and evapotranspiration over the basin. A trend of increasing mean annual discharge, which is apparent in both basins (Figs 8a and 8b), reflects the projected trend of increasing precipitation (Table 5). The trend from 2001 to 2098, based on the annual mean of the full 20-member ensemble, in BCHAR (2.1 m³ s⁻¹ y⁻¹) is more than twice that in BCHKL (0.9 m³ s⁻¹ y⁻¹) (see also Fig. 5). Because precipitation changes are uniformly distributed over the basin (not shown), the differences in annual discharge trends between the two sub-basins is attributed to differences in annual evapotranspiration. During the historical period, runoff is a higher proportion of precipitation (73%) in BCHAR than it is in BCHKL (61%) (Table 1), and this disparity continues into the future when evaporation in BCHAR tends to be more energy limited, whereas evaporation in BCHKL tends to be more water limited (see Section 4b). Consequently, a greater proportion of increased future precipitation is utilized for evaporation in the BCHKL than in the BCHAR. The differences in discharge changes between the A2, A1B, and B1 emissions scenarios, as reflected by the 30-year discharge means for each respective ensemble, become greater moving from the 2020s to the 2080s. In the 2050s, the greatest increase in mean annual discharge is projected under A1B in BCHAR and under B1 in BCHKL, with a 17% and 11% increase projected versus the 1970s, respectively. By the 2080s, projected increases are greatest under the A2 scenario at both locations, 27% in BCHAR and 17% in BCHKL. In the BCHAR, projected changes in annual discharge in the 2080s (as reflected by the 30-year scenario means) tend to lie outside the full ensemble range for the historical period, more so for the A2 scenario and less so for the B1 scenario. In the BCHKL, no scenarios tend to generate projections of annual discharge for the 2080s that lie outside the range of the historical ensemble.

Fig. 8 The spread of the mean annual discharge (grey) and 30-year means of mean annual discharge (horizontal black lines at 1961–90, 2011–40, 2041–70, 2071–98) are shown for 20 scenarios for a) the Columbia River at Keenlyside Dam (BCHAR) and b) the Kootenay River at Kootenay Canal (BCHKL). The means of the A1B, A2, and B1 scenarios (green, red, and light blue) are shown with solid lines, along with the 30-year means for the 2020s, 2050s, and 2080s.

Hamlet and Lettenmaier (1999) found that annual average streamflow increased for the 2020s, 2050s, and 2080s under the Hadley Centre Coupled Model verison 2 (HadCM2) GCM and decreased in the 2020s and 2050s under the Fourth Generation European Centre Hamburg Model (ECHAM4) for the Columbia River at the Dalles, Oregon. Annual runoff was projected to increase across Washington for the 2020s, 2050s, and 2080s for the A1B and B1 emissions scenarios (except A1B in the 2020s), although the magnitudes of the increases were much smaller than those found in this study because of smaller projected changes in annual precipitation (Elsner et al., 2010). Payne, Wood, Hamlet, Palmer, and Lettenmaier (2004) projected no change in annual runoff for the entire Columbia River basin upstream of the Dalles (because of a projection of stationary precipitation), but showed projections of increasing runoff during the mid- and end of the twenty-first century in the northern half of the UCRB (but decreasing runoff in the southern half of the UCRB).

2 monthly streamflow

Historical and future distributions of monthly discharge for four sub-basins of the UCRB is presented for the 1970s, 2020s, 2050s, and 2080s under the B1, A2, and A1B scenarios (Fig. 9). Distributions are constructed by pooling results for each scenario (i.e., samples sizes are 7 × 30 for the A1B and A2 scenarios and 6 × 30 for the B1 scenario for each month). Monthly discharge changes, based on differences between the median base and future consensus estimates for each scenario ensemble, are also summarized for each sub-basin for the 2080s in Table 6.

Table 6. Percentage change in median monthly streamflow for each sub-basin in the 2080s under each emissions scenario (A2, A1B, and B1). Historical values for 1961–90 (1970s) shown in italics.

		Oct	Nov	Dec	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep
BCHAR	1970s	553	361	194	121	90	104	275	982	2487	3274	2092	957
	A2	12	76	140	190	314	354	186	93	47	−2	−37	−24
	A1B	8	72	120	152	239	294	159	73	44	−6	−40	−29
	B1	7	45	74	89	157	205	123	61	29	3	−28	−21
COLAD	1970s	83	56	34	23	19	20	46	154	469	552	265	126
	A2	8	59	98	149	218	209	134	110	58	−32	−46	−36
	A1B	1	51	81	109	175	184	106	88	56	−31	−49	−36
	B1	−1	34	54	72	124	126	82	68	37	−13	−39	−25
BCHKL	1970s	301	297	222	173	148	187	487	1313	2592	1980	659	347
	A2	0	57	114	166	235	218	108	58	7	−54	−55	−48
	A1B	−1	56	102	139	181	184	89	47	10	−50	−55	−44
	B1	−4	43	62	76	132	132	69	40	9	−36	−44	−35
ELKAF	1970s	17	19	15	12	9	9	18	84	178	67	21	15
	A2	−14	31	84	148	250	289	218	80	−24	−63	−38	−41
	A1B	−11	26	77	124	208	243	177	66	−18	−57	−37	−36
	B1	−12	21	53	64	126	165	133	56	−14	−51	−30	−33

Fig. 9 Monthly discharge for four sub-basins. Solid lines indicate the median of the baseline (black), B1 (blue), A2 (green), and A1B (red) and dashed lines the 5th and 95th percentiles of the seven GCMs (six for B1) for a given emissions scenario, 30-year ensemble for each month. Light grey shows the range of the 25th to 75th percentiles and dark grey the 5th to 95th percentiles for the base period.

In all four sub-basins future monthly discharge tends to be larger than baseline discharge during winter and spring and smaller than baseline discharge from late summer through early fall under all three emissions scenarios, with the response becoming progressively more amplified through the 2020s, 2050s, and 2080s (Fig. 9). With few exceptions, differences between the baseline and future monthly distributions are significant in all sub-basins in all time periods under all scenarios for all months. Maximum relative increases are projected in February and March for all sub-basins, time periods, and emissions scenarios, reaching increases of up to 354% in BCHAR under the A2 emissions scenario in the 2080s (Table 6).

All four sub-basins respond qualitatively in the same way during the freshet months, wherein freshet runoff becomes progressively more advanced through the 2020s, 2050s, and 2080s (Fig. 9). However, some subtle differences between the sub-basins are revealed by their projected responses during this period. For instance, by the 2080s the peak discharge month occurs one month sooner in the COLAD, BCHAR and ELKAF sub-basins but not in the BCHKL sub-basin. Also, by the 2080s peak monthly discharge tends to increase in the BCHAR and COLAD sub-basins but not in the BCHKL sub-basin, and peak monthly discharge is projected to decrease in the ELKAF sub-basin. Nevertheless, even by the 2080s when mean annual temperatures are projected to warm by 4.4°C under the A2 scenario (Section 4a), all sub-basins maintain the characteristics of a nival regime. Subtle differences in response are also apparent for the summer and fall. Projected decreases in summer and early fall discharge tend to be higher in the southern portion of the UCRB (BCHKL and ELKAF, Fig. 9; see also Fig. 5). In the 2080s, projected summer decreases are estimated to be as high as 63% in July in ELKAF under the A2 emissions scenario (Table 6). In the ELKAF during the 2080s the period of decreased discharge is projected to occur from June through October, which is up to two months longer than that projected for the other sub-basins (Table 6).

The monthly streamflow response is not consistently different under one scenario over the others in any month in the 2020s or 2050s (Fig. 9). By the 2080s results for the different emissions scenarios begin to diverge noticeably (Fig. 9), and projected increases in winter and spring and projected decreases in summer and fall tend to be largest under the A2 and smallest under the B1 emissions scenarios (Fig. 9 and Table 6).

In snow-dominated regimes, such as the UCRB, projected changes in the seasonal distribution of streamflow predominantly reflect changes in the dynamics of snow accumulation and melt resulting from increasing temperatures. Given the strong agreement of increased temperature throughout the twenty-first century across a range of studies conducted throughout western North America (despite differences in methodology, including choice of GCMs and downscaling approach), projected changes in twenty-first century monthly runoff presented herein are broadly consistent with those projected for neighbouring watersheds, including the Fraser, Peace, and Athabasca rivers in British Columbia and Alberta (Kerkhoven & Gan, 2011; Schnorbus et al., 2012; Shrestha et al., 2012; Toth, Pietroniro, Conly, & Kouwen, 2006), and snow-dominated regimes in the Columbia basin in both Canada and the United States (Elsner et al., 2010; Hamlet & Lettenmaier, 1999; Wu et al., 2012).

3 Extremes

Low flows (7-day, 10-year return period) are projected to increase progressively for most GCM-driven ensemble members from the 1970s to the 2080s under the three scenarios in COLAD, BCHAR, ELKAF, and BCHKL, except for some ensemble members (Fig. 10). In these cold snow-dominated basins, low flows occur primarily in the winter during the 1970s when most precipitation remains in storage as snow and contributes little to streamflow over the winter (Fig. 7). With winter warming during the 2020s, 2050s, and 2080s, there is a tendency for more winter precipitation to fall as rain (Figs 3 and 7), producing runoff, which results in low flows of greater magnitude. HadGEM1 is one example where low flows are projected to increase in the 2020s and then decrease in the 2050s and 2080s for A1B and A2 in ELKAF and BCHKL, likely a result of dry antecedent conditions in fall resulting from high temperature and low precipitation in summer (see Figs 4-2 and 4-3 for the 2050s in Werner, 2011). In general, CGCM3.1 and MIROC3.2 driven ensemble members show the largest increase in low flows under all emissions scenarios likely because of larger precipitation increases projected in winter for these models (see Fig. 4-3 for the 2050s in Werner (2011)).

Fig. 10 Seven day, 10-year return period low flow for the Columbia River at Donald (COLAD), the Columbia River at Keenlyside (BCHAR), the Kootenay River at Kootenay Canal (BCHKL), and Elk River at Fernie (ELKAF) under the B1, A2, and A1B emissions scenario for the 1970s, 2020s, 2050s, and 2080s. Each box plot illustrates the median and inter-quartile range and the whiskers the upper and lower limits of the 30-year periods.

There is a wide range in projected future peak flows among individual GCMs; however, no one emissions scenario consistently shows a greater or lesser increase for any given sub-basin (Fig. 11). Future extreme peak-flow values (25-year return period) tend to be larger than the 1970s historical values for most sub-basins and scenarios, although the increase in peak discharge is not necessarily monotonic through the 2020s, 2050s, and 2080s (Fig. 11). An exception seems to be the ELKAF sub-basin where future 25-year peak flow magnitudes do not exhibit any clear tendency to be larger than historical values. The variability in projected changes in high streamflow extremes between climate projections reflects the subtle interplay between changes in temperature, seasonal snow storage, and precipitation. Peak flow events in the UCRB are most commonly attributed to snowmelt dynamics during the annual spring–summer freshet. In such cases, the magnitude of such peak-flow events would be governed by the intensity and duration of melt and the extent of the snow-covered area. Hence, the change in the magnitude of future events would depend on the interplay between opposing changes of increasing melt rates and decreasing extent and persistence of snow-covered area with increasing future temperatures. Rainfall can also play a role in the magnitude of peak-flow events. Rain-on-snow events (defined as events when precipitation occurs and snow depth decreases) that occur during the snowmelt season have been identified as the source of some of the largest local annual maximum flood events in the UCRB region (Schnorbus & Alila, 2004). Interior, high elevation sites in the western United States also experience rain-on-snow events, typically between October and May, although such events can occur as late as June as a result of a relatively long snow season, relatively large snowpack, and summer rainfall (McCabe, Hay, & Clark, 2007). Therefore, increased peak-flow magnitudes could also result from increased rainfall occurring along with snowmelt, either during the main snowmelt freshet period or during the fall and winter seasons but only in regions where the snowpack remains persistent in the future.

Fig. 11 As in Fig. 10 but for 25-year return period peak flow.

In the sub-basins of the UCRB during the 1970s, peak flows commonly occur in June or July in conjunction with the snowmelt freshet (Fig. 12). Peak-flow dates are projected to occur progressively earlier by a few weeks in all basins in the 2080s under the A2 and A1B emissions scenarios (Fig. 12). The shifts in the timing of the annual peak-flow events are attributed to earlier snowmelt with warmer temperatures and the possibility of increased rainfall in early summer (Fig. 7). Positive trends in the frequency of rain-on-snow events with warming temperatures have been found in the interior, high elevation areas of the western United States resulting from a persistent snowpack combined with precipitation falling as rain (McCabe et al., 2007). The coincident increase in rainfall and snowmelt found on a monthly time step in BCHAR (Figs. 6 and 7) would then also suggest an increased potential for fall or winter rain-on-snow conditions. Nevertheless, the evidence presented herein indicates that, by the end of the twenty-first century, annual maximum peak-flow events in three of the four examined sub-basins of the UCRB will continue to occur almost exclusively during the snowmelt freshet period. Although a few outlier simulations suggest that peak annual discharge may occasionally occur in the winter in BCHKL and BCHAR, this is not considered robust evidence for increased occurrence of fall/winter rain-on-snow events. Only ELKAF (and mainly for the A1B scenario) exhibits some evidence that future annual maximum peak-flow events may be generated more frequently by fall/winter rainfall or rain-on-snow events. Possible changes in the frequency of rain-on-snow events during the main snowmelt freshet period in all four sub-basins are ambiguous and cannot be diagnosed from Fig. 12.

Fig. 12 Annual maximum peak-flow magnitude versus day of peak flow for the Columbia River at Donald (COLAD), the Columbia River at Keenlyside (BCHAR), the Kootenay River at Kootenay Canal (BCHKL), and Elk River at Fernie (ELKAF) under the B1, A2, and A1B emissions scenario for the 1970s and 2080s.

5 Summary and conclusions

This study utilized a suite of seven GCMs driven by three emissions scenarios to project a wide range of plausible hydrologic futures for the 2020s, 2050s, and 2080s in the UCRB study area. Climate projections were statistically downscaled and used to force the VIC hydrologic model, which has been calibrated to 24 hydrometric sites throughout the UCRB. The hydrologic response over the UCRB was examined spatially, and streamflow projections were presented for four example sites within the study area.

Projected temperature, precipitation, SWE, evaporation, soil moisture, and runoff were evaluated over the UCRB under the A2 emissions scenarios out to the 2080s. Lower reaches were more susceptible to transitions from nival regimes to hybrid regimes, whereas some higher elevation areas saw increased snowpack moving from the 1970s to the 2080s. Evaporation became limited by lack of soil moisture in the southeastern part of the domain draining to the Kootenay River at the Kootenay Canal (BCHKL) project site. This suggests that in the future, water resources in the Kootenay River and its watershed, in the southeastern part of the UCRB, may become limited in the summer placing stress on forests and agriculture.

On an annual basis, mean discharge is projected to increase in the Columbia River at Keenlyside Dam (BCHAR) and BCHKL, with increases twice as great at BCHAR than at BCHKL. Response differs by emissions scenario in the 2020s and more so in the 2050s, but becomes quite dispersed in the 2080s, with the largest increases projected under the A2 emissions scenario. Seasonally, increased streamflow is projected for November through May consistently throughout the Columbia River. Decreases are projected in July through September, while signals can be mixed in June and October. Projections become greater in magnitude for each 30-year period, with larger increases in winter and spring and larger decreases in summer. In the 2050s, the largest increases are projected under the A1B emissions scenario and in the 2080s the largest winter/spring increases and summer decreases are projected under the A2 emissions scenario. These types of changes could challenge water management strategies by delivering water to reservoirs when water is already in excess and reducing water availability when resources are already limited and overallocated.

The 7-day, 10-year return period low flow is projected to increase for most sub-basins and GCM-driven ensemble members. The 7-day, 10-year return period low flows occur predominantly in winter, and their increases can be attributed to increases in the proportion of precipitation falling as rain and a larger area of the basin contributing to runoff as freezing levels rise. Future changes in the magnitude and timing of the 25-year peak-flow events results from the interplay between changes in temperature, seasonal snow storage, and precipitation. No one emissions scenario consistently shows a greater or lesser increase for any given sub-basin in the 25-year return period peak flows, although future extreme peak-flow values tend to be larger than the historical values for the 1970s for most sub-basins and scenarios. The changing temperature, snow/ice storage, and potential for rain-on-snow events also influence peak-flow date. The majority of peak-flow events will still tend to occur during the snowmelt freshet but will arrive, on average, approximately two weeks earlier, in June instead of July, under the A2 and A1B emissions scenarios by the 2080s. However, the occurrence of fall and winter rain-on-snow peak-flow events is projected to become more frequent in the Elk River at Fernie by the 2080s.

The hydrologic response to climate change differs by sub-regions in the UCRB. Effective management requires projected changes in streamflow not only for the UCRB as a whole but at project sites and other locations. Hydrologic changes have the largest amplitude under the A2 emissions scenario in the 2080s, which is an important period for the negotiation of international agreements and the management of heritage resources such as dams that help to produce a significant proportion of hydropower in BC.

Acknowledgements

We acknowledge the financial support of BC Hydro and the British Columbia Ministry of Environment. Katrina Bennett is gratefully acknowledged for her assistance in setting up the VIC model. Discussions with members of BC Hydro's Technical Advisory Committee and of the Climate Impacts Group, University of Washington improved the quality of this study. We appreciate the constructive feedback from two anonymous reviewers that helped to improve this manuscript.