We describe a state-of-the-art framework for projecting hydrologic impacts due to enhanced warming and amplified moisture fluxes in the subarctic environment under anthropogenic climate change. We projected future hydrologic changes based on phase 5 of the Coupled Model Intercomparison Project global climate model simulations using the Variable Infiltration Capacity hydrologic model and a multivariate bias correction/downscaling method for the Liard basin in subarctic northwestern Canada. Subsequently, the variable importance of key climatic controls on a set of hydrologic indicators was analyzed using the random forests statistical model. Results indicate that enhanced warming and wetness by the end of century would lead to pronounced declines in annual and monthly snow water equivalent (SWE) and earlier maximum SWE. Prominent changes in the streamflow regime include increased annual mean and minimum flows, earlier maximum flows, and either increased or decreased maximum flows depending on interactions between temperature, precipitation, and snow. Using the variable importance analysis, we find that precipitation exerts the primary control on maximum SWE and annual mean and maximum flows, and temperature has the main influence on timings of maximum SWE and flow, and minimum flow. Given these climatic controls, the changes in the hydrologic indicators become progressively larger under the scenarios of 1.5°, 2.0°, and 3.0°C global mean temperature increases above the preindustrial period. Hence, the framework presented in this study provides a detailed diagnosis of the hydrologic changes as well as controls and interactions of the climatic variables, which could be generalized for understanding regional scale changes in subarctic/nival basins.
Freshwater systems in Arctic and subarctic environments are changing rapidly in response to a warming climate and changing precipitation regime. Enhanced warming in the Arctic—which is almost twice as large as the global temperature average (Serreze and Barry 2011; Cohen et al. 2014)—coupled with amplified poleward moisture transport (Min et al. 2008; Zhang et al. 2013) is affecting different components of the hydroclimatic system (Larsen et al. 2014; Bring et al. 2016), including reductions in the duration and extent of snow cover (Liston and Hiemstra 2011; Derksen et al. 2015) and enhanced permafrost thaw (Romanovsky et al. 2010; Grosse et al. 2016). Intensification of the hydrologic regime is an expected manifestation of the enhanced warming and wetness in the region. For instance, increased precipitation is causing increased streamflow across pan-Arctic river basins (Peterson et al. 2002; Wu et al. 2005; Rawlins et al. 2010; Zhang et al. 2013), while peak streamflow is shifting earlier (Burn et al. 2010; Overeem and Syvitski 2010; Tan et al. 2011). Further, warmer winters result in a larger fraction of winter precipitation as rainfall, increased infiltration and subsurface water movement, which is leading to increased winter flows (Smith et al. 2007; Tananaev et al. 2016).
Future climate projections suggest continuations of enhanced warming and strong precipitation increase (Collins et al. 2013; Vihma et al. 2016). For instance, December–February and June–August mean air temperature over land across the pan-Arctic are projected to increase by 10°C (6°C) and 5°C (4°C), respectively, under phase 5 of the Coupled Model Intercomparison Project (CMIP5) representative concentration pathway RCP8.5 (RCP4.5) scenario by the end of the twenty-first century (CMIP5 multimodel mean for 2081–2100 relative to 1986–2005 averages; van Oldenborgh et al. 2013). Further, mean precipitation over land for October–March and April–September are projected to increase by 60% (35%) and 25% (20%), respectively, under the RCP8.5 (RCP4.5) scenario by the end of the twenty-first century (van Oldenborgh et al. 2013).
The ongoing hydrologic change in the Arctic/subarctic is expected to be amplified due to the enhanced future warming and wetness. Such alterations in the hydrologic regime, together with the loss of snowpack and permafrost, could have major implications on infrastructure and services in the regions (Larsen et al. 2014; Instanes et al. 2016). However, the regions have received limited focus in terms of basin-scale hydrologic impact studies. Typically, in other regions, sophisticated basin-scale hydrologic models forced with downscaled outputs from GCMs are used to project future hydrologic changes. For the Arctic and subarctic regions, however, there is still a need to rely on future projections based on global or pan-Arctic wide studies (e.g., Sperna Weiland et al. 2012; Koirala et al. 2014; Schewe et al. 2014; Bring et al. 2017). These large-scale studies certainty do have value, as their projections (e.g., increased annual flows and low flows, and earlier peak flows) provide a basic understanding of the hydrologic changes due to changes in the large-scale climate drivers. Yet, for river basin-scale assessments, spatial resolutions (≥0.5°) of these studies are often too coarse to represent the relevant physical processes and spatial heterogeneity, especially for the mountainous regions and small subbasins. For instance, Melsen et al. (2016) and Buitink et al. (2019) compared hydrologic models at different spatial resolutions for the Swiss Alps region, and found substantial improvement in the model performance and representation of intrabasin variability when using fine-resolution models. Further, in large-scale studies, meteorological forcing variables from GCMs are often used without downscaling, and model parameters are usually not calibrated to replicate basin-scale responses, which lead to large uncertainties in hydrologic projections (Krysanova et al. 2017; Hattermann et al. 2017). Some of these shortcomings have been addressed by recent basin-scale studies. For instance, Scheepers et al. (2018) employed an ensemble of statistically downscaled CMIP5 GCMs and a basin-scale hydrologic model for projecting changes in the Mackenzie River basin, western Canada. Although projected future changes—for example, earlier onset of spring snowmelt, higher winter and spring runoff, and lower summer runoff—are similar to results from large-scale studies, basin-scale studies are, arguably, more credible because climate model biases are accounted for explicitly, and snow and runoff processes are parameterized in a more realistic fashion.
Further, there have been recent innovations in the bias correction/downscaling methods used to obtain finer-scale, bias-adjusted representations of the GCM climate forcings. These include multivariate techniques, which, in contrast to commonly used univariate downscaling methods that only correct distributions of individual variables, also correct biases in intervariable dependence (Cannon 2018a). Given that the cold-regions hydrologic processes, such as snow accumulation and melt are highly dependent on temperature and precipitation interactions (Luce et al. 2014), the multivariate techniques by virtue of maintaining realistic relationships between climatic variables could potentially lead to improved hydrologic simulations. For instance, Meyer et al. (2019) demonstrated improvements in hydrologic simulations in terms of a more plausible rainfall–snowfall partitioning, and a better agreement with the observed snow water equivalent, glacier volume, and streamflow regime, when they compared results from a multivariate downscaling method with those from a univariate method. Thus, as suggested by Zscheischler et al. (2019), multivariate bias correction methods should be favored for impacts modeling when interactions between variables play an important role (such as for cold-regions hydrologic modeling) to ensure that the climatic variable dependencies are captured more realistically.
With respect to future hydrologic responses, changes in temperature and precipitation can affect different components of the hydrologic cycle. For instance, in snow-dominated regions, temperature changes mostly affect runoff timing while precipitation changes mostly affect runoff volume (Barnett et al. 2005, 2008). Further, the magnitude and seasonality of snowpack storage are sensitive to seasonal precipitation and temperature, with winter/spring temperature increases causing reductions that could be offset by precipitation increases (Mote 2006; Luce et al. 2014). However, the role of different climatic controls on key hydrologic indicators (e.g., annual flow, maximum and minimum flow and timings) remain to be analyzed in many basins. Such analyses could add value, not only by identifying the dominant climatic controls, but also by providing a general understanding of the systems behavior (e.g., control versus response) with potential usage beyond the studied basin. Input–output sensitivity can be investigated using variable importance techniques (Wei et al. 2015), which typically provide a measure of the relative importance of control variables (e.g., climatic drivers) on response variables (e.g., hydrologic indicators).
Further, with the signing of the 2015 Paris Agreement, which aims to limit global warming well below 2°C and to pursue efforts to limit it to 1.5°C above the preindustrial level (UNFCCC 2015) and subsequent release of the Intergovernmental Panel on Climate Change (IPCC) special report on global warming of 1.5°C (IPCC 2018), understanding local-scale impacts under different global warming thresholds has become a policy relevant issue. In this context, the results of VI analysis could be extended to link basin-scale temperature and precipitation changes at different global warming thresholds to hydrologic changes.
The study contributes to an improved understanding of hydrologic impacts of anthropogenic climate change in the subarctic environment by using a state-of-the-art modeling framework, consisting of a detailed hydrologic model, multivariate bias correction/downscaling method, and a statistical variable importance model. Specifically, in order to develop a better understanding of the river basin/subbasin scale future hydrologic responses due to enhanced warming and wetness in the subarctic environment, we employed the 1/16° grid resolution Variable Infiltration Capacity (VIC) hydrologic model (Liang et al. 1994; Hamman et al. 2018) for the Liard basin in northwestern Canada. We employed the N-dimensional multivariate bias correction algorithm (MBCn; Cannon 2018a) to downscale/bias-correct CMIP5 climate models (Taylor et al. 2012) under the RCP4.5 and RCP8.5 scenarios, and subsequently drive the VIC model. Further, we applied the random forests statistical model (Breiman 2001) to evaluate the dominant controls of temperature and precipitation changes on a set of hydrologic indicators. Additionally, we analyzed the interactions of climatic controls with hydrologic indicators under different global warming thresholds above the preindustrial period.
2. Study basin
The Liard is a large river system in the subarctic northwestern Canada and a major tributary of the Mackenzie River (Fig. 1). Covering about 16% of the Mackenzie basin with a drainage area of approximately 277 000 km2, the Liard River contributes about 25% of total annual discharge. The river and its tributaries flow through four Canadian provinces/territories (Yukon, British Columbia, Alberta, and Northwest Territories) before discharging into Mackenzie River at Fort Simpson. Major tributaries of the Liard River include the Fort Nelson River, South Nahanni River, Kechika River, Petitot River, Muskwa River, Dease River, and the Frances River.
Located in the Boreal and Taiga ecoregions in western Canada, the Liard River basin (LRB) is mostly in pristine state, hence this study provides an assessment of climate-induced hydrologic changes in a basin with minimum human impacts. The LRB is physiographically heterogeneous with three northwest–southeast trending Cordillera Mountains chains (Stikine Ranges, Selwyn, and Rocky Mountains) occupying more than 50% of the basin, and the rest in the interior plains. Elevation exceeds 3200 m in the mountains and ranges between 170 and 800 m in the southern plateaus and lowlands (Woo and Thorne 2006). Isolated to extensive discontinuous permafrost is present in parts of the basins (Heginbottom et al. 1995).
The region is dominated by westerly circulation, with topography exerting a major influence on atmospheric circulation, especially over the Cordillera Mountains causing a strong precipitation gradient (Szeto et al. 2008). Consequently, the southwestern mountainous part of the basin receives higher precipitation compared to the eastern leeward side (Woo and Thorne 2006). Mean annual precipitation, temperature, and runoff in the basin over the years 1979–2012 were about 570 mm, −2.0°C, and 290 mm, respectively [based on Pacific Northwest North American meteorological (PNWNAmet) gridded climate data (Werner et al. 2019) and Water Survey of Canada (WSC), Liard River at Mouth hydrometric station streamflow data]. Seasonally, the basin receives higher precipitation between April and September and lower precipitation between October and April when the basin mean temperature are below freezing (Fig. S1 in the online supplemental material). Runoff is characterized by a subarctic nival regime, in which snowmelt dominates generating annual high flows in combination with summer rainfall (Woo and Thorne 2006).
Recent studies suggest that there have been changes in the hydrologic response of the basin. For instance, Burn (2008) evaluated the spring peak streamflow timing for the years 1966–2005 and found significantly earlier shifts (at 10% significance level; p < 0.1) at several sites in the LRB. The winter baseflow contribution to the Liard River at Mouth streamflow was also found to have increased significantly (p < 0.1) over the years 1973–2007, while the changes in annual flows remained insignificant (St. Jacques and Sauchyn 2009). Further, Rood et al. (2016) found a significantly increasing long-term (1944–2013) (p < 0.05) trend in the Liard River at Fort Liard mean annual discharge.
Future climate projections for the LRB indicate continued enhanced warming and strong increases in precipitation, with hydrologic changes generally consistent with the historical trends. Specifically, Dibike et al. (2017) used an ensemble of six statistically downscaled CMIP5 GCMs (2071–2100 versus 1971–2000) and projected an about 3°–7°C (1°–4°C) increase in annual mean temperature and 20%–40% (10%–25%) increase in annual mean precipitation for the RCP8.5 (RCP4.5) scenario. They further analyzed the annual and June–August standardized precipitation and evapotranspiration indices, and found a gradual surplus in the annual water supply (expressed as precipitation minus potential evapotranspiration) and a progressively larger deficit in the summer water supply over the years 1950–2100. Hydrologic projections for the LRB are available in the studies by Thorne (2011) and Gosling et al. (2011), in which they evaluated future changes corresponding to 2°C prescribed warming for six GCMs and 1°–6°C warming for one GCM, with climate inputs derived by using a weather generator technique (Todd et al. 2011). Their hydrologic projections showed an earlier spring freshet, increased autumn to spring and annual discharge, and decreased summer discharge. Notably, neither future hydrologic projections with respect to CMIP5 RCP forcing scenarios and global warming thresholds nor evaluations of climate controls on hydrologic changes have been studied for the LRB.
3. Methods and data
a. GCM data and downscaling
We selected an ensemble of seven GCMs (summarized in Table S1) that participated in the CMIP5 experiment (Taylor et al. 2012), with two scenarios that correspond to high (RCP8.5) and moderate (RCP4.5) radiative forcings for each GCM. Five of the seven GCMs were used in the North American Coordinated Regional Climate Downscaling Experiment (Mearns et al. 2017), which favored climate models with a long, well-documented development history and record of credible simulations (Kotamarthi et al. 2016). Hence, uncertainties due GCM structure and greenhouse gas concentration and anthropogenic forcings, which are the most important sources of uncertainties in projecting hydrologic impacts of climate change (Najafi et al. 2011; Bennett et al. 2012; Hattermann et al. 2018), are considered in this study.
We employed a combination of spatial interpolation and the N-dimensional multivariate bias correction method (MBCn) (Cannon 2018a) to spatially disaggregate and bias correct the coarse-resolution GCM outputs to match the resolution (1/16°) and statistical characteristics of the gridded observations used to calibrate the hydrologic model. MBCn is a multivariate extension of quantile mapping that corrects biases in the marginal distributions of multiple climate model variables as well as the dependence between variables (i.e., the empirical copula). MBCn is an extension of an image processing technique (Pitié et al. 2007) that operates by iteratively (i) applying a random orthogonal rotation to both climate model and observational target datasets, (ii) correcting the marginal distributions via quantile mapping, and (iii) rotating datasets back to the original axes and checking convergence. Repeating these steps guarantees transfer of the climate model’s multivariate distribution to that of observations. Specifically, the marginal distributions and empirical copula in the historical calibration period will match observations. In the future projection period, projected changes in corrected quantiles are also constrained to match those of the raw climate model (Cannon et al. 2015). Extrapolation involves an approach called quantile delta mapping (Cannon et al. 2015) that preserves absolute changes in quantiles for temperature and relative changes in quantiles for precipitation. Further, the corrected empirical copula is free to evolve, and model-projected changes in multivariate dependence structure are retained by MBCn, although this depends on how much correction is needed to bring the historical empirical copula in line with observations. In terms of multivariate dependence, both features—correction of the historical copula and preservation of changes in dependence—can be important when considering future changes in multivariate extremes (Kirchmeier-Young et al. 2017; Zscheischler and Seneviratne 2017; Zscheischler et al. 2019), including those related to hydrologic processes in cold regions (Meyer et al. 2019). More details on the MBCn method are available in Cannon (2018a).
In this study, the MBCn R package (Cannon 2018b) was used to bias correct daily maximum and minimum temperature, daily precipitation amounts, and wind speed from GCMs over the LRB. First, GCM outputs were remapped onto the 1/16° observational grid using bilinear interpolation. Next, interpolated climate data were adjusted, on a grid point by grid point basis with respect to the PNWNAmet gridded climate data (Werner et al. 2019) using the 56-yr 1950–2005 period for calibration, with MBCn bias corrections applied over 3 × 19 years = 57-yr sliding windows to approximately match the length of the calibration period. Each time, the central 19 years were replaced, the window was slid 19 years, etc., until the end of the projection period was reached. Replacing the central 19 years helps avoid large discontinuities between edges of the full 57-yr sliding window and is computationally more efficient than, for example, replacing only the central year or other smaller block lengths. Within each 19-yr block—to ensure an unbiased seasonal cycle—bias corrections were applied over 3 × 30-day sliding day-of-year blocks, with the central 30 days replaced and the window slid 30 days at a time.
b. Hydrologic modeling
Given the mountainous terrain, fine-resolution hydrologic modeling is desirable for the LRB. Therefore, we employed the semidistributed macroscale VIC hydrologic model version 5.0.0 (Liang et al. 1994, 1996; Hamman et al. 2018) at a 1/16° spatial resolution for the LRB. Since its initial development in 1990s, the model has undergone a number of updates and refinements to better represent the cold climate processes, such as energy balance over snow and frozen ground (Cherkauer and Lettenmaier 2003; Andreadis et al. 2009). With these refinements, the model has been extensively applied in the snow dominated regions (e.g., Andreadis et al. 2009; Shrestha et al. 2012, 2014; Werner et al. 2013; Schnorbus et al. 2014), and over the pan-Arctic domain at a coarse resolution (Adam and Lettenmaier 2008; Tan et al. 2011; Hamman et al. 2016). Also, VIC has been widely used for assessing the hydrologic impacts of climate change (Hidalgo et al. 2009; Elsner et al. 2010; Shrestha et al. 2012; Werner et al. 2013; Schnorbus et al. 2014).
We set up the VIC model with geospatial data derived from: (i) 7.5-arc-s Global Multiresolution Terrain Elevation Data 2010 (Danielson and Gesch 2011), (ii) 250-m resolution land cover dataset from North American land change monitoring system (Latifovic et al. 2012), and (iii) 5-arc-min resolution soil classification and parameterization based on the data provided in the Global Soil Data Products CD-ROM (Global Soil Data Task 2014). The geospatial data were processed and gridded to match the 1/16° resolution of the VIC model setup, with elevation bands discretized at 200-m intervals. The LRB was subdivided into 28 subbasins based on the Water Survey of Canada hydrometric station network for the representation of spatially heterogeneous hydrologic responses. We used daily precipitation, maximum and minimum temperature, and wind speed from the PNWNAmet gridded climate data, which match the 1/16° resolution of the VIC model setup, for model calibration. The 3-hourly meteorological inputs of precipitation, maximum and minimum air temperature, wind speed, longwave radiation, shortwave radiation, atmospheric pressure, and vapor pressure data required to run VIC in energy balance mode were generated from the daily data using the Mountain Microclimate Simulation Model, version 4.3 (MTCLIM 4.3; Thornton and Running 1999) from VIC version 4.2.
We employed the Nondominated Sorting Genetic Algorithm (NSGA-II; Deb et al. 2002) for a multiobjective calibration of the VIC model parameters to match observed discharge characteristics, with three commonly used goodness-of-fit measures calculated at daily time steps used as objective functions: (i) Nash–Sutcliffe coefficient of efficiency (NSE), (ii) NSE of log-transformed discharge (LNSE), and (iii) volume bias (VB). Five standard runoff generation parameters were calibrated: variable infiltration curve parameter (Bi), fraction of maximum soil moisture where nonlinear baseflow occurs (Ws), maximum velocity of baseflow (Dsmax), fraction of Dsmax where nonlinear baseflow begins (Ds), and variation of saturated hydraulic conductivity with soil moisture (Expn). Ten years (1984–93) of observed discharge was used for model calibration and a further 10 years (1994–2003) was used for model validation, with two additional years (1982–83) used as a spinup period to exclude the effects of initial conditions from the hydrologic simulations. The NSGA-II optimization involved iterating a population size of 80 over 20 generations. From the multiobjective Pareto solutions, the model parameters with the best overall performance were selected using the fuzzy preference selection methodology (Shrestha and Rode 2008). Additionally, the VIC model snow simulations were compared with a number of snow pillow observations. After satisfactory model calibration, VIC was forced with the MBCn downscaled climate data for simulations corresponding to the historical and future GCM runs.
As previously stated, the LRB lies in the isolated to extensive discontinuous permafrost zone. Hence, the application of VIC with frozen ground/permafrost is the preferred means of simulating the hydrologic response. In practice, however, VIC model runs with frozen ground/permafrost turned on are computationally very expensive for large basins. We found in our tests that a VIC frozen ground run takes up to 1 h per grid cell for a 15-yr period. At this rate, model simulation for the LRB with 11 875 grids cells and 150 years of downscaled GCM data will take approximately 50 days to run on a 100-core processor cluster. Additionally, each subbasin needs to be calibrated, which involves running the model 1600 times. This makes it impractical to calibrate the model and run 14 GCMs with explicit frozen-ground processes. For this reason, model runs were limited to standard energy-balance mode. Thus, the lack of frozen ground in the hydrologic model leads to structural uncertainties, especially for subsurface-flow and low-flow simulation. However, the uncertainty in streamflow simulation is partly compensated by model calibration.
c. Variable importance analysis
Historical and future changes in temperature and precipitation, and their interactions, affect different components of the hydrologic cycle differently. Hence, we evaluated the response variables (i.e., hydrologic indicators) as a function of climate controls (seasonal temperature and precipitation) using the random forest (RF) ensemble machine learning method (Breiman 2001). The RFs are combination tree predictors, constructed with each tree using a random vector from bootstrap samples with the same distribution for all trees in the forest (Breiman 2001). We chose to use the RF model because of its simplicity and ability to account for interactions between variables, which is an important consideration in nival basins, for example, for snow accumulation and melt processes that are sensitive to temperature and precipitation interactions.
The RF model used in this study is based on (i) a large ensemble of regression trees (typically more than 500) in which a randomly chosen subset of predictors is used to find the best split at each node; (ii) a bootstrap sample of training cases, which increases the diversity between trees; (iii) growth of each tree to its maximum size without pruning; (iv) averaging predictions from the individual trees; and (v) using out-of-bag samples to estimate the generalization errors, thus eliminating the need for a separate cross validation. Averaging predictions from a large number of weakly accurate but diverse trees leads to an ensemble model that is robust against overfitting and performs well with little tuning (Breiman 2001).
The RF model provides a measure of sensitivities of the driving variables in terms of their relative variable importance (VI), as well as estimates of predictability of the response variables in terms of goodness-of-fit statistics. The algorithm estimates VI by considering the changes in prediction errors by randomly permuting the out-of-bag data for an individual variable while leaving all other variables unchanged (Liaw and Wiener 2002). The errors are computed on the out-of-bag data for each tree and after permuting a variable. The averaged differences (for all trees) normalized by the standard error gives the measure of importance for that particular variable. We used NSE for evaluating the model fit. Readers are referred to Breiman (2001) and Liaw and Wiener (2002) for further technical details. A number of studies (e.g., Prasad et al. 2006; Peters et al. 2007; Povak et al. 2014; Wang et al. 2015) illustrate applications of the RF model in hydrology and ecology.
In this study, the “caret” (Kuhn et al. 2018) interface to the “randomForest” R package (Liaw and Wiener 2002) was used. We trained a separate RF model for the six VIC-model-simulated hydrologic indicators: (i) maximum snow water equivalent (SWEmax), (ii) timing of SWEmax, (iii) annual mean flow (Qmean), (iv) annual maximum flow (Qmax), (v) timing of Qmax, and (vi) annual minimum flow (Qmin). The average SWE values and routed flows corresponding to the Liard River at Mouth station were used. The input variables consisted of seasonal temperature and precipitation (obtained from the MBCn downscaled climate data averaged over the entire basin and seasons) for the months: January–March (JFM), April–June (AMJ), July–September (JAS), and October–December (OND). These months were used instead of the commonly used meteorological seasons (e.g., December–February) because of better performance of the RF model during an initial model setup, likely due to the fact that the dominant processes in the basin align better with these months (e.g., snow accumulation in the basin starts in October). Hence, eight seasonal temperature and precipitation (JFM_T, AMJ_T, JAS_T, OND_T, JFM_P, AMJ_P, JAS_P, OND_P) variables were defined as RF inputs for a combined total of 1540 simulation years—seven historical forcing GCM simulations (1950–2005 for each), seven RCP4.5 GCM simulations (2006–2100 for each), and seven RCP8.5 GCM simulations (2006–2100 for each).
4. Results and discussion
a. Hydrologic model calibration/validation
We selected results for the Liard River at Mouth (Liard-M) hydrometric station to represent the aggregated response over the entire basin. Further, using the high-resolution VIC model outputs, we evaluated model performance and future projections for the South Nahanni River above Virginia Falls (S-Nahanni-VF) and Toad River above Nonda Creek (Toad-NC) tributaries (both Reference Hydrologic Basin Network stations; Brimley et al. 1999), and Liard River at Upper Crossing (Liard-UC) on the main stem (Fig. 1). The modeled flow results, with the NSE, LNSE, and Kling–Gupta efficiency (KGE) values greater than 0.70 (Table 1), indicate good model fits. Additionally, volume errors < 8% indicate a good ability of the model to simulate the water balance. The model is able to reproduce the observed streamflow dynamics reasonably well (Fig. S2), particularly for the annual major hydrologic event driven by the spring snowmelt. However, there are some discrepancies in the modeled results, for instance overprediction (underprediction) of the Toad-NC (S-Nahanni-VF) subbasin peak flows. A major source of uncertainty in streamflow simulation is the representativeness of the precipitation and temperature in the PNWMAmet dataset, especially the sparse station network in the high-latitude region (Mekis and Vincent 2011). Additionally, uncertainties arise from the quality of the observed discharge data, especially during ice-covered low-flow period (Hamilton 2008), hydrologic model structure (e.g., lack of frozen ground algorithm) and parameters (e.g., calibration of parameters based on discharge data alone).
Additionally, the comparison of VIC simulated SWE with observations at several locations in the basin (Fig. 1) generally shows a good representation of the seasonal high and low SWE dynamics (Fig. S3). However, there are discrepancies; for example, the mean absolute errors in the simulated values are 86, 52, 28, and 33 mm for Sikanni Lake, Wade Lake, Frances River, and Hyland River stations, respectively. Besides meteorological inputs, other factors contributing to the errors could include mismatch in the location and elevation of the observation station with the model grid and elevation band, respectively, and measurement errors (e.g., SWE calculation from snow depth). Nevertheless, given the general agreement in the magnitude and dynamics of the simulated SWE and streamflow with observations, the calibrated model was considered suitable for projecting future climate scenarios.
b. Temperature and precipitation changes
Figure 2 depicts the MBCn downscaled historical and future GCM projections of mean annual temperature and precipitation, together with changes in rainfall and snowfall. The projections illustrate progressive warming for both RCPs, with the ensemble median annual temperature (shown by solid lines) increasing from −2.4°C in 1976–2005 to about +0.9°C (RCP4.5) and +3.7°C (RCP8.5) by the end of century (2071–2100). The increases for the GCM ensemble range between +1.3° and +4.6°C under RCP4.5, and between +3.2° and +7.6°C under RCP8.5. Furthermore, the median annual temperature is projected to cross the 0°C threshold around 2040 (2050) for RCP8.5 (RCP4.5), with stabilization after 2060 under RCP4.5 and continued increases under RCP8.5. The temperature projections generally follow historical trends, which, for the years 1945–2012, are increasing significantly (p < 0.05) at a rate of about +0.3°C decade−1 (based on Mann–Kendall trend analysis with iterative prewhitening; Zhang and Zwiers 2004). Future projections indicate significant increasing trends for seven (six) out of seven ensemble members of RCP8.5 (RCP4.5) GCMs over the period 2006–2100.
Strong precipitation increases are projected in the basin under both RCPs, with annual increases for the GCMs ranging from +4% to +18% for RCP4.5 and from +9% to +30% for RCP8.5 (2080s versus 1976–2005). While the change in mean annual precipitation over the historical period is not significant, future projections indicate significant increases for seven (three) RCP8.5 (RCP4.5) GCM ensemble members over the period 2006–2100. Based on partitioning of rain and snow by VIC, precipitation increases are projected to mostly occur as rain (Fig. 2c), with significant increases for all seven ensemble members for both RCP4.5 and RCP8.5. Despite increased precipitation, snowfall is generally projected to decline (Fig. 2d), which is indicative of the temperature control on snowfall fraction and duration. This is also reflected in a larger snowfall decline for RCP8.5, which is statistically significant for five out of seven ensemble members. In the case of RCP4.5, the decline in snowfall is significant for a single ensemble member, with nonsignificant increases or decreases for three members each.
Following the IPCC Task Group on Scenarios for Climate Impact Assessment guideline (Carter et al. 2007), the 30-yr periods were considered for the evaluation of changes between the end of century (2071–2100) and historical (1976–2005) periods. Further, following a common practice (e.g., Jiménez Cisneros et al. 2014), the 2071–2100 period is referred to as 2080s. Compared to the historical period, the 2080s seasonal temperature and precipitation mostly show increases (Fig. S4, Table S2), which is consistent across GCMs and RCPs. Between the GCMs, HadGEM2-ES and GFDL-ESM2M tend to have the highest and lowest temperature increases, respectively. In the case of precipitation, the increases are mostly consistent among GCMs, RCPs, and seasons, with higher increases in OND and lower increases in JFM and JAS. However, there is a large variability in future changes for AMJ and OND, with CanESM2 RCP8.5 scenario increasing by +50%, and MPI-MR RCP4.5 and GLDL RCP4.5 showing no change for AMJ and OND, respectively.
c. Snow water equivalent changes
Given the potential impacts of temperature and precipitation increases on snow storage, we considered changes in monthly and annual SWE, as well as the magnitude and timing of SWEmax (Figs. 3, 4). Results are based on 30-yr means of each GCM, with the maximum, median, and minimum values from seven GCM ensemble members depicted. The comparison of monthly and annual SWE responses (2080s versus 1976–2005) revealed that despite the increase in mean annual and seasonal precipitation (Table S2), the mean annual SWE declines for all subbasins (Fig. 3). As expected, a larger decline in SWE occurs for the RCP8.5 scenario, corresponding to the larger decline in snowfall (Fig. 2d). Thus, under RCP8.5, disproportionately larger precipitation increases would be required to maintain a particular SWE level. In the case of mean monthly SWE, declines are projected for the Toad-NC subbasin and Liard-M basin, while some of the winter months could experience increases (especially for RCP4.5) for the Liard-UC and S-Nahanni-VF subbasins. The SWE values for the summer months are projected to decline for all subbasins. The Toad-NC subbasin becomes snow-free from July to September for the majority of RCP8.5 ensemble members. The summer snow storage loss also leads to considerably smaller SWE values from October to June for both RCPs compared to the historical period (Fig. 3a).
In contrast to these generally consistent declines in the monthly and annual SWE values, SWEmax values could either increase or decrease depending on the subbasin and RCP considered (Fig. 4a). Specifically, while there are declines for both future scenarios for the Toad-NC and Liard-M basins, there are increases for the majority of RCP4.5 and RCP8.5 ensemble members for the S-Nahanni-VF. In the case of Liard-UC, the future SWEmax values increase for the majority of RCP4.5 ensemble members and decreases for the majority of RCP8.5 ensemble members. A possible explanation for the differences is the north–south temperature gradient, with the warmer southern Toad-NC experiencing declines and colder northern S-Nahanni-VF experiencing increases. In the case of SWEmax timing, which is an indicator of snowmelt initiation, progressively earlier occurrences are projected under future warming scenarios. A closer examination of the RCP4.5 and RCP8.5 results revealed a generally consistent pattern of earlier timing for smaller SWEmax and later timing for larger SWEmax (not shown). This relationship between the magnitude and timing of SWEmax indicates the importance of temperature and precipitation interactions on snowpack accumulation and melt. Overall, the entire basin is projected to experience consistent declines in the monthly, annual, and maximum SWE in 2080s compared to 1976–2005. For instance, for the Liard-M basin, there are pronounced decreases in the projected SWEmean and SWEmax, and earlier SWEmax timings (Table S3).
d. Streamflow changes
The projected increases in monthly and annual mean flows by the end of century (2080s) compared to the baseline period for the subbasins of the LRB (Fig. 5) are consistent with the higher precipitation input to the basin (Fig. 2 and Fig. S4). Although the hydrologic regime of the basin is projected to remain nival, as the signature pattern of rapid rise in spring/summer flow following the snowmelt persists under both RCPs, the seasonality of future flow is projected to undergo prominent changes. For instance, enhanced flows for the fall–winter months are projected, attributable to the combined effect of the warmer and wetter future climate and consequent increases in rainfall fraction. This pattern is generally consistent, except for the northernmost and coldest (Table S2) subbasin of S-Nahanni-VF (Fig. 5c), likely due to the persistence of winter subfreezing conditions that inhibit rainfall events and snowmelt. On the other hand, future responses for the spring–summer months are similar across the LRB. Specifically, the months of April–May exhibit substantially higher flows from 1976 to 2005 to 2080s for all subbasins, a pattern that is consistent with the rapid depletion of snow storage for these months and RCPs (Fig. 4). For both winter and summer, flow increases are higher for RCP8.5 compared to RCP4.5. Flows in the month of June show continued increase for RCP4.5, and decline for RCP8.5 in response to an earlier depletion of snow storage. Continued declines in flows are projected for the months of July and August despite increased precipitation in these months (Fig. S4, Table S2), attributable to the earlier snowmelt-driven flow recession and higher rates of evapotranspiration loss. As a result, the seasonal flow ratios (expressed as October–March to April–September flow ratios) change progressively with the warmer climates (Table S4).
The effects of temperature and precipitation changes are also evident in streamflow extremes. Annual maximum flow Qmax, which occurs primarily as a result of snowmelt, is projected to increase in 2080s for most GCM ensemble members, and decline for some members (Fig. 6a). It is also interesting that while both Qmax and SWEmax increase for the majority of GCM ensemble members for the Liard-UC and S-Nahanni-VF subbasins, Qmax for RCP4.5 increase despite the decrease in SWEmax for the Toad-NC and Liard-M basins (Fig. 4a). Thus, there is no straightforward relationship between snow storage and maximum flow. This also suggests that mechanisms other than snowpack storage, such as rainfall and possibly rain-on-snow, play increasingly important roles in the generation of future Qmax. The increased contribution of rainfall could be a major factor for RCP8.5, as the increases in Qmax occur despite the greater snowpack declines (Fig. 3).
On the other hand, a generally consistent pattern of progressively earlier Qmax timing emerges with greater warming scenarios. Specifically, compared to the historical period, the median timings are projected to shift by about 10 days under RCP4.5 and a further 10 days under RCP8.5 (Fig. 6b). Further, Qmax dates lag SWEmax dates by more than 2 months (Fig. 4b). This is because the SWEmax date is an indicator of snowmelt initiation, while Qmax usually occurs after a depletion of the large fraction of snow storage, and is dependent on a number of factors, such as snowpack depth, melt rates, rainfall, and travel time for the meltwater from different parts of the basin to reach the basin outlet.
Minimum flow Qmin in the region occurs during the winter snow/river-ice-covered season, when subfreezing temperature inhibits rainfall events and prevents snowmelt, and baseflow from groundwater sustains low-flow discharge with the flow amount reflecting prewinter groundwater storage status (Woo and Thorne 2014). The results depict a consistent pattern of increase in winter low flows for all subbasins of the LRB, attributable to the increased fraction of winter rainfall and the subsurface flow generation under the warmer future climate. Consequently, there is a larger increase in low flow for the warmer RCP8.5 forcing. Overall, generally consistent patterns emerge from the future projections of monthly, annual and minimum hydrologic indicators. These include pronounced increases in Qmean and Qmin. In the case of Qmax, while the 2080s values could either increase or decrease depending on the choice of GCMs, their timings are projected to occur earlier for all GCMs (Table S3).
e. Variable importance of hydrologic projections
The VI percentage scores obtained from the RF models depict the relative controls of seasonal climate variables on hydrologic indicators (Fig. 7). The predictability scores from the RF model are given as NSE values. The sensitivity plots of these indicators with respect to temperature and precipitation changes (Fig. 8) are based on VI analyses, with multiple seasonal temperature and precipitation combined when any of the seasonal VI is higher than 10%.
The RF model results for SWEmax illustrate the dominance of OND_P and JFM_P with VI > 30%, while temperatures for these two seasons provide moderate influences (VI > 10%). The relationship between the variables is reinforced by a good fit (NSE = 0.76) of the RF model, hence there is a good predictability of SWEmax using seasonal temperature and precipitation. Figure 8a, which depicts the sensitivity of SWEmax in relationship to OND + JFM temperature and precipitation changes, illustrates the trade-off between increasing SWEmax with increasing precipitation and decreasing SWEmax with increasing temperature. Hence, a higher precipitation increase would be required to maintain a certain level of SWEmax as the basin gets progressively warmer. Under the RCP8.5 scenario, even the large winter precipitation increase is not able to offset the temperature driven decline in SWEmax, and the occurrences of mean historical SWEmax value (~210 mm) become rare. Overall, while the winter precipitation amount controls the SWEmax magnitude, winter temperature also exerts an important influence.
In contrast to SWEmax magnitude, winter temperatures have a larger influence on SWEmax timing. This is indicated by VI > 30% for JFM_T and AMJ_T, the two seasons when SWEmax typically occurs, although JFM_P also show some influence (VI > 10%). However, the RF model fit for this indicator is relatively low (NSE = 0.41), so the predictability of SWEmax timing is limited using seasonal temperature and precipitation. Further, although the sensitivity plot of SWEmax timing in relationship to JFM + AMJ temperature and precipitation changes show some evidence of temperature dependence (e.g., progressively earlier SWEmax with increasing warming Fig. 8b), the relationship is not clearly stratified, implying that seasonal precipitation and temperature have limited capability to explain the variability of SWEmax timing.
The VI for annual flow are >15% for precipitation in all four seasons, indicating that Qmean is predominantly controlled by precipitation. It is interesting to note that the JAS months have the highest precipitation amount (Table S2) but relatively low VI. This could be explained in terms of higher evapotranspiration losses in summer months and thus lower net water input to the basin. However, on an annual basis increased wetness tends to overshadow the effect of increased evaporative demand due to warming, which leads to increased annual flows. The model performance for Qmean is good (NSE = 0.78), hence, it is possible to provide good estimates of annual flows using seasonal temperature and precipitation. The precipitation dependence of Qmean can also be seen when the historical and projected future flows are considered with respect to annual temperature and precipitation changes (Fig. 8c). In this case, a clear stratified pattern emerges depicting higher flows with increasing wetness and vice versa.
The VI score greater than 20% for OND_P, JFM_P, and AMJ_P indicates that Qmax is primarily controlled by precipitation, while seasonal temperatures, with VI > 7%, exerts some influence. Thus, while the increased precipitation generally leads to increased Qmax, this could be offset by increased temperature and consequently reduced snowpack volume. The predictability of Qmax using seasonal temperature and precipitation is moderate (NSE = 0.50). Furthermore, while the sensitivity plot of Qmax in relation to OND + JFM + AMJ temperature and precipitation (Fig. 8d) show some evidence of precipitation control (e.g., increase in Qmax with precipitation increase), the lack of clear stratification is indicative of more complex interactions.
There is relatively high influence of temperature on Qmax timing, with VI > 14% for the months when the Qmax events occur (JFM_T and AMJ_T). Thus, similar to SWEmax timing, Qmax timing is mainly temperature controlled. However, the low RF model fit score (NSE = 0.35) suggests that seasonal temperature and precipitation do not sufficiently explain the variability of Qmax timing. This is further illustrated by the sensitivity plots of Qmax timing in relation to AMJ + JAS precipitation and temperature (Fig. 8e), where the pattern of early Qmax with increasing temperature is marked by the lack of clear stratification.
The primary controls for Qmin are JFM_T and OND_T, with the relative VI of 42% and 19%, respectively. The low flow in the region is mainly contributed by baseflow from groundwater (Woo and Thorne 2014) with an expectation of increased groundwater recharge and subsurface flow as temperature increases (Bense et al. 2009). This is especially the case for RCP8.5 that exhibits large increases in Qmin as the basin gets increasingly warmer (Fig. 6c). However, there is only a moderate predictability (NSE = 0.54) of Qmin using seasonal climate variables as predictors. This is further illustrated by the moderately stratified pattern of increasing Qmin with increasing warming in the sensitivity plot in relation to OND + JFM temperature and precipitation (Fig. 8f).
Overall, the VI analysis provides an effective means of evaluating climatic controls on hydrologic changes. Although, not a focus of this study, the results could be used to develop a computationally efficient statistical modeling framework, for example, for emulating hydrologic responses as a function of driving covariates (e.g., Luce et al. 2014; Schnorbus and Cannon 2014; Shrestha et al. 2017).
f. Global mean temperature, climatic controls, and hydrologic changes
With these insights from the VI analysis, we considered future changes in climatic controls and hydrologic indicators with respect to the policy relevant global mean temperature (GMT) changes. Figure 9 depicts the variability of hydrologic responses under the 1.0°, 1.5°, 2.0°, 2.5°, and 3.0°C GMT changes with respect to the preindustrial (PI) period of 1850–2000 following IPCC (2018). Accordingly, 1.0°C global warming corresponds to a period centered on 2017. It is to be noted that the periods of a GMT change differ between GCMs because each GCM has a different level of climate sensitivity and responds differently to the same radiative forcing scenario (e.g., Schleussner et al. 2016). In this study, we used the current/future 31-yr periods of GMT thresholds from Jeong et al. (2019), in which the periods for each GCM were calculated relative to 1986–2005 by assuming a 0.61°C GMT increase from 1850 to 1900 following IPCC (2013) (see Table S5). Some of the RCP4.5 GCMs do not reach more than 1.5°C GMT change, so differences between RCP4.5 and RCP8.5 responses could not be compared.
Figure 10 illustrates interactions of the seasonal/annual climatic controls with the hydrologic indicators under different GMT changes. Given the lack of observation data for the preindustrial period in the study region, the Liard basin temperature and precipitation changes were considered with respect to the common baseline (1976–2005) period. Thus, the temperature changes shown do not include an estimated 1.4°C warming from 1945 to 1976–2005 (based on temperature trend of +0.3°C decade−1, see section 4b) plus unknown warming from 1850–1900 to 1945. Nevertheless, as expected for the subarctic region, the basin-scale temperature increases are generally higher than the GMT changes. The results are summarized for the 31-yr means, and as in Fig. 8, multiple seasonal temperature and precipitation for the Liard basin were combined when any of the seasonal VI is higher than 10%.
Overall, the hydrologic responses tend to get more pronounced under higher GMT changes, while the basin-scale temperature and precipitation provide controls consistent with the results of VI analysis. Specifically, while SWEmax shows a general tendency of declines under higher GMT changes (Fig. 9a), the declines are partially moderated by the OND + JFM precipitation increases (Fig. 10a). In the case of SWEmax timing, there is a tendency to earlier occurrences under higher GMT changes (Fig. 9b), while the influence of basin-scale warming is evident in terms of earlier timing with higher temperature increases (Fig. 10b).
The ranges of Qmean tend to get wider under higher GMT changes, with the values for some GCMs exceeding those under 1.0°C (Fig. 9c). Consistent with the results of VI analysis, the Qmean increases are predominantly driven by basin-scale precipitation increases, with the higher precipitation increases under higher GMT changes leading to larger Qmean values (Fig. 10c). The patterns of Qmax responses are similar to Qmean, with some GCMs projecting larger extreme events under higher GMT changes (Fig. 9d). Such increases in Qmax occur when precipitation increases are high (e.g., under 2.5° and 3.0°C GMT changes), thus emphasizing the role of precipitation control (Fig. 10d). The Qmax timings respond similarly to SWEmax timings, with earlier shifts under higher GMT changes (Fig. 9e). The shifts are primarily influenced by temperature increases, while the precipitation declines seem to reinforce the shifts with the earliest Qmax occurring when precipitation decreases (Fig. 10e). Likewise, Qmin mostly increases with each 0.5°C GMT change (Fig. 9f). This signal in general, is reinforced by the basin-scale temperature increase, with higher the warming, larger the minimum flow (Fig. 10f).
Hence, in addition to the end of century projections, the results provide perspectives on the policy relevant GMT changes. Overall, under each 0.5°C global warming above the preindustrial period, the basin is projected to experience progressively enhanced warming and wetness, and more pronounced hydrologic changes. For instance, there are large differences in the ranges of Qmax, Qmax timing, and Qmin between 1.5° and 2.0°C GMT changes. Thus, differences in the severity of impacts, such as the timing and magnitude of spring flood, could be considerable if the GMT increase were to be limited to 1.5° versus 2.0°C. Comparing the projections under different GMT changes with the 2080s RCP8.5 scenario (Figs. 3–6), although the directions of changes are generally similar, the magnitudes are generally smaller. This is because the basin-scale warming under 3°C GMT change (Fig. 10) is less than RCP8.5 2080s ensemble (Table S2). On the other hand, hydrologic changes under a 3°C GMT increase generally match those of the RCP4.5 2080s ensemble due to the similar level of warming and precipitation increases.
5. Summary and conclusions
This study provides a state-of-the-art assessment of the hydrologic impacts of climate change due to enhanced warming and amplified moisture fluxes in the subarctic, with an application for the Liard River basin, western Canada. To this end, we employed the 1/16° grid resolution VIC hydrologic model, with the ensemble of CMIP5 climate forcings downscaled by using a multivariate method that explicitly takes into account the intervariable dependence. We evaluated the hydrologic projections by using the Random Forests model and identified the dominant climatic controls on a set of hydrologic indicators. Further, we analyzed the interactions of climatic controls with hydrologic indicators under 1.0°–3.0°C GMT changes. The lack of accounting of frozen soil/permafrost processes is the main limitation of this study and a development of efficient methods for simulating these processes need to be prioritized for future studies.
The projected changes for the end of century (2080s) in comparison to baseline period (1976–2005) indicate a substantially warmer and wetter climate, with higher temperature and precipitation increases for RCP8.5 GCM ensemble compared to RCP4.5. The projections for the hydrologic indicators reveal a range of changes, with more pronounced changes for the warmer RCP8.5 scenario. Specifically, despite the substantial increases in precipitation, snowfall is projected to decline leading to reductions in the annual and maximum SWE together with shifts to earlier maximum SWE. On the other hand, increased annual flows with larger contributions of rainfall-generated runoff are projected. Further, there are pronounced increases in the projected minimum flows, and earlier onsets of maximum flows. Results for maximum flow magnitudes are less consistent, with either increases or decreases depending on interactions of temperature and precipitation in GCM projections. The direction of changes under the GMT changes of 1°–3°C changes are similar to the 2080s projections, while the magnitudes under a 3°C GMT change generally match the RCP4.5 2080s ensemble.
The variable importance model captured the sensitivity of the hydrologic indicators to climatic variables and helped to provide physically plausible explanations for future hydrologic changes. Using the model, predictability is good for maximum SWE and mean flow, moderate for maximum and minimum flow, and low for timings of maximum SWE and flow. In terms of sensitivities, winter–spring precipitation is the primary control for maximum SWE, while January–June temperature has a greater influence in its timing. Consequently, the changes in maximum SWE reflect a trade-off between precipitation driven snow accumulation and temperature driven snowmelt initiation. On the other hand, the mean and maximum flows are also predominantly precipitation controlled with increasing wetness generally leading to higher discharge. Seasonal temperatures exert greater influences on maximum flow timing and minimum flows, apparently by driving snowmelt and contributing to increased subsurface runoff, respectively, as temperature increases. These hydrologic responses tend to get more pronounced under higher GMT changes, with a general pattern of smaller maximum SWE, earlier maximum SWE and maximum flow, and larger minimum flows for most GCMs; and larger mean and maximum flows for some GCMs. The basin-scale temperature and precipitation controls reinforce these signals, which are consistent with the results of VI analysis.
Hence, the framework presented this study goes beyond projecting hydrologic impacts of climate change. This study for a subarctic river basin is highly relevant because the region is experiencing enhanced warming and substantial precipitation increases, and detailed basin-scale hydrologic modeling studies for the region are scarce. The detailed diagnosis of the hydrologic changes, as well as controls and interactions of climatic variables on key hydrologic indicators, could be generalized for understanding regional scale changes in subarctic/nival basins. Additionally, by evaluating the hydrologic changes with respect to global mean temperature changes, a framework for understanding the policy relevant global versus regional warming was developed. Further, the method developed in this study could be used as a diagnostic tool for understanding future changes from different generations of GCMs and emissions scenarios, as well as a guidance to set up/improve computationally efficient statistical modeling approaches for simulating hydrologic changes.
We thank two anonymous reviewers for their comments, which led to an improved version of the manuscript. This paper is Canadian Crown Copyright.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JHM-D-18-0262.s1.