1. Introduction
Glaciers play an important role in a wide range of both human and ecological systems (O’Neel and others, Reference O’Neel2015; Immerzeel and others, Reference Immerzeel2020). Global observations show that glaciers have lost mass throughout the 20th and 21st centuries, and mass loss from glaciers in Alaska and northwest Canada (the Randolph Glacier Inventory (RGI) ‘Alaska region’; (RGI Consortium, 2017)) outpaced any other region (Larsen and others, Reference Larsen, Burgess, Arendt, O’Neel, Johnson and Kienholz2015; Zemp and others, Reference Zemp2019; Hugonnet and others, Reference Hugonnet2021). Global glacier mass loss will continue through the end of the century, and glaciers in the Alaska region are projected to provide the largest contribution to global sea level rise in the 21st century of any glacierized region on Earth (Rounce and others, Reference Rounce2023).
Repeated, systematic observations of glacier mass balance provide critical information for understanding how glaciers respond to a changing climate (Zemp and others, Reference Zemp, Hoelzle and Haeberli2009). However, fewer than 0.1% of glaciers in the Alaska region have consistent and long-term in situ observations (e.g. U.S. Geological Survey Benchmark Glacier Program and others, 2016). Geodetic remote sensing methods offer a complementary approach to collecting direct measurements, providing insight into glacier mass balance across local to global scales (Jourdain and others, Reference Jourdain2023; Zeller and others, Reference Zeller, McGrath, Sass, O’Neel, McNeil and Baker2023). However, geodetic mass balance approaches are affected by a variety of potential biases related to limitations on glacier area change datasets (Florentine and others, Reference Florentine, Sass, McNeil, Baker and O’Neel2023), different processing steps in the elevation data processing chain (Piermattei and others, Reference Piermattei2024) and the time intervals over which they can be applied (Huss, Reference Huss2013).
The line separating snow surfaces from ice or firn surfaces on a glacier is the snowline, and the average altitude of this boundary at any instant, particularly during the ablation season, is referred to as the transient snowline altitude (SLA) (Cogley and others, Reference Cogley2011). For glaciers with a distinct accumulation and ablation season (such as the Alaska region), the SLA at the end of the ablation season can be used as a proxy for the annual equilibrium line altitude (ELA) (Rabatel and others, Reference Rabatel2012). The area which remains snow covered at the end of the ablation season represents the glacier’s accumulation area, which can be used in conjunction with the total glacier area to calculate the glacier’s accumulation area ratio (AAR) (Cogley and others, Reference Cogley2011). Furthermore, this snow covered area on a glacier surface is readily observable from space-based platforms by differentiating the areas of a glacier which are covered in snow from those which are ice or firn, because snow has higher reflectance across visible and near-infrared wavelengths (Gao and Liu, Reference Gao and Liu2001).
The ELA and AAR provide insight into the relationship between local climate and glacier mass balance on short (annual) time scales as variations in the ELA and accumulation area respond to changing patterns and magnitudes of snow accumulation and ablation. Remotely sensed observations of these variables can be used for translating long-term geodetic mass balance measurements to annual time steps (Davaze and others, Reference Davaze, Rabatel, Dufour, Hugonnet and Arnaud2020), investigating glacier vulnerability and disequilibrium (McGrath and others, Reference McGrath, Sass, O’Neel, Arendt and Kienholz2017; Carturan and others, Reference Carturan, Rastner and Paul2020; Lorrey and others, Reference Lorrey2022) and calibrating glacier models (Geck and others, Reference Geck, Hock, Loso, Ostman and Dial2021).
Many previous studies have explored approaches to extracting the SLA from optical satellite imagery, with most focused on identifying the snow covered areas from Landsat imagery using methods that range from manual delineation to fully automated classification approaches (Rabatel and others, Reference Rabatel2017). These studies have focused on glaciers in, for example, the Tropics (Rabatel and others, Reference Rabatel2012), the European Alps (Rabatel and others, Reference Rabatel, Dedieu, Thibert, Letréguilly and Vincent2008; Rastner and others, Reference Rastner2019; Prieur and others, Reference Prieur, Rabatel, Thomas, Farup and Chanussot2022), Canadian Rockies (Jiskoot and others, Reference Jiskoot, Curran, Tessler and Shenton2009) and High Mountain Asia (Racoviteanu and others, Reference Racoviteanu, Rittger and Armstrong2019; Guo and others, Reference Guo, Geng, Shen, Wu, Chen and Wang2021) or a global sample of glaciers (Li and others, Reference Li, Wang and Wu2022), but very few have utilized the higher spatial and temporal resolution of more recently available Sentinel-2 imagery (Kavan and Haagmans, Reference Kavan and Haagmans2021).
These previous studies tend to be limited to a small sample size of glaciers (fewer than 250), focus on regional-scale SLAs rather than glacier-specific SLAs and/or have not made the derived datasets accessible for use in other studies. Additionally, few studies have explored SLA and snow cover extraction in the Alaska region, and those that have are focused on a small sample size or regional approximations of snowline variability from coarser-resolution MODIS imagery (Pelto, Reference Pelto2011; Mernild and others, Reference Mernild, Pelto, Malmros, Yde, Knudsen and Hanna2013; Shea and others, Reference Shea, Menounos, Moore and Tennant2013).
To evaluate the current state and future vulnerability of glaciers across the entire Alaska region, we create a high resolution inventory of 2018–22 end-of-summer accumulation areas and ELAs for over 3000 glaciers (
$ \gt 76,000\;\mathrm{km}^2$) using an automated classification of Sentinel-2 imagery (Zeller and others, Reference Zeller, McGrath, Sass, Florentine and Downs2024). Using this dataset, we assess the current state of Alaska glacier accumulation areas and investigate the climatic drivers of ELA variability. Lastly, we project how ELAs will change by the end of the century using an ensemble mean of 13 global climate models and estimate the resulting loss of glacier accumulation area across the region.
2. Study area
We derived accumulation areas and ELAs for glaciers defined by the RGI version 6.0 (RGI Consortium, 2017). All glaciers in the Alaska region larger than 
$2\,\mathrm{km}^2$ were evaluated (Fig. 1). This includes 3031 individual glaciers (11% of the 27,108 glaciers in the region) which together cover 
$76,569\,\mathrm{km}^2$ (
$ \gt 88\%$ of the total glaciated area). Glaciers smaller than 
$2\,\mathrm{km}^2$ were excluded from our analysis as their smaller size increases the relative effects of misclassification, and local terrain and topography provide more dominant controls on their patterns of snow accumulation and ablation (Florentine and others, Reference Florentine, Harper, Fagre, Moore and Peitzsch2018). Glaciers in the Brooks Range (of which 31 are larger than 
$2\,\mathrm{km}^2$) were excluded due to their short ablation seasons and limited mass turnover when compared to those at lower latitudes (Wendler and others, Reference Wendler, Fahl and Corbin1972).
Five second-order subregions (‘O2Regions’) of the first-order Alaska region, as defined by RGI Consortium 2017, encompass our study area: the Alaska Peninsula, Alaska Range, West Chugach Mountains, Saint Elias Mountains and North Coast Ranges (these O2Region boundaries are unchanged in the recently updated RGI version 7.0). We further divided these into 16 smaller subregions (‘O3Regions’), largely following the regions used by Kienholz and others Reference Kienholz, Herreid, Rich, Arendt, Hock and Burgess(2015). Glaciers are clustered mainly along the coast of the Gulf of Alaska in a maritime climate. Further from the coast, glaciers in the Alaska Range are generally smaller and experience a more interior climate with colder temperatures and less precipitation.
The study area in Alaska and northwest Canada. Glaciers included in the study are mapped, with colors corresponding to their O3Region. The inset legend indicates the color, name and labeled number of each O3Region. Red outlines and corresponding names indicate the O2Regions defined by the RGI (RGI Consortium, 2017). Inset globe and rectangle show the global context of the area. This map and all other figures are presented in the Alaska Albers equal-area projection (EPSG:3338).
3. Methods
The snow cover distribution on glacier surfaces was delineated in all 2018–22 Sentinel-2 Level-1C (top of atmosphere reflectance) multispectral imagery during the May–November months using a random forest classifier. After preliminary filtering and infilling of data gaps, each glacier’s SLA was extracted for each observation by identifying the elevation band at which the glacier surface transitioned from majority snow-free (ablation zone) to majority snow-covered (accumulation zone). The maximum SLA was identified from each annual time series of SLAs and taken to represent the ELA. This resulted in a 10 m resolution annual accumulation area product for each glacier each year, along with the annual ELA. The five annual accumulation area products for each glacier were then combined to produce a single average product, from which the 5 year average ELA was calculated. Each of these steps are described in detail in the following sections, and an overview workflow diagram in provided in the Supplementary material (Fig. S1).
3.1. Random forest classification
A random forest classifier (Breiman, Reference Breiman2001) was created to classify Sentinel-2 top of atmosphere reflectance images (Level-1C) of glaciers into six classes: snow, firn, ice, debris/rock, water and shadow. We created a training dataset in Google Earth Engine (Gorelick and others, Reference Gorelick, Hancher, Dixon, Ilyushchenko, Thau and Moore2017) from eight Sentinel-2 datatakes (the continuous image strip acquired by a single Sentinel-2 satellite, hereafter called images), with images selected to represent the full spatial and temporal range of the study area. Polygons of each class were manually digitized on each of these eight images. All image bands were resampled to 10 m resolution, and 20,000 pixels of each class were randomly chosen from each image and band values extracted, resulting in 16,0000 pixels of each class (96,0000 total).
In addition, a representative snow-on mosaic of the study region was created by combining all late-winter (between 19 February and 30 April) cloud-free Sentinel-2 imagery for 2018–22. Each image’s band reflectance was normalized relative to this snow-on mosaic, and these normalized band values were also used as inputs to the random forest classifier. This snow-on normalization was used as a way to provide a topographic correction for variable reflectances caused by variations in the solar incidence angle on the glacier’s surface (Fig. 2). The normalized difference water index (Gao, Reference Gao1996) and normalized difference snow index (Dozier, Reference Dozier1989) were calculated for each observation and included as features in the classification (Table S1).
An example of a Sentinel-2 image from which the training data were collected, showing Wright Glacier (RGI60-01.02602) (a), with example snow, firn and ice regions labeled in (b). Boxes in (a) indicate the extent of (b), (d) and (e). (c–e) show the process of normalizing imagery by the snow-on mosaic. (c) shows the snow-on mosaic of Wright Glacier, (d) shows a subset of the original image and (e) shows the effect of snow-on normalization. All images are displayed using the near-infrared, red and green bands.
The random forest classifier was trained using the scikit-learn RandomForestClassifier() implementation (Pedregosa and others, Reference Pedregosa2011) to classify pixels as one of the six classes (snow, firn, ice, debris/rock, water and shadow). A leave-one-out cross-validation approach was used during the training process to identify the optimal hyperparameters (namely, 
$n\_estimators$ and 
$max\_depth$) for the model to give the best results and to avoid overfitting to the training data. All pixels from seven of the eight images were used to train a random forest classifier with a set of hyperparameters, and the pixels from the eighth image were used to test the accuracy of the classifier. This was repeated eight times, using pixels from each of the training images as the test set, and the average accuracy across all eight ‘folds’ was used to calculate the accuracy of that specific hyperparameter combination. This was then repeated for many sets of hyperparameter combinations to identify the optimal combination that resulted in the most accurate classifier. The optimal hyperparameters found were 
$n\_estimators$=50 and 
$max\_depth$=15. A final random forest classifier was trained on pixels from all eight training images using these hyperparameters (as well as 
$min\_samples\_split$=500, 
$min\_samples\_leaf$=100 and 
$max\_features$=‘sqrt’, which are reasonable values for datasets of this size) and was imported to Google Earth Engine.
3.2. Snow cover identification
The random forest was used to classify all available Sentinel-2 Level-1C images of the study glaciers from 2018 to 2022 between May and November in Google Earth Engine (n = 1.4 million glacier images) at a 10 m resolution. Cloudy pixels were masked using the s2cloudless algorithm (Zupanc, Reference Zupanc2017) with a cloud probability threshold of 20% (pixels with cloud probabilities higher than this were masked). Areas with slopes steeper than 25∘ were also masked to avoid misclassification of these areas due to the high solar incidence angles. Terrain shadows were physically modeled and masked out using the Google Earth Engine TerrainShadow() function, with the Copernicus WorldDEM-30 (GLO-30 DEM) (European Space Agency and Airbus, 2022) and image-specific sun azimuth and incidence angles as inputs.
Snow cover maps of each glacier were made for each observation date from the image classification results, distinguishing between accumulation (snow), ablation (firn, ice, debris and water) and missing data (cloud and shadow) areas. A pixel-wise temporal filter was used to fill missing data and remove outliers in the surface classification by using the most frequently observed class in all other observations up to 7 days before and 3 days after. If there were no usable observations in that period, then the pixels were left as no-data. Frequent misclassification of debris cover as clouds by the s2cloudless algorithm (Zupanc, Reference Zupanc2017) was observed, so the debris cover product from Herreid and Pellicciotti Reference Herreid and Pellicciotti(2020) was used to fill in missing data gaps over these areas. Observations were then discarded if they had less than 30% coverage of the glacier surface before applying the filtering and infilling or less than 70% coverage after applying the filtering and infilling.
3.3. ELA and AAR extraction
A time-varying digital elevation model (DEM) of each glacier was constructed by combining the GLO-30 DEM (European Space Agency and Airbus, 2022) with elevation change products from Hugonnet and others Reference Hugonnet(2021), which provide spatially distributed elevation change rates (
$\mathrm{m\,a}^{-1}$) over glacierized areas globally at a 100 m pixel scale as 5 year temporal averages for four intervals from 2000 to 2020. The original GLO-30 DEM was taken to represent the 2013 glacier surface elevation, as it was constructed using data acquired during the 1 December 2010 to 31 January 2015 TanDEM-X mission (Rizzoli and others, Reference Rizzoli2017). X-band radar penetration into snow and ice may introduce uncertainty (less than the 4 m penetration depth expected for wet snow) (Millan and others, Reference Millan, Dehecq, Trouve, Gourmelen and Berthier2015; Rizzoli and others, Reference Rizzoli2017), but it is considered negligible for our purposes as we are most interested in relative differences in ELA across space and time rather than using the GLO-30 DEM for point measurements of surface elevation changes. DEMs for each year were constructed by applying the 2010–14 and 2015–19 glacier surface elevation change rates, with the 2015–19 product extended for the years 2020–22.
The SLA was extracted for each glacier’s daily snow-cover products using a modification of the approach described in Rastner and others Reference Rastner(2019). The glacier surface area was divided into elevation bands at 10 m increments. The fractional snow cover of the observed pixels within each band was calculated, and bands with fewer than 50 observed pixels were removed. Bands with at least a 50% snow cover fraction were labeled accumulation bands, and those with less than 50% were labeled ablation bands. The lowest point with at least five consecutive accumulation bands was then taken to represent the SLA of the glacier at that time step. If no SLA was identified, the process was repeated for four and then three consecutive bands. If there was still no SLA identified, then the SLA was assumed to be higher than the elevation range of the glacier, and the SLA is assigned the glacier’s maximum elevation. The glacier’s AAR was also calculated at each time step as the ratio between the observed snow-covered area and the total observed area of the glacier.
This resulted in a series of accumulation area maps of each glacier (referred to as ‘observations’ below) with an SLA and AAR associated with each. For each glacier each year, the observation with the maximum SLA was chosen to represent its annual accumulation area and ELA. But first, inaccurate observations were filtered out by identifying outliers in each glacier’s annual SLA and AAR time series. This works on the assumption that the trend in SLA and AAR over the course of a given melt season should be consistent (with SLA increasing and AAR decreasing). Widespread misclassification of the glacier surface in a single image (for instance, from atmospheric noise or clouds not being properly masked out) would result in the image’s derived SLA and AAR being inconsistent with the seasonal progression up to that point in time.
For each observation, a linear regression was fit to the SLAs from all other observations in the preceding 60 days, and the expected SLA for the given date was predicted. If the observation of interest had an SLA which was more than 200 m higher than predicted by this linear regression, then the observation was flagged as an outlier and removed. We tested various thresholds and found 200 m to optimize the trade-off between omitting erroneous outliers due to misclassification (e.g. clouds as snow or clouds as ice) without excluding actual changes in SLA based on visual inspection. The same process was applied using the AARs, and observations with AARs which were more than 0.15 lower than that predicted based on the linear regression were removed. Underestimations of SLA and over-estimations of AAR were never considered outliers, as these conditions are physically realistic effects of new snow accumulation. Similarly, observations from days after the observation of interest were not included in the linear regression because rapid decreases in SLA (and increases in AAR) at the onset of the accumulation season could result in these optimal end-of-summer observations being flagged as outliers.
After these filtering steps, the maximum SLA for each glacier in each year was found, and that observation was taken to represent the end-of-summer accumulation area product and annual ELA. We refer to these products as each glacier’s “annual” products. If multiple observations had the same SLA, then the one with the smallest AAR was chosen. For glaciers larger than 
$500\,\mathrm{km}^2$ (n = 22), additional manual selection was used to ensure that the optimal end-of-summer observation was chosen.
3.4. Final formatting
Each glacier’s annual products were combined to create an average accumulation area product (referred to as the ‘average’ products) by taking the most frequent observation of each pixel (defaulting to accumulation area if the number of years of accumulation and ablation were equal). The ELA was then recalculated for these average products to provide a baseline long-term (5 year average) ELA of each glacier. Remaining missing data gaps were then filled in each annual accumulation area and average accumulation area product. Missing data for elevations at or above the ELA were classified as accumulation areas, and those below the ELA were classified as ablation areas. The majority of these missing data pixels were located at higher elevation, north facing slopes which were in terrain shadows across the entire observation period, and were predominantly infilled as accumulation areas. Not accounting for these missing data gaps would lead to a systematic underestimation of accumulation area across the region. The AAR of each product was then recalculated, and these final infilled products were used for subsequent analyses. The RGI glacier outlines (derived predominantly from imagery over the 2004–10 period) and total glacier area were used for all analyses, as there are no time-varying glacier extent datasets available for the entire study region.
3.5. Climate data
To investigate the climatic controls on glacier ELAs and generate future ELA projections, we utilized ERA5-Land (Muñoz-Sabater, Reference Muñoz2019) reanalysis datasets to investigate the climatic controls on observed glacier ELAs and global climate models from CMIP6 (Eyring and others, Reference Eyring2016) to estimate end-of-century ELA changes. Average daily summer (June, July and August) temperature and total winter (January, February and March) precipitation of each of the 16 O3Regions were calculated from ERA5-Land daily values for 2018–22, and annual anomalies were calculated for each region in each year. A 13-model ensemble (Table S2), selected for its representativeness over North America (Mahony and others, Reference Mahony, Wang, Hamann and Cannon2022), was used for CMIP6 analyses. Ensemble-mean values for the SSP2-4.5 scenario, considered the standard ‘middle of the road’ scenario, were used (O’Neill and others, Reference O’Neill2016). Average summer temperature and total winter precipitation was found for the 2013–22 (modern) and 2090–99 (future) periods to calculate the projected end-of-century changes in each.
4. Results
4.1. Classification accuracy
The random forest classifier showed highly accurate results in the cross-validation testing. Throughout the eight cross-validation folds, the median classification accuracy and F1 score (which combines the precision and recall of the classifier into a single number) across all classes were 93.1% (with ranges of 67.4–99.8% and 66.7–99.8%, respectively). When considering just the performance in identifying snow, the median accuracy increases to 99.5% (ranging from 90.1% to 99.9%) and the median F1 score increases to 98.6% (ranging from 76.9% to 99.8%) (Fig. 3).
Confusion matrices showing the random forest classifier accuracy in classifying pixels within our training/validation dataset. A confusion matrix was created for each of the eight folds of the leave-one-out cross validation approach, and those eight were combined to a single confusion matrix (as shown here) by taking the mean value in each cell. The top shows the confusion matrix for all six surface classes which the classifier was trained to detect, and the bottom shows the same confusion matrix collapsed down to just the snow class and all others. Squares in the top plot are colored by the number of observations in each, and the text inside each true positive square indicates the accuracy of that class (i.e. snow was correctly classified 95% of the time, firn 59%, etc.).
The US Geological Survey (USGS) collects seasonal and annual glaciological mass balance observations at three ‘Benchmark Glaciers’ within our study area (U.S. Geological Survey Benchmark Glacier Program and others, 2016). Mass balance profiles are created from point observations of annual mass balance, and the ELA is calculated as the point at which the annual mass balance is zero (O’Neel and others, Reference O’Neel2019). We used the ELAs from these glaciers (Wolverine Glacier, Gulkana Glacier and Lemon Creek Glacier) to validate the accuracy of our remotely sensed ELA results (Fig. 4).
Comparison between glaciological ELAs and our derived ELAs (Remotely Sensed ELA) on three benchmark glaciers. (a) shows the ELA derived from each plotted against each other, and (b) shows the magnitude of difference between the two (remotely sensed minus in situ). Marker size in (a) indicates the number of days separating the observations, with smaller dots indicating larger time discrepancies (ranging from 5 to 59 days). Note that there are two points for Lemon Creek at 1500 m, indicating 2 years in which the ELA was above the elevation range of the glacier in both the remotely sensed and in situ datasets.
The mean absolute ELA error for all glaciers was 88 m (median absolute error of 44 m), with a mean bias of –4 m (median bias of –13 m). These errors were skewed by two outliers in the dataset on Gulkana and Lemon Creek glaciers which have large time discrepancies. The outlier on Gulkana was caused by significant cloud cover which was misidentified as firn and ice, and the outlier on Lemon Creek was caused by no cloud-free imagery being available in the later summer months. Removing these outliers results in a mean absolute error of 44 m and a bias of –18 m (indicating that the remotely sensed ELAs were on average 18 m lower than the glaciological ELAs).
4.2. Spatial patterns
The automated methods resulted in physically realistic patterns of glacier ELAs and AARs at both the individual glacier and regional scale. A clear relationship between continentality (the glaciers’ distance from the ocean) and ELAs was found, with glaciers with a geometric center nearer to the ocean having lower ELAs than more continental glaciers due to increased moisture availability and precipitation (Fig. 5). For glaciers within 100 km of the ocean, each additional kilometer from the ocean leads to an increase in ELA of 8.6 m (with an intercept of 1027 m at the ocean). However for those glaciers further than 100 km from the ocean, the trend is less clear.
The ELA of each glacier (indicated by the color), as calculated from the 5 year (2018–22) average accumulation area products. The inset scatterplot shows the relationship between distance-from-ocean and ELA. Light colored dots are all glaciers with observable ELAs, and darker purple dots are glaciers with observable ELAs that are larger than 10 km2. The red line indicates the definition of the coastline that was used.
The total accumulation area of all study glaciers is 
$31,334\,\mathrm{km}^2$ (calculated from the 5 year average products), with a corresponding AAR of 0.41 which is well below the AAR range of 0.54–0.64 needed for glaciers of this size to be at equilibrium with climate (Kern and László, Reference Kern and László2010). The AAR of individual subregions varies substantially across the study area, ranging from 0.26 (in the Aleutian Range) to 0.49 (in the Wrangell Mountains) (Fig. 6, Table S3). The glaciers with the highest AARs tend to be large tidewater glaciers (e.g. Hubbard Glacier, Yahtse Glacier Taku Glacier), while smaller glaciers on the periphery of their mountain ranges have the lowest AARs.
The AAR of each glacier (indicated by the color), as calculated from the 5 year average accumulation area products. Inset plots show the total accumulation and ablation area of each O3Region (bottom, left axis), as well as the total AAR (top, right axis), with numbers and colors corresponding to Figure 1.
Our observed AARs correlate well with coincident observations of glacier elevation change (Fig. 7). The average of the 2018 and 2019 AARs of each O3Region was compared to the average 2015–19 elevation change rate from Hugonnet and others Reference Hugonnet(2021). Regions with lower AARs (less accumulation area) exhibited greater elevation change (mass loss) than regions with higher AARs (
${R}^2=0.54$). We also note that regions with marine-terminating glaciers (regions 8–15) display greater surface elevation loss (i.e. greater mass loss) than other regions with similar AARs (Kochtitzky and others, Reference Kochtitzky2022). This is explainable by the increased glacier mass loss from frontal ablation of these glaciers (a process which is not captured by variations in accumulation area) in addition to surface mass balance (the process which variations in accumulation area do capture). Region 13 seems to not fit this pattern as well (Fig. 7), which may be attributable to the Glacier Bay tidewater glaciers undergoing rapid retreat in the early 20th century and now being more or less stable (McNabb and Hock, Reference McNabb and Hock2014).
Comparison between average AAR for 2018 and 2019 and the 2015–19 elevation change rate (Hugonnet and others, Reference Hugonnet2021) in each O3Region. Markers are colored and labeled according to Figure 1, with their size corresponding to the total glacier area in each. Points with black borders indicate regions that contain marine-terminating glaciers, while those with gray borders do not.
4.3. Temporal patterns
Substantial inter-annual variability in ELA was observed (Fig. 8). The median difference between the maximum and minimum ELA of individual glaciers was 470 m (interquartile range of 330–650 m) across the 5 year study period. The highest ELAs across the entire region were observed in 2019, with ELAs being 140 m higher than average across all glaciers, while the lowest ELAs were observed in 2021 at 90 m below the average. The particularly high ELAs in the North Coast Ranges in 2018 and 2019 match with observations on Taku Glacier of ELAs which were 150 m higher than any other year of record (1948–2022), caused by exceptionally warm summer temperatures even at high elevations (Pelto, Reference Pelto2019). The range in observed ELAs resulted in substantial variability in AAR across the study region, with the median difference between the maximum and minimum glacier AAR being 0.46 (interquartile range of 0.30–0.62), with the most pronounced variability found in the North Coast Ranges (Fig. 8, Table S4).
Interannual variability in the ELA of each glacier (a–e), presented as anomalies from each glacier’s average ELA in Figure 5. Red colors indicate higher-than average ELAs, and blue colors indicate lower-than average ELAs. (f) shows the variability in AAR for each glacier across the five years of observations.
5. Discussion
5.1. Methods’ limitations
The automated methods that are presented here allow for annual glacier AARs and ELAs to be efficiently extracted across large regions at high spatial resolution. However, as with any approach that is applied at scale, there are limitations. Although the 5 day return interval for Sentinel-2 provides imagery at a higher frequency than Landsat missions, capturing the true end-of-summer snow cover distribution is still dependent on the availability of suitable cloud-free imagery prior to the date of maximum SLA. In particulary cloudy years, this can result in offsets of more than a month between the true glaciological mass balance minimum and the timing of the best available image (Fig. 4).
The random forest classification approach which we implemented has high accuracy (Fig. 3), but misclassification is still possible particularly in steeper terrain or when cloud cover is not properly identified. Areas with slopes greater than 25 degrees (largely ice falls and glacier headwalls) were masked out in our analysis to limit misclassification of these areas. In rare cases, this may make it impossible to identify the true ELA on glaciers if the ELA falls within the elevation range of an extensive ice fall because too large of an area is masked out. Other machine learning approaches that are able to incorporate the wider spatial context of imagery, such as a convolution neural network, might be capable of improving the glacier surface facies classification. However, such an approach would require a much larger training dataset and be more computationally expensive.
Using the RGI glacier area dataset allows direct comparison between our findings and other studies and facilitates the use of our accumulation area, ELA and AAR products in future studies. RGI outlines in the Alaska region are derived from satellite imagery primarily from 2004–10, and glacier retreat in the intervening years means that non-glacierized areas are included in our analysis. This systematic over-estimation of glacier area causes the AARs to be smaller than if contemporaneous glacier extents were used, but the accumulation area and ELAs are not impacted. Roberts-Pierel and others Reference Roberts-Pierel, Kirchner, Kilbride and Kennedy(2022) analyzed glacier surface area changes over the majority (but not all) of our study region and found that total glacier area in 2020 was ∼94.5% of the 2006 glacier area (the approximate date of imagery from which the RGI glacier outlines were derived). This suggests that the modern-day AAR estimates which we present here may be over-estimated by ∼5.5%, and the true modern region-wide AAR would be 0.43.
Lastly, climate reanalysis products and future climate projections in mountainous regions have inherent limitations due to the scarcity of in situ data and complex topography (Zandler and others, Reference Zandler, Haag and Samimi2019). The end-of-century ELAs that we present here are meant to represent a reasonable estimate of future changes in these regions, rather than absolute projections of glacier-specific AAR and ELA changes.
5.2. Benchmark Glaciers and regional variability
One of the goals of the USGS Benchmark Glacier observations (U.S. Geological Survey Benchmark Glacier Program and others, 2016), and mass balance programs across the globe, is to provide detailed insight into the mass balances processes on the single glacier scale such that those observations can provide standards for comparison against other glaciers in the surrounding regions (Huss, Reference Huss2012; O’Neel and others, Reference O’Neel2019). The interannual variability in mass balance on benchmark glaciers is uniquely captured by seasonal in situ measurements, yet their geometry, location or hypsometry can result in short-term and long-term benchmark glacier mass balance changes that are atypical of the wider regional changes (Fountain and others, Reference Fountain, Hoffman, Granshaw and Riedel2009; Zemp and others, Reference Zemp, Hoelzle and Haeberli2009). The ELA and AAR products in this study provide an opportunity to investigate how well the interannual variability of site-specific benchmark glacier mass balance captures the wider regional mass balance by using the region’s AAR as a proxy for its annual mass balance.
To do so, we compare the annual mass balance of each of the three benchmark glaciers to the annual AAR of their surrounding regions (Fig. 9). The annual balances at Lemon Creek are strongly correlated (P < 0.01) with the AAR of the O2Region and O3Region, with 
${r}^2$ values greater than 0.96 (Fig. 9). The mass balance and regional AAR less well correlate at Gulkana Glacier (
${r}^2=0.73$ and P = 0.07 for the O3Region, 
${r}^2=0.58$ and P = 0.13 for the O2Region), while Wolverine Glacier is the least correlated (
${r}^2=0.33$ and 0.19, 
$P=0.32$ and 0.49).
The higher correlation that we find for Lemon Creek may be due to the larger range of annual balances (2.7 m w.e.) and AAR (0.45) for this glacier and region over the 5 years of observation, as compared to Wolverine (1.0 m w.e. and 0.26) and Gulkana (1.3 m w.e. and 0.27). Alternatively, it may be that there is greater intra-regional variability in mass balance in the regions surrounding Wolverine and Gulkana, or the 5 years may not provide a long enough temporal baseline for such an investigation.
5.3. Climate and ELAs
Winter precipitation and summer temperatures from ERA5-Land were used to investigate the climatic drivers of ELA variations at the regional scale (ELAs and climate variables were grouped by O3Region). Over the 5 year study period, ELAs were found to be positively correlated (P < 0.01) with summer temperatures (higher summer temperature led to higher ELAs) and negatively correlated (P < 0.01) with winter precipitation (increased winter precipitation led to lower ELAs) (Fig. 10). Across the entire Alaska region each 
$1 ^{\circ}\mathrm{C}$ increase in summer temperature was associated with a 62 m increase in ELA, while each 10% increase in winter precipitation was associated with a 24 m decrease in ELA. However, these relationships were not consistent across the entire study area. Coastal regions (O3Regions 8–16) and interior regions (O3Regions 3–7) showed contrasting patterns (Fig. 9). Winter precipitation was found to be a stronger predictor of ELA variability in coastal regions, while the linear relationship between ELA and summer temperature was not statistically significant (P = 0.13). Conversely, for the interior regions, summer temperature was a stronger predictor of ELA, while the linear relationship between ELA and winter precipitation was not statistically significant (P = 0.67).
A multiple linear regression model, using winter precipitation and summer temperature, was fit to all subregions. This model (adjusted 
${R}^2=0.46$) found that each 
$1 ^{\circ}\mathrm{C}$ increase in summer temperature resulted in an ELA rise of 85 m, and each 10% increase in winter precipitation resulted in ELA lowering 32 m.
Relationship between glaciological mass balances of the three Benchmark Glaciers (McNeil and others, Reference McNeil2016) and the AAR of their O2Region (pink dots) and O3Region (purple squares), with each point representing a single year. 
${r}^2$ values are shown for a linear regression of each. Higher 
${r}^2$ values indicate that the benchmark glacier mass balance variations are more representative of the annual AAR variations of their surrounding regions.
The relationship between climate and ELA variability across the study area. Y-axis on each subplots is the annual variability in ELA of each O3Region. X-axis of the left column is the ERA5-Land derived summer temperature (scaled as difference from mean), and the X-axis of the right column is the ERA5-Land derived winter precipitation (scaled by percent difference from mean). Points are colored according to the year which they represent. The top row includes all subregions in the study area. Middle row includes only the coastal subregions (8–16). Bottom row includes only the interior subregions (3–7).
Changes in summer temperature and winter precipitation derived from CMIP-6 data (Eyring and others, Reference Eyring2016) were used to predict the magnitude of end-of-century ELA changes. A single model was used, rather than a region-specific or glacier-specific approach, because the goal was to investigate the realistic magnitude of ELA changes, rather than provide absolute predictions. A more robust approach might include climate projection downscaling or including multiple SSP scenarios, but is beyond the scope of this study.
The CMIP6 ensemble shows that summer temperature increases will range from +1.7 to +3.2 degrees across the study region between the 2013–22 and 2090–99 periods. Winter precipitation is also projected to increase, ranging from +4% to +29% relative to modern precipitation. The resulting end-of-century projected ELA rise ranges from +86 to +249 meters relative to the present day (Fig. 11). The most extreme increases in ELA are projected to occur in the North Coast Ranges, due to the combination of the largest projected temperature increase and smallest increase in winter precipitation. These findings are similar to those of McGrath and others Reference McGrath, Sass, O’Neel, Arendt and Kienholz(2017) which utilized a similar approach but with a much smaller sample size of in situ observations. The regional AAR which we find in this study is substantially lower than the range of possible AARs which they considered when estimating modern glacier ELAs, highlighting the novelty and value which our glacier-specific observations of accumulation areas and ELAs provide.
Changes in summer temperature (left) and winter precipitation (center) between 2013–2022 average and 2090–2099 average from a 13-member ensemble of CMIP6 global climate models under the “middle of the road” scenario SSP2-4.5. The resulting change in ELA, as calculated from the multilinear regression described above, is shown on the right.
The projected ELA increases which we present here are increases in the long-term-average ELA. Individual years are likely to have more extreme changes (both positive and negative). The 2019 mass balance year highlights how anomalies of +300 m in some regions are realistic expectations for inter-annual ELA variability. Anomalies this large, on top of the projected increase in long-term average ELA, could result in individual melt seasons with ELAs 
$ \gt 500\, \mathrm{meters}$ above what we see today, with corresponding significant glacier mass losses in those years.
5.4. Glacier disequilibrium and vulnerability
We modeled the future changes in accumulation area for a range of ELA increases using a hypsometric approach which accounts for glacier area altitude distribution (Fig. 12). For each glacier, the modern ELA was incrementally increased and the resulting accumulation area and AAR were calculated by assuming that all surface area above the ELA would be in the accumulation zone and all surface area below the ELA would be in the ablation zone. The initial AAR conditions were calculated by applying the same hypsometric approach given the modern ELAs of each glacier, rather than using the average observed AAR, to remove the effects of snow cover below and exposed ice above the ELA. The ELA rise was calculated for each glacier by sampling the projected ELA rise shown in Fig. 11 using each glacier’s geometric center point (provided by the RGI ‘CenLon’ and ‘CenLat’).
(a) The end-of-century AAR for each glacier calculated by applying the ELA rise in Figure 11 to modern-day ELAs. (b) The percent of the modern-day accumulation area which would be lost with the projected ELA rise.
Across the entire region, total accumulation area is projected to decrease to 75% of its modern extent by the end of the century, resulting in an AAR of 0.33 given the modern glacier extents (compared to the modern hypsometric-derived AAR of 0.44) (Fig. 12). Substantial differences are found between individual regions, ranging from 12% to 73% loss of the modern accumulation area in the Eastern Chugach Mountains and Central Coast Mountains, respectively (Table S3), which match well with the patterns of end-of-century mass loss projected by Rounce and others Reference Rounce(2023).
A substantial number of glaciers are already entirely within the ablation zone in the present day, which we define as not having an observable ELA or having an AAR of 
$ \lt 10\%$ in the 5 year average products (Fig. 13). A total of 1009 glaciers in our study currently lack an accumulation zone, with a cumulative surface area of over 
$6400\,\mathrm{km}^2$ (8% of the studied area). Without an accumulation zone, these glaciers should be assumed to be entirely out of equilibrium with the modern climate and unable to persist in their present form (e.g.
Pelto, Reference Pelto2010). Under the projected end-of-century scenario of ELA rise, this will increase to 1944 glaciers and 
$16,500\,\mathrm{km}^2$ (22% of the studied area). The majority of these glaciers are small and on the periphery of their mountain ranges (Fig. 13). However, there are numerous large glaciers that are almost entirely ablation zone today or are likely to be by the end of the century, for example, Marvine/Hayden Glacier (
$445\,\mathrm{km}^2$), the eastern tributary of Bering Glacier (
$534\,\mathrm{km}^2$), Tana Glacier (
$515\,\mathrm{km}^2$), Miles Glacier (
$420\,\mathrm{km}^2$) and the glaciers within the Yakutat Icefield (
$ \gt 1000\,\mathrm{km}^2$).
Identification of glaciers with no modern (2018–22) accumulation area. Red indicates glaciers with no modern accumulation area. Orange indicates glaciers that would lose their accumulation area with the projected end-of-century ELA rise.
The end-of-century decreases in AAR will push glaciers even further out of equilibrium with climate than they are in the modern day. The present-day regional AAR of 0.41 is well below the range of values that would be needed for glaciers larger than 
$1\,\mathrm{km}^2$ to be at equilibrium with climate (0.54–0.64) (Kern and László, Reference Kern and László2010). Combining this range of AARs with the modern accumulation area of 
$31,000\,\mathrm{km}^2$ would suggest that the present-day accumulation area is capable of sustaining 49,000–
$58,000\,\mathrm{km}^2$ of total glacierized area, substantially less than the RGI-defined glacier extent of 
$76,000\,\mathrm{km}^2$ of our study glaciers. Under the projected end-of-century ELA rise scenario, the total sustainable glacier area drops to 39,000–
$46,000\,\mathrm{km}^2$, suggesting a potential loss of up to 49% of the current glacier surface area.
These projected glacier area losses will not occur immediately but rather will take decades to be realized, with this response time controlled largely by glacier geometry (Roe and Baker, Reference Roe and Baker2014). Furthermore, the modeled declines in AAR are based only on modern-day surface elevations and glacier area and do not account for future glacier surface elevation changes or surface area changes. In reality, the altitude-mass-balance feedback will increase the glacier surface area below the ELA and amplify the magnitude of glacier area loss, while the reduced glacier area will diminish AAR losses (e.g. Huss and others, Reference Huss, Hock, Bauder and Funk2012; Trüssel and others, Reference Trüssel, Motyka, Truffer and Larsen2013; Sass and others, Reference Sass, Loso, Geck, Thoms and Mcgrath2017).
Our projections of area change are in line with the most drastic scenarios considered by McGrath and others Reference McGrath, Sass, O’Neel, Arendt and Kienholz(2017) over these same glaciers (a 225 m increase in ELA under the RCP 8.5 climate projection), despite using less extreme climate forcings. This is due to the differences in the initial AAR and ELA values used in each study while using the same glacier geometry. McGrath and others Reference McGrath, Sass, O’Neel, Arendt and Kienholz(2017) assumed that glaciers are near equilibrium with modern climate (using AARs ranging from 0.5–0.75) in order to derive glacier-specific ELAs. In contrast, our projections of glacier area change include the committed area loss due to modern-day disequilibrium, highlighting the value of our glacier-specific observations to provide validation and to constrain model initialization in future studies.
6. Conclusions
In this study, we presented a detailed investigation of the modern and future state of glacier accumulation areas and ELAs for over 3000 glaciers which account for 88% of the glacier cover in the Alaska region. The automated methods which we developed using Sentinel-2 imagery could be applied in future years and to other glacierized regions. Furthermore, the derived products have potential uses in applications beyond what is presented here, such as for validating and constraining glacier models, to understand spatial variability in snow accumulation and ablation across a range of spatial and temporal scales, and as a supplement to in situ and geodetic products.
We found that the total AAR of these glaciers is 0.41, with an inter-annual range of 0.25–0.49. The observed AAR is substantially smaller than the steady-state AAR of 0.54–0.64, indicating that, on a regional scale, glaciers are out of equilibrium with modern climate. Furthermore, more than 1000 glaciers were found to have effectively no modern accumulation area. We observed substantial year-to-year variations in ELA and AAR, with an average glacier-specific ELA range of 470 m and AAR range of 0.46 across the 5 years.
We found that the temporal variability in ELAs can be explained by inter-annual climate variability. Summer temperatures and winter precipitation derived from ERA5-Land explained nearly 50% of regional-scale ELA variability, with winter precipitation being a stronger predictor of ELAs in coastal regions, while summer temperatures were more important for continental regions. An analysis of future climate projections indicated that ELAs will rise by an average of ∼170 m (ranging from 85–250 m) under a middle-of-the-road greenhouse gas emission scenario, with the largest changes expected in the North Coast Ranges. We investigated the vulnerability of modern glaciers to these changes and found that the projected ELA rise would result in a loss of 25% of the modern-day glacier accumulation area, while 1938 of the studied glaciers (22% of the investigated area) would be left with no functional accumulation area. In order to reach a new equilibrium with the projected end-of-century climate, up to 49% of the modern glacier surface area may be lost.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2024.65.
Data availability statement
The derived snow cover products and metrics are publicly available at https://doi.org/10.5066/P1QHST6F (Zeller and others, Reference Zeller, McGrath, Sass, Florentine and Downs2024). USGS Benchmark Glacier data are available at https://doi.org/10.5066/F7HD7SRF (U.S. Geological Survey Benchmark Glacier Program and others, 2016). Glacier elevation change products are available at https://doi.org/10.6096/13 (Hugonnet and others, Reference Hugonnet2021). All Sentinel-2 imagery, ERA5-Land reanalysis products, CMIP6 climate projections and the Copernicus WorldDEM-30 are available and were accessed via the Google Earth Engine Data Catalog at https://developers.google.com/earth-engine/datasets/catalog. All codes developed for this study are available through a public GitHub repository: https://github.com/lucas-zeller/Glacier-Snow-Lines.
Acknowledgements
We thank the two anonymous reviewers, Scientific Editor Lauren Vargo and Chief Editor Hester Jiskoot, for their feedback on this work. Figures were created with geographic data from Natural Earth. Free vector and raster map data were obtained from naturalearthdata.com. Color maps in some figures were created by Crameri Reference Crameri(2023). Any use of trade, firm or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Author contributions
Concept and design were done by L.Z., D.M., C.F., L.S. and J.D.; programming and formal analyses by L.Z.; writing by L.Z.; edit and review by L.Z., D.M., L.S., C.F. and J.D.; visualization by L.Z.; and funding acquisition and project administration by C.F., L.S. and D.M.
Funding statement
This work was funded by the U.S. Geological Survey Ecosystems Mission Area Climate R&D Program and the Alaska Climate Adaptation Science Center through USGS award G20AC00095.
Competing interests
The authors declare no competing interests.
Open access
