Diagnosis of future changes in hydrology for a Canadian Rockies headwater basin

Climate change is anticipated to impact the hydrology of the Saskatchewan River, which originates in the Canadian Rockies mountain range. To better understand the climate change impacts in the mountain headwaters of this basin, a physically based hydrological model was developed for this basin using the Cold Regions Hydrological Modelling platform (CRHM) for Marmot Creek Research Basin (∼9.4 km²), located in the Front Ranges of the Canadian Rockies. Marmot Creek is composed of ecozones ranging from montane forests to alpine tundra and alpine exposed rock and includes both large and small clearcuts. The model included blowing and intercepted snow redistribution, sublimation, energy-balance snowmelt, slope and canopy effects on melt, Penman–Monteith evapotranspiration, infiltration to frozen and unfrozen soils, hillslope hydrology, streamflow routing, and groundwater components and was parameterised without calibration from streamflow. Near-surface outputs from the 4 km Weather Research and Forecasting (WRF) model were bias-corrected using the quantile delta mapping method with respect to meteorological data from five stations located from low-elevation montane forests to alpine ridgetops and running over October 2005–September 2013. The bias-corrected WRF outputs during a current period (2005–2013) and a future pseudo global warming period (PGW, 2091–2099) were used to drive model simulations to assess changes in Marmot Creek's hydrology. Under a “business-as-usual” forcing scenario, Representative Concentration Pathway 8.5 (RCP8.5) in PGW, the basin will warm up by 4.7 ^∘C and receive 16 % more precipitation, which will lead to a 40 mm decline in seasonal peak snowpack, 84 mm decrease in snowmelt volume, 0.2 mm d⁻¹ slower melt rate, and 49 d shorter snow-cover duration. The alpine snow season will be shortened by almost 1.5 months, but at some lower elevations there will be large decreases in peak snowpack (∼45 %) in addition to a shorter snow season. Declines in the peak snowpack will be much greater in clearcuts than under mature forest canopies. In alpine and treeline ecozones, blowing snow transport and sublimation will be suppressed by higher-threshold wind speeds for transport, in forest ecozones, sublimation losses from intercepted snow will decrease due to faster unloading and drip, and throughout the basin, evapotranspiration will increase due to a longer snow-free season and more rainfall. Runoff will begin earlier in all ecozones, but, as a result of variability in surface and subsurface hydrology, forested and alpine ecozones will generate the greatest runoff volumetric increases, ranging from 12 % to 25 %, whereas the treeline ecozone will have a small (2 %) decrease in runoff volume due to decreased melt volumes from smaller snowdrifts. The shift in timing in streamflow will be notable, with 236 % higher flows in spring months and 12 % lower flows in summer and 13 % higher flows in early fall. Overall, Marmot Creek's annual streamflow discharge will increase by 18 % with PGW, without a change in its streamflow generation efficiency, despite its basin shifting from primarily snowmelt runoff towards rainfall-dominated runoff generation.

1 Introduction

The eastern slopes of the Canadian Rockies form the headwaters of the Saskatchewan River basin (SRB) and are a vital source of water supply to downstream users in the Prairie Provinces of Canada. These high mountain headwaters occupy only 12.6 % of the SRB's total drainage area but provide 87 % of streamflow, which is used downstream for agriculture, communities, and industries (Redmond, 1964). Climate in this region has been experiencing changes since the last century (Whitfield, 2014; DeBeer et al., 2016). Western Canada experienced a significant warming of 0.5 to 1.5 ^∘C from 1900 to 1998, with the greatest increases found for winter daily minimum temperatures (Zhang et al., 2000). Warming of the eastern slopes of the Canadian Rockies is greater than the regional average, with mean temperatures increasing by 2.6 ^∘C and winter minimum temperatures increasing by 3.6 ^∘C at middle elevations in Marmot Creek Research Basin (MCRB) since the early 1960s (Harder et al., 2015). With the warming air temperatures, the rainfall ratio (ratio of rainfall to total precipitation) is increasing as the fraction of precipitation as snowfall declines (Lapp et al., 2005; Shook and Pomeroy, 2012). In mountains of western North America this leads to decreases in the seasonal snowpack (Mote et al., 2005; Brown and Robinson, 2011) and consequently earlier spring runoff (Stewart et al., 2004). In contrast to ubiquitous warming temperatures, trends in precipitation volume are mixed for the Canadian Rockies, with some studies showing increasing trends of about 14 % over the period 1948–2012 (Vincent et al., 2015) and other studies finding neither trends nor change (Valeo et al., 2007; Harder et al., 2015). With further climate change expected in the future (IPCC, 2013), understanding the impacts of projected climate change on the hydrological cycle in mountain headwater basins is important for evidence-based water management of the SRB in the future.

Winter snow accumulation provides the greatest source of streamflow runoff in many mountain regions of the world (Grant and Kahan, 1974; Serreze et al., 1999), as snowmelt is the most important annual hydrological event (Gray and Male, 1981). Melt from the seasonal snowpack is the main contributor to streamflow in the eastern slopes of the Canadian Rockies (Kienzle et al., 2012; Pomeroy et al., 2012; Fang et al., 2013). Streamflow generation in mountain regions is highly variable and is controlled by many biophysical and hydrometeorological factors (Hunsaker et al., 2012; Zhang and Wei, 2014). Elevation affects both air temperature and precipitation, two of the most important drivers of snowpack variability in mountains (Lundquist and Cayan, 2007; Marks et al., 2013). Topographic features, such as slope/aspect and ecosystem features such as forest structure, are other important factors contributing to the heterogeneity of radiation, atmospheric energy, and wind flow in mountain environments (Föhn and Meister, 1983; Bernhardt et al., 2009; Marsh et al., 2012; Musselman and Pomeroy, 2017; MacDonald et al., 2018) and result in the high spatial variability of snow accumulation, melt patterns, evapotranspiration, and runoff in complex mountain terrain (Pomeroy et al., 2005; MacDonald et al., 2010; Ellis et al., 2013; Revuelto et al., 2014; Knowles et al., 2015; DeBeer and Pomeroy, 2017).

