Main

Recent global-scale estimates using simple models (for example, refs 3, 4, 5, 6) indicate that mountain glaciers could raise sea level by 0.39 m by 2100 (ref. 7). At regional-to-local scales efforts to project glacier mass changes have varied from models that apply glacier dynamics configured for single ice masses (for example, refs 8, 9) to those with greater geographical extent that rely on empirical scaling10,11, scaling in combination with a low-order treatment of ice dynamics12,13,14 or sub-grid parameterizations15. At these spatial scales the main effects of deglaciation are associated with changes in the hydrologic cycle16,17 and consequent impacts on water availability, aquatic habitat, hydroelectric power generation, recreation and tourism.

Projections of glacier surface mass balance (accumulation and ablation) can reveal the ultimate fate of glaciers, but they lack information on rates of change of thickness and extent. Glaciers individually respond to changes in the surface mass balance field and may survive an adverse climate by stabilizing at a higher elevation. This stabilization due to changes in glacier hypsometry (ice area altitude distribution) has been represented through scaling empiricisms in all current models of glacier evolution on regional and global scales5,6,8,15. A common feature of these models is that they lack a physics-based treatment of glacier dynamics. The central contribution of our study is thus to simulate the changes in ice thickness and extent over a large region using a high-resolution model of glacier dynamics, which yields year-to-year changes in ice area and volume for the entire study region.

Our study area is Alberta and British Columbia (BC) in western North America (Fig. 1), where glaciers account for an estimated area of 26,700 km2 (ref. 18) and volume of 2,980 km3 (ref. 19). The geographical scale is comparable to that of other glacierized mountain regions, such as South America (∼31,900 km2; ref. 7), the Himalaya and Karakoram (∼22,800 km2 and ∼18,000 km2; refs 16, 20) and Tien Shan (∼16,400 km2; ref. 21), where declining glacier melt will impact populations. In this study we project the evolution of regionally resolved glaciers from the present to 2100, using a regional glaciation model (RGM; Methods and Supplementary Sections 1 and 2): a high-resolution (200 m) glacier surface mass balance model coupled to a physics-based model of glacier dynamics. Projections are cast in the modelling framework of the Intergovernmental Panel on Climate Change’s Fifth Assessment Report (IPCC AR5; ref. 7) and use climate projections from six well-performing General Circulation Models (GCMs; refs 22, 23) forced by the four AR5 emissions scenarios. These scenarios, referred to as representative concentration pathways (RCPs), are labelled RCP2.6, RCP4.5, RCP6.0 and RCP8.5, where the numbers indicate the increase in radiative forcing (W m−2) by 2100 relative to pre-industrial values.

Figure 1: Study region and subregions in the Canadian Cordillera of western Canada.
figure 1

Present-day (2005) glacier extent is indicated in white. The yellow rectangle indicates the location of Columbia Reach drainage basin. Inset: Study region (black rectangle) within northwestern North America.

Full size image

Model skill can be assessed by comparing ice hypsometry from a 2005 inventory with modelled results from the same year (Fig. 2), as well as by comparing the number and area of observed and modelled ice masses (ref. 18 and Supplementary Table 3). The model is spun up from an ice-free state at year 0 to reach a quasi-steady state at 1901; subsequently, historical (1902–1979) and reanalysis (1980–2008) climate data are used (Methods and Supplementary Section 2.1). A comparison between spin-ups using steady and stochastic forcing confirms that by 1980 both procedures lead to the same result. Applying the observations of ice extent to a digital elevation model (DEM; Supplementary Section 1) together with an estimate of subglacial topography19 allows modelled and estimated volumes to be compared. The area comparisons show fractional errors of +17.8% for the Coast, −3.6% for the Interior, −2.9% for the Rockies and +14.1% for All. The fractional discrepancy between modelled and estimated ice volume (Supplementary Table 3) can be large (+60.7% for All), but the influence of these errors on normalized projections of area and volume (Fig. 3) is surprisingly small (Supplementary Section 6 and Figs 41 and 44).

Figure 2: Comparisons of observed and modelled ice hypsometry for reference year 2005.
figure 2

