An assessment of the disequilibrium of Alaska glaciers

The finite response time of alpine glaciers means that glaciers will be in a state of disequilibrium in the presence of a climate trend. Using a simple model of glacier dynamics, we use metrics of glacier geometry to evaluate the present-day disequilibrium for a population of 5200 alpine glaciers in Alaska. Our results indicate that glaciers throughout the region are in a severe state of disequilibrium. We estimate that the median glacier has only undergone 27 % of the retreat necessary to achieve equilibrium with the present-day climate. In general, glaciers with smaller areas have smaller response times, and so are closer to equilibrium than large glaciers. Because much of Alaska's glacier area is contained in a few large glaciers that are far from equilibrium, and because the rate of warming has increased in the last ∼ 50 years, the median equilibration weighted by area is only 16 %. Our estimates are sensitive to uncertainty in response time and to the shape of the warming trend. Uncertainty is greatest for intermediate glacier response times but is small for glaciers with the smallest and largest response times. Finally, we demonstrate that accounting for the increased rate of warming in the late-20th century is important for estimating glacier disequilibrium, whereas the shape of the warming trend in the early-20th century is less relevant. Our results imply substantial future glacier retreat is already guaranteed regardless of the trajectory of future warming.

1 Introduction

A glacier can be thought of as a reservoir of ice supplied by an input flux of snow accumulation and depleted by an output flux due to the ablation of ice (i.e., all processes that remove mass). Within the glacier, ice flows to redistribute mass from where it accumulates to where it is lost. As such, it takes time (typically decades to centuries) for a glacier's geometry to adjust to a change in climate. Therefore, in the face of a continuous, ongoing warming, glaciers are always in a state of disequilibrium – playing catch up to the evolving climate. Even if the climate were to stop changing today, there would be a period of continued glacier retreat reflecting the remaining adjustment towards a new equilibrium.

The lag of the change in glacier length behind a change in climate was an early target for modern theoretical studies of glacier response (e.g., Nye, 1960, 1965), and also for early numerical modelling studies (e.g., Budd and Jenssen, 1975; Oerlemans, 1986; Huybrechts et al., 1989). Jóhannesson et al. (1989) proposed a simple expression based on glacier geometry that yields response times of several decades for all but the largest alpine glaciers, and subsequent studies have supported this range. Such assessments of glacier response time underly the “very high confidence” of the Intergovernmental Panel on Climate Change (Fox-Kemper et al., 2021) that glacier retreat will continue in coming decades because glaciers are currently still responding to the warming that has already been observed.

The concept of “committed warming” was introduced into climate science as a metric of climate disequilibrium: how much future warming would there be if CO₂ levels were suddenly fixed at today's values (e.g., Hansen et al., 1985; Wetherald et al., 2001)? An alternative definition is the future warming from past human activity (e.g., Armour and Roe, 2011). The first definition is a more direct measure of the dynamical disequilibrium, and we use an equivalent version in this study. The concept of committed future change has been introduced to glaciology (e.g., Dyurgerov et al., 2009; Goldberg et al., 2015; Christian et al., 2018; Marzeion et al., 2018; Zekollari et al., 2020; Hartl et al., 2022): given that glaciers are out of equilibrium with the current climate, what is the “committed retreat” even if the climate were to stop changing today? An assessment of the current state of glacier disequilibrium is important for interpreting and attributing the cause of observed glacier retreat. It is also important for informing initialization strategies in numerical models projecting future glacier change, as initial conditions should be commensurate with the state of disequilibrium.

This study is an extension of Christian et al. (2018) and (2022b). Christian et al. (2018) used a hierarchy of numerical and analytical models to demonstrate that the key factors controlling committed retreat are the glacier response time, the strength of the climate trend, and the glacier geometry. Christian et al. (2022b) applied a simple analytical model to analyse disequilibrium for the glaciers of the Cascades Range in Washington State. Here we evaluate glacier disequilibrium for alpine glaciers in Alaska. Alaska glaciers were the largest regional contribution to sea-level rise from 1961–2016 (Zemp et al., 2019). They are projected to remain the largest regional glacier contribution to sea-level rise through 2100 (Rounce et al., 2023), underscoring the importance of understanding their response to climate change. We take advantage of new datasets and update our methods to better represent the much wider range of glacier geometries in Alaska. We also focus on how the shape of the anthropogenic warming trend impacts glacier disequilibrium. In Sect. 2 we provide a more detailed description of glacier disequilibrium and the metrics we use to characterize it. Our population of Alaska glaciers is described in Sect. 3, and in Sect. 4 our updated methods are applied to estimate the distribution of response times for our glacier population. In Sect. 5 we provide population distributions of the range of disequilibrium that we estimate, and for three different anthropogenic warming scenarios. Given the long response time of some Alaska glaciers, and given an acceleration of anthropogenic warming in recent decades, our method indicates a severe disequilibrium (and so a large committed retreat) for many Alaska glaciers. For example, we estimate that the median observed retreat of the largest glaciers (>250 km²) is only ∼ 15 % of the equilibrium retreat were climate to stop changing today. Finally, we also analyse how errors in our estimated input parameters might propagate to our answers. We conclude with a summary and a discussion of how our estimates might be refined in future work.

2 Glacier disequilibrium

The concepts of glacier disequilibrium and committed retreat are illustrated in Fig. 1 for an idealized industrial-era warming scenario. For the purpose of illustration, we use a simple model (described further below), but the concepts in Fig. 1 are general and not limited to a particular model. Consider the onset of a linear warming trend in 1880. For the sake of simplicity, we assume a constant prior climate and omit interannual variability. At the onset of the trend, we can define a temperature change, T^′(t) (Fig. 1a). One can calculate the glacier-length change, $L_{\mathrm{eq}}^{\prime }\left(t\right)$, that would be in equilibrium with T(t), which is shown by the dashed lines in Fig. 1b. The actual length change, L^′(t), (solid lines in Fig. 1b) is always less than $L_{\mathrm{eq}}^{\prime }\left(t\right)$, because it takes time for a glacier's geometry to adjust to a change in mass balance. Therefore, in the case of an ongoing warming trend, glacier length will always be in a state of disequilibrium with the climate. Let τ be the glacier's response time (defined more precisely below). The larger τ is, the greater will be the degree of disequilibrium (orange vs. blue lines in Fig. 1b). In our synthetic temperature scenario, we choose a year (2020, here) for the warming trend to cease. Disequilibrium will persist, albeit diminishing, as the glacier approaches a new equilibrium state (Fig. 1b).

