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

Domain: earth_science
Published: 2024-07-02
Editorial outcome: major_revision

Abstract
Globally, glaciers and ice fields 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.

Introduction
Globally, mountain glaciers and ice caps dominate loss of land ice1, and were responsible for 21% of observed global sea level rise from 1993–20171. The largest contributions (8 mm from 1961–2016) were from Alaskan and Western Canadian glaciers1,2, which hold a large volume of ice (46.4 ± 15.3 mm sea-level equivalent)3, and which will contribute to sea level rise for centuries to come4,5. Alaska will likely remain the largest regional contributor to 2100 AD6. Model projections indicate that, with emissions policies defined under the Paris Agreement, land ice will drive 25 cm (11, 40 cm at 5th, 95th percentiles) sea level rise by 2100 AD. Of this, 13 cm will come from glaciers, 25% of which will come from Alaska alone4. Under the higher-emissions scenario RCP4.5, around 30% of glacier ice will have disappeared from Alaska by the year 2100 (relative to 2015 AD), causing 17 ± 4 mm of global sea level rise7.

Alaskan icefields may be particularly vulnerable to accelerated melt as the climate warms. Firstly, a large ice volume in Alaska is hypsometrically top-heavy or plateau icefields (cf. 3, 8, Fig. 1) with a low-slope accumulation area, making them likely vulnerable to small changes in equilibrium line altitude (ELA) 8,9. Secondly, a decrease in height of a flatter icefield or ice cap will theoretically result in a decrease in surface mass balance (SMB) as the ice surface lowers into warmer air, which reinforces itself as an SMB-elevation positive feedback10. Thirdly, due to an absence of higher topography, flatter ice caps and icefields cannot retreat to higher elevations and find a new equilibrium11. Icefields therefore are likely to exhibit threshold behaviour10,12, with tipping points for rapid and irreversible recession. However, a paucity of long, multi-decadal empirical observations of icefield change impedes our ability to accurately constrain and evidence these processes and identify tipping points and thresholds of change. Understanding the physical processes that affect the different way icefields and ice caps respond to climate change is key to improving our ability to estimate icefield response to future climate change13–15.

Upskilling our ability to predict future glacier behaviour using model-data comparisons requires detailed observations for model parameterisation and calibration or weighting of projections16,17. In Alaska, the largest Holocene glacier expansion occurred during the “Little Ice Age” (LIA) 18, AD 1770-1850 AD18–20, with temperatures 1.3 °C19 to 2 °C20 cooler than modern (Supplementary Information Section 3). The geomorphological imprint of this expansion is clearly observable in Sentinel satellite imagery and the fine-resolution ArcticDEM product21. Temperature changes during this historical period were of a similar amplitude of forcing (±2 °C) to those projected in the next two centuries22. Datasets of glacier change through this historical period, with a measure of uncertainty, are required for understanding the behaviour of ice caps and icefields under a changing climate, and for upskilling projections through data-model comparisons4,14,16,23,24, but few exist with a long (more than a few decades) and high (decadal) temporal resolution. The aim of this study is therefore to quantify “Little Ice Age” to recent (2020 AD) glacier change at Juneau Icefield (Alaska and British Columbia, Figs. 1, 2a) to determine the rates of, and controls on, icefield-wide shrinkage and thinning.

Results
LIA glacial reconstruction
The reconstructed LIA icefield (Figs. 3, 4 and Supplementary Figs. 9–10) covered 5414.95 ± 111.76 km2; 29.5% larger than in 2019 AD (Table 1), with reconstructed ablation areas adding 95.26 ± 17.01 km3 compared with 2013 AD (Table 2). At this time, Juneau Icefield outlet glaciers extended down onto the plain in several valleys, with many forming large unconstrained piedmont lobes; for example, Hole-in-the-Wall and East and West Twin Glaciers, Mendenhall Glacier, Herbert Glacier, Llewellwyn Glacier and Eagle Glacier (Supplementary Fig. 10). At these piedmont glacier lobes, large concentric moraines represent a substantial period of stability and moraine building for the glacier. In many of these places, the moraines display multiple closely spaced ridge crests. There are some moraines visible between the innermost of these concentric piedmont glacier moraines and the 1948 glacier outlines, for example at Herbert Glacier (Supplementary Figs. 1h, 9), suggesting that the ice may have temporarily stabilized during recession into the valleys, potentially due to changing topographic constraints on the margin. Recessional moraines are also evident at Llewellwyn Glacier, but within the 1948 extent.

For the most part, over-deepenings inside the main moraine crests are infilled with proglacial lake water. Taku and Norris glaciers were the only two glaciers to reach the ocean and to calve directly into a marine fjord. East and West Twin Glacier and Hole-in-the-Wall Glacier drained into an ice-dammed lake, but the remaining piedmont glaciers were land-terminating. Kettle lakes occur frequently within the piedmont glacier moraine crests. The remainder of the outlet glaciers (e.g., Field, Gilkey and Meade glaciers) remained constrained as land-terminating valley glaciers at their most recent neoglacial maximum (Figs. 3, 4, Supplementary Fig. 10). These glaciers terminated on land and had a long, narrow, low-slope glacier snout. Again, as the glaciers have receded from their most recent maximum position, proglacial lakes have formed within the over-deepening, impeding observation of younger recessional moraines.

Taku Glacier shrank by 6.4% (49.95 km2) but was slightly thinner at the terminus during the LIA than today, potentially reflecting its different configuration and marine-terminating environment, compared with the thicker, terrestrially terminating glacier present today25,37. Recent recession now sees the formation of lagoons at the ice margin and the onset of a lacustrine or lagoonal terminus. At the LIA, Taku Glacier was coalescent with Norris Glacier in Taku Inlet28, calving into the fjord (Supplementary Fig. 11). It is likely that this marine-terminating glacier had a long, low-slope, low-elevation snout, as is typical for water-terminating glaciers. We suspect that the glacier was grounded, given the shallow fjord, though this is not constrained. Behind the terminus, we reconstruct an ice-dammed lake at 15 m elevation, following the clear break in slope between the forested mountainside and the low-slope boggy valley floor, as originally mapped by Lawrence, 195028. This is in line with ice-contact fans observed at Twin Glacier Lake moraines and at Wright Glacier terminus further up-valley (Supplementary Fig. 10), which provide additional geomorphological evidence for an ice-dammed lake at 15 m asl. However, little geomorphological evidence of shorelines is visible in the ArcticDEM or satellite imagery. Reconstructing the lake following the 15 m contour generates a lake 95.84 km2 in size.

Around the periphery of the main plateau, the LIA extent of valley glaciers and mountain glaciers is generally well constrained by substantial and clear moraines, enclosing ice-scoured bedrock or fluted glacial sedimentary surfaces. There are occasional recessional moraines within these larger moraines. For some glaciers, such as Denver Glacier, limited clear geomorphology visible on satellite imagery precluded reconstruction of the glacier extent, and so the 1948 extent is used as the earliest data source.

Reconstructed geometric equilibrium line altitudes (ELALIA) ranged from 1070 m (Taku Glacier, Lemon Creek Glacier) to 932 m (Mendenhall Glacier) and 1494 m (Llewellyn Glacier) (Supplementary Fig. 10). In comparison, the average ELA for Lemon Creek Glacier from 1998-2019 was 1327 m asl25,35. This represents an ELA lowering of ~150 m at Lemon Creek Glacier, relative to today25,35. Peripheral glaciers to the east of the icefield had an ELALIA higher than 1500 m, compared with under 1250 m in the west (Supplementary Fig. 10).

Overall changes since the “Little Ice Age”
Between the AD 1770 and AD 2019 inventories, 108 glaciers disappeared (Fig. 5), and only 70.53% of glacier area now remains relative to the Neoglacial maximum. 100% of glaciers mapped in 2019 had receded relative to their LIA position, with 47 new ice-contact proglacial lakes developing as glacier termini receded. The largest areal losses (52.5%) were from the Icefield’s outlet glaciers (which made up 64.25% of the area during the LIA). Of all the glacier area lost between 1770 and 2019 AD, 538.0 ± 32.4 km2 was from outlet glaciers. The main outlet glacier termini all receded by ~4 to 5 km, shrinking by up to 69 km2 (Supplementary Table 9).

The total cumulative ice loss from across Juneau Icefield is 315.3 ± 237.5 km3; this equates to a loss of 24.25% of the ice volume (from ref. 3) between the LIA maximum and 2020 AD. There has been thinning of 150 m at the terminus of Field Glacier, 220 m on Ogive Glacier, 150 m on Gilkey Glacier, and 200 m on Meade Glacier. Bacon Glacier has thinned by 190 m, and Tulsequah by 180 m.

Accelerating change through the 20th and 21st Centuries
Most glaciers (91%) receded substantially from their LIA extent between 1770 and 1948 AD (Fig. 5), during a period of warming (Fig. 6a). Total glacier area shrank by 12.18% at a rate of 0.07% a−1 (3.70 km2 a−1). This was accompanied by a conservative estimate of volume loss of 0.65 ± 0.92 km3 a−1 (Fig. 6b–g). Between 1948 and 1979 AD, 342 (30.6%) glaciers advanced, and overall rates of glacier recession were slow (0.12% a−1) (Fig. 6b–g; Table 1). This was likely caused by an increase in snowfall relative to the period before 1948 (Supplementary Table 2, Fig. 6f). Most advancing glaciers were either small mountain glaciers or glacierets, likely reflecting their faster response time. Despite these changes, a comparatively low rate of thinning persisted, at an icefield-wide average of 0.11 m a−1 (Fig. 6c). Rates of recession before 1948 remain well below those experienced by most outlet glaciers after the year 2005 (Supplementary Table 10).

Slower rates of areal loss and thinning continued until 1979, followed by an acceleration of ice loss during the latter 20th and start of the 21st centuries (Table 1). Fewer (8.5%) glaciers advanced between 1990 and 2005 (n = 95). Rates of glacier area loss rose, reaching 0.18% a−1 (8.33 km2 a−1, 1979–1990), and then sharply accelerated to 0.39% a−1 (17.54 km2 a−1, 1990–2005). From 1979–2000, icefield-wide volume loss reached 3.7 ± 1.6 km3 a−1 (Fig. 6d), also indicating a sharp acceleration from the preceding periods, with a mean thinning of −0.48 m a−1 across all glaciers (Fig. 6c). Glacier tongues thinned across the icefield (Fig. 7), with thinning reaching higher elevations up-glacier than before (Fig. 7). Thinning reached elevations of 1380 m asl on Mendenhall Glacier, 1490 m for Tulsequah Glacier, and the plateau at 1550 m for Meade Glacier.

This rate of ice loss continued into the 21st century. Between 2005–2015 AD, 99% of glaciers receded, and 10% of glaciers disappeared between 2005 and 2019 AD. From 2005 onwards, glacier areal loss accelerated sharply relative to previous time periods, reaching 0.60% a−1 (25.36 km2 a−1, 2005–2015) and then 0.96% a−1 (38.47 km2 a−1, 2015–2019) (Fig. 6b–g). Rates of recession were 1.5 times faster from 2015–2019 (Mean (M) = 5.88, Standard Deviation (SD) = 5.26) than from 2005–2015 (M = 2.48, SD = 2.11), t(1097), p < 0.001.

Icefield outlet glaciers experienced severe volume loss after 2010 (Fig. 6c, g; Table 2), with thinning across the icefield plateau at elevations of up to 1800 m asl (Fig. 7) and volume loss at a rate of 5.91 ± 0.8 km3 a−1 (2010–2020). The whole icefield loss equated to 5% of the total volume reconstructed in 2017–18 by ref. 3 from 2010–2020. Alongside this thinning we observe increased glacier fragmentation, with both terminus separation and glacier disconnections (sensu ref. 21), where glacier tributaries are detached from their accumulation areas, occurring increasingly frequently (Fig. 6h, i).

Overall, rates of glacier areal recession rose slowly from 1948–2005, followed by a sharp acceleration in rates of glacier area loss after 2005. Icefield-wide, rates of glacier area shrinkage were 7 times faster from 2015–2019 relative to 1948–1979, and 5 times faster relative to 1979–1990 (Table 1, Fig. 6b–g). There are three distinct periods in icefield volume loss, with relatively consistent rates from LIA-1979, followed by a rise from 1979–2010, and a doubling of rates of volume loss after 2010 (Fig. 6d). There has thus been a steep acceleration in the rate of ice mass loss across the icefields since the early twenty-first century.

Snowline and Albedo change
This rapid glacier recession was accompanied by ELAs rising to the plateau height in the period 2005–202025,35 (Fig. 6i), and exceptionally high snowlines on the plateau in 201838 and 2019 (Figs. 5, 6k, Supplementary Table 13; Supplementary Fig. 14). Late summer snowlines mapped in 2019 from Sentinel-2 imagery (Figs. 5, 6k) had a mean altitude of 1612 m asl (SD 162 m; n = 178). Snowlines on Taku Glacier were a mean altitude of 1445 m in 2019, compared with a measured ELA in 2019 of 1528 m35. The snowline here is therefore well above the southern plateau at 1200 m asl. No snow was observed in the Sentinel-2 imagery on Lemon Creek Glacier at the end of the ablation season in 2019, in line with the measured ELA of 2023 m that year35. Snowlines on Meade, Llewellyn and Field Glacier had a mean elevation of 1671 m (SD 119 m), well above the height of the northern plateau at 1500 m asl (Fig. 5). This icefield-wide trend has continued in subsequent years (Fig. 6k, Supplementary Tables 13, 14 and Supplementary Figs. 13, 14). These datasets show that, although snowline elevation was unusually high in 2019, snowlines are reaching the plateau frequently, with both snowlines and ELAs regularly occurring above the plateau height of 1200 m since 2005 and increasingly frequently in the last decade of the study period.

Alongside this trend of rising snowlines, we also observe a decrease in icefield wide average albedo (Fig. 6l, Supplementary Tables 15–17). When clipped to the 1990 glacier outlines, and thus taking into account terminus recession and icefield fragmentation as well as the impact of darkening snow and ice surfaces, the mean albedo from 1987–2009 was 0.81 ± 0.03 (95% confidence interval). From 2010 to 2023, this decreased to 0.67 ± 0.03. The albedo of the plateau above 1500 m (compared with the Copernicus DEM) also decreased from 0.92 ± 0.02 (1987–2009) to 0.78 ± 0.04 (2010–2023). The months of September and August in 2018 and 2019, which were characterised by exceptionally high snowlines, likewise had exceptionally low icefield-wide mean albedos of 0.44 to 0.48, respectively (LC08 sensor). In contrast, the lower icefield average snowline in August 2021 (Supplementary Table 16) is reflected with a higher albedo of 0.66 in the LC08 sensor.

When comparing the time periods 1987–2009 and 2010–2023, the icefield-wide albedo within the 2019 glacier outlines was lower from 2010–2013 (mean = 0.71, standard deviation = 0.14) than from 1987–2009 (mean = 0.84, standard deviation = 0.08); t(90) = 5.48, p < 0.05. This decrease is caused by more pixels with a lower albedo in summer, especially on the plateau; here, a regression for average albedo over time results in an r2 value of 0.23 (p < 0.05). Many bright pixels with high albedo however remain throughout the study period, especially at higher altitudes on the plateau.

Icefield fragmentation
The recession and thinning of Juneau Icefield glaciers has resulted in increasing icefield fragmentation. This includes the separation of glacier tongues in the valleys, and disconnections higher up the glacier, where bare rock appears within the glacier polygon as a result of thinning, usually over thin, steep and heavily crevassed ice such as in icefalls21. While 291 disconnections were first mapped in ref. 21, which were frequently observed in association with icefalls, here we calculate the timing of disconnection by comparing the mapped disconnections with the new glacier outlines, and analyse the thinning occurring at the site of the disconnections. The mean elevation of the disconnections has remained relatively constant over time (Fig. 6i), with a total mean elevation of 1299 m asl (1362 m in disconnections after 2005); this is controlled by the location of steep ice around the rim of the plateau21. This altitudinal range (1200–1400 m asl) was subjected to thinning in the period 2010–2020 (Fig. 7). For example, in a subsidiary accumulation basin on Gilkey Glacier, thinning of 1.4 m a−1 is observed at a height of 1340 m adjacent to a glacier disconnection. At Thiel Glacier, thinning of 1.7 m a−1 is observed at an altitude of 1034 m adjacent to the site of a pre-2019 disconnection.

