Accelerating glacier volume loss on Juneau Icefield driven by hypsometry and melt-accelerating feedbacks

Globally, glaciers and icefields contribute significantly to sea level rise. Here we show that ice loss from Juneau Icefield, a plateau icefield in Alaska, accelerated after 2005 AD. Rates of area shrinkage were 5 times faster from 2015–2019 than from 1979–1990. Glacier volume loss remained fairly consistent (0.65–1.01 km3 a−1) from 1770–1979 AD, rising to 3.08–3.72 km3 a−1 from 1979–2010, and then doubling after 2010 AD, reaching 5.91 ± 0.80 km3 a−1 (2010–2020). Thinning has become pervasive across the icefield plateau since 2005, accompanied by glacier recession and fragmentation. Rising equilibrium line altitudes and increasing ablation across the plateau has driven a series of hypsometrically controlled melt-accelerating feedbacks and resulted in the observed acceleration in mass loss. As glacier thinning on the plateau continues, a mass balance-elevation feedback is likely to inhibit future glacier regrowth, potentially pushing glaciers beyond a dynamic tipping point.


JUNEAU ICEFIELD GLACIOLOGY AND CLIMATE
3 THE "LITTLE ICE AGE" 3.1 Climate during the "Little Ice Age" In the Gulf of Alaska, the "Little Ice Age" (LIA) (1250 -1850 AD) is well recorded in a number of palaeoclimatic proxies, including dendrochronology 6,7 , coleopteran 8 , and ice cores 9 .In central Alaska, at the continental-climate Farewell Lake in central Alaska, the coldest part of the late Holocene culminated at 1700 AD 10 , with temperatures 1.7°C cooler than present.At Iceberg Lake, Bagley Icefield, the timing of maximum cooling is at 1650 AD and 1850 AD, when temperatures were 1°C cooler than modern 11 .In southwest Alaska, a coleopteran palaeoclimate record gives a summer temperature of 1.3°C cooler than the modern mean (1983-2016), with cooling centred on AD 1815 AD 8 .Ice core data from Mount Hunter indicate that temperatures are now 2°C warmer than the height of the LIA 9 .
3.2 Evidence for glaciation in Alaska during the 'Little Ice Age' Across the American Cordillera, there is widespread evidence of a Late Holocene Neoglaciation, often termed the "Little Ice Age" (LIA) [12][13][14] .This last Neoglacial advance left behind a distinctive geomorphic system of sharp-crested moraines, trimlines, ice-scoured bedrock and thicknesses of glacial sediment.These advances were frequently the most extensive of the last 10,000 years 14,15 .Most glaciers in Alaska below 1500 m asl have an uninterrupted history of continuous recession since this maximum 16 .
There is long history of research and substantial evidence of a readvance of glaciers around Juneau Icefield during the LIA [17][18][19][20][21][22] .A holistic geomorphic map was recently published 23 , but a reconstruction of icefield extent at this time has not yet been attempted.Around the Juneau Icefield, outlet glaciers advanced in the prior to the LIA, in the last millennium, but these advances were generally within "Little Ice Age" limits 14,24 , which dated to ca. 1750-1770 in this area 21,25 .The timing of glacier recession was largely synchronous across the icefield, with recession beginning for most glaciers between AD 1750 and 1785 17,25 .Local relative sea level data indicate that larger-than-present glaciers had stabilised by the mid-16 th century, with land first emerging as glaciers shrank between AD 1770 -1790 25 .Glaciers have been consistently shrinking since this maximum 18 , with the exception of Taku Glacier, the largest outlet glacier of Juneau Icefield, which has been shrinking since 2018 26 .
Taku Glacier (Supplementary figures 1C, 5) reached its maximum position three miles beyond the 1948 terminus, near Taku Point, in AD 1750 27 to AD 1755 17 , where it was confluent with Norris Glacier 17,28,29 .These dendrochronological ages agree with radiocarbon ages from nearby Loon Bog 27 .At this time, Taku Glacier blocked the end of Taku Inlet 18 , calving into the fjord, and dammed Taku River, causing an icedammed lake to form behind the glacier snout 17,27  In this work, we assume that the LIA maximum at Juneau Icefield occurred at AD 1770.However, due to the uncertainty in the timing of the LIA, and due to the presence of records in other locations indicating a later maximum 6,8,13,15 , we also calculate rates of recession at from a maximum at 1880 AD and for an earlier LIA maximum at AD 1675 15,32 .
4 SUPPLEMENTARY METHODS