Figure 1 An illustration of the concepts of glacier disequilibrium and committed retreat. (a) A warming trend applied from 1880 to 2020. (b) Solid lines show the change in glacier length, L^′, for two glacier response times: 25 years (orange) and 75 years (blue); dashed lines show the equilibrium glacier length, $L_{\mathrm{eq}}^{\prime }$, for the climate change up to that point in time. (c) The committed retreat, $\mathrm{|}{L}^{\prime }-L_{\mathrm{eq}}^{\prime }\mathrm{|}$. (d) $L^{\prime }/L_{\mathrm{eq}}^{\prime }$, which we define as the fractional equilibration (f_eq). Values near 0 indicate a glacier far from equilibrium with its climate, and values near 1 indicates a glacier that is near equilibrium with its climate. The shaded region shows the period of warming in all panels. These curves were calculated using Eq. (1), with β=100 (unitless) and $\mathrm{\Delta }{b}^{\prime }\left(t\right)=-\mathit{\mu }\mathrm{\Delta }{T}^{\prime }$, where the melt factor μ=0.65 m yr⁻¹ °C⁻¹.

In the context of a warming climate, disequilibrium can be characterized by the metric of committed retreat, $\left|L^{\prime }-L_{\mathrm{eq}}^{\prime }\right|$, which is the additional retreat after the warming trend ceases. Absolute-value brackets are used to emphasize that it is the magnitude of the difference between L^′ and $L_{\mathrm{eq}}^{\prime }$ we are interested in. Committed retreat is the future retreat that is already guaranteed to occur, even if there is no more warming (e.g., Dyurgerov et al., 2009; Mernild et al., 2013; Goldberg et al., 2015; Christian et al., 2018). At the onset of warming, glaciers are initially slow to adjust but, for a constant warming trend, increase their rate of retreat. In the linear model, the rate of retreat eventually matches the rate of change in $L_{\mathrm{eq}}^{\prime }$ (Fig. 1b), and so approaches a state of constant committed retreat (Fig. 1c, orange line). Over the 140-year trend in our synthetic example, the retreat rate of the smaller-τ glacier fully adjusts to the applied trend (albeit with a persistent lag in absolute length), while the larger-τ glacier is still in the early phase of adjustment and so experiences an accelerating retreat throughout the warming period.

Committed retreat is not the only metric of disequilibrium. Christian et al. (2018) defined the fractional equilibration as $f_{\mathrm{eq}}=L^{\prime }\left(t\right)/L_{\mathrm{eq}}^{\prime }\left(t\right)$. At any given moment in time, it is the ratio of the actual observed retreat to the glacier's eventual retreat if the climate were to stop warming immediately. It is related to the committed retreat, which is equal to $L_{\mathrm{eq}}^{\prime }\left(t\right)\left(\mathrm{1}-f_{\mathrm{eq}}\right)$. Fractional equilibration is a useful measure because it allows the committed retreat to be directly estimated from the observed retreat and, unlike committed retreat, it provides a direct comparison between glaciers with different τ. Moreover, because f_eq is a ratio, it is less sensitive to uncertainty in some relationships, such as that between temperature and melt, which affect both the transient and equilibrium response.

In this study we evaluate these aspects of glacier disequilibrium using the glacier model of Roe and Baker (2014; RB14). RB14 developed a model for glacier-length fluctuations linearized about some prescribed mean-state geometry, and that accurately emulates models of flowline ice dynamics (see Roe and O'Neal, 2009, RB14, Christian et al., 2018). In response to mass-balance anomalies, b^′, length fluctuations are governed by a linear, third-order differential equation:

$$\begin{array}{}\text{(1)}& {\left({\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathit{ϵ}\mathit{\tau }}}\right)}^{\mathrm{3}}{L}^{\prime }={\displaystyle \frac{\mathit{\beta }}{{\mathit{ϵ}}^{\mathrm{3}}{\mathit{\tau }}^{\mathrm{2}}}}{b}^{\prime }\left(t\right),\end{array}$$

where τ is the glacier response time, β is a dimensionless geometric constant that affects glacier sensitivity, and $\mathit{ϵ}=\mathrm{1}/\sqrt{\mathrm{3}}$. The two parameters in Eq. (1) are: $\mathit{\tau }=-H/b_{\mathrm{t}}$ (see Jóhannesson et al., 1989), where H is a characteristic thickness and b_t is the net (negative) mass balance at the terminus; and $\mathit{\beta }=A_{\mathrm{tot}}/wH$, where A_tot is glacier area and w is the width of the terminus zone. The linear assumption means that model parameters are taken to be constant over the length variations of interest. Linear models have been shown to closely emulate numerical flowline models through the course of several kilometers of terminus retreat (e.g., Christian et al., 2018), and are efficient tools for understanding glacier change (e.g., Oerlemans, 2005; Lüthi and Bauder, 2010; Roe et al., 2021; Christian et al., 2022b).

The equilibrium length response to a change in mass balance, Δb, is given by the solution to Eq. (1) when $\mathrm{d}{L}^{\prime }/\mathrm{d}t=\mathrm{0}$:

$$\begin{array}{}\text{(2)}& L_{\mathrm{eq}}^{\prime }=\mathit{\beta }\mathit{\tau }\mathrm{\Delta }b.\end{array}$$

We note that a glacier's transient response to a step-change in climate is sigmoidal rather than exponential in shape (e.g., RB14). It is important not to treat τ as an e-folding timescale. Doing so can significantly underestimate length disequilibrium relative to the true value (Christian et al., 2018).

For a linear climate trend starting at t=0, $b^{\prime }\left(t\right)=\dot{b}t$, where $\dot{b}$ is a constant, and the analytic solution to Eq. (1) is:

$$\begin{array}{}\text{(3)}& L^{\prime }\left(t\right)=\left[\mathrm{1}-{\displaystyle \frac{\mathrm{3}\mathit{ϵ}\mathit{\tau }}{t}}\left(\mathrm{1}-e^{-\frac{t}{\mathit{ϵ}\mathit{\tau }}}\right)+e^{-\frac{t}{\mathit{ϵ}\mathit{\tau }}}\left({\displaystyle \frac{t}{\mathrm{2}\mathit{ϵ}\mathit{\tau }}}+\mathrm{2}\right)\right]\mathit{\beta }\mathit{\tau }\dot{b}t\end{array}$$

(Roe et al., 2017). In the limit of t≫τ, committed retreat converges to $\left|L^{\prime }-L_{\mathrm{eq}}^{\prime }\right|=\mathrm{3}\mathit{ϵ}{\mathit{\tau }}^{\mathrm{2}}\mathit{\beta }\dot{b}$. From this formula we see that $\left|L^{\prime }-L_{\mathrm{eq}}^{\prime }\right|$ is greatest for larger τ, and for larger $\dot{b}$. For t>3τ, the disequilibrium is within 5 % of the asymptotic limit (e.g., Fig. 1c, orange line). In Sect. 4.3, we estimate that the area-weighted median τ=50 years for our population of Alaska glaciers, and so the implication is that many of Alaska's larger glaciers are currently far from the asymptotic disequilibrium limit, and are thus still in their early, accelerating phase of adjustment. For our linear climate trend, we have $L_{\mathrm{eq}}^{\prime }\left(t\right)=\mathit{\beta }\mathit{\tau }\dot{b}t$ which, with Eq. (3), gives:

