Surface mass balance modelling of the Juneau Icefield highlights the potential for rapid ice loss by the mid-21st century

Abstract

Plateau icefields are large stores of fresh water, preconditioned to enhanced mass loss due to their gently sloping accumulation areas. Accurate modelling of their mass balance is therefore crucial for sea-level rise projections. Here, we use the COupled Snowpack and Ice surface energy and mass-balance model in PYthon (COSIPY) to simulate historical and future mass balance of the Juneau Icefield, Alaska – a high elevation (>1200 m) plateau icefield. We force the model with dynamically downscaled climate simulations, for both past and future (RCP 8.5) conditions. The icefield's mass balance decreased from a mean of −0.22 ± 0.38 m w.e. a⁻¹ (1981–2019) to −1.52 ± 0.27 m w.e. a⁻¹ (2031–2060), with many glaciers shifting from positive to negative mass balances at the start of the 21st century. This mass loss is attributed to projected rising air temperatures and reduced snowfall, causing the equilibrium line altitude to rise and triggering albedo and melt-elevation feedbacks. These processes exacerbate melt, potentially leading to increased glacier disconnections at icefalls.

1. Introduction

Glaciers in Alaska are the largest source of released fresh water, having produced one-third of the total sea-level rise contribution from glaciers between 1961 and 2016 (Zemp and others, 2019). Although glaciers outside of the Antarctic and Greenland ice sheets contain a small proportion of global ice volume (~1%), they are disproportionally contributing to sea level rise (Marzeion and others, 2012; Hugonnet and others, 2021). Global-scale projections indicate that these glaciers will continue to disproportionally contribute to sea level rise, and that glaciers in Alaska will remain the largest suppliers of meltwater (Hock and others, 2019; Edwards and others, 2021).

In Alaska, the majority of glacial ice is held within its several icefields and ice caps situated along the mountainous coastal regions (Davies and others, 2022, 2024). Though global-scale projections show a linear loss of ice mass with temperature increases (Rounce and others, 2023), the distinctive hypsometry, or glacier area-elevation distribution, of plateau icefields like Juneau Icefield (Fig. 1), makes them more susceptible to rapid ice loss. This occurs through three key processes: (i) once the equilibrium line altitude (ELA) rises to the elevation of the plateau, small increases in ELA lead to a rapid decrease in accumulation area and widespread thinning (Bodvarsson, 1955; Bolibar and others, 2022); (ii) ELA rises cause the exposure of bare ice, which has a lower albedo than snow, promoting melting through an albedo-feedback (Johnson and Rupper, 2020); (iii) dynamic thinning across icefalls flowing over steep topography can lead to the disconnection of glacier tongues from the plateau, leaving them devoid of an accumulation area (Rippin and others, 2020; Davies and others, 2022). Thus, reliable forecasts of Alaskan glacier mass balance are critical for projections of global sea level rise.

Figure 1. Location of the Juneau Icefield and key glaciers (red outline) mentioned within the text. The low-slope interconnected high-elevation plateau area is highlighted from Davies and others (2022).