These disconnections are occurring at a uniform altitude and slope through time (Fig. 6i) but are occurring more frequently after 2005 AD as glacier ELAs increasingly intersect this altitudinal range (Fig. 6h; refs. 21, 35) and thinning occurs at higher elevations (Fig. 7). Glacier snowlines in 2019 AD were above the elevation of icefalls around the rim of the plateau21 (Fig. 5), resulting in net ablation at the height of many of the icefalls. Glacier tongues down-stream of glacier disconnections show increasing debris cover and thinning. For example, thinning of up to 9 m a−1 on Thiel Glacier tongue, with a glacier-wide mean of −3.33 m a−1, is far above the mean of 1.3 m a−1 for outlet glaciers (2010–2020, Fig. 6c).

Methods
Baseline glacier outlines were derived from the Randolph Glacier Inventory (RGI) v. 6.066,67. We used glacial geomorphological mapping (ref. 21 and reviewed in Supplementary Information Sections 3 and 5) and historical records28 to reconstruct LIA extent and then volume68,69. We then applied structure-from-motion algorithms in MicMac70,71 to 20th Century archival aerial photographs (1948 and 1979 AD) as well as topographic maps to reconstruct glacier area and surface elevation in 1948 and 1979 (see Supplementary Fig. 6 for data coverage). Satellite imagery was also used in 1990, 2005, 2015 and 2019 to altogether yield glacier outlines from 1770, 1948, 1979, 1990, 2005, 2015 and 2019 AD (Supplementary Fig. 7), and original icefield surface DEMs at ~1770, 1948, and 1979 AD. We combined these with datasets of ice surface elevation derived from structure-from-motion algorithms using ASTER imagery from 2000, 2010 and 202072.

Remote sensing datasets
Glaciers in 2019 AD, glacial structures, glacial geomorphology and lakes were used from ref. 21. These features were mapped in ESRI ArcGIS (projection WGS 84, UTM zone 8N), using 10 m resolution Sentinel-2 imagery (swath 290 km), and these data and methods are published in full in ref. 21; methods and datasets are summarised here for completeness. Late-summer satellite images with limited cloud cover and snow on the glacier ablation regions were chosen. Landsat imagery was used to map glacier extent in 2015 (Landsat 8), 2005 (Landsat 7 ETM+), 1990 (Landsat 4–5 MSS), and 1980–1982 (using a combination of archival aerial photographs and Landsat 3 MSS imagery where there was no coverage). Landsat imagery was acquired pre-registered to Universal Transverse Mercator (UTM) World Geodetic System 1984 ellipsoidal elevation (WGS84) zone 8N projection. Image resolution ranged from 4 m (aerial orthomosaics) to 10 m (Sentinel-2), and 30–80 m (Landsat).

Stereo aerial photographs from the 1979 Alaska High Altitude Photography (AHAP) survey were used to map glacier extent and structures from 11th–12th August 1979, especially in areas where significant glacier changes have occurred. This campaign acquired black and white photography at 1:120,000 scale resolution. Cloud cover is minimal and snow patches are limited in area and scope. These photographs were scanned by the USGS EROS Data Centre at 25 microns (1200 dpi).

Aerial stereo photographs from the 1948 aerial photography survey of Alaska were used to map glacier area. Photographs were acquired on 13th August 1948. USGS topographic maps dating from 1948 were used to supplement the 1948 aerial photographs, where the orthomosaic was poor, or where original aerial photographs were unavailable on the USGS server (Supplementary Fig. 6). This was mostly the case in Canada. Here, USGS maps (AK_Atlin_360559_1960_250000) were used, where the topographic maps were deemed sufficiently high in quality (see Supplementary Tables 3, 5), with a good correspondence between mapped lakes, rivers and topography. Where topographic maps showed a low correspondence, the 1979 AHAP mosaics were used, and, where these were also unavailable (in the eastern-most parts of the study area), Landsat imagery was used (Supplementary Fig. 6). Overall, 66.6% of the glacier area was mapped from the 1948 orthomosaics, and 27.8% of glacier area was mapped from topographic maps (Supplementary Table 3). This means that the glacier extent in 1948 is the best available estimate given the data availability. Annual rates of recession are calculated using the date of the imagery or data source used to map glacier extent. The glacier extent in 1948 is considered a minimum; use of later datasets results in some glaciers being underestimated due to continued glacier recession.

Topographic data were derived from the 2 m resolution ArcticDEM v3.0 release 773, except where tiles were missing, and ASTER GDEM74, ALOS DEM75 and novel DEMs derived from the aerial photographs using structure-from-motion methods. A derived composite DEM is mostly ArcticDEM, with data gaps being filled with a void-filled 30 m resolution ASTER GDEM v2.074. For this process, both DEMs were projected in UTM Zone 8N, and clipped to the same extent. All ArcticDEM ‘no data’ values were converted to 0. The ArcticDEM holes were then filled with the ASTER GDEM data (filled output upsampled to 2 m cell size to match ArcticDEM) using raster calculator. Comparison of the composite DEM with the original, unmodified GDEM raster yield differences for individual pixels of up to 15 m, which is within GDEM resolution. This likely stems from the resampling carried out as part of the final raster calculation. Hillshade and slope models were derived from the composite DEM mosaic to aid interpretation of glaciological features.

For some glaciers, historical maps and surveys prior to 1948 are available. Here, the maps provided in historical sources (e.g., ref. 28) were scanned, georeferenced, and the termini digitised and included in the available data. As only terminus positions are available, these are provided as polylines at the glacier terminus, annotated with the year (e.g., Fig. 2E, Supplementary Fig. 5).

Photogrammetry and generation of orthomosaics and DEMs
Aerial photographs from surveys flown in 1948 and 1979 are available in Earth Explorer (USGS). They are scanned images, without any geospatial data attached. The files are not georeferenced but the database includes an approximate footprint for some images. The surveys were undertaken in July–August 1948 and from a single flight in August 1979, during a period of minimal snow and cloud cover. We selected all available images from the Earth Explorer database (1446 from 1948 AD and 103 from 1979 AD) that intersected a polygon drawn around our study region, along with images that had no footprints in the database but were part of the same flight roll, and that were not damaged, torn or had holes.