$$\begin{array}{}\text{(4)}& f_{\mathrm{eq}}\left(t\right)=\mathrm{1}-{\displaystyle \frac{\mathrm{3}\mathit{ϵ}\mathit{\tau }}{t}}\left(\mathrm{1}-e^{-\frac{t}{\mathit{ϵ}\mathit{\tau }}}\right)+e^{-\frac{t}{\mathit{ϵ}\mathit{\tau }}}\left({\displaystyle \frac{t}{\mathrm{2}\mathit{ϵ}\mathit{\tau }}}+\mathrm{2}\right)\end{array}$$

We note that f_eq depends only on τ and the duration of the trend: because f_eq is a ratio, both β and $\dot{b}$ cancel. We note that the cancelation of β occurs even when f_eq is calculated for a trend that is not constant because of the linear nature of Eq. (1). We evaluate the impact of other trend shapes on f_eq in Sect. 5.2. The dependence of f_eq on only one parameter, τ, makes it useful for assessing disequilibrium across a population of glaciers because it reduces the number of uncertain inputs.

Figure 1d shows f_eq for our synthetic linear trend. It remains less than 0.1 for t<τ, and increases most rapidly between τ and 3τ, as the glacier converges to a constant value of $\left|L^{\prime }-L_{\mathrm{eq}}^{\prime }\right|$. Thereafter, f_eq continues to increase, but rather more slowly, as the absolute disequilibrium becomes a progressively smaller fraction of the total retreat. Figure 1 illustrates the importance of τ in determining future glacier behavior. The absolute retreats of our two glaciers are quite similar in 2020, but their different response times (τ=25, 75 years) mean there are large differences in committed retreat (∼ 3 km) and fractional equilibration (∼ 40 %): the glaciers are in different phases of their adjustment.

Our idealized climate scenario in Fig. 1 has neglected interannual climate variability – the year-to-year fluctuations in weather that occur even in a statistically constant climate. Natural interannual variability means that a glacier's length at any instant will rarely be in equilibrium with the long-term climate average. However, in this study we focus solely on the disequilibrium associated with the anthropogenic warming trend. One reason for this is that Roe et al. (2021) considered a wide range of synthetic, modelled, and reconstructed climate scenarios that included natural variability, and showed that post-industrial glacier disequilibrium is overwhelmingly associated with anthropogenic forcing. Secondly, glacier disequilibrium due to natural variability will fluctuate around zero.

3 Glacier inventory

We use the Randolph Glacier Inventory v6.0 (RGI Consortium, 2017), which provides metrics of glacier geometry, including length, maximum elevation, minimum elevation, area, and hypsometry. RGI calculates these statistics from glacier outlines, which reflect the glacier's state at a particular point in time. For RGI v6.0, outlines for Alaska were collected between 2004 and 2010. In this work we exclude: (i) tidewater glaciers, (ii) glaciers with reported areas less than 1 km², and (iii) glaciers with less than 250 m of elevation change. Tidewater glaciers have complex terminus dynamics that can complicate their response to climate trends. Smaller, flatter glaciers have lower driving stress, making their dynamics potentially less-well governed by the ice deformation our model was developed to emulate (e.g., Sanders et al., 2010). Additionally, century-scale changes for smaller glaciers would more likely violate our model's linear assumptions, though their response to climate trends is certainly deserving of study. Of the 27 108 glaciers in Alaska delineated by RGI, 5681 have areas greater than 1 km². Of the glacier area we exclude, 11 600 km² is from marine-terminating glaciers, 6400 km² is from glaciers with area less than 1 km², and 27 km² from the remaining glaciers that have less than 250 m of elevation change. Future work might further consider debris cover as an additional criterion for filtering the glacier inventory. However, after applying our stated selection criteria, 5169 glaciers remain for this analysis, which represents 77 % of Alaska's glacier area. Hereafter, our analysis refers to this subset of glaciers.

Figure 2a shows the geographic distribution of this subset, subdivided into the six glacierized second-order regions classified within the RGI. The histogram of glacier area shows that this population of glaciers is heavily skewed towards smaller glaciers (Fig. 2b). A small number of large glaciers contribute much of the total area: half the total area comes from the largest 1.5 % of the glaciers. In our analysis below we choose to present both probability density functions (PDFs) and cumulative distribution functions (CDFs) plotted as histograms. While the PDFs and the CDFs contain the same information, we find it useful to see both visual representations. We also present these distributions as both number-weighted (counting each glacier equally), and also area-weighted (weighting the contribution of each glacier by its area). Area-weighted distributions give greater emphasis to the larger glaciers in which most of the glacier mass resides, making it more relevant for an assessment of the aggregate mass of Alaska glaciers. For the PDF of a general glacier property, x, the area-weighted count, n^w, for the histogram bin with value x_i, is given by $n^{\mathrm{w}}\left(x_i\right)=\left(\mathrm{1}/A_{\mathrm{tot}}\right){\sum }_j{n}_j\left(x_i\right)A_j$, where the index j refers to all glaciers, with area A_j, that have value x_i, and A_tot is the total glacier area. The weights given to individual area bins can thus be seen in Fig. 2c, which shows how the total glacierized area is distributed across the range of area bins. The larger glaciers contribute strongly to the area-weighted distributions, even though their number is relatively few. In the Appendix, we provide more details of how our analyses vary with glacier size by defining five different area categories (Table A1), and across the six RGI second-order regions (Table A2).

Figure 2 (a) Map of glacier area in Alaska, for glaciers analysed in this study. The dot color shows the total glacierized area within each grid cell, as documented in the RGI v6 (RGI Consortium, 2017). Glacier area was aggregated using a hexagonal grid (displayed as circular markers for visual clarity) with each marker positioned at the center coordinates of its corresponding grid cell. Diamond markers indicate the location of long-term mass-balance records. The inset map shows the boundaries of second-order regions within Alaska as classified by the RGI: the Brooks Range (brown), Alaska Range (pink), Aleutian Range (red), Western Chugach Range (green), Wrangell and St. Elias Ranges (orange), and Coast Range (blue); (b) histogram of glacier numbers, by area; and (c) Distributions of total glacier area within each area bin, which indicates relative weights to each bin when results are area weighted.

4 Estimating glacier response time

4.1 Characteristic thickness, H

The ice thickness in the lower portion of the glacier is the most relevant for the dynamics governing advance and retreat (e.g., RB14). Christian et al. (2022b) estimates thickness for Cascadian glaciers using the method described in Haeberli and Hoelzle (1995) that assumes a uniform slab geometry and a critical shear stress. Given Alaska's larger glaciers, we here choose a method that can better represent thicknesses in the lower portion of the glacier.

