the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 20 May 2026
Lagged hydrological responses to glacier surface variability inferred from satellite-derived proxies in the Andes of central Chile
Abstract. Glaciers sustain downstream water resources by providing delayed meltwater inputs, but ongoing climate change is rapidly altering their dynamics and associated hydrological responses. Understanding how glacier variability propagates through hydrological systems remains a key challenge, particularly in data-scarce mountain regions. This study characterises glacier surface dynamics and quantifies their lagged influence on streamflow, and explores associated vegetation responses, across glacierised catchments in central Chile. Multi-decadal Landsat observations were used to derive physically based glacier proxies, including snow-covered fraction, albedo, land surface temperature (LST), and snow elevation. Principal component analysis was applied to decompose glacier behaviour into dominant modes of long-term change, short-term variability, and temporal persistence, while prewhitened cross-correlation analysis was used to assess lagged relationships with downstream responses. Physical consistency among glacier proxies was evaluated to ensure coherent representation of glacier surface conditions. Results reveal a pervasive signal of glacier degradation, with most glaciers exhibiting declining trends in snow-covered fraction and albedo alongside increasing LST. Variability shows strong scale dependence, with contrasting glacier-level and regional trends driven by the influence of large glaciers, whereas temporal persistence exhibits a consistent increase across scales. Lagged analyses show a transition from glacier-dominated headwaters, where glacier variability precedes streamflow, to downstream basins where hydrological responses become increasingly decoupled from cryospheric processes. A regression framework linking glacier-response modes to hydrological lag explains a substantial fraction of variability in streamflow timing (R² ~0.72). Vegetation responses are spatially heterogeneous and not significant at the system scale. These findings demonstrate the potential of satellite-based approaches to quantify the timing, persistence, and propagation of cryospheric signals in data-scarce mountain regions.
- Preprint
(13186 KB) - Metadata XML
-
Supplement
(3220 KB) - BibTeX
- EndNote
Status: final response (author comments only)
- RC1: 'Comment on egusphere-2026-2386', Anonymous Referee #1, 16 Jul 2026
-
RC2: 'Comment on egusphere-2026-2386', Anonymous Referee #2, 17 Jul 2026
Summary
The authors build a Landsat-based workflow in Google Earth Engine to derive four glacier-surface proxies (snow covered fraction (SCF), broadband albedo, land surface temperature (LST), and the median elevation of snow-covered pixels). They include 174 valley and mountain glaciers in the Valparaíso and Metropolitana regions of central Chile. Snow is classified with an adaptive Otsu threshold applied to a composite index calculated as snowScore = 0.7 x NDSI + 0.3 x VIS. The time series are deseasonalised, summarised by trend, rolling standard deviation ("variability") and lag-1 autocorrelation, and reduced to three PC1 indices. Prewhitened crosscorrelation functions (CCFs) are used to identify the lag of maximum correlation between catchment-aggregated glacier proxies and monthly streamflow (and NDVI). A multiple regression of that lag on the three PC1 scores yields R² ≈ 0.72.
The paper addresses an under-researched question, namely, how does a cryospheric signal propagate, in time, into downstream discharge? The GEE processing chain is thoughtfully constructed (topographic correction, adaptive thresholding, observation-corrected normalisation), the study region is well motivated, and the decomposition into trend/variability/memory modes is interesting.
However, in its current form I do not think the central results are established at the level of confidence the manuscript claims, and I have concerns about several of them being artefacts rather than signals. As a general comment, there is insufficient detail when it comes to the workflows creating the input remote sensing datasets, and insufficient detail regarding the hydrometric data. In addition, the statistical workflow has flaws that should be addressed prior to publication.
I believe that if these comments are addressed, this paper will be much stronger and an excellent addition to the published literature on this topic.
My comments are summarized below in Major and Minor Comments:
Major Comments- Static glacier outlines
The analysis is clipped to 2022 glacier outlines. Yet the analysis stretches back to the mid-1980s and the authors report area changes over time as well as other glacier statistics. I assume that the glaciers extended beyond the 2022 outlines, in the earliest parts of the timeseries. How is this workflow limitation currently being accounted for? This limitation could be overcome by time-series glacier outlines, or by clipping to glacier outlines from the earliest time step with the assumption that glaciers are not expanding beyond that polygon. - Water and rock classification.
Throughout the work, there is no mention of other land cover types in the glacier outlines, of particular note are: bare ground / rock and water. Water, in particular is often misclassified as snow/ice in glacier remote sensing. This is visible in the figures shown, proglacial lakes are clearly being mapped as snow/ice. - Issues with the trend analysis.
The reported trends are not defended against cross-sensor inconsistency (TM/ETM+/OLI band passes, cloud screening skill, SLC-off) and a roughly two-fold increase in observation density after 1999. The manuscript's only defences of multi-sensor consistency are (i) visual inspection of two glaciers (L265–271) and (ii) stability of the Otsu threshold (L272–276). L270–271 then states outright that the patterns "are not sensitive to sensor transitions or data gaps." That claim is not supported by anything presented, and it is doing an enormous amount of load-bearing work. The LST trend in particular (+0.168 K yr⁻¹, or ~+6.9 K over the record!) is not physically credible for a melting glacier surface and is, to me, a red flag for the whole proxy chain. The reported changes in "variability" and "persistence" are exactly what a change in sampling density would produce, and no test excludes this. Is the temperature being reported for the entire glacier polygon (i.e. including the exposed rock surfaces?). - The "physical consistency assessment" is not a validation.
The proxies are entangled by construction, so their mutual correlations are to be expected and do not provide an independent assessment of the results. There does not seem to be any external validation anywhere in the paper, despite regional benchmarks being available. - The validity of R2 = 0.72
The regression's response variable is the lag where each station's prewhitened glacier–streamflow cross-correlation happens to peak (the argmax of 25 noisy correlations from lags −12…+12), where even the winning r is only 0.15–0.25. Taking the maximum of 25 estimates is biased toward finding something, and under a true null the location of that peak is roughly uniform, so y itself is potentially mostly noise. The authors' F = 9.21 (p = 0.002) in Table 2 does guard against overfitting, so I don't dispute the R² on those grounds. The problem is that the F-test's validity rests on assumptions that are all violated here, each in the same optimistic direction: the ~15 catchments aren't independent (nested basins share the same glaciers, so the predictors are built from overlapping data and the effective n is smaller than 15); y is a derived argmax that can couple to variables like record length, catchment area and elevation (shorter records give noisier CCFs and more random peaks; the PC scores themselves correlate with glacier elevation and size in Fig. 9); and the whole y vector is the end of a long chain of defensible-but-arbitrary choices. So, the reported p-value is almost certainly too small, and nothing in the paper can currently distinguish R² = 0.72 from chance. Perhaps a model sensitivity test could help here: generate surrogate glacier timeseries that keep the autocorrelation but are independent of streamflow, e.g. via a fitted AR/ARMA model (cleanest here, since it reuses your rewidening step). Then run them through the entire pipeline — aggregation, prewhitening, cross-correlation, lag-picking, PCA, regression — a few thousand times and see how often you get an R² as high as 0.72 by pure chance. If almost never, the result is real and the paper is stronger for showing it; if it happens even 10–15% of the time, the main findings need reframing. - Snowline metrics
"Physically based glacier proxies" (L5, abstract). NDSI, VIS-brightness and an Otsu threshold are empirical/spectral, not physically based. Please soften this statement. "Snow elevation" is calculated as the median elevation of snow-covered pixels. This is uncommon in the literature and is not the snowline altitude. The median is confounded by hypsometry and by the snow-covered fraction. The field convention (Rastner et al., 2019; Racoviteanu et al., 2008 — and many more recent regional papers) is the transient snowline altitude. Note also that the reported trend (+1.64 m yr⁻¹ ⇒ ~+67 m over 41 years) is small relative to published SLA rises in the Andes, which is worth discussing. - Glacier volume and "water-equivalent-weighted" data
It is unclear how glacier volume and snow volumes were calculated and/or used. Volume and thickness are presented in the results and discussion, but not earlier in the paper. The aggregation of the results as “water-equivalent-weighted” data not sufficiently explained. The paper's own conclusion about scale dependence (L11–13, L464–466) is a direct consequence of this weighting choice, not an independent finding. Regional trends being "driven by the influence of large glaciers" is what a volume weighting does. Please present the area-weighted and unweighted alternatives alongside, and reframe this as a sensitivity of the aggregation rather than as a result. Weighting a snow-cover fraction by ice volume is an unusual choice. Area weighting would seem more defensible. Please justify. DGA inventory "volume" and "mean thickness" are almost certainly derived from area via volume–area scaling. If so, then the Area / Volume / Mean thickness rows of Fig. 9 are deterministic functions of one another, and presenting them as three separate structural controls is misleading. Please state how volume and thickness were obtained; if they are scaling-derived, collapse the rows and say so. - The vegetation component
The NDVI analysis returns a null result at every level (L365–369, L382–383, L404–410), and I think this is honest and correctly reported. But it occupies abstract real estate in the abstract, objectives, a figure, and several Discussion paragraphs, for a finding that amounts to "we looked and found nothing." Two concerns beyond that: PKU GIMMS NDVI is based on AVHRR (~9 km resolution), half-monthly, and ends in 2022. In narrow Andean catchments a 9 km pixel is insufficient; the product is also partly MODISfused/BPNN-modelled post-2000, so its pre/post-2000 interannual variability is not a clean AVHRR record. Given that the authors already have 30 m Landsat reflectance for the entire archive, using Landsat NDVI (or at minimum MODIS at 250 m) would be strictly better and essentially free. Please justify the choice or switch. - Hydrology data
The hydrometric data shared is barely explained with regards to: list of stations, are these stations correlated to eachother, basin size, other watershed characteristics that could be useful, which ones are nested? - Reproducibility does not meet TC requirements
L474–475 states the code "will be made available." For The Cryosphere, code and data must be archived (Zenodo or equivalent, with a DOI) and accessible at the time of review. A bare GitHub URL is not an archive. There is also no Data Availability section at all, which is mandatory. Please add one covering: the derived glacier-scale time series, the catchment aggregations, the station list, and the GEE scripts.
Minor comments
Abstract
L5–6: "physically based glacier proxies" should be “remotely sensed observations” or “remotely sensed proxies”
L11–13: Scale dependence is likely an artefact of the w.e. weighting choice (see major comments).
L15–16: "R² ∼ 0.72" — state n = 15 here. As written, this reads as a large-sample result.
L16: Consider whether the vegetation null belongs in the abstract at all, also see major comment regarding the vegetation analysis.
Introduction
L20–21: Owen et al. (2009) and Winkler et al. (2010) are unusual anchors for the claim about glaciers as sensitive Earth-system components. IPCC SROCC / Hock et al. (2019) / Zemp et al. (2019) would be more standard.
Study area
L80–83: Please give total glacierized area, the area distribution of the 174 glaciers (min/median/max), their elevation range, and — importantly — the percent glacierized of each study catchment. Without this the reader cannot judge whether a glacier signal in streamflow is even plausible at each station.
L83: Please identify the basins (Maipo and Aconcagua) and regions (Valparaiso and Metropolitana) on the map in Figure 1 as well as other important location features such as major centers for reference.
Datasets
L100–106: Static outlines. State the DGA classes retained and the minimum glacier size. Reconcile "174 glaciers" with the "Rocky and small glaciers" class in Fig. 5.
L108–111: State the Landsat collection and tier explicitly (Collection 2 Level-2? Tier 1 only?).
L109 vs L141 vs Table 1: "mid-1980s to 2025" (L109) vs "1984–2026" (L141) vs Table 1's last acquisition of 26-12-2025. Please correct.
L113–115: NASADEM has a ~2000 epoch. Over 41 years of thinning (tens of metres in places) this introduces a small but non-zero bias in the snow-elevation metric and in the illumination correction. Please note it.
L117–121: How many gauging stations? This is never stated anywhere in the paper. Provide a table (station ID, name, river, catchment area, % glacierised, record length, n overlapping months, selected lag, r at that lag, whether nested within another).
L120: "at least 240 overlapping monthly observations" — contiguous, or with gaps?
Preprocessing
L137: "<70 % valid glacier coverage" excluded. How does this impact the overall results? Is this based on literature or was a sensitivity analysis completed?
Snow quantification
L146: Can you comment on if the NDSI is good at differentiating snow vs ice vs water vs rock.
L150–152: The VIS term is "rescaled to the 0–1 range" — how, and over what domain? If min–max rescaled per scene or per glacier, the composite metric becomes scene relative and temporally non-comparable, which would undermine the entire time series.
L155: The 0.7/0.3 weights are asserted, not derived or referenced to previous research on the topic. Please provide a sensitivity analysis (e.g. weights from 1.0/0.0 to 0.5/0.5) showing the effect on SCF and, crucially, on the SCF trend.
L160–164: Over what histogram is Otsu applied — per glacier per date, per scene, or pooled? State it. Report the minimum pixel count and how Otsu behaves for the smallest glaciers.
L165–167: The fallback rule is ambiguous: when is 0.3 used and when is 0.45? Please also report the fraction of glacier-dates that used a fallback threshold, broken down by year and by sensor. If fallback frequency has a trend (it likely does, since it triggers under near-full or near-absent snow), that alone can generate a spurious SCF trend.
L169–171: Shadows are flagged but not excluded from snow classification. Cast shadow over snow gives low VIS but still-high NDSI; the composite metric is then pushed toward the threshold in a way that depends on solar geometry, which varies with acquisition date, latitude, and sensor overpass time. Please justify, and test the sensitivity of SCF trends to excluding shadowed pixels.
L179–181: Which narrowband-to-broadband coefficients, for which sensor? Liang (2001) was developed for ETM+/MODIS; Traversa et al. (2021) provide OLI coefficients. Applying one set across four sensors would introduce a step change. Any compensation for solar zenith angle?
L182–183: LST — see major comments.
L184–186: Median snow elevation is not typically referred to as snowline altitude (see major comments).
L185: Mention limitations of SRTM based DEMs on glacier elevation.
Deseasonalisation and trends
L214: The phrase "water-equivalent-weighted aggregation to account for glacier size" (L214–215,
223–225) appears repeatedly but is never defined. Weights = area × mean thickness from the inventory? Volume? Something else?L216–218: The 84-month window is asserted without justification. Please explain.
Lagged analysis
L222–226: Nested catchments share glaciers, i.e. pseudoreplication. Report the number of glaciers per catchment and the overlap matrix.
Results
L265–271: "These patterns are consistent across the full Landsat archive and are not sensitive to sensor transitions or data gaps." This is an unsupported claim resting on visual inspection of two glaciers. Either substantiate it quantitatively or delete it.
L293–296: The admission that the SCA trend "is sensitive to data completeness" needs to be turned into a table: trend vs. valid-area threshold (70/80/90/95/99 %). This is important because the reader currently cannot tell how much of the result is a QC choice.
L315–316: "except for LST at the regional scale" — this exception is passed over. Why does LST behave differently?
L352–355: "climatic forcing appears to dominate long-term glacier change" is inferred from the weakness of the correlations with glacier properties. Please soften this statement as there it is an argument in absence.
L379–381: How many basin clusters? If fewer than ~10, the cluster-robust result is not informative and should not be presented as reassurance. The authors half acknowledge this in the following clause; please state it plainly.
L384–387: The "trade-off between snowpack-mediated storage and event-scale dynamical variability" is a narrative fitted post hoc to two regression coefficients. It is not tested. Please label it as a hypothesis.
Discussion
L396–398: Direct contradiction of L240–241.
L419–420: Uncertainties are acknowledged in one sentence but nowhere quantified. There are no error bars on any glacier-month proxy value anywhere in the paper. Please propagate at least a first-order uncertainty (classification error × pixel count × valid fraction).
L425–441: The water-management implications run to ~17 lines and reach as far as ice stupas and snow fences. This is disproportionate to an evidence base of |r| ≈ 0.2 correlations in ~15 catchments. I would cut this by half and remove the adaptation technology paragraph.
L443–449: The limitations section is candid but incomplete. It should add: crosssensor consistency, time-varying observation density, static glacier outlines, multiple testing in the lag selection, pseudoreplication in the regression, the absence of external validation, and the absence of a precipitation control.
Back matter
L474–475: Code "will be made available" — must be archived with a DOI now. Missing: A Data Availability section (mandatory for TC).
L476–477: Author contributions are internally inconsistent: "I.F. and S.M. conducted the analysis... with guidance from I.F. and N.V." I.F. cannot guide I.F. Please correct.
L480: Ironically, in the reference to ChatGPT being used for language improvement, there are typos: “thepurpose” and "fullresponsibility"
Reference list
Many mistakes: "hess opinions", "the camels-cl dataset", "chile", "landsat", "modis", "ndsi", "andes", "greenland ice sheet", "köppen-geiger", "arima models", etc. In-text citation error: "Colin Cameron and Miller (2015)" (L262) should be "Cameron and Miller (2015)" — the author is A. Colin Cameron; the surname is Cameron. "Jpl, N. (2020)" should be "NASA JPL (2020)". Same in the Fig. 1 caption ("Jpl, 2020"). Many more were found. Please review carefully.Figures
Figure 1. Panel (d) label typo: "Anual rainfall" should be "Annual rainfall". Please add gauging-station labels/IDs to panel (b), and mark which catchments are nested.
Figure 3. Only two glaciers shown. Water is clearly being mapped as snow. Given the weight this figure carries (it is the paper's only classification evaluation), please show a larger, systematically-selected sample spanning glacier size, aspect, and debris cover — and, critically, include failure cases.
Figure 4. Caption typo: "dhashed" should be "dashed".
Figure 5. Legend entry "Rocky and small glaciers" is undefined in the methods. Font sizes are too small.
Figures 10 and 11. The CCF insets are unreadable at print size. More substantively: they show only six of the (~15?) stations, chosen without stated criteria. Show all stations, in a supplementary panel grid.
Table 1. Please add scenes per year (or a time-series plot of observation density in the supplement).
New table required. A station/catchment table (see comment on L117–121). Could be added to supplement.
Citation: https://doi.org/10.5194/egusphere-2026-2386-RC2 - Static glacier outlines
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 173 | 72 | 17 | 262 | 29 | 14 | 13 |
- HTML: 173
- PDF: 72
- XML: 17
- Total: 262
- Supplement: 29
- BibTeX: 14
- EndNote: 13
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Peer review for “Lagged hydrological responses to glacier surface variability inferred from satellite-derived proxies in the Andes of central Chile” by Fuentes et al, submitted to the Cryosphere
This manuscript presents a satellite-based analysis of glacier surface variability and its lagged relationships with downstream streamflow in central Chile. I found the central idea interesting and potentially innovative, particularly the effort to move beyond conventional trend analysis by examining variability, temporal persistence, and lagged hydrological responses. The manuscript is generally well structured, and the results provide interesting insights into spatial variability among glaciers and differences between glacier-scale and regional behaviour.
My main concern is the validation of the snow classification. The distinction between snow, firn, and exposed glacier ice is central to the analysis, but the current assessment relies largely on visual inspection and consistency among satellite-derived proxies. Some examples in Figure 3 do not appear consistent with the RGB imagery, which reduces my confidence in the snow-covered fraction and the subsequent analyses based on it. I am not suggesting that this necessarily changes the main results, but the classification accuracy and associated uncertainty need to be demonstrated more convincingly.
I also recommend tightening the manuscript’s scope. The vegetation analysis is presented from the beginning as secondary, produces limited results, and distracts from the stronger cryosphere–hydrology narrative. It should either be integrated more clearly into the study’s central questions or reduced substantially. Finally, the physical interpretation of the lag analysis should be more cautious: cross-correlations demonstrate temporal associations, but do not necessarily establish direct glacier contributions or hydrological decoupling. My review focuses primarily on the remote-sensing methodology, physical hydrological interpretation, and overall framing; I leave detailed assessment of the time-series statistics to reviewers with greater statistical expertise.
Overall, I think the study has strong potential, but some flaws need to be addressed before publication.
Line-by-line comments
Abstract: Vegetation is mentioned in the abstract but not in the title. More generally, its role in the manuscript is unclear. It is introduced as secondary from the outset and does not contribute strongly to the main narrative. I suggest either integrating vegetation more fully into the research questions or substantially reducing this component and maintaining a clearer cryosphere–hydrology focus.
Lines 80–97: More context on the regional hydroclimatic regime would be helpful. Please distinguish the lower-elevation winter rainfall-runoff response from high-elevation winter snow accumulation and subsequent spring–summer melt. The strong topographic influence on precipitation phase and water storage appears central to the lagged responses investigated here.
Line 96: Does “climatic gradient” refer to seasonality, elevation, or spatial variations in precipitation and temperature? This sentence does not flow naturally from the preceding discussion. Consider reorganising the section to move from the regional importance of water resources to climate, topography, and the resulting runoff regime.
Lines 104–106: Are debris-covered glaciers included? The manuscript states that rock glaciers were excluded, but the treatment of debris-covered glacier ice is unclear.
Lines 133–136: Please describe the preprocessing in sufficient detail to make it reproducible. What cloud and shadow masks, thermal filters, topographic corrections, and cross-sensor harmonisation procedures were applied?
Lines 137–139: Please report the number of glaciers here, in addition to the number of scenes. It would also be useful to report the typical number and temporal distribution of observations available per glacier.
Table 1: Was Sentinel-2 considered as an independent, higher-resolution dataset for validating the classification during the overlapping period?
Lines 150–167: How were the 0.7/0.3 weights, percentile limits, and fallback thresholds selected? The statement that these choices ensure robust snow detection requires supporting validation.
Lines 169–171: Shadowed pixels are retained in the snow classification. Please demonstrate that the method can reliably distinguish snow and ice under shadowed conditions.
Lines 174–177: Snow-covered fraction is calculated relative to the valid observed area. If clouds or sensor gaps obscure a particular elevation band, this could artificially increase or decrease the estimated fraction. How was this potential bias assessed? Why not use the full glacier outline as the denominator and retain only scenes with sufficiently complete coverage?
Line 194: Please report area using SI-compatible units, preferably km² rather than hectares.
Section 2.5: The physical-consistency analysis evaluates agreement among variables derived from the same satellite observations. This is useful as an internal diagnostic, but it does not independently validate the snow classification. Please include an assessment of classification reliability, particularly for snow, firn, clean ice, debris-covered ice, rock, and shadow.
Lines 195–211: Section 2.5 describes removal of a glacier-specific mean seasonal cycle, whereas Section 2.6 describes harmonic regression. Please clarify whether these are the same procedure or two separate deseasonalisation methods.
Lines 184–185 and elsewhere: This metric is the mean or median elevation of snow-covered pixels, not necessarily snowline elevation. Please use consistent terminology and interpret it cautiously, as it may also be influenced by glacier hypsometry and total snow extent.
Lines 214–225: Please define “water-equivalent-weighted” aggregation and explain how the weights were calculated.
Lines 228–240: Please state the lag-sign convention explicitly in the methods and relevant figure captions.
Figure 3: Some examples make me question the robustness of the snow classification. In the 23 January 2025 and 4 February 2009 scenes, parts of the lower glacier appear snow-free in the RGB images but are classified as snow. In the 3 December 1988 scene and several others, some apparently off-glacier areas also appear to be classified as snow. Please examine these discrepancies and provide a more rigorous validation of the classification.
Lines 265–270: Successfully retaining spatial coherence across Landsat 7 SLC-off gaps is a necessary baseline, but does not by itself demonstrate classification accuracy. I suggest focusing this section on evidence that the method correctly distinguishes snow from glacier ice and other surfaces.
Figure 4: The labels are very small. A seasonal pattern also appears to remain in the snow-elevation series. Please verify whether deseasonalisation was effective for this metric.
Lines 278–287: Is LST calculated across the entire glacier outline or only over snow-covered pixels? Please state this clearly in the results or caption.
Section 3.4: The scale dependence of variability is one of the most interesting and novel results. I suggest moving Figure S4 into the main manuscript, as the contrast between regional aggregation and individual-glacier behaviour deserves greater emphasis.
Figure 5: This figure contains useful information but is difficult to read. The labels and glaciers are very small, and the significant, non-significant, and excluded categories are difficult to distinguish. Consider enlarging the maps, zooming more closely to the analysed glaciers, and simplifying the colour scale and legend.
Figure 6: Panel A could potentially be moved to the supplement. Panel B should be described as an assessment of internal physical consistency rather than independent validation.
Figure 7 and supplementary figures: Please use consistent terminology for “memory” and “temporal persistence.” Currently, both terms are used for lag-1 autocorrelation.
Figure 8: Because the PCA is constructed from these metrics, strong correlations between the resulting PC scores and the original input variables are expected. Please clarify what additional result this figure provides. It may be more appropriate to report the PCA loadings directly and move this figure to the supplement.
Figures 7–9 and others: Figure titles do not need to be repeated within the figures; this information can remain in the captions.
Lines 352–355: This is a clear and effective synthesis of the PCA results.
Lines 357–364: I encourage the authors to distinguish more carefully between statistical association and physical causation. In particular, the interpretation that negative lags demonstrate direct glacier contributions, while positive lags indicate downstream decoupling, appears stronger than can be supported by cross-correlation alone. Please discuss other possible explanations, including shared climatic forcing and catchment storage.
Figure 10: The labels and point colours are difficult to read. Please enlarge the station symbols and plot them above the river network. The caption should also state clearly which glacier proxy was used to determine the mapped lag.
Lines 365–369: The vegetation results appear abruptly here. This reinforces my broader concern that the vegetation analysis is not well integrated into the main narrative.
Figure 12: This figure may not be necessary because the coefficients and their uncertainty are already reported in the text and Table 2.
Discussion: Subheadings would improve readability. I would also like to see more discussion of the novelty and advantages of the lag-analysis framework itself. What insight does it provide that conventional trend, correlation, or autocorrelation analyses would not?
Discussion: Please be cautious when interpreting statistical associations as evidence that glaciers “actively regulate” hydrological timing. The potential roles of shared climatic forcing and catchment storage should be discussed more explicitly.
Figure S2: Consider adding a horizontal dashed line at 273.15 K. Glacier and snow-surface temperatures substantially above the melting point would require explanation, including possible mixed pixels or LST-retrieval uncertainty.
Figure S3: The decline in short-term variability during approximately 2004–2024 is interesting. Could it be connected to the Central Chilean megadrought or another regional phenomenon? Some series also appear oscillatory rather than linear, so please consider whether a linear trend is an appropriate description.
Supplementary figures: Please improve the readability of all figures. Several labels are too small to read at normal viewing size.