Glaciers in western North America experience exceptional transient snowline rise over satellite era

Glacier polygons from the Randolph Glacier Inventory (RGI) Version 7 (Pfeffer et al 2014) serve as the basis for this analysis of Alaska (RGI Region 1) and Western Canada and USA (RGI Region 2) (figure 1). This dataset includes 46,239 glaciers, covering 101,229 km². Over 95% of the outlines originate from satellite imagery captured between 2004 and 2010, with the latest outlines from 2011. Although updated glacier outlines exist for British Columbia and Alberta (e.g. Bevington and Menounos 2022), they do not encompass our entire study area. To reduce edge effects and shadows, glaciers smaller than 2.0 km² are excluded from the study, along with 14 poorly digitized high elevation snowfields that lack ablation zones. These exclusions result in 4456 glaciers remaining, representing 84.5% (85,579 km²) of the total glaciated area.

Optical satellite imagery from 1984 to 2024 is used in this analysis, including five calibrated surface reflectance collections from the Landsat and Sentinel-2 programs with 20–30 m spatial resolution (table 1). Images acquired from July to September with <20% cloud cover in their metadata are kept. Each year, glaciers are, on average, imaged twice before the year 1999 and 12 times after 2019 (supplementary figure 1). Imagery is processed using the 'rgee' R package, which interfaces with Google Earth Engine (GEE) (Gorelick et al 2017, Aybar et al 2020). Cloud masks included in these collections perform poorly over glacier margins, and thus we include clouds in our image classification (Bevington and Menounos 2022).

Satellite imagery is classified using random forest (RF) models in GEE; these models are trained separately for each image collection (Gorelick et al 2017, Lu et al 2020, Rittger et al 2021). Training points are manually labeled into six classes: snow, glacier ice, cloud, water, bare ground, and vegetation. Surface reflectance bands (green, red, near infrared, and shortwave infrared) and seven spectral indices (e.g. NDVI, NDSI, see supplementary table 1) are included as predictors in the models. Training data are derived from three seasonal image mosaics per collection to capture spatial and temporal variability in surface reflectance (table 1). Each RF model was built using 900 training points (50 per class per image) randomly split into 70% training and 30% testing datasets. RF model parameters included 500 trees, five variables per split, a minimum leaf size of two, and a bag fraction of 0.5.

The TSL is calculated per image as the 5th elevation percentile of the snow within the glacier polygon which avoids biases from low-elevation snow patches (Shea et al 2013). The M_ATSL is calculated as the maximum TSL elevation observed per glacier in a given year. We use the 30 m Copernicus GLO-30 digital elevation model (DEM), and resample the imagery to align with the DEM grid (Rizzoli et al 2017). Images are excluded where <90% of a glacier polygon overlaps with valid imagery or >10% of pixels are classified as clouds. In some cases, partial snow accumulation areas may be used to calculate the TSL, but the majority of images cover the entire glacier (supplementary figure 2). In cases where two late summer partial images of the TSL exist, whichever has the higher TSL will be selected. Glaciers with fewer than one usable image per two-year interval are removed, excluding 257 glaciers, primarily in the Aleutian Islands (sub-region 01–01) and at the Yukon-Alaska border, leaving 4199 glaciers in the analysis (supplementary figure 3). When no snow is detected, the TSL is set to the glacier's maximum elevation, indicating snowline absence. Linear regression analysis per glacier allows us to quantify the rate of M_ATSL change over time, providing corresponding coefficients, standard errors and r² values (table 1). Area weighted sub-regional and regional statistics are also presented which scale the TSL values of each glacier by the glacier area.