We use the recently developed thickness dataset from Millan et al. (2022), which estimated thickness based on satellite observations of ice velocity. We choose Millan et al. (2022) for our main analyses because they report that their method exhibits less bias than the dataset of Farinotti et al. (2019) when compared to observed glacier thickness in Alaska. We report the impact of using other thickness datasets on our analyses in Sect. 5.3. We take H to be the mean thickness along the central flowline below the glacier's equilibrium-line altitude (ELA), with the flowline provided from the Open Global Glacier Model (OGGM; Maussion et al., 2019). The glacier's ELA is estimated using RGI hypsometry and an assumption that the accumulation-area ratio (AAR) is 0.6, which McGrath et al. (2017) found to be a plausible median for Alaska glaciers. Although the present-day AAR is estimated to be lower (e.g. Zeller et al., 2025), this is likely to be, at least in part, a transient consequence of the current disequilibrium. We analyse the sensitivity of our results to these assumptions in Sect. 5.3.

The PDF and CDF of H are presented in Fig. 3a and b, weighted by number (blue) and by area (orange). Weighted by number, the distribution of H is characterized by a median (and 90 % range) of 81 (42, 230) m. Weighted by area, the equivalent distribution has a median (and 90 % range) of 331 (70, 644) m, which shows the substantial skew towards larger glaciers. Table A1 shows a further breakdown of H by area category. In general, the smallest values of H come from small glaciers on steep slopes such as Williams Glacier (H=36 m) and Byron Glacier (H=56 m). The largest values of H are associated with valley glaciers that have termini on a shallow slope. Examples include Tana Glacier (H=700 m), Gilkey Glacier (H=473 m), and Kahiltna Glacier (H=448 m).

Figure 3 (a) PDF and (b) CDF of estimated characteristic thickness (H) using the method outlined in Sect. 4.1. Blue lines are number-weighted and orange lines are area-weighted. Lighter vertical lines indicate the medians of each distribution.

Glacier thickness remains a major source of uncertainty in glaciology. Farinotti et al. (2017, 2019, 2021) report substantial variations among different estimation methods, with an approximate standard error of 50 % when compared with observations (Farinotti et al., 2019). Nonetheless, for the purposes of evaluating a large network of glaciers, these geometry-based estimates capture the distribution of characteristic thicknesses well enough to be applied in our analysis.

4.2 Terminus balance rate, b_t

We calculate a glacier's terminus mass balance rate, b_t, by assuming a vertical mass-balance gradient, $\mathrm{d}b/\mathrm{d}z$, that is constant over the ablation area delineated by its estimated ELA. We estimate a representative $\mathrm{d}b/\mathrm{d}z$ for our population of glaciers from a compilation of point mass-balance observations of Alaska glaciers. We include glaciers with observations spanning a sufficient elevation range to evaluate a balance profile, with particular emphasis on sufficient data in the ablation zone.

The compilation includes point mass-balance records from glaciers monitored by the USGS Benchmark Glacier Project, the Juneau Icefield Research Program, and other publications on selected well-observed glaciers (see Table B1). The locations of these glaciers are marked in Fig. 2a. While Alaska has relatively few mass-balance observations relative to its number of glaciers, each of the state's glacierized regions from Fig. 2 are represented in our compilation. Figure 4 shows the vertical profile of mass balance of each glacier averaged over the available record, and demonstrates that our assumption of a linear balance profile is a good approximation over the ablation zone. On the basis of this compilation, we adopt a value for $\mathrm{d}b/\mathrm{d}z$ of −6.5 m yr⁻¹ km⁻¹ (black dashed line in Fig. 4) as a reasonable central estimate for the region. The calculated $\mathrm{d}b/\mathrm{d}z$ values and details about the data evaluated for each glacier can be found in Table B2.

Figure 4 Observations of the change in annual mass balance with elevation ($\mathrm{d}b/\mathrm{d}z$). Points are the average values of long-term measurement sites. Colored lines are the least-squares fit. The value of $\mathrm{d}b/\mathrm{d}z$ used for this analysis is shown by the dashed black line with a slope of 6.5 m yr⁻¹ km⁻¹. Data for McCall glacier is the mean profile combining intensive survey years 1969–1972 and 1993–1996 (Rabus and Echelmeyer, 1998). Marine-terminating glaciers are included here for context, though not in the rest of the analysis. Data sources are listed in Table B1.

Mass-balance profiles are known to differ by region due to regional differences in climatology, and among individual glaciers due to local topography and geometry (Oerlemans and Hoogendoorn, 1989; Kaser, 2001; Benn and Lehmkuhl, 2000; Larsen et al., 2015). However, we did not find statistically significant relationships between mass balance and continentality, slope, aspect, or latitude; and so do not adjust $\mathrm{d}b/\mathrm{d}z$ based on these factors. Although such relationships have been observed and have a physical basis (Machguth et al., 2006; McGrath et al., 2015; Olson and Rupper, 2019; McNeil et al., 2020; Florentine et al., 2020), quantifying them with confidence for all of Alaska is precluded by the small number of glaciers with mass-balance records. Further work on this topic is an important target of future research.

The PDF and CDF of calculated values for b_t are presented in Fig. 5a and b, weighted by number (blue) and by area (orange). Weighted by number, the distribution of b_t has a median (and 90 % range) of −2.4 (−6.0, −1.0) m yr⁻¹. Weighted by area, the equivalent distribution has a median (and 90 % range) of −5.7 (−10.2, −1.9) m yr⁻¹. The difference in shape between the area-weighted and unweighted CDFs results from larger glaciers typically spanning a larger elevation range.

Figure 5 (a) PDF and (b) CDF of estimated terminus balance rates (b_t). Blue lines are number-weighted and orange lines are area-weighted. Lighter vertical lines indicate the medians of each distribution. Like H, the area-weighted distribution is skewed towards larger (more negative) values.

The most negative values of b_t belong to glaciers with the largest elevation ranges, dz (Table A1). For example, Logan Glacier ($b_{\mathrm{t}}=-\mathrm{9.4}$ m yr⁻¹; dz=5145 m), Fairweather Glacier ($b_{\mathrm{t}}=-\mathrm{8.8}$ m yr⁻¹; dz=4675 m), and Matanuska Glacier ($b_{\mathrm{t}}=-\mathrm{8.7}$ m yr⁻¹; dz=3408 m). Glaciers with elevation ranges greater than 2500 m can be found in the Wrangell-St. Elias Mountains, Alaska Range, and Chugach Mountains. Large maritime glaciers from the icefields of the Coast Range in Southeast Alaska span elevation ranges less than 2500 m but extend down to sea level and have similarly large magnitudes of b_t (e.g., Meade Glacier: $b_{\mathrm{t}}=-\mathrm{8.8}$ m yr⁻¹; dz=2193 m). Glaciers with the least negative values of b_t tend to be found on the dry side of rain shadows (usually north-facing) and throughout the Brooks Range, where lower accumulation rates constrain a glacier's downslope extent. Examples include the north-facing Chamberlin Glacier ($b_{\mathrm{t}}=-\mathrm{2.5}$ m yr⁻¹) in the Brooks Range and Raven Glacier ($b_{\mathrm{t}}=-\mathrm{2.8}$ m yr⁻¹) in the western foothills of the Chugach. Less negative values of b_t also reflect glaciers which have already retreated upslope such as Ptarmigan Glacier ($b_{\mathrm{t}}=-\mathrm{1.6}$ m yr⁻¹) and Flute Glacier ($b_{\mathrm{t}}=-\mathrm{1.7}$ m yr⁻¹).

