Using GRACE and climate model simulations to predict mass loss of Alaskan glaciers through 2100

Abstract

Glaciers in Alaska are currently losing mass at a rate of ~−50 Gt a⁻¹, one of the largest ice loss rates of any regional collection of mountain glaciers on Earth. Existing projections of Alaska's future sea-level contributions tend to be divergent and are not tied directly to regional observations. Here we develop a simple, regional observation-based projection of Alaska's future sea-level contribution. We compute a time series of recent Alaska glacier mass variability using monthly GRACE gravity fields from August 2002 through December 2014. We also construct a three-parameter model of Alaska glacier mass variability based on monthly ERA-Interim snowfall and temperature fields. When these three model parameters are fitted to the GRACE time series, the model explains 94% of the variance of the GRACE data. Using these parameter values, we then apply the model to simulated fields of monthly temperature and snowfall from the Community Earth System Model, to obtain predictions of mass variations through 2100. We conclude that mass loss rates may increase between −80 and −110 Gt a⁻¹ by 2100, with a total sea-level rise contribution of 19 ± 4 mm during the 21st century.

1. INTRODUCTION

Global sea-level rise is likely to be one of the most socially disruptive consequences of future atmospheric warming. At present, about 1/3 to 1/2 of the observed rate of sea-level rise is caused by thermal expansion of the oceans, with the rest caused by mass loss from the world's ice-covered regions: Greenland, Antarctica, small ice caps and mountain glacier systems (Church and others, 2013). In the distant future, West Antarctica and Greenland will likely contribute the majority of the new sea-level rise (Velicogna, 2009; Joughin and others, 2014). For example, the Amundsen Sea portion of West Antarctica alone contains enough ice to raise sea level ~1.2 m (Rignot and others, 2014), whereas the entire volume of Alaska's glaciers, if melted, would only raise sea level ~0.04 m (Grinsted, 2013). Still, the relatively rapid response of mountain glacier systems to atmospheric warming implies that through 2100, contributions from glaciers are likely to be relatively important. Over the first decade of the 21st century, Alaskan glaciers contributed between 50 and 75 Gt a⁻¹ to global sea-level rise (e.g. Arendt and others, 2002; Jacob and others, 2012; Gardner and others, 2013; Luthcke and others, 2013; Sasgen and others, 2013; Larsen and others, 2015), equivalent to ~60–90% of the contribution from Antarctica over that same time period (e.g. Shepherd and others, 2012).

Thus, it is of interest to examine the evolution of Alaskan glaciers when trying to project sea-level rise over the next century. Several previous studies have addressed this problem. Radić and Hock (2011), Slangen and others (2012), Giesen and Oerlemans (2013), and Radic and others (2013) modeled all the world's glacier systems, not just Alaska, and projected the mass loss from those systems through the 21st century. Each of these studies constructed a parameterized mass-balance model for individual glaciers that allowed an estimate of how a glacier's volume will evolve with time from knowledge of atmospheric variables over that glacier. Parameters in these models were generally determined through calibration against a small sample of glacier mass-balance observations, world-wide. These models were then applied to all glaciers in every glaciated region outside Greenland and Antarctica, and were forced with global climate model output through (usually) to 2100.

Calibrating models with in situ observations, as the above studies have done, has some distinct advantages. However, this type of calibration is also heavily dependent on a very small sample size with known biases (Gardner and others, 2013). In Alaska, direct observations are made on only a handful of glaciers out of >27 000 glaciers statewide. Moreover, the glacier-to-glacier SD of mass balance in Alaska is approximately half the average regional mass balance (Larsen and others, 2015). Thus, calibrating models to sparse in situ observations brings with it the possibility of large biases.

A solution to this problem is to calibrate regional glacier models to regional observations rather than in situ data. Previously, Arendt and others (2009) took such an approach by calibrating two simplified mass-balance models to GRACE (Gravity Recovery and Climate Experiment) gravity observations from glaciers in Alaska's St. Elias region. They found that due to the coarseness of GRACE results, models incorporating elevation complexities could not be well constrained. They found that simpler models that do not explicitly incorporate elevation performed better, accounting for up to 60% of the variance in the GRACE time series. Their model, however, was forced by meteorological station data that are subject to local spatial variability in precipitation and temperature. Temperature data from climate re-analyses, however, has been found to correlate well the GRACE time series (Arendt and others, 2002).

Our approach aims to see if a simple mass-balance model, forced with regional re-analysis data, can more closely approximate Alaska-wide GRACE observations than a model forced with local meteorological station data. By considering GRACE results averaged over a larger region, we reduce possible problems caused by contamination from signals in glaciers outside the region of interest. Also, as in Arendt and others (2009), we use output from a hydrological model to subtract the effects of soil moisture and of seasonal snow mass. However, we do not remove model contributions of seasonal snow mass from within glaciated regions, under the assumption that the atmospheric snowfall fields we calibrate against contain contributions from snowfall both over glaciers and over unglaciated land within the glaciated regions.