For each of the ten subregions, pairs of bell-shaped curves show the distribution of ice area (normalized to the largest observed area) with elevation for the observed (green) and modelled (blue) ice extent.

Full size image
Figure 3: Projected changes for glaciers in the western Canadian study region.
figure 3

a–d, Area (left) and volume (right) are normalized to modelled values for the reference year 2005: projections for Coast glaciers (a), Interior glaciers (b), Rockies glaciers (c) and all glaciers (d). The mean (solid curves) and ±1σ limits (vertical hatching) for the multi-model GCM ensemble are plotted for four different emissions scenarios (RCP2.6, RCP4.5, RCP6.0, RCP8.5). NARR denotes the North American Regional Reanalysis.

Full size image

Large-scale projections of the model (Fig. 3) are summarized for the Coast, Interior, Rockies and All regions. Glaciers of the Coast region are most resistant to climate change. For these glaciers, depending on the scenario, the ensemble means indicate that, relative to 2005, 75 ± 10% of the 2005 ice area and 70 ± 10% of the volume will be lost by 2100 (here and throughout, error ranges are for model mean ±1σ). For the Interior and Rockies regions ice area and volume losses will exceed 90% of the 2005 amounts for all scenarios except RCP2.6. The resistance of Coast glaciers to climate-forced changes in area and volume is associated with subregions 1 (St Elias), 2 (Northern Coast) and 4 (Southern Coast), which have the highest present-day ice content (Supplementary Table 3). The remaining Coast subregions will experience total or near-total losses of ice area and volume. Comparison of area and volume projections for all GCMs and scenarios revealed that the MIROC-ESM GCM (Supplementary Table 1) most frequently represented the median member of the GCM ensemble. To examine the detailed spatio-temporal character of the model projection, we therefore identified six focus sites and extracted projection modelling time snapshots for the MIROC-ESM projections. The message that the magnitude of deglaciation will be significant does not differ from that of Fig. 3, but the visual consequences of ice loss are emphasized (Supplementary Figs 27–32).

The potential sea-level rise from glacier loss in the study region (6.3 ± 0.6 mm; ref. 19) is modest, but the hydrologic implications of the projected loss are substantial. As an example, the Columbia River, which flows from its headwaters in Interior BC to the Pacific Coast of Washington and Oregon, yields the largest hydroelectric production of any river in North America. The hydroelectric generating capacity of the Canadian headwaters exceeds 5 GW (ref. 24). In BC, the Columbia Reach drainage basin (Fig. 4a) is the most intensely glacierized part of the Columbia River Basin and thus the most susceptible to changes in ice cover. The main influence of glacier runoff in the Columbia River Basin is to maintain stream flow, which contributes water to the Mica Reservoir for hydroelectric power generation (Fig. 4a), and to regulate water temperature through summer months25,26. The Columbia River Treaty between Canada and the USA, signed in 1964 and renegotiable from 2024 onwards, provides a detailed framework for cooperation on hydroelectric power generation and flood control. Climate-forced changes in water availability will redistribute costs and benefits between the treaty partners.

Figure 4: Projected changes for glaciers in the Columbia Reach drainage basin within the Columbia River Basin of British Columbia.
figure 4

a,b, Observed (a) and modelled (b) glacier extent in the Columbia Reach drainage basin (yellow outline) for the reference year 2005. The Mica hydroelectric dam in the northwest quadrant is indicated by a light blue bar, with an arrow to indicate the water flow direction. c, Model projection for 2050 using MIROC-ESM with RCP2.6 forcing. The ice extents for RCP4.5, RCP6.0 and RCP8.5 (not shown) are virtually the same as for RCP2.6. d, Model projections for 2100 using MIROC-ESM with RCP2.6, RCP4.5, RCP6.0 and RCP8.5 forcings. e, Projected changes in ice area for the drainage basin. f, Projected changes in ice volume for the drainage basin. The graphs shown in e and f are for an ensemble of GCMs and four emissions scenarios; ensemble means are indicated by solid lines and the ±1σ limits by vertical hatching. The relative changes (right ordinate) are normalized to the modelled values for the 2005 reference year. g, Projected changes in mean annual meltwater input from glacier ice loss for the RCP2.6 scenario and an ensemble of GCMs. h, Projected changes in mean annual meltwater input from glacier ice loss for the RCP8.5 scenario and an ensemble of GCMs. The maps show the northern half of the drainage basin, which contains ∼85% of the present-day ice cover. The time series are the summed contributions for the entire basin.