4.3 Glacier response time, τ

We calculate τ for each glacier using our estimates of H and b_t. Figure 6 shows the PDF and CDF of τ. Weighted by number, τ has a median (and 90 % range) of 35 (14, 96) years. Weighted by area, the distribution has a median (and 90 % range) of 52 (19, 167) years. In general, larger area is associated with both larger H and larger b_t (Table A1). These have offsetting tendencies on τ, so the range of estimated τ is narrower than might be expected given the order-of-magnitude variations in both H and b_t. However, the fractional changes in Hare larger than the fractional changes in b_t, so the overall effect is a tendency for τ to increase with glacier area (Table A1). Although the variance among individual glaciers is high, a linear regression of log (Area) vs. log (τ) shows for every doubling of area, τ increases by a factor of 1.15. (Fig. C1; R²=0.09) for glaciers with areas greater than 5 km².

Figure 6 (a) PDF and (b) CDF of estimated glacier response times (τ). Blue lines are number-weighted and orange lines are area-weighted. Lighter vertical lines indicate the medians of each distribution. While there are many glaciers with short response times, they are generally small. The longer response times of larger glaciers means that over half of all Alaska glacier area we analysed is contained in glaciers with response times greater than 50 years.

In practice, glaciers with the shortest τ (<15 years) tend to be hanging glaciers and small valley glaciers, generally <10 km², where steep slopes result in thinner ice and higher values of b_t. Examples include the smaller glaciers of College Fjord (Vassar Glacier, τ=16 years; Holyoke Glacier, τ=10 years), Cantwell Glacier (τ=16 years), and Byron Glacier (τ=15 years). However, small glaciers that widen below their ELA or terminate on shallow slopes (e.g., cirque glaciers) can occasionally have τ exceeding 50 years (Laughton Glacier, τ=57 years; Anderson Glacier, τ=82 years; see also Barth et al., 2018). Glaciers with the longest response times (>100 years) tend to be large and also tend to terminate on shallow slopes. Some examples of glaciers with larger τ include Tweedsmuir Glacier (τ = 91 years), Tana Glacier (τ=132 years), and Brady Glacier (τ = 140 years).

We note that τ has been estimated using the modern, rather than the preindustrial, glacier geometry. For most glaciers this is a small effect relative to the other uncertainties that impact τ. The linear model closely emulates numerical flowline models through the course of several kilometers of terminus retreat (e.g., Christian et al., 2018) suggesting that the effective response time does not change rapidly with glacier state. The biggest errors are likely to be for the smallest glaciers that have undergone the largest fractional change in their geometry over the industrial era. We discuss how uncertainties in τ affect our results in Sect. 5.3.

5 Estimating glacier disequilibrium

5.1 Disequilibrium for a linear warming trend

We first estimate glacier disequilibrium for a linear warming trend. We choose the year 1880 as the onset of the trend, consistent with the IPCC (Eyring et al., 2021) and related prior work on glaciers (Roe et al., 2017, 2021; Christian et al., 2022b). We analyse glacier response through 2020 (i.e., t=140 years in Eq. 4). Figure 7 shows the PDF and CDF of fractional equilibration (f_eq). Recall that $f_{\mathrm{eq}}=L^{\prime }/L_{\mathrm{eq}}^{\prime }$, and that values near zero indicate a large disequilibrium, and values near one indicate near equilibrium. Weighted by number, the distribution of f_eq has a median (and 90 % range) of 0.57 (0.17, 0.83). A fractional equilibration of f_eq=0.57 means the median glacier still has 43 % of its total retreat ($L_{\mathrm{eq}}^{\prime }$) remaining to reach equilibrium with the present-day climate. Weighted by area, the distribution has a median (and 90 % range) of 0.41 (0.06, 0.76). The opposing skew in the PDFs between the number-weighted and area-weighted distributions mirrors that of glacier number and area (Fig. 2b, c). That is to say, smaller glaciers, of which there are many, tend to have small response times and are generally closer to equilibrium. Larger glaciers, which together constitute most of Alaska's glacierized area, are generally further from equilibrium, reflecting their larger response times (Table A1). For the largest category of glaciers we consider, the median fractional equilibration is f_eq=0.37 implying, for these assumptions, that a committed retreat (i.e., $\mathrm{|}{L}_{\mathrm{eq}}^{\prime }-L^{\prime }\mathrm{|}$) nearly twice the current retreat (L^′) is already built into the glacier's future response.

Figure 7 (a) PDF and (b) CDF of estimated fractional equilibration (f_eq) for a 140-year linear climate trend. Blue lines are number-weighted and orange lines are area-weighted. Low values of f_eq indicate more severe disequilibrium, high values indicate near-equilibrium states. Lighter vertical lines indicate the medians of each distribution.

5.2 Shape of the warming trend

Up to this point, we have approximated the industrial-era warming trend as linear so that f_eq has an analytical solution (Eq. 4). In this section, we evaluate the degree to which the shape of the anthropogenic climate forcing over time can influence glacier disequilibrium. We consider two alternatives to the linear warming trend, which represent alternative views on the shape of the anthropogenic warming.

The dotted grey line in Fig. 8a shows Alaska-wide summertime (JJAS) temperature anomalies since 1880, taken from the Berkeley Earth dataset (Rohde and Hausfather, 2020), and the solid grey line applies a 30-year low-pass filter. While there is substantial year-to-year variability, there has been an overall warming of ∼ 1.2 °C since 1880. This is a typical value for summertime trends seen elsewhere on land (e.g., Allen et al., 2018). The blue line is the linear warming scenario, for which we presented results in Fig. 7. The red line is the second warming scenario. It is characterized by a warming in the first half of the 20th century, a mid-century cooling, and a resumed warming in recent decades. This scenario incorporates what is sometimes deemed the “Early Twentieth-Century Warming” (hereinafter ETCW; Hegerl et al., 2018), which characterizes many, particularly high-latitude, temperature records. It is debated whether this shape is attributable to anthropogenic causes or reflects natural, internal climate variability (Haustein et al., 2019). The purple line is the third warming scenario. Its shape is an approximation of the Global Warming Index (hereinafter GWI), which is an estimate of the global anthropogenic radiative forcing (Haustein et al., 2017). The change in the rate of warming circa 1970 reflects an increase in the rate of industrial emissions. All three scenarios are standardized to start at 0 °C in 1880 and reach a value of 1.2 °C in 2020. We note that fractional equilibration (f_eq) here depends only on the shape of the trend and not its magnitude; because the bulk of the warming has occurred in the past 50 to 100 years, the results are not sensitive to the choice of 1880 as a starting point. In this study, we consider only temperature trends. Observed annual precipitation trends in Alaska have generally been small over the 20th century (McAfee et al., 2013; Bieniek et al., 2014), and can be difficult to quantify across the whole region due to sparse early data (McAfee et al., 2014; Ballinger et al., 2023).