In this study, we investigate climatic drivers of the M_ATSL since many note the influence of synoptic climatology on regional patterns of glacier mass balance (Moore et al 2009, Menounos et al 2019, Rounce et al 2023). In addition, many observe significant covariance between atmosphere-ocean teleconnections and the duration of both seasonal snow cover and the magnitude of streamflow from glacierized catchments (Bitz and Battisti 1999, Fleming et al 2006, Moore et al 2009, Bevington et al 2019, Gleason et al 2020). Climate drivers of M_ATSL are analyzed using ERA5-Land reanalysis data (Muñoz-Sabater et al 2021). Mean summer air temperature (MSAT; July–August) and mean spring snow water equivalent (MSSWE; May) are extracted from 9 km resolution datasets as per-glacier area-weighted averages and accessed via GEE. Atmosphere-ocean teleconnection indices, notably the Arctic Oscillation (AO), Oceanic Niño Index (ONI), Multivariate ENSO Index v2 (MEI), and Pacific Decadal Oscillation (PDO), are included and accessed using the 'rsoi' R package (Albers 2023).

The RGI polygons include terrain attributes computed using the Copernicus GLO-30 DEM (Pfeffer et al 2014). We process these attributes into i) average slope of the glacier (in degrees); ii) terrain aspect, which is separated into the unitless northness (sine of aspect in radians) and eastness (cosine); and iii) the hypsometric integral (HI), which is calculated from the mean, minimum and maximum elevation of the glacier (HI = mean–min/max–min). Terrain attributes are used to investigate regional patterns in the M_ATSL trend and in a multiple linear regression model to understand the influence of trends in climate on trends in M_ATSL while considering terrain at the level of the individual glacier.

A total of 1331, 1164, 693, 186 and 522 images were classified from the L5 TM, L7 ETM+, L8 OLI, L9 OLI, and S2 MSI, respectively. The accuracy of the five individual RF models exceeds 95% (table 1), with the greatest importance attributed to the near-infrared and red bands in all models (supplementary figure 4). The chosen classes were sufficiently unique to be mapped with confidence, and the surface reflectance values for each class in the training data were similar among sensors (supplementary figure 5). The snow class had high reflectance in the visible and near infrared bands and low in the shortwave bands, while glacier ice was similar with an attenuated visible signal (supplementary figure 5). Additional visual inspection confirmed that the models classified snow properly (e.g. figure 2). We analyzed the confusion matrix provided by each RF model and found that the water and glacier ice classes were most commonly confused and responsible for the 95%–99% accuracies (Stehman 1997).

Zoom In Zoom Out Reset image size

We provide two examples of individual glacier M_ATSL (figure 2) for Brintnell Glacier (RGI2000-v7.0-G-02-00562) in sub-region 02–01 (Mackenzie and Selwyn Mountains) and Bridge Glacier (RGI2000-v7.0-G-02-07266) in sub-region 02–02 (S Coast Mountains). The Brintnell Glacier M_ATSL has a positive trend of 4.5 ± 1.4 m a⁻¹ (r² = 0.2, p < 0.01). The maximum M_ATSLs occurred in 2023 and 2024, and were 2347 and 2346 m asl, respectively (approaching the maximum elevation of the glacier 2467 m asl). The snow covered area was less than 0.23 km² in both years (∼1% of the total glacier area of 22.1 km²), indicating an absent accumulation area. Similarly, the Bridge Glacier time series has a positive trend of 5.8 ± 1.5 m a⁻¹ (r² = 0.3, p< 0.01). Maximum M_ATSLs of 2281 and 2252 m asl were observed in 2023 and 2024, respectively. However, Bridge has a much higher maximum glacier elevation (2908 m asl) than Britnell Glacier, and as a result there was still a large snow covered area of 42.2 and 43.7 km² in 2023 and 2024, representing 51.2 and 53.1% of the total glacier area of 82.3 km².

The mean M_ATSL is lowest for coastal Alaskan regions (808 m asl in sub-region 01–03 and 1041 m asl in sub-region 01–04), and highest for South Rocky Mountains (sub-region 02–05, 3639 masl) (figure 1, table 2). The average M_ATSL over the entire study area rose by 2.7 ± 0.4 m a⁻¹ (r² = 0.3, p< 0.01) between 1984 and 2024. This rate increased fourfold from 2.1 ± 0.8 m a⁻¹ in 1984–2010 (r² = 0.1, p < 0.01) to 8.9 ± 1.7 m a⁻¹ in 2010–2024 (r² = 0.5, p < 0.01). M_ATSL in North Alaska (sub-region 01–01) had the greatest rate of M_ATSL rise in the 1984–2024 period, 5.6 ± 1.8 m a⁻¹ (r² = 0.2, p < 0.01), and the greatest rate from 2010–2024 was in the Cascade and Sierra Nevada (sub-region 02–04) with a rate of 15.5 ± 3.2 m a⁻¹ (r² = 0.6, p < 0.01). Following the year 2010, the greatest increases in M_ATSL rates of change were by factors of 6.4, 9.2 and 9.5 for sub-regions 02–04, 02–02 and 02–03, respectively (table 2).