Data sources, 1948 and 1979
Supplementary figures 6 and 7 below shows the source of the imagery used in the reconstruction of glacier area in 1948 and 1979.Aerial photographs were used wherever they were available.Where they were not available, such as over Canada, topographic maps for the main outlet glaciers from this era were used (Supplementary table 5).Some of the smaller glaciers peripheral to the icefield, especially in the northeast, were not viewable in the map or aerial photography, and so the earliest available imagery (in this case, the 1980s satellite imagery) was used.The date of the source imagery was recorded in all cases in attribute information and used when calculating rates of recession.The 1948 icefield reconstruction should therefore be considered a minimum, as a few glaciers were mapped using later imagery (Supplementary table 3), which is recorded in the attribute information.The date is taken into account when calculating annualised rates of change. Aerial

Assessment of Topographic Maps
Reconstructing glacier extent used USGS topographic maps where aerial photographs were not available, or where there were gaps in the orthomosaics.Data was taken directly from information on the maps.Topographic maps were downloaded from the USGS map viewer (https://ngmdb.usgs.gov/topoview/viewer/#4/52.19/-123.71)as georeferenced GEOTIFFS.Our qualitative quality assessment (Supplementary table 5) used comments directly from the maps and compared topographic maps to the Arctic DEM, ASTER GDEM and ESRI Basemap (usually GeoEye).The location of topography (valleys, hills, lakes, rivers, etc) was compared in the datasets.Maps with a high degree of correspondence to the high resolution remotely sensed datasets were rated as 'high confidence'; highly stylised glaciers and a low degree of correspondence between lakes, rivers, glaciers and topography were rated as 'low confidence'.Maps where the aerial photographs were unclear were also rated as 'low confidence'.The date of the glacier extent was taken to be the latest date of the aerial photographs used.Unless otherwise specified, 1948 aerial photographs were presumed to have been taken on 13.08.1948.5 SUPPLEMENTARY RESULTS

~1770
Short periods of faster rates of volume loss or readvances cannot be excluded due to a paucity of data prior to 1948, and Taku Glacier in particular is known to have receded behind the 1948 extent by 1890 AD (Supplementary figure 11, Supplementary table 12), followed by a readvance 17 .However, this behaviour is likely to be anomalous, and related to Taku Glacier's evolving terminal environment 26,33 .Temperature records from ice cores and tree rings 6,9 suggest uninterrupted, steady warming between the LIA and 1948, which would likely result in steady glacier recession, though temperatures reconstructed from varved lake sediments north of the study area suggest do some cooling from 1870-1900 (Figure 6a) 11 .Further, historical records from 1890 and 1908-1910 AD from Twin Glacier Lake (Supplementary table 12), Mendenhall (Figure 2), Eagle and Herbert Glacier (Figure 5) all indicate continued recession from the LIA to 1948 AD 17 .There are few moraines between the mapped LIA neoglacial maximum and the 1948 extent (Supplementary figure 9), which would suggest a steady recession, uninterrupted by readvances.
Icefield-wide volume loss from 1949 to 1979 AD was -1.06 ± 0.74 km 3 a -1 (Figure 6d).Thinning is concentrated in the glacier tongues, below the ELA, with elevation change on the icefield plateau close to zero (Figure 7).On Field Glacier, glacier thinning reaches an elevation of 1200 m asl, and 1246 m on Ogive Glacier.Thinning greater than 5 m is observed on Gilkey Glacier up to 1520 m asl, Meade Glacier to 1230 m asl, on Eagle Glacier to 1230 m asl, and on Herbert Glacier to 980 m asl.Rates of thinning (mean dh/dt across the glacier) were low (mean -0.27 m a -1 , Figure 6e), with the highest values from Ogive Glacier (areaaveraged elevation change of -1.58 m a -1 ) and nearby valley glaciers.Other rapidly thinning glaciers include Thiel (-1.09 ma -1 ) and Denver (-0.85 m a -1 ).Glacier area-averaged elevation change on Taku Glacier was very slightly positive in this time period.Taku Glacier terminus thickened and advanced 1,486 m between 1948 and 1979 as it began to fill in the forefield over-deepening with sediment 33 , with a total area growth of 8.22 km 2 .Shoaling sediment is visible in front of the terminus in the 1979 aerial photographs.
Icefield-wide volume loss from 1979-2000 reached -3.72 ± 1.57 km 3 a -1 , indicating an acceleration of thinning compared with the 1948-1979 time period.Median rates of area-averaged elevation change (dh/dt) for outlet glaciers reached -0.48 m a -1 , with valley glaciers at -0.72 m a -1 (Figure 6f, g).Thiel and Ogive were again the fastest-thinning glaciers in 1979-2000 (glacier-wide mean thinning of -2.43 and -2.03 m a -1 respectively).The terminus of Taku Glacier thickened by 20 m over this time period (1979-2000), though elevation change on the plateau was not observable.Taku Glacier advanced between 1979 and 1990 (a distance of just 140 m) (Supplementary tables 9, 12).Terminus thickening of up to 70 m allowed the Hole-In-The-Wall outlet to flow over a low col by the 1940s, forming a large piedmont glacier on the lowlands by 1960.Thinning was observed on the remainder of outlet glacier tongues across the rest of the icefield; the terminus of Meade and Field glaciers each thinned by up to 100 m.Gilkey Glacier thinned by up to 80 m, and similar values were observed on Mendenhall (Figure 7).
As a proportion of their area, the smallest glaciers shrank fastest from 2015-2019 (Figure 6d; Supplementary figure 12).Mountain glaciers shrank at total summed rate of 3.73 % a -1 from 2015-2019, and glacierets shrank at 7.76 % a -1 , whilst outlet glaciers shrank at just 0.40 % a -1 (Figure 6b; Supplementary figure 12).Supplementary figure 12 shows how, from 2015-2019 AD, glaciers peripheral to the main icefield shrank fastest relative to their area, whilst the main outlet glaciers had the lowest normalised rates of recession.A comparison of glaciers east and west of the ice divide and lacustrine versus land-terminating glaciers found no statistically significant difference in recession.Glaciers that have a significant proportion of debris (>10% debris cover by area) are receding more slowly; debris-rich glaciers are receding at a mean rate of 5.51 % a -1 compared with 3.91 % a -1 for debris-free glaciers (t(982) = 2.4, p = 0.00).However, glaciers with >10% debris cover are thinning faster (mean of -1.19 m a -1 versus 0.76 m a -1 for debris-free glaciers) (t(982) = 5.76, p = 0.00).This supports downwasting rather than areal recession for those glaciers with the most debris.
Icefield-wide volume loss from 2000-2010 AD was 30.79 ± 10.13 km 3 , and doubled to 59.09 ± 7.99 km 3 from 2010-2010 AD.This means that the icefield lost 5% of its volume in the period from 2010-2020 alone.Rates of volume loss are rapidly accelerating (Figure 6d, e, Table 2), with outlet glaciers exhibiting the fastest rates of thinning (dh/dt in Figure 6e), closely followed by valley glaciers.Thinning occurred across the plateau (Figure 7).Thinning at several outlet glacier termini is over 5 m a -1 (including Eagle, Gilkey, Field glaciers).Decreased thinning right at the snout of some glaciers, such as Tulsequah, indicates that these outlet glaciers are likely floating at their terminus, with thinning from below rather than the ice surface (Figure 7).11).
Mendenhall Glacier has receded across a large lake across the last 70 years, with substantial losses from the terminus as a result of increased calving due to terminus flotation 36,37 .The terminus is now partly grounded above lake-level, with a greatly reduced calving, no flotation of the terminus, and reduced frontal ablation potential from a much narrower lacustrine terminus.This has resulted in a slight slowing of the rate of terminus recession (Supplementary table 12).The rate of terminus recession at Herbert Glacier has also slowed in recent years, likely also due to recession into the narrow valley and away from the proglacial lake.
In contrast, Gilkey Glacier receded fastest from 2005-2015 (0.60 % a -1 ), with rates of recession slowing from 2015-2019.Field Glacier receded at 0.67 % a -1 from 2005-2015 and then recession also slowed, at 0.14 % a -1 from 2015-2019.Significant thinning (7 m a -1 ) is now occurring 1.1 km up from the glacier terminus (Figure 7), though the area at the terminus, where rifting and calving occurs, is not thinning, due to ice flotation in the proglacial lake.Gilkey and Field glaciers both receded back across an over-deepened forefield filled with a proglacial lake.In the 2019 imagery, abundant calving is visible in these lakes, again indicating substantial frontal ablation (cf. 23).A slowdown of recession for these glaciers may therefore indicate stretching and flotation of the ice as it thins over an over-deepened basin.In order to investigate whether differences in sensor contributed to this decrease in albedo, we calculated change in albedo over time using just Landsat 7 (LE07) data from all seasons (Supplementary table 18) and Landsat 5 and 7 data ( Supplementary figure 15c, d).For just Landsat 7 alone, the icefield wide albedo within the 1990 glacier outlines was also lower from 2010-2023 (mean = 0.81, standard deviation = 0.1) than from 1999-2009 (mean = 0.87, standard deviation = 0.02); t (131) = 2.76, p < 0.05.