Many studies have examined the impacts of climate change on snow accumulation, redistribution, snowmelt, evapotranspiration, soil moisture storage, and streamflow in mountain drainage basins through the simulations of hydrological models driven by future climate scenarios generated by downscaling of climate model outputs or perturbations of current meteorological observations (Kienzle et al., 2012; López-Moreno et al., 2014; Rasouli et al., 2015; Jepsen et al., 2016; Weber et al., 2016; Meißl et al., 2017). Physically based hydrological models are effective ways to analyse the hydrological response to climate change, as they can capture the complex hydrological processes governing streamflow generation and can be extrapolated beyond the hydrometeorological conditions under which they were developed. Empirical snowmelt modelling methods that use temperature-index techniques are inappropriate in cold mountain regions (Swanson, 1998) and generally do not perform well because of their lack of physical basis, need for calibration from sparse snowmelt observations, and neglect of sublimation contributions to ablation (Walter et al., 2005; Pomeroy et al., 2005, 2013). The Cold Regions Hydrological Modelling platform (CRHM; Pomeroy et al., 2007, 2016) offers a full suite of streamflow generation processes that commonly operate in the Canadian Rockies, such as wind redistribution of snow, snow avalanching, canopy snow and rain interception, sublimation, drip and unloading from forest canopies, infiltration to frozen and unfrozen soils, overland and detention flow, hillslope sub-surface water redistribution, and evapotranspiration. CRHM is an object-oriented, modular, and flexible platform for assembling physically based hydrological models. With CRHM, the user constructs a purpose-built model from a selection of possible basin spatial configurations, spatial resolutions, and physical process modules of varying degrees of physical complexity. Basin discretisation is performed via hydrological response units (HRUs) whose number and nature are selected based on the variability of basin attributes and the level of physical complexity chosen for the model. The user, in light of hydrological understanding, parameter availability, basin complexity, meteorological data availability, and the objective flux or state for prediction, selects physical complexity. A full description of CRHM is provided by Pomeroy et al. (2007). Physically based algorithms in CRHM have been developed from field process studies (Pomeroy et al., 2009; DeBeer and Pomeroy, 2010; Ellis et al., 2010; MacDonald et al., 2010; Harder and Pomeroy, 2013) and have been extensively evaluated using model falsification and multivariate evaluation in mountain headwater basins where models created using CRHM can be run successfully without calibration from streamflow (Fang et al., 2013; Pomeroy et al., 2013, 2016; Rasouli et al., 2015; Fang and Pomeroy, 2016).

Figure 1 Location and contour map of the Marmot Creek Research Basin (MCRB), showing hydrometeorological stations, a hydrometric station, a WRF grid centroid, and ecozones of the MCRB: alpine, treeline, upper forest, forest clearing blocks, forest circular clearings, and lower forest. Note that the WRF grid centroid represents the centre of the 4 km WRF grid, and the sizes and areas of circular clearings in Twin Creek are meant to represent a honeycomb pattern and are not to scale.

A recent application of the Weather Research and Forecasting (WRF) model provides 4 km simulation outputs for both current climate and a future climate scenario using dynamical downscaling from reanalysis data for large portions of western North America with perturbations from an ensemble of regional climate model (RCM) projections called “pseudo global warming” (PGW) as discussed by Liu et al. (2017) and Li et al. (2019). This remarkably high-resolution WRF application permits the representation of convective precipitation processes, resolves mountain topography, and so can capture variations in near-surface meteorology due to mesoscale orography such as found in mountains. The objectives of this paper are to use CRHM driven by WRF to (1) evaluate the ability to simulate snowpack and streamflow regimes in a Canadian Rockies headwater basin without calibration; and (2) diagnose the detailed changes in hydrology due to impending climate change for this headwater basin. By relying on physically based, uncalibrated simulations and dynamical downscaling, it is hoped that the approach introduces a highly robust method for evaluating the impacts of climate change on mountain hydrology.

Figure 2 HRU map and ecozones of Marmot Creek Research Basin. Note that the sizes and areas of circular clearings in Twin Creek are not to scale.

Table 1 Area and mean elevation, aspect, and slope for ecozones at the Marmot Creek Research Basin. Note that the aspect is in degree clockwise from north.

2 Methods

2.1 Study site

The study was conducted in the Marmot Creek Research Basin (MCRB) (50.95^∘ N, 115.15^∘ W) in the Kananaskis Valley, Alberta, Canada, located on the eastern slopes of the Canadian Rockies (Fig. 1). MCRB is a small headwater basin (9.4 km²) of the Kananaskis River which contributes to the Bow River and ultimately the Saskatchewan River basin. MCRB is composed of three upper sub-basins: Cabin Creek (2.35 km²), Middle Creek (2.94 km²), and Twin Creek (2.79 km²), which converge into the Confluence Sub-basin above the main stream gauge (1.32 km²). MCRB elevation ranges from 1590 m a.s.l. (above sea level) at the main gauged outlet to 2829 m at the summit of Mount Allan. Much of MCRB is covered by needleleaf forest; Engelmann spruce (Picea engelmanni) and subalpine fir (Abies lasiocarpa) are dominant in the upper part of the basin (1710 to 2277 m). The lower-elevation (1590 to 2015 m) forests are mainly Engelmann spruce and lodgepole pine (Pinus contorta var. Latifolia) with aspen woodland near the basin outlet (Kirby and Ogilvy, 1969). Alpine larch (Larix lyallii) and short tundra shrubs, tree islands, dwarf trees, and krummholz are found around the treeline at approximately 2016 to 2379 m. Exposed rock surfaces and talus predominate in the high alpine part of the basin (1956 to 2829 m). Forest management experiments conducted in the 1970s and 1980s left large clearcuts (1838 to 2062 m) in the Cabin Creek sub-basin and numerous small circular clearings (1762 to 2209 m) in the Twin Creek sub-basin (Golding and Swanson, 1986). Physiographic descriptions of these ecozones are shown in Table 1. Surficial soils are poorly developed mountain soils consisting principally of glaciofluvial, surficial till, and postglacial colluvium deposits (Beke, 1969). Relatively impermeable bedrock is found near the surface at the higher elevations and headwater areas, whilst the rest of the basin is covered by a deep layer of coarse and permeable soil allowing for rapid rainfall infiltration to subsurface layers overlying relatively impermeable shale (Jeffrey, 1965). The ecozones defined by vegetation, soils, bedrock, slope, aspect, and hydrography form the land-cover types that contribute to defining the HRUs by which CRHM discretises the basin and are shown in Fig. 2. Fang et al. (2013) describe the generation of these HRUs in MCRB in detail.

Continental air masses control the weather in the region, which has long and cold winters and cool and wet springs with a late spring/early summer precipitation maximum that can fall as rainfall or snowfall. Westerly warm and dry Chinook (föhn) winds lead to brief periods with air temperatures well above 0 ^∘C during the winter months. Annual precipitation ranges from 600 mm at lower elevations to more than 1100 mm at the higher elevations, of which approximately 70 % to 75 % occurs as snowfall, the percentage increasing with elevation (Storr, 1967). Mean monthly air temperatures range from 14 ^∘C observed at 1850 m in July to −10 ^∘C observed at 2450 m in January.

2.2 WRF model

2.2.1 Model overview