The highest M_ATSL values were in 2019 and 2023, which occurred in RGI Regions 1 and 2, respectively (figure 1). In Region 1, the 2019 median M_ATSL was 155 m above the long-term average, whereas in Region 2, the 2023 median M_ATSL anomaly was 148 m. In addition, the number of glaciers in each region with M_ATSL values greater than the 95th elevation percentile of the glacier (Z₉₅), a proxy for the complete loss of accumulation area, increased in Region 1 from 0.95 ± 0.46 (p > 0.05) to 1.74 ± 1.27 (p > 0.05) glaciers per year in 1984–2010 and 2010–2024 epochs, respectively. In Region 2, these rates changed from 0.35 ± 0.13 (p < 0.05) to 5.52 ± 2.37 (p < 0.05) glaciers per year. The years 2019 and 2023 saw 91 and 149 glaciers with M_ATSL > Z₉₅ in Regions 1 and 2, respectively.

We compare our results to studies available in our study area. Shea et al (2013) report a manually derived 2001 snowline for the Lillooet Icefield of 2102 m asl, and a 2009 snowline for the Wapta/Waputik Icefield of 2631 m asl. Our workflow detected the Lillooet and Wapta/Waputik snowlines to be 2078 and 2624, respectively. Shea et al (2013) used the 10th percentile of the snow covered area elevation to define snowline, whereas we use the 5th percentile, which explains why our numbers are lower. We also compare our findings to those reported by Pelto (2019) for Taku Glacier. In that study an exceptionally high M_ATSL value of 1425 m asl is reported for 10 October 2018, whereas we find a M_ATSL of 964 m asl on 16 September 2018, this is likely due to our study excluding the images from the month of October. For other years reported by Pelto (2019) our results agree with an average difference of ±46 m. These differences could be due to the difference in workflows, different DEM sources, interpretation of snow and firn, or due to the different date selected.

As expected, M_ATSL is inversely related to MSSWE in both Region 1 (−351 mm w.e.⁻¹, r² = 0.48) and Region 2 (−155 mm w.e.⁻¹, r² = 0.37), and positively correlated with the MSAT in both Region 1 (46 m °C⁻¹, r² = 0.54) and Region 2 (23 m °C⁻¹, r² = 0.30) (figure 3). In Region 1, 2019 had the 2nd highest MSAT (8.9°C) combined with the 2nd lowest MSSWE (5.03 m w.e.) in the same year, likely causing the very high M_ATSL values. Similarly in Region 2, both 2023 and 2024 had high MSAT (12.0 and 12.3°C) combined with the low MSSWE (0.97 and 1.04 m w.e.) in the same year.

Zoom In Zoom Out Reset image size

Strong positive and negative correlations are observed between M_ATSL and MSAT and between M_ATSL and MSSWE, respectively (figure 4). There are glaciers in Alaska and Yukon where the relation between MSSWE and M_ATSL is not significant, and there are only a few glaciers where the relationship is positive and significant (see sub-region 01–02 and 01–01). More regional variability exists when M_ATSL and atmosphere-ocean teleconnections are considered. M_ATSL is positively and negatively correlated with the ONI and PDO, although in many regions these were not significant. Resulting in higher M_ATSL values in El Niño years (positive ONI) and in cool-phase PDO years, with regional variability. There is a weak correlation between MEI and AO with M_ATSL and no apparent spatial pattern between these indices and M_ATSL exists (figure 4).

Zoom In Zoom Out Reset image size