Glacier snowlines
This result was replicated when including the Landsat 5 sensor as well, with lower albedos from 2010-2023 (mean = 0.81, standard deviation = 0.11) than from 1987-2009 (mean = 0.86, standard deviation = 0.11); t (211) = 3.28, p < 0.05.When including the longer time series provided by both Landsat 5 and 7, the difference was greatest in the summer season, though albedos decreased in all seasons (Supplementary figure 15d).Decreases in albedo in autumn and winter may reflect an observed shortening of the ablation season and lengthening of the accumulation season 41 .

Supplementary figure 1 .
Field photographs of Juneau Icefield, all from July 2022.A, B: the low-slope plateau accumulation area of Taku Glacier.C: The piedmont terminus of Taku Glacier, showing the shoaling moraines building up.D: Looking down Gilkey Glacier, with ogives on the tongue, and the icefall in the foreground.E: Fragmentation of a small glacier on the icefield.F: Avalanche Canyon, a deep, ice-scoured valley, with Gilkey Glacier in the background.G: recently exposed ice-scored bedrock at high elevations above Avalanche Canyon.H: Neoglacial terminal moraines in the Herbert Glacier forefield.All photographs credit BJD.Supplementary figure 4. ERA5 climate reanalysis data 5 for NW North America.Location of Juneau is shown with yellow star in panel a.A-d show the grid cell difference from 1990-2005 mean versus the 1950-1980 mean.E-h show the difference for the 2016-2020 mean versus the 1950-1980 mean. A. t2m air temperature at 2 metres above ground level.B. t850 temperature at 850 mb (approximately 1500 m asl, close to the Juneau Icefield plateau).C. vmid vertically integrated moisture difference.D. tp total precipitation.E. t2m air temperature at 2 metres above ground level.f. t850 temperature at 850 mb (approximately 1500 m asl, close to the Juneau Icefield plateau).g. vmid vertically integrated moisture difference.h.tp total precipitation.