The WRF model Version 3.4.1 was used in this modelling experiment. This version of WRF simulates convective weather systems at sub-synoptic scales and deals with mesoscale orography at a 4 km horizontal grid spacing for large portions of North America. Two 13-year experiments were conducted, consisting of a control (CTRL) simulation and a PGW simulation. The CTRL simulation is retrospective for the 2000–2013 period, with initial and boundary conditions from 6 h 0.703^∘ ERA-Interim reanalysis data (Dee et al., 2011). The PGW simulation is a 13-year (i.e. 2000–2013) simulation forced with the 6-hour ERA-Interim reanalysis data plus a climate perturbation. The climate perturbation was derived from 19-model ensemble mean change from the fifth phase of the Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012) under a “business-as-usual” forcing scenario: Representative Concentration Pathway 8.5 (RCP8.5; van Vuuren et al., 2011). The PGW simulation is equivalent to a future climate scenario for the 2086–2099 period. The perturbation for the WRF PGW simulation using the 19-model ensemble mean can reduce uncertainties in climate projections from global climate models (GCMs) due to inter-GCM variability (Li et al., 2019). A detailed description of the WRF model setup is provided by Li et al. (2019). For MCRB, the 4 km hourly WRF outputs from both CTRL and PGW simulations were extracted from the WRF grid shown in Fig. 1. The extracted variables include near-surface air temperature, vapour pressure, wind speed, precipitation, and shortwave irradiance; relative humidity is required by CRHM and was estimated with conversion equations using air temperature and vapour pressure (Tetens, 1930).

2.2.2 Bias correction

Although the 4 km WRF model allows direct use of microphysics and resolves mesoscale convection to provide a considerable level of spatial detail, it still produced biases in the near-surface meteorology for MCRB. The biases were caused by the complex mountain terrain in MCRB and other factors. Many of the topographic features such as alpine ridges, wind-exposed and wind-sheltered slopes, and valley bottoms influence the distribution of near-surface meteorological conditions at scales smaller than 4 km (Vionnet et al., 2015). To produce a more reliable driving meteorology dataset, the 4 km WRF outputs were bias-corrected using the quantile delta mapping (QDM) algorithm against data from five meteorological stations in MCRB (Cannon et al., 2015). The QDM algorithm corrects the systematic bias in quantiles of WRF outputs with respect to the observations and preserves model-projected relative changes in quantiles. First, the QDM algorithm was used to extract the climate change signal from projected future quantiles and then was used to detrend the series before reintroducing the trends in projected future quantiles. A transfer function that transformed the cumulative distributions of the model outputs to match the cumulative distributions of observed data was employed in the QDM algorithm to correct both historical and projected model outputs. More details of QDM are provided by Cannon et al. (2015). Observations of air temperature, relative humidity, wind speed, precipitation, and shortwave irradiance from the Centennial Ridge, Fisera Ridge, Vista View, Upper Clearing, and Upper Forest hydrometeorological stations were used to correct WRF model outputs using the QDM algorithm. These stations are shown in Fig. 1 and are described in several publications (DeBeer and Pomeroy, 2010; Ellis et al., 2010; MacDonald et al., 2010; Pomeroy et al., 2012). This bias correction downscaled the WRF grid to the locations of these stations to create a sub-grid WRF “virtual station” surface meteorology that was used to force CRHM simulations. The bias correction was performed for 8 water years (i.e. 1 October to 30 September) from 2005 to 2013 that had overlap of WRF CTRL outputs and MCRB observations. The WRF PGW outputs were similarly bias-corrected by preserving the model-projected relative changes in quantiles, resulting in 8 water years of corrected WRF PGW outputs from 1 October 2091 to 30 September 2099. Statistical indexes used to assess WRF CTRL outputs included the root mean square difference (RMSD) calculated as in Fang et al. (2013) and the mean absolute difference (MAD) computed as follows:

$$\begin{array}{}\text{(1)}& \mathrm{MAD}={\displaystyle \frac{\mathrm{1}}{n}}\sum \mathrm{|}{x}_{\mathrm{s}}-x_{\mathrm{o}}\mathrm{|},\end{array}$$

where n is the number of samples, and x_s and x_o are the observed and modelled values, respectively.

2.3 Hydrological model and simulations

CRHM (Pomeroy et al., 2007) was used to create a hydrological model for MCRB. For MCRB, a set of physically based modules was assembled to simulate the dominant hydrological processes by Pomeroy et al. (2012) and Fang et al. (2013), including wind redistribution of alpine snow, snow interception, sublimation, drip and unloading from forest canopies, sub-canopy radiation energetics, slope/aspect effects on radiation and wind flow, infiltration to frozen and unfrozen soils, overland flow, hillslope sub-surface water flow and storage, and evapotranspiration from forests, clearings, and alpine tundra. Recent updates were made to the evaporation and hillslope modules for better representation of runoff, including detention flow, on hillslopes and evapotranspiration from vegetation with seasonal variations in leaf area index and height, and the updated model was evaluated in the June 2013 flood when approximately 250 mm precipitation fell at MCRB during 17–24 June 2013 (Fang and Pomeroy, 2016; Pomeroy et al., 2016).

Hydrological model simulations were conducted driven by the bias-corrected WRF near-surface meteorological variables: air temperature, relative humidity, wind speed, precipitation, and shortwave irradiance for the CTRL and PGW periods, respectively. Simulations in both periods covered 8 water years: CTRL for the current period (i.e. 1 October 2005 to 30 September 2013) and PGW for the future period (i.e. 1 October 2091 to 30 September 2099). Model simulations of snow accumulation, spring snowmelt, and streamflow in the CTRL period were evaluated against observations of snow accumulation, snowmelt, and streamflow. Statistical indexes used to evaluate model simulations in the CTRL period were the Nash–Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970), RMSD, normalised RMSD (NRMSD), and model bias (MB) calculated as per Fang et al. (2013).

3 Results

3.1 WRF CTRL outputs

Near-surface hourly air temperature, relative humidity, wind speed, precipitation, and shortwave irradiance from observations, uncorrected WRF CTRL outputs, and bias-corrected WRF CTRL outputs were compared for the Centennial Ridge, Fisera Ridge, Vista View, Upper Clearing, and Upper Forest stations in MCRB, and the comparisons are shown in the quantile–quantile (Q–Q) plots (Fig. 3). The points in the Q–Q plots of observations and bias-corrected WRF outputs shown in Fig. 3b, d, f, h, and j are linearly distributed on the 1:1 line, while the points in the Q–Q plots of observations and uncorrected WRF outputs shown in Fig. 3a, c, e, g, and i do not appear to be linearly distributed. This suggests that the near-surface meteorological variables from observations and bias-corrected WRF CTRL outputs form the same distribution. The Q–Q plots for other stations show the same results and are provided as the Supplement. Table 2 shows MAD and RMSD indexes for the WRF CTRL outputs compared to observations. The uncorrected WRF CTRL relative humidity was converted using air temperature, specific humidity, and specific pressure of air outputs and had values higher than 100 % (Fig. 3c), indicating errors in these uncorrected WRF CTRL outputs. The MAD values for the bias-corrected WRF outputs were zero and were smaller than those for the uncorrected WRF outputs, suggesting there is no difference in the statistical distributions of observations and bias-corrected WRF CTRL outputs. The RMSD values for the bias-corrected WRF outputs were lower than those for the uncorrected WRF outputs, with two exceptions: the RMSD values were 4.88 m s⁻¹ and 0.6 mm for bias-corrected WRF wind speeds for the Centennial Ridge station and bias-corrected WRF precipitation for the Fisera Ridge station, respectively, and were slightly higher than the RMSD of 4.61 m s⁻¹ and 0.56 mm for the original WRF wind speed and precipitation. Despite that, the bias correction generally improved WRF outputs and reduced the difference from the observations.