Figure 8 Glacier disequilibrium for three warming scenarios: (a) Observed annual Alaska-wide summertime temperature anomalies (JJAS) from the Berkeley Earth dataset (dotted grey), and with a 30-year low-pass filter (zero-phase second-order Butterworth, solid grey). Three warming are considered: linear (blue), ETCW (red) and GWI (purple). Observations and scenarios are set at 0 °C in 1880; (b) number-weighted; and (c) area-weighted f_eq calculated at 2020 for each warming trend. The blue lines are the same as Fig. 7b. Lighter vertical lines indicate median values. The population medians and medians for different area categories defined by glacier area are reported in Table A1.

To assess f_eq for our three warming scenarios, we calculate L^′(t) by integrating Eq. (1) forward in time. We use our estimates of τ (Fig. 6), and assume that temperature anomalies, T^′(t), are related to mass balance by $b^{\prime }\left(t\right)=-\mathit{\mu }{T}^{\prime }\left(t\right)$, where μ is a constant melt factor (units of m yr⁻¹ °C⁻¹), and we have ${L^{\prime }}_{\mathrm{eq}}\left(t\right)=-\mathit{\beta }\mathit{\tau }\mathit{\mu }\cdot T^{\prime }\left(t\right)$. Finally, we have $f_{\mathrm{eq}}=L^{\prime }/L_{\mathrm{eq}}^{\prime }$. Note that when the ratio of lengths is taken, both μ and β cancel out, so our results do not depend on their values.

The effect of these three warming scenarios on f_eq were calculated at 2020 (Fig. 8b, c). Relative to the linear scenario, both the ETCW and GWI scenarios show a greater degree of disequilibrium; their greater rate of warming in recent decades means that glaciers are further from equilibrium with the current climate, thus lowering f_eq. For the ETCW scenario, weighted by number, the distribution of f_eq has a median (and 90 % range) of 0.29 (0.09, 0.59). Weighted by area, the distribution has a median (and 90 % range) of 0.20 (0.03, 0.48). For the GWI scenario, weighted by number, the distribution has a median (and 90 % range) of 0.27 (0.05, 0.62). Weighted by area, the distribution has a median (and 90 % range) of 0.16 (0.02, 0.51).

Our results demonstrate that the shape of the warming has a substantial impact on f_eq, for both number- and area-weighted distributions. The median values of f_eq are approximately halved in the ETCW and GWI scenarios compared to the linear scenario. Expressed another way: for the linear scenario, approximately 60 % of glaciers have f_eq >0.5 (i.e., are more than halfway equilibrated to the current climate); whereas for the ETCW and GWI scenarios, this falls to only 12 % of glaciers (Fig. 8b). The differences among scenarios are somewhat less for the area-weighted distributions because there is more weighting given to the larger τ glaciers, for which the decadal-scale details of the warming matter less. It is interesting to note that the ETCW and GWI scenarios have similar f_eq distributions (Fig. 8b, c). This is because we calculated f_eq for the year 2020, for which the previous 50 years of climate history are similar in both scenarios (Fig. 8a). In the ETCW scenario, f_eq calculated at 1960 can exceed one for glaciers whose retreat at that point exceeds the long-term equilibrium changes after a period of forced cooling.

5.3 Sensitivity of the results

The two factors that control fractional equilibration (f_eq) are the response time (τ) and the shape of the anthropogenic-warming scenario, and both have some uncertainty. In order to assess the relative importance of these uncertainties, we can plot f_eq (at 2020) as a function of τ and our three warming scenarios (Fig. 9). All three scenarios converge on the limits of f_eq→1 for small response times (τ→0 years) and of f_eq→0 for large response times (τ≳150 years). The largest difference between the linear and ETCW/GWI scenarios occurs for τ ∼ 50 years, for which the differences in recent warming have the biggest impact. This τ is close to the area-weighted median we found for our whole population (τ=52 years; Table A1). Figure 9 can be used to evaluate how uncertainties in τ lead to uncertainties in f_eq. The information is not available to perform a formal uncertainty analysis, so in previous work (Roe et al., 2017) we assigned broad uncertainty to our estimates of τ. We assumed a standard-deviation of $\mathit{\tau }/\mathrm{4}$, meaning the 95 % bounds (i.e., ±2σ) on the range for τ is equal to τ itself. In other words, the 95 % bounds for τ spans a factor of three ($\mathit{\tau }/\mathrm{2}$ to $\mathrm{3}\mathit{\tau }/\mathrm{2}$). We consider this to be a reasonable estimate for this study too. Applying this to Fig. 9 and the linear scenario, the largest uncertainty in f_eq arises for glaciers with τ ∼ 20–60 years. For such glaciers, uncertainty in τ and choosing a realistic warming scenario over the linear trend have comparable importance. When choosing between realistic scenarios however, the uncertainty is small compared to that of τ. Note that, for most values of τ, the symmetric uncertainty we assess around τ translates to asymmetric uncertainty around f_eq. For glaciers with long response times (e.g., τ≳150 years), uncertainties in f_eq are low because f_eq ∼ 0 across a wide range of τ regardless of warming scenario. For glaciers with small response times (τ→0 years), the absolute uncertainty in τ, and thus the uncertainty in f_eq, is lower than for the intermediate values of τ, despite the steeper slopes of the curves in Fig. 9 for smaller τ. In summary, for any individual glacier, especially those with τ of a few decades, the value of f_eq can have large uncertainty because of uncertainty in both τ and the true shape of the warming trend.

Figure 9 Fractional equilibration in 2020, as a function of response time for the three warming scenarios, using the method described in Sect. 5.2. For any given τ, we suggest a reasonable uncertainty range (at ±2σ) is $\mathit{\tau }/\mathrm{2}$ to $\mathrm{3}\mathit{\tau }/\mathrm{2}$. The shaded area and the colored bars on the right-hand side illustrate how, for a central estimate of τ=50  years, uncertainty in τ (grey bar) affects uncertainty in f_eq for the ETCW and linear scenarios.

We also tested alternative values of H and b_t in our estimation of τ. First, we repeated our analyses using the thickness dataset from Farinotti et al. (2019), which produces a lower median and narrower distribution of H as compared to Millan et al. (2022). Taking the GWI scenario as an example, using Farinotti et al. to estimate H gives a number-weighted median (and 90 % range) for f_eq of 0.35 (0.10, 0.72), compared to 0.27 (0.05, 0.62) for Millan et al. (2022). When weighted by area however, the difference in the resulting median and 90 % range between datasets is negligible. This contrast between weighting methods is mainly attributable to the largest glaciers having a greater estimated thickness in the Farinotti et al. dataset.

