Poleward shift of subtropical highs drives Patagonian glacier mass loss

Domain: earth_science
Published: 2025-04-23
Editorial outcome: minor_revision

Abstract
Patagonian glaciers have been rapidly losing mass in the last two decades, but the driving processes remain poorly known. Here we use two state-of-the-art regional climate models to reconstruct long-term (1940-2023) glacier surface mass balance (SMB), i.e., the difference between precipitation accumulation, surface runoff and sublimation, at about 5 km spatial resolution, further statistically downscaled to 500 m. High-resolution SMB agrees well with in-situ observations and, combined with solid ice discharge estimates, captures recent GRACE/GRACE-FO satellite mass change. Glacier mass loss coincides with a long-term SMB decline (−0.35 Gt yr−2), primarily driven by enhanced surface runoff (+0.47 Gt yr−2) and steady precipitation. We link these trends to a poleward shift of the subtropical highs favouring warm northwesterly air advections towards Patagonia (+0.14 °C dec−1 at 850 hPa). Since the 1940s, Patagonian glaciers have lost 1350 ± 449 Gt of ice, equivalent to 3.7 ± 1.2 mm of global mean sea-level rise.

Introduction
Situated in the Southern Andes, Patagonia hosts the largest glacier area within the southern hemisphere, excluding the Antarctic ice sheet1–3. Besides small mountain glaciers, Patagonia encompasses three major ice fields, namely the Northern (NPI; ~4000 km2)4 and Southern Patagonian Icefields (SPI; ~13,000 km2)5,6, and the Cordillera Darwin Icefield (CDI; ~2600 km2)7 located in the Tierra del Fuego archipelago (inset maps in Fig. 1a). Patagonian glaciers are on average 181 m thick2, with a maximum of 1400 m in deep glacial valleys3. Altogether, these glaciers hold 4772 ± 1242 Gt of ice8, eventually raising global sea level by ~13 mm if totally melted. Extrapolation of glaciological and geodetic observations suggests that Patagonia has experienced glacier mass loss since the early 1960s, contributing about 3.3 mm to global sea-level rise in the period 1961-20169. The latest assessment from the Glacier Mass Balance Intercomparison Exercise (GlaMBIE)8 reported rapid glacier mass loss of 26.5 ± 6.5 Gt yr−1 in the last two decades (2000-2023), equivalent to an average 0.07 mm sea-level rise each year. GlaMBIE estimates combine 13 mass change datasets including remote sensing (altimetry and gravimetry), geodetic and glaciological measurements. The three major ice fields are responsible for 83% of the total mass loss10, while the remainder is attributed to smaller neighbouring mountain glaciers. Storing only 3% of the Earth’s total ice volume outside polar ice sheets2, Patagonia is disproportionately responsible for ~10% of the glacier contribution to sea-level rise (2002-2016)11.

Mass change (MB) of Patagonian glaciers can be quantified by subtracting solid ice discharge (D), i.e., frontal ablation of calving icebergs, from the glacier surface mass balance (SMB), i.e., the difference between mass gained from precipitation accumulation and lost from meltwater runoff and sublimation. NPI and SPI ice fields generally discharge ice into the ocean on the western side of the Andes, while they typically flow into proglacial lakes on the eastern side1. Remotely sensed solid ice discharge estimates indicate that NPI and SPI accounted for 15.1 ± 1.1 Gt yr−1 on average in 1968–19991, which further increased to 24.1 ± 1.7 Gt yr−1 in 2000–201912. SPI glaciers contributed about 80% to the total solid ice discharge flux in 1968–19991, and up to 90% in 2000-201912. While remaining steady across most SPI glacier calving fronts, solid ice discharge has however experienced an overall decreasing trend since the 2000s, dominated by the continuous advance of a single glacier, namely Pio XI South12. This decreasing trend reveals that declining glacier SMB has been playing an increasing role in the recent Patagonian glacier mass loss13, which is expected to further accelerate in the future14,15. Previous NPI and SPI SMB reconstructions using positive degree day (PDD)16, regional climate models17, sometimes downscaled using fixed lapse rates to force simplified mass balance models13,18,19 (Methods), e.g., neglecting rain and meltwater retention in snow, suggested a positive but insignificant SMB trend in the period 1975-2015. This trend was ascribed to long-term snowfall increase combined with meltwater runoff decline. However, these results remain highly uncertain since low-resolution climate models are prone to overestimate precipitation over the steep Andes20, while they cannot accurately capture high meltwater runoff rates across narrow glaciers17. Furthermore, a positive SMB trend contradicts the recent mass loss increase observed in the last two decades10,21–25, highlighting large uncertainties in modelled SMB and a poor understanding of mass change drivers.

Here we use two state-of-the-art regional climate models, namely MAR version 3.14 (1940-2023) and RACMO version 2.3p2 (1979-2023)17,26, to reconstruct the contemporary SMB of Patagonian glaciers at 5 km and 5.5 km spatial resolution, respectively (Supplementary Fig. 1) (Methods). Both models are forced by the latest ERA5 climate reanalysis27 and statistically downscaled to a 500 m grid28,29 (Methods). In brief, statistical downscaling uses daily, local vertical gradients to correct SMB and components, i.e., total precipitation, snowfall, total sublimation, total melt and subsequent runoff, for elevation difference between the low-resolution surface topography prescribed in both models, and the high-resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) version 430 at 30 m resolution, down-sampled to a 500 m grid (Supplementary Fig. 2). Moderate Resolution Imaging Spectroradiometer (MODIS) albedo records at 500 m resolution are further used to correct surface melt and runoff for the albedo of exposed dark bare ice in summer. In line with previous studies17,20, RACMO strongly overestimates orographic-forced precipitation across the steep Andes. The statistically downscaled RACMO total precipitation is therefore adjusted and reduced by 36.5% over ice-covered grid-cells on a daily basis (Methods). This approach ensures that the precipitation adjustment has a dominating effect in regions where overestimates peak, notably across the steep and moist windside flank of the Andes, i.e., accumulation > 15 m water equivalent (w.e.) (Supplementary Fig. S1b). The justification for this adjustment is to eliminate the large positive bias relative to mass change records from the Gravity Recovery and Climate Experiment (GRACE) and Follow-On (GRACE-FO) missions24 (Supplementary Fig. 4b, c), when combined with solid ice discharge1,12. Such adjustment is not required for the statistically downscaled MAR product (Supplementary Fig. S4b). Hereafter, we use the high-resolution original MAR and precipitation-adjusted RACMO downscaled SMB products to identify the drivers of Patagonian glacier mass loss, and estimate their contribution to global sea-level rise.

Results
Mass change of Patagonian glaciers
Figure 1a shows long-term Patagonian glacier monthly cumulative mass change (MB = SMB-D) derived from MAR (green line; 1940-2023) and RACMO SMB (blue line; 1979-2023), further statistically downscaled to 500 m (Fig. 2c and Supplementary Fig. 3c), combined with published solid ice discharge estimates1,12. Since published estimates do not fully cover the period 1940-2023, we assume that the long-term solid ice discharge averages for 1968–19991 and 2000-201912 remain valid prior 1968 and post 2019, respectively. Therefore, we set solid ice discharge to 15.1 ± 1.1 Gt yr−1 in 1940-19991, with a step increase to 24.1 ± 1.7 Gt yr−1 in 2000–202312. Note that these published solid ice discharge estimates, i.e., annual averages for two periods, may not fully capture interannual and seasonal variabilities within their uncertainty range. Since 1940, Patagonian glaciers have been experiencing sustained mass loss (MB < 0), though briefly interrupted by a short mass gain episode (MB > 0) in the period 1945–1955 (Fig. 1a). This long-term negative trend is in line with geodetic mass change records spanning 1961–2016 (grey dots)9, with averaged observed and modelled mass loss of 26.2 ± 11.0 Gt yr−1 (Geodetic)9 and 22.3 ± 5.4 Gt yr−1 (MAR), respectively. Modelled mass loss from MAR (28.8 ± 5.4 Gt yr−1) and RACMO (29.9 ± 5.4 Gt yr−1) aligns well with GRACE/GRACE-FO satellite measurements (27.8 ± 11.0 Gt yr−1) (red line) in the last two decades (2002–2023), showing both a strong correlation (R2 = 0.98) and a low RMSE (~35 Gt) for both models (Supplementary Fig. 4c). Furthermore, our estimates are in excellent agreement with previously published mass change records including the latest GlaMBIE assessment8 (Supplementary Table 2), that reports mass loss of 26.5 ± 6.5 Gt yr−1 in 2000–2023 (orange dots) compared to 27.3 ± 5.4 Gt yr−1 in MAR and 28.5 ± 5.4 Gt yr−1 in RACMO, respectively. Using our two high-resolution products, we estimate that Patagonian glaciers have lost 1350 ± 499 Gt of ice since 1940, contributing 3.7 ± 1.2 mm to global sea-level rise.

Long-term surface mass balance trends
To explore the role of surface processes in the long-term mass loss, Fig. 1b shows annual cumulative SMB and components from both MAR (solid lines; 1940-2023) and adjusted RACMO models (bands; 1979-2023), statistically downscaled to 500 m (Fig. 2a–c and Supplementary Fig. 3a–c). For SMB evaluation, we use 74 stake measurements (1980-2019) collected at 38 in-situ sites across NPI (2 sites), SPI (26 sites), and CDI (19 sites) (yellow stars in the inset maps of Fig. 1a) (Methods). MAR and RACMO agree well with observations, i.e., R2 reaches 0.83 and 0.74, and RMSE ranges from 3.3 to 4.9 m w.e. respectively (Supplementary Fig. 5b, e). Cross-model comparison between each SMB component in the overlapping period 1979-2023 shows good agreement, with high correlation (0.92 ≤ R2 ≤ 0.94) and overall small model differences (1 to 7%) (Supplementary Fig. 6a–e) (Methods). Although glacier-integrated differences are small, low correlation (R2 = 0.32) is found for modelled retention and refreezing in firn (Supplementary Fig. 6f), i.e., the perennial compressed snow layer covering the glacier accumulation areas (SMB > 0 in the inset maps of Fig. 1a). We attribute these discrepancies to the different snow modules incorporated in MAR and RACMO (Methods). Supplementary Fig. 3d–f shows spatial differences between downscaled SMB components from MAR and RACMO at 500 m. Total precipitation shows large-scale spatial differences between the two models, i.e., with lower values in MAR on the western windside slope, and higher values on the eastern leeside (Supplementary Fig. 3d). This indicates a more pronounced foehn effect on the Patagonian Andes from prevailing westerly winds in RACMO, while MAR distributes precipitation more uniformly (Fig. 2a). This results in total precipitation being 8.0 Gt larger in MAR than in adjusted RACMO (Supplementary Fig. 6b, c and Supplementary Table 1). In addition, MAR has higher runoff than RACMO (5.6 Gt yr−1 in Supplementary Fig. 6d and Supplementary Table 1), notably across small mountain glaciers and at the ice fields margins (Supplementary Fig. 3e). This results from both higher meltwater production (1.7 Gt yr−1 in Supplementary Fig. 6e) and rainfall in MAR (1.3 Gt yr−1 in Supplementary Fig. 6c), combined with lower refreezing in firn (2.6 Gt yr−1 in Supplementary Fig. 6f). Although having a small contribution to SMB, we find that total sublimation is overall 2.1 Gt lower in MAR than in RACMO (Supplementary Table 1). Consequently, glacier-integrated SMB is almost equal in both models (0.3 Gt yr−1 difference in Supplementary Fig. 6a), and spatial differences are mostly driven by precipitation distribution patterns (Supplementary Fig. 3f). To ensure reliable surface accumulation, we compare modelled glacier-integrated total precipitation with corresponding estimates from the gridded meteorological PMET data set at 5 km spatial resolution (1980–2020)31 (Supplementary Fig. 7). We find excellent model agreement (R2 = 0.86 and 0.91) with low RMSE of 8.7 Gt yr−1 and 14.0 Gt yr−1, i.e., equivalent to 5% and 8% of the PMET total precipitation, in MAR and RACMO respectively.