In both regions terrain influences M_ATSL trends similarly, with the most pronounced M_ATSL trends occurring on south and east-facing aspects, steep slopes, on glaciers with lower HIs, and at lower elevations (supplementary figure 6). In a series of multiple linear regression models that predict M_ATSL per year, models that combine terrain with climate outperform models that use only climate (supplementary figure 7). In Region 1, the model r² values increase by a factor of 1.5 when terrain is considered, whereas in Region 2, the model is improved by a factor of 4. These annual models also offer insights into the importance of model covariates. By examining each year individually, we observed that aspect is becoming a less significant covariate over time, whereas glacier hypsometry is becoming increasingly important (supplementary figure 8).

Our workflow provides a means to automate per-glacier TSL mapping at the regional-to-global scale over the modern optical satellite era (last 40 years). Our method presents a significant improvement in speed and scalability over studies that use manual or threshold-based approaches to map TSL (Chan et al 2009, Pelto 2019, Racoviteanu et al 2019, Rastner et al 2019). The RF classification accuracy exceeds 95% (table 1) and multiple visual inspections and independent data validation corroborate the performance of the classifier (e.g. figure 2). Our workflow can be used for rapid and continual annual updates of TSL in the study area and using GEE this only takes a few hours to run. Expanding this work to other regions will require additional region-specific training data.

Our findings reveal a pronounced rise in M_ATSL in western North America, with an increased rise after 2010 (figure 1, table 2). Notable M_ATSL maxima occurred in 2019 (Alaska) and 2023 (Western Canada and USA), both of which reflect simultaneous low MSSWE and high MSAT for both regions (figure 3). The year 2024 was also an above normal M_ATSL year for Western Canada and USA. These findings align with the broader consensus of accelerated glacial retreat under the influence of global warming, corroborating the increased mass and area losses, and post-2010 acceleration reported in previous studies (Hugonnet et al 2021, Bevington and Menounos 2022, Roberts-Pierel et al 2022).

The ratio of the accumulation area to ablation area is likely one of the primary factors that helps explain the rapid acceleration of glacier mass loss (Criscitiello et al 2010, Rabatel et al 2012, Rastner et al 2019). The exceptional rise in M_ATSL not only signifies a reduction in the accumulation areas of glaciers but also a darkening of glaciers due to the darker optical reflective properties of firn and glacier ice with prolonged melt seasons (Elias et al 2021). Projections of deglaciation may currently underestimate the speed at which the early 21st century snowline is changing and, consequently, the rapid loss of the accumulation area (Clarke et al 2015, Rounce et al 2023).

Large area historic M_ATSL mapping provides improved observational datasets to calibrate and validate hydrologic models. These models often partition glacier runoff into flow from snow and ice surfaces with little data to validate these assumptions (Brandt et al 2022, Robinson et al 2024). In addition, near real-time snowline mapping could also improve short and medium-range hydrological flood forecasting (Tennant et al 2015, Fehlmann et al 2019).

This analysis of western North America reveals regional patterns in average M_ATSL (figure 1). For example, Alaska has lower elevation snowlines, higher MSSWE, and warmer MSAT than Western Canada and USA. Consequently, Region 1's M_ATSL sensitivity to changes in both MSAT and MSSWE is twice that of Region 2 (figure 3). This could be because glaciers are larger and less sheltered with flatter hypsometry. The regional differences in the M_ATSL trends also point to the complex interplay between local climatic conditions and global climate drivers. For instance, the variability in M_ATSL responses between different glacier regions may be influenced by local topographical features, continentality, and microclimates, which modulate local climate-glacier interactions (Fleming et al 2006). The M_ATSL sensitivity to climate is likely influenced by the elevation, slope, aspect and hypsometry of a glacier (Kienholz et al 2015, Davies et al 2024). Understanding the sensitivity of glaciers to changes in M_ATSL requires integrating terrain statistics with climate trends, as highlighted by the multiple linear regression model presented. These models emphasized the critical role of MSAT and MSSWE, and their interactions as a dominant control on M_ATSL. The results align with North American glaciological studies that link increased summer temperatures and reduced spring SWE to rising M_ATSLs, underscoring the influence of warming and snowpack decline on glacier mass balance (Huss and Hock 2018).