Figure 3 Quantile–quantile plots of observations and WRF CTRL outputs for the Fisera Ridge station in MCRB: (a) WRF CTRL and observed air temperature, (b) corrected WRF CTRL and observed air temperature, (c) WRF CTRL and observed relative humidity, (d) corrected WRF CTRL and observed relative humidity, (e) WRF CTRL and observed wind speed, (f) corrected WRF CTRL and observed wind speed, (g) WRF CTRL and observed incoming solar radiation, (h) corrected WRF CTRL and observed incoming solar radiation, (i) WRF CTRL and observed precipitation, and (j) corrected WRF CTRL and observed precipitation.

Table 2 Comparison of observations and WRF outputs for MCRB stations, with mean absolute difference (MAD) and root mean square difference (RMSD). Values are for bias-corrected WRF outputs; values inside parentheses are for WRF outputs without bias correction.

3.2 Hydrological model evaluations

Without changing any model parameters from previously published runs made using locally observed meteorological forcing, CRHM simulations of the snow water equivalent accumulation (SWE) using bias-corrected WRF CTRL near-surface outputs were compared to observed SWE for the sheltered, mid-elevation Upper Forest and Upper Clearing sites (Fig. 4a and b) and for the wind-blown, high-elevation Fisera Ridge site (Fig. 4c–f) for 2007–2013. The results demonstrate that the model forced with the bias-corrected WRF outputs was able to simulate SWE for both forest and alpine environments, with an exception in the Upper Forest site during the season of 2007–2008. In addition, Table 3 shows that there were large differences between the simulations and the observations for both forest and alpine sites during the season of 2011–2012, with RMSD ranging from 52.3 to 297.5 mm for spruce forest and lower alpine south-facing slope and MB ranging from −0.65 to −0.35 for spruce forest and larch forest treeline, respectively. For all six seasons, model simulations captured the general seasonal patterns of snow accumulation and ablation for these forest and alpine sites, with MB values of all seasons ranging from −0.43 for the forest clearing to 0.17 for the alpine ridge top (Table 3). All seasons' RMSD ranged from 46.5 mm for the mature spruce forest to 260.1 mm for the larch forest treeline, while the NRMSD ranged from 0.39 for upper alpine south slope to 0.84 for the mature spruce forest.

Figure 4 Comparison of the observed and simulated snow accumulation (SWE) for 2007–2013 at the sheltered, mid-elevation Upper Forest/Clearing and wind-blown, high-elevation Fisera Ridge sites in MCRB. (a) Mature spruce forest, (b) forest clearing, (c) ridge top, (d) upper alpine south-facing slope, (e) lower upper alpine south-facing slope, and (f) larch forest treeline.

Further model evaluation was conducted using the CRHM-simulated streamflow driven by bias-corrected WRF CTRL outputs and the observed outlet streamflow discharge measured by the Water Survey of Canada (WSC) gauge (05BF016) for 2005–2013 (Fig. 5). Again, no parameters were changed for these model runs from previously published CRHM simulations of Marmot Creek that were driven by observed meteorology. The WSC gauged streamflow from 1 May to 30 September during 2006–2012 and for part of 2013 before the flood. Streamflow during the 2013 flood was estimated by the University of Saskatchewan Centre for Hydrology (CH) using the best available information but with great uncertainty (Harder et al., 2015). After the flood, a streamflow gauge established near the same outlet by CH continued measurements for the remainder of 2013. The seasonal NSE values ranged from −0.33 in 2009 to 0.72 in 2012, and the overall eight-season NSE was 0.4 (Table 4), suggesting the model had some predictability for the temporal evolution of daily basin discharge in these eight seasons. The eight-season RMSD, NRMSD, and MB listed in Table 4 were 0.212 m³ s⁻¹, 0.79, and 0.001 for the predicted daily basin discharge, respectively, indicating relatively small differences between the simulated and observed Marmot Creek streamflow.

Table 3 Evaluation of simulated snow accumulation using the RMSD (mm SWE), normalised RMSD (NRMSD), and model bias (MB) at the Upper Forest/Clearing and Fisera Ridge sites, Marmot Creek Research Basin.

Table 4 Evaluation of simulated daily mean streamflow discharge for Marmot Creek using Nash–Sutcliffe efficiency (NSE), RMSD (m³ s⁻¹), NRMSD, and MB.

Figure 5 Comparison of the observed and simulated daily streamflow at the basin outlet over 2005–2013 for Marmot Creek.

3.3 Changes in WRF meteorology due to climate change

The WRF near-surface meteorological variables that had been bias-corrected with respect to MCRB stations were compared between CTRL and PGW. Figure 6 shows the comparisons of mean air temperature, relative humidity, wind speed, and shortwave irradiance and cumulative water year (WY) precipitation at the Centennial Ridge, Fisera Ridge, Vista View, Upper Clearing, and Upper Forest stations. There were consistent increases in mean air temperatures in PGW for all MCRB stations compared to those in CTRL. Mean temperatures in CTRL ranged from −1.8 ^∘C at Centennial Ridge to 1.4 ^∘C at Vista View; however, those in PGW ranged from 2.9 ^∘C at Centennial Ridge to 6.1 ^∘C at Vista View and so were about 4.7 ^∘C warmer for all stations (Fig. 6a). The mean relative humidity decreased slightly with PGW compared to that in CTRL, declining by 1.8 % for all stations (Fig. 6b), whilst the mean wind speed remained unchanged between CTRL and PGW (Fig. 6c). There were slight increases in the mean shortwave irradiance in PGW compared to CTRL, with mean shortwave irradiance increasing by 2.8 W m⁻² at both Centennial Ridge (i.e. from 141.4 to 144.2 W m⁻²) and Fisera Ridge (i.e. from 152.4 to 155.2 W m⁻²) and increasing by 2.1 W m⁻² at Upper Clearing (i.e. from 140.7 to 142.8 W m⁻²) (Fig. 6d). There was much more precipitation with PGW, and annual WY precipitation was 1287 and 882 mm with PGW at Fisera Ridge and Upper Clearing, respectively, about 147 and 150 mm more or 13 % and 20 % increases compared to the CTRL period at Fisera Ridge and Upper Clearing (Fig. 6e).

Figure 6 Boxplots of the bias-corrected WRF CTRL and PGW near-surface meteorology for MCRB station sites. (a) Air temperature, (b) relative humidity, (c) wind speed, (d) solar irradiance, and (e) precipitation. Note that the accumulation over the water year is used for precipitation, and the average value over the water year is presented for the other variables. The vertical line within the box is the median value of the 8 water year data, the box is the interquartile range (Q₁: 25 % to Q₃: 75 %) of the 8 water year data, the whiskers are the minimum and maximum, and the circle is the mean value of the 8 water year data.