For 1979, we used the calibration reports provided by NASA along with the fiducial markers to resample the images to a consistent size. The 1948 images do not have calibration reports but do have fiducial markers, so we used the average position of the fiducial marks to resample the images. We then applied contrast limit adaptive histogram equalization (CLAHE) to ensure the images had consistent brightness. Using MicMac71, we found tiepoints and calculated the relative camera orientations for each image and processed a relative orthophoto and DEM, using camera coordinates rather than geographic coordinates. The 1948 images were processed in four separate blocks due to the number of images involved. We then automatically registered the resultant orthophotos to a Sentinel-2 image to find ground control points. Elevation was determined using the Copernicus 30 m DEM (https://spacedata.copernicus.eu/collections/copernicus-digital-elevation-model). From these control points, we extracted the DEMs for each of the four blocks in 1948 and in 1979, using MicMac. For 1948 we co-registered to the Copernicus DEM and cropped the edges of the DEM to deal with edge effects, before stitching the DEMs together, using an established approach76. In post-processing, we masked the errors (largely in the accumulation area) and these are given as ‘no data’. Voids are filled using the local mean elevation change method77.

DEMs of difference were then produced for LIA-1948, 1948–1979, 1979–2000 AD. Volume change from 1979–2000 uses the STRM. DEMs for volume change analysis from 2000–2020 are derived from ref. 72, using ASTER imagery and ArcticDEM. DEMs were resampled prior to differencing to the extent of a Copernicus DEM over the study domain, which includes resampling to the 30 m pixel size of the Copernicus DEM. This is undertaken using bilinear resampling. The 1948 and 1979 DEMs are therefore upsampled, as they are higher resolution than the Copernicus DEM.

Glacier area determination and uncertainty
Initial ice divides and glacier outlines were downloaded from the Randolph Glacier Inventory (RGI) version 6.066,67, which has a census date of AD 200578. This historical inventory is our baseline from which recent changes are assessed. The RGI outlines were updated and removed since some were incorrect; they were misclassified snow cover or shadow, or inappropriately subdivided into glaciers. Some glaciers were missing, or were incorrectly delineated, or had overlapping outlines with another glacier polygon. These errors are clear when comparing outlines to the higher resolution, more recently acquired satellite imagery. RGI 6.0 glacier outlines were therefore manually edited to be consistent with Landsat 7 ETM+ imagery from 2005 (updated outlines published21, see Supplementary Table 6). Glacier outlines were then manually edited in ArcGIS according to observed changes around the glacier margin from the satellite imagery, aerial photograph orthomosaics, 1948 USGS topographic maps, both in the terminus and forefield, and at higher elevations, to reflect glacier extent in 2019, 2015, 2005, 1990, 1979, and 1948 AD. RGI IDs were maintained for each glacier polygon. Where glaciers had separated into separate polygons, the same RGI ID was used for both and the glacier was treated as a multipart polygon in ESRI ArcGIS. The high pixel resolution of the recent Sentinel-2 images and archival aerial photographs allowed a more detailed glaciological analysis and smaller minimum glacier and lake area for mapping than earlier inventories, which relied on coarser-resolution Landsat imagery78. The minimum glacier and lake size mappable was 0.001 km2. The date for each glacier is noted in the attribute information and is used to calculate annualised rates of change. Glaciers were categorised as glacierets, mountain glaciers, valley glaciers or outlet glaciers, following GLIMS guidelines79 and previous work21.

To generate glacier outlines for the latest neoglacial maximum (the LIA), the RGI v6.0 200567 outlines were manually edited to extend to the crests of LIA-associated moraines, bare ice-scoured rock and associated trimlines, using published geomorphological mapping21 and published chronological and geological data28,32,34,80,81, which are reviewed comprehensively in the Supplementary Information (e.g., Supplementary Fig. 5). The more subtle geomorphological evidence included stripped bedrock at elevations beneath surrounding vegetated hillslopes, and locally-depressed treelines. Glacier moraines at Herbert, Eagle and Twin Glacier Lake were ground truthed and field-checked during fieldwork in July 2022. Published ages were gleaned from the literature and included in the geodatabase. Radiocarbon ages were recalibrated using CALIB82. Calibrated mean ages (2σ) are presented as ‘cal. ka BP’. Morphostratigraphic principles83–85 were applied to differentiate between the most recent, LIA Neoglacial moraines, and older moraines deposited in the Younger Dryas or earlier Holocene glaciations. In most cases, Holocene neoglaciations are likely to have been within the limits of the LIA18, which makes confusion of older ice limits less likely. Where it was unclear which moraine related to the LIA, the innermost was chosen, to provide a conservative estimate of LIA area.

Uncertainty in glacier area was calculated following the methods outlined by ref. 86. Eleven representative glaciers, across the spectrum of glacier sizes and including typical examples of debris-covered and clean-ice glaciers, were digitised seven times, with a minimum of one day between each round of digitising (Supplementary Fig. 8). The mean, standard deviation and 95% confidence interval was then calculated for each glacier for each timeslice. There was a significant correlation between glacier size and confidence interval (Supplementary Table 7). A regression through the data points provided an equation that was used for size-class specific up-scaling to the full dataset66,86,87. As a comparison, we also calculated the standard error and standard error of the forecast for each timeslice. This gave a similar, but smaller, estimate of the uncertainty, and so the confidence interval, the more conservative estimate of uncertainty, was used.

Glacier recession, defined in terms of area change per annum (km2 a−1) and proportional rates of area change (% a−1), was calculated for each glacier. Note that where we refer to terminus retreat (i.e., a glacier length change) this is specifically stated. Annual rates of area change were calculated by dividing area change (km2) by the time between analyses for each glacier (time was taken from date of satellite or aerial photograph image acquisition). This allows results to be compared between different times and different glacier areas. For 1948 and 1979, where the same imagery is used for a few glaciers in each timeslice due to a lack of coverage, this means that rates of glacier recession are 0 and are removed from our analysis. However, all glaciers are included in the total summed areas given for each timeslice, and the total summed rates of change. Rates of glacier recession for individual glaciers (% a−1 and km2 a−1) are normally distributed. Mean rates are presented with the 95% confidence interval.

Snowlines
The late summer snowline from 2019 AD was mapped from the cloud-free Sentinel-2 imagery (Supplementary Table 6), using black and white band 4 imagery (in line with structural mapping conducted in ref. 21). These images, from 30th August to 7th September 2019, provided clear, cloud-free coverage of the entire icefield. The short timeframe between image acquisitions allowed temporal continuity and clear visualization of icefield snowlines. Snowlines were identified as a wavy, discontinuous boundary between bright white and duller grey; the higher albedo snow reflecting a brighter white colour.

Glacier disconnections
Glacier disconnections, where the tongue of a glacier had become separated from its accumulation area62,63,88, were mapped across the icefield published in ref. 21 without a date stamp. Individual disconnections were mapped that occurred prior to each timeslice. We classified glacier tongues as ‘disconnected’ if they were associated with one or more of these disconnections from the accumulation area. In the case of disconnected glaciers, multipart polylines were used to record glacier length along the same flowline. The timing of disconnection was calculated through comparison with the glacier area outlines, and the slope and elevation of the disconnection calculated through automated tools using tools within the GIS. Elevations were derived from ASTER GDEM74.

Surface albedo
We derived a stack of 299 Landsat 5–9 images in Google Earth Engine for the years 1986–2023, limited to WRS path/row 57/19 and 58/19, for all available scenes with less than 20% cloud cover on land (Supplementary Table 15). Clouds were masked using a cloud cover filter using the USGS QA_Mask, removing any pixels identified as cloud by the USGS algorithm. We calculated albedo for each pixel in the stack using the albedo equation developed for Landsat 789,90 bands blue, red, and infrared, and applied this across the TM, ETM+ and OLI sensors for Landsats 5–9. We then calculated average albedo for every scene within the icefield boundaries as of 1990, 2005 and 2019 AD (not including the peripheral ice masses). This captured the plateau and glacier tongues. Using the average albedos calculated within the 2019 outlines effectively analyses the darkening of snow and glacier ice only, whereas using the 1990 outlines captures glacier recession and fragmentation in addition to darkening snow and ice. We also calculated the average albedo in each scene for pixels above 1500 m in the Copernicus DEM (~2013 AD), this capturing just the plateau. Observations were limited to at least 60% of the plateau being covered and not masked out by cloud.

LIA volume reconstruction
Our LIA outlines were converted to lines with points along them at intervals of 100 m, and those points were used to extract elevations from the ALOS DEM75, which was selected for its seamless coverage. A surface to represent the LIA glacier ice surface below the ELA was interpolated between those points. Carrivick et al.68,69 have shown that there is a negligible effect of which DEM (resolution and date) is used in these types of analyses, because the points are preferentially located from the high-resolution 2 m ArcticDEM73 even if the elevation values are extracted from a coarser resolution DEM. However, Carrivick et al.68,69 have also shown that there is an effect on volume loss calculations of the choice of interpolation method and in this study, we followed their advice using a Natural Neighbour interpolation.

Due to an absence of moraines and trimlines above the equilibrium line altitude of glaciers, our mapping of a LIA glacier outline, our LIA ice surface and our analysis of elevation changes and hence volume loss only pertains to LIA ablation area. Following the method of Carrivick et al.68,69, LIA glacier-specific equilibrium line altitudes (ELALIA) were produced using the Area-Altitude Balance Ratio (AABR) method91. The AABR method was chosen to explicitly consider (spatially-distributed) glacier hypsometry. The glacier-specific ELA were subsequently merged with our mapped moraine crests and trimlines to create a single polygon for each glacier ablation area. We estimated the ELALIA using a BR of 1.88 for all glaciers as suggested to be typical of maritime glaciers92. Published work shows that the choice of balance ratio has little influence on glacier ELA93.

Volume change estimates presented in this study can only be considered as a proxy for mass balance if it is assumed that glacier surfaces above the ELA have not changed since the LIA. Thus, volume change using the methods in this work should be considered as minimum estimates because probably some mass loss has occurred above the ELA as they have risen in elevation due to atmospheric 20th Century warming. That mass loss above the ELA would likely be via surface lowering, rather than areal contraction, and so we would not detect it with our mapping.

Glacier volume change calculation
To determine volume changes (LIA-1948, 1948–1979, 1979–2000 AD), we used the LIA, 1948 and 1979 glacier outlines, using the outlines that correspond to the start date of the DEMs of difference for each time period. For each glacier, the outline was divided into elevation bands of either 50 m or, 10% of the glacier elevation range, whichever was smaller. In effect, glaciers with an elevation range of greater than 500 m were divided into 50 m elevation bands. The mean elevation change was calculated for each elevation band, removing outliers that were more than 3 normalised median absolute deviations (NMAD) away from the median of the band. We use the mean value for the elevation band if there is 25% or more coverage by area within the elevation band. If there is less than 25% coverage by area within the elevation band, then we use the icefield-wide average for that elevation band.

For each elevation band, this mean elevation change value was multiplied by the area of the band, to calculate volume change, per elevation band. Each band’s volume change was summed, and the sum of this is the total volume change for each glacier. For glaciers with no DEM coverage at all, we use the icefield-wide area-averaged thinning rate to scale up the total volume change for the region. This amounts to 44.6% of glacier area for LIA-1948, 16.8% for 1948–1979, 4.4% for 1979–2000, and 0% for 2000–2020 (Supplementary Table 8).

The sum of the total volume change for each glacier was then summed to calculate total volume change for each time period. Mean dh/dt per glacier was calculated by dividing the mean elevation change value by the time period. Mass change was not calculated. For LIA-1948, where pixels showed a positive elevation change, this was assumed to be erroneous, due to the substantial area change, these pixels were masked out. The LIA-1948 elevation change will be a minimum, because the LIA DEM only covers ablation areas.

We use the normalized median absolute deviation of off-glacier pixels to estimate the uncertainty in the DEMs of difference. To estimate the uncertainty in elevation change (σ dz) introduced by interpolating elevation bands with no data, we adapt the approach utilised in ref. 33 in Eq. (1):

σ dz = N z (1 − f) + M + N z f (1)

Where N z is the NMAD of off-glacier elevation differences calculated between the two DEMs, f is the fraction of the glacier that required interpolation, M is the standard deviation of all on-glacier pixels.

We estimate the total uncertainty in each glacier’s volume change as the quadratic sum of the elevation change and the area parts (Eq. 2).

σ dv = (σ dz A) 2 + (σ A dh) 2 (2)

Where A is the glacier area, σ A is the uncertainty in glacier area (following the methods above), and dh is the area-averaged elevation change on the glacier.

Uncertainty is higher in the LIA-1948 period due to the interpolation of the datasets over areas with no data. The LIA DEM is calculated only over the icefield ablation areas, and the 1948 DEM is likewise not available over Canada. This does not mean that the change is insignificant, but we emphasis that the uncertainty is large because of interpolation, not because of variability of the elevation differences.

Volume change, and uncertainty, from 2000–2010 and 2010–2020 AD was derived directly from ref. 72. Glacier outlines were taken as the RGI v6.0, as explained in the original publication. Glacier-averaged mean volume change and mean elevation change through time (mean dh/dt) was calculated following the method above. Thus they do not take into account area changes through this time period and should be considered a minimum value.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of the manuscript “Accelerating glacier volume loss on Janeau Icefield driven by non-linear processes” by Bethan Davies et al.\nThe research leaded by B. Davies, shows an interesting case study on Juneau Icefield in Alaska about the evolution of the glaciated area and volume of the glaciers in this region. Manuscript authors apply state of the art remote sensing techniques to delimitate glacier extent in different time periods. They are also able to reconstruct glacier surfaces to determine volume differences across distinct temporal windows. The manuscript reads well and shows reliable results for distinct time periods that without any doubt show the fast decline of Juneau Icefield in the last decade compared to previous period with data.\nMajor comments\nDespite I have enjoyed the read of the manuscript I consider that there are some major points to be addressed before considering for publication:\n• The manuscript claims that non-lineal processes drive the volume loss on Juneau Icefield but shits is not really shown in the manuscript. The database demonstrates that volume and area losses have accelerated in the last years and that many glaciers have disconnected or separated. On the contrary processes that control and will control glacier evolution are enumerated and related in a schematic (figure 6). These processes have not been measured or evaluated with the data processed and thus, from my understanding some sections of the manuscript must be changed considering this (including the title). See minor comment for lines 247-248.\n• The volume and glacier extent for 1948 and for the LIA maximum, have some uncertainties not included in the analysis. For example the glacier extent in 1948 is based in a 66.6% in images of this year, but the glacier area not observed in this year is inferred from topographic maps, and satellite images from 1970 to 1982. Similarly, LIA extent. Similarly, the uncertainty in area extent at the LIA maximum is of same order of magnitude that the extent in 2005, 2015, 2019, periods for which high spatial resolution satellite observations were available. Despite the LIA extent can be well identified from methods described from lines 250 to 271 of the supplementary information, an evaluation of the uncertainty of this extent must be included. Similarly the ELA reconstruction for the LIA might have an error assessment .\n• The readability of some figures must be improved. In some cases small circles avoid the identification of glacier extent in particular time periods. Similarly, some arguments are directly related with figures panels that do not have a correct identification doing difficult to follow the discussion. See minor comments for the figures.\nMinor comments:\nLine 18-19 an 34: In the second sentence of the abstract it is stated “The majority of the ice volume loss comes from Alaskan icefields,…”. I see a bit biased this previous sentence when in line 34-35 it is stated that ” with 13 cm coming from glaciers, a 25% of which will come from Alaska”. Definitely it is not the majority of volume loss. Change conveniently.\nLine 26-27: While it is true that that non-linear processes have to be incorporated into glacier evolution models; this is not a direct conclusion of the work. I would remove this last sentence of the abstract, and show a more general conclusion obtained on the frame of this study.\nLine 64: I have found quite long the supplementary information section about LIA ice extent and the climatic conditions observed in this period. I would encourage to reduce the extent or even completely remove these sections and just include appropriate references to this period. Include reference to research articles focusing in this time period and remove the reference to supplementary information.\nLine 84: Reconstruction of LIA icefield ablation areas: Provide a short explanation on how was obtained\nLine 84: Include uncertainty in the total ice volume estimation for LIA reconstruction “±XX km3”. Same for total cumulative ice loss (315.31km2±XX km2 )\nLine 110: The extent of 17.54km2 was an area decrease, right? It is not clear as in this sentence, it has not been computed as an area loss. In the next sentence it is clearly identified the 3.72 km as a loss. Please, along the manuscript, state clearly when loses are observed, eider with a “-” before the number, or clearing stating that a “decrease” or a “loss” has been observed.\nLine 118: M and SD are?? Mean and standard recession rates? Per year? Maybe I have missed the introduction of these acronyms, but I am not able to find them. Please, include them conveniently. Same for “p” and “t”\nLine 120: A volume loss of -5.91km3 might be volume gain. Please see comment for line 110 and change conveniently.\nLine 126: Rates of glacier shrinkage of area or volume? Please specify.\nLine 129: These doubling rates of volume loss after 2010 can be observed in Figure 4? I think they are, please, include which panel in figure 4 show it.\nLines153-168: The discussion is in some sections difficult to follow as figure 4 panels are not conveniently included.\nLine 180: Also include here this reference: https://doi.org/10.1038/s41558-018-0296-5\nLine 204: The debris-covered areas isolation from the climatic evolution might be very important. I think this can be longer discussed in these lines. Moreover, as glacier continue their area decrease, debris-covered areas will increase thus their impact will be more evident.\nLine 224-227: Include references to state where and which are the ice fields in Canada, Greenland, …\nLine 247-248: It is stated that “as our observation demonstrate, once the ELA rises above the plateau non-linear processes operate”. However the only reference to the ELA and the non linear processes in the manuscript is this: “The critical driver of this non-linear and accelerating glacier decline at Juneau Icefield is that the rising ELA now intersects the elevations of the icefield plateau (Figure 6).” Figure 6 is an schematic figure, showing processes, but it does not demonstrate with observations that non-lineal processes operate.\nFrom my understanding the main message of this manuscript (including the title) cannot be sustained with this figure, and thus must be conveniently changed. Moreover, the ELAs included in the extended data 1, do not include any temporal information. These are recent ELAs?? Please, clarify the ELA’s date computation in extended data 1\nOnline methods:\nLine 423: Cloud cover patches are limited in area, but which %?\nLine 425-436: The 28% of the glacier extent in 1948 is mapped from the topographic maps from which year?, the 3.7% from 1979 AHAP, the 1.67% from 1982 Landsat… (all data from Supplementary table 2). For me, this is not the true glacier extent in 1948 (more than 30 years in between some acquisitions). This must be clearly stated in the manuscript, and the glacier areas computed for 1948, include a higher uncertainty showing the not negligible % of the study area from a different (and later) observations. A good example where this must be shown is in figure 3 glacier areas.\nA good approach to determine if images from 1979, 1981 etc (or the topographic maps) can be exploited to identify glacier limits in 1948, is to analyze in areas with data in 1948 and later years, which is the deviations between these extents.\nLine 442 GDEM data is ASTER GDEM?, please clarify.\nLine 502: I would only talk about LIA maximum without including a specific year (i.e. 1770) For the maximum LIA extent I do now see relevant + or – 20 years. Same for supplementary table 6.\nLine 530: And what did authors do with images from 1980, 1981 and 1982 (supplementary table 2), which are also used for computing 1948 glacier extent?\nLine 579-580: From previous sentence this is obvious, remove the sentence: “In effect, glaciers…”\nLine 583-584: “Then we” appears twice.\nLine 587: Which is the % of glaciers with no DEM coverage?, please include this information here.\nLine 590: Mean dhdt should be dh/dt?? Same in line 609\nFigures:\nFigure 1: I encourage including here a map/satellite image showing the entire study area as this is the first figure of the manuscript and thus readers want to identify the main characteristics of the study site. I would remove panel D, which is not informative at this initial point of the manuscript. From this “complete” overview of Juneau Icefield I would delimitate with a rectangle where panels C and E are located.\nPanel A: Include boundaries for a better sea-land limit identification\nPanel E. Some contour lines (with transparency) might help to interpret the map. Glacier limits in 1931 and 1908-1910 are informative, but if I am right are not used in the manuscript, and thus are not needed. Disconnection separation marks are barely identified in panel E. I encourage to include here not only circles with distinct colors and use triangles, squares etc…for easier identification.\nFigure 2: The number of polygons (in brackets), is not informative of the number of glaciers, I would remove this information and change it by the total extent of the ice-dammed lakes, alluvial fin,… ELA (LIA) circles hamper the identification of the estimated LIA glacier extent. For me is more informative the LIA glacier extent (dark right) than the ELA estimation. I would remove all circles in A and include a “C pannel”( or include a new figure in the supplementary material), just including ELA (LIA) with features of different colors (with same color legend as these of the small circles).\nFigure 3: Again, the number of polygons in brackets is not informative. Include area extent instead the number of polygons. Regarding the volume rate of change, the error must be also included. Disconnection and separation circles are barely seen in the figures. I encourage manuscript authors to change symbols and colors for a proper identification.\nFigure 4: Panel letters are not included (a), b), c)…) what makes difficult figure interpretation. Figure 4 caption: Summer air temperatures are computed for which time period?, please, include it.\nFor panel b) is stated “The timing of an early, mid and late LIA is shown.” Where is this shown? I think this is an interpretation of the graph, at least nothing is marked in the graphs. If this sentence is an interpretation I would include it in the text not in the figure caption.\nFor the bottom left panel, please include a legend for the red line and the green triangle (mean and median I guess).\nJuneau airport data are used in the manuscript but it is not located in any figure or the distance to the study area is included. Please, provide, at least the distance from the airport to the ice field and the elevation of the weather station in the text/caption of this figure.\nExtended data 1: Please include the dates for which the ELA was computed. Hypsometric indexes are used in the manuscript? I guess these values are not used. Please, if so, remove them in this figure.\nSupplementary information:\nLine 57 and 60. Please include precipitation in mm (International system units)\nLine 76 There is a typo in the year. I guess it is not 2948.\nLine 156 : Where is Taku Point?? Show it or include a refence here.\nSupplementary figure 1: 0-100m and 100-200 m, ice thicknesses representation can be improved with more classes (0-50, 50-100….). This would show the areas with extreme glacier retreat in the forthcoming years.\nSupplementary Figure 3: Include the location of the study area in NW North America (at least a point with the location or an small rectangle).\nSupplementary Figure 5: I might remove the 5 upper glaciers in this figure as it is difficult to identify the distinct lines. From my understanding is better to show the methodology in 2-3 glaciers with good examples than showing it with 10 examples barely observed.", "reviewer_2": "This paper reconstructs ice volume and areal changes for the Juneau Icefield (JIF), a plateau icefield with outlet glaciers, and adjacent glaciers from the 19th century to 2019, using multiple data sets, including geomorphological and topographic maps, air photos and satellite imagery. The most noteworthy finding is that the rate of ice loss from the JIF has not been linear and has accelerated sharply since 2005 and that the acceleration is due, not to an acceleration in temperature increases, but is due to the JIF hypsometry, as the ELA for the icefield has now risen to the plateau surface. This means that ablation zone is now rapidly expanding over a wide area and the rate of ice loss will continue to increase even if temperature increases stay low.\nThese findings are of significance to glaciologists and climate and sea-level modellers for a number of reasons. Most importantly, the study suggests that current model projections of ice loss and resulting sea level rise from Alaskan and other ice fields may be underestimates, as they are based on a linear rate of ice loss. This has important implications for the rate of sea-level change in the next 100 years, as Alaskan glaciers are an important contributor in the short term.\nSecondly, the paper provides a model for the mechanisms of ice field retreat which is applicable to both modern and past ice fields globally. The paper builds on previous studies of Juneau ice field and is a significant new contribution to the literature. The use of sequential DEMs to assess ice loss is not new, (see, for example Berthier et al., 2010 and 2018), but this approach has been previously applied to changes in this century only. The generation of a sequential DEMs covering the period from 1770-2020 is an important achievement, and is a result of combining multiple datasets, including geomorphological and topographical maps, air photos and satellite imagery, to produce the time series.\nThe methodology used here is sound and the resulting data is of high quality; possible errors are discussed in detail in the methods and supplementary data expands on some of these data sets. There is sufficient detail provided to reproduce the work. However, I do have some issues with the conclusions drawn from the work presented. This is due more to the writing style and presentation of data within the main text than to any inaccuracies in approach. Nevertheless, some rearrangement/minor rewriting is necessary to highlight particular aspects of the data and make it clear when the authors are referring to their own or other datasets. In particular, I think the following needs addressing:\n1. A key finding of the paper is that the rising ELA of the Juneau Icefield now intersects the elevations of the icefield plateau; this is stated without preamble on line 171. However, the only clear ELA data presented in the results section is for the LIA. Instead, evidence is presented for thinning across the icefield between 2010-2020, plus a large increase in disconnections in this period. However, this data is buried within a paragraph (lines 114-124) and it took me some time to figure out that this must be the supporting evidence for the statement on line 171, and that this is clearly different to previous periods. The authors should clarify if this is the case. If it is not, then they need to make clear what the evidence is for the ELA rising to the plateau height.\n2. Line 178 – It is not clear here whether the dramatic reduction in icefield accumulation area refers to your data or the predicted changes. I assumed on first reading it was the predicted change, as the authors have not calculated a change in accumulation area volume, so cannot say it is dramatic as yet? If this does refer to observations presented in the paper, the authors need to make it clear clear why specifically theyu think the accumulation area has decreased dramatically, as this would require retreat of the ELA beyond the edge and toward the interior of the plateau.\n3. Lines 178- 196, I think, refer to the model of ice recession presented in Fig. 6, but the absence of an initial linking statement to this effect means that it reads as if the work presented here directly measure the changes implied, e.g. snow cover duration. This fails to distinguish between the authors’ own work and that of others adequately and needs addressing.\n4. I was surprised that there is no reference to the two papers by Berthier et al. (2010 and 2018) on estimating recent rates of change for Alaskan including JIF (2010) and JIF specifically (2018). The 2018 paper is mentioned in passing, but both papers show an earlier attempt to estimate ice loss using DEMs, although at a lower resolution, with much less coverage and for a much shorter time. In particular, it is worth considering a comparison of the authors’ work to the 2018 Dh/Dt work on specific glaciers, to see if there is a significant change. This need not be included in the main paper, but should be in the extended methods?\nThere are other minor issues/typos within the figure and text to address:\nLine 58: “way a..” – remove a\nFigure 1: I didn’t find this particularly helpful, other than the map showing the location of the JIF. Consider merging 1A with Figure 3 (but see below re Figs 2-3) and moving the rest to the supplementary/extended sections.\nLine 83: I would extend this to say why Lemon Creek is used specifically – presumably because of the excellent data set? But not everyone knows that as a matter of course.\nLine 88 – worth saying “underwent dramatic changes *in area* from the LIA..” – I had to reread the paragraph to figure that out. I recommend checking throughout that you have been precise in specifying area/volume etc. – it’s quite confusing otherwise.\nFigure 2A: the ELA circles almost entirely obscure the smaller glaciers so it is hard to see the LIA extent vs the 2019 extent. You could make the circles smaller, but a better alternative might be to change the order of the figures – put Figure 3 before Figure 2 as it shows the areal extent clearly. There is also room to add a v small map showing the location of the JIF on Fig. 3.\nFigure 2b: The figure caption vs the legend is confusing. The legend implies all the outlet glaciers except for Taku have increased in height, but I think the inverse applies? i.e. the figure caption is correct and this shows the amount of lowering, rather than the amount of change.\nLine 100 – also refer to your Figure 3.\nLine 107: grammar\nLine 115: grammar – should be ”…areal loss accelerated sharply…”\nLine 121: you need to distinguish more clearly between the “icefield plateau” (above) and the icefield (this line) – suggest “the whole icefield”. I also got confused in places between glaciers (presumably including outlets) and glaciers (not including outlets possibly?) – maybe worth saying “all glaciers” when you mean the lot.\nLines 119-122: key points buried in the middle of a paragraph – need to be highlighted….\nLine 125: this reads as a contradictory statement – i.e. “rates of areal recession rise steadily from 1948-2005” compared to lines 107 and 128, where you mention the increase in rate from 1979-2010. Needs clarifying.\nLine 130: - do you mean ice loss? Or have glaciers completely disappeared? – if the latter, more info needed!\nFigure 3: Love this, but the disconnection and separation dots are a little hard to see – brighter colours needed?\nFigure 4: It was a pain relating the captions to the graph. I suggest you rearrange the caption and lettering, going across instead of down, i.e.\na b\nc d\ne f etc.\nLine 155: I strongly recommend making your overall point at the end of the first sentence, i.e. temperatures important, but the relationship between temperature and volume/area loss is non-linear. Otherwise this is a very disjointed read. You can still make it at the end and suggest a threshold response.\nLine 171: see my point above – you haven’t drawn the links back to the data.\nLine 176: add “… observational data *presented here*…” – all very abstracted otherwise.\nLine 178: see main comments above\nFigure 6: This is ok, but the photograph caption is not accurate – that’s not the JIF, just a small area of it.\nLine 204: refer to supplementary evidence and reference 17.\nLine 234: grammar…\nLine 274: I think it would be better if 17 was not in superscript – it forms part of the main sentence, it’s not the equivalent of an asterix. The same replies to references 46 and17 below, and references in the extended data where the reference number is part of the main sentence.\nLine 466: you mention that 1979 was in 4 blocks below – include here with 1948?\nLines 587: how many are there?\nLine 590: This could be compared to the Berthier et al. (2018) example to see if your DEMs make a difference in calculating area/volume?", "reviewer_3": "The authors present a compilation of observations and reconstructions of changes of the Juneau Icefield that is probably unique for a glacierized mountain chain of this scale in its spatial and temporal scope, and in its comprehensiveness. Similar densities of observations are found, if at all, in regions like the Alps or Scandinavia, which arguably are less relevant than Alaska to many glaciological research questions. These other data sets have not been created because of their urgency, but because of convenience and traditions. The authors therefore deserve a lot of praise for their effort to improve the knowledge in a region that is much more relevant and much harder to investigate, and they have obviously succeeded in connecting the different qualities, quantities and characteristics of their data to build a panoramic view of the evolution of the Juneau Icefield since the Little Ice Age.\nThere are therefore very strong arguments that this work should be published; it would be of high value to the glaciological research community (and probably beyond). However, I see two fundamental problems in the manuscript that are not related to the data, but to their interpretation and the framing of the entire manuscript.\n(1) The smaller of these problems concerns the discussion of the relevance of the findings with respect to glacier models and projections of glacier evolution. The manuscript claims that “Current analyses of future glacier behaviour predict a linear future contribution of global (8) and Alaskan ice masses to sea level rise (5,6,8), and a linear rate of volume loss for Juneau Icefield (9).”, and that “the simplified physics deployed by many of the models in Glacier MIP […] accounts poorly for variation across different hypsometries (11,14)”. The manuscript argues that the nonlinearities supposedly found in the data (see below for more on this) cannot be captured by the current generation of models. I strongly disagree with this assessment. The given references (5, 6, 8, 9) show nonlinearities in the mass change rates over the 21st century, such that invoking “such linear projections” does not adequately reflect their characteristics. It is true that regional and global scale projections “hide” nonlinear behavior of individual glaciers in the aggregated numbers typically shown. This does not mean, however, that this kind of behavior is not represented at the scale of individual glaciers (and it even shows in the aggregate, as stated above). In particular, the well-known feedbacks between glacier hypsometry/geometry and mass balance referred to in the manuscript are well represented in state-of-the-art global glacier models as used in references 6, 8, and 9. Some other nonlinearities are possible in some, but not all of the global-scale models (e.g., concerning disconnections). There are additional nonlinearities which are not captured well by any of this type of model (the manuscript mentions some), these are addressed, e.g., in ref. 14. But those nonlinearities are not the ones the manuscript claims to be affecting the Juneau Icefield.\nThis issue does not affect the center of the manuscript, and because of this, I see this as the “smaller” of the two fundamental problems. I see it as fundamental since it affects a main motivation of the manuscript, being prominently given in the abstract (and repeatedly mentioned in the manuscript). It should be possible to reframe the manuscript so that it is more adequately referring to the state-of-the-art in glacier modeling.\n(2) The bigger of the two problems concerns the claim that nonlinear processes are driving the mass loss acceleration of the Juneau Icefield (see title).\nThe manuscript has not convinced me that there is in fact a nonlinear response of the icefield observable in the data. Essentially, this conclusion would have to rest on the scatter plot in Fig. 4j. Here, a nonlinear response of the rates of area and thickness change to temperature change would have to be visible. With only 4 (thickness) and 6 (area) point available, and quite a bit of scatter (which might be explained by precipitation anomalies), the argument for nonlinearity is weak. The statement in lines 165 – 168 is therefore not justified; neither the part concerning the present manuscript, nor concerning the modeling studies (see above). Arguing the other way around (“nonlinear responses found in projections are not yet robustly detected in observations/reconstructions”) would be a similar misstatement, but arguably not worse than what is said in the manuscript. As an alternative to showing a nonlinear response of the entire icefield (Fig. 4j), it would be possible to demonstrate that nonlinear effects are affecting the glaciers by showing that the ELA has crossed into the plateau region (i.e., by looking at ELA changes and how they relate to the hypsometry of the glacier), or by investigating the mass balance profile and ice flow velocities in the region of glacier disconnections (see explanation of feedbacks around Fig. 6). Unfortunately, this is not done in the manuscript.\nIn this context, the section “Hypsometric controls of icefield behavior” is very problematic: the relevant mechanisms are well-known, but they are also explained (and referenced) well in the manuscript. However, the manuscript does not show that a single one of them is actually affecting the Juneau Icefield in their data. Remarkably, they use an edited photograph (Fig. 6) for illustrating “how it might look in a few decades”. That is exactly the point: in this manuscript, I don’t find any evidence that these feedbacks are active in present day. However, the manuscript very boldly claims that it is the case in the title, the abstract, scattered in the text, and in the final, conclusive paragraph.\nThe problem of a missing demonstration of nonlinearity in the data may arise from a confusion of “nonlinearity” and “acceleration”. E.g., an acceleration in mass change can be a linear response to an increased temperature anomaly. The feedbacks explained in the second and third paragraph of the manuscript (and Fig. 6) would lead to an acceleration in mass change without additional changes in the temperature anomaly. Therefore, I would suggest that the authors carefully decide whether they actually want to claim (and then have to demonstrate) a nonlinear response of the Juneau Icefield to climate change (which, with the data on hand, seems to be very difficult), or whether they wish to make a case for an impressive acceleration of glacier change in Alaska in the recent decades, which would be a very strong case.\nEither way, very substantial work on the text (and potentially, the data analysis) would be needed, so that I would probably consider the manuscript as “new”. This would imply a suggestion to reject the manuscript. However, as I am convinced that the work definitely should be published eventually, I would also agree with considering this as a case of “major revisions”. Expecting the manuscript to go through another round of revisions (potentially in a different format, pending the editor’s decision), I have refrained from specific/minor comments on all parts of the manuscript except the main text, where I hope also these specific comments will be helpful (also, I don’t see any major issues in the methods).\nSpecific/minor suggestions for the main text:\nL18: „The majority…“ would imply > 50 %, which is not the case according to Hugonnet et al. (https://doi.org/10.1038/s41586-021-03436-z); perhaps change to: “The largest regional contribution…”\nL19: “…projected to become even more dominant…” I’m not sure I would agree, and would expect that this depends on the scenario: in a high-emission scenario, e.g., I would expect that regions that have more ice left (peripheries of Greenland and Antarctica, the Canadian Arctic) become more dominant over time.\nL22: “…with…” since this refers to the acceleration (but no acceleration value is given), perhaps change to: “…leading to…”\nL35: suggest to change to: “Under the higher-emission scenario RCP 4.5…”\nFig. 1 is a bit disorganized: e.g., I can’t find any ice-dammed lakes, lakes, or sandurs on the map (or photographs), and referring to north- and south-flowing glaciers in panel d, without including a north arrow, is not ideal (perhaps label the glaciers in panel d directly?). I would also expect that the individual glaciers mentioned in the text (e.g., lines 83-87) appear in the map, or vice versa. Please include a reference to (e.g.) Fig. 2a in the relevant sections of the text.\nL82-90: please indicate the uncertainties of the assessed changes in ELA, volume, and thickness. These uncertainties are important for understanding the significance of the changes.\nFig. 2: I appreciate the difficulty to put the rich data set on the map(s), and unfortunately, I don’t have a good solution. But it is problematic that (i) the circles of the ELA estimates cover a large fraction of the trimline and moraine markers and (ii) moraine polygons are shown gray on gray background.\nL118: what do “M” and “SD” mean? What units are used here?\nL129: It has not been shown that the acceleration is really exponential; I would suggest to state this more carefully (also in the abstract and other places in the manuscript).\nFig. 3: The tables in the lower left should either really be formatted as tables (outside of the figure) or left out (they are already shown in Fig. 4). Regarding number of glaciers: how are separations treated? Do they increase the number of glaciers?\nFig. 4a: please explain red line here, not just in panel f. b: Where is the timing of an early, mid and late LIA shown? Last sentence of the caption: what does this refer to? Panel f does not show reduced precipitation during 2005 – 2015?\nFig. 5: legend and caption: please correct from dhdt to dh/dt. More generally, I would suggest to move this figure to the supplement; it is not that central to the arguments in the manuscript."}, "author_response": "Reviewer #1 (Remarks to the Author):\nReview of the manuscript “Accelerating glacier volume loss on Juneau Icefield driven by non-linear processes” by Bethan Davies et al.\nThe research leaded by B. Davies, shows an interesting case study on Juneau Icefield in Alaska about the evolution of the glaciated area and volume of the glaciers in this region. Manuscript authors apply state of the art remote sensing techniques to delimitate glacier extent in different time periods. They are also able to reconstruct glacier surfaces to determine volume differences across distinct temporal windows. The manuscript reads well and shows reliable results for distinct time periods that without any doubt show the fast decline of Juneau Icefield in the last decade compared to previous period with data.\nWe are delighted that the reviewer found the manuscript enjoyable and reliable. Thank you for the constructive and helpful comments.\nMajor comments\nDespite I have enjoyed the read of the manuscript I consider that there are some major points to be addressed before considering for publication:\n• The manuscript claims that non-lineal processes drive the volume loss on Juneau Icefield but this is not really shown in the manuscript. The database demonstrates that volume and area losses have accelerated in the last years and that many glaciers have disconnected or separated. On the contrary processes that control and will control glacier evolution are enumerated and related in a schematic (figure 6). These processes have not been measured or evaluated with the data processed and thus, from my understanding some sections of the manuscript must be changed considering this (including the title). See minor comment for lines 247-248.\nThere is a substantial data (McNeil et al., 2020; McNeil et al., 2019; O'Neel et al., 2019; Pelto, 2019) that shows the rising ELAs and transient snowlines in Juneau Icefield, which now reach the plateau. This was rather buried in the Supplementary Information, but has now been moved to the main text. We have emphasised this on figure 6i by showing the height of the plateau with respect to the height of the ELAs, Secondly, we now include an analysis of the late summer snowline on Juneau Icefield in 2019 AD from sentinel 2 imagery. This late summer snowline clearly intersects the plateau and demonstrates the significance of the icefield hypsometry in sensitivity to rising snowlines. We have noted this in the Results section:\n“This rapid glacier recession was accompanied by ELAs rising to the plateau height in the period 2005-2020(McNeil et al., 2020; McNeil et al., 2019) (Figure 6i), and exceptionally high snowlines on the plateau in 2018 and 2019(Pelto, 2019) (Figure 5). Late summer snowlines mapped in 2019 from Sentinel 2 imagery (Figure 5) had a mean altitude of 1600 m asl (SD 161 m; n=163). Snowlines on Taku Glacier were a mean altitude of 1444.7 m in 2019; this compares with a measured ELA in 2019 of 1528 m(McNeil et al., 2019). The snowline here is therefore well above the southern plateau at 1200 m asl. No snow was observed on Lemon Creek Glacier, in line with the measured ELA of 2023 m that year(McNeil et al., 2016). Snowlines on Meade, Llewellyn and Field Glacier had a mean height of 1671 m (SD 119 m), well above the height of the northern plateau at 1500 m asl (Figure 5).”\nThis rising ELA then results in a number of processes related to albedo, thinning on the plateau, shrinking accumulation area, disconnections and so on, which are enumerated in the Figure 6 (now Figure 8). We have edited the Discussion, highlighting the rise in the ELA and its intersection with the icefield plateau.\nFinally, in line with suggestions from Reviewer #3, we have reframed the abstract and discussion on conclusions around the accelerated mass change observed driven by the ELA and glacier thinning now intersecting with the icefield plateau, which drives a number of hypsometrically driven melt-accelerating feedbacks (e.g. related to albedo, thinning, accumulation area ration, fragmentation). Therefore we have strengthened the discussion with reference to these processes and hope that this is now addressed.\n• The volume and glacier extent for 1948 and for the LIA maximum, have some uncertainties not included in the analysis. For example the glacier extent in 1948 is based in a 66.6% in images of this year, but the glacier area not observed in this year is inferred from topographic maps, and satellite images from 1970 to 1982. Similarly, LIA extent.\nWe have clarified the data sources and analysis in the text and in the response to reviewers below. All uncertainties are carefully included in the manuscript. For 1948 in particular, careful analysis of the topographic maps used is included in the supplementary information, including the year used for the imagery. This date (of either the year of publication or the year if the imagery used) is included in the attribute information for this dataset – all of which is now available in GLIMS viewer (https://www.glims.org/maps/glims).\nThe date of the imagery used is used to calculate annual rates of recession. For those small number of glaciers (12.5% of glaciers) where the earliest available date is the 1980s Landsat imagery, this is included in the attribute information as well. Since glacier shrinkage has occurred consistently across the icefield over this time period, this means that the 1948 glacier area is a minimum. This is detailed in the Methods and in the Supplementary Information.\nSimilarly, the uncertainty in area extent at the LIA maximum is of same order of magnitude that the extent in 2005, 2015, 2019, periods for which high spatial resolution satellite observations were available. Despite the LIA extent can be well identified from methods described from lines 250 to 271 of the supplementary information, an evaluation of the uncertainty of this extent must be included.\nWe have conducted a thorough analysis of uncertainty in glacier area extent for each timeslice (supplementary table 6). The moraines, trimlines and areas of ice scour are easily identified and these data have been published in detail already (Davies et al., 2022). The use of Sentinel satellite imagery for these geomorphic landforms, as opposed to the lower resolution Landsat imagery available in the 1990, 2005 and 2015 timeslices, allows a high degree of confidence in the mapped geomorphology. These geomorphic details are briefly summarised in Supplementary figure 4.\nThe detailed geomorphological mapping allows glacier outlines at the LIA to be drawn with precision and the detailed assessment of uncertainty in Supplementary figure 6 follows best practice guidelines (Paul et al., 2017). There is limited scope for greatly differentiated LIA glacier extent because glacier tongues largely remain constrained in topographic valleys, with a very clear geomorphic imprint (and with few marine terminating glaciers). Therefore the uncertainty in glacier extent as presented in the supplementary information at the LIA is robust.\nSimilarly the ELA reconstruction for the LIA might have an error assessment.\nUncertainty in ELA reconstruction for palaeoglaciers derives from choice of balance ratio (Rea, 2009). Analysis by our team members on the influence of different balance ratios on glacier ELA shows that ELA is only weakly sensitive to choice of balance ratio (Carrivick et al., 2023). Added to manuscript:\n“Published work shows that the choice of balance ratio has little influence on glacier ELA (Carrivick et al., 2023).”\nIn detail, Carrivick et al. (2023) their supplementary information Figure SI_8 shows that varying BRs between 1.2 and 2.6 [noting Oien et al., (2022) recommend 1.56 for paleoglaciers and Rea et al., suggest 2.24 for arctic glaciers] produces mean ELAs varying between 1040 and 950 m. asl.\n• The readability of some figures must be improved. In some cases small circles avoid the identification of glacier extent in particular time periods. Similarly, some arguments are directly related with figures panels that do not have a correct identification doing difficult to follow the discussion. See minor comments for the figures.\nFigures have been improved as suggested. Panels are all labelled.\nMinor comments:\nWe have made all of these amendments as suggested, which has improved the clarity of the manuscript. Thank you for these suggestions.\nLine 18-19 an 34: In the second sentence of the abstract it is stated “The majority of the ice volume loss comes from Alaskan icefields,…”. I see a bit biased this previous sentence when in line 34-35 it is stated that ” with 13 cm coming from glaciers, a 25% of which will come from Alaska”. Definitely it is not the majority of volume loss. Change conveniently.\nAmended: “. The glaciers experiencing the most ice loss are the Alaskan icefields”\nLine 26-27: While it is true that that non-linear processes have to be incorporated into glacier evolution models; this is not a direct conclusion of the work. I would remove this last sentence of the abstract, and show a more general conclusion obtained on the frame of this study.\nThis is the important implication and wider significance of the work. It is one of our key findings and is supported by our data, as well as published mass balance observations across the icefield. Therefore, it should be included as it is the key finding. We have moved the glaciological and mass balance observations from the Supplementary Information to the main text, which now supports this argument more clearly.\nLine 64: I have found quite long the supplementary information section about LIA ice extent and the climatic conditions observed in this period. I would encourage to reduce the extent or even completely remove these sections and just include appropriate references to this period. Include reference to research articles focusing in this time period and remove the reference to supplementary information.\nWe argue that the climate and evidence for glacier extent in the LIA is a crucial part of the paper, and a comprehensive recent review of this is not available in the literature. This data is needed as it explains why the glaciers have been dated to a particular timeframe. There is some complexity in the behaviour of some glaciers, which requires detailed explanation. However, we have shortened and condensed this section of the Supplementary Information, and removed unnecessary detail.\nLine 84: Reconstruction of LIA icefield ablation areas: Provide a short explanation on how was obtained.\nThis is the area below the reconstructed ELA and is explained in line 559 in the Online Methods (“Reconstruction of ELAs and ablation areas).\nLine 84: Include uncertainty in the total ice volume estimation for LIA reconstruction “±XX km3”. Same for total cumulative ice loss (315.31km2±XX km2 ).\nAmended. The full details of uncertainties are additionally provided in extended data 4.\nLine 110: The extent of 17.54km2 was an area decrease, right? It is not clear as in this sentence, it has not been computed as an area loss. In the next sentence it is clearly identified the 3.72 km as a loss. Please, along the manuscript, state clearly when loses are observed, eider with a “-” before the number, or clearing stating that a “decrease” or a “loss” has been observed.\nWe try to be clear in the manuscript by being consistent in use of signs. Talking about glacier shrinkage and then giving a “-“ sign can be confusing and a double negative. This sentence is talking about rates of glacier area loss. We have been clear throughout the manuscript and careful to note that these are rates of glacier area or volume loss.\nLine 118: M and SD are?? Mean and standard recession rates? Per year? Maybe I have missed the introduction of these acronyms, but I am not able to find them. Please, include them conveniently. Same for “p” and “t”\nM is Mean, SD is Standard Deviation. They are commonly used terms. Defined now in manuscript:\n“Rates of recession were 1.5 times faster from 2015-2019 (Mean (M) =5.88, Standard Deviation (SD)=5.26)”\nt is the T statistic, p is the p value. These are given in the usual APA statistics format.\nLine 120: A volume loss of -5.91km3 might be volume gain. Please see comment for line 110 and change conveniently. Amended as asked\nLine 126: Rates of glacier shrinkage of area or volume? Please specify. Amended.\nLine 129: These doubling rates of volume loss after 2010 can be observed in Figure 4? I think they are, please, include which panel in figure 4 show it. Amended.\nLines153-168: The discussion is in some sections difficult to follow as figure 4 panels are not conveniently included. Amended as asked\nLine 180: Also include here this reference: https://doi.org/10.1038/s41558-018-0296-5. Added as asked.\nLine 204: The debris-covered areas isolation from the climatic evolution might be very important. I think this can be longer discussed in these lines. Moreover, as glacier continue their area decrease, debris-covered areas will increase thus their impact will be more evident. Noted in the manuscript.\nLine 224-227: Include references to state where and which are the ice fields in Canada, Greenland, …\nThere isn’t a single reference that compiles these. These are icefields named in maps largely. We have added references where available or relevant.\nLine 247-248: It is stated that “as our observation demonstrate, once the ELA rises above the plateau non-linear processes operate”. However the only reference to the ELA and the non linear processes in the manuscript is this: “The critical driver of this non-linear and accelerating glacier decline at Juneau Icefield is that the rising ELA now intersects the elevations of the icefield plateau (Figure 6).” Figure 6 is an schematic figure, showing processes, but it does not demonstrate with observations that non-lineal processes operate.\nAmended to “once ELAs rise and thinning occurs on the plateau”. However, we have ELA reconstructions from published work (McNeil et al., 2020; O'Neel et al., 2019; Pelto, 2019). These ELAs are shown in figure 4i, and the height of the plateau (seen in figure 6) is now described in the main text. The ELAs and glaciology were previously outlined in detail in the SI. We have moved this information to the main manuscript, which aids clarity considerably, as it is much clearer that the rising ELAs are driving the thinning observed on the plateau. Secondly, we have added data on transient snowlines from 2019, following on from published work (McNeil et al., 2020; Pelto, 2019), showing that the snowlines are reaching the plateau now.\nThe schematic Figure 6 (now Figure 8 due to restructuring) illustrates the processes that are driving the behaviour observed on Juneau Icefield, and is supported by reference to other published work and observation of the rising ELAs.\nFrom my understanding the main message of this manuscript (including the title) cannot be sustained with this figure, and thus must be conveniently changed.\nWe have updated the main message of the manuscript and changed the title in line with this.\nMoreover, the ELAs included in the extended data 1, do not include any temporal information. These are recent ELAs?? Please, clarify the ELA’s date computation in extended data 1.\nThis figure (now Figure 1) provides the most recent available mean ELAs – references are provided in the figure caption. The dates for the ELAs are added to the figure.\nOnline methods:\nLine 423: Cloud cover patches are limited in area, but which %?\nOverlapping imagery is used to ensure that cloud cover does not impede mapping. Images listed in the supplementary information.\nLine 425-436: The 28% of the glacier extent in 1948 is mapped from the topographic maps from which year?, the 3.7% from 1979 AHAP, the 1.67% from 1982 Landsat… (all data from Supplementary table 2). For me, this is not the true glacier extent in 1948 (more than 30 years in between some acquisitions). This must be clearly stated in the manuscript, and the glacier areas computed for 1948, include a higher uncertainty showing the not negligible % of the study area from a different (and later) observations. A good example where this must be shown is in figure 3 glacier areas.\nThis is shown clearly in Supplementary Figure 6. As is shown in this figure, the 1948 aerial photo mosaic covers most of the glacier termini in the west, north and south of the icefield. With the exception of Llewellwyn Glacier, the icefield outlet glaciers are well covered, and it is the small peripheral glaciers to the east of the icefield that mostly lack coverage from aerial photos. The outlet glaciers not covered by the photo mosaic are covered by the topographic maps. These smaller glaciers mapped from Landsat imagery account for only a very small proportion of glacier area, and therefore uncertainties here are unlikely to affect the conclusions of the paper. For the photo mosaic, larger patches not covered remain over the accumulation area. Where there is coverage there is little difference observed in the accumulation area, and uncertainties or differences here are unlikely to affect the conclusions of the manuscript.\nA good approach to determine if images from 1979, 1981 etc (or the topographic maps) can be exploited to identify glacier limits in 1948, is to analyze in areas with data in 1948 and later years, which is the deviations between these extents.\nFirstly, we note in the manuscript that the date of image acquisition is provided in the attribute information for each glacier, and the rate of recession provided is annualised using this date.\nSecondly, we note in the manuscript that the area for 1948 is a minimum, as a small amount of glacier area is calculated using later imagery. Therefore the area in 1948 is the best available estimate, and is clearly acknowledged to be a minimum. We have added text to the Methods emphasising this.\nThirdly, for the 117 glaciers mapped from topographic maps (supplementary table 3), detail on the date of the topographic maps is provided in Supplementary Table 5; this includes the date of the map published and the number of glaciers used. The published date of the map is used as the date for the glacier attribute information if clear information on source data is not provided.\nAdded to the Methods in the manuscript:\n“This means that the glacier extent in 1948 is the best available estimate given the data availability. Annual rates of recession are calculated using the date of the imagery or data source used to map glacier extent. The glacier extent in 1948 is considered a minimum; use of later datasets results in some glaciers being underestimated due to continued glacier recession.”\nLine 442 GDEM data is ASTER GDEM?, please clarify. Yes, amended.\nLine 502: I would only talk about LIA maximum without including a specific year (i.e. 1770) For the maximum LIA extent I do now see relevant + or – 20 years. Same for supplementary table 6. Amended\nLine 530: And what did authors do with images from 1980, 1981 and 1982 (supplementary table 2), which are also used for computing 1948 glacier extent?\nAs stated, if the same imagery was used, then the glaciers are removed from analysis of rates of recession.\nLine 579-580: From previous sentence this is obvious, remove the sentence: “In effect, glaciers…”\nYes, this may be obvious, but we retain for clarity.\nLine 583-584: “Then we” appears twice. Amended\nLine 587: Which is the % of glaciers with no DEM coverage?, please include this information here.\nThis data has been added to the Supplementary Information (Supplementary Table 8). Added to the text:\n“This amounts to 44.6% of glacier area for LIA-1948, 16.8% for 1948-1979,4.4% for 1979-2000, and 0% for 2000-2020 (Supplementary table 8).”\nLine 590: Mean dhdt should be dh/dt?? Same in line 609 Amended throughout\nFigures:\nFigure 1: I encourage including here a map/satellite image showing the entire study area as this is the first figure of the manuscript and thus readers want to identify the main characteristics of the study site. I would remove panel D, which is not informative at this initial point of the manuscript. From this “complete” overview of Juneau Icefield I would delimitate with a rectangle where panels C and E are located.\nThis figure is useful as it highlights the different data sources used in the publication. However, we have edited it as requested, enlarging Panel A to show in more detail the icefield, and including the coastline and ocean more clearly.\nPanel A: Include boundaries for a better sea-land limit identification. Amended as requested.\nPanel E. Some contour lines (with transparency) might help to interpret the map. Amended as requested.\nGlacier limits in 1931 and 1908-1910 are informative, but if I am right are not used in the manuscript, and thus are not needed.\nThey are used to test and evidence continued glacier recession from LIA-1948. This is discussed in detail in the supplementary results.\nDisconnection separation marks are barely identified in panel E. I encourage to include here not only circles with distinct colors and use triangles, squares etc…for easier identification.\nThe symbology here is used to be consistent with later figures. The symbology is clearly discernible on my machine and printed version. Darker colours have been used to aid clarity.\nFigure 2: The number of polygons (in brackets), is not informative of the number of glaciers, I would remove this information and change it by the total extent of the ice-dammed lakes, alluvial fin,…\nThe number of glaciers is important; it highlights graphically the number of glaciers disappearing between the two observational periods. We have added the area of key landforms (glaciers, lakes, fans) as requested.\nELA (LIA) circles hamper the identification of the estimated LIA glacier extent. For me is more informative the LIA glacier extent (dark right) than the ELA estimation. I would remove all circles in A and include a “C pannel”( or include a new figure in the supplementary material), just including ELA (LIA) with features of different colors (with same color legend as these of the small circles). We have removed the ELAs from the figure in the main text and added them to a figure in the SI, and referred back to this.\nFigure 3: Again, the number of polygons in brackets is not informative. Include area extent instead the number of polygons. For clarity, the area is given in the figure, in the bottom right. There isn’t enough space in the top left of the figure.\nRegarding the volume rate of change, the error must be also included. Note that this is area change, not volume change. Following best practice we do not provide an uncertainty here.\nDisconnection and separation circles are barely seen in the figures. I encourage manuscript authors to change symbols and colours for a proper identification. These are an important part of the story; in any case the disconnections are provided as shapefiles if needed for further clarity.\nFigure 4: Panel letters are not included (a), b), c)…) what makes difficult figure interpretation. Panel letters were included in the y axis for each graph. However we have added them as titles for each panel to aid clarity.\nFigure 4 caption: Summer air temperatures are computed for which time period?, please, include it. The time period is shown on the x axis of the graph.\nFor panel b) is stated “The timing of an early, mid and late LIA is shown.” Where is this shown? I think this is an interpretation of the graph, at least nothing is marked in the graphs. If this sentence is an interpretation I would include it in the text not in the figure caption. Apologies, this was an error – this refers to panel d. Amended.\nFor the bottom left panel, please include a legend for the red line and the green triangle (mean and median I guess). This is already in the caption for graph e. Clarified in the caption.\nJuneau airport data are used in the manuscript but it is not located in any figure or the distance to the study area is included. Please, provide, at least the distance from the airport to the ice field and the elevation of the weather station in the text/caption of this figure. Juneau Airport is added to the updated Figure 5.\nExtended data 1: Please include the dates for which the ELA was computed. These data have been added to the figure. Hypsometric indexes are used in the manuscript? I guess these values are not used. Please, if so, remove them in this figure. The hypsometric index is referred to in the discussion (bottom heavy versus top heavy icefields).\nSupplementary information:\nLine 57 and 60. Please include precipitation in mm (International system units). Amended.\nLine 76 There is a typo in the year. I guess it is not 2948. Amended\nLine 156 : Where is Taku Point?? Show it or include a refence here. We have added a new figure with these placenames all clearly marked.\nSupplementary figure 1: 0-100m and 100-200 m, ice thicknesses representation can be improved with more classes (0-50, 50-100….). This would show the areas with extreme glacier retreat in the forthcoming years.\nAmended as requested. Symbology changed to show this more clearly.\nSupplementary Figure 3: Include the location of the study area in NW North America (at least a point with the location or an small rectangle). Amended as requested.\nSupplementary Figure 5: I might remove the 5 upper glaciers in this figure as it is difficult to identify the distinct lines. From my understanding is better to show the methodology in 2-3 glaciers with good examples than showing it with 10 examples barely observed. This figure clearly shows this in the larger Willison, Echo and other glaciers; all glaciers are included to show the diversity of glacier sizes and locations used.\nReviewer #2 (Remarks to the Author):\nThis paper reconstructs ice volume and areal changes for the Juneau Icefield (JIF), a plateau icefield with outlet glaciers, and adjacent glaciers from the 19th century to 2019, using multiple data sets, including geomorphological and topographic maps, air photos and satellite imagery. The most noteworthy finding is that the rate of ice loss from the JIF has not been linear and has accelerated sharply since 2005 and that the acceleration is due, not to an acceleration in temperature increases, but is due to the JIF hypsometry, as the ELA for the icefield has now risen to the plateau surface. This means that ablation zone is now rapidly expanding over a wide area and the rate of ice loss will continue to increase even if temperature increases stay low.\nThese findings are of significance to glaciologists and climate and sea-level modellers for a number of reasons. Most importantly, the study suggests that current model projections of ice loss and resulting sea level rise from Alaskan and other ice fields may be underestimates, as they are based on a linear rate of ice loss. This has important implications for the rate of sea-level change in the next 100 years, as Alaskan glaciers are an important contributor in the short term.\nWe are delighted by these positive comments and that the reviewer recognises the significance of the work.\nSecondly, the paper provides a model for the mechanisms of ice field retreat which is applicable to both modern and past ice fields globally. The paper builds on previous studies of Juneau ice field and is a significant new contribution to the literature.\nThank you for these generous comments.\nThe use of sequential DEMs to assess ice loss is not new, (see, for example Berthier et al., 2010 and 2018), but this approach has been previously applied to changes in this century only. The generation of a sequential DEMs covering the period from 1770-2020 is an important achievement, and is a result of combining multiple datasets, including geomorphological and topographical maps, air photos and satellite imagery, to produce the time series.\nWe are delighted that the reviewer recognises achievement here. Thank you for these generous comments.\nThe methodology used here is sound and the resulting data is of high quality; possible errors are discussed in detail in the methods and supplementary data expands on some of these data sets.\nWe are delighted that the reviewer recognises our method as sound and the data as high quality.\nThere is sufficient detail provided to reproduce the work.\nThank you for this comment.\nHowever, I do have some issues with the conclusions drawn from the work presented. This is due more to the writing style and presentation of data within the main text than to any inaccuracies in approach. Nevertheless, some rearrangement/minor rewriting is necessary to highlight particular aspects of the data and make it clear when the authors are referring to their own or other datasets. In particular, I think the following needs addressing:\n1. A key finding of the paper is that the rising ELA of the Juneau Icefield now intersects the elevations of the icefield plateau; this is stated without preamble on line 171. However, the only clear ELA data presented in the results section is for the LIA. Instead, evidence is presented for thinning across the icefield between 2010-2020, plus a large increase in disconnections in this period. However, this data is buried within a paragraph (lines 114-124) and it took me some time to figure out that this must be the supporting evidence for the statement on line 171, and that this is clearly different to previous periods. The authors should clarify if this is the case. If it is not, then they need to make clear what the evidence is for the ELA rising to the plateau height.\nThe evidence for the ELA rise is from the detailed literature and long mass balance records from the Icefield (McNeil et al., 2020; McNeil et al., 2019; O'Neel et al., 2019). Cumulative mass balance data is shown in Supplementary figure 2a and rising ELAs can be seen in Figure 6i. These details were reviewed in the Supplementary information, but they are an important part of the story and key evidence. We have therefore moved these details, with data on the rising ELAs and comparison with the plateau height, to the introduction of the manuscript. We have also added the mapped snowline in 2019, with comparison to the plateau, in figure 5. We have added further detail on the newly mapped snowlines to the Supplementary Information. This should clarify the ELA rise and support the interpretations more clearly.\n2. Line 178 – It is not clear here whether the dramatic reduction in icefield accumulation area refers to your data or the predicted changes. I assumed on first reading it was the predicted change, as the authors have not calculated a change in accumulation area volume, so cannot say it is dramatic as yet? If this does refer to observations presented in the paper, the authors need to make it clear why specifically they think the accumulation area has decreased dramatically, as this would require retreat of the ELA beyond the edge and toward the interior of the plateau.\nThis is from the mass balance modelling in the paragraph above. Citation added for clarity.\n3. Lines 178- 196, I think, refer to the model of ice recession presented in Fig. 6, but the absence of an initial linking statement to this effect means that it reads as if the work presented here directly measure the changes implied, e.g. snow cover duration. This fails to distinguish between the authors’ own work and that of others adequately and needs addressing.\nWe have added a link to the figure to aid clarity.\n4. I was surprised that there is no reference to the two papers by Berthier et al. (2010 and 2018) on estimating recent rates of change for Alaskan including JIF (2010) and JIF specifically (2018). The 2018 paper is mentioned in passing, but both papers show an earlier attempt to estimate ice loss using DEMs, although at a lower resolution, with much less coverage and for a much shorter time. In particular, it is worth considering a comparison of the authors’ work to the 2018 Dh/Dt work on specific glaciers, to see if there is a significant change. This need not be included in the main paper, but should be in the extended methods?\nThese are excellent papers, and the lack of referencing was only in an attempt of brevity. Direct comparison with Berthier et al 2010 is challenging because this work only covered a small fraction of Juneau Icefield, and because the data are not publicly available for interrogation.\nThe comparison would rather be between the Hugonnet data (2000-2019) and the Berthier results (2000-2016), which is more of a comparison between the gaussian process registration used in Hugonnet et al and the linear regression used in Berthier et al. A different comparison would be with the 2010 paper that used DEM differencing, but challenging without having access to the original data. It is also not likely to be particularly illuminating due to difference in the methods and timeframes used.\nWe have added a reference to both papers in the ‘Glacier mass balance’ section at the paper’s introduction.\nThere are other minor issues/typos within the figure and text to address:\nWe have made all of these amendments as suggested, which has improved the clarity of the manuscript. Thank you for these suggestions.\nLine 58: “way a..” – remove a Amended\nFigure 1: I didn’t find this particularly helpful, other than the map showing the location of the JIF. Consider merging 1A with Figure 3 (but see below re Figs 2-3) and moving the rest to the supplementary/extended sections. Figure 1 has been updated.\nLine 83: I would extend this to say why Lemon Creek is used specifically – presumably because of the excellent data set? But not everyone knows that as a matter of course. Because ELAs have been measured here (see the text moved to the introduction).\nLine 88 – worth saying “underwent dramatic changes *in area* from the LIA..” – I had to reread the paragraph to figure that out. I recommend checking throughout that you have been precise in specifying area/volume etc. – it’s quite confusing otherwise. Clarified in text.\nFigure 2A: the ELA circles almost entirely obscure the smaller glaciers so it is hard to see the LIA extent vs the 2019 extent. You could make the circles smaller, but a better alternative might be to change the order of the figures – put Figure 3 before Figure 2 as it shows the areal extent clearly. There is also room to add a v small map showing the location of the JIF on Fig. 3. This figure (now figure 4) has been updated and amended, so it is clearer. The reconstructed LIA ELAs have been moved to the supplementary information. The location is shown clearly in figures 1 and 2.\nFigure 2b: The figure caption vs the legend is confusing. The legend implies all the outlet glaciers except for Taku have increased in height, but I think the inverse applies? i.e. the figure caption is correct and this shows the amount of lowering, rather than the amount of change. figure amended as requested.\nLine 100 – also refer to your Figure 3. Amended as requested\nLine 107: grammar amended\nLine 115: grammar – should be ”…areal loss accelerated sharply…” Amended as requested\nLine 121: you need to distinguish more clearly between the “icefield plateau” (above) and the icefield (this line) – suggest “the whole icefield”. I also got confused in places between glaciers (presumably including outlets) and glaciers (not including outlets possibly?) – maybe worth saying “all glaciers” when you mean the lot. Amended as requested\nLines 119-122: key points buried in the middle of a paragraph – need to be highlighted….Have inserted paragraph break so these key points are emphasised in a new paragraph.\nLine 125: this reads as a contradictory statement – i.e. “rates of areal recession rise steadily from 1948-2005” compared to lines 107 and 128, where you mention the increase in rate from 1979-2010. Needs clarifying. Clarified in text.\nLine 130: - do you mean ice loss? Or have glaciers completely disappeared? – if the latter, more info needed! 108 glaciers have disappeared; but this is referring to ice mass loss rather than loss of individual glaciers. Clarified in the text.\nFigure 3: Love this, but the disconnection and separation dots are a little hard to see – brighter colours needed? Amended as suggested.\nFigure 4: It was a pain relating the captions to the graph. I suggest you rearrange the caption and lettering, going across instead of down, i.e. We go in this order so that the time lines are the same on the up-down axes. Left column is 1650-2020 AD, right column is 1940-2020 AD. This makes the datasets easier to compare. We have added titles to the panels (a, b, c, d etc), which makes it easier to see and identify the panels.\na b\nc d\ne f etc.\nLine 155: I strongly recommend making your overall point at the end of the first sentence, i.e. temperatures important, but the relationship between temperature and volume/area loss is non-linear. Otherwise this is a very disjointed read. You can still make it at the end and suggest a threshold response. We have re-ordered the paragraph as suggested.\nLine 171: see my point above – you haven’t drawn the links back to the data. Amended\nLine 176: add “… observational data *presented here*…” – all very abstracted otherwise. Amended as suggested\nLine 178: see main comments above Amended as suggested\nFigure 6: This is ok, but the photograph caption is not accurate – that’s not the JIF, just a small area of it. Noted. This is Vaughan Lewis Icefall. Noted in the caption.\nLine 204: refer to supplementary evidence and reference 17. Reference to Davies et al 2022 included.\nLine 234: grammar… Amended\nLine 274: I think it would be better if 17 was not in superscript – it forms part of the main sentence, it’s not the equivalent of an asterix. The same replies to references 46 and17 below, and references in the extended data where the reference number is part of the main sentence. The citation is included now in the sentence.\nLine 466: you mention that 1979 was in 4 blocks below – include here with 1948? 1979 was in one block and 1948 in 4 blocks, which is outlined in the manuscript.\nLines 587: how many are there? This information has been added in the Supplementary Information (supplementary table 8) and the text.\nLine 590: This could be compared to the Berthier et al. (2018) example to see if your DEMs make a difference in calculating area/volume? As noted above, any difference would be down to the difference in regression approach, which isn’t the point of the study undertaken here.\nReviewer #3 (Remarks to the Author):\nThe authors present a compilation of observations and reconstructions of changes of the Juneau Icefield that is probably unique for a glacierized mountain chain of this scale in its spatial and temporal scope, and in its comprehensiveness.\nThank you for this positive and supportive comment.\nSimilar densities of observations are found, if at all, in regions like the Alps or Scandinavia, which arguably are less relevant than Alaska to many glaciological research questions. These other data sets have not been created because of their urgency, but because of convenience and traditions. The authors therefore deserve a lot of praise for their effort to improve the knowledge in a region that is much more relevant and much harder to investigate, and they have obviously succeeded in connecting the different qualities, quantities and characteristics of their data to build a panoramic view of the evolution of the Juneau Icefield since the Little Ice Age.\nThank you for these very positive and generous comments and we are delighted that the reviewer acknowledges the value of the dataset.\nThere are therefore very strong arguments that this work should be published; it would be of high value to the glaciological research community (and probably beyond).\nThank you for this positive and generous comment.\nHowever, I see two fundamental problems in the manuscript that are not related to the data, but to their interpretation and the framing of the entire manuscript.\n(1) The smaller of these problems concerns the discussion of the relevance of the findings with respect to glacier models and projections of glacier evolution. The manuscript claims that “Current analyses of future glacier behaviour predict a linear future contribution of global (8) and Alaskan ice masses to sea level rise (5,6,8), and a linear rate of volume loss for Juneau Icefield (9).”, and that “the simplified physics deployed by many of the models in Glacier MIP […] accounts poorly for variation across different hypsometries (11,14)”. The manuscript argues that the nonlinearities supposedly found in the data (see below for more on this) cannot be captured by the current generation of models. I strongly disagree with this assessment. The given references (5, 6, 8, 9) show nonlinearities in the mass change rates over the 21st century, such that invoking “such linear projections” does not adequately reflect their characteristics. It is true that regional and global scale projections “hide” nonlinear behavior of individual glaciers in the aggregated numbers typically shown. This does not mean, however, that this kind of behavior is not represented at the scale of individual glaciers (and it even shows in the aggregate, as stated above). In particular, the well-known feedbacks between glacier hypsometry/geometry and mass balance referred to in the manuscript are well represented in state-of-the-art global glacier models as used in references 6, 8, and 9. Some other nonlinearities are possible in some, but not all of the global-scale models (e.g., concerning disconnections). There are additional nonlinearities which are not captured well by any of this type of model (the manuscript mentions some), these are addressed, e.g., in ref. 14. But those nonlinearities are not the ones the manuscript claims to be affecting the Juneau Icefield.\nWe have refocused the manuscript’s motivation and conclusions to focus more on the empirical observations and have removed focus on the limitations of current numerical modelling attempts. While existing numerical simulations forecast a linear rise in rates of volume loss (cf. Rounce et al. abstract, “mass loss is linearly related to temperature increase”), even using PISM (Ziemen et al., 2016), and OGGM, which can simulate disconnections and more of the full processes, and while many of the other models in GlacierMIP are vastly simpler, some relying on just processes like volume-area scaling (Marzeion et al., 2020), we argue that the story presented by the impressive acceleration of Juneau Icefield is alone sufficiently interesting. We therefore simplify our introduction and discussion by focusing less on model limitations and more on our own empirical data.\nThis issue does not affect the center of the manuscript, and because of this, I see this as the “smaller” of the two fundamental problems. I see it as fundamental since it affects a main motivation of the manuscript, being prominently given in the abstract (and repeatedly mentioned in the manuscript). It should be possible to reframe the manuscript so that it is more adequately referring to the state-of-the-art in glacier modeling.\nWe hope that since we have removed this from the top and tail of the manuscript and leave the motivation as requiring longer term studies to understand in more detail the trends in area and volume change. we have therefore reframed the manuscript.\n(2) The bigger of the two problems concerns the claim that nonlinear processes are driving the mass loss acceleration of the Juneau Icefield (see title).\nThe manuscript has not convinced me that there is in fact a nonlinear response of the icefield observable in the data. Essentially, this conclusion would have to rest on the scatter plot in Fig. 4j. Here, a nonlinear response of the rates of area and thickness change to temperature change would have to be visible. With only 4 (thickness) and 6 (area) point available, and quite a bit of scatter (which might be explained by precipitation anomalies), the argument for nonlinearity is weak. The statement in lines 165 – 168 is therefore not justified; neither the part concerning the present manuscript, nor concerning the modeling studies (see above). Arguing the other way around (“nonlinear responses found in projections are not yet robustly detected in observations/reconstructions”) would be a similar misstatement, but arguably not worse than what is said in the manuscript.\nJust to note that the scatter is caused by precipitation anomalies.\nAs an alternative to showing a nonlinear response of the entire icefield (Fig. 4j), it would be possible to demonstrate that nonlinear effects are affecting the glaciers by showing that the ELA has crossed into the plateau region (i.e., by looking at ELA changes and how they relate to the hypsometry of the glacier), or by investigating the mass balance profile and ice flow velocities in the region of glacier disconnections (see explanation of feedbacks around Fig. 6). Unfortunately, this is not done in the manuscript.\nThe non-linear effects are indeed driven by the ELA crossing the plateau, which we had not adequately demonstrated in the previous version of the manuscript. The review of the ELA was rather buried in the Supplementary Information. We have moved the description of the vast previous work on Juneau Icefield ELAs (McNeil et al., 2020; McNeil et al., 2019; McNeil et al., 2016) and snowlines (Pelto, 2019) to the main text. Further to this, we have built on previous late-summer snowline mapping (Pelto, 2019) by analysing the 2019 late-summer snowlines; see new table in Supplementary Information and addition to figure 5. These ELA changes show the rising ELAs and their relationship to glacier hypsometry, with observed summer snowlines high on the plateau in 2019.\nIn terms of mass balance profiles at the site of disconnections, our data do show increased thinning at the mean height of the disconnections. Firstly, figure 6i shows that the rising ELAs have reached the mean elevation of the disconnections, resulting in increased numbers of observations, as seen in figure 6h. secondly, the elevation at which thinning occurs rises through the observational period (as shown in figure 7). We have strengthened the discussion of these results under ‘Icefield fragmentation’ in the Results.\nWe hope that this has strengthened our identification of the non-linear processes affecting the icefield that are outlined in Figure 8.\nIn this context, the section “Hypsometric controls of icefield behavior” is very problematic: the relevant mechanisms are well-known, but they are also explained (and referenced) well in the manuscript. However, the manuscript does not show that a single one of them is actually affecting the Juneau Icefield in their data. Remarkably, they use an edited photograph (Fig. 6) for illustrating “how it might look in a few decades”. That is exactly the point: in this manuscript, I don’t find any evidence that these feedbacks are active in present day. However, the manuscript very boldly claims that it is the case in the title, the abstract, scattered in the text, and in the final, conclusive paragraph.\nWe argue here that the hypsometry of the icefield is driving the acceleration of thinning across the plateau. As suggested by the reviewer above, here we again strengthen the argument by referring back to the observations of ELA and snowlines on the plateau, resulting in the processes and feedbacks identified in figure 8 (previously figure 6). This includes the reduced accumulation area and changed albedo as a result of the loss of high-albedo snow on the plateau – this is mapped in figure 5. Thinning and SMB feedback is then an inevitable response to the thinning already observed on the plateau. We have developed this section of the discussion to emphasise the evidence available for each of the processes named in the blue box in figure 8 (the schematic).\nWe also show increasing icefield fragmentation as glaciers disconnect, which increases bare rock as the icefield fragments, driving albedo changes, increasing debris cover, and changing ice flow. We also show increased thinning on glacier tongues. Our data therefore does show that these processes are affecting Juneau Icefield.\nThe problem of a missing demonstration of nonlinearity in the data may arise from a confusion of “nonlinearity” and “acceleration”. E.g., an acceleration in mass change can be a linear response to an increased temperature anomaly. The feedbacks explained in the second and third paragraph of the manuscript (and Fig. 6) would lead to an acceleration in mass change without additional changes in the temperature anomaly.\nAnd this is indeed what we have shown here.\nTherefore, I would suggest that the authors carefully decide whether they actually want to claim (and then have to demonstrate) a nonlinear response of the Juneau Icefield to climate change (which, with the data on hand, seems to be very difficult), or whether they wish to make a case for an impressive acceleration of glacier change in Alaska in the recent decades, which would be a very strong case.\nThank you for this constructive comment. We agree that our work shows an acceleration in mass change due to a steadily increasing temperature forcing. The hypsometrically controlled melt-accelerating feedbacks identified in the discussion have driven this impressive acceleration of melt. These feedbacks contribute to a hysteresis in the system, that will impede glacier regrowth in event of cooling. We have updated and reframed the abstract, discussion and Figure 8 (the schematic) to better reflect these findings, and focused on acceleration rather than ‘nonlinearity’. We therefore focus on making a very strong case for the impressive acceleration of glacier change in Alaska in recent decades.\nEither way, very substantial work on the text (and potentially, the data analysis) would be needed, so that I would probably consider the manuscript as “new”. This would imply a suggestion to reject the manuscript. However, as I am convinced that the work definitely should be published eventually, I would also agree with considering this as a case of “major revisions”. Expecting the manuscript to go through another round of revisions (potentially in a different format, pending the editor’s decision), I have refrained from specific/minor comments on all parts of the manuscript except the main text, where I hope also these specific comments will be helpful (also, I don’t see any major issues in the methods).\nThank you for these constructive comments, which we feel have improved the manuscript. We have reframed the discussion and abstract to focus on the observed accelerated mass change driven by hypsometrically controlled feedbacks, and have provided substantial new data and analysis on the ELAs and snowlines which are now intersecting the plateau.\nSpecific/minor suggestions for the main text:\nL18: „The majority…“ would imply > 50 %, which is not the case according to Hugonnet et al. (https://doi.org/10.1038/s41586-021-03436-z); perhaps change to: “The largest regional contribution…”\nAmended.\nL19: “…projected to become even more dominant…” I’m not sure I would agree, and would expect that this depends on the scenario: in a high-emission scenario, e.g., I would expect that regions that have more ice left (peripheries of Greenland and Antarctica, the Canadian Arctic) become more dominant over time.\nAlaska is the largest regional contributor to global mean sea level from 2015-2100 (Rounce et al., 2023). The Greenland periphery, Antarctic etc, contribute all less. Citation added to the introduction (abstracts are unreferenced).\nL22: “…with…” since this refers to the acceleration (but no acceleration value is given), perhaps change to: “…leading to…” Amended as suggested.\nL35: suggest to change to: “Under the higher-emission scenario RCP 4.5…” Amended as suggested.\nFig. 1 is a bit disorganized: e.g., I can’t find any ice-dammed lakes, lakes, or sandurs on the map (or photographs), features not visible have been removed from the legend.\nand referring to north- and south-flowing glaciers in panel d, without including a north arrow, is not ideal (perhaps label the glaciers in panel d directly?). North arrow added to figure. All maps are on the same projection.\nI would also expect that the individual glaciers mentioned in the text (e.g., lines 83-87) appear in the map, or vice versa. Please include a reference to (e.g.) Fig. 2a in the relevant sections of the text. We have added key glacier names to Figure 1 to aid clarity in the manuscript.\nL82-90: please indicate the uncertainties of the assessed changes in ELA, volume, and thickness. These uncertainties are important for understanding the significance of the changes. Amended as requested.\nFig. 2: I appreciate the difficulty to put the rich data set on the map(s), and unfortunately, I don’t have a good solution. But it is problematic that (i) the circles of the ELA estimates cover a large fraction of the trimline and moraine markers and (ii) moraine polygons are shown gray on gray background. In line with other reviewers, we have redesigned this figure and moved the ELA data to the SI. Constraints on numbers of display items make showing all the data challenging.\nL118: what do “M” and “SD” mean? What units are used here? Mean and standard deviation; clarified in the text.\nL129: It has not been shown that the acceleration is really exponential; I would suggest to state this more carefully (also in the abstract and other places in the manuscript). Amended to ‘steep’ throughout\nFig. 3: The tables in the lower left should either really be formatted as tables (outside of the figure) or left out (they are already shown in Fig. 4). Removed from figure as suggested. They’re also shown in the table moved from the SI to the main text.\nRegarding number of glaciers: how are separations treated? Do they increase the number of glaciers? In line with inventory best practice, they are retained as a single glacier but are a multipart polygon. If just the terminus is separated but they still share an accumulation area, then they remain as one single polygon.\nFig. 4a: please explain red line here, not just in panel f. Amended as requested. b: Where is the timing of an early, mid and late LIA shown? This was an error - this should refer to panel d. Amended. Last sentence of the caption: what does this refer to? Panel f does not show reduced precipitation during 2005 – 2015?\nRemoved from caption.\nFig. 5: legend and caption: please correct from dhdt to dh/dt. Amended as requested. More generally, I would suggest to move this figure to the supplement; it is not that central to the arguments in the manuscript. We argue that this figure is crucial as it demonstrates some of the key arguments in the text; this includes the thinning observed across the plateau, and also the relationship between the rising elevation of the thinning and the increasing frequency of observations of glacier disconnections."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Review of the manuscript “Accelerating glacier volume loss on Janeau Icefield driven by hypsometry and melt-accelerating feedbacks” by Bethan Davies et al.\nI appreciate the changes included in the manuscript by the authors and I think their work has now been improved. Nonetheless from my point of view there is still work to do. And I recommend major review accomplishing these important points (the first of them not demonstrated in the manuscript).\n1. Despite I agree with the hysteresis behavior of glaciers and that Juneau icefield regrowth is challenged due to the hypsometrical changes of the icefield, I consider that all the feedbacks have only been named (yes, with suitable references) but not shown with observations and data in the manuscript. This way, I consider it is not worth to claim about which are the main drivers of the acceleration. The authors clearly show the acceleration of the area and volume loss but do not analyze what is causing this acceleration (the only data shown in this regard is temperature rise). Such strong assertion (what causes the acceleration), needs to be verified with data (albedo trends, mean slope evolution, debris cover extent trend, ice flow evolution…).Otherwise, the title and the main conclusions must be consequently changed. Somehow, this is related with Reviewer 3 remarks. I completely agree about what he/she points out: “I don’t find any evidence that these feedbacks are active in present day”. From my perspective, the response of the authors and the changes included in the text do not answer this lack. No albedo change is shown, no debris cover increase is shown, no iceflow change is shown through the study period. I consider the manuscript must be conveniently changed in this regard before publications. For example, lines 433-435 should be removed and also some sections in the discussion. The named feedbacks can be described or presented as an hypothesis, but not as a trend observed.in this sense, I agree on the way of presenting them in lines 376 and 377:”…feedbacks CAN lead to an increase ….”. Overall, the melt accelerating feedbacks can be discussed, but one schematic figure (figure 8) do not justify these are driving (and will drive) Juneau Icefield volume loss acceleration and that now a “tipping point is reached controlled by a number of feedbacks” (stated in the abstract) as these have not been shown in the article.\n2. One potential way to show if a tipping point has been reached could be to show the ELA evolution in the last years. In lines 269-276, the manuscript shows the ELA observed in 2019. Why the ELA mean elevation is not updated until 2023? Sentinel images are available until yesterday (literally) and this will definitely show if the ELA is rising in the last years and area above/below the plateaus mean elevation is being reduced. Please, update this information to firmly show (or not) that this accelerating trend is being observed nowadays. I consider this is a relevant point, update the entire database of volume and area loses would be time consuming but Sentinel derived ELAs (potentially 2022 and 2023) is much faster to derive and will help to see the recent evolution of Juneau icefield.\nI still have a list of minor comments to be addressed:\nLine 247. Volume loss of 0.65±0.92 km2a-1 means that, a gain of 0.27 km2a-1 can be a true value during the study period. Such a big error compared to the mean value, is related to the method used to derive it or is it related to the interannual variability of glaciers mass balance? I think it is worth to comment this point.\nFigure 1. Please include a coast limits with a black line. Change the color to represent the ocean/water, as in low elevation areas it is not possible to differentiate land and ocean. A small location map (of north America which is easily recognizable) is needed in Figure 1 (it can be included in the upper right area of the upper panel). I think it is not needed to include all area vs elevations histograms for Alaskan ice fields. The cumulative area vs cumulative elevation is informative enough. I would only include the Juneau area elevation histograms.\nFigure 2: Lower right map. Symbols are difficult to see. Make them bigger and change colors (specifically disconnection in panel E have same color as the LIA extent and this is confusing).\nFigure 3: Almost impossible to identify (in printed version and online version) disconnection and separations with the symbols and colors used. Maybe there are not disconnections or separations in these areas, if this is the case, do not include these symbols in the legend please. Change conveniently.\nFigure 6 a): The temperature anomaly is computed for which reference time period? Please include it in the text\nFigure 4: answer to reviewer 1 comment on figure 2: Yes number of glaciers is important but I guess features are not in all cases glaciers. Clearly indicate that the number of features (in brackets) is not glacier number. I appreciate that the area has been included in the new version, this is informative.\nFigure 4: In the LIA area and glacier extent is not needed to include decimals if the error is three others of magnitude larger than the decimals included. Remove them please. Same for decimals in lines 223, 226, 259, Table 1 area changes, table 2 total volume loss, etc.… The 50 km scale bar is confusing. Edit it including smaller divisions (as those included in figure 5 or 7).\nFigure 4 caption. For clarity, include the period for which the ice surface lowering has been computed.\nFigure 5: Same comment of Figure 3. Disconnection symbols are difficult to identify. I would change colors. Also, I encourage to make more evident the difference between 1200 and 1500 m plateau elevations. In both printed version and in screen version, is barely observed the transition.\nFigure 7: Include the number of disconnections registered in each time period please.\nSupplementary figure 4: Despite in the reviewers’ response it is said that the location of the study area is included, it is not. Please mark in all ERA5 reanalysis maps where Juneau icefield is located. Also include a legend for it.\nResponse to points 2 and 3 (Lines 178 and 178-196 respectively) of reviewer 2: Please include more details on how these points have been addressed, also including the changes done in the text and why you consider is appropriate.", "reviewer_2": "This is my second time reviewing this article and I am happy that all my suggestions for improvement have been incorporated. As I stated previously, this is a significant piece of work, based on extensive, thoroughly analysed datasets, which support the conclusion and claims. The rearrangement and rewording of the paper has highlighted the significant key findings.\nA couple of typos:\nLine 110: should be \"in the same time period\" or \"at the same time\"\nLine 137/8: a 95% confidence interval, rather than the 95%...?\nLine 226: \"... the icefield is 315.\" rather than \"of\"", "reviewer_3": "I would like to thank the authors for having done a great job revising the manuscript following (or, where appropriate, arguing against) the reviewers’ comments, and I am very happy with the overall quality of the manuscript! Now even more than before I am convinced that these results should be published, and I see less issues with how some of the results are presented and discussed. The manuscript also reads very well and the figures have been improved substantially, so that I have hardly any minor suggestions. But I still see two points that I would strongly suggest to be changed before publication.\nThe first concerns the claim of a tipping point being reached. The manuscript makes now a much better case arguing that the tipping point is in reach (i.e., that ELAs are approaching the plateau elevations, and are crossing them in individual years), but some of the statements still are too bold. Mostly, I’m concerned with the last sentence of the abstract (“This icefield has now reached a tipping point controlled by a number of feedbacks, with hysteresis inhibiting glacier regrowth.”). This is a categorical statement, and while I would agree that it seems very likely that this threshold will be crossed (and in fact might be, and interannual variability clearly already pushes the ELA above threshold elevations), I don’t think the claim that the hysteresis is already affecting the SMB is acting right now, is substantiated as clearly in the manuscript as the abstract makes believe. In fact, within the main text, the statements are a lot more careful (“…a decrease in height of a flatter icefield or ice cap will theoretically result in a decrease in surface mass balance (SMB) as the ice surface lowers into warmer air, which reinforces itself as an SMB-elevation positive feedback”", "reviewer_or": "the mass balance-elevation feedback will inhibit glacier regrowth (cf.15), likely forming an irreversible tipping point in glacier dynamics", "reviewer_point": "Please check the legend of Fig 3, there is a double entry (sandur) and I’m not sure all the features in the legend actually occur in the maps (I don’t see disconnections and separations)."}, "author_response": "Response to reviewers\nReviewer #1 (Remarks to the Author):\nReview of the manuscript “Accelerang glacier volume loss on Janeau Icefield driven by hypsometry and melt-accelerang feedbacks” by Bethan Davies et al.\nI appreciate the changes included in the manuscript by the authors and I think their work has now been improved. Nonetheless from my point of view there is sll work to do. And I recommend major review accomplishing these important points (the first of them not demonstrated in the manuscript).\n1. Despite I agree with the hysteresis behavior of glaciers and that Juneau icefield regrowth is challenged due to the hypsometrical changes of the icefield, I consider that all the feedbacks have only been named (yes, with suitable references) but not shown with observaons and data in the manuscript. This way, I consider it is not worth to claim about which are the main drivers of the acceleraon. The authors clearly show the acceleraon of the area and volume loss but do not analyze what is causing this acceleraon (the only data shown in this regard is temperature rise). Such strong asseron (what causes the acceleraon), needs to be verified with data (albedo trends, mean slope evoluon, debris cover extent trend, ice flow evoluon…).Otherwise, the tle and the main conclusions must be consequently changed. Somehow, this is related with Reviewer 3 remarks. I completely agree about what he/she points out: “I don’t find any evidence that these feedbacks are acve in present day”. From my perspecve, the response of the authors and the changes included in the text do not answer this lack. No albedo change is shown, no debris cover increase is shown, no iceflow change is shown through the study period. I consider the manuscript must be conveniently changed in this regard before publicaons. For example, lines 433-435 should be removed and also some secons in the discussion. The named feedbacks can be described or presented as an hypothesis, but not as a trend observed.in this sense, I agree on the way of presenng them in lines 376 and 377:”…feedbacks CAN lead to an increase ….”. Overall, the melt accelerang feedbacks can be discussed, but one schemac figure (figure 8) do not jusfy these are driving (and will drive) Juneau Icefield volume loss acceleraon and that now a “pping point is reached controlled by a number of feedbacks” (stated in the abstract) as these have not been shown in the arcle.\nAs also suggested by Reviewer 3, we have rephrased the abstract (please see comment to reviewer 3 below) regarding the pping point in the abstract, in order to be more clear about what the arcle data shows:\n“Rising equilibrium line altudes and increasing ablaon across the plateau has driven a series of hypsometrically controlled melt-accelerang feedbacks and resulted in the observed acceleraon in mass loss. As glacier thinning on the plateau connues, a mass balance-elevaon feedback is likely to inhibit future glacier regrowth, potenally pushing glaciers beyond a dynamic pping point.”\nWe have undertaken new work to show both snowline and albedo change across mulple decades, and present this new data in the Supplementary informaon. We have also brought the snowline and albedo data together under a new heading in the Results secon, to emphasise the new data and to differenate it from the glacier area and volume change. We believe that this strengthens our argument in the manuscript.\nIn the Results and from published datasets, we have therefore more clearly demonstrated that the ELA and snowline are rising to intersect with the plateau, the albedo is decreasing, and document the fragmentaon of glacier ice and the diminishing ice cover due to glacier recession.\nIn the discussion, we are careful to, as the reviewer suggests, present the feedbacks as a hypothesis, as we agree that we have not demonstrated albedo change (bold text denotes the changes):\n“We hypothesise that the dramac reducon in icefield accumulaon area is decreasing overall icefield albedo. As snow cover duraon and extent is reduced (cf. Figure 5), darker rock, firn and ice surfaces are increasingly exposed, which absorb more solar radiaon, enhancing the ablaon rate. This albedo feedback (cf.1,2 and illustrated in Figure 8) is likely to then contribute to further icefield thinning. This effect is compounded by deposion of black carbon3 and dust4, which is set to increase in future, due to a rise in tourism and wildfires5.\nThinning on low-slope icefield plateaux, as observed in Figure 7, is likely to drive an elevaon-SMB feedback by lowering the icefield surface and driving increased melt in warmer air temperatures (due to the lapse-rate induced rise in air temperature6-9) (Figure 8). This response is fundamentally unstable, as thinning drives increasingly negave mass balances, even in the absence of further climate forcing9. In some places across the Juneau plateau, the ice is over 600 m thick10 (Supplementary figure 2), so the ice-elevaon feedback has a large range over which to operate. This topographic feedback can lead to an increased frequency of extreme negave mass balance rates, especially when the climate data diverges from the mean6. Once iniated, long-term thinning is inevitable, because even in the event of climate stabilisaon, the mass balance-elevaon feedback will inhibit glacier regrowth (cf.11), likely forming an irreversible pping point in glacier dynamics.”\nWith regards to debris cover, a stascal relaonship between glacier disconnecon and glacier debris cover was demonstrated in a previous publicaon by these authors, which is cited here. We have added the citaon more clearly when needed. We have also been more careful to cite classic papers demonstrang the influence of debris on glacier melt (e.g., 12) (changes in bold):\n“Thin debris accumulates on the ice surface, darkening the surface and thereby enhancing melt. In some places, debris accumulates to the point that the glacier becomes debris-covered (>50% of the tongue is covered in thin debris13). Thicker debris cover on glacier ice encourages in situ down- wasng but slows rates of glacier area loss due to insulang impact of debris cover on surface melt12. The proporon of debris-covered glacier area to clean-ice glacier area will likely increase as glacier shrinkage connues, increasing the relave importance of debris cover.” … “We hypothesise that increased bare rock exposure and the loss of ice within and around glaciers will act to darken the region, reducing regional albedo, which may contribute to further regional warming.”\nWe therefore have, as the reviewer suggests, taken care to ensure that the feedbacks are discussed and presented as a hypothesis, and not an observaon.\nUnder the secon headed ‘Wider Implicaons’, we have amended lines 433-435 to remove ‘as our observaons demonstrate’. These feedbacks are well understood and are cited in the arcle extensively so this sentence is jusfied.\n2. One potenal way to show if a pping point has been reached could be to show the ELA evoluon in the last years. In lines 269-276, the manuscript shows the ELA observed in 2019. Why the ELA mean elevaon is not updated unl 2023? Sennel images are available unl yesterday (literally) and this will definitely show if the ELA is rising in the last years and area above/below the plateaus mean elevaon is being reduced. Please, update this informaon to firmly show (or not) that this accelerang trend is being observed nowadays. I consider this is a relevant point, update the enre database of volume and area loses would be me consuming but Sennel derived ELAs (potenally 2022 and 2023) is much faster to derive and will help to see the recent evoluon of Juneau icefield.\nThere is already an extensive and long-term monitoring dataset over Juneau icefield, which shows the rising ELAs which are intersecng the plateau. This is shown for example in figure 6i – here, the measured ELAs for two glaciers and the height of the plateau is emphasised. You can clearly see the rising ELAs intersecng the plateau. Snowline elevaons are a proxy for these more detailed datasets. There are also long-term datasets measuring transient snowlines on Taku Glacier from Pelto 2019 14 which support this analysis.\nIn order to thoroughly address this point, we mapped new end of summer snowlines from 2020-2023, meaning that there is a new icefield wide dataset from 2019-2023. These data are shown in detail in the revised Supplementary Informaon secon 5.6, with updated tables and maps showing the snowline elevaon. We also updated the table of snowline informaon in the Source Data. These new datasets show that snowlines have connued to reach the plateau in recent years, though 2019 was an unusually high year. We have supported these new observaons by comparing these icefield wide snowline datasets to the long record of glacier mass balance recorded by the USGS at Lemon Creek Glacier and Taku Glacier, with new text and two new figures added to the Supplementary Informaon:\n“Supplementary figure 14a a shows the elevaons of the icefield wide mapped end-of-summer snowlines in comparison with mapped snowlines from Taku Glacier and Lemon Creek Glacier from 1946-2023 (Supplementary figure 14b), using an updated dataset from Pelto (2019)14, and measured equilibrium line altude for these two glaciers over the same me period, using data derived from the USGS Benchmark Glacier Program15; see also reference publicaons16,17 (Supplementary figure 14c). These datasets show that, although the year 2019 was an unusually high year for snowline elevaon, snowlines are reaching the plateau frequently, with both snowlines and ELAs occurring above the plateau height of 1200 m since 2005.”\nWe have also updated the main manuscript to note this important point, with snowlines and ELAs occurring on the plateau regularly over the last 15-20 years.\nWe believe that this change has strengthened the argument in the manuscript and does indeed show that the icefield has reached a threshold, with increasingly high snowlines driving change.\nI sll have a list of minor comments to be addressed:\nLine 247. Volume loss of 0.65±0.92 km2a-1 means that, a gain of 0.27 km2a-1 can be a true value during the study period. Such a big error compared to the mean value, is related to the method used to derive it or is it related to the interannual variability of glaciers mass balance? I think it is worth to comment this point. The large uncertainty here is due to the interpolaon of the datasets over areas with no data. The LIA DEM is calculated only over the icefield ablaon areas, and the 1948 DEM is likewise not available over Canada. This does not mean that the change is insignificant, but we emphasise that the uncertainty is large because of interpolaon, not because of variability of the elevaon differences.\nWe have added a note to the Methods to clarify this point, as writen above.\nFigure 1. Please include a coast limits with a black line. The complexity of the coastline makes this very difficult to show at this scale; it ends up dominang the figure. Instead, we have added a colour to the ocean, which makes the land stand out much more.\nChange the color to represent the ocean/water, as in low elevaon areas it is not possible to differenate land and ocean. See above comment. We have coloured the ocean to beter differenate the water and land.\nA small locaon map (of north America which is easily recognizable) is needed in Figure 1 (it can be included in the upper right area of the upper panel). Amended as requested.\nI think it is not needed to include all area vs elevaons histograms for Alaskan ice fields. The cumulave area vs cumulave elevaon is informave enough. I would only include the Juneau area elevaon histograms.\nWe think that these histograms illustrate the argument well – the cumulave graphs alone are difficult to interpret.\nFigure 2: Lower right map. Symbols are difficult to see. Make them bigger and change colors (specifically disconnecon in panel E have same color as the LIA extent and this is confusing). The symbols were enlarged to improve their visibility. The colour is different to the LIA and this is now clearer; they are a unique colour.\nFigure 3: Almost impossible to idenfy (in printed version and online version) disconnecon and separaons with the symbols and colors used. Maybe there are not disconnecons or separaons in these areas, if this is the case, do not include these symbols in the legend please. Change conveniently. Thanks for picking up the error – the disconnecons were not layered appropriately in the GIS table of contents. This is now corrected. The symbols were also enlarged to improve their visibility.\nFigure 6 a): The temperature anomaly is computed for which reference me period? Please include it in the text. The anomaly is computed against the year 1000-1998 mean, as done in the original publicaon (see y axis for the panel a). In panel f, the temperature anomaly is computed against the 1986-2005 mean, as is commonly used for the IPCC. This is now noted addionally in the capon.\nFigure 4: answer to reviewer 1 comment on figure 2: Yes number of glaciers is important but I guess features are not in all cases glaciers. Clearly indicate that the number of features (in brackets) is not glacier number. Noted in the legend.\nI appreciate that the area has been included in the new version, this is informave. Thank you.\nFigure 4: In the LIA area and glacier extent is not needed to include decimals if the error is three others of magnitude larger than the decimals included. Remove them please. Same for decimals in lines 223, 226, 259, Table 1 area changes, table 2 total volume loss, etc.… Amended to give results to 1 d.p.\nTable 1 is now given to 1 d.p. for total glacier area, mean glacier area and change in area.\nHowever, in Table 2 I think that for mean area-averaged dh/dt, two decimal places are needed due to the smaller numbers in the confidence interval.\nThe 50 km scale bar is confusing. Edit it including smaller divisions (as those included in figure 5 or 7).\nAmended as requested.\nFigure 4 capon. For clarity, include the period for which the ice surface lowering has been computed.\nAmended as requested.\nFigure 5: Same comment of Figure 3. Disconnecon symbols are difficult to idenfy. I would change colors.\nThe colours are unique – no other colours are really available to enhance the change. Instead I have enlarged the disconnecon symbols to make them more visible.\nAlso, I encourage to make more evident the difference between 1200 and 1500 m plateau elevaons. In both printed version and in screen version, is barely observed the transion. We have darkened the plateau colours to enhance visibility.\nFigure 7: Include the number of disconnecons registered in each me period please. Amended as requested. However, note that as the date of disconnecon was calculated using the satellite imagery, the mings do not exactly match.\nSupplementary figure 4: Despite in the reviewers’ response it is said that the locaon of the study area is included, it is not. Please mark in all ERA5 reanalysis maps where Juneau icefield is located. Also include a legend for it. The locaon of Juneau is already shown in this figure, with a yellow star in panel A. This is noted in the capon. To emphasise this further we have moved the note on this to the second sentence of the capon and added a further label to panel A.\nResponse to points 2 and 3 (Lines 178 and 178-196 respecvely) of reviewer 2: Please include more details on how these points have been addressed, also including the changes done in the text and why you consider is appropriate.\nWe welcome the opportunity to provide more informaon here.\nLines have now changed as extensive informaon was moved from the Supplementary Informaon to the main manuscript. These points refer to the third and fourth paragraphs of the Discussion, under ‘hypsometric controls on icefield behaviour’.\nPoint 2 in the first round of reviews quesoned whether the dramac reducon in icefield accumulaon area referred to our data or to predicted changes. This is because in the first submission of the manuscript we referred back to another paper by our team that looks at the future mass balance of the icefield (Ing et al., The Cryosphere). In the second submission of the manuscript, for clarity about what our data explicitly and aer receiving the reviewers’ comments, the reference to Ing et al has been removed.\nTherefore, in Point 2, where we are talking about the mass balance changes on the icefield, we are now more clearly only discussing our own data, rather than the predicted changes from Ing et al. The second paragraph under ‘Hypsometric controls’ now more clearly refers back to the published mass balance work from McNeil et al., 2020 and uses this to contextualise and understand the processes acng on the icefield.\nThe Reviewer 2 point 3 comment is related, and discusses again the need to clarify if this is referring to our own data or to the mass balance modelling work. We emphasise that we are referring to our own data by linking back to the figure that shows the 2019 AD snowline.\nReviewer #2 (Remarks to the Author):\nThis is my second me reviewing this arcle and I am happy that all my suggesons for improvement have been incorporated. As I stated previously, this is a significant piece of work, based on extensive, thoroughly analysed datasets, which support the conclusion and claims. The rearrangement and rewording of the paper has highlighted the significant key findings.\nWe thank the reviewer for their generosity and me in reviewing the manuscript and for their construcve comments.\nA couple of typos:\nLine 110: should be \"in the same me period\" or \"at the same me\". Amended as requested.\nLine 137/8: a 95% confidence interval, rather than the 95%...? Amended as requested.\nLine 226: \"... the icefield is 315.\" rather than \"of\". Amended as requested.\nReviewer #3 (Remarks to the Author):\nI would like to thank the authors for having done a great job revising the manuscript following (or, where appropriate, arguing against) the reviewers’ comments, and I am very happy with the overall quality of the manuscript! Now even more than before I am convinced that these results should be published, and I see less issues with how some of the results are presented and discussed. The manuscript also reads very well and the figures have been improved substanally, so that I have hardly any minor suggesons. But I sll see two points that I would strongly suggest to be changed before publicaon.\nWe thank the reviewer for their generous and construcve comments and are delighted that they see the high quality of the manuscript.\nThe first concerns the claim of a pping point being reached. The manuscript makes now a much beter case arguing that the pping point is in reach (i.e., that ELAs are approaching the plateau elevaons, and are crossing them in individual years), but some of the statements sll are too bold. Mostly, I’m concerned with the last sentence of the abstract (“This icefield has now reached a pping point controlled by a number of feedbacks, with hysteresis inhibing glacier regrowth.”). This is a categorical statement, and while I would agree that it seems very likely that this threshold will be crossed (and in fact might be, and interannual variability clearly already pushes the ELA above threshold elevaons), I don’t think the claim that the hysteresis is already affecng the SMB is acng right now, is substanated as clearly in the manuscript as the abstract makes believe. In fact, within the main text, the statements are a lot more careful (“…a decrease in height of a flater icefield or ice cap will theorecally result in a decrease in surface mass balance (SMB) as the ice surface lowers into warmer air, which reinforces itself as an SMB-elevaon posive feedback”"}, {"round": 3, "reviewer_comments": {"reviewer_1": "After reading the response to referees letter and the updated version of the manuscript, I would like to congratulate manuscript authors as I consider now the work ready for publication. The results have high interest for the scientific community and are conveniently supported by the data and the analysis. Moreover, some changes (i.e. new figure 8, albedo trend analysis in the supplementary material...) are highly suitable to show the main conclusions. Moreover, I can say that I have enjoyed the read of the final manuscript. Thanks for your effort and your time on rising this article.", "reviewer_3": "I would like to thank the authors for their effort to better express the somewhat hypothetical nature of their findings regarding the importance of the different feedbacks for the acceleration of mass loss. Their responses and edits to the manuscript (including the new design of Fig. 8) are adequate and I now find the manuscript of suitable quality to be published. I have no remaining questions or comments."}, "author_response": ""}]

Decision letter