Full size image

The observed (Fig. 4a) and the modelled (Fig. 4b) Columbia Reach ice cover for 2005 broadly agree with the observed and modelled area (733 km2 and 666 km2) and volume (44 km3 and 36 km3). Widespread glacier loss occurs by 2050 (Fig. 4c), with near-total ice disappearance by 2100 (Fig. 4d) for the MIROC-ESM GCM and the RCP8.5 scenario. Using ice area and volume projection results from the six GCMs (Supplementary Table 1) and four AR5 scenarios, we present time series for the averages and ±1σ ranges for the multi-model ensemble (Fig. 4e, f). Until mid-century (∼2050), the fate of all glaciers in this area is virtually independent of the emission scenario and climate model used for the projections. By the end of the century, however, the ensemble averages range from ∼70% (RCP2.6) to ∼95% (RCP8.5) reduction of both area and volume relative to 2005 values. The rate of change in ice volume (Supplementary Section 4.5) for each GCM and the RCP2.6 (Fig. 4g) and RCP8.5 (Fig. 4h) scenarios yields the projected changes in meltwater input from glacier ice loss. The graphs clearly show the effect of an unsustainable ‘deglaciation discharge dividend’ (ref. 27). For the majority of GCMs and scenarios, runoff from glacier wastage is characterized as a well-defined peak in meltwater discharge, having a typical amplitude of ∼15 m3 s−1, roughly 3% of the ∼500 m3 s−1 annual average discharge of Columbia River at the Mica Dam, followed by decades of declining flow. Our simulated deglaciation discharge (Fig. 4g, h) supplements the annual cycle of glacier storage and melt, and corresponds to annual average rates of mass loss that would mainly occur in summer and early fall. Meltwater discharge projections for all GCMs and scenarios and for all regions and subregions, as well as for the Canadian Columbia River Basin, are also calculated (Supplementary Sections 4.5 and 4.7). Most of these runs indicate a clear peak in discharge, in marked contrast to recent runoff projections28 obtained using a model without ice dynamics that show no peak for western Canada.

Uncertainty in the RGM projections results from the uncertainty concerning which emissions pathway will be followed, from the range of GCM projections of future climate and from shortcomings of the surface mass balance model, the estimated subglacial topography and the ice dynamics model. Uncertainties associated with GCM projections are examined in the AR5 (ref. 7), but for mountainous regions additional uncertainty is contributed by orographic effects on precipitation and temperature29 (Supplementary Section 2.2). Uncertainties associated with the mass balance model are dominated by the ensemble variability within each scenario rather than by limitations of the model. The error contribution of the ice dynamics model is small relative to that from the surface mass balance treatment (Supplementary Section 9). Our sensitivity analysis (Supplementary Tables 5–10) indicates that projections of ice area and volume have low sensitivity to reasonable parameter assignments for the ice dynamics model but high sensitivity to parameters of the surface mass balance model. Ice dynamics calculations increase the computational demands but do not greatly complicate the projection methodology. Projections performed with and without ice dynamics have been compared (Supplementary Section 7); when ice dynamics are neglected the RGM systematically underestimates ice volume loss at 2100 (47 ± 12% loss without dynamics and 60 ± 10% with dynamics for ‘All’ and RCP2.6; Supplementary Fig. 48). Comparisons of observed and modelled ice extents offer a powerful approach to refining surface mass balance fields, but one that is possible only using models that include dynamics. The effects of introducing a mass balance bias correction were carefully examined (Supplementary Section 6). The main challenge remains that of improving the surface mass balance treatment.