Figure 7 Boxplots of the simulated annual water balance variables for WRF CTRL and PGW for the alpine (AE), treeline (TE), upper forest (UF), forest clearing blocks (FCB), forest circular clearing north-facing (FCCNF), forest circular clearing south-facing (FCCSF), and lower forest (LF) ecozones. (a) Rainfall, (b) snowfall, (c) actual evaporation, (d) sublimation, (e) blowing snow transport, (f) surface runoff, (g) subsurface flow, (h) groundwater flow, and (i) storage change. The vertical line within the box is the median value of the 8 water year data, the box is the interquartile range (Q₁: 25 % to Q₃: 75 %) of the 8 water year data, the whiskers are the minimum and maximum, and the circle is the mean value of the 8 water year data.

Table 5 Simulated basin-scale mean annual air temperature, relative humidity, wind speed, mean annual water balance fluxes, mean annual radiation fluxes to snow cover, mean seasonal cumulative snowmelt, peak snow accumulation, snow-cover duration, and snowmelt rate for WRF CTRL and PGW.

Figure 8 Simulated mean annual snow accumulation (SWE) for WRF CTRL and PGW. (a) Alpine, (b) treeline, (c) upper forest, (d) forest clearing blocks, (e) forest circular clearing north-facing, (f) forest circular clearing south-facing, and (g) lower forest ecozones, and (h) MCRB. The line represents the mean and the shadow represents the standard deviation of the 8 water year SWE.

3.4 Changes in water balance variables

The simulated annual water balance variables for all ecozones were compared between the CTRL and PGW periods. Compared to CTRL, rainfall increased and snowfall decreased with PGW for all ecozones (Fig. 7a and b). The increase in annual rainfall ranged from 201 mm in the lower forest ecozone to 328 mm in the alpine ecozone, whilst the decrease in annual snowfall ranged from 63 mm in the lower forest ecozone to 168 mm in the alpine ecozone. Averaged over the whole basin, annual rainfall rose by 268 mm and annual snowfall declined by 112 mm, and there was a 156 mm or 16 % increase in total precipitation with PGW (Table 5). Actual evapotranspiration is the sum of evaporation from soil, intercepted rain, and open water, and transpiration from plants; this increased with PGW for all ecozones because of more rainfall (Fig. 7c). The increase in annual actual evapotranspiration ranged from 59 mm at the alpine ecozone to 179 mm at the upper forest ecozone, with a 124 mm increase for the whole basin shown in Table 5. Sublimation is the total flux of snow sublimated from the surface snowpack and during blowing snow and forest canopy interception processes and declined with PGW for all ecozones as a result of decreased snowfall and the impact of warmer air temperatures in limiting blowing snow occurrence and increasing unloading of intercepted snow from forest canopies; the decrease in annual sublimation ranged from 6 mm at the forest circular clearing north-facing ecozone to 58 mm at the alpine ecozone (Fig. 7d), with a reduction of 40 mm for the entire basin (Table 5). For the sheltered and sparsely vegetated forest clearing ecozones, neither blowing snow nor forest interception processes occur, so there was only a change in surface snowpack sublimation for these ecozones. For the alpine and treeline ecozones, blowing snow is a very important process in controlling seasonal snow accumulation, and the alpine ecozone is the source area for blowing snow (i.e. negative value for blowing snow transport), whilst the treeline ecozone accumulates blowing snow in deep snowdrifts (i.e. positive value for blowing snow transport). Figure 7e shows that blowing snow transport was suppressed with PGW, resulting in a smaller annual blowing snow transport loss from the alpine ecozone of 24 mm and a smaller annual blowing snow transport gain to the treeline ecozone of 53 mm. Blowing snow does not occur below the treeline in MCRB, and so no changes in blowing snow transport occurred in subalpine ecozones. Annual surface runoff increased by 250 mm with PGW for the alpine ecozone (Fig. 7f) and by 87 mm for the entire basin (Table 5). In other, subalpine, ecozones, flow predominantly occurred in the subsurface and did not show change in the surface runoff with PGW. Annual subsurface flow from alpine and treeline ecozones decreased by 61 and 24 mm, respectively, with PGW. Although there were increases in subsurface flow with PGW for other ecozones, ranging from 8 mm from the lower forest ecozone to 56 mm from the forest circular clearing north-facing ecozone (Fig. 7g), on average over the basin, there was a 12 mm reduction in subsurface flow with PGW (Table 5). In contrast, groundwater flow stayed relatively constant with PGW for all ecozones (Fig. 7h). There were increases in total subsurface storage in soil and groundwater with PGW for some ecozones, ranging from 9 mm in the upper forest ecozone to 27 mm in the forest circular clearing south-facing ecozone, whilst the subsurface storage dropped, ranging from 2 mm at the lower forest ecozone to 40 mm at the forest clearing blocks ecozone with PGW (Fig. 7i). Subsurface storage for the entire basin declined by 12 mm, from 416 mm under CTRL (i.e. 45 % saturation) to 404 mm with PGW (i.e. 43 % saturation), as shown in Table 5.

Figure 9 Boxplots of the simulated water year's (a) cumulative snowmelt, (b) peak snow accumulation (SWE), (c) snow-cover duration, (d) incoming solar radiation to snow, (e) longwave irradiance to snow, and (f) net radiation to snow for WRF CTRL and PGW for the alpine (AE), treeline (TE), upper forest (UF), forest clearing blocks (FCB), forest circular clearing north-facing (FCCNF), forest circular clearing south-facing (FCCSF), and lower forest (LF) ecozones. The vertical line within the box is the median value of the 8 water year data, the box is the interquartile range (Q₁: 25 % to Q₃: 75 %) of the 8 water year data, the whiskers are the minimum and maximum, and the circle is the mean value of the 8 water year data.

3.5 Changes in snow regime

The simulated snowpack accumulation (SWE) for all ecozones and the entire basin was compared between CTRL and PGW periods. Figure 8 illustrates the annual time series of SWE for CTRL and PGW simulations and demonstrates the impacts of PGW on seasonal SWE for different ecozones. For all ecozones, the peak SWE occurred earlier and was much lower with PGW, with declines in peak SWE ranging from 7 mm at the upper forest ecozone to 166 mm at the treeline ecozone. In the alpine ecozone, the cold mid-winter snowpack was not impacted by PGW and had similar to higher levels by early April compared to the CTRL, but shortly after alpine SWE ablated rapidly and disappeared 26 d earlier with PGW. For the subalpine ecozones, the seasonal snowpack declined substantially and throughout the season with PGW. The date of seasonal snowpack depletion advanced from early August to late June with PGW at the treeline ecozone and from mid-June to late May in the other ecozones. The mean snowmelt rate was estimated by dividing mean annual peak SWE by the number of days from peak SWE to snowpack disappearance, and with PGW it was lower for the treeline and forest clearing ecozones, ranging from 0.9 mm d⁻¹ in the forest clearing blocks ecozone (i.e. from 1.4 mm d⁻¹ in CTRL to 0.5 mm d⁻¹ in PGW) to 1.6 mm d⁻¹ in the forest circular clearing north-facing ecozone (i.e. from 2.9 mm d⁻¹ in CTRL to 1.3 mm d⁻¹ in PGW), whilst the melt rate was slightly higher for the alpine ecozone (i.e. from 1.9 mm d⁻¹ in CTRL to 2.0 mm d⁻¹ in PGW) and remained unchanged for the forest ecozones (i.e. 0.6 mm d⁻¹ at upper forest and 0.5 mm d⁻¹ at lower forest in both CTRL and PGW). For the entire basin, there was a very small decline in the melt rate from 1.3 mm d⁻¹ in CTRL to 1.1 mm d⁻¹ in PGW (Table 5).