Terrain factors also play a crucial role, as indicated by the significant negative slope coefficient and the negative HI coefficient, suggesting that flatter glaciers are more sensitive to M_ATSL changes (supplementary figure 5). The positive slope coefficient indicates that the steepness of a glacier enhances its response to M_ATSL changes, likely through dynamic effects, while the negative hypsometric coefficient highlights the vulnerability of flatter glaciers to climatic shifts because of their limited accumulation areas and less protective elevation gradients. Together, these findings reflect the interplay of dynamic and climatic controls on glacier sensitivity. Steep glaciers may respond dynamically to changes, but flatter glaciers are more prone to losing mass relative to their total area when the M_ATSL rises due to their geometry and elevation distribution. This dual influence aligns with observations of glacier behavior across diverse regions in North America, where both steep, dynamically active glaciers and flatter, low-elevation glaciers display varying degrees of sensitivity to climate forcing (Clarke et al 2015).

Several factors continue to challenge the accuracy and reliability of automated TSL mapping. For example, the small number of historical images limits our ability to detect the M_ATSL accurately every year (supplementary figure 1). It is likely that we occasionally underestimate M_ATSL when less images are present, particularly for older epochs with fewer images. We conducted a sensitivity analysis on the number of images included per year and found that the overall M_ATSL trend stabilizes when we include eight images per glacier per year or more (supplementary figure 9). This question should be looked at in more depth in future work. We also tested the impact of the cloud cover and valid pixel thresholds used in our image filtering and found that small changes in the thresholds do not substantially change the results (supplementary figure 10). Our analysis often has less than eight scenes per year, which suggests that we may not be capturing the M_ATSL every year. To test this, we estimate regional M_ATSL using multiple linear regression analysis with MAST and MSSWE (supplementary figure 11). Our M_ATSL observations align with the model; however, in 2023 and 2024, Region 2 shows higher-than-expected M_ATSL. The under prediction of regression model could signal enhanced melt due to darkening of snow and firn surfaces (e.g. impurity loading or increased snow grain wetness/growth). Combining workflows with higher frequency datasets (e.g. MODIS, Planet, etc.) would increase the resolution of detecting the M_ATSL, but these datasets do not cover the entire period of study. We also investigate the day-of-year of the detected M_ATSL over time and do not observe any trend in the timing of the M_ATSL (supplementary figure 12).

Another limitation is that machine learning workflows require large amounts of training data. In this study, we manually digitized 900 training points per image collection and used five image collections. However, more training data that covers more diverse snow facies (e.g. dirty ice, snow algae, etc.) may lead to an even better model performance. In some cases, landslides occurred over glacier surfaces, volcanoes erupted onto glaciers in the Aleutian Islands, and debris emerged from within the glacier over the time series. These site-specific challenges require more attention in future work and may be improved using sematic segmentation or computer vision algorithms (Shankar et al 2024).

Additionally, we based our analysis on the publicly available RGI Version 7 glacier polygons. Therefore, our workflow inherits any errors in that dataset. We manually removed glaciers smaller than 2 km² to avoid edge effects and complex shadows, we also removed glaciers that were digitized entirely in high elevation accumulation areas as those polygons have no ablation area and the TSL cannot be detected (discussed earlier). As such, we only included 4199 glaciers. We assessed the sensitivity of the workflow to the sample size and found that overall M_ATSL trends over time remained within 15% with as few as 900 randomly selected glaciers in the study area (supplementary figure 9).

Lastly, we computed M_ATSL using a static elevation model over time that does not reflect the surface elevation change of the glacier over the period of study. The lowering of glacier surfaces over time (e.g. Hugonnet et al 2021) may have caused an overestimation of our findings because we did not use a time-dynamic DEM. Annually resolved time-dynamic DEMs would improve future work of mapping M_ATSL over time. However, using a static DEM could offer a clearer signal of climate change impacts (Østrem 1973, Østrem and Haakensen 1999).