We thank both reviewers for their thorough and constructive comments. The manuscript has been substantially revised in response. Major changes include:

(1) the ArcPy-based DEM reprojection, alignment, and differencing workflow has been replaced with a fully open-source Python pipeline using rasterio and rioxarray — DEM co-registration via demcoreg was already open-source in the submitted version and is unchanged;

(2) a 1975 glacier inventory (DGA, 2011; the closest available to the 1954 IGM DEM baseline) has been incorporated as the reference outline for all 1954-based epochs, extending the area change time series and improving the accuracy of historic mass balance calculations;

(3) a Pléiades 2025 DEM, built by DGA from Pléiades-1A stereoscopic imagery acquired on 17 April 2025 and processed using the Ames Stereo Pipeline (ASP, v3.5.0; Beyer et al., 2018), was obtained under the Chilean Transparency Act and added to the analysis, providing full-complex mass balance estimates for 2000–2025 and complementing the UAV surveys which cover 11 glaciers in the Cerro Blanco and Las Termas subcomplexes;

(4) a bias-corrected 30 m ASTER DEM for 2018 was obtained from R. Hugonnet and incorporated into the mass balance time series as an intermediate epoch; a 2008 scene was also evaluated but excluded after quality filtering due to near-zero median elevation differences over glacier terrain indicative of snow contamination;

(5) uncertainty propagation has been revised following Hugonnet et al. (2022) and Rolstad et al. (2009) using the xDEM Python package, replacing the previous NMAD-only approach; and

(6) the E-Divisive changepoint analysis has been removed and replaced with Sen’s slope and Mann–Kendall tests for monotonic trend detection, following Reviewer 2’s concern about the appropriateness of distributional changepoint methods for gradual climate trends.

All responses below are numbered to correspond to the reviewers’ original comments. Manuscript line numbers refer to the revised version.

RC1.1

I think some care needs to be taken with the vertical datum used by each of your two DEMs, especially the 1954 topographic map DEM. The SRTM that you are using is provided as the height above the EGM96 geoid, while the 1954 map is most likely the same as what is reported by Farías-Barahona et al. (2020), orthometric height (m a.s.l.) measured from tide gauges along the coast. I think there’s a good chance of large differences between the vertical datum of the topographic maps and the EGM96 geoid, which would likely show up as differences (errors) on the off-glacier areas of your DEMs, even after co-registration. It might be comparatively minor, but I think it’s at least worth mentioning in the data description. That said, I think you might also want to mask your “source” DEM for the co-registration - it looks as though there are areas outside of the DEM that are still being included for the co-registration, which may also have an impact here.

Response:

We thank the reviewer for raising this point. Source DEMs were already being masked during co-registration; we have clarified this in §3.4. We have also conducted a thorough review of the vertical datums across all DEMs. The SRTM DEM is provided as heights above the EGM96 geoid (orthometric). The 1954 IGM DEM is based on Chilean topographic maps using the PSAD56 datum, which approximates orthometric heights at this latitude and is consistent with EGM96 to within approximately 1 m at 37°S — a difference we consider negligible relative to the overall uncertainty of the IGM DEM (stable-terrain NMAD of ~26 m). The UAV-derived DEMs and Pléiades 2025 DEM were generated in WGS84 ellipsoidal heights; we now apply a pixel-wise EGM96 geoid undulation correction (N ≈ 22–24 m at 37°S, 71°W) using the pyproj library prior to co-registration and differencing, converting all ellipsoidal DEMs to orthometric heights consistent with SRTM. The ASTER DEM provided by R. Hugonnet was co-registered to TanDEM-X (ellipsoidal) and has also been geoid-corrected using the same procedure. These corrections and their spatial variation across the study area are now described in §3.3.1, §3.3.3, and §3.3.4.

RC1.2