We derive a parameterized model of glacier mass variability by constructing a simple accumulation/melt model that is independent of elevation and forced by monthly ERA-Interim (ERAI) (Dee and others, 2011) snowfall and temperature fields. We calibrate our model against observations that cover the entire glaciated region, using regional estimates of mass loss determined using monthly GRACE gravity fields (Jacob and others, 2012) from August 2002 through December 2014. The model involves only three adjustable parameters, which we tune so that the model results best match the GRACE estimates.

Despite its simplicity, this model accounts for 94% of the variance in the GRACE data. The model replicates monthly variations including a secular decrease, seasonal variations and interannual fluctuations. Given the quality of this fit, we experiment with forcing the model with Community Earth System Model (CESM; Hurrell and others, 2013) monthly projections of future temperature and precipitation, to predict how Alaskan glacier volumes might change through 2100, assuming plausible future greenhouse gas emission levels.

2. THE GRACE DATA

The GRACE satellite mission, launched in the spring of 2002 by NASA and Deutsches Zentrum für Luft- und Raumfahrt (the German space agency), has been providing users with monthly, global estimates of the Earth's gravity field in the form of spherical harmonic coefficients (e.g. Tapley and others, 2004; data available at http://podaac.jpl.nasa.gov). These gravity fields have been used to study a variety of geophysical processes that involve changes in the Earth's mass distribution (e.g. Wouters and others, 2014). One such process is the change in the snow mass and ice in the world's mountain glacier systems. Here, we estimate the mass change in Alaskan glaciers, following the methods described by Jacob and others (2012) and summarized briefly here.

We use monthly gravity fields provided by the Center for Space Research (CSR) at the University of Texas, in the form of spherical harmonic coefficients. We replace the monthly GRACE values of C ₂₀ (the zonal, degree-2 spherical harmonic coefficient of the geopotential) with the more accurate values provided by satellite laser ranging (Cheng and Tapley, 2004), and we include degree-1 terms using the method described by Swenson and others (2008). We remove contributions from glacial isostatic adjustment (GIA), by subtracting the GIA model of A and others (2013).

We solve for changes in snow/ice mass over Alaska's glaciated regions, by covering the region with ‘mascons’ (small, arbitrarily defined areas of the Earth) and fitting mass values for each mascon to the GRACE gravity fields. Each mascon is defined as the union of points on a 0.5° latitude/longitude grid. Figure 1 shows the points we use to define glacier mascons in Alaska. There are seven such mascons, composed of a total of 216 0.5° grid points (covering a total area of 3.3 km × 10⁵ km) (Jacob and others, 2012). We sum the time series for all seven Alaskan mascons to obtain the total Alaska time series shown by the black curve in Figure 2. The best-fitting trend of that time series is −52 ± 4 Gt a⁻¹. This rate has overlapping error bars with Luthcke and others (2013) and Arendt and others (2002) and is somewhat lower than Larsen and others (2015). The range of results is due to a combination of differences in spatial extent, time period considered and differences in GRACE processing and terrestrial water storage models (Jacob and others, 2012; Luthcke and others, 2013).

Fig. 1. Locations of the seven Alaskan glacier mascons, shown in different colors, and the half-degree latitude/longitude points that make up those mascons.

Fig. 2. Monthly results for the total mass of the mascons, from August 2002 through December 2014. GRACE results are shown in black. Results from the best-fitting accumulation minus melting model are shown in orange. The difference, shown in blue, has been offset from 0 to make it easier to distinguish.

We also generated GRACE solutions where we separated the Figure 1 Alaskan grid points into subsets of points and solved for mass variations of each of those subsets. This allowed us to focus on smaller glaciated regions. The difficulty with this approach was partly that the solutions became less accurate, and partly that the resulting subregions were close enough to one another that GRACE could not cleanly separate one subregion from another, and so mass variations in one subregion contaminated the GRACE estimate of mass variations in other subregions. Moreover, while glacier mass balance has been found to vary substantially glacier-to-glacier, variations between regions are small (Larsen and others, 2015). For these reasons, we elected to concentrate on the total Alaskan glaciated region.

Before solving for the mass variations, we subtract (from the GRACE data) monthly water storage estimates predicted by the Community Land Model (CLM4.5) land surface model (Oleson and others, 2013), to minimize the contamination from non-glacier water storage signals. Specifically, we subtract CLM4.5 estimates of soil moisture, vegetation and river water storage at every grid point. Over every CLM4.5 grid point that does not correspond to a glacier mascon point, we further remove CLM4.5 estimates of snow water storage. Thus, the GRACE mascon time series (Fig. 2) includes contributions from changes in snow and ice on glaciers within those mascons, as well as contributions from snow variability over unglaciated land within the mascons.

The removal of CLM4.5 is done in the spherical harmonic domain. Namely, we expand the relevant water storage components of CLM4.5 into harmonics for each monthly CLM4.5 field, and we subtract the resulting monthly CLM4.5 coefficient of each harmonic from the monthly GRACE value of that harmonic. We then fit the monthly mascon solutions to the resulting GRACE-minus-CLM4.5 harmonics.

3. THE MASS CHANGE MODEL