Among the largest of the Alaskan icefields is Juneau Icefield (Fig. 1). Situated within the Coast Mountain Range, and on the border between Alaska and British Columbia. Juneau Icefield consists of a high elevation (>1200 m) plateau which is drained by 40 outlet glaciers (Sprenke and others, 1999; Davies and others, 2022). The majority of the outlet glaciers terminate terrestrially, though several terminate in proglacial lakes (e.g. Llewellyn, Tulsequah, Mendenhall and Gilkey glaciers), While the largest outlet, Taku Glacier, is a tidewater glacier (Davies and others, 2022). Detailed observational assessments of the mass balance of the Juneau Icefield have been conducted. The most comprehensive of these comes from the Juneau Icefield Research Program, which has collected direct field measurements of accumulation and ablation at Lemon Creek and Taku glaciers annually since 1947 (LaChapelle, 1954; Wilson, 1959; Miller, 1975; O'Neel and others, 2019; McNeil and others, 2020). Geodetic methods have also provided estimates of the surface mass balance (SMB) of Juneau Icefield since 1948 (Larsen and others, 2007, 2015; Berthier and others, 2010, 2018; Larsen, 2010; Melkonian and others, 2014). Of these, Larsen and others (2015) found that all glaciers had a negative mass balance between ~1994 and 2013, except for Taku Glacier, which subsequent measurements have shown is now receding (McNeil and others, 2020). Additionally, Davies and others (2024) reported a doubling of mass loss from 1979 to 2010 to 2010–2020, driven by widespread thinning across the high-elevation plateau that triggered a series of ice-elevation feedbacks.

Although observations of Juneau Icefield are rich, modelling icefield-wide mass balance remains challenging. The mountainous terrain and maritime climate of the Juneau Icefield makes modelling the climate in the region difficult as it requires high-resolution climate models to capture orographic processes (Bozkurt and others, 2019). Global-scale projections of mass balance, of course, include the Juneau Icefield within their domain (e.g. Hock and others, 2019; Rounce and others, 2023) but the scale of these experiments means that they use simplified mass balance models. Examples include positive-degree day (e.g. Marzeion and others, 2012; Huss and Hock, 2015) and simplified energy balance models (e.g. Giesen and Oerlemans, 2013), that neglect the full energy balance of glacier surfaces. Additionally, these models commonly rely on coarse resolution global climate models as input. Although often statistically downscaled to higher resolutions, such approaches underrepresent orographic processes (Roth and others, 2018). For the Juneau Icefield, this means that detailed spatial patterns of mass change are unobtainable from global-scale models. The first study to specifically model the mass balance of the Juneau Icefield is that of Ziemen and others (2016) that used a dynamically downscaled climate model (CCSM dynamically downscaled to 20 km resolution) as input to a positive degree-day model. These authors found that a fit to mass balance and ELA observations was only obtainable by manually introducing a precipitation gradient across the icefield; a choice designed to compensate for the lack of orographic processes producing a rain shadow effect, but one that lacks physical realism (Roth and others, 2018). Ziemen and others (2016) also found that statistically downscaled climate data also provided a poor fit to observations, due to simplistic elevation dependent downscaling methods for precipitation, despite being a higher resolution (2 km). Similarly, statistically downscaled climate data (MERRA-2 downscaled to 200 m) was used in Young and others (2021), the second and most recent study to model the mass balance of the icefield. Despite not accounting for orographic processes, their simulations provided a good fit to the geodetic observations of Berthier and others (2018), against which their model was calibrated.

In this study, we aim to improve the SMB modelling of the Juneau Icefield. We utilise high-resolution, dynamically downscaled climate model simulations for Southeast Alaska from Lader and others (2020), and employ a more complex mass-balance model than has previously been applied to the Juneau Icefield – the COupled Snowpack and Ice surface energy and mass-balance model in PYthon (COSIPY; Sauter and others, 2020). The long glaciological mass balance time series from the Juneau Icefield Research Program (e.g. O'Neel and others, 2019; McNeil and others, 2020) provides an ideal opportunity to optimise the model, and assess the accuracy of historical simulations. We then perform projections of the near future (2031–2060, RCP8.5) mass balance of the Juneau Icefield and use these simulations to evaluate drivers of future change.

2. Methods

Our overall approach is summarised in Figure 2, with subsequent sections providing more detail on each stage. We input datasets (Section 2.2) pertaining to the topography, glaciology and climate (Lader and others, 2020) of the Juneau Icefield into the COSIPY surface energy and mass balance model (Section 2.1).

Figure 2. Overview of the workflow of this study. The schematic highlights which of the three climate model outputs from Lader and others (2020) are used in each experiment. Acronyms presented in this figure are defined in Section 2.2.

For the climate, we utilise three models for our mass balance simulations. A reanalysis model: the Climate Forecast System Reanalysis (CFSR; Saha and others, 2010) which utilises observations to generate a ‘best estimate’ of the past climate. And two ‘free-running’ global climate models: Geophysical Fluid Dynamics Laboratory Climate Model version 3 (GFDL-CM3; Donner and others, 2011) and Community Climate System Model version 4 (NCAR-CCSM4; Gent and others, 2011) which are forced by a time series of atmospheric greenhouse gas concentrations, solar variability and volcanic activity, but are unconstrained by observations. These simulations should capture the part of climate variability that is due to these external forcings, hence should produce general climatic patterns and trends, but not necessarily the exact timing of individual events. The reanalysis and climate datasets were projected onto a single grid system, and the climate datasets were additionally bias-corrected (Section 2.2).

The reanalysis climate data were used to drive a perturbed parameter ensemble of simulations. Comparison of these reanalysis-driven simulations to mass-balance observations was used as the basis of model optimisation (Section 2.3). The best-fitting parameter combination was used to drive mass-balance simulations from historical climate simulations, to evaluate whether the global climate models adequately conform to the broad climatic patterns and trends across the Juneau Icefield, before they were utilised in projections. As such, we produced three sets of simulations of the reference SMB of the Juneau Icefield; an evaluation simulation spanning 1981–2019 which utilised reanalysis data, a historical simulation utilising global climate model simulations (1981–2010), and projections of future SMB under RCP 8.5 between 2031 and 2060.

2.1. COSIPY

COSIPY, an updated version of COSIMA (COupled Snowpack and Ice surface energy and mass-balance model in MAtlab), is an open-source mass-balance model that provides a flexible framework for modelling snow and glacier mass changes (Huintjes and others, 2015; Sauter and others, 2020). Since its release it has been used for various applications; from simulating the behaviour of Schiaparelli Glacier, Chile, during the Little Ice Age (Weidemann and others, 2020) to investigating the sensitivity of Halji Glacier in the Himalaya (Arndt and others, 2021). Since COSIPY is a 1-D model, it neglects any lateral and basal processes such as snowdrift or lateral mass and energy fluxes. As such, COSIPY represents a middle ground in complexity when compared to other SMB models. Here, we include only a brief description of the model. For a full description of COSIPY, see Sauter and others (2020).

COSIPY consists of two coupled models: a surface energy balance model and a multilayer subsurface snow and ice model. The model assumes full conservation of mass and energy in the snowpack, with the total mass balance consisting of the SMB and internal mass balance. The SMB is calculated through the sum of surface accumulation (snowfall and deposition) and ablation (surface melt and sublimation) processes. The calculated surface melt is used as an input to the subsurface model which accounts for many subsurface processes including meltwater percolation, retention, refreezing and subsurface melt, all of which are resolved in a vertical layer structure.

Snowfall at each gridpoint is calculated using a logistic transfer function combining precipitation, temperature and relative humidity (Hantel and others, 2000). Fresh snow is added in the model only if snowfall exceeds 0.001 m day⁻¹. Snow albedo (a_snow) is parameterised using the approach by Oerlemans and Knap (1998), whereby

(1)$$a_{snow} = a_{\,firn\;} + ( {a_s-\;a_{\,firn\;}} ) \exp \left({\displaystyle{s \over {\tau^\ast }}} \right), \;$$

where a_firn and a_s are the albedo values for firn and fresh snow as set by the user. τ* is the albedo timescale and defines the time taken for the snow albedo to drop from fresh snow to firn, s defines the number of days since the last snowfall. The model also considers the effect of snowpack thickness on albedo. After falling, snow settles and compacts during metamorphism, causing the density of the underlying snowpack to increase, impacting thermal conductivity and liquid water content. The rate of density change of each snow layer is calculated using the method described in Boone (2002), where the densification happens through the overburden of pressure (Essery and others, 2013). At the surface, turbulent fluxes are parameterised based on a flux-gradient similarity theory. A linear roughness length relationship with time is used to calculate the turbulent fluxes using the method from Mölg and others (2012).

COSIPY was initialised with the default parameterisations, with the model parameters optimised via the procedure described in Section 2.3.

2.2. Input data to COSIPY

Two types of input data are required by COSIPY: static data and climate data. Static data do not change during a model simulation. These consist of a shapefile defining a glacier mask and a digital elevation model (DEM). For our region of interest, we use the glacier outlines of Juneau Icefield in 2019 from Davies and others (2022). The DEM used in this study is derived from the Shuttle Radar Topography Mission, which has a global coverage of 1-arc (~30 m). This was aggregated to the coarser resolution of 600 m, which formed the model grid. From the aggregated DEM, other input grids of surface slope, aspect, hill shading and an ice mask were calculated using the Geospatial Data Abstraction Library (GDAL). COSIPY does not account for the evolution of an ice surface or margin during a model simulation, hence, the SMB presented here should be considered as reference SMB, according to the definitions of Cogley and others (2011).

To force COSIPY, we utilise reanalysis and global climate model data that were dynamically downscaled to an hourly temporal and 4 km spatial resolution using the Weather and Research Forecasting model version 4 (Skamarock and others, 2019) by Lader and others (2020). At this resolution, the model is able to accurately resolve the effect of mountainous topography and reduce the number of parameterisations needed to model processes such as cumulus convection (Lader and others, 2020).

The original data used by Lader and others (2020) for dynamical downscaling consists of output from three models. One reanalysis model (CFSR; Saha and others, 2010), is used here for evaluation simulations of past SMB between 1981 and 2019. This was chosen by Lader and others (2020) as a target for dynamical downscaling as it is one of the top performing reanalysis models for southeast Alaska (Lader and others, 2016). Furthermore, output from the dynamically downscaled CFSR model shows good agreement with observations from the automatic weather station at Juneau airport, albeit with a slight positive precipitation bias because of the large range of topography within the gridcell of Juneau airport (Lader and others, 2020). Two dynamically downscaled global climate models (GFDL-CM3, Donner and others, 2011, and NCAR-CCSM4, Gent and others, 2011), provide both historical runs for 1981–2010, and projections for the period 2031–2060. These two models were chosen by Lader and others (2020) as they routinely rank in the top five of all CMIP5 models for Alaska in evaluations against near-surface temperature, precipitation and sea level pressure data (Walsh and others, 2018). Usefully, the two climate models also have different climate sensitivities – NCAR-CCSM4 has a low climate sensitivity, GFDL-CM3 has a high climate sensitivity (Flato and others, 2013). Climate sensitivity refers to the degree to which global temperatures respond to changes in greenhouse gas concentrations, with high sensitivity indicating greater temperature changes. Henceforth, we refer to GFDL-CM3 as GFDL and NCAR-CCSM4 as CCSM. The use of these specific model datasets is outlined in Figure 2.

The data obtained from Lader and others (2020) provides dynamically downscaled projections of climate under the RCP8.5 emissions scenario. Recent studies have suggested scenario RCP8.5 is unlikely, and cautioned its overuse (Hausfather and Peters, 2020; Burgess and others, 2021). However, as also highlighted by Lader and others (2020), since our projections are focussed only on the short-term future (to 2060), the choice of climate scenario alone may be inconsequential due to the atmospheric residence time of the greenhouse gases and aerosols already emitted, with many models in CMIP5 predicting similar increases in global air temperatures for all RCP emissions until 2060 (Overland and others, 2014).

Despite the dynamically downscaled products from Lader and others (2020) having a spatial resolution of 4 km, many glaciers of Juneau Icefield have a width of less than 1 km, creating the need for further downscaling or interpolation. For 2 m air temperature (T ₂) and surface pressure (p), nearest-neighbour interpolation was used to select the climate data gridpoint of each variable closest to the COSIPY gridpoint before applying the chosen downscaling method. T ₂ was downscaled using a constant lapse rate of 5 K km⁻¹, chosen through analysis of several AWS sites in McNeil and others (2020). Surface pressure was downscaled using the barometric formula. The precipitation, incoming longwave and shortwave radiation were all bi-linearly interpolated onto the COSIPY model grid. Wind speed was interpolated onto the model grid in the same manner, after being derived from wind speed components at 10 m, assuming a logarithmic wind profile (Arndt and others, 2021). The incoming shortwave radiation was additionally adjusted after being interpolated to the model grid by a radiation model based on Wohlfahrt and others (2016). The model corrects the incoming shortwave radiation by taking into account the effect of shadowing from nearby terrain at each grid-point throughout the year.

The final step in preparing the climate input was to conduct bias correction of the downscaled global climate model data (CCSM and GFDL). Since these models are only forced by radiative and aerosol forcings, they do not replicate past interannual variability well, but should capture longer-term trends. Here, we used empirical quantile mapping based on the method used in Amengual and others (2012). Empirical quantile mapping adjusts the cumulative distribution function of the future projected climate data (CCSM and GFDL) by adding the mean regime shift and deviation of the climate models to the quantiles of the observed dataset (CFSR reanalysis). Values of certain variables that are outside of the range of the observational dataset are corrected using constant extrapolation. Additionally, negative precipitation and radiation values are set to zero and relative humidity is bounded from 0 to 100. For simplicity, the wet day frequency and drizzle threshold of the bias-corrected climate data was not changed, and all variables were adjusted additively. Table S2 displays the mean and decadal trends of precipitation, snowfall, air temperature and incoming shortwave radiation for the raw and bias-corrected GFDL and CCSM data. This bias correction procedure preserved trends in the projection datasets, While shifting the model output towards the observed climatology and removing outliers. For both models, bias correction decreased future precipitation (rain and snowfall). For air temperature, the mean annual increased for CCSM, but decreased for GFDL. These adjustments reflect the opposing climate sensitivities of the two models, which were evident in the historical simulations when comparing the means between the global climate models and the reanalysis model (see Table S2).

2.3. Model optimisation

To optimise the COSIPY model, we use a calibrated glacier-wide mass balance dataset provided by the Juneau Icefield Research Program (Pelto and others, 2013; McNeil and others, 2020). This calibrated dataset is derived from midsummer accumulation and ablation observations, geodetic mass balance estimates, glacier hypsometry, and transient snowline observations to provide a comprehensive reanalysis of the mass balance of Taku Glacier from 1946 to 2023 and Lemon Creek Glacier from 1953 to 2023.The associated uncertainty with this mass balance is ± 0.45 m w.e. a⁻¹ (McNeil and others, 2020).

Initial parameter values for model optimisation were selected based on Blau and others (2021), with the albedo ranges adjusted not to overlap with each other. A Latin Hypercube was used to generate a random sample of 100 values for each of the 10 parameters, with the sampling evenly distributed across the entire parameter space. COSIPY was then run in an ensemble of 100 model simulations for Taku Glacier (which accounts for 19% of the total icefield area) from 2003 to 2008, forced with the CFSR reanalysis data (Fig. S1). The optimisation process was not conducted on Lemon Creek Glacier due to its smaller size and its location on the moister southwest side of the icefield. The model simulations were then scored against the observations of SMB from the Juneau Icefield Research Program using the root mean square error (RMSE) and the coefficient of determination (R ²). The optimal model parameters were selected from the simulation with the lowest RMSE (0.24 m w.e. a⁻¹). All the top 10 scoring runs captured the interannual variability of SMB well (Fig. S2) and had a RMSE of less than 0.54 m w.e. a⁻¹. Optimised model parameters are summarised in Table 1.

Table 1. Optimised model parameters

Note that the initial testing range is adapted from those tested in Blau and others (2021).

To evaluate the impact of each parameter on the modelled average annual SMB, a sensitivity analysis using the 100 simulations of Taku Glacier and a Random Forest regression model was conducted. The most influential parameter was the multiplication factor for total precipitation (Fig. S3). The same set of 100 simulations was used to estimate a conservative uncertainty related to the COSIPY parameters. We calculated the mean absolute error between the observed mass balance from the Juneau Icefield Research Program and the top 10 performing model optimisation runs, based on RMSE, for Taku Glacier. The mean absolute error was found to be ± 0.38 m w.e. a⁻¹. We applied this uncertainty measure to our historical and projected icefield-wide model simulations, assuming that the parametric uncertainty from the 100-model optimisation runs for Taku Glacier is representative of the model performance for the entire icefield. Furthermore, we presume this uncertainty remains consistent when using different climate models, which have been bias-corrected against the CFSR reanalysis model used in the optimisation simulations. While this approach is simplistic, conducting a full Monte Carlo uncertainty analysis for the entire icefield is computationally prohibitive.

3. Results

3.1. Comparison of modelled and observed mass balance data for Taku and Lemon Creek glaciers

The optimised parameters were used to conduct a longer model simulation for both Taku and Lemon Creek glaciers (1981–2019), to investigate whether the model would produce realistic results outside of the optimisation period (Fig. 3). The model matches the Juneau Icefield Research Program observations for Taku Glacier in and outside the optimisation period, with the COSIPY simulations matching the observed interannual variation (Fig. 3a). The model appears to fit less well for the Lemon Creek Glacier, likely owing to the glaciers' relatively small area, which is less than three grid cells of the dynamically downscaled climate data. A slight positive bias was noted for the simulations of both glaciers possibly due to the model using a shapefile of the glaciers 2019 extent and missing grid points in the ablation areas.

Figure 3. Comparison between the optimised COSIPY SMB (red) and the calibrated Juneau Icefield Research Program SMB (blue) for (a) Taku Glacier and (b) Lemon Creek Glacier, between 1981 and 2019 using the CFSR reanalysis data. The bottom left box on each subplot displays the RMSE and R2 between the COSIPY model and the Juneau Icefield Research Program SMB data.

3.2. Evaluation simulations of surface mass balance from reanalysis (1981–2019)

Across the Juneau Icefield, the evaluation simulation from CFSR displayed a non-statistically significant (Mann–Kendall test, p-value = 0.24) warming trend of + 0.02°C decade⁻¹, with a mean annual air temperature of −3.38°C (Table S2). Elevation was the strongest control upon temperature, with a lapse rate of 6 K km⁻¹ positioning the lowest mean annual air temperatures in the centre of the plateau (≃8°C) and the highest at low elevation tongues of outlet glaciers (Fig. 4a). The mean annual total precipitation across the icefield is 3065 mm a⁻¹, with a maximum occurring in the southwest and a minimum occurring at the tongue of Llewellyn Glacier (Fig. 4b). Snowfall follows a similar spatial pattern, the difference being a more concentrated maximum in the southwest portion of the icefield, and little to no snowfall at the tongues of outlet glaciers due to higher temperatures (Fig. 4c). Both precipitation and snowfall have a non-statistically significant (Mann–Kendall test, p-value = 0.07 and 0.06, respectively) negative trend during the historic period (Table S2). The mean annual incoming shortwave radiation across the icefield is 110 W m⁻² and shows a strong gradient, with lower values in the southwest (95 W m⁻²) and higher values in the northeast (~130 W m⁻²) (Fig. 4d).

Figure 4. Climate averages for the Juneau Icefield from CFSR-COSIPY simulations. (a) Mean daily air temperature at 2 m, (b) mean total annual precipitation, (c) mean annual snowfall, (d) mean incoming shortwave radiation.

The SMB of the whole icefield for the historic period is summarised in Figure 5a. The reanalysis-driven simulation (CFSR-COSIPY) showed an average annual SMB of −0.22 ± 0.38 m w.e. a⁻¹, with a statistically significant negative trend (Mann–Kendall test, p-value = 0.001) of −0.02 m w.e. a⁻². Within the accumulation area of the plateau, mass balance is often greater than 2 m w.e. a⁻¹, while at the low elevation ablation zone of the glaciers the mass balance is below −4 m w.e. a⁻¹. The model simulates that most of the select glaciers in Fig. 5b (Mendenhall, Lemon Creek, Taku, Llewellyn, Field, Willison, Meade and Tulsequah glaciers) had a positive SMB at the start of the evaluation period. An exception to this is the Gilkey Glacier, which had a negative SMB throughout the majority of the simulations (Fig. 5b). The annual SMB of the icefield began to turn negative ~ the late 1980s, as does the annual SMB of the Mendenhall, Gilkey, Field, Meade and Lemon Creek glaciers.

Figure 5. Surface mass balance estimates from CFSR-COSIPY during the historic period (1981–2019). (a) The mean annual surface mass balance. (b) The annual cumulative surface mass balance of select glaciers across Juneau Icefield.

The average icefield-wide ELA was ~1360 m, with individual glacier ELAs ranging between 1150 to 1640 m. The average accumulation area ratio for this period is 56%, and covers the majority of the plateau. The spatial pattern of simulated ELA using the CFSR simulations, compares favourably with ELA observations (Ziemen and others, 2016; Fig. S6). The lowest ELAs occurred at the southerly glaciers (e.g. Taku, Mendenhall and Tulsequah glaciers), because of the higher accumulation rates at these glaciers. Across the nine major glaciers identified in Figure 5b, an average positive trend in the ELA of 6.25 m a⁻¹ was observed. This trend was statistically significant (revealed by a Mann–Kendall test) for all glaciers except Field Glacier. For several years in the past simulations, Lemon Creek Glacier's ELA was incalculable, as there was no accumulation area and all grid points had a negative annual SMB.

Table 2 compares previously derived mass balance estimates from the Juneau Icefield, using both geodetic (Larsen and others, 2007; Berthier and others, 2010, 2018; Melkonian and others, 2014) and modelling studies (Ziemen and others, 2016; Young and others, 2021), against the CFSR simulation to provide context for where our estimates fit within the broader range of studies. The CFSR specific mass balance is consistent with other studies, albeit with a slight positive bias, showing a negative icefield-wide mass balance for the 20th century. In all cases, the CFSR results fall within the observational uncertainties of other studies and our own estimated uncertainty. Per-glacier simulations also show a slight positive bias when compared to mass balance estimates for nine glaciers derived from laser altimetry in Larsen and others (2015) (Table S1). These differences between geodetic estimates and the CFSR simulation are likely due to disparities in the icefield area, with our results use a fixed area throughout the simulations, and our simulations not accounting for mass loss from calving processes. The differences between modelling studies and the CFSR likely reflect the use of different mass-balance models, approaches to model calibration, and methods for downscaling climate data.

Table 2. Comparison of the specific mass balance results from this study (CFSR-COSIPY) with previous estimates

Note the different study periods for each estimate. The area of the icefield is not listed for Ziemen and others (2016) due to its evolution with modelled ice dynamics. Additionally, the estimate from Melkonian and others (2014) is likely an underestimate due to biases from radar penetration depths in different seasons, as noted by Berthier and others (2018). The CFSR-COSIPY estimate has been subset to match the time period of other studies.

3.3. Historical simulations of surface mass balance from climate models (1981–2010)

The evaluation simulation presented in Section 3.2 was driven by climate reanalysis data (CFSR) that is constrained by observations, and thus represent our most realistic output. However, to obtain projections under future climate, observationally unconstrained simulations are required, forced by scenarios of future atmospheric greenhouse gas concentrations. To ascertain whether the two global climate models studied here, GFDL and CCSM, were able to simulate realistic patterns of climate pertinent to mass balance processes, we first conducted simulations in COSIPY driven by the modelled historical climates from these two models covering the period 1981–2010. These historical climate simulations were forced by past atmospheric greenhouse gas concentrations, solar variability, and volcanic forcing.

The means and decadal trends of total annual precipitation, total annual snowfall, annual air temperature, and annual incoming shortwave radiation for the GFDL, CCSM and CFSR simulations are summarised in Table S2. Spatial averages for GFDL-COSIPY and CCSM are depicted in Figs S4 and S5, respectively. The historical simulations from GFDL and CCSM exhibit higher mean annual total precipitation (3450–3630 mm a⁻¹) and snowfall (2.01–2.05 m w.e. a⁻¹) compared to CFSR (3070 mm a⁻¹ for precipitation and 1.80 m w.e. a⁻¹ for snowfall). Mean annual incoming shortwave radiation was relatively similar among the models, with values of 109, 102 and 110 W m⁻¹ for GFDL, CCSM and CFSR, respectively. The GFDL model, which has higher climate sensitivity, shows more pronounced decadal trends, with statistically significant (as revealed by a Mann–Kendall test) increases in temperature, and decreases in annual precipitation and snowfall over the period 1980–2010. In contrast, the lower-sensitivity CCSM model only shows a significant decrease in snowfall during the same period (Table S2). This led to the two models producing a spread of simulated SMB over the historical period (Fig. 6). As expected, the climate models do not match the interannual variability of the evaluation simulations; however, the average annual SMB for the historical period from each model (GFDL, CCSM and CFSR) overlaps within their estimated uncertainties. If the mean of the two models (GFDL and CCSM) is considered, the SMB is within ± 0.03 m w.e. a⁻¹ of the CFSR evaluation simulation (Fig. 6). Furthermore, the spatial distribution of SMB from the mean of the two model runs, matches closely with the CFSR-driven simulations (Fig. S7). This gives confidence in the two models for providing projections of future SMB.

Figure 6. Time series of SMB across Juneau Icefield for 1981 to 2010 for the CFSR-COSIPY evaluation simulation (black), and the two global climate models (GFDL-COSIPY – purple, and CCSM-COSIPY – green). The mean SMB of GFDL-COSIPY and CCSM-COSIPY is shown in red. The average annual SMB of the whole period for each model is noted in the legend (m w.e. a⁻¹).

3.4. Future surface mass balance of the Juneau Icefield (2031–2060, RCP 8.5)

The change in modelled climate between the means of the historical (1981–2010) and future bias-corrected (2031–2060) simulations of each model is shown in Figure 7. Both models project an increase in temperature everywhere on the icefield, with a Mann–Kendall test of the annual averages showing a statistically significant trend for the future period at most grid points, except for at some cells on glacier tongues. The GFDL simulation predicts a mean increase of 3.50°C, While CCSM predicts a mean increase of 1.01°C. There are spatial differences between the two models. The GFDL simulation shows the largest temperature increases at the tongues of glaciers in the north and northeast, most notably at Llewellyn Glacier (Fig. 7a), While the CCSM indicates that the tongues of glaciers in the south and west will have the largest temperature increases, most notably at Taku Glacier (Fig. 7e). Both models project a small decrease in annual total precipitation almost everywhere on the icefield (Figs 7b and f). This is specific to the icefield domain, as precipitation across the rest of southeast Alaska is projected to have a positive trend under both models due to a deepening of the Aleutian Low (Gan and others, 2017; Lader and others, 2020). Both models also project a reduction in snowfall across Juneau Icefield (Figs 7d and h). In both cases, the change in snowfall is greater than the total change in precipitation, indicating that a smaller proportion of precipitation will fall as snow. Little change is projected to occur in incoming shortwave radiation, though both models show a slight increase at the tongues of Llewellyn Glacier and Tulsequah Glacier.

Figure 7. Changes in the future climate averages for the Juneau Icefield domain between the future (2031–2060; RCP8.5) and historic (1981–2010) means. (a) Change in daily 2 m air temperature, (b) percentage change in total annual precipitation, (c) percentage change in daily incoming shortwave radiation and (d) percentage change in annual snowfall for the GFDL-COSIPY data. (e–h) same as (a–d) but for the CCSM-COSIPY data. Additional text boxes on each plot show the average change for the whole icefield domain.

The projected future (2031–2060) SMB under the RCP8.5 emissions scenario is shown in Figure 8. Both models display an increasingly negative annual SMB despite their different climate sensitivities (Fig. 8). These simulations compare well with the RCP8.5 glacier mass loss projections for Alaska from Hock and others (2019). This similarity occurs despite our simulations lacking a consideration of changing glacier hypsometry due to ice flow. The mean annual SMB of the two models of −1.52 ± 0.27 m w.e. a⁻¹ (2031–2060), represents a decrease of 1.41 m w.e. a⁻¹ when compared to the mean annual SMB of the historical simulations (1981–2010). The mean trend in the annual SMB for the average of the two models over 2031–2060 was −0.064 m w.e. a⁻², which was found to be statistically significant (Mann–Kendall test, p = 0.0001). When considering the mean of the two models, there are multiple years towards the end of the simulation where no grid points on the icefield have a positive SMB (not shown). The largest decreases in SMB are projected to be in the accumulation areas on the icefield plateau, particularly in the southwest catchment zones of Norris, Mendenhall and Taku glaciers (Fig. 8c). The predictions suggest that glaciers in the south (e.g. Mendenhall, Taku and Lemon Creek glaciers) will have a SMB ~ 1.0 m w.e. a⁻¹ lower than elsewhere on the icefield. The glaciers with the lowest SMB by the end of the simulations are Mendenhall and Lemon Creek glaciers.

Figure 8. Projections of future SMB of Juneau Icefield (2031–2060, RCP8.5). (a) Time series of annual SMB across Juneau Icefield for 2031–2060 for the two global climate models (GFDL-COSIPY – purple, and CCSM-COSIPY – green). The mean annual SMB of GFDL-COSIPY and CCSM-COSIPY is shown in red. The average annual SMB of the whole period for each model is noted in the legend (m w.e. a⁻¹). Regional RCP8.5 projections of specific mass balance for Alaska from Hock and others (2019) are displayed in grey. (b) The spatial distribution of the mean annual SMB from the mean of GFDL-COSIPY and CCSM-COSIPY. (c) The change in mean annual SMB compared to the GFDL-COSIPY and CCSM-COSIPY simulated mean from the historic period.

The surface mass balances shown in Figure 8, averaged from both models, indicate an increase in ELA of ~370 m by 2060, with the highest increases (~400 m) seen at Field, Willison and Meade glaciers in the north. The projections reveal that the ELA of Lemon Creek Glacier remains above the elevation of its highest gridpoint. The accumulation area ratio of glaciers across the icefield was found to decrease from 54% in the mean of the GFDL and CCSM historic simulations to 18% in projections (2031–2060), and to only 6% when considering the final decade of the simulations (2051–2060). The most rapid rise in ELA occurs during the first 10 years of the simulation, with the icefield wide accumulation area ratio becoming ~10%. By 2060, accumulation is confined to the highest elevations of the icefield.

Coupled with these changes in elevation and accumulation areas is a notable shift in the melt season duration. Historically, the melt season predominantly spanned from early May to late August, based on the mean of the GFDL and CCSM data from 1981 to 2010. However, future projections suggest an earlier onset starting in late April and extending into mid-September, lengthening the melt season by approximately 24 days (GFDL-CCSM mean, 2031–2060) (Fig. 10b). If comparing to the last 5 years of the future projection (2056–2060), the melt season increases by 38 days when compared to the historic simulations. This prolonged melt season is accompanied by a 39% increase in the total annual surface melt (GFDL-CCSM mean) from the past to the future simulations.

4. Discussion

Our projected mass balance simulations utilise the RCP8.5 emissions scenario. As outlined in Section 2.2, although this high-end emissions scenario is now considered unlikely (Hausfather and Peters, 2020; Burgess and others, 2021), the atmospheric residence times of greenhouse gas emissions render the choice of scenario largely inconsequential for the timeframe of our study (up to 2060) (Overland and other, 2014). This assertion is supported by global mass balance projections from Hock and others (2019), which indicate that the differences in glacier mass loss between the RCP2.6, 4.5 and 8.5 emissions scenarios for Alaska are minimal until approximately 2060 (Hock and others, 2019, Fig. 4). Similar findings were also reported in regional projections of ice volume loss of western Canada by Clarke and others (2015), with small variations between RCP2.6, 4.5, 6.0 and 8.5 scenarios. This context allows us to postulate on the potential fate of the Juneau Icefield, which is already committed to significant mass loss, while also exploring a range of possible responses of the icefield to climate change based on the two contrasting climate sensitivities of the global climate models used (Flato and others, 2013).

The following discussion is presented with an additional important context. Our model simulations assume a fixed glacier geometry and do not account for the impacts of ice dynamics, such as glacier thinning and retreat. To a certain extent these two processes have opposing effects on SMB, potentially reducing their overall impact on our projections (Schaefer and others, 2013). To fully explore the impact of these processes on SMB, it would be necessary to couple COSIPY with an ice-flow model, which is out of the scope of this study. However, the results of our SMB modelling allow us to postulate on the mechanisms of mass loss and the potential fate of Juneau Icefield.

4.1. Drivers of change across the Juneau Icefield

To investigate the potential drivers of the projected decrease in SMB (Fig. 8), we divided the annual SMB from the future simulations into ablation and accumulation components (Fig. 9a). Both models show a strong negative correlation between ablation and annual SMB (Fig. 9b), indicating that interannual variability and the negative trend in SMB are strongly linked to the variability of ablation processes. This increase in ablation is primarily associated with rising air temperatures (Fig. 9c), with a multi-model (GFDL and CCSM) mean summer temperature increase of 3.03°C for 2030–2060 compared to 1980–2010. The stronger correlation between mean annual air temperature and SMB in future projections (Fig. 9c), compared to the past, suggests the increasing influence of rising air temperatures on the negative SMB of the icefield.

Figure 9. (a) Time series of the GFDL-COSIPY and CCSM-COSIPY future simulations, separating SMB, ablation, and accumulation. Note that the sign of ablation has been swapped, so that a decrease on the graph represents more ablation. (b–g) Scatter plots showing the relationships between different variables and the ablation and accumulation across the icefield for the CFSR-COSIPY evaluation and GFDL-COSIPY and CCSM-COSIPY projections. The Pearson correlation coefficient (r) is shown for each.

There is also a strong positive correlation between accumulation and annual SMB (Fig. 9e), with higher accumulation years generally leading to more positive SMB (Fig. 9a). In the GFDL projections, unlike the historical evaluation simulations (CFSR), high-accumulation years did not mitigate ablation by increasing the mean albedo of the icefield (Fig. 9a) (Moore and others, 2009; Schaefer and others, 2013). This lower correlation between accumulation and mean icefield albedo in future GFDL projections (Fig. 9g), compared to the past, suggests that the increasing confinement of snowfall to higher elevations limits its ability to enhance the icefield's overall albedo, even in years of high snowfall.

Rising air temperatures also affect the ratio between solid (e.g. snow) and liquid (e.g. rain) precipitation, with a shift in the ratio towards rain projected for Alaska (McAfee and others, 2014). Historical simulations indicate only a weak negative correlation between air temperature and accumulation (Fig. 9f). However, future projections show a moderate negative correlation, with higher air temperatures leading to more precipitation to fall as rain rather than snow (Figs 10e, d). This shift in precipitation phase is most pronounced during autumn, a period when the Juneau Icefield typically receives the greatest snowfall (Fig. 10d). Historically, snowfall occurs year-round and starts to increase in August, peaking in October (Fig. 10d). This autumn peak is due to the thermal contrast between the ocean and the land being at its highest around October, resulting in a deeper Aleutian Low, more frequent storms and large influxes of moisture (Wendler and others, 2016). Similar to the findings of McGrath and others (2017), our future projections suggest that daily snowfall will not increase until mid-to-late September, with an autumn peak noticeable only for the period 2031–2035 (Fig. 10d). This occurs despite an intensification of the Aleutian Low in both models, due to the increasing thermal contrast between the ocean and land (Gan and others, 2017).

Figure 10. The historical and projected changing seasonal cycle across Juneau icefield from the mean of the GFDL-COSIPY and CCSM-COSIPY simulations. (a) Average daily air temperature, (b) total daily surface melt, (c) average daily albedo, (d) total daily snowfall and (e) total daily total precipitation. A 30-day running mean has been applied to all fields.

The decreasing SMB of the Juneau Icefield is primarily driven by rising air temperatures, which increase ablation and reduce the amount of precipitation falling as snow, thereby decreasing accumulation. We hypothesise that these drivers are further amplified by feedback mechanisms such as the ice-albedo and ice-elevation feedbacks, which exacerbate melt. The temperature-driven ice-albedo feedback is illustrated in Figure 10, where increased temperatures (Fig. 10a) lead to higher surface melt rates (Fig. 10b), exposing more low-albedo ice surfaces to melting (Fig. 10c). This, in turn, increases melting due to the increased absorption of incoming solar radiation. This ice-albedo feedback mechanism was demonstrated in Johnson and Rupper (2020) to be responsible for up to 80% of the resultant melt caused by a 1°C rise in air temperature. Over recent decades, decreasing albedo has already been observed to significantly contribute to the ice volume loss of the Juneau Icefield (Davies and others, 2024). With the multi-model mean summer air temperature projected to rise by 3.03°C from 1980–2010 to 2030–2060, the importance of this mechanism is expected to increase significantly. This is evidenced by a stronger correlation between ablation and albedo in future simulations compared to historical simulation (Fig. 9d). Figure 10c further illustrates the albedo decrease, with each successive 5-year period from the future simulations, approaching the bare ice albedo of 0.45 in COSIPY.

The GFDL-COSIPY and CCSM-COSIPY projections indicate a significant rise in the ELA across the plateau, shrinking the mean accumulation area of the icefield to just 3.6% for the last 5 years of the future simulations. This leads to extensive thinning on the plateau, which in turn lowers the surface elevation and exposes it to warmer air temperatures and will likely increase the rain–snow ratio (O'Neel and others, 2015), thus exacerbating the thinning process and triggering ice-elevation feedback (Bodvarsson, 1955). This can prompt a nonlinear response to mass balance, with the icefield initially responding steadily to ELA rises, but then nonlinear once the ELA reaches the low slope plateau where minor increases dramatically decrease the accumulation area and cause widespread thinning (McGrath and others, 2017; Hock and Huss, 2021; Bolibar and others, 2022; Davies and others, 2024). Recent observations have already documented this thinning of the high-elevation plateau (Davies and others, 2024), and our projections suggest that this trend will persist. Eventually, the plateau surface may cease to serve as a source of snow and ice for outlet glaciers. In certain areas, the ice thickness on the plateau exceeds 600 m (Millan and others, 2022), indicating that this nonlinear mass balance feedback can operate over a substantial elevation range.

4.2. Potential response of glaciers across the Juneau Icefield

The top-heavy hypsometry of the Juneau Icefield with its high-elevation plateau which feeds multiple outlet glaciers through steep icefalls (Fig. 1) will likely influence the response of individual glaciers to the projected decrease in SMB (Furbish and Andrews, 1984; Jiskoot and others, 2009). These icefalls, which flow steeply over topographic steps, are prone to thinning and disconnection (Rippin and others, 2020; Davies and others, 2022). Once disconnected these outlet glaciers will lose their accumulation area, causing a large decrease in their mass balance.

The average annual SMB at the icefall connections of the 13 outlet glaciers to the main high-elevation plateau (Davies and others, 2022) was −1.31 m w.e. a⁻¹ (2031–2060, GFDL-COSIPY and CCSM-COSIPY mean). Of further concern are the East and West Twin glaciers, with icefalls already observed to be thinning and narrowing in 2019 (Hugonnet and others, 2021; Davies and others, 2022). Our simulations of future SMB indicate that the mean annual SMB at these icefalls is −3.20 and −2.20 m w.e. a⁻¹ for the East and West Twin glacier icefalls, respectively (Fig. 11b). These high thinning rates suggest a glacier disconnection is likely if temperatures rise as projected. Furthermore, our SMB model projects widespread icefield thinning, potentially creating new icefall locations over steep topographic steps (Rippin and others, 2020; Davies and others, 2022). Thus, glacier disconnection and rapid down-wasting may become increasingly common in the future.

Figure 11. Potential future and ongoing changes to outlet glaciers of Juneau Icefield. (a) and (b) Mean annual SMB for 2050–2060 from the GFDL and CCSM simulations superimposed and interpolated on Google-Earth imagery, for key icefalls on (a) Herbert Glacier and (b) West and East Twin glaciers. Here, the annual mass balance is projected to be highly negative. These areas are thus likely locations of future glacier detachment. (c) Photo of West Twin Glacier icefall (July 2022), a location already undergoing thinning. (d) The terminus of Mendenhall Glacier (July 2022). Note how the front has already partially receded from its proglacial lake.

Isolated mountain glaciers, such as Lemon Creek Glacier, are not connected to the main high-elevation plateau of the icefield. These glaciers face similar challenges as those that have lost their accumulation areas through disconnection at icefalls. Without access to the large catchment area of the main plateau, which can provide additional ice to slow their retreat, these isolated glaciers are prone increased ablation rates. While Bolibar and others (2022) suggest that the retreat of these glaciers to higher elevations might mitigate losses due to future warming, our historical simulations indicate that Lemon Creek may already have an ELA near or above its highest altitude. In such cases, climatic forces likely surpass any mitigating effects from topographic feedbacks from Lemon Creek and other smaller mountain glaciers in the Juneau Icefield region.

Conversely, tidewater glaciers are likely to respond differently under our SMB projections. Historically, the formerly tidewater Taku Glacier has been an outlier, as it was the only glacier to have advanced during the 20th century (Post and Motyka, 1995; Pelto and others, 2013; McNeil and others, 2020). However, between 2013 and 2018, Taku Glacier began to retreat (McNeil and others, 2020). As the annual SMB of Taku becomes increasingly negative, as shown in both of our future projections, the glacier's apparent resilience to climate change observed in the 20th century will cease. This will likely prompt a retreat into its 40 km long over-deepened basin (Nolan and others, 1995) and the eventual re-initiation of calving (McNeil and others, 2020). The retreat of the icefield's outlet glaciers through increased ablation of their tongues is likely to alter their terminus type and subsequently affect ice dynamics and frontal ablation. For outlet glaciers with low-gradient or over-deepened beds, such as Tulsequah, Gilkey, Meade and Field (Davies and others, 2022), this will likely lead to the formation of proglacial lakes, increasing frontal ablation through calving (Motyka and others, 2003; Boyce and others, 2007; Davies and others, 2022). In contrast, other outlet glaciers with steeper slopes, like East Twin Glacier, may retreat beyond the water bodies they flow into. The beginning of this transition can already be observed at Mendenhall Glacier, where retreat has resulted in half of the terminus being situated on land above the proglacial lake (Fig. 11d).

5. Conclusions

High-resolution dynamically downscaled climate models (Lader and others, 2020) have enabled us to simulate the historical and projected future SMB of Juneau Icefield using the COSIPY model. Tuning of the model to the rich empirical record collected by the Juneau Icefield Research Project enabled us to accurately simulate the pattern of past changes and provides confidence in future projections. This highlights the value of such long-term monitoring programmes. Under RCP8.5 projections, our modelling suggests that the mass balance across Juneau Icefield is set to become increasingly negative in the middle of the 21st century, with a dramatic rise in ELA and reduction of the accumulation area. For the period 2031–2060, a multi-model mean SMB of −1.52 ± 0.27 m w.e. a⁻¹ is projected. This is attributed to an icefield-wide increase in air temperature, which causes increased snowmelt and a higher percentage of precipitation to fall as rain rather than snow. Reduction of snow cover in the model leads to longer and more extensive exposure of lower albedo ice, leading to an albedo-induced melt feedback. Negative mass balances are likely to spread across the plateau, and at icefalls ice thinning is likely to promote glacier disconnections. These stark projections of future mass balance are likely to lead to numerous feedbacks which augment future ice losses from the Juneau Icefield. Future work should consider ice-flow feedbacks, and the timescales over which similar processes are likely to occur on other plateau-icefields and ice caps globally.