Changes in the seasonal snowmelt, peak SWE, snow-cover duration, and radiation fluxes to snow cover were also compared between CTRL and PGW. Figure 9a shows that cumulative snowmelt volume decreased with PGW for all ecozones, and seasonal snowmelt at the treeline ecozone declined the most (215 mm), with decreases elsewhere ranging from 32 mm at the upper forest to 113 mm at the alpine ecozone. Peak SWE declined for all ecozones with PGW, with the largest decrease in the treeline ecozone (149 mm) and the lowest at the upper and lower forest ecozones (11 mm), as shown in Fig. 9b. The duration of seasonal snow cover declined for all ecozones with PGW, ranging from 31 d at the forest circular clearing north-facing ecozone to 49 d at the treeline ecozone (Fig. 9c). Table 5 shows that the basin-wide snowmelt volume, peak SWE, and duration of seasonal snow cover decreased by 84 mm, 40 mm, and 49 d, respectively. The seasonal net radiation to snow cover increased slightly with PGW for all ecozones, ranging from an additional 2 W m⁻² at the lower forest to 4 W m⁻² at the other ecozones (Fig. 9f). The increases in the net radiation to snow cover were due to higher longwave irradiance with PGW for all ecozones, ranging from an additional 10 W m⁻² to the lower forest ecozone to 17 W m⁻² to the treeline and forest clearing ecozones (Fig. 9e), whilst the seasonal solar irradiance to snow cover declined for all ecozones with PGW, ranging from 2 W m⁻² less at the upper forest ecozone to 17 W m⁻² less at the forest clearing blocks ecozone (Fig. 9d). Basin-wide, seasonal solar irradiance to snow cover decreased by 11 W m⁻² and longwave irradiance to snow cover increased by 14 W m⁻², resulting in an increase of 4 W m⁻² in annual net radiation to snow cover with PGW (Table 5).

Figure 10 Simulated mean annual (a) Marmot Creek daily streamflow discharge and (b) cumulative discharge for the WRF CTRL and PGW periods. The line represents the mean and the shadow represents the standard deviation of the streamflow discharges over the 8 water years.

Figure 11 Change between WRF CTRL and PGW periods in the simulated mean Marmot Creek monthly streamflow discharges during March to October for the 8 water years.

3.6 Changes in streamflow

Simulated daily streamflow discharge was compared between CTRL and PGW in Fig. 10a, which shows the annual time series of Marmot Creek discharge for the CTRL and PGW periods. The results suggest that basin discharge will remain very similar over winter through mid-March. However, the average onset of spring freshet advanced forward 45 d, from 8 May for CTRL to 24 March with PGW, with the centre of flow volume occurring 12 d earlier, from 22 June for CTRL to 10 June with PGW. The peak basin discharge was 1.13 and 1.01 m³ s⁻¹ in CTRL and with PGW, respectively, both on 21 June. Daily streamflow discharge during the recession after the date of peak discharge through to late August declined with PGW compared to the CTRL period, because of higher evaporation losses with PGW (Fig. 7c). Despite the lower peak and faster recession, the cumulative annual discharge volume increased with PGW by 18 % from 3973 dam³ in CTRL to 4683 dam³ with PGW (Fig. 10b). The change in monthly discharge was calculated by subtracting the monthly discharge with PGW from that in the CTRL period. Figure 11 shows noticeable increases in monthly discharge with PGW for March to May and September, with PGW increases ranging from 0.02 m³ s⁻¹ in March to 0.27 m³ s⁻¹ in May compared to CTRL monthly discharges of 0.01 m³ s⁻¹ in March and 0.13 m³ s⁻¹ in May. In contrast, monthly discharge with PGW declined notably in June and July by 0.03 and 0.09 m³ s⁻¹ from CTRL values of 0.69 and 0.29 m³ s⁻¹, respectively. Monthly discharge in PGW decreased by only 0.01 m³ s⁻¹ in August from the CTRL value of 0.13 m³ s⁻¹ and was a similarly low value (0.06 m³ s⁻¹) in October. Monthly fractional changes in discharge with PGW ranged from a 573 % increase in April to a 33 % decrease in July; seasonal fractional changes in discharge ranged from a 236 % increase in spring (March to May), a 12 % decline in summer (June to August), to a 13 % increase in fall (September to October).

Figure 12 Simulated mean annual daily runoff for WRF CTRL and PGW. (a) Alpine, (b) treeline, (c) upper forest, (d) forest clearing blocks, (e) forest circular clearing north-facing, (f) forest circular clearing south-facing, and (g) lower forest ecozones. The line represents the mean and the shadow represents the standard deviation of runoff over the 8 water years.

Figure 13 Simulated mean annual cumulative runoff for WRF CTRL and PGW. (a) Alpine, (b) treeline, (c) upper forest, (d) forest clearing blocks, (e) forest circular clearing north-facing, (f) forest circular clearing south-facing, and (g) lower forest ecozones. The line represents the mean and the shadow represents the standard deviation of the 8 water year runoff.

The simulated daily runoff fluxes (surface, sub-surface, and groundwater runoff) and annual runoff volumes were plotted for all ecozones in MCRB to examine changes between CTRL and PGW. Figures 12 and 13 consistently show minimal change in winter months and an advance in the onset of the spring freshet for all ecozones with PGW. The annual peak runoff from the alpine ecozone decreased from 25.6 mm d⁻¹ in CTRL to 23.2 mm d⁻¹ with PGW, both occurring on 20 June (Fig. 12a), whilst the greatest decline in annual peak runoff occurred from the treeline ecozone, from 27.8 mm d⁻¹ in CTRL to 19.5 mm d⁻¹ with PGW, on 21 and 17 June, respectively (Fig. 12b). There were moderate declines in annual peak runoff with PGW from other ecozones, ranging from 0.3 mm d⁻¹ less at the lower forest ecozone to 2.0 mm d⁻¹ less at the forest circular clearing north-facing ecozone. The change in dates of annual peak runoff with PGW ranged from no change at the forest clearing north-facing ecozone to 9 d later at the forest clearing south-facing ecozone (Fig. 12c–g). There were moderate increases in annual runoff volumes from the forest clearing and lower forest ecozones with PGW, ranging from a 12 dam³ increase from the lower forest ecozone to a 17 dam³ increase from the forest clearing blocks ecozone (Fig. 13d–g). The annual runoff volume from the upper forest ecozone increased from 258 dam³ in CTRL to 316 dam³ with PGW (Fig. 13c). The alpine and treeline ecozones are the primary sources for Marmot Creek discharge; here the annual runoff volume increased 25 %, from 2457 dam³ in CTRL to 3065 dam³ with PGW, from the alpine ecozone, but decreased 2 % from 1007 dam³ in CTRL to 986 dam³ with PGW from the treeline ecozone (Fig. 13a and b).