The GRACE mass variability shown in Figure 2 is the result of accumulated snowfall, which adds mass, and of ablation – principally through melting – which removes mass. Tidewater calving in Alaska, we consider small enough to exclude as they contribute to only ~7% of the total mass loss from Alaska and calving is not expected to grow in the next century (Larsen and others, 2015). We also exclude a term to represent meltwater refreezing as GRACE is insensitive to internal vertical redistribution of mass. In addition, refreezing contributions are absorbed into our numerical solution for the degree-day factor (DDF), described below.

Temperature index models are commonly used to estimate melt at discrete points on a glacier or snow covered surface (Hock, 2003). Below 0°C, such models assume no melting and above 0°C melting are assumed to be linearly proportional to the number of degrees above zero where the proportionality constant is referred to as the DDF (Hock, 2003). Because GRACE does not recover mass values at individual points, but provides values averaged over an entire region using a traditional point, implementation of a temperature index model is ineffective (Arendt and others, 2009). Thus we design a simple, ad hoc model based on this proportionality of melting to the degree-day sum, but that represents averages over the entire glaciated region. We combine a temperature index model with a model for accumulation to develop a simple model for Alaska glacier mass variability. The model involves three adjustable parameters, and we find the numerical values of those parameters that give the best fit to the GRACE data. Appendix A provides a list of all variables used in our model.

For this model we use monthly, gridded, ERAI reanalysis fields of precipitation rate and temperature, available from the European Centre for Medium-Range Weather Forecasts (ECMWF). ERAI precipitation is partitioned into the sum of snowfall rate and liquid rainfall rate fields using a temperature-based algorithm from CLM4.5 (Oleson and others, 2013).