In addition to the hydrologic implications of reduced late-summer surface flows in the Columbia Basin, our projected changes of ice cover in western Canada have broader ramifications for aquatic ecosystems, agriculture, forestry, alpine tourism, water quality and resource development. Our study departs from previous large-scale projections of glacier response to climate change by including the contribution of ice dynamics and by operating at high spatial resolution despite the large size of the study region. By archiving the deglaciation projections for all the GCMs and scenarios for the entire study region (http://www.unbc.ca/research/supplementary-data-unbc-publications) we open the possibility for a wide range of local- and regional-scale impact studies. With appropriate modifications of the climate projections and surface mass balance models, the present work provides a template for an assessment of future mass change in Earth’s other glacierized mountain regions.

Methods

The regional glaciation model (RGM) combines a surface mass balance model that quantifies mass fluxes (accumulation and ablation) at the glacier and land surfaces with an ice dynamics model to simulate ice flow (Supplementary Figs 2 and 3). Both models rely on gridded representations of the present-day surface elevation and ice extent. Data from the 2000 Shuttle Radar Topography Mission (available at http://srtm.csi.cgiar.org) were reprojected and resampled at a resolution of 200 m to obtain a DEM of surface topography; glacier outlines from around 2005 (available at http://www.glims.org/RGI) were used to derive a rasterized and co-registered ice mask. From these and other data, a DEM of the hidden subglacial topography was estimated using an optimization method (Supplementary Section 1.1).

The surface mass balance model calculates ablation by a distributed temperature-index model, accumulation using a temperature threshold to differentiate snow from rain precipitation and refreezing following a parameterized thermodynamic approach (Supplementary Section 2.3). To calibrate the model we force it with downscaled monthly climate fields derived from the North American Regional Reanalysis for 1980–2008 (NARR; ref. 30). NARR temperature and precipitation fields are downscaled from a horizontal resolution of 32 km to 200 m using methods described elsewhere29. Model calibration consists of tuning the model parameters to minimize a misfit between modelled and observed glacier mass balances derived from all available in-situ and geodetic measurements in the region within the period 1980–2008 (Supplementary Section 2.3). To obtain past and future forcings we take the downscaled NARR monthly averaged temperature and precipitation fields for 1980–2008 to represent the baseline of these fields and superimpose the anomalies from the CRU data set for 1902–1979, and from an ensemble of GCM outputs for 2009–2100 (Supplementary Section 2.5). Forced by these climate fields, the mass balance model yields annual mass balance fields at a resolution of 200 m for 1902–2100 which are then coupled with the glacier dynamics model. The coupling takes into account the feedback mechanisms between glacier mass balance and glacier geometry changes (for example, positive feedback between glacier thinning and mass balance; negative feedback between glacier shrinking and mass balance), while the dynamics model redistributes the mass through ice flow.

The glacier dynamics component of the RGM assumes the shallow-ice approximation and isothermal ice. The model is 2.5D (two-dimensional vertically integrated) with a grid spacing of 200 m. Evolution equations for surface elevation are approximated as finite-difference expressions and solved using a super-implicit numerical scheme and flux-limiters (Supplementary Section 1.2). The ice dynamics model is spun up from an ice-free state at year 0 by randomly selecting the modelled annual mass balance for each year in the range 1902–1931 and applying this forcing from year 0 to 1901, by which time a quasi-steady state (‘stochastic equilibrium’) is attained. Our simulations start from this modelled 1901 ice configuration, progress through the twentieth century to 2009, and are then projected from 2010 to 2100. Model skill is assessed by comparing the observed ice hypsometry from a 2005 inventory with modelled results for the same year. To improve the model skill, we apply a bias to the surface mass balance field, following similar bias-correction strategies as in other projection studies (for example, refs 3, 6). The bias is additive, spatially varying and constant in time (Supplementary Section 2.7). RGM simulations forced with biased and unbiased mass balance fields exhibit differences in the total area and volume of modelled ice, but similar population statistics for glaciers modelled with and without the bias (Supplementary Section 6). Projection results for 2009–2100 for all GCMs and scenarios of our study are archived at http://www.unbc.ca/research/supplementary-data-unbc-publications.

Code availability.

The code is not available.