The relationship between rainfall ratio (RR) and runoff efficiency (RE) was examined for all ecozones and the Marmot Creek basin for the CTRL and PGW periods. The RR is the cumulative rainfall divided by cumulative precipitation over a WY, and RE is the cumulative runoff (surface, subsurface, and groundwater runoff) divided by cumulative precipitation for a WY. A RR > 0.5 indicates a rainfall-dominated precipitation regime, and a RR < 0.5 indicates a snowfall-dominated precipitation regime. The RE describes how well the basin or ecozone converts precipitation volumes to runoff volumes, and it normally varies between 1 and 0. Figure 14 illustrates the changes between CTRL and PGW in mean RR and RE values for the ecozones and the basin. The RR in CTRL was 0.43 for the alpine ecozone, meaning it is snowfall-dominated, and for the treeline ecozone it was 0.51, close to an equal snowfall and rainfall precipitation regime. For other ecozones, the RR values in CTRL were between 0.6 and 0.62, indicating rainfall dominance. With PGW, the RR increased to range from 0.61 at the alpine ecozone to 0.76 at the lower forest ecozone. For the basin, the RR rose from 0.52 in CTRL (i.e. close to equal snowfall and rainfall) to 0.68 with PGW (i.e. rainfall-dominated). The RE stayed relatively unchanged for the forest ecozones, ranging from 0.1 for the lower forest ecozone to 0.11 for the upper forest ecozone in both CTRL and PGW periods. For the forest clearing ecozones, RE values dropped by 0.02. The RE changed substantially in the alpine and treeline ecozones, dropping from 1.04 in CTRL to 0.91 with PGW for the treeline ecozone but increasing from 0.62 in CTRL to 0.69 in PGW for the alpine ecozone. The value near 1 in CTRL for treeline ecozones refers to melt of late-lying snow patches. For the entire basin, the RE increased by only 0.01 with PGW, from 0.44 in CTRL to 0.45 in PGW despite the basin shifting towards domination by rainfall runoff.

Figure 14 Mean rainfall ratio and runoff efficiency for WRF CTRL and PGW for alpine (AE), treeline (TE), upper forest (UF), forest clearing blocks (FCB), forest circular clearing north-facing (FCCNF), forest circular clearing south-facing (FCCSF), and lower forest (LF) ecozones and the entire Marmot Creek basin.

4 Discussion

The results suggest that the bias-corrected 4 km resolution WRF CTRL near-surface meteorological variables were comparable to detailed station observations of a well-instrumented mountain basin, and the hydrological simulations driven by these variables could make acceptable predictions of seasonal snow accumulation, snowmelt, and streamflow that compared well to field observations. WRF can be run at higher resolutions than 4 km, which may provide even more realistic precipitation patterns, especially for extreme events (Tao et al., 2016; Li et al., 2017). However, these higher-resolution WRF simulations require greater computational capacity, and so the 4 km WRF simulations are currently considered the best options to assess future climate change for large regions (Liu et al., 2017; Prein et al., 2017; Warscher et al., 2019).

For the business-as-usual scenario (RCP8.5) with PGW, MCRB on average warmed up by 4.7 ^∘C, and this degree of warming is comparable to other findings for this region by the end of the 21st century (Nogués-Bravo et al., 2007; Kienzle et al., 2012). As air warms, the water-holding capacity of the atmosphere increases based on the Clausius–Clapeyron equation, resulting in 13 % to 20 % greater precipitation to MCRB in the projected PGW period. This combination of a warmer and wetter future climate produces complex hydrological responses. MCRB is composed of land covers that vary with elevation (ecozones), for which there are different dominant hydrological processes in each ecozone, leading to variable hydrological responses to the warming climate. Annual precipitation shifted towards more rainfall and less snowfall for all ecozones, with larger changes to alpine and treeline ecozones, and MCRB shifted from a balance of snowfall and rainfall towards rainfall domination with PGW. As a result of decreased snowfall, the seasonal snowpack dropped substantially for all ecozones with PGW, even though sublimation losses from blowing snow in the alpine ecozone and intercepted snow in the forested ecozones also decreased. The treeline ecozone developed large snow drifts and so sustained the highest seasonal snow accumulation in MCRB and suffered the greatest loss due to reduced snowfall and diminished blowing snow transport. These results are consistent with modelling results for warmer climates in other mountain basins where the model has a detailed representation of snow processes (López-Moreno et al., 2013; Pomeroy et al., 2015; Marty et al., 2017). In contrast to the ubiquitous slower snowmelt rates predicted by Musselman et al. (2017), snowmelt rates exhibited a mixed response to PGW, with lower melt rates for the treeline and forest clearings ecozones, a slightly higher melt rate for the alpine ecozone, unchanged melt rates for the forest ecozones, and a slightly lower melt rate over the whole basin. These melt rate responses are due to the complex interaction amongst changes in precipitation, air temperature, albedo decay, and radiation fluxes to snow cover with PGW. On the other hand, all ecozones at MCRB experienced earlier depletion of the seasonal snowpack and a shorter snow-cover duration with PGW; as a result, the onset of snowmelt runoff and timing of Marmot Creek streamflow shifted earlier towards spring. Marmot Creek streamflow declined during summertime with PGW because of higher evaporation loss and earlier water availability. The shift in snowmelt runoff, snow-cover depletion, and decline in summer discharge have been found for many snow-dominated mountain basins (Stewart et al., 2004; Barnett et al., 2005; Kienzle et al., 2012; Jepsen et al., 2016) and suggest that the water management of mountain basins will be severely impacted (Rood et al., 2005; Chen et al., 2006; Bongio et al., 2016). The results here suggest an overall increase in discharge volume, which contrasts with earlier studies suggesting declines in Canadian Rockies streamflow in the 21st century (Rood et al., 2005; Tanzeeba and Gan, 2012). However, earlier studies used much coarser climate models to inform hydrological predictions and model outputs were not verified against current high-altitude observations.

It is interesting to note the variability of runoff response from different ecozones in MCRB to PGW. Despite an earlier onset of runoff for all ecozones, the annual runoff volume increased only very moderately for the forest clearings and lower forest ecozones and moderately for the upper forest ecozone with PGW. Increases in precipitation, particularly rainfall, with PGW were almost completely consumed by increases in evaporation from these densely vegetated ecozones. At higher elevations where most runoff is produced, the results show opposing responses between the alpine and treeline ecozones. For the alpine ecozone, annual runoff volume increased with PGW; this is caused by a combination of limited available subsurface storage and increases in rainfall that overwhelmed the increased evaporation loss. Increases in surface runoff more than compensated for the slightly reduced subsurface runoff from this ecozone with PGW. In contrast, the annual runoff volume decreased for the treeline ecozone, caused by the following factors: a large reduction in seasonal snow accumulation due to decreased snowfall and suppressed blowing snow transport, decreased runoff from reduced snowmelt, higher evaporation in the longer snow-free period that more than compensated for increased rainfall, and more available subsurface storage, subsequently leading to less flow from subsurface runoff. However, the decline in runoff from the treeline ecozone was compensated for by an increase in runoff from the alpine ecozone, which contributed to an overall increase in annual basin streamflow discharge of 18 % for a 16 % increase in precipitation with PGW. The alpine and treeline ecozones comprise about 45 % of MCRB and receive more than half of the basin's annual precipitation, so their runoff response to PGW dominates changes in basin streamflow. Interestingly, MCRB has shown hydrological resilience to past changes in climate, forest disturbance, and extreme events (Harder et al., 2015); it also exhibits some resilience to the projected future climate in respect to its RE. That is, despite the basin shifting towards rainfall-runoff dominance with PGW, the RE values are almost unchanged.