Let SF (t) be the resulting snowfall rate, as a function of time, averaged over the entire GRACE mascon region (covered by the x's in Fig. 1). Let T (t) be the ERAI temperature values, similarly averaged over the mascon region. We assume the snowfall accumulation rate averaged over the mascons is K ₀SF (t), where the factor K ₀ is introduced to account for any bias between the true snowfall rate and the ERAI/CESM rate, and will be estimated when fitting to the GRACE time series. K ₀ may differ from 1.0 because of a bias in the ERAI precipitation, or a bias in how the CESM model separates that precipitation into rain and snowfall.

We approximate the melting rate, averaged over the mascons using a temperature index method (Hock, 2003) generalized to each mascon. We assign no melt in a mascon when the ERAI model surface temperature, (T (t)) is below a reference temperature, T ₀. For T (t) >T ₀, the melting rate is linearly proportional to the degrees above T ₀, and is scaled by a DDF to have the form DDF (T(t) − T ₀).

The application of a temperature index model to an entire mascon carries key differences from the typical implementation (Hock, 2003). In particular, the DDF and T ₀ values will depend on the elevation distribution within each mascon (which influences the spatial variability of melt and precipitation within the mascon), and ERAI biases. Because the ERAI temperature represents the regional average temperature, T ₀ does not simply represent freezing; rather it may be offset from 0°C due to differences between the average temperature of the glaciated areas within the mascon and that within the ERAI pixel. For example, within each mascon the snow is concentrated at higher elevations, where temperatures are lower than the mascon mean temperature. At these elevations, melt will begin after the mascon average temperature climbs above 0°C, which will correspond to an ERAI temperature i.e. >0°C. We assume the values of the parameters K ₀, T ₀ and DDF are the same for glaciers as for snow over unglacierized land within the mascons.

Let S _gl be the total mascon area that is covered by glaciers, and let S _ugl be the total mascon area that is covered by unglacierized land. We use hypsometry (area-vs-elevation) data for every Alaskan glacier (Larsen and others, 2015), and sum them over all glaciers that lie within the mascons shown in Figure 1, to obtain an estimate of the total glaciated area lying in each 30 m elevation band from sea level up to 6165 m. This hypsometry is used to estimate the change in K ₀ and T ₀ as the regional glacier mass changes (Appendix C). Integrating over all elevations, we find that S _gl  = 7.4 km × 10⁴ km: ~22% of the total mascon area, implying that the unglaciated portions of those mascons cover ~78% of the total mascon area.

Let M _gl (t) and M _ugl (t) be the total mass over glaciated and unglaciated land relative to the mass at time 0. Let C _gl and C _ugl be the accumulation rates integrated over glaciated and unglaciated land, and let A _gl and A _ugl be the melting rates integrated over glaciated and unglaciated land. Then, we assume:

(1) $$\displaystyle{{{\rm d}{M_{gl}}} \over {{\rm d}t}}(t) = {C_{gl}}(t) - {A_{gl}}(t)$$

and

(2) $$\displaystyle{{{M_{ugl}}} \over {{\rm d}t}}(t) = {C_{ugl}}(t) - {A_{ugl}}(t)$$

where

(3) $${C_{gl}}(t) = {K_0}{\rm SF(}t{\rm )}{S_{gl}}$$

(4) $${C_{ugl}}(t) = {K_0}{\rm SF(}t{\rm )}{S_{ugl}}$$

and

(5) $${A_{gl}}(t) = \left\{ {\matrix{ 0 & {{\rm if}\quad T(t) \lt {T_0}} \cr {{\rm DDF}(T(t) - {T_0}){S_{gl}}\; \; \;} & {{\rm if}\quad T\;(t) \gt {T_0}} \cr}} \right\}$$

(6) $${A_{ugl}}(t) = \left\{ {\matrix{ 0 & {{\rm if}\quad T(t) \lt {T_0}} \cr {{\rm DDF}(T(t) - {T_0}){S_{gl}}\; \; \;} & {{\rm if}\quad T\;(t) \gt {T_0}} \cr}} \right\}.$$

Each of the functions (dM _gl /dt)(t), (dM _ugl /dt)(t), C _gl (t), C _ugl (t), A _gl (t) and A _ugl (t) represent a time rate of change of mass. But GRACE mass results, Mass (t), represent the total mass at time t. That total mass is obtained by temporally integrating the rate of change. So

(7) $$(t) = \mathop \int \limits_0^t \left( {\displaystyle{{{\rm d}{M_{gl}}} \over {{\rm d}t}}(\tau ) + \displaystyle{{{\rm d}{M_{ugl}}} \over {{\rm d}t}}(\tau )} \right){\rm d}\tau + {M^0}$$

where M ⁰ is the mass at time t = 0 (and cannot be determined using either this simple model or the GRACE data). For $M_{gl}^0 $ we use Grinsted's (2013) value:

$$M_{gl}^0 = V_{gl}^0 \times {\rho _{ice}} = 1.6 \times {10^4}\;{\rm k}{{\rm m}^3}$$

Up to this point, there is no difference between the models over glaciers and unglaciated land. But we now impose the additional constraint that melting over unglaciated land vanishes once all the accumulated snow from the previous season has melted. Thus, at times t when

$$\int_0^t {({C_{ugl}}(\tau) - {A_{ugl}}(\tau)){\rm d}\tau \lt 0}$$

(which happens every summer), we change A _ugl (t) so that

$$\int_0^t {({C_{ugl}}(\tau ) - {A_{ugl}}(\tau)){\rm d}\tau = 0}$$

We impose no such constraint on A _gl (τ). As a result, the unglaciated contribution to GRACE (t) i.e. $\int_0^t {({C_{ugl}}(\tau ) - {A_{ugl}}(\tau )){\rm d}\tau + M_{ugl}^0}$ has seasonal and other short-period variations, but exhibits no secular trend. In principle, this need not be true. It could be, for example, that if there is a winter with an unusually large snowfall, or a summer where the temperatures are unusually low, some of the snow that fell on unglaciated land might not have melted by the end of the following summer. This can lead to long-period variability, but in practice it did not occur in our solution below. Instead, we find that the unglaciated contribution to GRACE falls to zero every summer, thus exhibits no secular trend. Instead, the entire trend comes from the glacier contribution:

(8) $${M_{gl}}(t) = \int_0^t {({C_{gl}}(\tau ) - {A_{gl}}(\tau )){\rm d}\tau + M_{gl}^0} $$

The glaciers can also contribute seasonal and short-period variability.

4. SOLVING FOR THE MODEL PARAMETERS

To find the best-fitting values of DDF, T ₀ and K ₀ we perform a grid search, where we vary the values of those three parameters over a wide range, and find the set of values that minimizes the variance

(9) $${\rm variance} = \mathop \sum \limits_t {(GRACE(t) - M(t))^2},$$

where GRACE (t) is the GRACE time series, temporal averages have been removed from both GRACE (t) and Mass (t), and the sum over t in (9) represents the sum over the monthly values. We find:

(10) $${\matrix{{{\rm T}_0 = 2.2^{\circ} C.} & {1 - \sigma \;{\rm interval:}\; [2.1,\;2.3]} \hfill \cr {\quad \quad \quad \quad {\rm DDF} = 2.3\;{\rm mm}\;^{\circ} {{\rm C}^{ - 1}}\;{{\rm d}^{ - 1}}.} & {1 - \sigma \;{\rm interval:}\; [2.2,\;2.4]} \hfill \cr {{K_0} = 0.69\quad } & {1 - \sigma \;{\rm interval:}\; [0.66,\;0.74].} \hfill \cr }}$$

The 1 − σ uncertainty intervals in (Eqn (10)) are computed by adding synthetic Gaussian noise to the GRACE time series, and repeating the fit. We do this for 500 different synthetic Gaussian noise time series, one at a time, and we look to see what range of parameter values spans 68% of the 500 sets of results (68% corresponds to a 1-σ confidence interval). When constructing the synthetic noise, we have to adopt a value for the variance of that noise. We find (see the following paragraph) that by using the preferred parameter values given in (Eqn (10)), we are able to explain 94% of the variance in the GRACE data. We thus assume that the variance of the Gaussian distribution is the 6% (=100 − 94%) variance of the GRACE-minus-model residuals. By using that 6% variance, rather than the variance of the errors in the GRACE data, we are trying to account at least partially not only for the GRACE errors, but also for deficiencies in the simple model itself.

Figure 2 shows the results for M (t) (in orange) when these best-fitting values are used in Eqns (3)–(7). The agreement with the GRACE data is excellent considering that there are only three free parameters and the model parameterizes all elevation-related complexities. The model explains 94% of the variance of the GRACE results, and provides a good match to the secular trend and the seasonal terms, and captures much of the interannual variability. The most noticeable discrepancy between the two time series occurs in the first half of the year, when the model tends to overestimate peak accumulation relative to GRACE.

The difference (Fig. 2, in blue) between GRACE and the model shows that there is virtually no residual trend. Instead, the residuals consist almost entirely of seasonal variability (the model slightly overestimates the annual cycle) and short-period scatter. Since it is the glacier contribution to the model that determines the secular trend (Fig. 3), the implication is that the values of DDF, T ₀ and K ₀ obtained from the least squares fit are probably more representative of the glaciers than of snow over unglaciated land. We hypothesize that the values of those parameters for unglaciated land are close enough to those for glaciers that, by using the same parameter values for the unglaciated land, we get good agreement with GRACE at seasonal periods, in addition to obtaining a good match for the secular trend. We have tried, incidentally, including six parameters in the fit: DDF, T ₀ and K ₀ for unglaciated land, in addition to DDF, T ₀ and K ₀ for glaciers, but the correlations between these two sets of parameters are high enough that the uncertainties of the parameter solutions are much larger than those given in Eqn (10), and so we have elected to proceed with the three parameter solution.

Fig. 3. Results, from the best-fitting model, for mass variability on glaciers (orange) and over unglaciated land (blue).

5. PROJECTIONS

We apply this same model of glacier mass variability to simulations of future snowfall and temperature variations, to predict the future rate of change of Alaskan glacier mass. Projections of future precipitation and temperature were obtained from the CESM1 Large Ensemble Community Project (Kay and others, 2015). The CESM-LE ensemble includes 30 core members covering 1920–2100, using a fully-coupled model CESM1 (CAM5). The spatial resolution of each simulation is 1° latitude/longitude. For the period 1920–2005, historical radiative forcing due to solar variability, volcanic aerosol emissions and greenhouse gas concentrations was used, and for the period 2006–2100 the Representative Concentration Pathway (RCP) 8.5 radiative forcing scenario was used.

Ensemble members use identical model configurations, except that the initial atmospheric state is perturbed. The CESM-LE ensemble therefore provides estimates of both internal climate variability and the response to external forcings.

We apply the analysis described below to each of the 30 core model projections. The models are identified within the CESM1 Project by assigning them numbers from 001 to 030. In what follows, we first show detailed results for model 001. Later, we will summarize the results for all models.

We start with the model output for the simulated monthly gridded global temperature and snowfall rate fields, and average those fields over all the grid points that correspond to the GRACE Alaskan mascons shown in Figure 1. The resulting averages are shown in Figure 4 (for ensemble member 001). Note that, starting shortly after 2000, the snowfall rate decreases and the temperature increases, with the trends accelerating after about 2050. Figure 5 compares the simulated snowfall rate and temperature results for 2003–2014, with the corresponding results for the ERAI data. The two sets of results are in remarkable agreement, considering that the simulated results did not use real precipitation or temperature data. One implication of this agreement is that the values for K ₀, T ₀ and DDF we obtained by fitting the ERAI data to the GRACE results are likely to be relevant when using the simulated fields.

Fig. 4. Simulation of the monthly snowfall rate and the temperature averaged over the Alaskan glacier mascons, for 1920–2100, from the CESM.

Fig. 5. Results from the CESM atmospheric simulation compared with those from the observation-based ERAI fields, for the snowfall rate and temperature averaged over the Alaskan glacier mascons, between January 2003 and December 2014.

We use the simple model of glacier mass balance described by Eqns (3), (5) and (8), with the parameter values given in Eqn (10), together with the simulated 1920–2100 snowfall rate and temperature values shown in Figure 4, to predict the Alaskan glacier mass through 2100. The time derivative of the glacier mass, dM/dt is shown as the blue curve in Figure 6. Note that dM/dt is of the order of −30 to −40 Gt a⁻¹ until just before 2000, decreases slightly to ~−50 Gt a⁻¹ from 2000 to 2010 (consistent with the August 2002–December 2014 trend of −52 Gt a⁻¹ in the GRACE time series shown in Fig. 2), and then decreases dramatically after that, reaching −220 Gt a⁻¹ by 2100 and showing no signs of leveling off. The undulations in the mass loss rate post 2010 are all due to natural variations projected in the climate simulation used in this case. The model results prior to the GRACE era are at least partially consistent with those derived from observations. Berthier and others (2010), for example, used sequential DEMs to conclude that Alaskan glaciers lost mass at a rate of −42 Gt a⁻¹ between 1962 and 2006.

Fig. 6. The rate of change of total Alaskan glacier mass, from 1920 through 2100, based on the simulated atmospheric data and the model parameters inferred from fitting the ERAI snowfall and temperature fields to the GRACE data. The blue curve shows results where neither the total glaciated area, nor the model parameters T ₀ and K ₀, are allowed to vary as the glacier volume changes. The other curves show results when various combinations of those quantities are allowed to vary. The most realistic result is the black curve, in which all three quantities are allowed to vary. These results assume a lapse rate of −0.98°C (100 m)⁻¹, and a precipitation gradient of 0.08 (100 m)⁻¹, when computing the effects of the variations in T ₀ and K ₀, respectively.

The application of this simple glacier mass model does not account for the fact that the numerical values of the parameters S _gl , K ₀ and T ₀, which we are using to predict future mass variability, are likely to change with time as the total glacier mass decreases. For example, as the glaciers melt, their total area, S _gl is likely to decrease. S _gl is used in Eqns (3) and (5), and so a change in S _gl impacts the predictions of future mass loss. Basically, the smaller the area over which the snowfall or melting rates are integrated, the smaller the total amount of snowfall or melted mass, respectively.

In addition, because the reduction in glacier mass will come preferentially from lower elevations where the temperatures are higher, the average glacier elevation will increase, thus decreasing the average temperature over the glaciers and increasing the average snowfall rate. This will impact the values we use for K ₀ and T ₀; we assume it does not lead to changes in DDF.

In Appendix B, we describe how we compute the change in area caused by a change in glacier mass, and how we incorporate the effects of that change in area back into our predictions of future mass change. We use a simpler approach than used by others (e.g. Radić and Hock, 2011) since we are not modeling on the scale of individual glaciers. The fundamental result is (B17), which says that if $M_{gl}^1 $ is the mass change computed using Eqn (8) under the (incorrect) assumption that the total glacier area, S _gl , does not change with time, and M _gl is the mass change computed after allowing that area to change as the total glacier volume changes, then

(11) $$\Delta {M_{gl}}(t) = M_{gl}^0 \left( {\left[ {1 + 0.264\displaystyle{{\Delta M_{gl}^1 (t)} \over {M_{gl}^0}}} \right] - 1} \right)$$

where $M_{gl}^0 $ is the total initial mass of all the glaciers, which we take to equal to 1.5 × 10⁴ Gt (see (B18)).

The orange curve in Figure 6 (marked as ‘with area feedback’) shows results for dM _gl /dt. For those results, we keep the total glacier area, S _gl , constant until 2010, when we begin applying Eqn (11). This causes the orange and blue curves in Figure 6 to agree up until 2010. Our rationale for doing that is simply that the main objective in this section is to predict future mass loss, and not to speculate on mass loss in the past. Note that the inclusion of the area feedback reduces future mass loss rates considerably. By 2100, the predicted mass change rate is now of the order of −130 to −140 Gt a⁻¹, and appears to have begun to stabilize. We note, however, that this rate is still 2–3 times larger than the present-day mass change rate of ~−50 Gt a⁻¹.

Our method of estimating the impact of the changing glacier geometry on K ₀ and T ₀ is described in Appendix C. Our approach involves determining the increase in mean glacier elevation at each time step, and then using estimates of the atmospheric lapse rate, lr, and the fractional precipitation gradient, d _prec (notation taken from Radić and Hock, 2011), to estimate the resulting changes in K ₀ and T ₀. Those modified values of K ₀ and T ₀ are then used in Eqns (3), (5) and (8) to compute the total mass at each time step.

Note that there is an inconsistency between this approach and the inversion of GRACE data for K ₀, T ₀ and DDF described above. When using the GRACE data, we inverted for just one set of parameters, K ₀, T ₀ and DDF, which we assumed to describe both glaciated and unglaciated land. Therefore the values we obtained are presumably some sort of average of the values for glaciers and the values for unglaciated land. But we are now describing a method of propagating into the future, where we modify those parameters to reflect the change in elevation of the glaciers alone. Presumably there would also be changes in the parameters over unglaciated land, and so to be consistent with the approach we took for the inversion, we should adopt changes in K ₀ and T ₀ that continue to be consistent with averaging the parameters for glaciers and unglaciated land. As described at the end of Section 4, however, we have concluded that the parameter values we obtain from the GRACE inversion are values that mostly describe the glaciers alone; and that the impact of unglaciated land on those values is secondary.

Our method of estimating the values of K ₀ and T ₀ as the glacier area changes, requires estimates of lr and d _prec . For lr, we use the results of Gardner and others (2009), who looked at how the temperature varies with elevation along the surface of four glaciers in the Canadian Arctic. They found lapse rates that varied between lr = −0.31°C (100 m)⁻¹ and lr = −0.64°C (100 m)⁻¹, but that if a single lapse rate is required, lr = −0.49°C New Roman (100 m)⁻¹ would be most appropriate. Although Gardner and others' results were based on observations of glaciers outside of Alaska, we have opted to use lr = −0.49°C (100 m)⁻¹ as our default value, and to use −0.31 and −0.64°C (100 m)⁻¹ as upper and lower bounds on lr when estimating uncertainties. This range of values, incidentally, is similar to Radić and Hock's (2011) conclusion that lr = −0.44 ± 0.18°C (100 m)⁻¹; a range they inferred from tuning a glacier mass model to direct observations from 36 glaciers world-wide.

For d _prec we use Radić and Hock's (2011) result: d _prec  = 0.08 ± 0.05 (100 m)⁻¹. This result came from the same tuning process they employed when deriving their temperature lapse rate. We use d _prec  = 0.08 (100 m)⁻¹ as our default value, and 0.03 and 0.13 (100 m⁻¹) as the limits on d _prec when estimating uncertainties.

We do not modify K ₀ or T ₀ prior to 2010. Figure 6 shows the resulting predictions for the future mass loss rates. To show the impact of allowing K ₀ and/or T ₀ to vary, we include each case as a separate curve in Figure 6. In each case, we include the effects of changes in total glacier area after 2010, as described in Eqn (11). The preferred results are, of course, those when the modifications to both K ₀ and T ₀ are included (black curve).

Note that the changes in both K ₀ and T ₀ cause a reduction in the mass loss rate (both the purple and green curve are less negative than the orange). This occurs because an increase in the mean glacier elevation (due to the mass loss being preferentially concentrated at lower elevations) implies that the mean glacier temperatures decrease, which leads to decreased melting and increased snowfall rates, both of which make dM/dt less negative.

When feedbacks from both K ₀ and T ₀ are included (black curve, Fig. 6), the values of dM/dt decrease more modestly after 2010, reaching ~−90 Gt a⁻¹ by 2100. The rate appears to be starting to level off by then, though it may still be decreasing somewhat. The 2100 rate in this case is still substantially more negative than the 2010 rate of ~−50 Gt a⁻¹. When we integrate the rates shown by the black curve in Figure 6 between 2000 and 2100, we find a total mass loss from Alaskan glaciers of ~6800 Gt, corresponding to 19 mm of global sea-level rise.

Uncertainties. There are several sources of uncertainty on our mass change estimates. These include the following.

1. (1) Errors in the model we use to reproduce the GRACE data. These could include errors in the values of the three parameters used in that model, and also the fact that the model formulation itself is ad hoc, and so cannot be expected to provide a perfect match, no matter what parameter values are adopted. Our model to the GRACE data, for example, still leaves 6% of the time series variance unexplained. Our 1-σ uncertainties for the parameters in Eqn (10) attempt to take into account both the errors in the parameter solutions, and the imperfect nature of the model. Figure 7 shows what happens when the projected mass changes are computed using, not only the default values shown in Eqn (10), but also the 1-σ limits on those parameters. As can be seen, the spread in the results is small: of the order of only 1.3 mm of 21st-century sea-level rise between these curves.

2. (2) Errors in our assumption of how the area-vs-elevation dependence changes as the area decreases. This is discussed in Appendix C (step 2). For our default method, employed in finding the results shown in Figure 6, we use the ‘compression approach’. Figure 8 shows the difference in the mass change projections between using that approach and the ‘truncation approach’. Both these approaches are ad hoc, and unlikely to be correct. But we are assuming that the difference between the results provides some idea of the possible errors caused by using either one of them. Either way the errors appear small, of the order of only 0.4 mm of 21st-century sea-level rise.

3. (3) Errors in the values we have adopted for lr and d _prec . We have recomputed our dM/dt results, using the range of values for lr and d _prec described above, when finding the impact of the decreasing glacier area on K ₀ and T ₀. Specifically, we compute results for lr = −0.31 and −0.64°C (100 m)⁻¹ (instead of the default value lr = −0.49°C (100 m)⁻¹), and for d _prec  = 0.03 and 0.13 (100 m)⁻¹ (instead of the default value d _prec  = 0.08 (100 m)⁻¹). Figure 9 shows results for dM/dt that use those endpoints for these parameters. The black curve in Figure 9 is the same as the black curves in Figures 6–8. The total glacier mass change between 2000 and 2100 for this set of parameter values, ranges from ~5800 Gt (corresponding to a sea-level rise of 16 mm) for lr = −0.63°C (100 m)⁻¹ and d _prec  = 0.13 (100 m)⁻¹, to ~8100 Gt (corresponding to a sea-level rise of 23 mm) for lr = −0.31°C (100 m)⁻¹ and d _prec  = 0.03 (100 m)⁻¹. These results imply that the uncertainty due to errors in our assumed values of lr and d _prec is far larger than the uncertainties caused by errors in the model used to reproduce the GRACE results, and by our imperfect assumption of how the area-vs-elevation curve changes as the glacier area decreases.

4. (4) Uncertainties in the CESM-LE projections of future atmospheric temperatures and snowfall rates. To help quantify the effects of these uncertainties, we have applied our formulation to all 30 core CESM-LE ensemble members, each run from 1920 through 2100. Results for the mass loss rates as functions of time are shown in Figure 10, computed using the model parameters shown in Eqn (10), and assuming the default values lr = −0.44 ± 0.18°C (100 m)⁻¹ and d _prec  = 0.08 (100 m)⁻¹. The black curve in Figure 10 is the same as the black curves in Figures 6–9. The total 2000–2100 mass loss projected from these models varies from 18.2 to 20.6 mm of global sea-level rise. This is considerably smaller than the range of 16–23 mm we obtain using different values of lr and d _prec (point 3).

Fig. 7. The predicted rate of mass change computed using the 1-σ limits of the model parameters given in Eqn (10). The effects of K ₀ and T ₀ feedback are computed using values of the lapse rate and the precipitation gradient: lr = −0.49°C (100 m)⁻¹ and d _prec  = 0.08 (100 m)⁻¹. The black curve shows the same results as the black curve in Figure 6.

Fig. 8. Projections of the mass loss rate for two assumptions of how the hypsometry changes as the glaciated area decreases. The initial values of the model parameters, DDF, K ₀ and T ₀, come from the fit of the model to the GRACE time series, and are given in Eqn (10) The effects of K ₀ and T ₀ feedback are computed using values of the lapse rate and the precipitation gradient: lr = −0.49°C (100 m)⁻¹ and d _prec  = 0.08 (100 m)⁻¹. The black curve shows the same results as the black curve in Figure 6.

Fig. 9. The predicted rate of mass change, computed using different values of the lapse rate and the precipitation gradient when estimating the K ₀ and T ₀ feedback. The initial values of the model parameters, DDF, K ₀ and T ₀, come from the fit of the model to the GRACE time series, and are given in Eqn (10). The black curve shows the same results as the black curve in Figure 6.

Fig. 10. The projected rate of mass change, computed for 30 different core ensemble models for 1920–2100. The black curve is for the first ensemble member (i.e. model 001) and is the same as the black curves in Figure 6–9.

It is apparent from these results that lr and d _prec are the greatest sources of uncertainty in the projections. Adding, in quadrature, the bounds on the 21st-century contributions to sea level from the uncertainty sources (1)–(4), we conclude that the contribution of Alaskan glaciers to global sea level rise during the 21st century, is likely to be 19 ± 4 mm.

6. SUMMARY

We use monthly GRACE gravity field solutions to recover a time series of month-to-month surface mass variations summed over all glaciated areas in southern Alaska from August 2002 through December 2014. We find that the average mass loss rate, over this region and during this time span, was −52 ± 4 Gt a⁻¹. We construct a simple model of ice/snow mass variability for this glaciated region that uses monthly ERAI atmospheric temperature and snowfall fields, along with three adjustable parameters, to estimate monthly variations in accumulation and melting. We adjust the values of those three parameters so that the model output best matches the GRACE results and find that the resulting model can explain 94% of the variance in the GRACE time series.

We then apply that same simple, three-parameter model to 30 simulations of atmospheric temperature and snowfall through 2100, to compute projections of mass loss from Alaskan glaciers. When applying that simple model to the simulated atmospheric data, we use the same values of the three parameters that we inferred from the fit to GRACE for projections through 2010. After 2010 we adjust those values as time progresses, by accounting for how the values of those parameters would change as the total glacier volume and area decrease. Our imperfect knowledge of how those parameters change is the source of greatest uncertainty in our final estimates of future glacier mass loss.

We conclude that the total contribution to global sea-level rise from Alaskan glaciers will be 19 ± 4 mm between 2001 and 2100. This corresponds to a total mass loss of ~7000 Gt, which is nearly 50% of the present-day total Alaskan glacier mass. Note, by comparison Radić and Hock (2011) estimated an Alaskan glacier mass loss during the same time period equivalent to 26 ± 7 mm of sea-level rise (using the A1B emissions scenario). Radić and others, (2014) provided two estimates: one equivalent to 18 ± 7 mm of sea-level rise (RCP4.5), and the other to 25 ± 8 mm of sea-level rise (RCP8.5). Slangen and others (2012) projected the 21st-century sea-level rise to be 27 mm SLE (results as quoted by Giesen and Oerlemans (2013)) estimated 15 mm of sea-level rise between 2012 and 2099 (A1B scenario).

The fact that we use GRACE mass change estimates to calibrate our model is both a strength and a weakness when comparing our approach with the approach used in earlier studies. By using GRACE results, we are able to calibrate our model against estimates of the total mass loss from all Alaskan glaciers, rather than against just a few tens of glaciers, few of which are actually in Alaska. However, because of the limited spatial resolution of GRACE, our mass change model is more ad hoc than models constructed to explain mass loss from individual glaciers. Nonetheless, this model comes remarkably close to other more physically-based projections of Alaska 2100 mass loss.

DEDICATION

We dedicate this paper to John Wahr who passed away before its acceptance. His limitless curiosity and selfless pursuit of understanding the natural world were an inspiration to us all. John's extraordinary contributions to GRACE and the Earth science community have forever changed our field. He will be missed.

Fig. 11. Area-vs-elevation results, based on summing the hypsometry curves of all glaciers that lie within the Alaskan mascons shown in Figure 1. The black curve shows the observed results for the present day. The orange and blue curves show how we have modified those results to reflect future reductions in the total area of the Alaskan glaciers.