Supplementary table 2. Climatic data for each period from the Juneau meteorological station 1 . Summer temperature anomaly is compared with the 1986-2005 AD mean.
. By 1794 AD, Taku Glacier had already receded across Taku River 27 and Taku Inlet was ice-free for most of the 19 th century18.It was surveyed in 1890 AD by the United States Coast and Geodetic Survey, and was well inland of the 1948 position then, calving into Taku Inlet 17 .Taku Glacier advanced to a new position at the mouth of the valley by 1948, and then was then slowly advancing for most of the twentieth century26,30 , advancing 7.3km between 1890 and the early 1990s 27 .Calving into Taku Inlet decreased in the early 20 th century, as a moraine shoal rose above sea level.Calving had largely ceased by 1952.The glacier is now separated from tidewater by terminal moraines and outwash deposits 27 .Taku Glacier has only recently began to recede again, with a negative mass balance recorded from AD 2013 onwards 26 .Herbert Glacier reached its LIA maximum, where it formed a substantial six-meter high moraine, in AD 1765 17 (Supplementary figure1h).Significant recession (~2 miles) had occurred by 1948 17 .Eagle Glacier moraines were loosely dated to AD 1785-1787 by dendrochronology published in 1950 17 , with recession of ca.0.7 miles noted from this position in 1909-1910 20 and further recession of 0.7 miles by 1948 (ibid).Both these glaciers developed ice-marginal lakes during 20 th century recession16, leading to increased rates of recession due to calving.A similar history was noted at Mendenhall Glacier (near the town of Juneau), where the glacier in 1948 had receded from its LIA maximum (AD 1769) by two miles.Lemon Creek Glacier also formed LIA moraines dated to c. AD 1750 by dendrochronology 31 .These moraines demarcate the Holocene maximum extent of this glacier, 2.5 km downstream of the 1958 position.