Using long-term MAR time series (Fig. 1b), we find that Patagonian glacier SMB (Fig. 2c) remained close to equilibrium at 1.6 ± 0.3 Gt yr−1 on average (1940-2023), the result of total precipitation (164.7 ± 8.0 Gt yr−1) (Fig. 2a), i.e., partitioned in snowfall (79%) and rainfall (21%), approximately balancing surface runoff (163.0 ± 5.6 Gt yr−1) (Fig. 2b) (Supplementary Table 1). At the same time, SMB experiences large interannual variability (SD = 25.7 Gt yr−1), primarily driven by runoff fluctuations (SD = 21.1 Gt yr−1) being 40% larger than those of total precipitation (SD = 15.2 Gt yr−1). Large interannual variability is illustrated by peak high and low SMB reaching +59.4 Gt in year 1948, and -52.2 Gt in year 2016 (Fig. 1b). In these years, SMB is two to three-folds larger than the corresponding solid ice discharge flux (D = 15.1 to 24.1 Gt)1,12, indicating that SMB fluctuations primarily drive the variability and trend of Patagonian glacial mass change. In the period 1940–2023, we find that Patagonian glaciers have undergone sustained, and significant SMB decline (-0.35 ± 0.23 Gt yr−2, p < 0.01) (Fig. 1b). Our product suggests that long-term SMB decline is driven by insignificant total precipitation increase (0.12 ± 0.13 Gt yr−2, p > 0.05) notably in the glacier interior accumulation zones (Fig. 2d), combined with significantly enhanced surface runoff (0.47 ± 0.16 Gt yr−2, p < 0.01) across low-lying glacier ablation zones (Fig. 2e). Interestingly, the insignificant precipitation trend stems from steady snowfall accumulation (0.03 ± 0.12 Gt yr−2, p > 0.05) with significant rainfall increase (0.09 ± 0.03 Gt yr−2, p < 0.01). The increase in runoff coincides with a significant rise in annual mean glacier near-surface air temperature (T2m) (0.10 ± 0.03 oC dec−1, p < 0.01), hence triggering enhanced meltwater production (0.42 ± 0.14 Gt yr−2, p < 0.01). Atmospheric warming is in line with two gridded meteorological datasets at 5 km spatial resolution, namely CR2MET (1960-2021)32 and PMET (1980-2020)31 (Fig. 1c). We conclude that, following long-term atmospheric warming (1940-2023), enhanced runoff at the glacier margins (Fig. 2e) exceeds the small precipitation increase inland (Fig. 2d), in turn significantly reducing SMB (Fig. 2f).

Drivers of long-term runoff increase
Patagonian glacier SMB strongly correlates with surface runoff in both MAR and RACMO (r = 0.80) (Fig. 3a), while the signal is much weaker for total precipitation (r = 0.61) (Fig. 3b). Our products thus indicate that long-term Patagonian glacier SMB responds predominantly to runoff rather than precipitation changes (Fig. 2d–f). At the same time, runoff is highly sensitive to glacier near-surface temperature anomalies with respect to 1960-1989 (r = 0.81) (Fig. 3c), indicating that long-term atmospheric warming triggers runoff increase. To investigate the link between atmospheric warming and runoff, Fig. 4a–c shows long-term ablation zone extent (SMB < 0) as a fraction of the total glacier area; the fraction of total precipitation falling as rain; and the glacier refreezing capacity, i.e., the fraction of rain and meltwater effectively retained or refrozen within the firn layer. The latter three quantities are first correlated to near-surface temperature anomalies (Fig. 4d–f) and then to surface runoff (Fig. 4g–i). MAR and RACMO agree well in the overlapping period 1979-2023, with strong correlation and small model differences for the ablation zone extent (R2 = 0.83, difference = -2.4%), rainfall fraction (R2 = 0.96, difference = -0.2%), and firn refreezing capacity (R2 = 0.70, difference = -0.8%) (Fig. 4a–c).

Using long-term MAR time series, we find that the glacier ablation zones have significantly expanded since 1940 (+0.9 ± 0.5% dec−1, p < 0.01), following near-surface atmospheric warming (r = 0.70, Fig. 4d). This stems from a significant long-term increase of the equilibrium line altitude (ELA, SMB = 0) by 8.9 ± 5.5 m dec−1 (p < 0.01). As a result, the ELA migrated upward from ~1220 m a.s.l. to ~1300 m a.s.l. between 1940 and 2023, further expanding the ablation zone area. This is further supported by the fact that MAR captures, within uncertainties, the ELA of NPI and SPI derived from previous observational and model estimates (Supplementary Table 3). Ablation zone expansion triggers runoff increase (r = 0.78, Fig. 4g) as enhanced melt eventually drains on top of impermeable bare ice exposed at the surface in summer. Furthermore, the contribution of rainfall to total precipitation has significantly increased (+0.4 ± 0.2% dec−1, p < 0.01), the result of atmospheric warming favouring precipitation in the liquid phase (r = 0.67, Fig. 4e), which preferentially runs off on top of bare ice surfaces (r = 0.84, Fig. 4h). In the interior accumulation zone, firn retains roughly 18% of surface melt and rain in its pore space on average (1940-2023). Since 1940, firn refreezing capacity has significantly declined (-0.3 ± 0.1% dec−1, p < 0.01) (Fig. 4c) following atmospheric warming (r = 0.70, Fig. 4f). This results in increased surface runoff across the glacier interior accumulation zone (r = 0.86, Fig. 4i) through i) the progressive removal of the firn layer, ii) a reduced firn replenishment as rainfall increases at the expense of snowfall, and iii) the depletion of available firn pore space following enhanced surface melt and rain storage. Through these combined processes, near-surface atmospheric warming drives enhanced surface runoff not only in low-lying ablation zones, but also, though to a smaller degree, in the higher interior accumulation zones (Fig. 2e).

Poleward shift of subtropical highs drives glacier mass loss
Atmospheric warming strongly correlates throughout the air column from the near-surface (T2m) up to the 850 hPa level (T850) (r = 0.95) (Fig. 3d), suggesting a large-scale circulation change origin besides the effect of global warming. Figure 5a shows long-term spatial correlation between T850 derived from the ERA5 reanalysis27 and glacier-integrated surface runoff from MAR at 500 m (1940-2023). Strong links are generally found across the South Pacific and Atlantic Oceans, notably nearby the Southern Andes (r > 0.75, yellow contour in Fig. 5a and Supplementary Fig. 8a). Similar patterns hold for RACMO glacier runoff in the period 1979-2023 (Supplementary Fig. 8c). To explore the drivers of this atmospheric warming, Fig. 5b maps long-term trends in T850 (°C dec−1; background colour), geopotential height Z850 (m dec−1; grey contours), and wind speed and direction (arrows) derived from ERA5 in 1940-2023. The Southern Pacific Subtropical Gyre, i.e., also called Saint Helena High (SHH in Fig. 5b; 30–45oS/130–160oW) experiences enhanced Z850 (> 2.5 m dec−1), indicating a long-term strengthening. This is supported by stronger anticlockwise atmospheric circulation, and associated synoptic warming through enhanced air subsidence (> 0.2 oC dec−1). Superimposed on this, the SHH has undergone a long-term poleward shift since 1940, highlighted by a pronounced southeastward curvature of the geopotential contours (80-110oW). Patagonian glaciers experience increased atmospheric pressure (1.01 ± 0.81 m dec−1, p < 0.05) (Fig. 5c) that drives enhanced northwesterly warm subtropical air advection (Fig. 5b). This triggers significant atmospheric warming near the Southern Andes (0.14 ± 0.03 oC dec−1, p < 0.01) (Fig. 5d), that is 17% larger than the global average (0.12 ± 0.01 oC dec−1, p < 0.01). We conclude that this atmospheric warming is primarily responsible for the long-term runoff increase (r = 0.83) (Fig. 5d).

Methods
Surface mass balance and mass change records
For SMB evaluation, we use 74 stake measurements spanning the period 1980-2019, and collected at 38 locations across NPI (2 sites), SPI (26 sites) and CDI (19 sites) (inset maps in Fig. 1a). SMB records were estimated from stakes and firn cores drilled on: i) San Rafael (1984)37 and Nef (1996)38 glaciers in NPI; ii) Perito Moreno (1980–198539 and 1996-200140,41), De Los Tres (1995-1996)42, Chico (1994–2001)43, Cerro Gorra Blanca (1995–2001)44,45, Tyndall (1998–2000)46, and Pio XI (2000–2006)45 glaciers in SPI; and iii) Monte Sarmiento Massif (2013-2019)47 in CDI. We discarded data compiling less than 3 months of monitoring. For a meaningful comparison, we cumulated downscaled SMB at 500 m from MAR and RACMO over the overlapping months of measurements (Supplementary Fig. 5).
We compare monthly GRACE and GRACE-FO mass change records (and uncertainties)24 for the period 2002-2023 with input-output method estimates (MB = SMB-D) (Fig. 1a and Supplementary Fig. 4). Annual mean solid ice discharge (D) was set to 15.1 ± 1.1 Gt yr−1 in 1940-19991 and 24.1 ± 1.7 Gt yr−1 in 2000–202312,48, assuming that these long-term estimates for 1968–1999 and 2000–2019 remain valid pre-1968 and post-2019 respectively, and were equally distributed within each month of the year. In future work, these solid ice discharge estimates could be improved by exploiting historical imagery to cover the pre-satellite era49, hence capturing potentially large long-term variability. Monthly SMB is derived from MAR and RACMO spatially integrated over glacier-covered grid-cells, i.e., ignoring seasonal snow and melt outside glaciers. While this may affect the amplitude of mass change, i.e., with peak snow/melt in winter/summer, it does not affect its trend since snow cannot accumulate outside glacier areas over the years.

Gridded meteorological records
We use two gridded meteorological datasets at 5 km, i.e., CR2MET (1960-2021)32 and PMET (1980-2020)31, to evaluate downscaled glacier-integrated total precipitation and snowfall from MAR and RACMO at 500 m spatial resolution (Supplementary Fig. 7). While the meteorological datasets correlate well in the overlapping period 1980-2020 (R2 = 0.71), total precipitation is on average 32% (55.8 Gt yr−1) larger in PMET than in CR2MET. Total precipitation in MAR and RACMO show excellent agreement with PMET (R2 = 0.86 and 0.91; RMSE = 8.7 Gt yr−1 and 14.0 Gt yr−1 respectively), while being 44% to 37% larger than in CR2MET (51.6 Gt yr−1 and 45.4 Gt yr−1 respectively). We attribute the difference between the two meteorological datasets to the fact that CR2MET may only capture solid precipitation. This is supported by the relatively low RMSE obtained when comparing downscaled snowfall from MAR (15.3 Gt yr−1) and RACMO (10.6 Gt yr−1) to the CR2MET data set, i.e., 13% and 9% of the total precipitation in CR2MET, respectively.
We also compare anomalies in annual mean near-surface air temperature (T2m), spatially averaged over glacier-covered grid-cells, from MAR and RACMO with those of CR2MET (1960-2021) and PMET (1980-2020). For a meaningful comparison, anomalies are estimated with respect to the period 1960–1989, or the overlapping period in each data set, i.e., 1979–1989 for RACMO and 1980–1989 for PMET (Fig. 1c). All four datasets show excellent agreement with 0.74 < R2 < 0.89, and 0.14 °C < RMSE < 0.28 °C.

