The rapidly warming Arctic is experiencing permafrost degradation and shrub expansion. Future climate projections show a clear increase in mean annual temperature and increasing precipitation in the Arctic; however, the impact of these changes on hydrological cycling in Arctic headwater basins is poorly understood. This study investigates the impact of climate change, as represented by simulations using a high-resolution atmospheric model under a pseudo-global-warming configuration, and projected changes in vegetation, using a spatially distributed and physically based Arctic hydrological model, on a small headwater basin at the tundra–taiga transition in northwestern Canada. Climate projections under the RCP8.5 emission scenario show a 6.1°C warming, a 38% increase in annual precipitation, and a 19 W m−2 increase in all-wave annual irradiance over the twenty-first century. Hydrological modeling results suggest a shift in hydrological processes with maximum peak snow accumulation increasing by 70%, snow-cover duration shortening by 26 days, active layer deepening by 0.25 m, evapotranspiration increasing by 18%, and sublimation decreasing by 9%. This results in an intensification of the hydrological regime by doubling discharge volume, a 130% increase in spring runoff, and earlier and larger peak streamflow. Most hydrological changes were found to be driven by climate change; however, increasing vegetation cover and density reduced blowing snow redistribution and sublimation, and increased evaporation from intercepted rainfall. This study provides the first detailed investigation of projected changes in climate and vegetation on the hydrology of an Arctic headwater basin, and so it is expected to help inform larger-scale climate impact studies in the Arctic.
Recent changes in the Arctic region climate (Wanishsakpong et al. 2016; Whitfield et al. 2004), vegetation (Xu et al. 2013), and other environmental functions (Hinzman et al. 2005) motivate investigation of the future hydrology of the Arctic. Changes in streamflow discharge (Mendoza et al. 2015; Arheimer and Lindström 2015; Gelfan et al. 2017), permafrost thaw (Woo et al. 2007), subsurface water storage and flow (Walvoord et al. 2012), and snow accumulation and cover (Liston and Hiemstra 2011) are of great interest under scenarios of changing climate. An improved understanding of hydrological change is necessary to more effectively adapt and mitigate the potential impacts of climate change; however, given the complexity of the environment and the uncertainty associated with climate projection, this has represented a great scientific challenge.
Climate projections from global circulation models (GCMs) agree on a warmer and wetter Arctic by the end of the century (Kattsov et al. 2005). Because of problems in GCM representations of regional or local surface weather, higher-resolution (tens of kilometers) regional climate models (RCMs) are used to dynamically downscale GCMs; however, large-scale RCMs still fail to represent surface weather, particularly precipitation in areas with complex topography, deep convection, or extreme events (Prein et al. 2015), requiring further downscaling. Downscaling techniques are classified into statistical and dynamical approaches, the former using empirical relationships between observed and simulated climate and the latter requiring the implementation of a high-resolution climate model (Maraun et al. 2010). Dynamical downscaling produces physically connected weather variables (Fowler et al. 2007), which is critical for cold regions (Pomeroy et al. 2015a), as the lack of physical realism in relationships among driving meteorological variables restricts the implementation of statistical downscaling in hydrological studies using physically based cold-region hydrological models. Recent studies have argued that dynamical downscaling using convection permitting models (CPMs; spatial resolution < 4 km) is required to properly represent changes in extreme precipitation events (Kendon et al. 2017; Prein et al. 2015) that are hydrologically important.
Shrub expansion and densification have been well documented in the Arctic (Lantz et al. 2013; Myers-Smith et al. 2011). However, there are no published studies investigating changes in the forest structure (i.e., density, height, and extension) over northwestern Canada. Payette and Filion (1985) found that white spruce tree lines in northern Quebec have not substantially changes over the past centuries, whereas Suarez et al. (1999) found that the tundra–taiga tree line in Alaska advanced northward between 80 and 100 m north over the last 200 years. Gamache and Payette (2004) studied black spruce height near the Arctic tree line in eastern Canada and found that height growth has not significantly changed. Trends in greening and browning have been studied in Canada and Alaska using remote sensing, showing an spatially heterogeneous response but a clear greening trend in northwestern Canada (Ju and Masek 2016), which is likely driven by reported shrub expansion and densification. Zhang et al. (2013) investigated Arctic vegetation projections under future climate conditions and showed an overall shrubification of the tundra; however, virtually no change in the tundra–taiga transition in the northwestern Canadian Arctic was found.
Arctic hydrological processes needed to calculate basin hydrology below or at the tree line include snow accumulation and melt; sublimation and unloading of intercepted snow from forest canopy; blowing snow redistribution and sublimation; evapotranspiration; infiltration into frozen and unfrozen soils; ground freeze and thaw; water flow through snowpack, surface, and subsurface flow; and groundwater and streamflow routing (Kane et al. 1991; Pomeroy et al. 2008). These processes have been included in the spatially distributed and physically based Arctic Hydrology Model (AHM) developed and verified by Krogh et al. (2017) using the Cold Regions Hydrological Modeling (CRHM) platform (Pomeroy et al. 2007).
The purpose of this study is to investigate the effect of future climate and vegetation changes on the hydrological processes of a small headwater Arctic tree-line basin underlain by continuous permafrost. These steps were followed to pursue this goal: 1) climate projections from a high-resolution climate model under a convection-permitting configuration were compared to surface observations and used to force the CRHM-AHM under historical and future conditions; 2) vegetation projections based on observed rates of changes were used to parameterize land cover for the CRHM-AHM; 3) climate and hydrological projections were analyzed and discussed; and 4) a sensitivity analysis examined the impact of vegetation change on the basin water balance.
2. Study site
Havikpak Creek basin (HPC; Fig. 1), located in the Northwest Territories, Canada, was selected as it has a history of hydrological process studies and hydrological modeling applications and is located near the tundra–taiga transition, where changes in vegetation are anticipated and may impact the hydrology. HPC is a small headwaters basin (16.4 km2) that is underlain by continuous permafrost and covered primarily by taiga forest (>50%), with large areas of open tundra, shrubs, wetlands, and open water. HPC mean annual precipitation and temperature from 1980 to 2010 are 327 mm and −8.2°C, respectively (Krogh et al. 2017), resulting in relatively dry conditions and long winters (October–April). Soils at HPC are characterized by an upper layer of permeable organic peat, composed of decomposed vegetation, lichen, and moss, followed by a lower peat layer over a relatively impermeable mineral soil layer. A detailed description of HPC, including soil characteristics, meteorology, and other characteristics, are presented and discussed by Krogh et al. (2017).
a. Automated weather stations
An hourly and long-term meteorological time series for Havikpak Creek was reconstructed by Krogh and Pomeroy (2018) for the period 1960–2016, including precipitation, temperature, wind speed, relative humidity, and shortwave and longwave irradiance, based on a combination of in situ meteorological observations, the Adjusted and Homogenized Canadian Climate Dataset (AHCCD; Mekis and Vincent 2011) at the Inuvik station (ID 2202578), and the ERA-Interim (Dee et al. 2011) and ERA-40 atmospheric reanalyses (Uppala et al. 2005). The automated weather stations used were the Inuvik Climate, the Inuvik Airport, and the Inuvik Upper Air (Fig. 1), all maintained by Environment and Climate Change Canada (ECCC). For details about the time series reconstruction, readers are referred to Krogh and Pomeroy (2018).
Daily streamflow records at the Havikpak Creek station (ID 10LC017) have been measured by the ECCC Water Survey of Canada (WSC) since 1995. The hydrometric station is located downstream from the Havikpak Creek crossing with the Dempster Highway (Fig. 1). Arctic stream gauging is challenging, particularly in small creeks due to the presence of ice and snow in the cross section during the spring snowmelt runoff, in which the annual peak streamflow discharge and the majority of the annual discharge volume typically occurs. These problems are acknowledged in the metadata provided by ECCC through the Environment Canada Data Explorer.
c. Atmospheric model: The WRF Model
1) Historical simulations
The Weather Research and Forecasting (WRF; Skamarock et al. 2008) Model is a state-of-the-art numerical weather prediction (NWP) and atmospheric modeling system developed by a group of U.S. government agencies led by the National Center for Atmospheric Research (NCAR). WRF (version 3.4.1) was run at a convection-permitting resolution of 4 km over western Canada for the period 2000–13 (Li et al. 2016). The initial and lateral boundary conditions used were 6-hourly time series from the ERA-Interim reanalysis at a 0.7° spatial resolution, and surface states and fluxes from the Noah-MP (Niu et al. 2011) land surface scheme. Main benefits of these runs are 1) the large extent of the spatial domain, 2) the decadal period, 3) the high-resolution topographic representation, and 4) the convection-permitting configuration. The last two have been shown to greatly improve summer and winter precipitation representation, as opposed to climate models using cumulus parameterizations with lower spatial resolution (Brisson et al. 2016; Fosser et al. 2015; Prein et al. 2015; Rasmussen et al. 2011, 2014), providing more robust precipitation projections under future climate scenarios (Kendon et al. 2014, 2017). The outputs from this run used in this study are 2D hourly time series at the ground surface of precipitation, air temperature, wind speed, specific humidity, and shortwave and longwave irradiance.
2) Future weather simulation
The pseudo-global-warming (PGW) approach (Schär et al. 1996) was used to produce future weather simulations. The PGW approach adds a mean monthly perturbation to the ERA-Interim initial and lateral boundary conditions to the period 2000–13, which are calculated as the difference between the 25-yr monthly values between the 1975–99 and 2075–99 periods from the CMIP5 model intercomparison experiment (Taylor et al. 2012) under the RCP8.5 gas concentration scenario (Riahi et al. 2011), resulting in a time series associated with the 2086–99 period. Two benefits of the PGW approach are the reduction of uncertainty caused by the interannual variability and the reduction of model bias contained in the GCM projections (Kawase et al. 2008), and the main disadvantage is that it does not allow for future interannual variability. Previous studies have used PGW to quantify hydrological changes under future climate scenarios (Ma et al. 2010; Mendoza et al. 2016).
Figure 2 presents the modeling flowchart that summarizes the methodology used in this study. The first modeling stage comprises historical hydrological simulations using bias-corrected simulated weather from WRF. Historical runs are validated against daily observed streamflow discharge. Future simulations use the bias-corrected simulated weather from WRF-PGW with vegetation projections based on observed rates of growth. Finally, the sensitivity of the future mass balance to projected changes in vegetation is performed.
a. WRF bias correction
WRF weather was bias corrected to generate forcing data that are as representative as possible of observed records. This study uses the univariate quantile mapping correction from the R package MBC (https://CRAN.R-project.org/package=MBC), which uses the quantile delta-mapping algorithm described by Cannon et al. (2015). Quantile mapping correction was performed to hourly precipitation, air temperature, wind speed, and water vapor pressure, using the observed weather time series presented by Krogh and Pomeroy (2018). Air temperature was divided into two periods: 1) spring (April–May) and 2) nonspring (June–March) for quantile mapping correction in order to compensate for the larger WRF cold bias that was during spring. Precipitation was also divided into two periods: 1) winter and 2) summer using the 0°C mean temperature threshold, which was important to properly represent cumulative snowfall and extreme rainfall events. Wind speed correction was performed to the full period. Water vapor pressure correction was required to properly calculate relative humidity following the Clausius–Clapeyron relationship. Quantile correction was performed to the entire period, and then relative humidity was calculated based on corrected air temperature, water vapor pressure, and the Buck formula (Buck 1981) to compute the saturated water vapor with respect to water and ice. Relative humidity was not allowed to exceed 100%.
b. Hydrological model
The hydrological model used in this study is the AHM (Krogh et al. 2017) developed with the CRHM platform (Pomeroy et al. 2007). The CRHM-AHM is a physically based and spatially distributed hydrological model that includes the following key Arctic physical processes: blowing snow redistribution and sublimation, snowmelt energy balance, sublimation/evaporation of canopy intercepted snowfall/rainfall, soil moisture storage and flow, evapotranspiration, infiltration into frozen and unfrozen ground, flow through organic terrain and snowpack, ground freeze and thaw, surface runoff, and streamflow routing. The model uses hydrological response units (HRUs; Flügel 1995) to spatially discretize HPC based on land cover classes, elevation, and topographic features such as gullies and snowdrifts. This discretization resulted in 11 HRUs: upper and lower tundra, upper and lower sparse shrubs, upper and lower gully/drift, close shrubs, taiga forest, forest, wetland, and open water. CRHM-AHM is run forced by hourly precipitation, air temperature, wind speed, relative humidity, and shortwave and longwave irradiance. Details about model parameterization, calibration of few subsurface and surface hydraulic and storage parameters, validations against daily streamflow, snow accumulation and active layer thickness, and sensitivity analysis to key parameters are presented by Krogh et al. (2017). An updated streamflow validation is presented by Krogh and Pomeroy (2018).
1) Historical modeling
The CRHM-AHM was run for the period January 2001–October 2013 using the corrected weather time series from the four closest WRF centroids to HPC (Fig. 1), for which centroid 1 provided the best simulation, and therefore, it was used for all the analyses. Note that there were small differences among weather from the different centroids, resulting in very small differences between the four streamflow simulations, and, therefore, the impact of the centroid selection is expected to be minimal for the analysis and discussion presented in this study.
Shrub extension and density characteristics were taken as the average for the years 2001–13 based on the extrapolated rates presented by Krogh and Pomeroy (2018), built on observations by Pomeroy and Marsh (1997) and observed shrub changes from Lantz et al. (2013). The vegetation characteristics used for the historical modeling are presented in Table 1. Subsurface and surface hydraulic and storage parameters calibrated by Krogh et al. (2017) were used for this simulation; nevertheless, the uncertainty introduced by this new weather time series was investigated by performing another automated calibration with the dynamically; dimensioned search algorithm (DDS; Tolson and Shoemaker 2007).
2) Future modeling
Future hydrological modeling includes both climate and vegetation projections. Climate projections are those from the bias-corrected WRF-PGW for the equivalent period of 2087–99. The year 2086 was used as a spinup period for the climate model. Vegetation projections assume that the observed rates of changing shrub cover and density from Lantz et al. (2013) remain in the future, as no other projections are available. To include the “new” sparse shrubs in the CRHM-AHM, two new HRUs were added: the upper and lower “new” sparse shrubs, resulting in a model with 13 HRUs. Projected shrub area, stem density, and leaf area index (LAI) for the projected HRUs are presented in Table 1. Vegetation height of the new sparse shrubs was estimated to be 0.8 m, which corresponds to roughly half of the estimated average height of current sparse shrubs at HPC (Krogh et al. 2017) and reflects that these plants are colonizing previous tundra-covered surfaces. There are no quantitative projections of forest height, density, or extension available, and, therefore, the forested HRUs were held constant in the future model configuration.
c. Mean change analysis
Projected mean annual changes between historical and future simulations for several water fluxes and state variables were tested using the nonparametric Wilcoxon rank test (Wilcoxon 1945) with a significance threshold of p ≤ 0.05. The Wilcoxon test assumes that the two series are independent and continuous.
a. WRF validation
Figure 3a presents a comparison of annual cumulative precipitation between observed, raw, and bias-corrected WRF historical data. Raw WRF mean bias was only 8 mm yr−1; nevertheless, this bias was removed after applying the quantile mapping correction of daily precipitation, which also improved the representation of high-precipitation events (Fig. 3b) that can produce important rainfall–runoff responses. Figure 3c presents the Q–Q plot between observed and bias-corrected WRF precipitation, in which the high-precipitation events are well represented. Despite the good representation of mean precipitation and its quantiles, bias-corrected WRF has mixed performance representing the variability of annual precipitation. Figure 3d compares monthly observed, raw, and bias-corrected WRF air temperature. An overall cold bias of 2.7°C was found for the raw WRF; however, seasonal temperatures and quantiles (Fig. 3e) were well represented by raw WRF. After bias correction of hourly temperature, WRF air temperatures agreed well with observations and the cold bias was removed (Fig. 3f). Ten-meter hourly wind speed (not shown) simulated by WRF overestimated surface observations by an average of 1.5 m s−1; this was removed after bias correction. Calculated relative humidity from WRF bias-corrected water vapor pressure and air temperature represented observed mean and quantiles well (not shown).
Table 2 presents a comparison between observed and bias-corrected WRF daily precipitation for wet and dry spells. A wet (dry) spell is defined as any period with at least three consecutive days with precipitation above (below) 0.1 mm day−1. Observed mean annual number of dry and wet spells is 28 and 14, respectively, whereas simulated mean annual number of dry and wet spells is 30 and 13, respectively. The mean annual length of dry spells is 8.9 and 8.1 days for observation and simulations, respectively, whereas the mean length of wet spells is 5.4 and 5.5 days for observations and simulations, respectively.
b. Projected changes in climate
Figure 4 presents mean daily time series of projected changes using bias-corrected WRF simulations of air temperature, cumulative precipitation, wind speed, and relative humidity, and raw WRF simulations of shortwave and longwave irradiance; mean annual changes are also included in bold when they are statistically significant (p ≤ 0.05). Table 3 presents seasonal and annual changes for the same variables presented in Fig. 4 and the all-wave irradiance. A significant warming of 6.1°C in mean annual air temperature is projected by the WRF-PGW, particularly during winter and spring, for which an increase of 6.8° and 7.1°C, respectively, is indicated. The date at which temperature reaches 0°C is projected to occur roughly 2 weeks earlier in spring, whereas the date at which temperature drops below freezing is delayed by 16 days in fall. Cumulative precipitation significantly increased by 117 mm yr−1 or 38% with respect to the historical period, with the largest seasonal increase in winter (October–April, 94.3 mm). A negligible increase in mean annual wind speed of about 0.1 m s−1 is projected. The simulated mean annual increase in relative humidity is about 4%, with the largest increases in winter of about 10%. Shortwave irradiance decreased annually by 2.1 W m−2; however, a small increase of roughly 2.5 W m−2 was found in summer, whereas longwave irradiance significantly increased by 21.2 W m−2 annually, with the largest seasonal increase in winter of 22.8 W m−2, which is partially explained by the increase in air temperature and specific humidity and decrease in shortwave irradiance and hence atmospheric transmittance. The total change in all-wave irradiance is an annual increase of 19 W m−2, which is consistent throughout the year except in spring, during which a small decrease of 3.4 W m−2 was projected.
c. Hydrological model performance
The CRHM-AHM performance using the parameters from Krogh et al. (2017) and the corrected WRF meteorology resulted in a Nash–Sutcliffe efficiency (NSE) and mean bias of 0.44 and −16%, respectively, for the period 2002–12 using daily streamflow. These results are consistent with those presented by Krogh et al. (2017) and Krogh and Pomeroy (2018), suggesting a robust representation of Havikpak Creek hydrology. Nevertheless, to assess the sensitivity of the model to this new forcing data, the few surface and subsurface storage and flow parameters calibrated by Krogh et al. (2017) were recalibrated using corrected WRF data and the NSE as objective function. Five-hundred model iterations using DDS for the period 2002–08 were performed, and 2009–12 for validation. The result of the parameter recalibration and the one performed by Krogh et al. is presented in the appendix (Table A1). NSE and mean bias for the entire 2002–12 period using corrected WRF were 0.45 and −18%, respectively, suggesting a marginal improvement in NSE and a deterioration of mean bias. Parameters from the recalibration show substantial differences with some of those presented by Krogh et al. (2017); this is partially explained by the conceptual nature of some of the few calibrated parameters (e.g., surface and subsurface lag and routing), which does not guarantee a consistency between calibration experiments.
Figure 5 shows a streamflow comparison between observed (blue), simulated using parameters from Krogh et al. (2017) (red), and simulated using recalibrated parameters (green). Annual streamflow discharge (Fig. 5a) is underestimated in five years; the largest difference in 2008 was likely due to an underestimation of precipitation in late 2007 (Fig. 3) reducing snow accumulation and subsequent snowmelt runoff. Simulated mean monthly streamflow (Fig. 5b) represents some of the observed seasonality, in particular the rise in May streamflow, but underestimated the rise in the end of the summer streamflow. Flow duration curves (Fig. 5c) show a good overall agreement with some errors in representing high flows. Cumulative streamflow shows that simulations rise slightly earlier than observations (Fig. 5d), and correctly simulate the rate of increase in rising discharge, but underestimate spring discharge volumes. Overall, the CRHM-AHM showed consistent streamflow response with those presented by Krogh et al. (2017); however, a larger bias was found, likely due to the imperfect interannual variability of corrected WRF precipitation. Differences between CRHM-AHM simulations using the few calibrated parameters from Krogh et al. (2017) or the ones recalibrated for this study showed only marginal differences; therefore, the parameter values presented by Krogh et al. (2017) are used as they produce a slight smaller mean bias and are presumed to be more realistic, as they were derived using observed weather.
d. Changes in snow accumulation and cover
Figure 6 presents historical and future mean and standard deviation of daily snow water equivalent (SWE) for the HPC basin and four representative land cover classes. Basin-scale peak SWE is projected to increase in 80 mm and occur at the same day of the year (19 April). Snow-cover duration is projected to shorten by 26 days, as there is a 15-day delay in the initiation of snow accumulation and 11 days earlier snow-cover depletion (Table 4). Because of the larger peak SWE and the earlier snow-cover depletion date, the snow ablation rate increased significantly from 1.8 to 3.5 mm day−1. An increase in peak SWE for most land covers of 67–83 mm was found (Fig. 6), mostly due to the increasing snowfall and warmer conditions that dampened redistribution or interception and hence sublimation. Peak SWE in the gully/drift decreased by 164 mm due to suppression of blowing snow redistribution due to warmer temperatures and increasing shrub density. The smaller increase in peak SWE in the dense forest compared to the taiga forest is due to the dense forest’s higher interception capacity, which allows for higher sublimation losses. SWE interannual variability, represented by the standard deviation in Fig. 6, does not change substantially in the future for most land covers; however, in the gully/drift, future variability is reduced to a more temporally consistent and lower snow accumulation due to increased shrubs and warmer temperatures diminishing blowing snow redistribution to the drift.
e. Changes in active layer thickness
Figure 7 shows changes associated with the active layer thickness (ALT) for selected land cover and the basin average. The basin average ALT is projected to increase by 0.25 m, thawing from approximately 0.96 down to 1.21 m in the future. Spatial heterogeneities among land cover types were found in the increased ALT, ranging from roughly 0.2 to 0.35 m between the historical and future scenarios. The 26-day increase in the snow-free season and increased air temperature can explain the deepening ALT depth (Table 4). Ground thaw initiation started roughly 5 days earlier in future simulations, shifting from 22 May to 17 May for the historical and future scenarios, respectively, consistently with the earlier depletion of snow cover. Ground freeze initiation is delayed by 17 days in the future, driven by the 2-week-later initiation of snow cover. Average thawing rates increased in the future from 0.78 to 0.83 cm day−1, driven by the longer thawing season and the warmer air temperatures. ALT interannual variability for individual land covers or over the basin does not change substantially in the future from current conditions.
f. Changes in the mass balance
Figure 8 shows mean annual cumulative fluxes for the historical and future scenarios at the basin scale. Historically, snowfall was the largest precipitation component with 172 mm (57% of mean annual precipitation), whereas rainfall contributed 133 mm (43%). Future partitioning between rainfall and snowfall shows the largest increase, 63 mm, going toward rainfall, whereas snowfall increased by only 54 mm. Despite this, snowfall remains the largest precipitation component in the future with 53% of the mean annual precipitation. Figure 8c shows that ET increased by 27 mm in the future, which, along with the 117-mm increase in precipitation, results in a reduction of the evaporation ratio (ET/precipitation) from 0.49 to 0.42. Increased ET is the result of the warmer and wetter conditions, and the larger all-wave irradiance (Table 3). Slightly higher mean sublimation rates from intercepted snowfall were found; however, cumulative intercepted sublimation dropped by 1 mm (Fig. 8d), which is not a substantial change, but when compared with the 54 mm increase in snowfall it results in a decrease from 11% to 8% of the mean annual snowfall. Decreasing total sublimation from intercepted snowfall is explained by the shortening of the snowfall season permitting snow interception on the canopy and warming air temperatures that induce more rapid and earlier unloading of canopy snow. Total blowing snow sublimation (Fig. 8e) increased by 2 mm due to the increased snowfall, which to some degree overcame the impact shrub expansion and densification had in restricting blowing snow redistribution. Daily rates of sublimation at the snow surface remain virtually the same (Fig. 8f); however, cumulative surface sublimation decreased by 6 mm due to the shortening of the snow-cover season. The cumulative sublimation decreased by 5 mm, with a substantial drop in the sublimation ratio to snowfall from 34% to 23% (Fig. 8g). Mean annual soil moisture increased by roughly 7 mm (Fig. 8h), which is attributed to the increased soil storage capacity associated with deeper ALT (Fig. 7) and the increase in precipitation that was not matched by a proportionate increase in ET or sublimation. Larger rates of soil recharge were found at the beginning of the summer; however, these were somewhat compensated by the faster soil moisture depletion projected by midsummer, resulting in virtually the same minimum soil moisture. The same minimum soil moisture in historical and future simulations is likely due to the storage in deeper layers of the soil with low permeability (mineral soil) that can hold moisture for longer periods of time. A significant increase in streamflow discharge volume, by roughly 100 mm, was projected to be driven by increasing surface and subsurface runoff (Fig. 8i), suggesting that most of the increased precipitation (117 mm) is translated into streamflow. The interannual variability of the main hydrological fluxes, represented by the standard deviation in Fig. 8, does not change except for an increase in the interannual variability in blowing snow sublimation. Years with the smaller blowing snow sublimation are less frequent in the future; however, years with higher snow accumulation have substantially larger sublimation losses as more snow is redistributed by vegetation and wind with increasing snowfall.
g. Changes in the hydrological regime
Mean monthly streamflow (Fig. 9a) shows a significant peak in spring runoff, which increased by roughly 130%, from 0.28 to 0.65 m3 s−1, and a much smaller increase in fall flows. The flow duration curve shows (Fig. 9b) an increase in daily streamflow discharge for most of the exceedance probabilities, particularly for low exceedance probabilities. For example, the 1%, 5%, and 10% exceedance probability, associated with return periods of 100, 20, and 10 years, increased by 0.6, 0.3, and 0.25 m3 s−1, respectively. Mean streamflow discharge (Fig. 9c) initiates about a week earlier in the future, consistently with earlier snow depletion, whereas the end of the streamflow discharge is delayed by 6 days. The mean annual peak flow (Fig. 9d) changed from 0.9 to 1.6 m3 s−1, whereas the date at which it occurs advanced by a week from 22 May to 15 May. Runoff ratio increased by 45% from 0.33 to 0.48.
h. Mass balance sensitivity to projected vegetation changes
To explore uncertainty and assess the mass balance sensitivity to the projected changes in vegetation characteristics, 1296 vegetation scenarios were created (Table 5). These scenarios assumed that vegetation height in those HRUs with previously existing shrubs remain the same, as there is no evidence of increasing shrub height in the region. Scenarios of shrub expansion range from a slight increase in shrub cover to a complete shrubification of the tundra. Changes in forest LAI are included to assess the potential impact of increasing forest density on the mass balance.
The left panel in Fig. 10 shows that the overall sensitivity of the water balance to vegetation change is relatively small for most mass fluxes. The hydrological flux with the largest absolute sensitivity is streamflow (corresponding to the letter i on the x axis), ranging from 186 to 206 mm yr−1 (runoff ratio between 0.47 and 0.52), followed by total sublimation (letter f) and sublimation of intercepted snowfall (letter d) ranging from 49 to 64 mm yr−1 and from 16 to 26 mm yr−1, respectively. The right panel in Fig. 10 shows the relative sensitivity of vegetation projections with respect to each mass flux as calculated with the vegetation projection presented in section 4b(2). The largest relative sensitivity is associated with evaporation from intercepted rainfall (a; 66%–135%), as the forest covers the majority of the basin. The second one is associated with blowing snow sublimation (c; 73%–124%) and is due to the effect of shrub expansion in blowing snow. The third one is from sublimation of intercepted snowfall (d; 77%–121%) due to the direct relationship between changes in forest LAI and canopy interception capacity. The sensitivity to these two sublimation terms (letters c and d) drives the sensitivity to total sublimation (f; 87% and 112%). The other mass fluxes show a lower sensitivity (<±10%).
Historical climate simulations using WRF at 4 km generally represented precipitation well at Inuvik with a small mean bias of 8 mm yr−1 or 3%, and a good agreement in the number and mean length of dry and wet spells (Table 2); however, the largest summer rainfall events were somewhat underestimated. These errors were greatly reduced with the quantile mapping bias correction (Fig. 3). A cold bias in simulated temperature of 2.7°C was found and corrected, which has also been found in other Arctic studies using high-resolution atmospheric models. Cai et al. (2018) dynamically downscaled ERA-Interim using polar WRF (Hines et al. 2011) at a 10-km spatial resolution over Alaska and found a mean cold bias of 1.4°C.
Climate projections from WRF under the PGW configuration show significantly warmer (6.1°C) and wetter (117 mm yr−1 or 39%) conditions, small changes in wind speed and relative humidity, and an increase in all-wave irradiance (19 W m−2). These results are consistent with, but slightly larger than estimations from the polar WRF forced by the CESM1.0 model under the RCP8.5 scenario in northern Alaska from Cai et al. (2018) (5°C and 25% in mean annual air temperature and precipitation, respectively). Overland et al. (2014) investigated seasonal average changes in Artic surface air temperature projections using CMIP5 ensembles mean under the RCP8.5 and RCP4.5 scenarios and found an increase of roughly 5° and 3°C, respectively. There are inherent uncertainties in climate projection due to CO2 emission scenarios, model structure and parameters, and intrinsic internal variability (Hodson et al. 2013); however, to include all the sources of uncertainties is computationally very expensive (Liu et al. 2017), particularly over large domains.
Large hydrological changes under changing climate and vegetation were projected to the end of the century. An increase of 80 mm or 70% in peak SWE, a shortening in 26 days or 9% of the snow-covered period, and a doubling in the snow ablation rate from 1.8 to 3.5 mm day−1 are projected. These changes are driven by the increase in snowfall by 54 mm yr−1 or 31% and the warmer temperatures. Callaghan et al. (2011) analyzed changes in Arctic SWE and snow-cover duration using GCMs from CMIP3 for the 2049–60 period and found a smaller increase in peak SWE of up to 15% for most of the Arctic and a larger decrease in snow-cover duration of 10%–20%. However, large-scale GCMs cannot properly resolve snow accumulation processes in environments with complex topography or where snow redistribution is important. López-Moreno et al. (2013), Musselman et al. (2017), Pomeroy et al. (2015b), and Rasouli et al. (2014) found decreasing snowmelt rates with climate warming in the western United States, Spain, and southwestern and subarctic Canada, in contrast to the accelerating melt rates found farther north in this Arctic study location. The results suggest that the impact of warming on melt rates cannot be generalized and a detailed analysis that includes the snow processes driving snow accumulation and melt needs to be considered.
A shorter snow-cover season produces earlier exposure of bare ground and later ground freeze, resulting in a longer ground thaw season and a 0.25-m-thicker active layer. Woo et al. (2007) projected a similar 0.3-m increase in the active layer thickness under the A2 scenario by 2100 for two Arctic sites with a 0.2-m peat cover. Total sublimation was projected to slightly decrease by 5 mm or 9%, mostly driven by the shortening of the snow-covered period and warmer temperatures dampening the sublimation from canopy-intercepted snowfall. The combined effect of these projected changes resulted in a significant increase of 100 mm or 100% in streamflow volume, mostly due to the doubling of spring runoff (Fig. 9). Note that the 117-mm increase in annual precipitation and the 6°C warming did not produce a particularly large increase in ET (26 mm or 18%); in fact, it decreased the evaporative index from 0.49 to 0.42. Similarly, the sublimation ratio dropped from 34% to 23%. Capturing this particular hydrological behavior is critical in cold regions, where the hydrological regime is dominated by snowmelt, as misrepresenting the timing and magnitude of precipitation can produce a completely different hydrological response (e.g., larger ET and less streamflow if precipitation increase shifts to the summer); this encourages the use of dynamical downscaling approaches. Having the appropriate hydrological processes represented realistically in the model is also important—the drop-in ratios of sublimation and ET were critical to the increase in runoff ratio and required including canopy resistance, interception of rainfall and snowfall, and blowing snow transport processes in the model.
Krogh and Pomeroy (2018) analyzed historical change at Havikpak Creek using the CRHM-AHM over 1960–2016, showing some discrepancies between the historical trends and the future projections presented in this study. Decreasing historical annual trends in ET, soil moisture and blowing snow sublimation, and a negative changepoint in streamflow volume oppose the projections presented in this study. These discrepancies are due to the large increase in precipitation projected for the future (39%) as opposed to the declining historical precipitation found by Krogh and Pomeroy (2018). Nevertheless, similarities also exist. Earlier snow-cover depletion and peak streamflow, delayed ground freeze and snow-cover initiation, and thicker active layer were found in both the historical change analysis and this study. These processes are mostly controlled by changes in air temperature, which explains consistencies between historical trends and future projections.
To the best of the authors’ knowledge, there are only two studies that have investigated the impact of climate change in Arctic headwater basins. Hinzman and Kane (1992) in Imnavait Creek, Alaska, used the HBV model and defined three scenarios of climate change by increasing observed temperature in 4°C and precipitation by 0% and ±15%; however, these simulations were performed over a 1-yr period, limiting the generalization of the results. Pohl et al. (2007) in Trail Valley Creek, Northwest Territories, Canada, used WATFLOOD and a similar approach to generate climate change scenarios for 2050 and 2080; however, mean weather changes were informed by two GCMs and two emission scenarios (A2 and B2). Pohl et al. (2007) showed results that are in the same direction as those presented in this study, such as increasing runoff volume and evaporation, and earlier spring runoff and peak streamflow. Importantly, snow redistribution and sublimation processes were not included despite their demonstrated importance in this environment (Liston et al. 2002; Pomeroy and Li 2000; Pomeroy et al. 1997), and no permafrost calculations or changes in vegetation were included. Those shortcomings have been addressed in this study, producing simulations that are expected to be more robust under future conditions. Limitations in this study arise from the use of a single future climate, due to computational costs of the PGW approach. Future studies should aim to include a larger number of future climate projections to incorporate some of the uncertainty produced by different emission scenarios and model structure. The physical consistency between climate variables provided by dynamically downscaled, high-resolution atmospheric models provided an advantage in driving physically based hydrological models such as the one used in this study; this advantage should not be abandoned in pursuit of representing uncertainty from climate models by adopting statistical downscaling.
This study presented the implementation of a spatially distributed and physically based Arctic hydrological model (CRHM-AHM) forced with historical and future dynamically downscaled weather from a high-resolution (4 km) atmospheric model (WRF) under a pseudo-global-warming configuration, as well as projections of vegetation changes. Future climate at HPC showed much warmer (6.1°C) and wetter (39%) conditions, which produced several hydrological changes, including an intensification of the hydrological regime by increasing spring runoff by 130%, producing earlier (1 week) and larger (77%) peak streamflow, increasing peak snow accumulation (70%), shortening of the snow-cover duration (26 days), increasing the melt rate (94%), thickening of the active layer thickness (0.25 m), increasing ET (18%), and decreasing sublimation (9%). These projected changes are strongly conditioned to projected changes in climate; however, projected increase in shrub cover and density also play an important role in annual ET from intercepted rainfall and blowing snow redistribution and sublimation, as revealed by the sensitivity analysis. Overall, the hydrological processes shifted from controls exerted by cold regions processes toward summer processes, though the future basin remains a snow- and permafrost-dominated cold-region basin despite the impacts of climate change.
The high-resolution WRF run provided physically consistent weather time series that did not require further downscaling and drove realistic hydrological model responses that did not require model calibration when compared with observed streamflow discharge at Havikpak Creek. This suggests that atmospheric modeling at 4-km resolution provides suitable driving weather for small-scale hydrological modeling in environments with relatively low topographic gradients, such as the one found in the delta of the Mackenzie River. Therefore, future climate from WRF under the PGW configuration is expected to be robust and suitable for hydrological applications; however, the main limitation of this approach remains the large computational and human resources required, currently restricting the number of future projections to one simulation. This study provides the first detailed investigation of projected changes in climate and vegetation on the hydrology and hydrological processes of an Arctic headwater basin, which is expected to help inform other Arctic climate impact studies.
Financial support for this study was provided by CONYCIT under the BECAS CHILE scholarship program, Global Water Futures, Changing Cold Regions Network, NSERC Discovery Grants, Canada Research Chairs, and Yukon Environment. The authors thank Tom Brown of the Centre for Hydrology for providing technical support with the CRHM platform. Lucia Scaff and Zhenhua Li from the University of Saskatchewan provided valuable support to process WRF outputs.