figure 6. Sources of imagery used in A) 1948 and B) 1979 timeslices. Note that some data gaps in the USA imagery meant that holes were filled using the Topomaps. C) The different USGS Topomaps used. See Supplementary table 2 for more information. D) the 1948 ariel photograph mosaics over the USA. Aerial photographs are available from USGS Earth Explorer (see 'Data Availability'). New orthomosaics are available in Mendeley Data (see 'Data availability'). Historical topographic maps are available from the USGS map viewer (https://ngmdb.usgs.gov/topoview/viewer/#4/52.19/-123.71) as georeferenced GEOTIFFS. AHAP: Alaska High Altitude Photography. Supplementary figure 7. Imagery used in the reconstruction of glacier area in each timeslice. AHAP: Alaska High Altitude Photography. The source imagery for each single glacier polygon is provided in the attribute information of the shapefile, and summarised below in Supplementary table 3. Supplementary table 3. Source imagery for 1948 glacier outlines. AHAP: Alaska High Altitude Photography. Source imagery (1948 glacier outlines) Number of glaciers Area of glaciers (km 2 ) % Glacier area
photograph imagery was more widely available in 1979, with Landsat 3 satellite data used for gapfilling where needed (Supplementary table 4).A list of all satellite imagery used is available in Supplementary table 6.The date of the source imagery is recorded in the attribute information of each glacier polygon, and this is accounted for (number of days) when calculating annualised rates of change.

Supplementary table 6. Sources of satellite imagery used in this study
Analysis of uncertainty in glacier areaMethods of deriving an uncertainty estimate of glacier area are explained in detail in Methods, and are illustrated below in Supplementary Figure8and in Supplementary Table7.

Supplementary table 9. Examples of outlet glacier length and area change between the LIA and 2019 AD. LIA: "Little Ice Age".
In total, 91.04% of glaciers receded from AD 1770-1948 AD, though this is a minimum estimate as it was not possible to accurately reconstruct LIA glacier extent for all glaciers.Sparse historical records 17 suggest relatively slow but consistent recession from the LIA to ~1900, (Supplementary table 9), with accelerating recession from the early 1900s to 1948 and then onwards into the late twentieth century.

table 13. Summary dataset of icefield-wide snowline elevation for end of summer snowlines from 2019 to 2023.
Glacier snowline data were manually mapped from end-of-summer cloud-free Sentinel 2 satellite imagery (black and white, band 4 images) annually from 2019-2023 AD (Supplementary figure13; Supplementary table13).For each glacier, the mean elevation of multiple snowlines from one image is calculated.Mean end of summer snowline elevations for key named glaciers are shown in Supplementary table 14 and the full dataset is provided in the Source Data and in the shapefile of glacier snowlines in the Supplementary Data.The location of snowlines with respect to glacier outlines and the plateau is shown in Supplementary figure13.