I don’t think the approach taken in §3.5, to convert the off-glacier NMAD value to m w.e. and use this as the uncertainty for the mass balance is correct. Because the mass balance is computed as an aggregation of the elevation change, the uncertainty of the volume change needs to be calculated by propagating the uncertainty of a single elevation change pixel through the spatial aggregation — see, for example, Hugonnet et al. 2022 or Rolstad et al. 2009. Because you are using python for some aspects of this, you might want to look at the xDEM python package (https://xdem.readthedocs.io), which includes tools for estimating the DEM uncertainty and propagating that uncertainty to volume changes.

Response:

We agree with this comment and thank the reviewer for the specific methodological guidance. The previous approach was incorrect and has been fully revised. We now implement a two-step uncertainty propagation pipeline following Hugonnet et al. (2022) and Rolstad et al. (2009) using the xDEM Python package.

In the first step, we infer a spatially variable per-pixel elevation error map σ(x,y) from stable off-glacier terrain using the xDEM function infer_heteroscedasticity_from_stable. This function models error as a function of terrain slope and maximum curvature via a two-dimensional binning approach, capturing the well-documented dependence of DEM error on terrain ruggedness (Hugonnet et al., 2022). Where insufficient stable terrain pixels were available for heteroscedastic modelling (e.g. the Cerro Blanco UAV DEM, which has only ~1,200 stable pixels), a spatially uniform error equal to the stable-terrain NMAD was used as a fallback.

In the second step, the spatial correlation structure of elevation errors was characterized by fitting a spherical variogram model to the empirical variogram of off-glacier elevation differences, computed using the Dowd estimator. The uncertainty in mean glacier elevation change σ_d̅h was then propagated following Rolstad et al. (2009, Eq. 14), which accounts for spatial autocorrelation: σ_d̅h = σ_pixel × √(A_corr / A_glacier), where A_corr = πL²/2 is the area of spatial correlation and L is the variogram range. Mass balance uncertainty is then derived as σ_MB = σ_d̅h × ρ_ice/ρ_water / Δt. This approach substantially reduces uncertainty estimates compared to the previous NMAD-only method. For example, for Glaciar Nevado (Pléiades 2000–2025), uncertainty decreased from ±0.23 to ±0.10 m w.e. yr⁻¹, yielding a signal-to-noise ratio of 7.2. The revised methodology is described in full in §3.5 and the code is available in Notebook 03 of the updated GitHub repository.

RC1.3

Were you able to digitize the glacier outlines from the topographic map? I know that not every topographic map will include glaciers as a distinct feature, but given that other studies such as Farías-Barahona et al. (2020) were able to include the glacier outlines digitized from similar maps, I figure it’s at least worth asking about. If glaciers were not identified on the map used, this should be mentioned in the appropriate section as an explanation for why it is not included. If it is possible to include the 1954 outlines, I think it’s worth re-doing at least part of the analysis by comparing the mass changes computed with these outlines.

Response:

Glacier outlines are not clearly delineated as a distinct feature on the 1954 IGM topographic map, precluding direct digitization of 1954 glacier extents. We have added a sentence to §3.3.1 explaining this. However, we have obtained a 1975 glacier inventory for the Nevados de Chillán complex digitized from Landsat imagery, which provides the earliest available outlines and substantially improves the analysis. The 1975 inventory is now used as the reference outline for all 1954-based epochs, replacing the 2000 outlines used in the submitted version. The 1975 inventory records a total complex area of 16.02 km², compared to 2.91 km² in 2000, representing a loss of 81.8%. Mass balance for 1954-based epochs is now calculated using a time-weighted mean area approach following Cogley et al. (2011), with the 1975 inventory serving as the start-of-period area estimate (the closest available outline to the midpoint of the 1954–2000 observation period). The use of the 1975 rather than 1954 outlines, and its implications for mass balance estimates, is discussed in §3.4 and §3.5.

We note that 1954-based mass balance estimates for the Las Termas subcomplex yield slightly positive but statistically insignificant values (+0.07 ± 0.47 m w.e. yr⁻¹ for SRTM−IGM; +0.06 ± 0.30 m w.e. yr⁻¹ for Pléiades−IGM). We attribute this to the geometric mismatch between the large 1975 Las Termas polygons (4.99 km²) and the much smaller current glacier footprint (0.65 km²), which results in inclusion of recently deglaciated terrain with near-zero elevation change in the mean area calculation. These estimates are flagged as unreliable in §4.3 and excluded from trend analysis.

Minor and Specific Comments

1. 101–102

I would not characterize these as “images” - the IGM DEM you describe is based on topographic maps which were generated from many images, but what you describe is not using those images. Similarly, the SRTM DEM is not an “image”, but rather a DEM.

Response:

Corrected throughout. The IGM 1954 dataset is now described as a topographic map-derived DEM and SRTM is described as a DEM, not an image. Table 1 has been updated accordingly.

1. 145

Why would manual flying impact the accuracy of the DEM?

Response:

Manual flight paths result in suboptimal image overlap, departing from the recommended ~70–80% forward and side overlap, which reduces the number of common tie points available for Structure-from-Motion photogrammetric reconstruction and decreases DEM accuracy in areas of reduced overlap. This is reflected in the Metashape processing statistics: 15 images failed to align for Cerro Blanco (versus zero for Las Termas), and the RMS reprojection error was higher for Cerro Blanco (0.26 px versus 0.22 px), consistent with reduced tie-point density in manually flown sections. The higher stable-terrain NMAD of the Cerro Blanco UAV DEM (~12 m) relative to the Las Termas DEM (~5 m) after co-registration is consistent with this explanation. We have added this clarification to §3.1 and §3.2.

1. 162

Should this be two files, if both points and lines were digitized?

Response:

The reviewer is correct. The digitization produced two shapefiles: one containing spot height points and one containing contour polylines. The contour shapefile was converted to a point file by extracting vertices. The resulting vertices shapefile was then merged with the spot heights shapefile into a single combined point file for interpolation. The manuscript has been corrected to reflect this two-stage workflow in §3.3.1.

1. 163

What is the density/spacing of the vertices for the polylines? Is this consistent throughout the area digitized?

Response:

The mean distance between adjacent vertices is 36.4 ± 21.5 m (standard deviation), computed as the mean nearest-neighbor distance for contour vertices converted to points. Spacing is therefore variable rather than perfectly uniform, reflecting the curvature of contour lines — denser in areas of complex topography and sparser on gentle slopes. This information has been added to §3.3.1.

1. 165

Why IDW, rather than a method like Topo to Raster?

Response:

Topo to Raster implements the ANUDEM algorithm, which uses the river network to enforce hydrological consistency in the interpolated surface. While ANUDEM is generally superior to IDW for hydrological modelling applications, its primary advantage (the river network correction) is not appropriate for our use case. The streams in the glaciated zone of Nevados de Chillán vary substantially in shape and position as the glaciers retreat, meaning that the river network assumed by ANUDEM would not be representative of the 1954 landscape being reconstructed. IDW is therefore appropriate here and is consistent with the approach used by Farías-Barahona et al. (2020) for similar topographic map-derived DEMs in the Chilean Andes. A brief justification has been added to §3.3.1.

1. 171–179

I don’t think this is a correct characterization of the approach taken by both Dussaillant et al. (2019) or Hugonnet et al. (2021). Both of those approaches used many ASTER DEMs, but not to create a single composite DEM. There are issues with ASTER DEMs over steep (north-facing) slopes or low-contrast areas such as snow and ice, but unknown/unquantified signal penetration in the SRTM DEM also poses an issue. I am also not sure that characterizing the SRTM C-band DEM as “comparing images taken seconds apart” is correct, given that it is a mosaic of acquisitions over an 11-day period (Rabus et al., 2003).

Response:

We apologize for the inaccurate characterization and have revised §3.3.2 and §3.3.3 accordingly. The text previously confused individual ASTER stereo DEMs with the ASTER GDEM composite product; we have corrected this distinction. Both Dussaillant et al. (2019) and Hugonnet et al. (2021) used individual ASTER stereo acquisitions rather than a composite, and we now describe this correctly. The description of SRTM as “comparing images taken seconds apart” has been removed and replaced with an accurate description of the 11-day acquisition period (Rabus et al., 2003). We have also revised the characterisation of ASTER limitations: ASTER is known to suffer from signal attenuation and reduced vertical accuracy over steep north-facing slopes due to shadowing effects, and over low-contrast snow and ice surfaces (Berthier et al., 2018). We have additionally acknowledged the potential for unquantified C-band signal penetration in SRTM. SRTM was acquired in February 2000 at the end of the austral summer ablation season, when glacier surfaces in the ablation zone consist predominantly of wet ice into which C-band penetration is negligible (<1 m; Rignot et al., 2001). Any residual penetration into firn in the accumulation zone would cause a systematic underestimation of mass loss, and our estimates should therefore be considered conservative in this respect. This is consistent with the approach of Farías-Barahona et al. (2020) for nearby Chilean glaciers.

1. 186

Neither Shean et al. 2023 nor Berthier et al. 2024 are included in your references.

Response:

Both references have been added to the reference list.

1. 200, 201, 258

Use MB rather than Mass–Balance; include a citation for the chosen snow/ice density value; Glaciar Nevado.

Response:

Corrected throughout: the non-standard hyphenated form “Mass-Balance” has been replaced with “MB”; Huss (2013) is now cited for the ice density value of 850 kg m⁻³; “Glaciar Nevado” was corrected.

Fig. 3

Showing this as a stretched raster, rather than classified, would make it easier to see the local variation in elevation difference. I would also show the elevation difference of off-glacier areas, though this could be done with some transparency.

Response:

We have updated Fig. 3 to use a stretched (continuous) color ramp rather than a classified scheme, which better conveys local variability in elevation change within individual glaciers. Regarding showing off-glacier elevation differences: we are cautious about displaying off-glacier elevation change in this volcanic setting, as apparent elevation changes in ice-free terrain may reflect volcanic uplift, subsidence, or other non-glaciological processes rather than measurement error. To avoid ambiguity, we have retained the glacier-masked display but reduced the opacity of the off-glacier basemap to provide topographic context. We note this limitation in the figure caption.

1. 290

What is the spatial scale over which the differences are autocorrelated?

Response:

Variogram ranges (correlation lengths) are now reported for each DEM pair in the revised §3.5. Ranges vary substantially by DEM type: 882 m for the IGM 1954 (reflecting large-scale systematic errors in digitized topographic maps), 13 m for the ASTER 2018 DEM (indicating largely uncorrelated pixel-level errors), 505–637 m for the UAV DEMs, and 187 m for Pléiades 2025. These values are physically plausible and consistent with previously reported correlation lengths for similar DEM types (Hugonnet et al., 2022).

Fig. 4

Change the labels in the legend to be a line, rather than the outline of a rectangle, using matplotlib.lines.Line2D.

Response:

Updated. All legend handles in Fig. 4 and other figures using glacier outline overlays (Figs. 1 and 3) now use matplotlib.lines.Line2D, consistent with the reviewer’s suggestion.

1. 300–304

How does the % area change compare to the area-averaged elevation difference for glaciers in different sub-complexes? Did you also compare the area change to the area-averaged elevation/mass change for each glacier?

Response:

We have added a per-glacier analysis comparing area change (2000–2019) to mean elevation change (2000–2025, Pléiades) for each individual glacier in §4.2. Spearman correlation between area change (%) and mean elevation change rate was ρ = −0.018 (p = 0.94, n = 20 glaciers), indicating no significant relationship between area retreat and surface lowering rate at the individual glacier scale. This decoupling is discussed in §4.2 and likely reflects the influence of volcanic geothermal heating on basal melt, which drives ice loss independently of surface area.

1. 360 and Figs 6–9

Is this meant to be January and February? In Fig. 6-9, it seems as though the trends are only shown if they are significant, and there is no trend shown in Fig. 7a-c. Additionally, include the p-value. Show vertical changepoints on all plots and include the year in the legend. Also show the slope pre- and post-change.

Response:

Following the comments of Reviewer 2 (see RC2 response below), we have removed the E-Divisive changepoint analysis from the hydroclimate section entirely and replaced it with Sen’s slope trend lines and Mann–Kendall significance tests. Segmented regression was applied to test for acceleration in trends; it did not converge on a statistically robust breakpoint for any series, confirming that hydroclimatic change is better characterized as monotonic rather than abrupt. Pre- and post-2000 slopes are now shown on all panels where trends are significant, and p-values are reported in the legend for all panels. The analysis has been extended to include March (DJFM rather than DJF) as the most snow-free and glacier melt-relevant month for the study area. Non-significant panels are now clearly labelled as such.

1. 400

Why would the single year show a larger magnitude than multiple years?

Response:

We have added an explanation to §5.1. The substantially larger magnitude of the DGA’s single-year estimates relative to our multi-decadal averages reflects the well-documented sensitivity of annual mass balance to inter-annual variability in snowfall and temperature: a particularly warm or dry year may produce several times the long-term average annual loss, while an anomalously cold or snowy year may approach near-zero balance. Multi-year geodetic averages smooth this variability and better represent the sustained glacier response to climate forcing. This limitation of point-in-time estimates is now explicitly acknowledged.

1. 407–420

I think you should also compare to Hugonnet et al. (2021) per-glacier estimates and WGMS annual mass change estimates for the same glaciers.

Response:

We have added a comparison to the regional mass balance estimates from Hugonnet et al. (2021) in §5.1. We note that Hugonnet et al. (2021) aggregates estimates at the regional and subregional level and does not provide individual glacier-scale values for small glaciers of this size; per-glacier comparisons are therefore not possible from that dataset. Our whole-complex 2000–2025 estimate (−0.51 ± 0.17 m w.e. yr⁻¹) is discussed in relation to the Hugonnet et al. (2021) regional estimate of −0.67 ± 0.15 m w.e. yr⁻¹ for 2000–2019 in §5.1. Furthermore, as noted in our response to RC2.3, we obtained the bias-corrected 2018 ASTER DEM from R. Hugonnet directly, providing a direct methodological connection to the Hugonnet et al. (2021) dataset and strong independent validation of our 2000–2018 estimates.

1. 418

You could also discuss the value of longer-term local studies like this one — these are comparatively rare worldwide, but especially in the Andes.

Response:

We have added a paragraph to §5.1 emphasizing the scientific value of long-term, high-resolution local studies in the context of Andean glaciology, where such records remain rare despite the importance of glacier meltwater for downstream communities.

1. 433–434

Why the disagreement [with previous studies on summertime precipitation]?

Response:

The discrepancy between our finding of no significant trend in summertime (DJFM) precipitation and the regional decreases reported by Garreaud et al. (2020) and Rubio-Álvarez & McPhee (2010) likely reflects differences in spatial scale and record length. Both referenced studies analyzed regional reanalysis or multi-station datasets spanning broader domains and longer periods, whereas our analysis is based on two local stations with records from 1980. At the station scale, DJFM precipitation exhibits substantially higher interannual variability than annual totals, reducing statistical power to detect trends over ~40-year records. We have added a sentence to §5.1 acknowledging that our local finding does not contradict regional-scale drying, which is primarily driven by winter and spring precipitation reductions rather than summer precipitation (Garreaud et al., 2020).

1. 448

Can you report the surface area change here?

Response:

We have added glacier surface area values to §4.2: the complex declined from 16.02 km² in 1975 to 2.91 km² in 2000 and 1.90 km² in 2019, representing an 88% reduction over the study period.

1. 518

I have to think this would be somewhat negligible compared to the uncertainties in the SRTM, which are typically reported as being at least 5 m or so.

Response:

We agree with the reviewer’s assessment. The potential contribution of volcanic uplift to apparent ice surface lowering is negligible relative to the dominant SRTM uncertainty (±5–16 m; Rodriguez et al., 2006), and geodetic uplift rates at active volcanic systems rarely exceed a few centimeters per year during non-eruptive periods. We have simplified the text at §5.2 to acknowledge this source of uncertainty briefly without overstating its magnitude.

We thank both reviewers for their thorough and constructive comments. The manuscript has been substantially revised in response. Major changes include:

(1) the ArcPy-based DEM reprojection, alignment, and differencing workflow has been replaced with a fully open-source Python pipeline using rasterio and rioxarray — DEM co-registration via demcoreg was already open-source in the submitted version and is unchanged;

(2) a 1975 glacier inventory (DGA, 2011; the closest available to the 1954 IGM DEM baseline) has been incorporated as the reference outline for all 1954-based epochs, extending the area change time series and improving the accuracy of historic mass balance calculations;

(3) a Pléiades 2025 DEM, built by DGA from Pléiades-1A stereoscopic imagery acquired on 17 April 2025 and processed using the Ames Stereo Pipeline (ASP, v3.5.0; Beyer et al., 2018), was obtained under the Chilean Transparency Act and added to the analysis, providing full-complex mass balance estimates for 2000–2025 and complementing the UAV surveys which cover 11 glaciers in the Cerro Blanco and Las Termas subcomplexes;

(4) a bias-corrected 30 m ASTER DEM for 2018 was obtained from R. Hugonnet and incorporated into the mass balance time series as an intermediate epoch; a 2008 scene was also evaluated but excluded after quality filtering due to near-zero median elevation differences over glacier terrain indicative of snow contamination;

(5) uncertainty propagation has been revised following Hugonnet et al. (2022) and Rolstad et al. (2009) using the xDEM Python package, replacing the previous NMAD-only approach; and

(6) the E-Divisive changepoint analysis has been removed and replaced with Sen’s slope and Mann–Kendall tests for monotonic trend detection, following Reviewer 2’s concern about the appropriateness of distributional changepoint methods for gradual climate trends.

All responses below are numbered to correspond to the reviewers’ original comments. Manuscript line numbers refer to the revised version.

RC2.1

The large number of datasets used in this study means that presenting a coherent story about the analysis is not an easy task. I think that the authors could benefit from understanding what their analysis can and cannot say. Most notably, care should be given to how the polygons (from the year 2000 or 2019) are used in analysis of 1954 datasets – I note, for instance, that the 2000 polygons were used to derive stable terrain, but this would inevitably include ice melted from 1954-2000 (and thus, the terrain would not be ‘stable’). As glacier area change also cannot be included in the analysis, is it more accurate to report the overall figures as ice elevation change rather than mass balance?

Response:

We thank the reviewer for this important methodological concern. To address it, we have incorporated the 1975 glacier inventory as an additional mask during co-registration, constructing a conservative union of the 1975, 2000, and 2019 glacier outlines (§3.4). This substantially reduces the inclusion of formerly glaciated terrain as “stable” ground during co-registration of the 1954 IGM DEM. We acknowledge that terrain glaciated in 1954 but deglaciated by 1975 remains unmasked; however, given the magnitude of area loss between 1975 and 2000 (81.8%), we consider the 1975 inventory to be a reasonable conservative proxy. The use of the 1975 inventory also allows us to report results as specific mass balance using a time-weighted mean area approach following Cogley et al. (2011), described in §3.5, rather than elevation change alone, which we consider more physically meaningful. We are open to the reviewer’s feedback on whether this addition is considered sufficient to justify reporting mass balance rather than elevation change.

RC2.2

I’m also not sure that the change point analysis is the appropriate statistical test for a time series where you are not expecting to see a system change, but rather a linear trend. Care with regard to this analysis would improve confidence in the study.

Response:

We agree with the reviewer and have removed the E-Divisive changepoint analysis entirely. As the reviewer notes, distributional changepoint methods assume sudden regime shifts, which is not appropriate for gradual anthropogenically-driven climate trends. Hydroclimatic trends are now characterized solely using Sen’s slope estimator and Mann–Kendall significance tests, which are appropriate for monotonic trends in non-normally distributed time series. To confirm the absence of a statistically robust breakpoint, we additionally applied segmented (piecewise linear) regression; this did not converge on a significant breakpoint for any series, confirming that the observed changes are better described as continuous monotonic trends than abrupt regime shifts. Pre- and post-2000 Sen’s slope estimates are reported to enable comparison with the glacier mass balance record, which has its baseline at 2000.

RC2.3

There are also opportunities that could improve the study. The analysis to 2024 includes 11 glaciers due to practical considerations with the UAV survey. However, you could make use of secondary datasets – e.g. Hugonnet et al. (2021) – to report elevation thinning to 2020 for all glaciers. I pose this suggestion only as an idea for the authors to consider; the study is complete without this extra analysis.

Response:

We thank the reviewer for this suggestion. We contacted R. Hugonnet directly and he generously shared his bias-corrected 30 m ASTER DEMs for our study tile. After quality filtering (stereo correlation threshold ≥60 following Hugonnet et al., 2021, spatial outlier removal vs SRTM, snow bias assessment, and a minimum 30% glacier coverage threshold), one scene was retained: 1 March 2018. A March 2008 scene was also evaluated but excluded due to near-zero median elevation differences over glacier terrain indicative of residual snow bias. The 2018 ASTER result (−0.51 ± 0.01 m w.e. yr⁻¹ for the whole complex) matches our Pléiades 2025 result almost exactly, providing strong independent validation of both datasets. We additionally include a Pléiades 2025 DEM obtained from DGA under the Chilean Transparency Act, built from Pléiades-1A stereoscopic imagery acquired on 17 April 2025 and processed using the Ames Stereo Pipeline (ASP, v3.5.0; Beyer et al., 2018), which covers the full complex at 0.73 m ground sampling distance and extends the full-complex record to 2025. Together these additions substantially increase temporal resolution between the 2000 SRTM baseline and the 2024 UAV surveys.

Specific Comments

1. 21–22

Report change from 1954-2000 and 2000-2024 to demonstrate that there is a clear acceleration between the two time periods, rather than the overlapping window.

Response:

Revised. The abstract and introduction now report non-overlapping periods (1954–2000 and 2000–2024/2025) to clearly demonstrate acceleration in mass loss, rather than the overlapping window used previously.

1. 36

Suggested rephrase to “The relative contribution of glacier melt is…”

Response:

Agreed, sentence rephrased accordingly.

1. 39

What do “small” and “significant” mean in this context? Quantifying will help to sell the importance of glaciers in this region.

Response:

We have replaced the qualitative descriptors with quantitative values in the revised text, reporting specific glacier area, meltwater contribution estimates, and downstream population exposure where data are available.

1. 42, 49, 50, 73, 81

“glaciated” > “glacierized”; repetition of “Fig. 1B”; provide values in km².

Response:

“Glacierized” is used throughout. The repeated “Fig. 1B” reference has been deleted. Glacier areas are now reported in km² at the relevant locations.

Figure 2

Why not also produce mass balance calculations for 2000–2024 for the glaciers in the UAV survey? Note: this Fig references 18 glaciers in the UAV survey, but in the abstract this was 11 glaciers?

Response:

Mass balance for 2000–2024 was computed for the UAV-surveyed glaciers and is reported in §4.3 and Fig. 5; Figure 2 is a methodological flowchart and did not previously include this step explicitly, which was an omission now corrected. The reference to 18 glaciers in Figure 2 was a labelling error from an earlier version of the analysis. The UAV survey covered 11 glaciers: 1 in the Cerro Blanco subcomplex (Glaciar Nevado, CL108116004) and 10 in the Las Termas subcomplex, based on the DGA 2000 glacier inventory. This is consistent with the abstract and has been corrected in Figure 2 and throughout the manuscript.

1. 157

Did you perform any vertical transformation or georeferencing to ensure that the PSAD56 datum was comparable to the vertical SRTM and UAV-based datums?

Response:

Yes — see response to RC1.1 above. We now apply pixel-wise EGM96 geoid corrections to all ellipsoidal DEMs (UAV, Pléiades, ASTER) prior to co-registration and differencing. The PSAD56 datum of the IGM DEM approximates EGM96 to within ~1 m at 37°S, which we consider negligible given the IGM DEM’s overall uncertainty. This is described in §3.3.1, §3.3.3, and §3.3.4.

1. 175

While true, have you considered and accounted for the ice penetration bias that might be present in SRTM data (see Gardelle et al. 2017)?

Response:

See response to RC1 l. 171–179 above. C-band penetration into wet ice is negligible (<1 m; Rignot et al., 2001), and SRTM was acquired during the late austral summer ablation season (February 2000) when glacier surfaces are predominantly wet ice. Any residual penetration into firn would cause a conservative (underestimated) mass loss signal. This limitation is now acknowledged in §3.3.2.

1. 189, 195

Is there a risk that by using polygons from the year 2000 to co-register a DEM from 1954, that you are unintentionally adding in unstable terrain? Do you have glacier boundaries from 1954?

Response:

See RC2.1 response above. We do not have 1954 glacier boundaries (outlines are not identifiable as a distinct feature on the IGM topographic map, as noted in §3.3.1) but now incorporate the 1975 inventory to substantially reduce this issue. The 1975 outlines are used in the glacier union mask for co-registration (§3.4) and as the reference area for mass balance calculation in 1954-based epochs (§3.5).

1. 211

Where do these control points come from?

Response:

No ground control points (GCPs) were used. Image georeferencing relied entirely on the onboard GPS-RTK positioning system of the DJI Mavic 3 Enterprise (RTK accuracy: 1 cm + 1 ppm horizontal, 1.5 cm + 1 ppm vertical), as described in §3.1. SfM processing was performed in Agisoft Metashape Professional (v.2.1.2). For Las Termas, 1,046 of 1,046 images were successfully aligned (RMS reprojection error: 0.22 px; coordinate precision: 4.31 cm). For Cerro Blanco, 1,095 of 1,110 images were aligned (RMS reprojection error: 0.26 px; coordinate precision: 2.50 cm). These processing statistics have been added to §3.2, and the manuscript text has been clarified to note that RTK-based positioning was used in lieu of GCPs.

1. 221

It is not necessary to divide by the number of years for the error estimates, as you aren’t assessing change over time of the stable terrain, but rather identifying what the difference of this (non-changing) terrain is and using that as the error.

Response:

Agreed. The previous approach of annualizing the NMAD was incorrect and has been removed. We now report stable-terrain NMAD as a direct measure of elevation error, and propagate this through the spatial aggregation using the Rolstad et al. (2009) equation as described in RC1.2. The error is not divided by the number of years.

1. 231

Why not use the average temperatures from the dataset rather than manually deriving from minimum and maximum values?

Response:

We have updated the temperature analysis to use daily mean temperatures computed directly from the dataset, rather than averaging minimum and maximum values. The revised results are consistent with the submitted version, confirming that this change does not materially affect our conclusions. The methods text has been updated accordingly in §3.6.

1. 287

Topographic variables*. Which statistical test did you use here?

Response:

Corrected to “topographic variables” throughout. We now report both Spearman rank correlation (primary, as our elevation change distributions are non-normal) and Pearson correlation coefficients for all glacier-level and pixel-level topographic analyses. The statistical tests are now explicitly described in the new §3.7 (Geospatial and Topographic Analysis) and results are reported in §4.2.

Fig. 5

Is the bar missing for Las Termas in the 2000 outlines? Better to plot 1954–2000, then 2000–2024, then 1954–2024.

Response:

The missing bar was a plotting error and has been corrected. The figure has been reordered to show 1954–2000, 2000–2024/2025, and 1954–2024/2025 sequentially, which more clearly demonstrates the acceleration in mass loss between the earlier and later periods.

1. 409, 425, 430

Also Hugonnet et al. (2019); can you describe (quantify) these trends and similarities to your data?; topographic variables.

Response:

Hugonnet et al. (2021) is already cited and the regional comparison is quantified in §5.1; Hugonnet et al. (2019) has been added to the reference list and cited alongside the other regional studies in §5.1. Trends and similarities with our mass balance results are now quantified in §5.1. “Geographic” has been corrected to “topographic” throughout.

1. 445–454

This section about peak water having passed could be strengthened with further reference to previous studies that have quantified this effect in the region – e.g. McCarthy et al. (2022) or Ragettli et al. (2016).

Response:

We have added references to McCarthy et al. (2022) and Ragettli et al. (2016), which provide regional context for peak water in central Chile, and have expanded the discussion of peak water in §5.1 to situate our streamflow findings within this broader literature.

1. 465

Is it not possible to identify this from UAV images?

Response:

Distinguishing debris-covered ice from recently deglaciated terrain in UAV imagery is challenging at this site due to similar spectral and textural properties. We are developing machine learning approaches for automated classification in an upcoming study, and have noted this in §5.1.

1. 475

Given that glacier area change isn’t reflected in the calculations, is it more accurate to report glacier thinning estimates than mass balance change estimates?

Response:

See RC2.1 response above. The incorporation of the 1975 inventory and the use of a time-weighted mean area approach following Cogley et al. (2011), described in §3.5, substantially addresses this concern. Glacier area change is now explicitly accounted for in the mass balance calculation for all epochs, and we consider reporting as specific mass balance (m w.e. yr⁻¹) to be appropriate given these revisions.