The hydrological model set up in CRHM provided a physically and process based analysis of climate change effects on the hydrological changes at MCRB by including the relevant hydrological processes for streamflow generation. These processes were parameterised from local research results without calibration from streamflow (Pomeroy et al., 2009, 2016; Ellis et al., 2010; MacDonald et al., 2010; Harder and Pomeroy, 2013; Fang et al., 2013; DeBeer and Pomeroy, 2017). This is a strength of the modelling approach in this study. However, the CRHM model does have weaknesses in simulating groundwater systems and groundwater–surface interactions in the hillslope module given the insufficient information available to describe and fully parameterise all of the processes that may influence this system (Fang et al., 2013). The parameterisation in the hillslope module assumes surface tension in the pore spaces is the dominant factor, while in mountain basins, coarse-textured and unconsolidated materials over underlying bedrock or impeding layers could be found in subsurface and groundwater layers. This may result in preferential flow regimes influencing the groundwater–surface interactions in ways that are incompletely understood and described in this model (Sidle et al., 2000; McClymont et al., 2011). There is an ongoing effort in CRHM to improve the parameterisation for groundwater–surface interactions. Thus, some caution should be taken when interpreting changes in groundwater storage and flow.

In this study, the assessment of future changes in hydrology at MCRB held land-cover and soil properties static between the CTRL and PGW periods. Caution should therefore be taken in interpreting this study's results, as land-cover and soil properties could also change in the warmer climate (Rasouli et al., 2019). Modelling studies in the southern Rocky Mountains, USA, and the eastern European Alps examine future changes in basin hydrology and water availability by coupling climate change and land-cover disturbance and suggest that the feedbacks of land-cover disturbance with a warmer climate should be considered when assessing future impacts on mountain hydrology (Buma and Livneh, 2015; McDowell et al., 2016; Bennett et al., 2018; Strasser et al., 2019). However, land-cover changes such as forest expansion into the alpine ecozone under a warmer climate are less likely in the Canadian Rockies than elsewhere because of the limitations imposed by geological, microclimate, and geomorphological processes on tree establishment above current treelines (Macias-Fauria and Johnson, 2013). In the Canadian Rockies, greater collaboration amongst hydrologists, ecologists, micrometeorologists, and soil scientists is warranted to examine the regional changes in vegetation and soil for the future climate and to come up with more meaningful scenarios of land cover and soil disturbance under changing climate for a more comprehensive assessment of future mountain hydrology.

5 Conclusions

A physically based hydrological model was set up using the CRHM platform and was forced by bias-corrected 4 km WRF near-surface meteorology outputs in CTRL (i.e. 1 October 2005 to 30 September 2013) and PGW (i.e. 1 October 2091 to 30 September 2099) periods to assess the future changes in hydrology at Marmot Creek Research Basin. Model simulations using the bias-corrected WRF outputs over the CTRL period provided reasonable predictions for snow accumulation and snowmelt for forest and alpine sites compared to the field observations. Model streamflow simulations were also evaluated against observed streamflow discharge, which showed adequate predictability. Compared to the CTRL period, Marmot Creek on average warmed up by 4.7 ^∘C and received 16 % more precipitation for the business-as-usual climate warming scenario; as a result, the rainfall ratio rose above 0.5, with a 268 mm increase in rainfall and 112 mm decrease in snowfall. However, changes in basin hydrology were more complex than those in precipitation, as the basin is composed of seven ecozones ranging from forest to alpine tundra and spanning a large elevational gradient. Under climate change, all ecozones developed lower snowpacks, with decreases in peak SWE ranging from 11 mm for the upper and lower forest ecozones to 149 mm for the treeline ecozone, shorter snow-cover duration ranging from 31 d for the forest circular clearing north-facing ecozone to 49 d for the treeline ecozone, and total seasonal snowmelt volume decreasing from 32 mm for the upper forest ecozone to 215 mm for the treeline ecozone. These changes were impacted by several cold region hydrological processes, including suppressed blowing snow transport and sublimation for the alpine and treeline ecozones and reduced sublimation from smaller canopy snow interception for forest ecozones. With climate change, snowmelt rates declined by 1.1 mm d⁻¹ for treeline ecozones and by 0.9 to 1.6 mm d⁻¹ for forest clearing ecozones, but increased by 0.1 mm d⁻¹ for the alpine ecozone and remained unchanged for forest ecozones, showing that slower snowmelt is not ubiquitous with warming. Additionally, evaporation losses increased with climate change from 59 mm for the alpine ecozone to 179 mm for the upper forest ecozone, with a 124 mm greater evaporation loss for the whole basin. Runoff responses to climate change were even more variable in the different Marmot Creek ecozones. Runoff began earlier for all ecozones with climate change, and the forested and alpine ecozones generated annual runoff volumes that were from 12 % to 25 % larger, whereas runoff from the treeline ecozone declined by 2 %. For the whole basin, the declining runoff from the treeline ecozone was more than compensated by higher runoff from the other ecozones, resulting in an 18 % increase in basin streamflow discharge with climate change. Given the 16 % increase in precipitation, the runoff efficiency for Marmot Creek is virtually unchanged under climate change despite the hydrological shift from snowmelt runoff towards rainfall-dominated runoff.

Appendix A: Abbreviation list

a.s.l.	above sea level
CH	Centre for Hydrology
CMIP5	Coupled Model Intercomparison Project fifth phase
CRHM	Cold Regions Hydrological Modelling platform
CTRL	control from WRF
dam3	cubic decametre
GCMs	Global climate models
HRUs	hydrological response units
IPCC	Intergovernmental Panel on Climate Change
MAD	mean absolute difference
MB	model bias
MCRB	Marmot Creek Research Basin
NRMSD	normalised RMSD
NSE	Nash–Sutcliffe efficiency
PGW	pseudo global warming from WRF
Q–Q	quantile–quantile
QDM	quantile delta mapping
RCM	regional climate model
RCP8.5	Representative Concentration Pathway 8.5
RE	runoff efficiency
RMSD	root mean square difference
RR	rainfall ratio
SRB	Saskatchewan River basin
SWE	snow water equivalent
WRF	Weather Research and Forecasting
WSC	Water Survey of Canada
WY	water year