1. Introduction
Projections of glacier mass change estimate that 26–41% of global glacier mass will be lost by 2100, relative to 2015, for global mean temperature increases of 1.5–4°C, relative to preindustrial levels (Rounce and others, Reference Rounce2023). Glaciers in Alaska and northwest Canada (hereafter ‘Alaska’) are currently losing 0.85 m w.e. a−1 (Zemp and others, Reference Zemp2019) and are projected to lose 50–75% of their current mass by 2100 (Rounce and others, Reference Rounce2023). Alaska consequently is the greatest contributor to present-day and future glacier mass loss of all the glacierized regions outside of the ice sheets. Refining estimates of the timing and magnitude of glacier mass loss is thus important for sea-level rise and more localized impacts such as changes to water resources, ecosystems, and culture. Refining estimates of future glacier change is important for long-term planning and hazard assessment efforts across the region, and particularly within Alaska’s national park units, which contain approximately half of the region’s glacierized area (42 774 km2 out of 86 708 km2; Loso and others, Reference Loso, Fahnestock, Larsen and Tober2026).
The calibration and validation of glacier evolution models used for projections rely on glacier change observations. Cloud computing along with recent satellite missions now provide regional- to global-scale estimates of glacier mass change (e.g. The GlaMBIE Team, 2025), though uncertainties remain high at the individual glacier scale (Dussaillant and others, Reference Dussaillant2025). These various observational datasets provide unique opportunities and challenges for calibrating models. For example, elevation change estimates from the digital elevation model (DEM) differencing or altimetry must be aggregated to glacier-wide scales, such that the climatic mass balance (i.e. the sum of accumulation, melt, and refreeze) equals the total mass balance, or the models must properly account for the climatic mass balance, flux divergence, and density variations. Spatial and temporal resolutions are also important considerations. DEM differencing currently provides decadal estimates of elevation change at a spatial resolution of 100 m (Hugonnet and others, Reference Hugonnet2021), whereas the highest temporal resolution altimetry products provide monthly changes but at a spatial resolution of 500 m (Jakob and Gourmelen, Reference Jakob and Gourmelen2023). In contrast, glacier evolution models are often run at a monthly timestep using elevation bins (Zekollari and others, Reference Zekollari2025).
In situ measurements of the climatic mass balance are ideal for calibrating and validating large-scale glacier evolution models; however, long-term (>5 years) measurements of seasonal or annual mass change are only available for ∼300 of the world’s 274000 glaciers (WGMS, 2025). These measurements are also disproportionately concentrated on relatively small glaciers in lower-latitude, more accessible regions; hence, only 2% of these long-term in situ mass-balance measurements are from glaciers in Alaska. Targeted campaigns, such as NASA’s Operation IceBridge, help to address this gap by measuring seasonal and annual changes in surface elevation over multi-year/-decade periods at high spatial resolution (Larsen and others, Reference Larsen, Burgess, Arendt, O’Neel, Johnson and Kienholz2015; MacGregor and others, Reference MacGregor2021). However, these data present similar challenges of requiring models to accurately estimate the climatic mass balance and flux divergence to meaningfully leverage these data for calibration, a central focus of this work.
Most current state-of-the-art large-scale glacier evolution models (Zekollari and others, Reference Zekollari2024, Reference Zekollari2025) calibrate each glacier independently using a single multi-decadal glacier mass change estimate from DEM differencing (Hugonnet and others, Reference Hugonnet2021). The use of one observation to calibrate multiple parameters leads to equifinality issues (Rounce and others, Reference Rounce, Khurana, Short, Hock, Shean and Brinkerhoff2020; Schuster and others, Reference Schuster, Rounce and Maussion2023), which may be reduced by leveraging additional data. Unfortunately, despite recent advances in open-source glacier evolution modeling (Maussion and others, Reference Maussion2019; Rounce and others, Reference Rounce2023), current models are not set up to integrate multiple data types into their calibration scheme (Morlighem and Goldberg, Reference Morlighem, Goldberg, Ismail-zadeh, Castelli, Jones and Sanchez2023). Given the growing abundance of observational datasets, there is an urgent need to develop flexible calibration frameworks that assimilate diverse observations to improve our predictive capabilities of how glaciers will evolve over the next century.
Here we develop a new flexible calibration framework for a large-scale glacier evolution model that leverages spatially distributed elevation change measurements derived from airborne laser altimetry between 1994 and 2021 in addition to glacier-wide mass-balance estimates. Detailed analyses are used to quantify how the integration of the spatially distributed data reduces model parameter uncertainties for the region and improves hindcast simulations of glacier surface elevation change. We then project glacier changes from 2015 to 2100 for the 185 surveyed glaciers (representing 56% of Alaska’s regional glacier area) for various future climate scenarios, and aggregate projections for Alaska’s major glacierized national parks to support park management and planning.
2. Data and methods
2.1. Surface elevation change data
We use spatially distributed surface elevation change observations of glaciers in Alaska (Loso and others, Reference Loso, Fahnestock, Larsen and Tober2026), derived from laser altimetry profiles through NASA’s Operation IceBridge (Larsen and others, Reference Larsen, Burgess, Arendt, O’Neel, Johnson and Kienholz2015; Larsen, Reference Larsen2020) and prior flights (Echelmeyer and others, Reference Echelmeyer1996; Arendt and others, Reference Arendt, Echelmeyer, Harrison, Lingle and Valentine2002; Johnson and others, Reference Johnson, Larsen, Murphy, Arendt and Zirnheld2013), to calibrate our glacier evolution model (see Section 2.2). Altimetry profiles were acquired from 1994 to 2021 for 205 glaciers, typically along the glacier centerline, and repeated for each glacier at irregular intervals for a total of 2–15 (median 5) surveys in the 27 year span. Time series of surface elevation changes from these data are estimated through a novel differencing approach that accounts for spatial variations associated with potentially non-coincident flight lines (Loso and others, Reference Loso, Fahnestock, Larsen and Tober2026). Specifically, offsets of altimetry data from a reference DEM were first calculated for each survey and aggregated into 10 m elevation bins, using Randolph Glacier Inventory version 7.0 (RGI Consortium, 2023) glacier outlines. These binned offsets from Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026) comprise the fundamental measurements of elevation change we use in this paper.
Binned surface elevation offsets were mapped to RGI 6.0 (RGI Consortium, 2017) binned offsets, since not all of the input data products (e.g. glacier-wide mass-balance rates and debris cover products) are readily available as RGI 7.0 inputs to our glacier evolution model (see Section 2.2). Direct mapping of binned offsets from RGI 7.0 to RGI 6.0 is possible for the glaciers included in this study because no changes were made to their outlines between versions. For each glacier, binned offsets from survey pairs separated by only an integer number of years (
$ \pm $1 month) were differenced and used to reduce seasonal variations such as snow accumulation, melt, or firn densification (Fig. 1). We obtain surface elevation change estimates for 185 glaciers in Alaska and adjacent Yukon, Canada. Despite these representing fewer than 1% of the total number of glaciers in the region, these mostly large glaciers account for 56% of the glacierized area (Pfeffer and others, Reference Pfeffer2014). We note that Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026) process glacier-wide mass-balance rates for 205 glaciers—20 more glaciers than the elevation change dataset we use here for calibration. This discrepancy is caused by the integer-year differencing threshold applied here, in contrast to the seasonal differencing threshold used in Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026).
Example of the elevation change measurements used for model calibration for Baird Glacier, Alaska. (a) Basemap showing altimetry flight lines from Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026) and inset map showing the glacier location. RGI 6.0 glacier outlines, shown in white, highlight Baird Glacier. (b) Derived elevation changes between consecutive airborne laser altimetry surveys separated by an integer number of years.
Figure 1 Long description
Map with a grayscale basemap and a long, narrow glacier-shaped outline. Multiple colored lines run along the glacier length, with a legend listing date-pair labels: 199606, 199906, 199909, 200006, 200009, 200106, 200109, 200206, 201006, 201106, 201206, 201306, 202006. A scale bar labeled 5 kilometer appears near the lower right of the map. Line graph with multiple colored curves and a legend listing date-pair labels: 199606 colon 199906, 199906 colon 200006, 199909 colon 200009, 200006 colon 200106, 200009 colon 200109, 200106 colon 200206, 200109 colon 201006, 200206 colon 201006, 201006 colon 201106, 201106 colon 201206, 201206 colon 201306, 201306 colon 202006. The x-axis is labeled Elevation (m), with tick labels at 0, 1000 and 2000. The left y-axis is labeled Elevation change (m), with tick labels at 0, minus 20, minus 40 and minus 60. A secondary right y-axis is labeled Area (km superscript 2), with tick labels at 0, 5 and 10. A gray filled area plot is drawn against the right y-axis, spanning the elevation range and rising to a broad peak near the middle of the x-axis before decreasing toward higher elevations. The colored curves generally lie near 0 at low elevations and diverge downward at higher elevations, with the lowest curve reaching near minus 60 by around 2000 on the x-axis.
2.2. Glacier evolution model
We model glacier changes using the Python Glacier Evolution Model (PyGEM; Rounce and others, Reference Rounce2023). The model estimates the mass balance of every glacier using ∼10 m elevation bins and a monthly time step forced with air temperature and precipitation data from ERA5 (Hersbach and others, Reference Hersbach2020) for historical periods and an ensemble of general circulation models (GCMs) and shared socioeconomic pathways (SSPs) for future periods (Eyring and others, Reference Eyring2016). Glacier melt is computed using a degree-day model, accumulation from a temperature threshold, and refreezing from the mean annual temperature. Sub-debris melt enhancement factors account for the enhanced or reduced melting associated with debris thickness (Rounce and others, Reference Rounce2021). A frontal ablation parameterization is used for marine-terminating glaciers (Rounce and others, Reference Rounce2023). PyGEM is coupled to the Open Global Glacier Model (OGGM), which is used to estimate the glacier dynamics annually using a one-dimensional flowline model (OGGM; Maussion and others, Reference Maussion2019) based on the shallow-ice approximation (Hutter, Reference Hutter1983). OGGM is also used to estimate each glacier’s initial ice thickness using a previously calibrated ice viscosity parameter for the region from Rounce and others (Reference Rounce2023) that was calibrated to be consistent with regional multi-model ice thickness estimates (Farinotti and others, Reference Farinotti2019). Details regarding the model calibration framework and simulations are described below.
2.2.1. Surface elevation change
By definition, the elevation change measurements comprise changes due to the climatic mass balance and flux divergence. While PyGEM models the monthly climatic mass balance, the flux divergence is estimated annually based on the shallow-ice approximation (Maussion and others, Reference Maussion2019). We thus assume that the flux divergence is constant throughout the year and compute the total mass balance in each elevation bin as the climatic mass balance minus the flux divergence. Changes in this binned total mass balance are then computed for a specified time span and converted to elevation change using a model density factor (see Section 2.2.3).
2.2.2. Dynamical spin-up
IceBridge altimetry observations demonstrate widespread thinning across the ablation zone of glaciers in Alaska. When compared with these observations, previous model estimates (Rounce and others, Reference Rounce2023) revealed too little thinning—and in some cases thickening—at the terminus, resulting from model initialization under an assumption of equilibrium (Fig. S1). Since IceBridge data capture the transient state of glacier elevation change, we adopted a dynamical spin-up routine from OGGM (Aguayo and others, Reference Aguayo2024) to represent this transient state. This spin-up routine estimates a glacier’s past state by modeling its dynamical evolution from a specified date in the past and matching the glacier area from RGI 6.0 at its reference year, which ranges from 2000 to 2011 for our glaciers.
Integrating dynamical spin-up that accounts for differences in model parameters directly into the Bayesian inference framework is currently computationally prohibitive. We thus estimate the transient state of the glacier at the time of the first observation (2000 or earlier if IceBridge data are available) by performing an initial dynamical spin-up that is based on a single set of model parameters from Rounce and others (Reference Rounce2023), which used a simplified calibration scheme to only match the glacier-wide geodetic mass balance from Hugonnet and others (Reference Hugonnet2021) using Gaussian process mass-balance emulators. This transient dynamical state was then used to initialize the dynamical runs within the Bayesian inference framework (see Section 2.2.3). Marine-terminating glaciers also include a frontal ablation parameterization using the frontal ablation model parameter from Rounce and others (Reference Rounce2023). The duration of the spin-up period required such that modeled elevation change aligns with observations was found to vary for glaciers in this study, presumably due to differences in glacier response times (Zekollari and others, Reference Zekollari, Huss and Farinotti2020). Thus, an initial optimization routine that identifies the best spin-up duration ranging from 0 to 60 years was implemented that minimizes the mean absolute error between modeled and observed surface elevation change data. The maximum spin-up period was limited to 60 years due to the use of 1940–2021 ERA5 climate data.
2.2.3. Bayesian inference calibration framework
With PyGEMv1.1.0 presented here, we expand PyGEM’s calibration capabilities by developing a custom and flexible Metropolis–Hastings Markov Chain Monte Carlo (MCMC) sampling framework to leverage glacier-wide geodetic mass-balance data and binned elevation change data. In this framework, at each step of the MCMC chain, PyGEM’s three mass-balance parameters—temperature bias (
${T_{{\text{bias}}}})$, precipitation factor (
${k_{\text{p}}})$, and degree-day factor of snow (
${f_{{\text{snow}}}}$)—as well as the accumulation and ablation area densities (
${\rho _{{\text{acc}}}}$ and 
${\rho _{{\text{abl}}}}$) are randomly perturbed by a factor of their prior standard deviations.
As changes in the observed surface elevation may result from a combination of accumulation and ablation of ice, firn, and snow, modeled density values were used to convert binned modeled mass changes (Section 2.2.1) to elevation changes. We simplify these changes assuming separate densities for the accumulation area and the ablation area based on estimates of the equilibrium line altitude.
The log posterior probability is then calculated according to Bayes’ theorem as the sum of the log prior probability, the normal log likelihood, and the log potential function:
\begin{equation}{\text{log}}\,P{\text{(}}\theta \left| {\hat d{\text{)}}\, = \,{\text{log}}\,P{\text{(}}\theta {\text{)}}\, + \,{\text{log}}\,P{\text{(}}\hat d} \right|\theta {\text{)}}\, + \,{\text{log}}\,\phi \,{\text{(}}\theta {\text{)}}\end{equation} For each glacier, the initial flowline state was uniformly taken from the dynamical spin-up at the earliest year of available calibration data. As described by Rounce and others (Reference Rounce, Khurana, Short, Hock, Shean and Brinkerhoff2020, Reference Rounce2023), the joint prior probability distribution for both the precipitation factor and the temperature bias is approximated for each RGI second-order subregion, while the log potential function (
$\,log\,\phi \,{\text{(}}\theta {\text{)}}$) ensures that nonphysical model parameters, such as an accumulation area density factor greater than that of the ablation area, are rejected. The log likelihood assumes that the differences between the model’s predictions and the observed data follow a normal distribution. We compute the mean normal log likelihood across observation types, assigning equal weight to the binned elevation-change data and the scalar glacier-wide mass balance. This prevents the higher dimensionality of the elevation-change observations from dominating the joint likelihood and causing the sampler to preferentially fit elevation changes at the expense of the glacier-wide mass balance.
The prior distribution for the accumulation and ablation area densities is assumed to be normally distributed with mean and standard deviations of 900 
$ \pm $ 17 and 600 
$ \pm $ 60 kg m−3 (Huss, Reference Huss2013), respectively. Accumulation and ablation areas were separated based on the minimum equilibrium line altitude over 2000–2019 using OGGM’s default mass-balance model (Maussion and others, Reference Maussion2019), thereby avoiding any circularity issues with the equilibrium line altitude changing as the model parameters vary between steps within the Markov chains.
Three MCMC calibration chains with differing initial parameter states were produced for each glacier for which surface elevation change data were available to explore sampling performance. One chain was initialized from the mean value of each parameter’s prescribed prior distribution, while the other two chains were initialized from the tails of the 95% confidence intervals. We note that performing calibration of surface elevation change requires running the full dynamical model for each step in the calibration chain and thus results in a >150× slowdown compared to calibrating against the glacier-wide mass-balance rate alone.
2.2.4. Hindcast and future model simulations
Following dynamical spin-up and model calibration against both glacier-wide mass balance from Hugonnet and others (Reference Hugonnet2021) and binned elevation change observations from Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026), hindcast simulations were performed from 1994 to 2021 using ERA5 reanalysis data (Hersbach and others, Reference Hersbach2020) to allow for comparison against IceBridge elevation change data. An ensemble of 12 GCMs forced by four SSPs (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5; Table S1) was then used to simulate the mass and area change of each glacier from 2000 to 2100. PyGEM’s default 20 year (2000–19) delta change bias adjustment option (Weathers and others, Reference Weathers, Rounce, Fasullo and Maussion2025) was used to account for differences between the GCMs and ERA5 over the calibration period. For each of the 185 glaciers in this study, we ran 50 simulations for each GCM/SSP ensemble member by randomly sampling the posterior distributions from model calibration, thereby enabling us to quantify uncertainty associated with the model forcing and model parameters.
3. Results
3.1. Optimized dynamical spin-up
The duration of the dynamical spin-up strongly impacts the flux divergence, particularly at the terminus of the glacier and in earlier model years. An insufficient spin-up period often leads to excessively high terminus emergence velocities and consequent thickening, whereas an overly long spin-up period often results in low emergence velocities and overthinning of the terminus (Fig. S1). The duration of dynamical spin-up required to provide sensible terminus flux divergence estimates and allow for calibration against observed elevation change data varies in a heterogeneous manner across the region. No correlation is observed between glacier location, area, slope, aspect, elevation range, or type (land- vs marine-terminating) and the optimal spin-up period (Fig. S2). While the median spin-up period is 20 years (interquartile range = 39 years), the distribution is bimodal, with modes at 10 and 60 years.
3.2. Calibration performance
The performance of the MCMC sampling scheme was evaluated using the Gelman–Rubin statistic and effective sample size. An analysis of varying chain lengths revealed that a chain length >20 000 steps at a step size of 0.5 relative to the prior standard deviation led to robust sampling and chain convergence (Fig. 2).
MCMC sampling performance analysis. Top row shows the effective sample size (
${n_{{\text{eff}}}}$) for each parameter, while the bottom row shows the Gelman–Rubin convergence statistic (
$\hat R$). Dashed horizontal lines show desirable metric thresholds (
${n_{{\text{e}}\,{\text{ff}}}}\, \geq \,\,{\text{100}},\hat R\,\, \leq {\text{1}}{\text{.1}}$). Columns show the temperature bias (
${T_{{\text{bias}}}}$), precipitation factor (
${K_{p}}$), degree-day factor of snow (
${f_{{\text{snow}}}}$), ablation area density (
${\rho _{{\text{abl}}}}$), and accumulation area density (
${\rho _{{\text{acc}}}}$), respectively.
Figure 2 Long description
Diagnostic line plots evaluating MCMC sampling performance for five parameters: temperature bias, precipitation factor, degree-day factor of snow, ablation area density and accumulation area density. Each parameter is represented in two rows. The top row shows the effective sample size (n subscript eff) with the vertical axis ranging from 0 to 2000 and the horizontal axis representing steps in thousands, from 0 to 30. Dashed horizontal lines indicate desirable metric thresholds. The solid line represents the median, with shaded intervals showing uncertainty. The n subscript eff generally increases with steps, crossing the threshold around 20,000 steps for most parameters. The bottom row displays the Gelman-Rubin convergence statistic (R hat), with the vertical axis ranging from 1.0 to 1.2. The horizontal axis is the same as the top row. The R hat values decrease toward 1.0, indicating convergence, with most parameters reaching below the threshold by 20,000 steps. The plots collectively show sampling efficiency and convergence over time, with the temperature bias parameter converging fastest. All panels share similar styling with solid lines and shaded intervals, providing a consistent visual representation across parameters.
The posterior distributions of the temperature bias, precipitation factor, and the degree-day factor of snow were compared to those produced for the model being calibrated against glacier-wide mass balance alone (PyGEMv0.2.0; Fig. 3). The per-glacier median temperature bias decreased by a mean of 0.1°C, while the precipitation factor increased by a mean of 0.2, and the degree-day factor of snow increased by a mean of 0.3 mm w.e. d−1 °C−1. Additionally, the per-glacier confidence interval of each model parameter decreased, most notably for the degree-day factor of snow where the mean decrease was 0.8 mm w.e. d−1 °C−1.
Change in the model posterior distributions for each parameter between PyGEMv1.1.0 (this study) and PyGEMv0.2.0 (previous version; Rounce and others, Reference Rounce2023). The left, middle, and right columns refer to the temperature bias (°C), precipitation factor (–), and degree-day factor of snow (mm w.e. d−1 °C−1), respectively. The top row shows the prior (black) and posterior distributions for PyGEMv0.2.0 (blue) and PyGEMv1.1.0 (orange). The middle row shows a comparison between the per-glacier median values (med), while the bottom row shows a comparison between the per-glacier 95% confidence intervals (CI). For the middle and bottom rows, the upper left corner of each panel shows the mean difference.
Figure 3 Long description
The top row contains three density plots comparing prior and posterior distributions for different parameters. The first plot shows temperature bias with density on the vertical axis and T subscript bias on the horizontal axis, ranging from negative 5 to positive 5. The second plot displays the precipitation factor, K subscript p, ranging from 0 to 10. The third plot illustrates the degree-day factor of snow, f subscript snow, ranging from 0 to 10. Each plot includes lines for prior, Post v0.2.0 and Post v1.1.0 distributions. The middle row presents scatter plots comparing per-glacier median values for each parameter between PyGEMv0.2.0 and PyGEMv1.1.0. The first scatter plot shows med(T subscript bias, v1.1.0) versus med(T subscript bias, v0.2.0) with a mean difference of negative 0.1. The second scatter plot displays med(K subscript p, v1.1.0) versus med(K subscript p, v0.2.0) with a mean difference of 0.2. The third scatter plot shows med(f subscript snow, v1.1.0) versus med(f subscript snow, v0.2.0) with a mean difference of 0.3. The bottom row contains scatter plots comparing per-glacier 95 percent confidence intervals for each parameter. The first plot shows CI(T subscript bias, v1.1.0) versus CI(T subscript bias, v0.2.0) with a mean difference of negative 0.2. The second plot displays CI(K subscript p, v1.1.0) versus CI(K subscript p, v0.2.0) with a mean difference of 0.3. The third plot illustrates CI(f subscript snow, v1.1.0) versus CI(f subscript snow, v0.2.0) with a mean difference of negative 0.8.
3.3. Hindcast elevation change
Following dynamical spin-up and model calibration, 1994–2021 hindcast simulations were run to evaluate modeled elevation change rates against observations. We calculate the mean absolute error between observed and modeled elevation change for each glacier and compare it to the same mean absolute error for the model being calibrated against only glacier-wide mass balance (PyGEMv0.2.0; as in Rounce and others, Reference Rounce2023) to quantify the changes in model performance. Overall, the mean absolute error of the binned surface elevation change was reduced from a mean (by glacier count) of 1.24 
$ \pm $ 0.63 to 1.06 
$ \pm $ 0.63 m a−1, with improvement varying considerably by glacier (Fig. 4). The mean absolute error in the binned surface elevation change decreased for 139 glaciers, while it increased for 34 of the glaciers (Figs 4a and S3). The simulated glacier-wide mass balance from 2000 to 2019 stayed consistent at approximately −0.96 
$ \pm $ 0.68 m w.e. a−1.
Comparison between observed and modeled elevation change rates for the new model calibration framework (PyGEMv1.1.0; orange) and old version (PyGEMv0.2.0; blue). (a) Per-glacier change in the average mean absolute error (MAE), with cooler colors indicating a reduction, as well as histograms showing the overall distribution of the (top) MAE for each model version and (bottom) the change in MAE between versions. (b–d) Examples of the performance compared to observations (black) and their associated uncertainties (gray shading) for Gulkana, Kahiltna, and Kennicott Glaciers, respectively. The locations of these three glaciers are indicated in panel (a), along with the change in their MAE.
Figure 4 Long description
The composite image includes a map and graphs comparing elevation change rates for glaciers in Alaska using PyGEMv1.1.0 (orange) and PyGEMv0.2.0 (blue) models. The map shows per-glacier change in mean absolute error (MAE) with cooler colors indicating reduction. Locations of Gulkana, Kahiltna and Kennicott glaciers are marked with their elevation change rates: Gulkana (-0.24 m per year), Kahiltna (-0.05 m per year), Kennicott (-0.35 m per year). Insets display histograms of MAE distribution for each model version and change in MAE between versions. Graphs b, c and d show elevation change over time for Gulkana, Kahiltna and Kennicott glaciers, respectively, with observations in black and uncertainties in grey shading. Each graph is divided into time periods: 2000-2005, 2013-2016, 2016-2017 for Gulkana; 2010-2013, 2013-2016 for Kahiltna; 2000-2013, 2007-2018 for Kennicott. Elevation is plotted on the x-axis in meters and elevation change on the y-axis in meters per year. The graphs illustrate model performance compared to observations, highlighting differences in elevation change rates and associated uncertainties.
While the overall quantitative misfit broadly shows improvement, we highlight representative examples comparing these changes in metrics with the available surface elevation change data for Gulkana Glacier, Kahiltna Glacier, and Kennicott Glacier, chosen in part for their long-term monitoring records from the United States Geological Survey and the National Park Service. While Gulkana Glacier shows only a moderate improvement in model fit of 0.24 m a−1, we observe clear improvement at all elevations, particularly for the earliest time spans (2000–16; Fig. 4b). Though Kahiltna Glacier demonstrates a subtle improvement of 0.05 m a−1, we observe some improvement in model fit at both ends of the glacier (Fig. 4c). For Kennicott Glacier, PyGEMv0.2.0 hindcast results agree poorly with observations below ∼2000 m elevation (Fig. 4d). In fact, the modeled elevation change is negatively correlated with elevation and implies terminus thickening, whereas observations show elevation change rates decreasing with elevation. In contrast, PyGEMv1.1.0 results align with the observations and provide an overall improvement in the misfit of 0.35 m a−1. Glaciers that result in an increased mean absolute error are largely those for which PyGEMv0.2.0 hindcast elevation changes align with observations; hence, any increases are fairly minimal (Fig. S4).
3.4. Regional projections of glacier change
For each SSP, we present the ensemble (N = 12) median (and associated 95% confidence interval) annual glacier mass and area for each of the 185 glaciers in this study. Glaciers of the Chugach Mountains (Knik, Matanuska, Harvard) and Wrangell–St. Elias Mountains (Copper, Nabesna, Bering, Yahtse, Guyot, Kaskawulsh, Klutlan, Donjek, Logan, Walsh, Chitna, Seward and Hubbard) generally lose the least mass across SSP scenarios, each retaining at least 50%, and up to 94% of their mass at 2100, relative to 2015, for SSP2-4.5 (Figs 5 and S5–S7). Mass losses increase for greater emissions scenarios, with Yale, Bering, and Hubbard being the only glaciers with at least 50% of their mass remaining by the end of the century for the highest emissions scenario (SSP5-8.5).
Projections of glacier mass loss through 2100, relative to 2015, showing (a) each glacier for SSP2-4.5, and temporal variations in (b) mass and (c) area for all four emissions scenarios for all 185 glaciers as well as the glaciers aggregated by National Park units (dashed boundaries). For the time series, the solid lines show the ensemble median, while box-and-whisker plots show the ensemble interquartile range and 95% confidence interval at year 2100. The total glacierized area and percent thereof included in this study are indicated for each region at the top of panel (b).
Figure 5 Long description
The image includes a map and multiple graphs projecting glacier mass and area changes through 2100. The map highlights regions like Denali, Lake Clark, Kenai Fjords, Wrangell-St. Elias and Glacier Bay, with color coding indicating mass at 2100 relative to 2015 under SSP2-4.5. The legend shows a range from 0.0 to 1.0, representing the fraction of mass retained. Graphs are divided into two sections: mass (B) and area (C) for each region. The x-axis shows years from 2000 to 2100 and the y-axis shows mass and area relative to 2015. Each graph includes lines for SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5 scenarios, with ensemble median lines and box-and-whisker plots indicating interquartile ranges and 95 percent confidence intervals at 2100. For mass, regions like RGI Region 1 and Glacier Bay show significant losses under higher emission scenarios, with SSP5-8.5 showing the steepest decline. Area graphs show similar trends, with area loss accelerating over time. The graphs highlight that higher emission scenarios lead to greater mass and area loss, with SSP1-2.6 retaining more mass and area compared to SSP5-8.5. Visual differentiation includes distinct line styles for each scenario and box-and-whisker plots for variability at 2100. The description provides a comprehensive understanding of the projected changes in glacier mass and area across different scenarios and regions.
Projections were aggregated to assess glacier change across Alaska and within five of the state’s largely glacierized National Park units, with the glacier cover in each unit from Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026). Altogether, the 185 glaciers included in this study represent 56% of the total regional glacier area and are projected to lose 51 
$ \pm $ 22% (6648 
$ \pm $ 2858 Gt) of their mass and 29 
$ \pm $ 19% (13 700 
$ \pm $ 8881 km2) of their area by the end of the century for SSP2-4.5, relative to 2015 (Fig. 5). For the highest emissions scenario (SSP5-8.5), the model estimates the glaciers will lose 67 
$ \pm $ 27% (8760 
$ \pm $ 3593 Gt) of their mass and 48 
$ \pm $ 27% (22 196 
$ \pm $ 12 454 km2) of their area. Projected glacier changes vary across Alaska’s national parks. Wrangell–St. Elias National Park and Preserve, with the most glacier cover and the largest glaciers, is projected to lose between 28 
$ \pm $ 16% of its glacier mass for SSP1-2.6 and up to 51 
$ \pm $ 34% for SSP5-8.5. Conversely, glaciers in Lake Clark and Kenai Fjords are generally smaller and are projected to lose 83 
$ \pm $ 20% and 59 
$ \pm $ 21% of their total mass, relative to 2015, for SSP1-2.6, respectively, as opposed to 99 
$ \pm $ 6% and 86 
$ \pm $ 19% for SSP5-8.5.
For several of the most visited and well-studied glaciers, we redistribute the flowline area and ice volume to a two-dimensional (2-D) grid to better visualize the changes (Fig. 6). Tanaina Glacier in Lake Clark National Park and Preserve is projected to lose between 79 
$ \pm $ 29% and 98 
$ \pm $ 7% of its mass by 2100, relative to 2015, depending on the emissions scenario, with ice remaining only at elevations greater than 1300 m (Fig. S8). We project Gulkana Glacier, one of the United States Geological Survey’s benchmark glaciers with mass-balance measurements dating back to 1966 (Florentine and McKeon, Reference Florentine and McKeon2022), will lose between 92 
$ \pm $ 19% and 100 
$ \pm $ 5% of its mass by 2100, relative to 2015, with very little high-elevation ice remaining only for SSP1-2.6 and SSP2-4.5, and being entirely deglaciated for SSP3-7.0 and SSP5-8.5 (Fig. S9). We estimate that Exit Glacier in Kenai Fjords National Park will lose between 97 
$ \pm $ 11% and 100 
$ \pm $ 4% of its mass, with the glacier’s lower tongue retreating to the Harding Icefield before mid-century, regardless of the emissions scenario (Fig. S10). Kahiltna Glacier in Denali National Park and Preserve and the combined Kennicott–Root glaciers in Wrangell–St. Elias National Park and Preserve are projected to lose ∼45% of their mass under SSP1-2.6 and 70% for SSP5-8.5, accompanied by substantial retreat (Figs S11 and S12). Note that our projections estimate ∼10–20% more mass loss for the Kennicott–Root glacier complex than estimated by Wells and others (Reference Wells2025). Margerie Glacier in Glacier Bay National Park and Preserve is projected to lose between 28 
$ \pm $ 25% and 63 
$ \pm $ 41% of its mass by the end of the century, retreating between ∼5 and 20 km (Fig. S13).
Projections of glacier thickness (H) and mass for Tanaina, Kahiltna, Exit, Gulkana, Kennicott–Root, and Margerie Glaciers, respectively, from top to bottom. The three-dimensional images show the gridded glacier thickness in meters in 2015 and 2100 for SSP2-4.5 and SSP5-8.5. The time series show the total annual mass, relative to the year 2015, for SSP2-4.5 (green) and SSP5-8.5 (orange), where the solid line is the ensemble median, and the shaded area represents the 95% confidence intervals. RGI 6.0 glacier area is noted to the right of each color bar.
Figure 6 Long description
Projections of glacier thickness and mass for Tanaina, Kahiltna, Exit, Gulkana, Kennicott-Root and Margerie Glaciers. Each row shows three 3D images of glacier thickness in meters for 2015, 2100 under SSP2-4.5 and SSP5-8.5 scenarios. The thickness is indicated by a color gradient, with a scale bar on the left showing the range in meters. Adjacent to these images, time series graphs display the total annual mass relative to 2015 for SSP2-4.5 (green) and SSP5-8.5 (orange). The solid line represents the ensemble median and the shaded area indicates the 95 percent confidence intervals. The graphs show a decline in mass over time, with SSP5-8.5 generally showing a steeper decline compared to SSP2-4.5. The RGI 6.0 glacier area is noted to the right of each colorbar, providing additional context for the glacier′s size and changes over time.
4. Discussion
4.1. Expanded model calibration: advances and limitations
The Bayesian calibration framework presented herein demonstrates the ability to leverage spatially distributed elevation change measurements for PyGEM calibration. While we showcase this capability using altimetry-derived elevation change data for 185 glaciers in Alaska, the framework is well suited to leverage other spatially distributed elevation change datasets from DEM differencing (e.g. Dehecq and others, Reference Dehecq2020; Berthier and others, Reference Berthier2024) and altimetry (e.g. Sochor and others, Reference Sochor, Seehaus, Braun, Sochor, Seehaus and Braun2021; Wang and others, Reference Wang, Yi and Sun2021; Jakob and Gourmelen, Reference Jakob and Gourmelen2023; Hassan and others, Reference Hassan2025). As the number of available elevation change datasets continues to grow, this framework is uniquely suited to ingest these data along with their associated uncertainties, enabling differences among datasets to be handled within the calibration process rather than requiring modelers to select only a single product or determine how to aggregate different products. For this to be effective though, the uncertainty associated with the data, as well as its assumed distribution (e.g. normal), must be accurately characterized.
Calibrating against spatially distributed elevation change data shifts the posterior distributions of all model parameters compared to calibrating against glacier-wide mass-balance rates alone. Moreover, it improves the confidence intervals for the temperature bias, precipitation factor, and degree-day factor of snow on average by 5%, 16%, and 18%, respectively. These improvements translate to reduced uncertainty in glacier change projections from the model parameters, although the overall uncertainty in mass loss projections remains the same as it is dominated by the uncertainty associated with the climatic forcing. The PyGEMv1.1.0 projections do estimate slightly less mass loss at 2100, relative to 2015, compared to PyGEMv0.2.0, ranging from 3.1% less for SSP1-2.6 to 1.3% less for SSP5-8.5 (Fig. S14), although these differences are attributed to both the differences in model calibration as well as model spin-up. Nonetheless, despite the improved model calibration and spin-up framework, persistent model-data elevation change mismatches remain for some glaciers, particularly large and marine-terminating glaciers (Fig. S3). For some marine-terminating glaciers, the model predicts terminus thickening where altimetry observations indicate thinning (e.g. Yale Glacier), while for others, thinning is overestimated (e.g. Hubbard Glacier; Fig. S4). Improving the representation of frontal ablation (e.g. Malles and others, Reference Malles, Maussion, Ultee, Kochtitzky, Copland and Marzeion2023) may improve model performance for marine-terminating glaciers. Large glaciers with many tributaries and complex geometries also pose challenges for the elevation-band flowline model, in which a glacier’s 2-D shape is reduced to a bulk-average flowline (Huss and Farinotti, Reference Huss and Farinotti2012; Maussion and others, Reference Maussion2019), limiting the model’s ability to reproduce spatial patterns of elevation change. Coupling PyGEM to higher-order dynamical ice models (e.g. Brinkerhoff, Reference Brinkerhoff2023; Jouvet and Cordonnier, Reference Jouvet and Cordonnier2023; Brinkerhoff and others, Reference Brinkerhoff2025) may alleviate this by improving our ability to capture the dynamical component of thickness or elevation change.
In addition to structural model limitations, the spatial elevation change observations used for calibration are often characterized by substantial noise, even after filtering and smoothing. Many glaciers in Alaska are also known to surge episodically (Post, 1969), which is expressed in the elevation change record presented by Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026; e.g. Bering Glacier). While the MCMC calibration framework seeks to reduce misfit in the mean binned surface elevation change, performance can be degraded where observations are dominated by noise or characterized by surging, which is not represented in the model.
Despite these challenges, this calibration framework is highly flexible and is not restricted to elevation change data, allowing integration of other surrogate measurements of glacier change, such as snowlines (e.g. Bevington and Menounos, Reference Bevington and Menounos2025) and melt extents (Wells and others, Reference Wells, Rounce and Fahnestock2026), as these data become increasingly available. While many of these emerging products remain uncertain, lack extensive in situ validation, or are not yet broadly available, their integration within PyGEM’s Bayesian framework provides a pathway for future multi-constraint calibration and/or synthesis efforts. Future work should thus seek to leverage these diverse datasets and expand this framework globally to refine hindcasts and projections.
4.2. Hazards and park management challenges associated with glacier change
Our projections reveal continued rapid deglaciation across Alaska, with mass loss ranging from 43% to 67% in total for the glaciers represented in this study. Because approximately half of Alaska’s glacierized land area falls within national parks, these future glacier changes will impose a number of park management-related challenges, including hazard mitigation, wildlife and landscape management, and the maintenance of park infrastructure and glacier access trails in dynamically evolving glacial and postglacial landscapes.
As marine-terminating glaciers, such as those in Kenai Fjords, Glacier Bay, and Wrangell–St. Elias, continue to retreat, their termini will become grounded and cease calving, impacting local marine habitats and the tourism industry. This is demonstrated by projections of Margerie Glacier, which foretell the glacier’s terminus retreat by at least ∼5 km by the end of the century (Fig. 6). These projected changes build on the extensive fjord expansion already documented across the region in response to historical retreat (e.g. McNabb and Hock, Reference McNabb and Hock2014; Carlson and others, Reference Carlson2017).
Projections indicate continued thinning and retreat of land- and lake-terminating glaciers across the region throughout the 21st century, which will lead to the expansion and formation of proglacial lakes, extending the regional pattern of lake growth observed in recent decades (e.g. Shugar and others, Reference Shugar2020; Field and others, Reference Field, Armstrong and Huss2021; Black and Kurtz, Reference Black and Kurtz2023). For example, at Kennicott Glacier, where two proglacial lakes have grown to a cumulative area of 1.6 km2 in response to terminus stagnation and retreat (Petersen and others, Reference Petersen2025; Wells and others, Reference Wells2025), terminus retreat is projected to exceed 10 km by 2100, which—based on best estimates of bed topography (Wells and others, Reference Wells2025)—will drive further proglacial lake growth.
Continued glacier thinning across the region may increase the potential for ice-dam failure at glacier-marginal lakes, which can trigger catastrophic outburst floods and pose major downstream hazards. For example, future glacier retreat may further prime Bear Glacier in Kenai Fjords National Park to produce floods that endanger recreationalists at Bear Lake (Wilcox and others, Reference Wilcox, Wade and Evans2014), while those from Kennicott Glacier’s marginal Hidden Creek Lake threaten residents, visitors, and transportation infrastructure in Wrangell–St. Elias National Park and Preserve (Rickman and Rosenkrans, Reference Rickman and Rosenkrans1997). Historical events across Alaska demonstrate the significance of such hazards (e.g. Post and Mayo, Reference Post and Mayo1971; Kienholz and others, Reference Kienholz2020); however, these threats will eventually diminish once glaciers retreat so much that the ice dams are no longer present. Conversely, glacier thinning and retreat may also make the region more prone to landslides that have the ability to generate tsunamis and could affect tourism (Shugar and others, Reference Shugar2026).
As Kahiltna Glacier continues to thin, crevasse fields may become more routinely exposed at high-elevation glacier airstrips, posing increasing challenges for ferrying mountaineers to Denali National Park and Preserve’s Mount McKinley. Our projections indicate Exit Glacier will become unviewable from Kenai Fjords National Park’s Glacier Overlook Loop Trail by mid-century and approach near-complete deglaciation by the end of the century. Providing park visitors direct access to Root Glacier in Wrangell–St. Elias National Park and Preserve will also be more challenging and may become unfeasible within the next several decades as the glacier continues its dramatic thinning and retreat. Expected changes such as these may adversely impact visitor satisfaction to Alaska’s national parks.
While the glaciers included in this study do not encompass all glacierized area within Alaska’s national parks or the region as a whole, they account for a sizeable share, including over half of Alaska’s total glacierized area and roughly two-thirds to three-quarters within Denali, Kenai Fjords, Wrangell–St. Elias, and Glacier Bay National Park and Preserves. We thus consider these glaciers representative of the broader regional trends, which together foretell widespread glacier loss across Alaska. Even under the lowest emissions scenario (SSP1-2.6), 43 
$ \pm $ 17% of the regional glacier mass and 24 
$ \pm $ 12% of the regional glacier area represented here will be lost by 2100, relative to 2015. With the exception of Wrangell–St. Elias, which we project will lose up to 24 
$ \pm $ 24% of its total glacier area, Alaska’s national parks may lose between ∼30 and 90% of their glacier area, depending on the emissions scenario, resulting in widespread landscape changes.
5. Conclusion
We present a novel calibration framework for the open-source PyGEM, in which the model’s Bayesian inference calibration is expanded to leverage spatially distributed elevation change data. A dynamical spin-up routine was implemented to capture the transient state of glacier change evident in the elevation change data. This routine functions as a short-term glacier reanalysis, reconstructing the glacier state over the recent past by modeling its evolution from a specified time such that its area matches the observed area in a reference year. We showcase our expanded calibration functionality using the comprehensive record of airborne altimetry-derived elevation change observations presented by Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026) for Alaska.
Calibrating against spatially distributed elevation change data results in a subtle shift of the posterior distributions of PyGEM’s model parameters and tightens the confidence interval of each parameter by up to 18%, while improving the alignment of hindcast simulations to elevation change observations, with the mean absolute error decreasing from 1.23 
$ \pm $ 0.63 to 1.06 
$ \pm $ 0.63 m a−1. The model results presented here indicate widespread glacier retreat across Alaska, likely resulting in significant landscape changes, and potentially increasing future hazards. For the 185 glaciers in this study, projections estimate glaciers will lose from 43 
$ \pm $ 17% to 67 
$ \pm $ 27% of their aggregated mass by 2100, relative to 2015. We discuss some of the significant impacts these changes may have on Alaskan national parks. Aggregated mass losses vary by park unit, with the potential for nearly complete deglaciation of Lake Clark National Park and Preserve, while between 28 
$ \pm $ 16% and 51 
$ \pm $ 34% of glacier mass is projected to remain for Wrangell–St. Elias National Park and Preserve. Future efforts should leverage additional spatially distributed elevation change measurements for the remaining glaciers in Alaska and elsewhere to refine model calibrations and improve projections.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10186.
Data availability statement
PyGEM version 1.1.0 is publicly available online at 10.5281/zenodo.20278458. The elevation change offsets used for model calibration are provided by Loso and others (Reference Loso, Fahnestock, Larsen and Tober2026). Model calibration, projections, and code supplementing this manuscript are publicly available at 10.1184/R1/32319708. 3D glacier visualizations were produced using the open-source Glacier:3D-Viz tool (10.5281/zenodo.18338218).
Acknowledgements
B.S.T. and D.R.R. were primarily supported by the National Park Service Inventory and Monitoring Program’s Southeast Alaska Network through Department of the Interior (DOI) award P22AC02208, with additional support from National Aeronautics and Space Administration (NASA) awards 80NSSC24K1530 and 80NSSC24K1575. Thanks to M. Fahnestock, C. Larsen, and others at the University of Alaska Fairbanks who collected and analyzed the glacier surface elevation change measurements utilized in this study. The authors are grateful for the constructive feedback provided by the scientific editor (H. Jiskoot) and two anonymous reviewers.
Open access