MAR: modèle atmosphérique régional
The Modèle Atmosphérique Régional version 3.14 (MARv3.14) implements the dynamical core developed by Gallée and Schayes (1994)50 and the physics discussed by Fettweis et al. (2005)51. MAR has been used and thoroughly evaluated for glaciated areas, including the Greenland ice sheet52, glaciers and ice caps of the Arctic region53, and the Antarctic ice sheet54. The model incorporates a dedicated Soil Ice Snow Vegetation Atmosphere Transfer (SISVAT) module55 that is specifically adapted for snow and ice processes56,57. Surface melt, percolation and retention in firn, and subsequent surface runoff are simulated in a 21-layer snowpack. Snow albedo is computed in the CROCUS sub-module58 based on snow grain properties (size, sphericity, dendricity), snow densification, solar zenith angle, cloud optical thickness, and surface meltwater ponding. Here, MAR is run at 5 km spatial resolution and forced at its lateral (seven pixels) and upper (stratosphere) atmospheric boundaries by the latest ERA5 climate reanalysis27 (1940-2023). Forcing fields include temperature, pressure, specific humidity, wind speed and direction prescribed at the 24 model atmospheric levels. Sea surface temperature is prescribed by ERA5 reanalysis on a 6-hourly basis. Firn is initialised with a former MARv3.13 10-year spin-up simulation at 10 km spatial resolution using initial snowpack densities of 500 kg m−3 and 920 kg m−3 above and below the equilibrium line at ~1200 m a.s.l., respectively. Before conducting the 1940-2023 MARv3.14 simulation at 5 km, an additional 4-year spin-up is carried out with snowpack initialisation from the former MARv3.13 run at 10 km. Ice albedo is fixed at 0.55 as a constant in space and time. Surface topography and ice mask are fixed and derived from the 1’ resolution digital elevation model ETOPO0159 and the European Space Agency (ESA) Climate Change Initiative (CCI) Land Cover User Tool version 3.1060 at 1 km resolution, both down-sampled to 5 km.
Compared to the high-resolution glacier outlines from the Randolph Glacier Inventory version 6 (RGI 6.0)61, MAR at 5 km generally overestimates ice extent in low-lying ablation areas, notably in SPI (Supplementary Fig. 1a). Surface topography in glacier areas is overall 609 m lower in MAR than in the high resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM)30 (Supplementary Fig. 2b), with small outlet glaciers being generally too high, and mountain divides and promontories being too low. As a result, glacier hypsometry, i.e., the area-elevation distribution (Supplementary Fig. 9a), is overestimated at low elevations in MAR (green line) relative to SRTM DEM (black line), notably in ablation areas below 1200 m a.s.l., and vice-versa further inland (Supplementary Fig. 9b). Using 74 in-situ SMB measurements for model evaluation, we find that MAR overestimates SMB across low-lying outlet glaciers, and underestimates SMB in inland accumulation zones, where erroneous ablation conditions are captured (Supplementary Fig. 5a). Poor agreement with in-situ observations (R2 = 0.05) and large RMSE (6.1 m w.e) indicate that MAR typically underestimate SMB across Patagonian glaciers and ice fields. This is supported by an overall mass loss overestimate when comparing MAR-derived mass change (SMB-D) with GRACE/GRACE-FO records (RMSE = 606 Gt) (green and red lines in Supplementary Fig. 4a).

RACMO: regional atmospheric climate model
The Regional Atmospheric Climate Model version 2.3p2 (RACMO2.3p2) incorporates the dynamical core of the High Resolution Limited Area Model (HIRLAM)62 and the physics package cycle CY33r1 of the European Centre for Medium-Range Weather Forecasts-Integrated Forecast System (ECMWF-IFS)63. The model is specifically adapted to represent surface processes of polar ice sheets and ice caps, including the Greenland ice sheet26, Canadian Arctic64, Svalbard65, Iceland66, Patagonia17 and Antarctica67. The model has a 40-layer snow module simulating melt, percolation and retention into firn and subsequent surface runoff68. The model represents dry-snow densification69, drifting snow erosion70, and snow albedo based on grain size, cloud optical thickness, solar zenith angle, and impurity content71. RACMO2.3p2 at 5.5 km spatial resolution is forced by ERA5 reanalysis27 (1979–2023) within a 24-grid-cell-wide relaxation zone at the lateral model boundaries. Forcing consists of temperature, pressure, specific humidity, wind speed and direction being prescribed at the 40 model atmospheric levels every 3 hours. Upper atmospheric relaxation is active72. Sea surface temperature is prescribed from the ERA5 reanalysis on a 3-hourly basis. Firn is initialised on 1st January 1979 by prescribing an initial snow depth (4 m), temperature (-10 °C), and density profile (300–900 kg m−3), followed by a 5-year spin-up simulation (1979-1983). The spin-up ensures that the snowpack properties reach equilibrium with atmospheric conditions. The presented 1979–2023 simulation is then branched from the 5-year spin-up. Ice albedo is prescribed as a constant field in space and time at 0.55. Surface topography and ice mask are fixed and derived from the 30” resolution digital elevation model GTOPO3073 and the Global Land Cover Characteristics63 at 1 km resolution, both down-sampled to 5.5 km.
RACMO at 5.5 km spatial resolution resolves the three major ice-fields relatively well (Supplementary Fig. 1a), but does not capture the smaller neighbouring mountain glaciers outlined in the high-resolution RGI 6.0 product61. Relative to the SRTM DEM30, RACMO shows patterns of surface elevation bias similar to those of MAR, though being smaller on average (529 m) (Supplementary Fig. 2c). As opposed to MAR, glacier hypsometry in RACMO at 5.5 km is underestimated below 600 m a.s.l. (Supplementary Fig. 9d), as low-lying glaciers are generally not captured. This is supported by model evaluation using 74 in-situ measurements (Supplementary Fig. 5c), showing an overall SMB overestimate in the ablation zone. RACMO aligns better with in-situ observations (R2 = 0.55) than MAR (R2 = 0.05), although with similar RMSE (6.3 m w.e.) and an over three-fold larger positive bias (3 m w.e.). Comparison to GRACE/GRACE-FO mass loss records confirms a SMB overestimate in RACMO, with an erroneous and persistent mass gain since 2002 (blue line in Supplementary Fig. 4a), in line with previous studies17,20. In addition, RACMO-derived mass change RMSE (1227 Gt) is two-fold larger than that of MAR (606 Gt) (Supplementary Fig. 4a).

Statistical downscaling
MAR (1940-2023) and RACMO (1979–2023) SMB components are statistically downscaled from the native model resolution of 5 km and 5.5 km, respectively, to a 500 m ice mask and topography derived from RGI6.0 (1950–2010)61 and the high-resolution SRTM DEM version 4 (2000)30. The downscaling procedure corrects individual SMB components for elevation on the 500 m topography, using daily-specific vertical gradients estimated on the native model grids. SMB components (X in Eq. (1)) include total precipitation (PR), total sublimation (SU), total melt (ME), and runoff (RU). Drifting snow erosion (ER) is accounted for in RACMO. Vertical gradients are estimated as linear regressions using at least six grid-cells, i.e., the current one and at least five (up to eight) adjacent pixels. To obtain realistic local estimates, the regression slope (a) is applied to the current grid-cell to compute an intercept (b), i.e., value at sea level. These two regression coefficients are bi-linearly interpolated from the low-resolution model grids onto the high resolution one, and applied to the SRTM DEM surface topography at 500 m (h) as,
X 500m = a 500m x h 500m + b 500m (1)
Melt and runoff are further corrected for surface albedo in regions exposing dark bare ice (albedo < 0.55) in summer, which are unresolved in MAR and RACMO. To that end, we use a 500 m MODIS 16-day product averaged for the period 2000–2023. MODIS bare ice albedo is estimated as the mean of the 5% lowest surface albedo recorded each year, averaged for the period 2000-2023. In addition, daily snowfall fraction (SFfrac), i.e., the fraction of snowfall on total precipitation, is statistically downscaled to 500 m following Huai et al. (2022)74. Snowfall is thus estimated as,
SF 500m = PR 500m x SF frac 500m (2)
Rainfall is estimated as a residual as,
RA 500m = PR 500m - SF 500m (3)
SMB is reconstructed using individual components statistically downscaled to 500 m as,
SMB 500m = PR 500m - RU 500m - SU 500m - ER 500m (4)
Note that the drifting snow erosion flux (ER) is not accounted for in MAR. Refreezing (RF) is estimated as a residual following,
RF 500m = ME 500m + RA 500m - RU 500m (5)

Spatial refinement and adjustment
In MAR, statistical downscaling results in a 6% increase in total precipitation (9.4 Gt yr−1 for 1940-2023) combined with a 19% decrease in surface runoff (39.4 Gt yr−1), in turn enhancing SMB by 48.0 Gt yr−1 (Supplementary Table 1). Total precipitation mostly increases around the hypsometry peak of Patagonian glaciers situated between 1000-1600 m a.s.l. (Supplementary Fig. 9a–c), where MAR at 5 km substantially underestimates surface elevation and glacier area. This is in excellent agreement with PMET-derived vertical precipitation profile (cyan line in Supplementary Fig. 9c). Surface runoff decreases below the hypsometry peak, where MAR at 5 km largely overestimates glacier area, and increases to a smaller extent inland above the hypsometry peak, where glacier area is underestimated (Supplementary Fig. 9a–c). Compared to in-situ SMB records, statistical downscaling significantly improves upon the native MAR product at 5 km, with higher correlation (R2 = 0.83) and essentially halved RMSE (3.3 m w.e.) (Supplementary Fig. 5a–b). Statistical downscaling also improves MAR-derived mass change agreement with GRACE/GRACE-FO, i.e., the RMSE is almost eliminated (36 Gt) and reduced by 94% relative to the native product at 5 km (green lines in Supplementary Fig. 4a, b). In view of the good agreement with in-situ measurements, PMET precipitation accumulation, and GRACE/GRACE-FO mass change, we deem that no further adjustment is required for downscaled MAR.
In RACMO, runoff increases by 92% (79.7 Gt yr−1 for 1979-2023) through statistical downscaling. This is almost balanced by a 40% increase in total precipitation (72.3 Gt yr−1), hence only reducing SMB by 7% (6.8 Gt yr−1) relative to the native product at 5.5 km. RACMO generally underestimates glacier area at all elevations (Supplementary Fig. 9d), notably below and around the Patagonian glacier hypsometry peak (1000-1600 m a.s.l.). The largest increases in total precipitation and runoff are thus found in low-lying (respectively elevated) regions where outlet (respectively mountain) glaciers were not resolved at 5.5 km (Supplementary Fig. 9d, e). As mentioned in previous studies17,20, RACMO at 5.5 km exaggerates orographic-forced precipitation across the Andes. This process is amplified through statistical downscaling as demonstrated by large precipitation overestimates relative to the PMET dataset (cyan line in Supplementary Fig. 9f). This is supported by large SMB overestimate compared to in-situ measurements in the accumulation zone (Supplementary Fig. 5d). Nonetheless, we find that statistical downscaling improves upon the native product at 5.5 km, with higher correlation (R2 = 0.79) and lower RMSE (5.3 m w.e.) (Supplementary Fig. 5c, d). Comparing downscaled RACMO-derived mass change with GRACE/GRACE-FO shows similar positive bias, correlation and RMSE to the native product at 5.5 km (Supplementary Fig. 4a, b), the result of almost unchanged SMB as precipitation increase compensates for enhanced runoff. In line with remote sensing, the amplitude of seasonal mass change has increased, i.e., with higher winter accumulation and summertime ablation. Reducing total precipitation by 36.5% (92.0 Gt yr−1 in Supplementary Table 1) across all glacier grid-cells on a daily basis eliminates the remaining positive bias and decreases mass change RMSE to 31 Gt relative to remote sensing (blue line in Supplementary Fig. 4c). Adjusted precipitation also results in improved vertical profile agreement between downscaled RACMO and PMET products (cyan line in Supplementary Fig. 9g), better aligns with glacier-integrated MAR and PMET precipitation estimates (Supplementary Fig. 7a), and further decreases the SMB RMSE to 4.9 m w.e. compared to in-situ records (Supplementary Fig. 5e).

Spatial integration
To estimate mass fluxes in Gt, we spatially integrate SMB and its individual components over ice-covered grid-cells. To that end, we use RGI6.061 glacier outlines to derive a fractional ice mask with grid-cell values ranging from 0 (ice-free) to 1 (fully ice-covered). The use of a fractional mask (0 to 1 in Supplementary Fig. 10a) instead of a binary one (0 or 1 in Supplementary Fig. 10b) ensures that large glacier area overestimate is prevented, notably for small mountain glaciers covering less than 0.25 km2 (i.e., area of a 500 m grid-cell). Here, we estimate glacier-integrated SMB and individual components (M in Gt) as,
M = sum(m i,j x A i,j x I i,j) x 1e-6 (6)
where mi,j is the SMB (or individual components) expressed in kg m−2 (or mm w.e.) and estimated for grid-cell i,j; Ai,j is the area of grid-cell i,j set to 0.25 km2; Ii,j is the fractional ice mask of grid-cell i,j.