Secondly, we tested the impact of altering the AAR in our calculation of b_t. Zeller et al. (2025) provide AARs for a large population of Alaska glaciers based on satellite imagery of end-of-summer snow cover. They report an average AAR of 0.4 for all glacier area in Alaska over the period 2018–2022. Using this dataset, we examine two variations on our analysis: (i) applying the average AAR of 0.4 uniformly to match our standard analysis, and (ii) using the AARs Zeller et al. derived for individual glaciers. Detailed results and summarized values of b_t, τ, and f_eq are presented in Appendix D (Table D1), using the GWI warming scenario as an example.

In our model, a smaller AAR yields a higher ELA, resulting in more negative b_t and thus a smaller τ. For small glaciers with small magnitudes of b_t, which tend to also have smaller H, even a modest change to b_t causes a proportionally larger change in τ. This corresponds to greater sensitivity for small τ to alternative values of AAR. Conversely, because large glaciers tend to have large τ (Fig. C1), we find that our estimates of area-weighted f_eq remain relatively insensitive to AAR assumptions. Overall, these sensitivity studies show our main findings remain consistent to alternative dataset choices and the analyses can readily be updated as new datasets and better observations are obtained.

6 Summary & Discussion

In this study we have estimated the state of disequilibrium for a large population of Alaska glaciers. We excluded small glaciers and tidewater glaciers because our model of glacier dynamics is less applicable to such systems, leaving us with a population of approximately 5600 glaciers representing 79 % of the region's glacier area. There are considerable observational uncertainties in glacier thickness and mass balance, which means that estimating disequilibrium for any single glacier involves substantial uncertainty. By analysing a whole population of glaciers, we aim to provide a region-wide picture, without depending on the uncertain details of one particular setting.

Building on previous work (Christian et al., 2022b), we define a metric of disequilibrium that uses a simple linear model of ice dynamics; and uses characteristic ice thickness, H, and terminus mass balance, b_t, to estimate the glacier response time, τ. A key advantage of our metric is that, because it is a ratio, other geometric parameters cancel out. One disadvantage of using our linear model is that parameters are fixed, so the predicted retreat does not account for settings where the glacier geometry (i.e., width, bed slope) varies strongly upslope of the current margin. More sophisticated and detailed numerical modelling might address such settings. However, given the substantial uncertainties in input parameters and initial conditions, and given our focus on region-wide disequilibrium, we do not think that a more complex modelling system would necessarily provide better assessments. Ice thickness and mass-balance gradients still must be estimated in the absence of direct observations, and biases in these estimates will still propagate into the results of a model with more sophisticated ice dynamics.

We provide results as statistical distributions for our selected population of glaciers, both number weighted (i.e., by count), and area weighted. Table A1 presents our results in separate categories defined by glacier area. Table A2 shows the same results but instead categorized by RGI second-order region. A reader can get a sense of the geographic variations of our results by referencing Tables A1 and A2, and the distribution of areas in the different glacierized second-order regions shown in Fig. 2. Table A2 shows how the varying glacier geometries within second-order regions aggregate to produce distributions of response times and disequilibrium. However, for the full Alaska region and the second-order regions, our estimated distributions of H, b_t, and τ appear reasonable relative to published estimates for individual glaciers. We defined fractional equilibration, f_eq, as the ratio of current retreat to the eventual equilibrium retreat if the climate were to immediately stop changing. Our analyses identify two key factors affecting f_eq: glacier response time and the shape of the anthropogenic forcing, and we also demonstrate that there is interplay between them. Small-τ glaciers (τ →0 years) are in near-equilibrium with the current climate (f_eq →1), and so do not depend sensitively on climate history. Large-τglaciers (τ≳100 years) are far from equilibrium (f_eq → 0), and so respond slowly enough that decadal-scale details of past climate trajectory also do not matter. However, for intermediate-τ glaciers (representing the majority of glaciers in our population), f_eq varies sensitively and depends on both the shape of the forcing and on τ.

The linear warming trend has an analytic solution for f_eq, which is attractive as a first approximation and provides insight into the physical factors affecting how f_eq evolves over time. However, it is likely that the ETCW and GWI scenarios are better representations of the true anthropogenic forcing, both of which imply that glaciers are further from equilibrium than in the linear scenario. These two scenarios have similar-enough climate histories over the past five decades that the values of f_eq are also similar (Figs. 8, 9). Taking just the GWI scenario, when weighted by number, we estimate the distribution of f_eq has a median (and 90 % range) of 0.27 (0.05, 0.62). Weighted by area, the distribution has a median (and 90 % range) of 0.16 (0.02, 0.51). These results indicate that, overall, Alaska glaciers are in a state of dramatic disequilibrium due to anthropogenic climate change, with substantial continued retreat already guaranteed. This degree of disequilibrium also implies that most of the glacier mass in Alaska is in the early, still accelerating, phase of its response to climate change.

For hydrologic-resource planning, or for other local purposes, the state of disequilibrium of one specific glacier may be of importance. Our algorithm is useful for providing an initial estimate, but it should be considered only as a rough approximation. A more comprehensive assessment might be made for individual glaciers using more of the available information. For instance, if the vertical profile of the mass balance is known then, for a given catchment geometry, the glacier extent that would be in equilibrium with that profile can be estimated (e.g., Rasmussen and Conway, 2001). One challenge is that mass balance remains poorly observed in most locations. Moreover, due to large interannual variability, it can take many years of observations to accurately assess the long-term position of the ELA, which is all the harder in the face of long-term trends. For specific locations, it is often possible to assess the accuracy of the glacier-thickness estimates using radar measurements. Given the severe disequilibrium indicated in our results, it seems important to assess the impact of varying thickness on glacier simulations. For individual glaciers, it may be that more complex ice-dynamics modelling (2D or 3D) over the realistic catchment geometry offers improved estimates of disequilibrium, although that should always be evaluated relative to the other uncertainties, such as thickness, mass balance, and climate history. Finally, there is an increasing availability and use of remotely sensed imagery (e.g., Maurer et al., 2019; Shean et al., 2020; Friedl et al., 2021; Jakob and Gourmelen, 2023; Knuth et al., 2023). In the face of warming, glacier adjustment takes the form of overlapping sequential stages of (i) thinning, (ii) reduced fluxes, and (iii) retreat (e.g., RB14), all of which can now be observed with remote sensing. It may be such observations can be combined to estimate the dynamical phase of a retreat.

Our algorithm offers a straightforward and efficient method to assess glacier disequilibrium at regional scales. One robust result is that the wide range of glacier response times implies a broad distribution in the degree of glacier disequilibrium within any region. Constraining current disequilibrium is important for the numerical modelling of glacier populations. In particular, care must be taken to properly represent the initial conditions of glacier simulations. It should also be recognized that uncertainty in climate history and glacier response time means significant uncertainty is attendant on disequilibrium estimates (especially for individual glaciers), and likely points to the need for ensemble modelling techniques to provide probabilistic projections of future glacier evolution (e.g., Christian et al., 2022a). Finally, our assessment of Alaska glaciers implies a severe degree of glacier disequilibrium with the current climate. It seems important for policy makers, resource managers, and the general public to appreciate how much retreat is already baked into the future evolution of glaciers.

Appendix A

Table A1 shows summary statistics for our population of glaciers categorized into bins by area. The median value within each bin is reported along with the 90 % range. The median values for the magnitudes of H, b_t, and τ increase with glacier area. However, the 90 % ranges are wide and often overlap the ranges of adjacent area categories, reflecting the wide range of individual glacier geometries in each. As shown by Fig. 2b, c, there are fewer glaciers with large areas. The geometry for these few largest glaciers is much more influential to our area-weighted estimates, highlighting their individual characterization as a key target for improving estimates of disequilibrium.

Table A2 shows the same summary statistics for our population of glaciers as Table A1, but instead categorized by RGI second-order regions, to illustrate geographic variation in glacier characteristics across Alaska. Differences in glacier geometry among subregions reflect varying topographic and climatic conditions. Maritime subregions (Coast Range, Western Chugach Range) generally support some of the largest, thickest glaciers (note the 90 % ranges, rather than medians) and also more negative terminus balance rates compared to continental subregions (Brooks Range, parts of the Alaska Range). The Brooks Range, though having the fewest glaciers in our analysed population (76), exhibits distinctive characteristics with the thickest ice relative to glacier size (median H=108 m) and the longest response times (median τ=54 years). The Coast Range, in a more temperate maritime climate, ranks second regionally in both median ice thickness and response time (H=92 m, τ=39 years). It is important to bear in mind the differences in glacier count and area between these regions (as well as the spread in their 90 % ranges), but the comparisons reveal how different climate and topographic settings can produce relatively similar distributions of response times, and by extension, disequilibrium states.

Table A1 Median calculated values binned by glacier area. Parentheses are the 90 % range. Note that, due to the shape of the distribution of glacier areas in Fig. 2b, the areas of glaciers within each bin are skewed towards the lower bound. The summary statistics for the population are shown in the bottom rows. These values are equivalent to the number-weighted statistics reported in the main text. For comparison, area-weighted statistics for the population are included in the bottom row.

Table A2 Median calculated values binned by RGI second-order region. Parentheses are the 90 % range.

Appendix B

Table B1 presents the Alaska mass-balance data we compiled to produce a representative region-wide value for $\mathrm{d}b/\mathrm{d}z$ in the ablation zone (i.e., Fig. 4). The estimated $\mathrm{d}b/\mathrm{d}z$ values for each glacier and a summary of the data evaluated is shown in Table B1. Because there are few glaciers with long-term records, we include shorter records where there are sufficient measurements in the ablation zone. For glaciers with long records, we find the shape of the balance profile and magnitude of $\mathrm{d}b/\mathrm{d}z$ to be similar among different years.

Table B1 Summary of mass balance data evaluated for Fig. 4. The first and second columns give the slope ($\mathrm{d}b/\mathrm{d}z$) and standard error of the least-squares fit of annual point mass balance on elevation. Mass-balance data were detrended for time. The remaining columns show: the range of years where mass-balance observations were evaluated, the total number of point observations in each regression, the range of elevations among observations, and the glacier's total elevation range for reference.

Appendix C

For valley glacier geometries, intuition suggests that a glacier with greater area is typically also thicker and able to sustain a terminus extent with greater ablation. These tendencies of H and b_t offset each other in our distribution of estimated τ and result in a smaller spread than would be expected for independent variables (which can be seen from the spreads of H, b_t, and τ in the different size categories in Table A1). H and b_t are indeed somewhat anticorrelated for glaciers in our population ($r=-\mathrm{0.52}$). We find a geometric relationship between τ and glacier area in our population for glaciers larger than 5 km². For a linear regression of log (τ) on log (Area), every doubling of area corresponds to a predicted τ that increases by a factor of 1.15. Although it only explains 9 % of the variance in log-log space, the large number of data points mean the relationship is highly robust (Fig. 9). When using the Farinotti et al. (2019) thickness dataset to estimate τ, we find a larger scaling factor of 1.22, with the relationship explaining 22 % of the variance (not shown). For τ derived from the Farinotti et al. dataset, glaciers smaller than 5 km² show a weak positive correlation with area, yielding a scaling factor of 1.08. The relationship explains less than 1 % of the variance but is statistically significant and distinct from that of larger glaciers. No corresponding relationship is observed when using the Millan et al. (2022) dataset.

Figure C1 Log-log relationship between glacier area and response time, τ using the Millan et al. (2022) dataset. Lines show the least-squares linear fit to the log-transformed data for glaciers with area <5 km² (red), and ≥5 km² (black). The black line has the equation: τ=25 years × $(A/A_{\mathrm{0}}{)}^{\mathrm{0.15}}$, where A₀=1 km². A doubling in glacier area is associated with τ increasing by a factor of 1.15 (SE = 0.014; R²=0.09). The red line shows the linear fit for glaciers with area <5 km², where there is no equivalent correlation. Points are scaled by area for visual clarity.

Appendix D

Table D1 shows summary statistics from the alternative AAR cases mentioned in Sect. 5.3, using the GWI warming scenario as an example. The results from the standard analysis are included for comparison, denoted as the AAR_0.6 case here. The median (and 90 % range) of the affected variables are reported for each case, weighted both by number and by area. The difference between the AAR_0.6 and AAR_0.4 cases represents the sensitivity to the AAR value chosen. The difference between the AAR_0.4 and AAR_Z cases demonstrates the impact of assuming a single AAR for all glaciers.

Relative to the standard analysis, the AAR_0.4 case corresponds to a number-weighted (and area-weighted) increase in ELA of 100 (250) m. Despite sharing a median AAR, the number-weighted (and area-weighted) median ELA is an additional 150 (20) m higher in the AAR_Z case than in the AAR_0.4 case, suggesting individual AARs may be particularly relevant for representing small glaciers. For our population, using individual AARs yields significantly more negative estimates of b_t than assuming an equivalent uniform AAR. The size of this effect is comparable to the difference between the AAR_0.4 case and the standard AAR_0.6 case. When binned by glacier area (as in Table A1), the difference in b_t among cases is notably consistent in magnitude. For glaciers with small magnitudes of b_t, which tend to also have smaller H, even a modest change to b_t causes a proportionally larger change in τ. This sensitivity decreases for glaciers with more negative b_t. Figure 9 shows that sensitivity in f_eq is amplified for small τ, and damped for large τ Because large glaciers have long response times, we find that our estimates of area-weighted f_eq are more insensitive to our estimated AAR.

Table D1 Median calculated values for versions of our analysis using alternative AARs. Parentheses are the 90 % range. Values for each variable are given weighted both by number and by area. Results from the standard analysis, denoted as the AAR_0.6 case, are included for comparison with the two alternative cases. The AAR_0.6 case is identical to the main analysis. The AAR_0.4 case applies a uniform AAR of 0.4 for all glaciers. The AAR_Z cases uses individual glacier AARs provided by Zeller et al. (2025) based on satellite imagery of end-of-summer snow cover.