Uncertainty estimates
For glacier-integrated SMB (components), uncertainties are estimated as the annual model difference (MAR minus RACMO) in glacier-integrated SMB (components), statistically downscaled to 500 m and averaged for the overlapping period 1979–2023. Here RACMO estimates include the adjustment in total precipitation (-36.5%). These uncertainties are shown in Supplementary Fig. 6 (Diff) and listed in Supplementary Table 1 (Difference).
For glacier mass change (MB = SMB-D), we sum the uncertainty in D estimated at 1.7 Gt yr−1 in Minowa et al. (2021)12, with that of SMB estimated as the RMSE of the model difference (MAR minus RACMO; 7.3 Gt yr−1 in Supplementary Fig. 6a), further halved to 3.7 Gt yr−1 to yield an equally distributed positive and negative uncertainty around our SMB estimates. We thus obtain a model mass change uncertainty of 5.4 Gt yr−1 as listed in Supplementary Table 2. For cumulative mass change and contribution to sea-level rise, we cumulate monthly mass change uncertainty in time (0.45 Gt per month), from January 1940 to December 2023 (Fig. 1a). Sea-level rise uncertainties are derived from the latter monthly cumulative mass change uncertainties, assuming that 362.5 Gt of ice raises global sea level by 1 mm, as in Cogley et al. (2011)75.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a very interesting manuscript aiming to resolve some of the crucial questions we still have about Patagonian glaciers.\nWe have known for a while that glaciers in the region are retreating and shrinking in recent decades, process modulated by\ncalving but utterly triggered by climatic changes. However, these climatic changes are still poorly known due to the lack of in\nsitu measured long term climatic series, the availability of few stations closer to the main glaciers and the difficulties of\nmeasuring mass balance and other variables on the ground. In this sense, most of our understanding about these climatic\nchanges are based upon gridded low resolution data that without proper comparison with field data for validation purposes\nare having high uncertainties. Recent papers using gridded data are indicating some controversial results, since they\nobtained positive surface mass balances, but, especially the calving glaciers of the main icefields are having very strong\ntotal mass losses due to dynamic responses, indicating that something is not fitting very well and require further research.\nAre these glacier dynamic changes responses having large delays to climate changes? If the surface mass balance is\npositive, even slightly positive, are these dynamic retreats going to stop soon? The work presented by Noël et al addresses\nsome of these issues, therefore, it is very valuable. Despite the above, after reading the text carefully I found some points that\ndeserve to be clarified or resolved to fully accept their findings. Below some of my main points:\nOne first question is the study period. The authors claimed to have reconstructed surface mass balance between 1940 and\n2023. That’s certainly new and useful, but I’m wondering how they estimated D (solid ice discharge) from periods when no\nmaps with accurate elevation data are available. They quoted Rignot, Rivera an Casassa 2003 for D values between 1940\nand 1999, however, in this study were used the older available, accurate enough, elevation data from 1968 (Argentina) and\n1975 (Chile). Both were compared to SRTM (2000). These old DEMs are not very accurate but in front of the strong\nelevation changes, the delta H results were significant. From where they obtained D, needed for calculating MB in the period\npre-1968 in Argentina and pre-1975 in Chile? The showed uncertainties are also very small, almost invisible in this series,\nespecially at the beginning. That’s surprising, since older data with no validation neither control is expected to have higher\nuncertainties. Then, Minowa et al is quoted for D between 2000 and 2023, but their work was comparing 2000 and 2019.\nWhat you did to obtain D between 2019 and 2023?\nIn line 17 it is quoted Aniya et al., 1996 for the SPI, that’s fine as an historical and important reference, but there are very\nstrong changes between his work and more recent estimates. I suggest using more recent inventories for glacier areas, for\nexample De Angelis 2014 https://doi.org/10.3189/2014JoG13J127\nIn line 18 it is quoted Temme et al 2023 for the total area of Cordillera Darwin, however, this paper is about Monte Sarmiento\nand doesn’t estimate total areas. Maybe is better to quote the inventory done by Bown et al., 2014\nhttps://link.springer.com/chapter/10.1007/978-3-540-79818-7_28\nIn line 38 they said “SPI floating glaciers”. As far as we know, no calving Patagonian glaciers are floating (only occasionally\nor ephemeral conditions as described by Minowa et al 2023 https://doi.org/10.1016/j.epsl.2021.116811). Some of these\nglaciers are close to flotation and buoyancy forces explain strong calving, but there is not floating ice since as soon as these\ntemperate glaciers approach flotation (for example in deep fjords) the fronts break away.\nIn line 40 they quoted Minowa et al when stating that since the 2000s the SPI ice discharge experiences a decreasing trend,\nhowever, this trend is missing if you exclude Pio XI. In his paper, Minowa et al realized that frontal ablation of the whole SPI\nwithout Pio XI is relatively stable (see his Figure 3a green line). This is explained by the almost continuous advance of Pio\nXI glacier, the biggest single basin in the SPI, and not necessarily by a global positive SMB. In this sense, I’m not totally sure\nthat the positive surface mass balance is explaining the reduction of the frontal ablation of the whole SPI. Another big issue\nis why Pio XI is advancing, and SMB is probably an important factor, but not the only one, since reduced ablation at the front\nby sedimentation is also playing a role (Muto et al., https://doi.org/10.3389/feart.2020.576044).\nLine 46-47 the SMB uncertainties are high because the used climatic models have low resolution, but at least in some of the\nmodels (ref 14, 19 and 20), the data were downscaled to 450 m. Why your model downscaling to 500 m will decrease the\nuncertainties, I mean, from this point of view?\nLine 64-65. The correction of 12.5 and 32.5% of MAR and RACMO is the most critical point in the text. Without these\ncorrections, RACMO shows a strong positive trend. After reading the methods I think is better to add to the main text a better\nexplanation of these corrections. By the way, correcting by a constant value in every cell in the domain is quite arguable. I\nunderstand that the final result compares very well with GRACE/GRACE-FO, but the reasons explained for the under or over\nestimations in these models are not homogenously distributed on space. A better justification is needed.\nLine 166-167. The ablation zones have expanded as you said, but the glaciers as a whole have shrinkage significantly,\nmainly due to the ablation zones/fronts calving away. Did you account for this area reductions? Which area are you using for\nthis calculation? Only from 1986 to year 2016 the Patagonia glaciers lost an area of near 2200 km2 (Meier et al., 2018\nhttps://doi.org/10.3389/feart.2018.00062), almost all within ablation areas. What about the ELA elevation changes and/or the\nAAR? A better understanding of the ablation zone expansion as you concluded, needs a better discussion about total area\nreductions, ELA and AAR changes. I suggest to address these concerns when discussing your results.\nLine 248-251. This asseveration is quite strong as it is. Assuming no significant fluctuation in frontal ablation is a BIG\nassumption. We know that calving glacier are dynamically responding with very anomalous responses. It is just a matter of\nseeing Pio XI’s continuous advances. I rather prefer to see a discussion about dynamic readjustments or delayed\nresponses, before predicting the disappearing of the huge volume of ice in the icefields.\nFigure 1. Maybe you can add a reference for the conversion between mass change and 1 mm of sea level? I don’t see the\nRACMO coloured bands in b. After increasing the ZOOM to a very maximum, I only got a grey colour underneath the SMB\nblack line. Is this the band? Invisible for to me…Maybe leaving only MAR in here and showing RACMO in the\nsupplementary info? Add a line to legend indicating what are the stars (sites with ¬mass balance measurements). I know\nyou said so in the text, but I think is better to see the figure with all the info in the legend.\nFigure 3. I understand that only yellow polygon derived from MAR was used here but why not other variable with high r?.\nWhy 3c was obtained with N=84 and not 129 as in most of the other correlations in this and other figures (e.g. figure 4)? It is\nbecause T2m is averaged over glacier areas? If so, why in 3b N=129 and not 84? Please explain the “averaged over\nglacier” approach.\n(In the high resolution version of Figure 3 (page 27 of the merged file) there are many purple and yellow points…I\npresume these are just printing errors, or have any meaning? I don’t see these “anomalous” points in Figure 3 in page 22. The same\nis happening in Figure 4 in page 28)\nFigure 5. I don’t see the RACMO bands in d. The yellow polygon is the MAR r>0.75 correlation area shown in Supple 8a,\nright? Why using this polygon and not the RACMO derived >0.75? or the one correlating Z850 and pp? Seems to be\narbitrary the selection of this specific one. Please add a phrase explaining why this and not any other one was used.\nSupplementary Figure 4. In the legend the time frame is different from the time scale of these figures. I suggest explaining\nboth (one is the time used for the whole series but here only using data from 2002 and 2023 for comparison purposes)\nSupplementary 6. Why N=45? I thought that you can compared the whole grid (N=129).\nSupplementary Figure 7: I don’t see RACMO coloured bands either. Are there? The alignment between total pp of MAR and\nPMET is clear, and also between MAR snowfall and CR2Met, but I don’t see the rest of the relationships indicated in the\nlegend. Maybe you can add the R2 between them? Using max fluxes in primary “Y” is a bit confusing as you are only\ndiscussing precipitations in this figure and I think m we/yr is enough. Maybe saying something about the conversion from m\nwe/yr to GT/yr? or moving GT/yr to the secondary Y?\nSupplementary 8. The yellow line for RACMO?\nFormat details:\nRACMO and RACMO2 are used alternatively, maybe use just one acronyms to avoid misunderstandings.\nIn Line 444 you said “cyan in Supplementary Fig.9e”… but I think you are referring to Figure 9h.", "reviewer_2": "This is a comprehensive and well-executed study on the surface mass balance (SMB) and long-term mass loss of\nPatagonian glaciers, which are key contributors to global sea-level rise. The integration of high-resolution regional climate\nmodels and statistical downscaling provides robust framework for analyzing SMB over an extended period (1940–2023).\nThe authors effectively combine these results with GRACE/GRACE-FO satellite data and solid ice discharge estimates,\noffering valuable insights into the drivers of glacier mass loss. The study is particularly notable for identifying the\npoleward shift of subtropical highs and its implications for warm air advection, demonstrating that long-term fluctuations in\nsurface runoff are a more significant driver of SMB trends than precipitation.\nThe paper is very clearly written and well-presented, especially considering the complexity of the analysis, and I see no\nerrors or omissions in the detail of the paper. I only have minor suggestions exclusively pertaining to figures:\n-Size may be an issue, but it would be helpful if Fig 1a had a scalebar\n-Fig 2 is also missing a scalebar. What is the background hillshade image the data is overlaid on?\n-Several of the plots are not colorblind compatible, particularly with regards to red/green colorblindness in Fig 1a, 1c, 2a and\n2b."}, "author_response": "We would like to thank the Editor for securing two reviews, and for inviting us to revise our manuscript. Below you\nwill find in blue our point-by-point responses to the Editor’s and Reviewers’ comment, while changes applied to the\nrevised manuscript are displayed in red. The manuscript has been thoroughly and substantially revised based on the\nEditor and Reviewers’ comments. Major changes include the followings:\n1. Based on a constructive comment from Reviewer #1 on glacier retreat, we realized that our integration ice mask\nwas overestimated for small mountain glaciers in the initial submission. We addressed this in the revised\nmanuscript by creating a fractional ice mask ensuring that glacier area is not overestimated anymore.\n2. We re-conducted statistical downscaling for both MAR and RACMO using this new fractional ice mask. As a result,\nwe found that MAR does not require runoff adjustment anymore to align with GRACE/GRACE-FO mass change;\nRACMO still requires total precipitation adjustment, but the correction has decreased relative to the initial\nsubmission. These technical changes did not affect the main findings of our manuscript.\n3. We clarified in a new Methods Section “Product strengths and limitations” that our regional climate models do\nnot account for glacier retreat and thinning as they are not coupled to a dedicated ice dynamics module. Such\ncoupling has, to our knowledge, never been achieved in the literature for small Patagonian glaciers and is beyond\nthe scope of our study. As an alternative, we quantified the SMB response to glacier retreat by comparing glacier-\nintegrated SMB (components) over glacier outlines of the Randolph Glacier Inventory version 6 (RGI 6.0, 1950-2010)\nwith the latest version 7 (RGI 7.0, 1950-2021), i.e., smaller by roughly 1,500 km2. We find that changes\nremain statistically insignificant. This comparison was only made possible through the creation of fractional ice\nmasks for both RGI 6.0 and 7.0.\n4. We now elaborate on the strengths of our study relative to previous work in a new Methods Section “Product\nstrengths and limitations”. Strengths include i) the use of polar regional climate models combined with a\nsophisticated statistical downscaling technique, ii) thorough model evaluation using a comprehensive set of\nindependent observational measurements, i.e., now including the latest GlaMBIE assessment.\n5. We upgraded our GRACE/GRACE-FO timeseries with an improved new release extending to December 2023. We\nalso included mass change estimates from the latest GlaMBIE assessment (2000-2023). We find excellent model\nagreement with this new and independent data set (orange dots in Fig. 1a, Supplementary Table 2).\nThese corrections, additional data, and clarifications improved the robustness of our results, without altering our\nmain messages: i) the long-term trend in Patagonian glacier mass loss is dominated by increased surface runoff; ii) a\npoleward shift in subtropical highs primarily drives runoff increase in Patagonia.\nAs mentioned by reviewer #1, we acknowledge that observational data sets in Patagonia are sparse and come with\nlarge uncertainties, but we deem that this should not be a constraint to improve our knowledge on glacier mass\nchange and its driving processes. Here, we made substantial efforts to compile a comprehensive set of independent\nobservational measurements, including the latest GRACE/GRACE-FO records and GlaMBIE assessment, enabling a\nthorough model evaluation that ensures solid results. We also present two independent high-resolution data sets\nderived from two state-of-the-art regional climate models, statistically downscaled to 500 m. These products not\nonly show excellent agreement with each other (Supplementary Fig. 6), but also with the extensive observational\ndata sets we compiled. To our knowledge, such detailed model evaluation has never been achieved in the literature,\nand the generally good model agreement demonstrates the robustness and strength of our downscaled products.\nThis is clarified in the new Methods Section “Product strengths and limitations” as: “The resulting statistically\ndownscaled products were thoroughly evaluated with a comprehensive set of independent observations, including\ni) averaged solid ice discharge estimates for 1968-1999 and 2000-2019, ii) previously published geodetic, glaciological\nand remote sensing mass change measurements, complemented with GRACE/GRACE-FO records and the latest\nGlaMBIE assessment (Supplementary Table 2), iii) point in situ measurements spanning the three main icefields over\nthe period 1980-2019 (Supplementary Fig. 5), iv) gridded precipitation and temperature observations from two\nmeteorological products CR2MET (1960-2021) and PMET (1980-2020) (Supplementary Fig. 7). Such detailed\nevaluation has never been achieved in the literature and yielded generally good model agreement, demonstrating\nthe robustness and strength of our statistical downscaling approach and associated high-resolution products.”\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nAuthors’ general comment. We would like to thank Reviewer #1 for their interest in our study, and for their thorough\nand constructive comments. Based on your remarks on glacier retreat, we noticed that our two downscaled data sets\nsuffered from an ice mask extent overestimate. In the original submission, our ice mask was binary (0 or 1 in\nSupplementary Fig. 10b), i.e., it did not account for the ice fraction within each grid-cell (i.e., ranging from 0 to 1 in\nSupplementary Fig. 10a). Using a binary mask proved to largely overestimate the total area of small mountain glaciers\noutside the main icefields (i.e., ~20.000 km2 previously relative to ~10.000 km2 now). This overestimate was less\npronounced for the more extensive icefields, and we therefore did not notice this problem.\nSupplementary Figure 10: Fractional and binary ice masks. a Fractional ice mask, i.e., ranging from 0 (ice free grid-\ncell) to 1 (fully ice-covered grid-cell), and b binary ice mask, i.e., set to 0 (ice free) or 1 (ice-covered), of Patagonian\nglaciers at 500 m spatial resolution. The black outline shows the RGI 6.0 glacier outline. The inset map locates the\nselected area displayed in a and b (red contour). The background image of each subfigure is derived from Google\nEarth.\nIn view of this, we relaunched our statistical downscaling technique for both MAR and RACMO. We figured out that\nthe 12.5% runoff overestimate in downscaled MAR in the original submission largely resulted from too large\nmountain glaciers area in our binary mask (Supplementary Fig. 10b). Applying our new fractional mask\n(Supplementary Fig. 10a) to downscale MAR almost eliminates the runoff overestimate previously found relative to\nGRACE/GRACE-FO mass change records. As opposed to our initial submission, here, downscaled MAR no more\nrequires runoff adjustment when using a fractional ice mask. For downscaled RACMO, the use of a fractional mask\ndoes not eliminate the precipitation overestimate found in the original submission, notably because peak\nprecipitation overestimates occur inland across steep Andes divides that are generally not affected by the fractional\nice mask. Therefore, a precipitation adjustment (36.5%) is still required to eliminate the positive bias compared to\nGRACE/GRACE-FO. This adjustment has slightly decreased in the updated RACMO downscaled product, i.e., from a\nprecipitation reduction of 98.4 Gt yr-1 to 92.0 Gt yr-1. We find that overall, the use of a fractional mask has improved\nthe robustness of our results, and did not affect the main messages of our study: 1) the long-term trend in Patagonian\nglacier mass loss is driven by increased surface runoff (r = 0.80 in Fig. 3a), rather than precipitation changes (r = 0.61\nin Fig. 3b); 2) a poleward shift in subtropical highs primarily drives the latter runoff increase in Patagonia. We have\ncorrected all Figures, numbers and Tables across the revised manuscript and supplements accordingly. We now clarify\nhow integrated fluxes in Gt were estimated using our new fractional mask in the revised Methods Section “Spatial\nintegration” as: “To estimate mass fluxes in Gt, we spatially integrate SMB and its individual components over ice-\ncovered grid-cells. To that end, we use RGI 6.0 glacier outlines to derive a fractional ice mask with grid-cell values\nranging from 0 (ice free) to 1 (fully ice-covered). The use of a fractional mask (0 to 1 in Supplementary Fig. 10a)\ninstead of a binary one (0 or 1 in Supplementary Fig. 10b) ensures that large glacier area overestimate is prevented,\nnotably for small mountain glaciers covering less than 0.25 km2 (i.e., area of a 500 m grid-cell). Here we estimate\nglacier-integrated SMB and individual components (M in Gt) as,\nM = ∑ m * A * I * 1e-6\ni,j i,j i,j\nwhere m is the SMB (or individual components) expressed in kg m-2 (or mm w.e.) and estimated for grid-cell i,j; A\ni,j i,j\nis the area of grid-cell i,j set to 0.25 km2; I is the fractional ice mask of grid-cell i,j.”\ni,j\nThis is a very interesting manuscript aiming to resolve some of the crucial questions we still have about Patagonian\nglaciers. We have known for a while that glaciers in the region are retreating and shrinking in recent decades, process\nmodulated by calving but utterly triggered by climatic changes. However, these climatic changes are still poorly\nknown due to the lack of in situ measured long term climatic series, the availability of few stations closer to the main\nglaciers and the difficulties of measuring mass balance and other variables on the ground. In this sense, most of our\nunderstanding about these climatic changes are based upon gridded low resolution data that without proper\ncomparison with field data for validation purposes are having high uncertainties.\nThank you, we acknowledge that observational data sets in Patagonia are sparse and come with large uncertainties,\nbut we deem that this should not be a constraint to improve our knowledge on glacier mass change and its driving\nprocesses. Here, we made substantial efforts to compile a comprehensive set of independent observational\nmeasurements, enabling a thorough model evaluation that ensures solid results. To our knowledge, such detailed\nmodel evaluation has never been achieved in the literature, and the generally good model agreement demonstrates\nthe robustness and strength of our downscaled products. This is now clarified in a new Methods Section “Product\nstrengths and limitations”. Please, see our response to the Editor for additional detail.\nNote that our analysis now includes a new and improved GRACE/GRACE-FO release extending to December 2023, as\nwell mass change estimates from the latest GlaMBIE assessment (2000-2023) (orange dots in Fig. 1a, Supplementary\nTable 2). This is now stressed in the text as: “Modelled mass loss from MAR (28.8 ± 5.4 Gt yr-1) and RACMO (29.9 ±\n5.4 Gt yr-1) aligns well with GRACE/GRACE-FO satellite measurements (27.8 ± 11.0 Gt yr-1) (red line) in the last two\ndecades (2002-2023), showing both a strong correlation (R2 = 0.94) and a low RMSE (~50 Gt) for both models\n(Supplementary Fig. 4c). Furthermore, our estimates are in excellent agreement with previously published mass\nchange records including the latest GlaMBIE assessment (orange dots in Fig. 1a), that reports mass loss of 26.5 ± 6.5\nGt yr-1 in 2000-2023 compared to 27.3 ± 5.4 Gt yr-1 in MAR and 28.5 ± 5.4 Gt yr-1 in RACMO (Supplementary Table 2),\nrespectively.”\nRecent papers using gridded data are indicating some controversial results, since they obtained positive surface mass\nbalances, but, especially the calving glaciers of the main icefields are having very strong total mass losses due to\ndynamic responses, indicating that something is not fitting very well and require further research. Are these glacier\ndynamic changes responses having large delays to climate changes? If the surface mass balance is positive, even\nslightly positive, are these dynamic retreats going to stop soon?\nTo date, state-of-the art regional climate models do not capture glacier dynamics (i.e., glacier retreat and thinning)\nbut instead prescribe a fixed present-day glacier geometry (surface elevation and glacier extent) to reconstruct\ncontemporary SMB time series. Coupling regional climate models to a dedicated ice dynamics model is beyond the\nscope of our study, and, to our knowledge, has never been achieved for the small Patagonian glaciers and icefields\nin the literature. This is now clarified in a new Methods Section “Product strengths and limitations”.\nWe also wish to stress that our study focus is to 1) demonstrate, through detailed model evaluation, that the long-\nterm trend in Patagonian icefields mass loss is driven by increased surface runoff (r = 0.80 in Fig. 3a), rather than\nprecipitation changes (r = 0.61 in Fig. 3b); 2) identify a poleward shift in subtropical highs as main driver for the latter\nrunoff increase in Patagonia. We deem that our detailed analyses and thorough model evaluations provide solid and\nnovel insights on the mechanism driving long-term Patagonian glacier mass loss.\nThe work presented by Noël et al addresses some of these issues, therefore, it is very valuable. Despite the above,\nafter reading the text carefully I found some points that deserve to be clarified or resolved to fully accept their\nfindings. Below some of my main points: Thank you. Please, find below our point-by-point responses to your\ncomments.\n1a) One first question is the study period. The authors claimed to have reconstructed surface mass balance between\n1940 and 2023. That’s certainly new and useful, but I’m wondering how they estimated D (solid ice discharge) from\nperiods when no maps with accurate elevation data are available. They quoted Rignot, Rivera an Casassa 2003 for D\nvalues between 1968-1999, however, in this study were used the older available, accurate enough, elevation\ndata from 1968 (Argentina) and 1975 (Chile). Both were compared to SRTM (2000). These old DEMs are not very\naccurate but in front of the strong elevation changes, the delta H results were significant. From where they obtained\nD, needed for calculating MB in the period pre-1968 in Argentina and pre-1975 in Chile?\nTo clarify, we do not estimate solid ice discharge ourselves, as our regional climate models are not coupled to an ice\ndynamics model. Instead, we use previously published estimates from Rignot et al. (2003) for the period 1968-1999,\nand Minowa et al. (2021) for the period 2000-2019. Since no data are available prior-1968 and post-2019, we assume\nthat the latter two long-term D estimates remain valid for the extended periods 1940-1999 and 2000-2023,\nrespectively. This is now clarified as: “Since published estimates do not fully cover the period 1940-2023, we assume\nthat long-term solid ice discharge averages for 1968-1999 and 2000-2019 remain valid prior 1968 and post 2019,\nrespectively. Therefore, we set solid ice discharge to 15.1 ± 1.1 Gt yr-1 in 1940-1999, with a step increase to 24.1 ±\n1.7 Gt yr-1 in 2000-2023.” We also clarified this in the caption of Fig. 1a and added a vertical dashed line to highlight\nthe transition between both D estimates as: “Mass change is estimated as the difference between downscaled\nsurface mass balance (SMB) from MAR (green line, 1940-2023) and RACMO (blue line, 1979-2023) at 500 m\nresolution, and solid ice discharge estimates from 1968-1999 and 2000-2019 assuming that these remain valid pre-\n1968 (1940-1999) and post-2019 (2000-2023), respectively. The vertical dashed line marks the transition between\nthese two solid ice discharge products.”\n1b) The showed uncertainties are also very small, almost invisible in this series, especially at the beginning. That’s\nsurprising, since older data with no validation neither control is expected to have higher uncertainties.\nThank you, mass change uncertainty is estimated as the sum of both D and SMB uncertainties. D uncertainty is\nestimated as the largest uncertainty value extracted from the literature, i.e., 1.7 Gt yr-1 for 2000-2019 in Minowa et\nal. (2021). In the original submission, SMB uncertainty was estimated as the mean difference between adjusted MAR\nand RACMO downscaled SMB, spatially integrated over the glacier area for the overlapping period 1979-2023 (0.7\nGt yr-1 in our former Supplementary Table 1). We acknowledge that the latter SMB uncertainty was too conservative.\nInstead, we now estimate the RMSE of the model difference (MAR minus RACMO; 7.3 Gt yr-1) for the overlapping\nperiod 1979-2023, yielding a full error range. The RMSE is further halved (3.7 Gt yr-1) to obtain an equally distributed\npositive and negative uncertainty around our SMB estimates. Thus, annual mass change uncertainty reaches 5.4 Gt\nyr-1, or 0.45 Gt per month (i.e., 5.4/12). Monthly mass change uncertainty is cumulated in time until December 2023.\nThis is clarified in the new Methods Section “Uncertainty estimates” as: “For glacier mass change (MB = SMB - D),\nwe sum the uncertainty in D estimated at 1.7 Gt yr-1 in Minowa et al. (2021), with that of SMB estimated as the RMSE\nof the model difference (MAR minus RACMO; 7.3 Gt yr-1 in Supplementary Fig. 6a), further halved to 3.7 Gt yr-1 to\nyield an equally distributed positive and negative uncertainty around our SMB estimates. We thus obtain a model\nmass change uncertainty of 5.4 Gt yr-1 as listed in Supplementary Table 2. For cumulative mass change and\ncontribution to sea-level rise, we cumulate monthly mass change uncertainty in time (0.45 Gt per month), from\nJanuary 1940 to December 2023 (Fig. 1a). Sea-level rise uncertainties are derived from the latter monthly cumulative\nmass change uncertainties, assuming that 362.5 Gt of ice raises global sea level by 1 mm, as in Cogley et al. (2011).”\n1c) Then, Minowa et al is quoted for D between 2000 and 2023, but their work was comparing 2000 and 2019.\nWhat you did to obtain D between 2019 and 2023?\nHere we assumed that averaged D estimate for 2000-2019 (24.1 ± 1.7 Gt yr-1) remains valid post 2019. This is now\nclarified in the revised paper. Please see also our response to comment 1a).\n2) In line 17 it is quoted Aniya et al., 1996 for the SPI, that’s fine as an historical and important reference, but there\nare very strong changes between his work and more recent estimates. I suggest using more recent inventories for\nglacier areas, for example De Angelis 2014 https://doi.org/10.3189/2014JoG13J127\nThank you for pointing this out, we have added the reference accordingly.\n3) In line 18 it is quoted Temme et al 2023 for the total area of Cordillera Darwin, however, this paper is about Monte\nSarmiento and doesn’t estimate total areas. Maybe is better to quote the inventory done by Bown et al.,\n2014 https://link.springer.com/chapter/10.1007/978-3-540-79818-7_28\nThank you again, we have updated the reference accordingly.\n4) In line 38 they said “SPI floating glaciers”. As far as we know, no calving Patagonian glaciers are floating (only\noccasionally or ephemeral conditions as described by Minowa et al 2023\nhttps://doi.org/10.1016/j.epsl.2021.116811). Some of these glaciers are close to flotation and buoyancy forces\nexplain strong calving, but there is not floating ice since as soon as these temperate glaciers approach flotation (for\nexample in deep fjords) the fronts break away.\nThank you, we reformulated as: “SPI glaciers contributed about 80% …”.\n5) In line 40 they quoted Minowa et al when stating that since the 2000s the SPI ice discharge experiences a\ndecreasing trend, however, this trend is missing if you exclude Pio XI. In his paper, Minowa et al realized that frontal\nablation of the whole SPI without Pio XI is relatively stable (see his Figure 3a green line). This is explained by the\nalmost continuous advance of Pio XI glacier, the biggest single basin in the SPI, and not necessarily by a global positive\nSMB. In this sense, I’m not totally sure that the positive surface mass balance is explaining the reduction of the frontal\nablation of the whole SPI. Another big issue is why Pio XI is advancing, and SMB is probably an important factor, but\nnot the only one, since reduced ablation at the front by sedimentation is also playing a role (Muto et\nal., https://doi.org/10.3389/feart.2020.576044).\nHere we examine the whole of SPI when quantifying Patagonian glacier mass change, i.e., including Pio XI South,\nso we deem it relevant to highlight that the trend since the 2000s is negative. However, we now stress that discarding\nPio XI South, solid ice discharge would have remained stable. This does not affect our statement that steady or\ndeclining solid ice discharge in SPI demonstrates the important role of SMB processes in enhancing mass loss. This is\nnow clarified as: “While remaining steady across most SPI glacier calving fronts, solid ice discharge has however\nexperienced an overall decreasing trend since the 2000s, dominated by the continuous advance of a single glacier,\nnamely Pio XI South.”\n6) Line 46-47 the SMB uncertainties are high because the used climatic models have low resolution, but at least in\nsome of the models (ref 14, 19 and 20), the data were downscaled to 450 m. Why your model downscaling to 500 m\nwill decrease the uncertainties, I mean, from this point of view?\nThank you for this very interesting point. Spatial resolution is not the only factor that impacts SMB uncertainty. Other\nfactors matter, including the quality and resolution of the forcing climate reanalysis, the ability of regional climate\nmodels to capture relevant surface mass balance processes (notably firn processes), the assumptions pertaining to\nthe statistical downscaling approach used. For instance, previous studies were based on:\n1. Older, low-resolution climate reanalysis forcing, i.e., NCEP-NCAR at 2.5º (~250 km) and ERA-Interim at 0.75º (~\n75 km) resolution. Here we use the latest ERA5 reanalysis at 31 km spatial resolution to force MAR and RACMO.\n2. Non-polar regional climate models (WRF and RegCM4) that do not account for firn processes, i.e., rain and\nmeltwater retention/refreezing in firn and subsequent runoff, as opposed to MAR and RACMO.\n3. Simplified statistical downscaling approaches based on spatially and temporally fixed meteorological lapse rates.\nIn contrast, we downscale individual SMB components by locally estimating daily specific gradients in MAR and\nRACMO. We also use MODIS satellite records of bare ice albedo at 500 m spatial resolution to correct downscaled\nmelt and runoff in dark ablation zones.\n4. Simplified mass balance models that estimate ablation from downscaled temperature, solar radiation, and\nprescribed snow/ice albedo. In contrast, polar regional climate models like MAR and RACMO i) simulate the\nsurface albedo of snow, ii) quantify melt by explicitly solving the full surface energy budget at the glacier surface,\ni.e., including radiative and turbulent fluxes, iii) capture rain and meltwater retention/refreezing in firn to\naccurately estimate surface runoff.\nThese differences are now briefly introduced in the Introduction as “Previous NPI and SPI SMB reconstructions using\npositive degree day (PDD), regional climate models, sometimes downscaled using fixed lapse rates to force simplified\nmass balance models (Methods), e.g., neglecting rain and meltwater retention in snow, suggested a positive but\ninsignificant SMB trend in the period 1975-2015.”, and discussed in more detail in a new Methods Section “Product\nstrengths and limitations” as: “The SMB of Patagonian icefields was previously estimated on a 450 m grid. In contrast\nto the current study, these previous SMB estimates i) were temporally limited to 1975-2015 (41 years), ii) spatially\nrestricted to NPI and SPI icefields, relied on non-polar regional climate models that iii) were forced by former, low-\nresolution climate reanalyses including NCEP-NCAR (~250 km) and ERA-Interim (~75 km), and iv) do not capture snow\nprocesses such as rain and meltwater retention/refreezing in firn, and subsequent runoff. In addition, these studies\nassumed spatially and temporally fixed lapse rates to statistically downscale modelled meteorological data to 450 m,\nincluding precipitation, temperature and solar radiation. These data at 450 m are then used to force simplified mass\nbalance models that derive accumulation from downscaled precipitation (i.e., neglecting condensation/sublimation),\nand simulate ablation using downscaled temperature, solar radiation and prescribed snow and ice albedo (i.e.,\nneglecting rain and melt retention in firn). In contrast, MAR and RACMO are forced by the latest ERA5 reanalysis (31\nkm), quantify melt by solving the full glacier surface energy budget, i.e., including both radiative and turbulent fluxes,\nand simulate the albedo of snow, firn retention processes and runoff using multi-layer snow modules. Individual\nmodelled SMB components are further statistically downscaled to 500 m using local (spatially varying), daily specific\ngradients. In addition, we use MODIS satellite records at 500 m to correct downscaled melt and runoff fluxes for the\nalbedo of dark bare ice exposed in summer across ablation zone areas. The resulting statistically downscaled products\nwere thoroughly evaluated using a comprehensive set of independent observations, including i) averaged solid ice\ndischarge estimates for 1968-1999 and 2000-2019, ii) previously published geodetic, glaciological and remote sensing\nmass change measurements, complemented with GRACE/GRACE-FO records and the latest GlaMBIE assessment\n(Supplementary Table 2), iii) point in-situ measurements spanning the three main icefields over the period 1980-\n2019 (Supplementary Fig. 5), iv) gridded precipitation and temperature observations from two meteorological\nproducts CR2MET (1960-2021) and PMET (1980-2020) (Supplementary Fig. 7). Such detailed evaluation has never\nbeen achieved in the literature and yielded generally good model agreement, demonstrating the robustness and\nstrength of our statistical downscaling approach and associated high-resolution products.”\n7a) Line 64-65. The correction of 12.5 and 32.5% of MAR and RACMO is the most critical point in the text. Without\nthese corrections, RACMO shows a strong positive trend. After reading the methods I think is better to add to the\nmain text a better explanation of these corrections.\nThank you, we figured out that our ice mask area in the original submission was overestimated, notably over small\nmountain glaciers due to the use of a binary (0 or 1 as soon as a grid-cell is partly ice-covered) instead of a continuous\nfractional (0 to 1) ice mask. This area overestimate largely explained the runoff overestimate in downscaled MAR in\nthe original submission. Using our new fractional ice mask, we find that downscaled MAR does not require further\nadjustment to align with GRACE/GRACE-FO mass change. However, an adjustment in RACMO total precipitation\n(36.5%) is still required after applying our new fractional ice mask. Please see our Authors’ general comment above\nfor additional detail.\n7b) By the way, correcting by a constant value in every cell in the domain is quite arguable. I understand that the\nfinal result compares very well with GRACE/GRACE-FO, but the reasons explained for the under or over estimations in\nthese models are not homogenously distributed on space. A better justification is needed.\nWe fully agree with Reviewer #1, and wish to stress that our precipitation adjustment is not constant in space and\ntime, as a 36.5% decrease of daily precipitation across the moist western windside of the SPI divide (> 15 m w.e. yr-\n1 in Fig. S3a) accounts for more mass reduction than in the drier leeside eastern regions (< 5 m w.e. yr-1 in Fig. S3a).\nIn other words, our RACMO adjustment has a dominating effect on peak (overestimated) precipitation values across\nthe wet windside of Patagonian divides. We therefore deem that our approach is best suited to locally eliminate the\nremaining mass change bias compared to GRACE/GRACE-FO. This is now clarified in the main text as follows: “In line\nwith previous studies, RACMO strongly overestimates orographic-forced precipitation across the steep Andes. The\nstatistically downscaled RACMO total precipitation is therefore adjusted and reduced by 36.5% over ice-covered grid-\ncells on a daily basis (Methods). This approach ensures that the precipitation adjustment has a dominating effect in\nregions where overestimates peak, notably across the steep and moist windside flank of the Andes, i.e., accumulation\n> 15 m water equivalent (w.e.) (Supplementary Fig. S1b). The justification for this adjustment is to eliminate the large\npositive bias relative to mass change records from the Gravity Recovery and Climate Experiment (GRACE) and Follow-\nOn (GRACE-FO) missions (Supplementary Fig. 4b,c), when combined with solid ice discharge. Such adjustment is not\nrequired for the statistically downscaled MAR product (Supplementary Fig. S4b).” In addition, adjusting RACMO\nprecipitation not only results in excellent agreement with observed PMET and modelled MAR time series in the\noverlapping period 1980-2020 (Figs S6b and c, and S7), but also better aligns with the vertical precipitation profile\nderived from PMET meteorological data (Figs S9g and h), demonstrating an improved precipitation distribution upon\nthe original downscaled RACMO product. These improved agreements support the robustness of our approach. This\nis clarified in the Methods as: “Adjusted precipitation also results in improved vertical profile agreement between\ndownscaled RACMO and PMET products (cyan line in Supplementary Fig. 9f), better aligns with glacier-integrated\nMAR and PMET precipitation estimates (Supplementary Fig. 7a), and further decreases the SMB RMSE to 4.9 m w.e.\ncompared to in-situ records (Supplementary Fig. 5e).”\n8a) Line 166-167. The ablation zones have expanded as you said, but the glaciers as a whole have shrinkage\nsignificantly, mainly due to the ablation zones/fronts calving away. Did you account for this area reductions? Which\narea are you using for this calculation? Only from 1986 to year 2016 the Patagonia glaciers lost an area of near 2200\nkm2 (Meier et al., 2018 https://doi.org/10.3389/feart.2018.00062), almost all within ablation areas.\nOur regional climate models do not account for ice dynamics (glacier retreat and thinning), instead we use a recent\nfixed glacier area and surface topography from the Randolph Glacier Inventory version 6 and the SRTM digital\nelevation model version 4, respectively. This approach is commonly used in the literature, and to our knowledge, no\nprevious study has accounted for glacier area and thickness changes to estimate SMB. To explore the response of\nintegrated SMB to glacier retreat, we estimated changes in glacier-integrated SMB and main components (runoff and\nprecipitation) when using the latest RGI 7.0 ice mask, i.e., smaller than the RGI 6.0 ice mask used in our estimates.\nWe discuss these points in the new Methods Section “Product strengths and limitations” as: “Patagonian glaciers and\nicefields have been retreating since the Little Ice Age, with a glacier extent reduction of roughly 2,200 km2 between\n1986 and 2016. However, our two regional climate models and statistical downscaling approach do not account for\nice dynamics processes driving glacier retreat and thinning. Instead, we use fixed surface elevation and glacier extent\nfrom the SRTM DEM (2000) and RGI 6.0 (1950-2010), yielding conservative present-day SMB estimates. However,\nglacier retreat can affect the spatial distribution of SMB components, notably in marginal ablation zones. For\ninstance, area-integrated meltwater runoff may decrease as ablation zone extent shrinks. In contrast, glacier thinning\nmay trigger pronounced melt-elevation feedback, i.e., enhanced meltwater runoff rates as temperature increases\nwith surface lowering. Quantifying these competing mechanisms requires coupling regional climate models to a\ndedicated ice dynamics module, as was previously done for the Greenland ice sheet. To date, such coupling has never\nbeen achieved for the small Patagonian glaciers and icefields. As an alternative, we explored the response of SMB\ncomponents to glacier retreat by quantifying the difference obtained when spatially integrating downscaled MAR\nSMB components over the reduced glacier extent of the latest RGI version 7 (RGI 7.0, 1950-2021). Relative to RGI 6.0\n(1950-2010), glacier area in RGI 7.0 declines by roughly 1,500 km2. We find that using the RGI 7.0 ice mask, glacier-\nintegrated total precipitation and runoff decrease by 3.7 Gt yr-1 (2.2%) and 13.0 Gt yr-1 (8.0%) on average for 1940-\n2023, hence resulting in a 6.4 Gt yr-1 SMB increase. These changes are insignificant as they remain lower than one\nSD of the long-term SMB components averages (1940-2023; Supplementary Fig. 11), i.e., SD = 15.2 Gt yr-1, 21.1 Gt\nyr-1, and 25.7 Gt yr-1 for total precipitation, runoff and SMB, respectively. Note that these changes do not account for\nmelt-elevation feedback resulting from associated glacier thinning, which could significantly increase surface runoff.”\nSupplementary Figure 11: SMB response to glacier retreat. Scatterplot comparing a downscaled MAR SMB spatially\nintegrated over RGI 6.0 and RGI 7.0 fractional ice masks for the period 1940-2023. b and c same as a but for total\nprecipitation and runoff, respectively. Relevant statistics including number of records (N), the slope (a) and intercept\n(b) of the regression line (y = ax + b), coordination coefficient (R2), mean model differences, i.e., as an absolute and\nrelative (%) value, and standard deviation (SD) of downscaled MAR spatially integrated over RGI 6.0.\n8b) What about the ELA elevation changes and/or the AAR? A better understanding of the ablation zone expansion\nas you concluded, needs a better discussion about total area reductions, ELA and AAR changes. I suggest to address\nthese concerns when discussing your results.\nFirst, we wish to stress that the ablation zone fraction presented here corresponds to 1 – AAR,\n• Ablation zone fraction = ablation zone area / total glacier area = 1 - AAR\nwith AAR being the Accumulation Area Ratio or:\n• Accumulation Area Ratio = accumulation zone area / total glacier area\nAs our manuscript examines the role of surface runoff in the ablation zone, we deem that discussing the ablation\nzone fraction is more appropriate.\nWe estimated changes in Equilibrium Line Altitude (ELA, SMB = 0) in both MAR and RACMO, showing a significant\nincreasing trend since 1940 (8.9 m per decade; p < 0.01), hence expanding the ablation zone higher up inland.\nELA in MAR and RACMO align well with a mean difference of 43 m for the overlapping period 1979-2023. Using MAR\nlong-term trend, the ELA moved from ~1220 m a.s.l. to ~1300 m a.s.l. between 1940 and 2023 (see Fig. below). For\nconciseness, the Figure below was not added to the manuscript. Instead, we now discuss these ELA changes in the\nmain text as: “Using long-term MAR timeseries, we find that the glacier ablation zones have significantly expanded\nsince 1940 (+0.9 ± 0.5% dec-1, p < 0.01), following near-surface atmospheric warming (r = 0.70, Fig. 4d). This is\nsupported by a significant long-term increase of the equilibrium line altitude (ELA, SMB = 0) by 8.9 ± 5.5 m dec-1 (p <\n0.01). As a result, the ELA migrated upward from ~1220 m a.s.l. to ~1300 m a.s.l. between 1940 and 2023, further\nexpanding the ablation zone area.”\nAs mentioned in comment 8a), we cannot estimate glacier area reduction as regional climate models, used here and\nin previous studies, do not account for ice dynamics (glacier retreat and thinning). This limitation is now discussed in\na dedicated section “Product strengths and limitations”.\n9a) Line 248-251. This asseveration is quite strong as it is. Assuming no significant fluctuation in frontal ablation is a\nBIG assumption. We know that calving glacier are dynamically responding with very anomalous responses. It is just\na matter of seeing Pio XI’s continuous advances.\nHere we give insight on mass loss to be expected, combining our current MAR and RACMO SMB estimates in the\noverlapping period 1979-2023 with published solid ice discharge estimates. We have reformulated our statement\nas: “At the current rate of mass loss, i.e., 24.2 ± 5.4 Gt yr−1 in MAR and 24.5 ± 5.4 Gt yr−1 in RACMO for 1979-2023\n(Fig. 1a), Patagonian glaciers could completely melt away within the next 220 years. These estimates should be\ninterpreted as conservative since none of MAR, RACMO, or the statistical downscaling procedure account for glacier\nretreat and thinning, nor for the evolution of calving processes.” Note that we now assume that 362.5 Gt of ice raises\nsea-level by 1 mm, as in Cogley et al. (2011).\n9b) I rather prefer to see a discussion about dynamic readjustments or delayed responses, before predicting the\ndisappearing of the huge volume of ice in the icefields.\nPlease, see our response to comments 8a) and 9a).\n10a) Figure 1. Maybe you can add a reference for the conversion between mass change and 1 mm of sea level?\nThank you, we added a reference to Cogley et al. (2011) and clarified as: “Mass change is converted to global sea-\nlevel rise equivalent assuming that 362.5 Gt of ice raises sea-level by 1 mm.”\n10b) I don’t see the RACMO coloured bands in b. After increasing the ZOOM to a very maximum, I only got a grey\ncolour underneath the SMB black line. Is this the band? Invisible for to me…Maybe leaving only MAR in here and\nshowing RACMO in the supplementary info?\nFor completeness and consistency with Fig. 1a and c, we prefer showing RACMO values as bands in Fig. 1b.\n10c) Add a line to legend indicating what are the stars (sites with ¬mass balance measurements). I know you said so\nin the text, but I think is better to see the figure with all the info in the legend. Thank you, this is done.\n11) Figure 3. I understand that only yellow polygon derived from MAR was used here but why not other variable with\nhigh r?\nFor completeness we now show the correlation between glacier-integrated total precipitation and SMB (Fig. 3b), to\nhighlight that SMB fluctuations are primarily driven by runoff changes (r = 0.80), and to a lesser extent by\nprecipitation changes (r = 0.61). Figures 3c and d thus focus on relevant variables that explain the significant runoff\nincrease (0.47 Gt yr-2, p < 0.01), that we identify as main driver of the significant SMB declining trend (-0.35 Gt yr-2,\np < 0.01), since precipitation remain overall steady (0.12 Gt yr-2, p > 0.05). Figures 3c and d demonstrate that the\nrunoff increase is linked to upper atmospheric temperature (T850) and subsequent near-surface temperature\nfluctuations. To clarify, T850 in Fig. 3d is extracted from the yellow contours, i.e., where runoff strongly correlates\nwith T850, shown in Supplementary Fig. 8a and c for MAR and RACMO, respectively. Furthermore, the yellow\ncontours conveniently outline a relevant region that encompasses the whole of Patagonia. This is ow clarified in the\nFigure caption.\n12a) Why 3c was obtained with N=84 and not 129 as in most of the other correlations in this and other figures (e.g.\nfigure 4)? It is because T2m is averaged over glacier areas? If so, why in 3b N=129 and not 84?\nThank you very much for pointing this out, this was a mistake from our side. Figure 3c shows 129 dots, i.e., 84 years\nof MAR (1940-2023) and 45 years of RACMO (1979-2023) data. This is now corrected in the revised manuscript. We\nalso clarified this point in the figure caption as: “Linear regression and relevant statistics, i.e., number of records (N),\nslope (a), intercept (b), and correlation (r), include both MAR (84 years) and RACMO (45 years) data sets.”\n12b) Please explain the “averaged over glacier” approach.\nNear-surface temperature is spatially averaged over glacier covered grid-cells, i.e., spatially summed and divided by\nthe number of glacier-covered grid-cells with ice fraction > 0. This is clarified as: “… are spatially averaged over glacier-\ncovered grid-cells (i.e., ice fraction > 0)” in the revised Figure caption.\n12c) In the high resolution version of Figure 3 (page 27 of the merged file) there are many purple and yellow points…I\npresume these are just printing errors, or have any meaning? I don’t see these “anomalous” points in Figure 3 in\npage 22. The same is happening in Figure 4 in page 28).\nThank you for reporting this. The high-resolution images uploaded to the Nature Communications server do not show\nthese artifacts. This must be due to compression or printing issues.\n13a) Figure 5. I don’t see the RACMO bands in d.\nWe acknowledge that the band is sometimes difficult to see, as the difference between MAR and RACMO remains\nsmall. For completeness, we prefer keeping the bands in the revised manuscript.\n13b) The yellow polygon is the MAR r>0.75 correlation area shown in Supple 8a, right? Why using this polygon and\nnot the RACMO derived >0.75?\nIn the original submission, we used the same outline for both MAR and RACMO obtained in Fig. 5a. The reason was\nthat spatial correlation maps between runoff and T850 were almost identical in MAR and RACMO. We agree that this\nwas inconsistent, and we now separately extract T850 for MAR and RACMO based on their respective yellow outline\n(r > 0.75) shown in revised Figure S8a and c. Figure 5 has been updated accordingly and the caption now mentions:\n“d Timeseries of T850 anomalies relative to 1960-1989 derived from ERA5 reanalysis (solid black line; 1940-2023)\nand retrieved from the yellow contour extracted for MAR and RACMO in Supplementary Fig. 8a and c, respectively.”\n13c) or the one correlating Z850 and pp? Seems to be arbitrary the selection of this specific one. Please add a phrase\nexplaining why this and not any other one was used.\nIn Fig. 5c, we deem appropriate to extract Z850 from a representative region that encompasses Patagonia, to show\nthat glaciers have experienced increased atmospheric pressure following a poleward shift in subtropical highs. In\ncontrast, we deem that using the region where precipitation and Z850 strongly correlate, i.e., centered on the ocean\nand not on land, is irrelevant. We clarified this in the caption of Fig. 5 as: “c Timeseries of Z850 anomalies relative to\n1960-1989 derived from ERA5 reanalysis (1940-2023) and extracted from a relevant region that encompasses the\nwhole of Patagonia (yellow outline in a-b).”\n15) Supplementary Figure 4. In the legend the time frame is different from the time scale of these figures. I suggest\nexplaining both (one is the time used for the whole series but here only using data from 2002 and 2023 for\ncomparison purposes).\nThank you for spotting this. We reformulated the Figure caption to stress that only 2002-2023 is shown: “Monthly\ncumulative mass change (MB = SMB - D) of Patagonian glaciers for 2002-2023 estimated as the difference between\nmodelled surface mass balance (SMB) from MAR (green line) and RACMO (blue line)”.\n16) Supplementary 6. Why N=45? I thought that you can compared the whole grid (N=129).\nUnfortunately, we did not run RACMO back to 1940, and thus compare MAR and RACMO for the overlapping period\n1979-2023 (i.e., N = 45 years). This is now clarified in the caption as: “Cross model correlation of a SMB, b snowfall,\nc rainfall, d runoff, e total melt, and f retention and refreezing, between MAR and RACMO statistically downscaled to\n500 m resolution for the overlapping period 1979-2023 (45 years).”\n17a) Supplementary Figure 7: I don’t see RACMO coloured bands either. Are there? The alignment between total pp\nof MAR and PMET is clear, and also between MAR snowfall and CR2Met, but I don’t see the rest of the relationships\nindicated in the legend. Maybe you can add the R2 between them?\nWe agree that the RACMO bands can sometimes be difficult to see as the differences between MAR and RACMO\nremain overall small, but we prefer showing them for completeness. As suggested, we added R2 estimates in\nSupplementary Fig. 7a and modified the caption as: “Coefficient of determination (R2) between the different data\nsets are also listed.”\n17b) Using max fluxes in primary “Y” is a bit confusing as you are only discussing precipitations in this figure and I\nthink m we/yr is enough.\nThank you, we reformulated the primary Y-axis as: “Precipitation fluxes (Gt yr-1)”. As our manuscript extensively\ndiscusses mass fluxes in Gt yr-1, we decided to keep the primary axis unit as is.\n17c) Maybe saying something about the conversion from m we/yr to GT/yr? or moving GT/yr to the secondary Y?\nWe clarified how glacier-integrated SMB (and components) are estimated in the revised Section “Spatial integration”\nas: “To estimate mass fluxes in Gt, we spatially integrate SMB and its individual components over ice-covered grid-\ncells. To that end, we use RGI 6.059 glacier outlines to derive a fractional ice mask with grid-cell values ranging from\n0 (ice free) to 1 (fully ice-covered). The use of a fractional mask (0 to 1 in Supplementary Fig. 10a) instead of a binary\none (0 or 1 in Supplementary Fig. 10b) ensures that large glacier area overestimate is prevented, notably for small\nmountain glaciers covering less than 0.25 km2 (i.e., area of a 500 m grid-cell). Here we estimate glacier-integrated\nSMB and individual components (M in Gt) as,\nM = ∑ m * A * I * 1e-6\ni,j i,j i,j\nwhere m is the SMB (or individual components) expressed in kg m-2 (or mm w.e.) and estimated for grid-cell i,j; A\ni,j i,j\nis the area of grid-cell i,j set to 0.25 km2; I is the fractional ice mask of grid-cell i,j.”\ni,j\n18) Supplementary 8. The yellow line for RACMO?\nThe yellow outline is now shown for both MAR and RACMO. These areas are used to extract T850 time series for\nMAR and RACMO in Figs. 3c and 5d, respectively. See also our response to comment 13c).\nFormat details:\n19) RACMO and RACMO2 are used alternatively, maybe use just one acronyms to avoid misunderstandings.\nThank you, for consistency we now use “RACMO” across the main manuscript and supplements.\n20) In Line 444 you said “cyan in Supplementary Fig.9e”… but I think you are referring to Figure 9h.\nThank you for spotting this. The figure and caption have been modified accordingly.\nAdditional references:\n• The GlaMBIE Team. Community estimate of global glacier mass changes from 2000 to 2023. Nature (2025).\nhttps://doi.org/10.1038/s41586-024-08545-z\n• Cogley, J. et al. Glossary of glacier mass balance and related terms (UNESCO-IHP, IHP-VII Technical documents\nin hydrology 86, IACS contribution 2, 2011). http://unesdoc.unesco.org/images/0019/001925/192525E.pdf\nReviewer #2 (Remarks to the Author):\nThis is a comprehensive and well-executed study on the surface mass balance (SMB) and long-term mass loss of\nPatagonian glaciers, which are key contributors to global sea-level rise. The integration of high-resolution regional\nclimate models and statistical downscaling provides robust framework for analyzing SMB over an extended period\n(1940–2023). The authors effectively combine these results with GRACE/GRACE-FO satellite data and solid ice\ndischarge estimates, offering valuable insights into the drivers of glacier mass loss. The study is particularly notable\nfor identifying the poleward shift of subtropical highs and its implications for warm air advection, demonstrating that\nlong-term fluctuations in surface runoff are a more significant driver of SMB trends than precipitation.\nThe paper is very clearly written and well-presented, especially considering the complexity of the analysis, and I see\nno errors or omissions in the detail of the paper. I only have minor suggestions exclusively pertaining to figures:\nThank you very much for your enthusiasm and strong support of our study. Please, find below our point-by-point\nresponse to your comments.\n1) Size may be an issue, but it would be helpful if Fig 1a had a scalebar. Fig 2 is also missing a scalebar. What is the\nbackground hillshade image the data is overlaid on?\nWe added a scalebar to each map in the manuscript and supplements. The hillshade is based on the SRTM DEM\ndown-sampled to 250 m spatial resolution. This is now clarified in captions where appropriate as: “The background\nhillshade is derived from the SRTM DEM down-sampled to 250 m spatial resolution.”\n2) Several of the plots are not colorblind compatible, particularly with regards to red/green colorblindness in Fig 1a,\n1c, 2a and 2b.\nTo remain consistent across the manuscript e.g., MAR data displayed in green, RACMO in blue, observations in red\nand orange, etc…, we prefer keeping Figs. 1a-c as are. The same holds for Figs. 2a,b which use a color palette\nconsistent with our previous publications in Nature Communications e.g., Fig. 3 in Noël et al. (2022) or Fig. 4 in Noël\net al. (2023). If necessary, we can modify these colors, but this may strongly impact the display of all figures and may\nnot necessarily improve readability. We leave this decision to the discretion of the Editor.\nAdditional references:\n• Noël, B., Lenaerts, J.T.M., Lipscomb, W.H. et al. Peak refreezing in the Greenland firn layer under future\nwarming scenarios. Nat Commun 13, 6870 (2022). https://doi.org/10.1038/s41467-022-34524-x\n• Noël, B., van Wessem, J.M., Wouters, B. et al. Higher Antarctic ice sheet accumulation and surface melt rates\nrevealed at 2 km resolution. Nat Commun 14, 7949 (2023). https://doi.org/10.1038/s41467-023-43584-6"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have reviewed your new version of the manuscript and the answers to my previous report, as well as the responses to the second reviewer. I think you have considered and appreciated my comments as a positive support to your work. My concerns were in detail addressed and I'm satisfied with all of them, however, I still have a couple of points to be mentioned:\nRegarding the use of a constant value (and negative) for the period pre-1968 in Argentina or pre-1975 in Chile is very arguable. I understand that you don't have many more options but considering that the MB MAR shows a positive trend in the oldest period, I suggest discussing a bit further the use of these constants. Maybe it is enough to say that this assumption could be improved in such or such a way. For example, using TRIMETROGON oblique images for reconstruct surface topographies or using historical glacier outlines for reconstructing altitudes from old works. Structure-from-motion is widely used at present (see https://doi.org/10.1016/j.geomorph.2017.01.008) or the use of outlines is an old method still valid (https://doi.org/10.3189/S0260305500001178).\nELA: I suggest adding a couple of references in the discussion section about the ELA trends you are obtaining for discussing differences or coincidences (e.g. https://doi.org/10.1038/s41598-021-95725-w or older papers estimating future changes like https://doi.org/10.1016/j.gloplacha.2006.11.025).\nAlso I suggest changing the following text:\n“……This is supported by a significant long-term increase of the altitude of the\nequilibrium line altitude (ELA, SMB = 0) by 8.9 +- 5.5 m dec−1 (p < 0.01). As a result,\nthe ELA migrated upward from ∼1220 m a.s.l. to ∼1300 m a.s.l. between 1940 and\n2023, further expanding the ablation zone area.”"}, "author_response": "Response letter\nWe would like to thank the Editor for inviting us to revise our manuscript, and Reviewer #1 for his positive comments.\nBelow you will find in blue our point-by-point responses to the Reviewer’s comments, while changes applied to the\nrevised manuscript are displayed in red.\nReviewer #1 (Remarks to the Author):\nI have reviewed your new version of the manuscript and the answers to my previous report, as well as the responses\nto the second reviewer. I think you have considered and appreciated my comments as a positive support to your\nwork. My concerns were in detail addressed and I'm satisfied with all of them, however, I still have a couple of points\nto be mentioned: Thank you, you can find our responses below.\nRegarding the use of a constant value (and negative) for the period pre-1968 in Argentina or pre-1975 in Chile is very\narguable. I understand that you don't have many more options but considering that the MB MAR shows a positive trend in\nthe oldest period, I suggest discussing a bit further the use of these constants. Maybe it is enough to say that this\nassumption could be improved in such or such a way. For example, using TRIMETROGON oblique images for\nreconstruct surface topographies or using historical glacier outlines for reconstructing altitudes from old works.\nStructure-from-motion is widely used at present (see https://doi.org/10.1016/j.geomorph.2017.01.008) or the use\nof outlines is an old method still valid (https://doi.org/10.3189/S0260305500001178).\nAs previously mentioned, we do not quantify solid ice discharge but rather use published data to conduct a\nmeaningful comparison with GRACE/GRACE-FO mass change. We acknowledge that published estimates, i.e., annual\naverages for two periods, may not fully capture solid ice discharge seasonal and interannual variabilities within their\nuncertainty range. This is now clarified as: “Note that these published solid ice discharge estimates, i.e., annual\naverages for two periods, may not fully capture interannual and seasonal variabilities within their uncertainty range.”\nWe deem that elaborating on how solid ice discharge estimates could be improved is technical, i.e., beyond the scope\nof this study, and may divert readers from our core findings. However, we agree that solid ice discharge estimates\ncould be improved in future work by extending available times series to cover the pre-satellite era, to better capture\n(long-term) mass change variability. This is clarified as: “In future work, these solid ice discharge estimates could be\nimproved by exploiting historical imagery to cover the pre-satellite era77,78, hence capturing potentially large long-\nterm variability.\nELA: I suggest adding a couple of references in the discussion section about the ELA trends you are obtaining for\ndiscussing differences or coincidences (e.g. https://doi.org/10.1038/s41598-021-95725-w or older papers estimating\nfuture changes like https://doi.org/10.1016/j.gloplacha.2006.11.025). Also I suggest changing the following text:\n“……This is supported by a significant long-term increase of the altitude of the equilibrium line altitude (ELA, SMB =\n0) by 8.9 +- 5.5 m dec−1 (p < 0.01). As a result, the ELA migrated upward from ∼1220 m a.s.l. to ∼1300 m a.s.l.\nbetween 1940 and 2023, further expanding the ablation zone area.”\nSupplementary Table 3 now compares MAR and RACMO equilibrium line altitude (ELA) estimates for NPI and SPI\nwith remote sensing and model estimates presented in Table 1 of Bravo et al. (2024). We find that MAR ELAs agree\nwithin uncertainties with previous work; RACMO overestimates ELAs by 96 and 151 m compared to MAR in NPI and\nSPI, respectively (i.e., see model uncertainty in Supplementary Table 3). This supports our finding that RACMO has a\nlarger ablation zone area than MAR, as shown e.g. in Fig. 4a covering the three icefields and neighboring mountain\nglaciers. This is now discussed as: “This stems from a significant long-term increase of the equilibrium line altitude\n(ELA, SMB = 0) by 8.9 ± 5.5 m dec−1 (p < 0.01). As a result, the ELA migrated upward from∼1220 m a.s.l. to∼1300 m\na.s.l. between 1940 and 2023, further expanding the ablation zone area. This is further supported by the fact that\nMAR captures, within uncertainties, the ELA of NPI and SPI derived from previous observational and model estimates\n(Supplementary Table 3).”"}]

Decision letter