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
-
AC1: 'Reply on RC1', Ignacio Fuentes, 11 Aug 2026
Comment 1: 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.
Response: We sincerely thank the reviewer for the careful evaluation of our manuscript and for the positive assessment of its novelty, structure, and overall contribution. We appreciate the recognition of our effort to move beyond conventional glacier trend analyses by investigating variability, temporal persistence, and lagged hydrological responses. We have carefully considered all of the reviewer's comments and believe that the revisions outlined below will substantially strengthen the manuscript.
Comment 2: 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.
Response: We thank the reviewer for this comment and fully agree that the original manuscript did not provide a sufficiently rigorous assessment of the snow-classification accuracy. We agree that the reliability and uncertainty of the snow-covered fraction estimates should be demonstrated more convincingly, as these variables underpin all subsequent analyses.
In the revised manuscript, we will incorporate an independent quantitative validation of the snow-classification procedure using Sentinel-2 imagery as a higher-resolution reference dataset. Validation of long-term glacier snow-cover products is inherently challenging because no independent dataset directly measures snow-covered fraction over the full Landsat archive. Existing in situ observations (e.g., snow water equivalent or snow depth) quantify snow accumulation rather than snow extent, while regional snow products such as MODIS have substantially coarser spatial resolution and are themselves derived from spectral snow indices, making them unsuitable as an independent reference for many of the glaciers analysed in this study.
To address this limitation, we will adopt a validation strategy based on independently interpreted Sentinel-2 imagery. Seven representative glaciers spanning a range of glacier morphologies and elevations will be evaluated using three Landsat acquisition dates each, resulting in 21 independent validation cases. Only Landsat–Sentinel image pairs acquired within one day of each other will be retained, and Sentinel-2 cloud-contaminated pixels will be excluded using the Scene Classification Layer (SCL) before manual interpretation.
Approximately 200 manually interpreted snow/no-snow reference points will be generated for each validation case. To provide a conservative assessment, reference points will be sampled within the bounding rectangle of a 500 m buffer surrounding each glacier rather than being restricted to glacier surfaces alone, thereby including surrounding terrain (e.g., rock, debris, terrain shadows and occasional water bodies) that makes the evaluation more challenging than the subsequent glacier-only analyses.
The revised manuscript will include a new subsection describing the validation methodology, representative validation examples, confusion matrices, pooled and case-by-case classification metrics (overall accuracy, precision, recall, F1 score, balanced accuracy and Cohen's kappa), and a supplementary table summarising the performance of all validation cases. We believe these additions will provide a substantially more rigorous assessment of classification reliability and directly address the reviewer's concern.
Comment 3: 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.
Response: We thank the reviewer for this thoughtful and constructive comment. We agree that the manuscript benefits from a clearer focus on the cryosphere–hydrology interactions, and we will revise both the framing and the interpretation accordingly. We agree that the original manuscript gave excessive weight to the vegetation analysis relative to the primary cryosphere–hydrology objectives. In the revised manuscript, we will reorganise the Results and Discussion to clearly establish glacier–streamflow interactions as the central narrative, while presenting vegetation responses as complementary evidence of downstream propagation of glacier-related signals.
Second, we agree that the interpretation of the lag analyses should be more cautious. Throughout the revised manuscript, we will replace causal wording with terminology that more accurately reflects the observational nature of the analysis. In particular, we will avoid statements implying direct glacier contributions or hydrological decoupling based solely on cross-correlation analyses, and instead describe the results as lagged temporal associations that are physically consistent with, but do not by themselves demonstrate, underlying hydrological processes.Finally, we will revise the Introduction, Results and Discussion to strengthen the overall cryosphere–hydrology narrative and ensure that the central scientific questions remain focused on how glacier surface variability propagates through downstream hydrological systems. We believe these revisions will improve the overall coherence of the manuscript while preserving the broader perspective provided by the vegetation analysis.
Line-by-line comments
Comment 4: 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.
Response: We agree that the role of the vegetation analysis required clearer definition within the overall manuscript. Rather than expanding this component into a co-equal objective, we chose to strengthen the cryosphere–hydrology focus of the study while retaining vegetation as a complementary analysis that explores whether glacier-related signals extend beyond streamflow into downstream ecosystem dynamics.
In the revised manuscript, we will clarify this hierarchy throughout the Introduction, Objectives, Results, Discussion and Conclusions. The central scientific focus will remain on the propagation of glacier variability to downstream hydrological systems, while the vegetation analysis will be explicitly presented as a secondary line of evidence. The Results and Discussion will be reorganised accordingly, placing greater emphasis on the glacier–streamflow relationships and interpreting the vegetation results primarily in the context of the spatial limits of glacier signal propagation.
Although vegetation responses exhibited significant lagged associations in some catchments, they were considerably more heterogeneous than streamflow responses and were not significantly associated with the dominant glacier-response modes identified in the PCA analysis. Rather than weakening the manuscript, we believe this contrast provides useful insight into the extent to which glacier variability propagates through different components of the coupled cryosphere–hydrology–ecosystem system.
We believe these revisions will provide a clearer narrative centred on cryosphere–hydrology interactions while retaining the vegetation analysis as a complementary contribution that helps define the downstream extent of glacier-related signals.
Comment 5: 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.
Response: We thank the reviewer for this helpful suggestion and agree that the original description of the regional hydroclimatic setting did not sufficiently emphasise the elevation-dependent partitioning of precipitation and runoff generation.
In the revised manuscript, we will expand the description of the study area to explicitly distinguish the contrasting hydrological behaviour of low- and high-elevation catchments. Specifically, we will clarify that winter precipitation at lower elevations falls predominantly as rainfall, producing relatively rapid runoff responses. In contrast, precipitation at higher elevations is largely stored as seasonal snowpack, with additional storage in glacier ice, delaying water release until the spring and summer melt seasons. We will further emphasise that this elevation-dependent partitioning between rainfall, snow accumulation and meltwater release provides the physical basis for the lagged hydrological responses investigated throughout the study.
We believe these additions will strengthen the physical motivation of the study and provide clearer context for interpreting the spatial patterns of lagged glacier–streamflow relationships.
Comment 6: 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.
Response: We agree that the original text did not clearly specify what was meant by "climatic gradients" and that the description of the study area can be organised more logically.
In the revised manuscript, we will reorganise this section to progress from the regional importance of Andean water resources to the climatic and topographic setting, followed by the elevation-dependent partitioning of precipitation into rainfall and seasonal snow storage and the resulting mixed runoff regime. We will also replace the more general expression "climatic gradients" with a more explicit description of the spatial variability in elevation, temperature and precipitation that underpins the hydrological behaviour of the study region.
Comment 7: 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.
Response: We thank the reviewer for pointing out this ambiguity. Debris-covered glacier ice was included in the analysis, whereas rock glaciers and glaciarets were excluded. The original wording did not clearly distinguish between these two cases. To clarify this point, we will revise the Methods section to explicitly state that glacier polygons were retained irrespective of surface debris cover, while rock glaciers were excluded because they represent a distinct cryospheric landform with different geomorphological and hydrological characteristics. This revision removes the ambiguity regarding the treatment of debris-covered glaciers.
Comment 8: 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?
Response: We will substantially expand the preprocessing description in the Methods to improve reproducibility. The revised manuscript now specifies the Landsat Collection 2 Level-2 products used, the application of USGS scale factors, the QA_PIXEL cloud and cloud-shadow masking procedure, the glacier-specific coverage threshold, and the hybrid Minnaert/C topographic correction. We also clarify that no additional empirical cross-sensor harmonisation was applied to the reflectance-derived glacier proxies because cross-sensor consistency was explicitly evaluated using temporally matched observations, which showed high agreement and no evidence of systematic offsets in snow-covered fraction, glacier albedo or snow elevation (new Figure S1). During this evaluation, we identified a progressive warm bias in Landsat 7 thermal observations after 2020, consistent with the satellite's late-stage orbital drift. Consequently, Landsat 7 land surface temperature observations acquired from 2021 onward were excluded from all LST analyses. These additions improve the transparency and reproducibility of the processing workflow while providing quantitative justification for the preprocessing choices.
Comment 9: 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.
Response: We have revised the Methods to report not only the number of retained Landsat scenes but also the number of glaciers included in the analysis and the temporal sampling available for each glacier. Specifically, we now state that the dataset comprises repeated observations for 174 glaciers and report the median number of valid observations per glacier together with the interquartile range and overall range after quality filtering. To further illustrate the temporal support of the dataset, we have added a new supplementary figure (new Figure S2) showing the distribution of valid observations per glacier across the study period. These additions provide a clearer description of the temporal density and robustness of the glacier time series used throughout the analyses.
Comment 10: Table 1: Was Sentinel-2 considered as an independent, higher-resolution dataset for validating the classification during the overlapping period?
Response: Following this recommendation, we incorporated an independent validation of the Landsat snow classification using Sentinel-2 Level-2A imagery during the overlapping observation period. A new subsection entitled "Classification validation" will be added to the Materials and Methods (Section 2.4.4), describing the validation procedure. Seven glaciers representing a range of glacier characteristics were evaluated using three Landsat acquisition dates each (21 independent validation cases). For each case, the closest Sentinel-2 acquisition (maximum one-day separation) was selected, and approximately 200 snow/non-snow reference points were manually interpreted. Overall, 3,599 independent reference points were analysed. The corresponding validation results will be incorporated into a new Results subsection, including new Figure 3 and Table 2, demonstrating excellent classification performance (overall accuracy = 0.98; Cohen's κ = 0.96).
Comment 11: 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.
Response: We thank the reviewer for raising this point. We agree that the original manuscript did not provide sufficient justification for the 0.7/0.3 weighting coefficients, percentile limits, and fallback thresholds used in the SnowScore and adaptive-threshold procedure. These values were initially selected empirically during algorithm development, with greater weight assigned to NDSI because of its established sensitivity to snow, while retaining a smaller contribution from visible reflectance to incorporate the high reflectance of snow in the visible wavelengths. However, they were not formally optimised against an independent validation dataset in the original analysis.
In the revised manuscript, we will therefore use the newly developed Sentinel-2 validation dataset (21 Landsat–Sentinel-2 image pairs and 3,599 manually interpreted reference points) to evaluate the sensitivity of classification performance to reasonable variations in these parameters. This analysis will allow us to determine whether the selected parameterisation is supported by independent validation and whether classification performance is robust to alternative weighting and threshold choices. The final parameter values and corresponding sensitivity results will be reported explicitly in the revised Methods and Supplementary Materials.
Comment 12: 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.
Response: We agree that the treatment of topographic shadow requires explicit evaluation. Shadowed pixels were retained in the snow-classification procedure because masking them would systematically remove portions of steep, high-relief glacier surfaces and substantially reduce temporal coverage. Snow classification is performed after topographic correction of the optical reflectance data, whereas shadowed pixels are excluded from the glacier- and snow-albedo calculations. We acknowledge, however, that the original manuscript did not demonstrate classification performance specifically under shadowed conditions.
In the revised manuscript, we will therefore use the new Sentinel-2 validation dataset to evaluate classification performance separately for illuminated and topographically shadowed Landsat pixels. Validation points will be stratified using the terrain-illumination/shadow mask already generated by the processing workflow, and confusion matrices and classification metrics will be calculated independently for both illumination classes. This will allow us to quantify directly whether retaining shadowed pixels introduces a systematic loss of snow-classification performance. The Methods will also be revised to distinguish clearly between the treatment of shadows in the snow classification and in the albedo calculations.
Comment 13: 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?
Response: We agree that normalising snow-covered area by the valid observed area can lead to biased fractions when missing pixels are spatially structured, for example when clouds or sensor gaps preferentially obscure particular elevation bands. In the revised analysis, we will therefore evaluate snow-covered fraction relative to the full glacier outline and retain only scenes meeting a minimum usable-coverage criterion. We will additionally compare this formulation with the original valid-area-normalised fraction across a range of coverage thresholds to quantify the sensitivity of the resulting time series to the denominator choice. This change avoids interpreting spatially incomplete observations as fully representative of the glacier surface while preserving temporal coverage through an explicit scene-quality criterion.
Comment 14: Line 194: Please report area using SI-compatible units, preferably km² rather than hectares.
Response: We modified hectares to m2
Comment 15: 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.
Response: We agree that the physical-consistency analysis is an internal diagnostic and does not constitute an independent validation of the snow classification. We have therefore added a separate validation based on near-contemporaneous Sentinel-2 imagery, as described in the new "Classification validation" subsection. The classification developed in this study is intentionally binary, distinguishing snow-covered from non-snow-covered surfaces. It is not designed to resolve firn, clean ice, debris-covered ice, rock, and shadow as separate thematic classes. Accordingly, the independent validation was formulated consistently with the classification objective, using manually interpreted snow/no-snow reference points across 21 validation cases and 3,599 reference observations. The validation domain was deliberately extended beyond glacier polygons to include surrounding rock, debris, terrain shadows, and occasional water bodies, thereby testing the method under spectrally challenging conditions. We will revise Section 2.5 to make clear that the physical-consistency analysis is complementary to, rather than a substitute for, the independent classification validation. We agree that explicitly separating snow, firn, clean ice, debris-covered ice, and rock would be valuable for studies focused on glacier-surface facies, but such a classification is beyond the scope of the present analysis and will be discussed in the discussion section.
Comment 16: 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.
Response: We thank the reviewer for identifying this ambiguity. The two passages were intended to describe the same deseasonalisation procedure, but the wording in Section 2.5 was imprecise. In the analysis, deseasonalisation was performed consistently using harmonic regression with two annual harmonics for both the physical-consistency assessment and the trend/robustness analyses. The glacier-specific seasonal component estimated by the harmonic model was removed while retaining the longer-term trend and interannual variability. We will revise Section 2.5 to explicitly refer to the harmonic-regression procedure described in Section 2.6 and remove the misleading reference to subtraction of a mean seasonal cycle.
Comment 17: 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.
Response: We thank the reviewer for this clarification. We agree that the metric used in this study represents the median (or mean) elevation of snow-covered pixels rather than the snowline elevation itself. The original terminology was therefore imprecise. In the revised manuscript, we have adopted the more accurate terminology "median (mean) snow-covered elevation" throughout the manuscript and explicitly state that these metrics describe the vertical distribution of snow-covered areas across each glacier. We also clarify that they are influenced by both snow extent and glacier hypsometry and therefore should not be interpreted as direct estimates of snowline or equilibrium-line altitude.
Comment 18: Lines 214–225: Please define “water-equivalent-weighted” aggregation and explain how the weights were calculated.
Response: We thank the reviewer for pointing out that this terminology was insufficiently defined. In the revised manuscript, we will explicitly describe the weighting procedure. Glacier-scale proxy values were aggregated using glacier water-equivalent volume estimates (EQ_AGUAKM3) obtained from the Chilean glacier inventory as weighting factors. Specifically, regional and catchment-scale values were calculated as weighted means, where each glacier contributed in proportion to its estimated water-equivalent volume. This approach gives greater influence to larger glaciers, which contribute disproportionately to regional glacier storage and potential meltwater production. The corresponding equation and data source will be added to the Methods.
Comment 19: Lines 228–240: Please state the lag-sign convention explicitly in the methods and relevant figure captions.
Response: We agree that the lag-sign convention should be stated explicitly. In the revised manuscript, we will clarify in the Methods that negative lags indicate glacier variability preceding the downstream response (streamflow or vegetation), whereas positive lags indicate the opposite relationship. The same convention will also be stated explicitly in the captions of all figures presenting lag maps or cross-correlation results to avoid ambiguity.
Comment 20: 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.
Response: We thank the reviewer for these observations. We agree that some of the illustrative examples show commission errors outside glacier boundaries or local discrepancies between the RGB imagery and the binary snow segmentation. These examples highlight the limitations of any automated snow-segmentation approach and are consistent with the non-zero commission and omission errors reported in the independent validation.
It is important to note that the objective of the proposed method is to quantify snow conditions over glacier surfaces rather than to produce a general land-cover classification of the surrounding landscape. The validation was intentionally designed to be substantially more demanding than the operational application by evaluating the segmentation over the entire bounding rectangle of a 500 m glacier buffer, thereby including terrain shadows, rock outcrops, water bodies and other non-glacier surfaces. Consequently, the reported validation statistics represent a conservative assessment of algorithm performance.
To address this concern, the revised manuscript will include an independent validation against manually interpreted Sentinel-2 imagery comprising 21 glacier-date combinations and 3,599 reference points. This analysis yielded a pooled overall accuracy of 98%, an F1-score of 0.98 and a Cohen's κ of 0.96, while also explicitly documenting the residual commission and omission errors. We have clarified that old Figure 3 is intended to illustrate the temporal consistency of the workflow across the Landsat archive rather than to serve as the primary validation of the segmentation.
Comment 21: 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.
Response: We agree with the reviewer that visual consistency across the Landsat archive, including successful handling of Landsat 7 SLC-off scenes, does not constitute an independent assessment of classification accuracy. In the revised manuscript, classification accuracy is now evaluated separately using an independent Sentinel-2 validation dataset comprising 21 glacier-date combinations and 3,599 manually interpreted reference points, yielding a pooled overall accuracy of 98%, an F1 score of 0.98 and Cohen's kappa of 0.96 (Section 3.1, Figure 2, Table 2).
Accordingly, we have revised the text describing old Figure 3 to clarify that its purpose is to illustrate the temporal consistency and operational behaviour of the adaptive SnowScore segmentation across the Landsat archive rather than to validate classification accuracy. Evidence supporting the ability of the method to distinguish snow from surrounding surfaces is now provided by the independent Sentinel-2 validation.
Comment 22: 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.
Response: We thank the reviewer for this observation. We agree that the labels in Figure 4 are too small and will increase their size in the revised figure. We will also explicitly reassess the effectiveness of the deseasonalisation procedure for the snow-covered elevation metric. In particular, we will examine the monthly distribution and residual annual autocorrelation of the harmonically deseasonalised series and revise the implementation where necessary to ensure that the harmonic regression uses actual acquisition dates rather than assuming regularly spaced observations. The revised figure and analysis will therefore verify that the reported long-term snow-elevation trend is not driven by residual seasonal structure.
Comment 23: 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.
Response: We appreciate the reviewer's suggestion. The manuscript has been clarified to explicitly state that land surface temperature (LST) was calculated as the median temperature of all valid observed glacier pixels within each glacier outline, rather than only snow-covered pixels. In contrast, snow albedo was calculated exclusively over illuminated snow-covered pixels. This distinction is now stated explicitly in the Methods when describing glacier proxies and clarified in the corresponding Results and figure captions.
Comment 24: 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.
Response: We appreciate this suggestion and agree that the contrast between regional aggregation and glacier-specific behaviour is an important aspect of the study. Figure S4 has therefore been moved from the Supplementary Material to the main manuscript, and the corresponding Results section has been expanded to explicitly discuss the scale dependence of glacier variability and the differences between regional and individual-glacier responses.
Comment 25: 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.
Response: We appreciate the reviewer's suggestion. Figure 5 will be redesigned to improve readability by combining both glacier regions into a single map for each variable, substantially enlarging the mapped glaciers. The figure layout will be simplified, with a single histogram per variable, larger labels, and simplified legends, making both the spatial patterns and trend distributions easier to interpret.
Comment 26: 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.
Response: We agree that Panel B represents an assessment of the internal physical consistency among independently derived glacier proxies rather than an independent validation of the snow-classification procedure. The manuscript has been revised accordingly to consistently describe these analyses as assessments of physical consistency. In addition, Panel A will be moved to the Supplementary Materials to improve the overall presentation of Figure 6.
Comment 27: Figure 7 and supplementary figures: Please use consistent terminology for “memory” and “temporal persistence.” Currently, both terms are used for lag-1 autocorrelation.
Response: We thank the reviewer for this suggestion. We agree that the terminology was inconsistent. Throughout the revised manuscript, we will consistently use temporal persistence, quantified as lag-1 autocorrelation (AC1), to describe this property of glacier dynamics. References to "memory" will be replaced by "temporal persistence" (or "temporal persistence mode" in the PCA analysis), and figure titles, captions, and text will be revised accordingly.
Comment 28: 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.
Response: We appreciate this observation and agree that strong correlations between the principal component scores and the original variables are expected because the PCA was constructed from these metrics. The purpose of the figure was not to validate the PCA, but to facilitate physical interpretation of the principal components by illustrating which glacier proxies dominate each mode and their direction of association. To avoid redundancy, we will remove this figure from the main manuscript and move it to Supplementary materials and instead report the interpretation of the PCA loadings directly in the Results section.
Comment 29: Figures 7–9 and others: Figure titles do not need to be repeated within the figures; this information can remain in the captions.
Response: Titles from those figures will be removed in the revised version.
Comment 30: Lines 352–355: This is a clear and effective synthesis of the PCA results.
Response: We thank the reviewer for this positive assessment. We are pleased that this synthesis clearly conveys the interpretation of the PCA results.
Comment 31: 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.
Response: We agree that cross-correlation alone cannot establish causal relationships. The revised manuscript now distinguishes more clearly between the observed lag structure and its physical interpretation. Specifically, we now describe negative and positive lags as temporal ordering between glacier variability and streamflow rather than evidence of direct glacier control, and we explicitly acknowledge that the observed relationships may also arise from shared climatic forcing, catchment storage, routing effects, and other hydrological processes. The interpretation has therefore been revised to emphasise consistency with these mechanisms rather than causal attribution.
Comment 32: 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.
Response: We thank the reviewer for these helpful suggestions. The revised version of this figure will be improved by enlarging the gauging-station symbols, plotting them above the river network, increasing the symbol outline thickness, and improving the overall visual contrast. In addition, the caption and colour-bar label will be revised to state explicitly that the mapped quantity is the lag corresponding to the maximum statistically significant prewhitened cross-correlation between glacier snow-covered fraction and streamflow. The lag-sign convention will also be stated explicitly in the caption.
Comment 33: 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.
Response: We agree that the transition to the vegetation analysis was abrupt in the original manuscript. The revised text now introduces the vegetation results by explicitly motivating them as an extension of the lag analysis from hydrological to ecological responses. We also clarify that the vegetation analysis is necessarily restricted to vegetated catchments, resulting in reduced spatial coverage relative to the streamflow analysis. This revised structure better integrates the vegetation results into the overall narrative.
Comment 34: Figure 12: This figure may not be necessary because the coefficients and their uncertainty are already reported in the text and Table 2.
Response: We appreciate the reviewer's suggestion. Although the regression coefficients are reported numerically in Table 2, we will retain Figure 12 because it provides a concise visual comparison of the relative magnitude, direction, and uncertainty of the standardised regression coefficients, facilitating interpretation of the respective contributions of the three PCA modes. We will revise the caption to clarify this purpose.
Comment 35: 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?
Response: We appreciate this suggestion and agree that the Discussion benefits from clearer thematic organisation. We will introduce descriptive subheadings to improve readability. We will also expand the discussion of the lag-analysis framework to clarify its methodological contribution. In particular, we will explain that, unlike conventional trend, correlation, or autocorrelation analyses, the proposed framework explicitly quantifies the temporal ordering and response delays between glacier variability and downstream hydrological and ecological responses while accounting for serial dependence. This temporal perspective provides additional insight into how glacier changes propagate through mountain catchments and better motivates the novelty of the proposed approach.
Comment 36: 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.
Response: In the revised manuscript, we will adopt more cautious language to distinguish statistical association from physical causation. Expressions implying direct regulation (e.g., "glaciers actively regulate hydrological timing") will be replaced by wording indicating statistical association or consistency with plausible physical mechanisms. We will also expand the Discussion to explicitly acknowledge that the observed lag relationships may reflect the combined influence of glacier variability, shared climatic forcing, and catchment storage (e.g., snowpack, firn, groundwater, lakes and reservoirs).
Comment 37: 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.
Response: We will add a horizontal reference line at the melting point (273.15 K) in old Figure S2 to facilitate interpretation of the LST time series. We will also clarify in the Methods and the figure caption that glacier LST represents the median Landsat-derived radiometric surface temperature across the entire glacier polygon rather than the temperature of snow or clean glacier ice alone. Consequently, values exceeding 273.15 K may occur because glacier outlines include debris-covered ice, exposed rock, mixed pixels and partially snow-free surfaces, as well as uncertainties inherent to satellite-based thermal retrievals.
Comment 38: 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.
Response: We thank the reviewer for this observation. During the revision process, the glacier time series were recalculated following improvements to the masking procedure and quality-filtering criteria introduced elsewhere in the revised manuscript. These refinements affected the rolling variability and temporal persistence estimates, and the previously apparent decline in short-term variability around 2004–2024 is no longer observed in the updated analysis. The revised Figure S3 therefore reflects the improved dataset. We agree that glacier variability exhibits non-linear temporal behaviour, and we will clarify in the manuscript that the reported linear trends are intended as first-order summaries of long-term tendencies rather than complete descriptions of temporal evolution.
Comment 39: Supplementary figures: Please improve the readability of all figures. Several labels are too small to read at normal viewing size.
Response: We will revise all supplementary figures to improve readability by increasing font sizes, axis labels, tick labels, legends, and annotation text where appropriate. We will also adjust panel spacing and figure dimensions to improve clarity at normal viewing size.
Citation: https://doi.org/10.5194/egusphere-2026-2386-AC1
-
AC1: 'Reply on RC1', Ignacio Fuentes, 11 Aug 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 -
AC2: 'Reply on RC2', Ignacio Fuentes, 16 Aug 2026
Comment 1: 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.
Response: We thank the reviewer for the careful and detailed assessment of our manuscript and for recognising the novelty of the study, the remote-sensing workflow, and the glacier-response framework. We appreciate the reviewer's thoughtful suggestions regarding the description of the remote-sensing processing, hydrometric data, statistical analyses, and interpretation of the results. In the revised manuscript, we will substantially expand the methodological descriptions, incorporate an independent Sentinel-2 validation of the snow-classification procedure, clarify the statistical workflow and terminology, improve figure readability, and revise the discussion to distinguish more carefully between statistical association and physical interpretation. Detailed responses to each of the reviewer's specific comments are provided below.
My comments are summarized below in Major and Minor Comments:
Comment 2: 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.
Response: Our analysis uses the most recent national glacier inventory to define a consistent analysis domain across the complete Landsat archive. Consequently, the results describe temporal changes within the glacier area that remained present in the reference inventory, rather than reconstructing historical changes in glacier extent. We agree that historical glacier outlines would provide a more complete representation of long-term glacier evolution. However, consistent time-varying glacier inventories are not available for the study region over the full 1984–2026 period. Actually, the only previous version is from 2014. Conversely, adopting the earliest glacier outlines throughout the analysis would progressively include areas that subsequently became ice-free terrain, introducing biases into glacier albedo, LST and snow-cover estimates because these metrics would increasingly represent exposed bedrock rather than glacier surfaces. We will therefore retain a fixed reference inventory to ensure consistent spatial sampling through time and will add a discussion of this limitation to the manuscript.
Comment 3: 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.
Response: The objective of our classification is to distinguish snow-covered from non-snow surfaces rather than to perform a complete land-cover classification. Consequently, bare rock, glacier ice, debris-covered ice, water bodies and vegetation are all treated as non-snow surfaces. We agree that optically bright water bodies may occasionally be classified as snow, particularly near glacier margins. To evaluate the practical impact of these commission errors, we substantially expanded the manuscript by including an independent validation using Sentinel-2 imagery (21 validation cases; 3,599 manually interpreted reference points). The validation was intentionally performed within the bounding rectangle of a 500 m glacier buffer, thereby including surrounding rock, debris, terrain shadows and occasional water bodies. This represents a more demanding test than validation restricted to glacier surfaces alone. The results demonstrated high classification performance (overall accuracy = 0.98, F1 = 0.98, κ = 0.96), while the remaining commission errors were found to occur primarily along snow–rock transition zones and isolated water bodies. We have clarified these points in both the Methods and Results sections and will include this as a limitation suggesting the inclusion of additional water indices or the use of ANDSI for improving snow segmentation.
Comment 4: 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?).
Response: We thank the reviewer for this comment. We agree that demonstrating the temporal consistency of the Landsat record is essential for interpreting long-term glacier trends. In the revised manuscript we have substantially expanded this assessment. Specifically, we now present a quantitative cross-sensor consistency analysis comparing nearly contemporaneous Landsat 5, 7, 8 and 9 observations for all glacier proxies. This analysis demonstrates high consistency for snow-covered fraction, glacier albedo and snow-covered elevation, supporting the use of the Collection-2 Level-2 products without additional empirical cross-sensor harmonisation. For LST, the analysis identified a systematic warm bias affecting Landsat 7 observations during its final years of operation. Consequently, Landsat 7 thermal observations acquired after 2020 were excluded from the LST analyses, while all other glacier proxies retained the full Landsat record.
We also clarify that glacier LST represents the median Landsat-derived radiometric surface temperature across all valid pixels within each glacier polygon rather than the temperature of snow or glacier ice alone. Accordingly, increasing LST reflects the combined effects of changing snow cover, increasing exposure of debris-covered ice, and changes in surface energy balance, and should not be interpreted as warming of melting glacier ice.
Finally, although the number of available Landsat observations increased over time, all glacier proxy series were aggregated to monthly resolution prior to trend, variability and persistence analyses. Consequently, the analysed time series have a uniform temporal sampling interval, substantially reducing the influence of changes in acquisition frequency on the derived metrics.
Comment 5: 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.
Response: We agree that the physical-consistency analysis is an internal diagnostic and does not constitute an independent validation of the snow classification. We have therefore added a separate validation based on near-contemporaneous Sentinel-2 imagery, as described in the new Classification validation subsection. The classification developed in this study is intentionally binary, distinguishing snow-covered from non-snow-covered surfaces. Accordingly, the independent validation was formulated consistently with the classification objective, using manually interpreted snow/no-snow reference points across 21 validation cases and 3,599 reference observations. The validation domain was deliberately extended beyond glacier polygons to include surrounding rock, debris, terrain shadows, and occasional water bodies, thereby testing the method under spectrally challenging conditions. We will revise Section 2.5 to make clear that the physical-consistency analysis is complementary to, rather than a substitute for, the independent classification validation.
Comment 6: 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.
Response: We thank the reviewer for this important observation. We agree that the conventional regression -test does not fully account for the fact that the response variable is itself a derived statistic obtained by selecting the lag corresponding to the maximum absolute cross-correlation over multiple candidate lags, nor for the dependence introduced by nested catchments and shared glacier information. Consequently, the nominal regression -value may overstate the evidence for an association.
To address this concern, we will add a surrogate-based null-model analysis that explicitly reproduces the principal analytical steps leading to the observed regression. Surrogate glacier-proxy time series will be generated under a null hypothesis of no glacier–streamflow coupling while preserving their temporal dependence. Each surrogate realisation will then be propagated through the same preprocessing, prewhitening, cross-correlation and lag-selection procedure used for the observations, followed by the corresponding PCA/regression analysis. Repeating this procedure across a large number of surrogate realisations will provide an empirical null distribution of model , against which the of the revised observed model will be evaluated.
This analysis will allow us to quantify directly how frequently a model of comparable explanatory power could arise from the lag-selection procedure under the null hypothesis. We will report the empirical probability of obtaining an at least as large as the observed value and revise the interpretation of the glacier–streamflow relationship accordingly. We will also retain the catchment-overlap diagnostics and explicitly acknowledge the non-independence associated with nested catchments.
Comment 7: 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.
Response: We thank the reviewer for these observations. We agree that the phrase "physically based glacier proxies" was overly strong and have replaced it throughout the manuscript with terminology such as "complementary glacier-surface proxies" or " remotely sensed glacier-surface proxies" as appropriate. We also agree that the median elevation of snow-covered pixels should not be interpreted as transient snowline altitude. In the revised manuscript, we have clarified the definition of this metric throughout the Methods and Results and expanded the Discussion to explain why trends in median snow-covered elevation are expected to differ from published changes in transient snowline altitude. Specifically, we emphasise that the metric integrates snow extent and glacier hypsometry and therefore provides a complementary descriptor of snow-cover distribution rather than a direct estimate of snowline migration.
Comment 8: 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.
Response: We appreciate this suggestion and agree that glacier-area weighting is more appropriate and transparent for the regional aggregation of the glacier proxies analysed here, particularly for snow-covered fraction. Accordingly, we will replace the previous water-equivalent weighting by glacier-area weighting throughout the analysis. Glacier area was obtained directly from the DGA Public Glacier Inventory (IPG2022_v2) and corresponds to the mapped glacier extent reported in the inventory. We will revise the Methods accordingly and soften the interpretation of regional aggregation in the Discussion.
Comment 9. 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.
Response: We thank the reviewer for this thoughtful comment. We agree that the vegetation analysis produced weaker and less consistent relationships than the streamflow analysis, and we have revised the manuscript accordingly by reducing its prominence in the Abstract, Results and Discussion. We nevertheless decided to retain the analysis because negative results provide useful information regarding the limited explanatory power of glacier variability alone for basin-scale vegetation dynamics. We have also expanded the Methods to justify the choice of the PKU GIMMS NDVI dataset. The primary motivation was to obtain a temporally consistent vegetation record spanning the complete Landsat analysis period (1980s–2022), thereby enabling direct comparison with the glacier proxy time series over multiple decades. Although higher-spatial-resolution products such as Landsat or MODIS are available, the objective of this study was to characterise catchment-scale vegetation dynamics rather than local vegetation patterns. At this spatial scale, the additional detail provided by Landsat imagery is largely averaged during catchment aggregation while substantially increasing data volume and computational requirements. We now explicitly acknowledge that the coarse spatial resolution of GIMMS may contribute to the weaker vegetation relationships observed, particularly in narrow Andean catchments, and identify higher-resolution vegetation products as an important avenue for future work focused on local glacier–vegetation interactions.
Comment 10. 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?
Response: We thank the reviewer for this suggestion. We have substantially expanded the description of the hydrometric dataset. The revised manuscript now reports the 17 gauging stations used in the analysis and provides their basin areas, glacierised fractions, number of contributing glaciers, mean glacier elevation, streamflow record length and completeness, and nesting relationships in a new Supplementary Table. Catchments range from approximately 110 to 14,915 km² and include both small glacier-influenced headwaters and large downstream basins. We also now explicitly state that several stations are hydrologically nested and therefore not statistically independent, because downstream discharge integrates upstream flows.
Comment 11: 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.
Response: We thank the reviewer for highlighting this issue. We have revised the manuscript and associated repository to improve the reproducibility and accessibility of the study. The complete workflow used to derive the glacier-scale time series from Google Earth Engine has now been made publicly available through the Github repository, together with the derived glacier-scale time series, catchment-level aggregations, and the station/catchment information used in the analyses. The associated derived datasets have been archived in Zenodo [10.5281/zenodo.21908291], providing a persistent and versioned record corresponding to the revised manuscript. A development version of the code is also maintained on GitHub.
We have additionally added a dedicated Data availability section to the manuscript describing the sources of the original remote-sensing and ancillary datasets and providing access to the derived datasets and code required to reproduce the analyses. The revised text also distinguishes between third-party datasets, which remain available from their original providers, and the derived products generated in this study.
Comment12. Minor comments. Abstract. L5–6: "physically based glacier proxies" should be “remotely sensed observations” or “remotely sensed proxies”
Response: Thanks for the suggestion. We agree and have replaced “physically based glacier proxies” in the Abstract with “remotely sensed glacier-surface proxies” consistent with the terminology adopted throughout the revised manuscript.
Comment 13: L11–13: Scale dependence is likely an artefact of the w.e. weighting choice (see major comments).
Response: We thank the reviewer for raising this concern. We agree that weighting regional glacier proxies by estimated water equivalent could disproportionately emphasise glaciers with greater estimated ice volume and thereby influence the apparent contrast between regional and glacier-level behaviour. We will therefore revise the regional aggregation throughout the manuscript to use glacier surface area from the Chilean Public Glacier Inventory as the weighting factor. This provides a more direct spatial representation of the remotely sensed surface proxies and avoids introducing estimated ice volume into their aggregation. After this change, differences between regional area-weighted and individual-glacier variability trends remain for some proxies, indicating that the observed scale dependence is not solely an artefact of the previous water-equivalent weighting. We will nevertheless revise the Abstract and corresponding Results and Discussion to describe this pattern more cautiously as scale-dependent behaviour rather than attributing it generally to the influence of large glaciers.
Comment 14: L15–16: "R² ∼ 0.72" — state n = 15 here. As written, this reads as a large-sample result.
Response: We agree and thank the reviewer for this suggestion. We will report the sample size (n=15) alongside the model R² in the Abstract to make the limited sample size explicit and avoid implying that this estimate is based on a larger sample. The regression analysis will also be revised following the methodological changes introduced during revision, and the reported R² and associated statistics will be updated accordingly.
Comment 15: L16: Consider whether the vegetation null belongs in the abstract at all, also see major comment regarding the vegetation analysis.
Response: We thank the reviewer for this suggestion. We agree that the relevance of the vegetation result in the Abstract depends on the scope and robustness of the revised vegetation analysis. In response to the reviewer’s major comment, we will reassess this component of the analysis and revise its interpretation accordingly. We will then reconsider whether the vegetation result warrants inclusion in the Abstract; if retained, its wording will be revised to accurately reflect the scope and strength of the evidence.
Comment 16: 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.
Response: We agree with the reviewer. The original references are not the most appropriate anchors for this broad statement regarding glacier sensitivity to climatic variability and long-term warming. We have therefore replaced them with more widely established assessments and observational evidence, including Hock et al. (2019) and Zemp et al. (2019).
Comment 17: 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.
Response: We thank the reviewer for this suggestion. We agree that this information is important for evaluating the physical context and plausibility of the glacier–streamflow relationships. We have expanded the Study Area section to describe the glacier sample and catchment glacierisation more explicitly. From the Chilean Public Glacier Inventory, our analysis retained 177 glacier units classified as valley or mountain glaciers, excluding rock glaciers and other glacier classes not considered in the analysis. These glaciers cover a total area of 346.3 km2 with individual areas ranging from 0.26 to 24.47 km2 (median = 0.78 km2) and elevations spanning 2602–6049 m a.s.l. We also now report that glacier cover across the 17 gauged catchments ranges from 0.13% to 6.27%. Catchment-specific glacier cover, number of glaciers, mean glacier elevation, basin area, and streamflow record characteristics are provided in Table S1.
Comment 18: 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.
Response: We thank the reviewer for this suggestion. Figure 1 has been revised to provide additional geographic context. The Maipo and Aconcagua basins and the Valparaíso and Metropolitana administrative regions are now explicitly identified in the map, and major urban centres have been added as geographic reference points. The legend has also been revised to more clearly distinguish administrative boundaries, selected catchments, and the glacier subset used in the analysis.
Comment 19: 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.
Response: We thank the reviewer for identifying this ambiguity. We have revised the Methods to explicitly state that the analysis retained DGA glacier units classified as valley or mountain glaciers, while rock glaciers, glaciarets, and other glacier classes were excluded. No minimum glacier-area threshold was imposed; the smallest retained glacier unit has an inventory area of 0.26 km$^2$. The resulting dataset comprises 177 glacier units. Four glacier identifiers contain separate A and B inventory units, which were retained independently because they correspond to distinct polygons in the DGA inventory. We have also revised Figure 5 by replacing the previous “Rocky and small glaciers” category with “Other glacier classes (excluded),” clarifying that these polygons were excluded according to glacier classification rather than glacier size. We note that the 5-ha criterion described for the physical-consistency assessment refers only to the minimum valid observed area within individual monthly observations and is not a glacier-size filter.
Comment 20: L108–111: State the Landsat collection and tier explicitly (Collection 2 Level-2? Tier 1 only?).
Response: We thank the reviewer for noting this omission. We have revised the Methods to explicitly identify the Landsat products used in the analysis as Collection 2, Tier 1, Level-2 products from Landsat 5 TM, Landsat 7 ETM+, and Landsat 8--9 OLI/TIRS. The revised text also clarifies that the analysis used the corresponding surface reflectance and surface temperature products.
Comment 21: 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.
Response: We thank the reviewer for identifying this inconsistency. The reference to 1984--2026 was a typographical error. The Landsat dataset extends from 14 October 1984 to 26 December 2025, as reported in Table 1. We have corrected the text and standardised the study period as 1984-2025 throughout the manuscript.
Comment 22: 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.
Response: We thank the reviewer for highlighting this limitation. We agree that NASADEM represents topography near the 2000 SRTM acquisition epoch and therefore does not capture glacier-surface elevation changes occurring throughout the 1984--2025 study period. We have now explicitly acknowledged this limitation in the Methods and Discussion. In particular, glacier thinning may introduce a non-zero bias in snow-covered elevation estimates and, to a lesser extent, in terrain-based shadow and illumination corrections. We therefore clarify that temporal changes in snow-covered elevation are derived relative to a static reference topography and should be interpreted accordingly.
Comment 23: 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).
Response: We thank the reviewer for identifying that the number of gauging stations was not explicitly stated in the text. We have revised the Methods to clarify that 17 gauging stations met the selection criteria and were retained for the streamflow analysis. Station-specific information is provided in Table S1, including station ID, catchment area, glacier cover, number of glacier units, mean glacier elevation, streamflow record and completeness, and whether the catchment is nested within another analysed catchment. Regarding select lag, r at lag, these will be added to the results section.
Comment 24: L120: "at least 240 overlapping monthly observations" — contiguous, or with gaps?
Response: We thank the reviewer for requesting this clarification. The 240 overlapping monthly observations were not required to be contiguous; isolated gaps in either the streamflow or glacier time series were permitted. We have revised the Methods to make this explicit. In the revised version, we will also correct the station-selection implementation so that the temporal overlap criterion will be based on the number of valid paired monthly observations rather than the total calendar span between the first and last overlapping observations. Missing months are retained explicitly in the monthly time axis and are not treated as consecutive observations.
Comment 25: 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?
Response: In the revised version we will add a sensitivity test to better communicate the glacier coverage filter decision as an offset between data availability and results consistency.
Comment 26: Snow quantification. L146: Can you comment on if the NDSI is good at differentiating snow vs ice vs water vs rock.
Response: We thank the reviewer for raising this point. NDSI is effective at distinguishing snow from many snow-free surfaces because snow exhibits high visible reflectance and strong absorption in the SWIR, but it does not provide perfect discrimination among all surface classes. Our newly added validation analysis confirms this limitation. In particular, we found that water was the principal non-snow class occasionally confused with snow, whereas discrimination from exposed rock was generally more robust. We have therefore clarified that the method should not be interpreted as a general-purpose classifier capable of uniquely separating snow, ice, water, and rock.
Importantly, the validation was deliberately conducted over a broader area than that used in the glacier analysis. Validation samples were obtained within the rectangular bounds of a 500-m buffer around each glacier, thereby including surrounding water bodies, exposed rock, and other non-glacier surfaces. In contrast, snow-covered fraction was calculated only within the glacier boundaries. Consequently, the occurrence of water pixels in the validation dataset is substantially greater than expected within the operational classification domain, and the observed snow–water confusion likely represents a relatively limited source of uncertainty for the glacier-scale snow fractions. Nevertheless, because supraglacial water or water included within the adopted glacier masks cannot be entirely excluded, we now acknowledge this limitation explicitly in the revised manuscript in the discussion and modified the NDSI introduction in methodology to clarify this: “NDSI exploits the high reflectance of snow in the visible spectrum and its strong absorption in the SWIR and is widely used for mapping snow-covered surfaces in mountain environments (Salomonson and Appel, 2004). However, NDSI does not uniquely discriminate among all surface types, with potential spectral confusion between snow and other surfaces, particularly water and glacier ice under some conditions (Mohammadi et al., 2023)”.
Comment 27: 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.
Response: We thank the reviewer for identifying this ambiguity. The VIS term was not min–max rescaled independently for each scene or glacier. VIS was calculated as the mean of the topographically corrected blue, green, and red reflectances and was subsequently linearly transformed using fixed reflectance bounds of 0.05 and 0.60, corresponding to scaled values of 0 and 1, respectively. Thus, (VIS_{scaled}=(VIS-0.05)/(0.60-0.05)), with the same transformation applied to all Landsat observations and glaciers. Consequently, a given VIS reflectance corresponds to the same scaled value throughout the time series, and the VIS component of the SnowScore is not scene-relative. We agree that the original wording ("rescaled to the 0–1 range") did not adequately describe this procedure and could be interpreted as scene-wise or glacier-wise min–max normalisation. We have therefore revised the Methods to explicitly state the scaling procedure and the fixed bounds used.
Comment 28: 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.
Response: We thank the reviewer for raising this point. We agree that the original manuscript did not provide sufficient justification for the 0.7/0.3 weighting coefficients, percentile limits, and fallback thresholds used in the SnowScore and adaptive-threshold procedure. These values were initially selected empirically during algorithm development, with greater weight assigned to NDSI because of its established sensitivity to snow, while retaining a smaller contribution from visible reflectance to incorporate the high reflectance of snow in the visible wavelengths. However, they were not formally optimised against an independent validation dataset in the original analysis.
In the revised manuscript, we will therefore use the newly developed Sentinel-2 validation dataset (21 Landsat–Sentinel-2 image pairs and 3,599 manually interpreted reference points) to evaluate the sensitivity of classification performance to reasonable variations in these parameters. This analysis will allow us to determine whether the selected parameterisation is supported by independent validation and whether classification performance is robust to alternative weighting and threshold choices. The final parameter values and corresponding sensitivity results will be reported explicitly in the revised Methods and Supplementary Materials.
Comment 29: 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.
Response: We thank the reviewer for requesting this clarification. Otsu thresholding was performed independently for each glacier and observation date and was not based on histograms pooled across dates, glaciers, or complete Landsat scenes. Specifically, for each glacier-date observation, the (snowScore) histogram was constructed using valid pixels within a 500-m buffer surrounding the glacier. The buffered domain was used both to increase the number of pixels available for threshold estimation and to improve representation of snow and non-snow surfaces in the distribution used for segmentation. This is particularly relevant for small or strongly snow-covered glaciers, for which restricting the histogram to the glacier polygon could result in very few pixels or a distribution strongly dominated by a single surface class. The resulting threshold was then applied within the glacier boundary to derive the snow-related glacier proxies.
Otsu was implemented as part of an adaptive thresholding procedure. When the fraction of pixels with (snowScore>0.5) was between 0.10 and 0.90, the Otsu threshold was calculated from the corresponding buffered glacier-date histogram and constrained between its 10th and 90th percentiles. When this fraction exceeded 0.90 or was below 0.10, fixed thresholds of 0.45 and 0.30, respectively, were used instead. Thus, the procedure does not rely on Otsu segmentation when the threshold-estimation domain is strongly dominated by one class.
We agree that reporting the number of pixels used for threshold estimation would provide an additional useful diagnostic. Although this quantity was not retained in the original processing output, the 500-m buffered domain was specifically used to avoid deriving thresholds from the limited number of pixels contained within the smallest glacier polygons. We have revised the Methods to explicitly state the histogram domain and the behaviour of the adaptive thresholding procedure for class-dominated observations.
Comment 30: 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.
Response: We thank the reviewer for identifying the ambiguity in the fallback procedure and for raising the possibility that temporal changes in fallback use could influence the inferred snow-covered fraction (SCF) trend. We have clarified the decision rule in the revised Methods. For each glacier-date observation, the fraction of pixels within the threshold-estimation domain with (snowScore>0.5) was evaluated. When this fraction exceeded 0.90, indicating a predominantly high-snowScore distribution, a fixed threshold of 0.45 was applied. When the fraction was below 0.10, indicating a predominantly low-snowScore distribution, a fixed threshold of 0.30 was applied. For intermediate cases, the percentile-constrained Otsu threshold was used.
We additionally quantified fallback occurrence through time. Across the complete dataset, 48.5% of glacier-date observations used a fallback threshold, comprising 32.5% using the high-snowScore threshold 0.45 and 16.0% using the low-snowScore threshold (0.30). Overall fallback frequency did not exhibit a significant temporal trend (−0.240 percentage points yr-1, R2=0.054, p=0.139). However, the two fallback conditions showed opposing temporal changes: use of the 0.30 threshold increased by 0.545 percentage points yr-1 (R2=0.570, p<0.001), whereas use of the 0.45 threshold decreased by 0.786 percentage points yr-1 (R2=0.369, p<0.001). Because fallback selection is determined directly from the observed snowScore distribution, these opposing changes indicate a shift from predominantly high- toward predominantly low-snowScore conditions through time.
To test the reviewer's concern that this changing fallback composition could itself generate the inferred SCF trend, we repeated the regional trend analysis after excluding all observations for which either fallback threshold was used and retaining only classifications based on Otsu-derived thresholds. The negative SCF trend persisted in this independent subset (−0.00252 yr-1; 95% CI: −0.00385 to −0.00118 yr-1, HAC p=0.00022). Thus, the declining SCF signal is also present when fixed fallback thresholds are completely excluded, indicating that the long-term decline is not solely generated by temporal changes in fallback use.
We have clarified the fallback decision rule in the Methods and added the fallback-frequency and Otsu-only sensitivity analyses to the revised manuscript/Supplement.
Comment 31: 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.
Response: We agree that residual illumination effects in cast-shadow areas may influence the VIS component of the composite snow metric even though reflectance is topographically corrected before calculation of NDSI, VIS, and snowScore. Our original motivation for retaining shadowed pixels was to avoid systematically excluding recurrently shaded portions of steep glacier terrain, while albedo was restricted to illuminated pixels because of its direct dependence on reflectance magnitude. We have clarified this rationale in the Methods. To assess whether retention of shadowed pixels influenced the inferred SCF trends, in the revised version we will repeat the SCF analysis after excluding DEM-derived shadowed pixels from both the snow-covered and valid glacier areas while retaining the original classification thresholds.
Comment 32: 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?
Response: We thank the reviewer for requesting this clarification. Broadband albedo was estimated using the narrowband-to-broadband conversion of \citet{liang2001narrowband}, applied to the corresponding spectral bands of each Landsat sensor. The applicability of this formulation to OLI observations over snow- and ice-covered surfaces is supported by \citet{traversa2021landsat}, and this has now been clarified in the Methods. We also note that potential discontinuities among Landsat sensors were explicitly evaluated through our cross-sensor comparison using temporally matched observations (now included as Supplementary materials). For glacier albedo, median inter-sensor biases were +0.0001 for Landsat 5–7, −0.0039 for Landsat 7–8, and −0.0002 for Landsat 8–9, providing no evidence of an appreciable systematic step change in the albedo series. Accordingly, no additional empirical sensor harmonisation was applied. Regarding solar zenith angle, no additional solar-zenith or BRDF correction was applied specifically to the broadband albedo calculation. However, the surface reflectance entering the albedo calculation had previously undergone DEM-based topographic correction incorporating scene-specific solar geometry. We have clarified the sensor applicability of the albedo formulation in the revised Methods.
Comment 33: L182–183: LST — see major comments.
Response: We thank the reviewer. This issue is addressed in detail in our response to Major Comment 4. In response to that concern, we conducted a quantitative cross-sensor consistency assessment for all glacier proxies. While SCF, glacier albedo, and snow-covered elevation showed strong agreement among overlapping Landsat sensors, the analysis identified a systematic warm bias in Landsat 7 LST observations during the final years of the mission. Landsat 7 thermal observations acquired after 2020 were therefore excluded from the revised LST analyses.
We have also clarified the definition and interpretation of LST in the Methods. LST represents the median Landsat Collection 2 Level-2 surface temperature across all valid pixels within the glacier outline and is not restricted to snow-covered pixels. It therefore characterises the thermal state of the glacier surface as a whole, including changes associated with snow-cover loss and increasing exposure of debris-covered ice and other non-snow surfaces, rather than the temperature of melting glacier ice alone. The revised Methods now explicitly state this definition and the Landsat 7 exclusion.
Comment 34: L184–186: Median snow elevation is not typically referred to as snowline altitude (see major comments).
Response: We agree and thank the reviewer for highlighting this distinction. As detailed in our response to Major Comment 7, the metric calculated in this study is the median (and mean) elevation of snow-covered pixels and should not be interpreted as transient snowline altitude. We have revised the terminology throughout the manuscript accordingly and now refer to this metric as snow-covered elevation rather than snowline altitude. The Methods have also been revised to clarify that this metric describes the vertical distribution of snow-covered area across the glacier and is influenced jointly by snow extent and glacier hypsometry; it therefore represents a descriptive snow-cover proxy rather than a direct estimate of transient snowline or equilibrium-line altitude.
Comment 35: L185: Mention limitations of SRTM based DEMs on glacier elevation.
Response: We agree. The revised manuscript now explicitly acknowledges this limitation in the Satellite imagery and topography subsection. NASADEM is a reprocessed SRTM-derived DEM representing topography near the 2000 acquisition epoch and therefore does not capture glacier-surface elevation changes occurring over the 1984–2025 study period. We have further clarified that temporal changes in our snow-covered elevation metrics consequently represent changes in the spatial distribution of snow over a fixed reference topography rather than changes in the glacier surface elevation itself.
Comment 36: 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?
Response: The term “water-equivalent-weighted” was removed from the revised manuscript as the reviewer suggested using an aggregation based on glacier surface area, not glacier thickness, volume, or water-equivalent storage. Specifically, each glacier was weighted by its mapped area reported in the Chilean Public Glacier Inventory, and the regional value at each time step was calculated as the weighted mean across glaciers with valid observations at that time. We have now explicitly defined this weighting scheme and included the corresponding equation in the Methods. The same glacier-area weighting was used when aggregating glacier proxy series to catchment scale for the lagged analyses.
Comment 37: L216–218: The 84-month window is asserted without justification. Please explain.
Response: We thank the reviewer for requesting clarification. The 84-month window corresponds to seven years and was selected as a compromise between statistical stability and temporal resolution. Shorter windows provide fewer monthly observations for estimating SD and, particularly, lag-1 autocorrelation, making these metrics more sensitive to individual years and missing observations, whereas substantially longer windows would increasingly smooth multi-year changes and reduce the number of effective temporal positions available over the 1984–2025 record. A seven-year window provides up to 84 monthly observations and spans multiple annual cycles while retaining sufficient temporal resolution to examine changes in variability and persistence. We have added this rationale to the revised Methods.
Comment 38: Lagged analysis. L222–226: Nested catchments share glaciers, i.e. pseudoreplication. Report the number of glaciers per catchment and the overlap matrix.
Response: We thank the reviewer for highlighting this issue. We agree that the nested structure of several analysed catchments results in partial non-independence among the catchment-scale glacier signals because upstream glaciers may also contribute to the glacier ensemble of downstream catchments. We have therefore made this structure explicit in the revised manuscript. Supplementary Table S1 now reports the number of glaciers associated with each of the 17 gauged catchments, together with their glacier cover and other catchment characteristics. In addition, we have added Figure S7 showing the pairwise overlap in glacier membership among the analysed catchments, separately for the Aconcagua and Maipo basins. The figure reports the number of glaciers shared by each pair of catchments and makes the nested structure of the glacier ensembles explicit.
We have also revised Section 2.7 to clarify that lag estimates from nested catchments are not statistically independent and should therefore be interpreted as spatially distributed catchment-level diagnostics rather than as independent replicates. The purpose of the station-wise CCF analysis is to examine how the timing and strength of glacier–streamflow/vegetation associations vary through the drainage network, rather than to treat the 17 catchments as independent experimental units.
Comment 39: 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.
Response: We agree with the reviewer that the original statement was too strong and could not be supported by visual inspection of two representative glaciers alone. We have therefore removed the statement that the classification patterns are “not sensitive to sensor transitions or data gaps” and substantially revised this section.
In the revised manuscript, Figure 4 is explicitly presented only as a qualitative illustration of the behaviour of the segmentation approach under contrasting acquisition conditions and sensors. Classification performance is now evaluated quantitatively using an independent Sentinel-2 validation comprising 3,599 manually interpreted reference points across 21 glacier-date cases. This yielded a pooled overall accuracy of 0.98, snow F_1 score of 0.98, and Cohen's kappa of 0.96, with case-level validation statistics reported in the Supplementary Material. We have also added quantitative diagnostics of the temporal behaviour of the adaptive thresholding procedure, including the frequency of Otsu and fallback threshold use through time and across sensors.
More broadly, cross-sensor consistency of the derived glacier proxies is now evaluated quantitatively using nearly contemporaneous Landsat observations, as described in our response to Comment 4. We have revised the Results accordingly so that the multi-temporal examples are no longer used to support an archive-wide claim of sensor insensitivity.
Comment 40: 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.
Response: We thank the reviewer for highlighting the sensitivity of the original snow-covered-area analysis to the valid-area threshold. We agree that presenting a trend obtained only after imposing a ≥99% valid-area criterion did not adequately demonstrate the robustness of this result.
Rather than retaining this threshold-dependent analysis, we substantially revised the treatment of incomplete spatial coverage. We now evaluate snow-covered fraction (SCF) using both the complete inventory glacier area and the valid observed glacier area as alternative denominators across minimum usable-coverage thresholds of 50%, 70%, 80%, 90%, and 95%. This sensitivity analysis demonstrates that normalisation by valid observed area is substantially less dependent on spatial completeness, while agreement between both formulations progressively increases at higher coverage thresholds. At the same time, increasingly restrictive coverage criteria substantially reduce temporal sampling. Based on this trade-off, we adopted a minimum usable-coverage threshold of 80% and calculate SCF relative to the valid observed glacier area throughout the revised analysis (Figure S3).
As an additional consistency check, we compared total snow-covered glacier area estimated directly from classified Landsat pixels with an independent inventory-based proxy calculated as SCF × inventory glacier area, using the revised quality-controlled dataset. The two estimates were nearly identical (n=1363, Pearson r = 0.999, Spearman rho=1), with a regression of inventory-proxy SCA = pixel-based SCA + 0.556 km². This indicates that, under the revised coverage criterion, estimates of total snow-covered area are not strongly dependent on the choice between direct pixel summation and inventory-area scaling.
The previous statement concerning a snow-covered-area trend obtained using ≥99% valid glacier coverage has therefore been removed, and the revised manuscript no longer relies on this threshold-dependent result.
Comment 41: L315–316: "except for LST at the regional scale" — this exception is passed over. Why does LST behave differently?
Response: We thank the reviewer for highlighting this exception. Following the methodological revisions introduced elsewhere in the manuscript, including the revised glacier-area weighting, updated quality filtering, and correction of the Landsat 7 LST record, the temporal-persistence analysis was recalculated. The previously reported regional LST exception is no longer present in the revised results. We have therefore removed the original statement and updated the Results to reflect the recalculated AC1 trends.
Comment 42: 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.
Response: We agree with the reviewer that the previous statement was too strong. The weak or moderate relationships between glacier characteristics and observed proxy changes do not, by themselves, demonstrate that climatic forcing dominates long-term glacier change. We have therefore removed this inference. In addition, because the PCA and glacier-property analyses were recalculated following the methodological revisions, this section has been rewritten to reflect the revised multivariate modes and their relationships with glacier characteristics. The revised interpretation now distinguishes between coherent regional-scale behaviour and glacier-specific heterogeneity without attributing the unexplained component specifically to climatic forcing. We also clarify that the analysed glacier characteristics account for only part of the observed between-glacier variability and that the remaining variability may reflect environmental forcing, glacier-specific processes, and observational effects not explicitly represented by these predictors.
Comment 43: 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.
Response: We agree. The 17 gauged catchments analysed in this study belong to only two major basin systems (Aconcagua and Maipo), and we therefore agree that basin-level cluster-robust standard errors do not provide a reliable basis for inference. The previous wording could incorrectly imply that this sensitivity analysis provided statistical reassurance despite the very small number of clusters. We have therefore removed the cluster-robust results and the corresponding statement from the revised manuscript. Instead, we now explicitly characterise the dependence arising from the nested catchment structure by reporting the number of glaciers associated with each gauged catchment (Table S1) and a pairwise glacier-membership overlap matrix for the Aconcagua and Maipo catchments (Figure S7). We also acknowledge that lagged relationships from nested catchments are not statistically independent and interpret spatial patterns across gauges accordingly.
Comment 44: 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.
Response: We agree with the reviewer. The previous wording interpreted the regression coefficients too mechanistically and presented the proposed balance between snowpack storage and shorter-term variability as though it had been directly tested. The analysis establishes associations between multivariate modes of glacier-surface variability and downstream response timing, but it does not identify the hydrological mechanisms responsible for those associations. We have therefore revised this interpretation and now explicitly present the proposed mechanism as a hypothesis consistent with the observed relationships rather than as a demonstrated trade-off. We have also softened the interpretation of the non-significant vegetation results, which are now described as an absence of detectable statistical association rather than evidence that NDVI dynamics are not controlled by glacier variability.
Comment 45: Discussion. L396–398: Direct contradiction of L240–241.
Response: We agree with the reviewer. The wording in the original Discussion was inconsistent with the inferential limits explicitly stated in the Methods. In particular, the expressions “exert a strong control” and “actively regulate” implied a causal relationship that cannot be established from the cross-correlation and regression analyses used here. We have therefore removed this causal language throughout the Discussion. The revised text describes the relationships between glacier-surface variability modes and downstream response timing as statistical associations and avoid expressing causal relationships.
Comment 46: 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).
Response: We agree that the original manuscript acknowledged uncertainty without quantifying its magnitude. We have therefore added a first-order propagation of snow-classification uncertainty to SCF estimates using the independent Sentinel-2 validation dataset. Sensitivity and specificity were estimated by case-level bootstrap resampling of the 21 validation cases, thereby preserving dependence among reference points acquired under common glacier-date conditions. These distributions were propagated to 104,656 glacier-date SCF observations. The resulting median absolute difference between observed and classification-corrected SCF was 0.0031, while the median half-width of the propagated 95% uncertainty interval was 0.0049 SCF units (95th percentile = 0.0218). The full uncertainty analysis is now presented in another Supplementary Figure.
We evaluated incomplete spatial coverage separately because it represents a distinct source of uncertainty from snow/no-snow classification and cannot reasonably be treated as independent pixel-level classification error. This sensitivity analysis is presented in Figure S3. We have also revised the Discussion to distinguish these quantified sources of uncertainty from additional uncertainties associated with spatial resolution, glacier delineation, the static DEM, and residual sensor-related effects, which are not included in the propagated classification intervals.
Comment 47: 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.
Response: We agree that the water-management implications in the previous version extended beyond the evidence provided by our analysis. We have substantially shortened this section and removed the paragraph discussing specific adaptation technologies, including snow-retention measures and artificial ice storage. We have also revised the remaining text to frame the management implications more cautiously, focusing on the potential relevance of changes in hydrological timing and buffering for water-resource reliability rather than implying that our analysis directly evaluates specific management responses.
Comment 48: 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.
Response: We agree and have expanded the limitations section to explicitly address these issues. The revised manuscript now acknowledges residual cross-sensor effects and time-varying observation density, despite the quantitative sensitivity analyses introduced in the revision; the use of static glacier outlines and topography; non-independence arising from shared glaciers among nested catchments; statistical selection associated with identifying maximum correlations across multiple lags; and the absence of an explicit precipitation control. We have also clarified the scope of external validation. Snow classification is now independently validated against Sentinel-2 observations, whereas the inferred glacier–hydrological relationships lack independent external validation and are therefore interpreted as statistical associations rather than causal effects.
Comment 49: Back matter. L474–475: Code "will be made available" — must be archived with a DOI now. Missing: A Data Availability section (mandatory for TC).
Response: We agree. As described in our response to Comment 11, the code and associated analysis materials have now been archived in a permanent repository with a DOI, and the corresponding repository information has been added to the revised manuscript. We have also added the mandatory Data Availability section, specifying the sources and accessibility of the datasets used in the study.
Comment 50: 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.
Response: Thank you for identifying this inconsistency. The Author Contributions statement has been corrected to remove the erroneous reference to I.F. providing guidance to himself and to clarify the respective author contributions.
Comment 51: L480: Ironically, in the reference to ChatGPT being used for language improvement, there are typos: “thepurpose” and "fullresponsibility"
Response: This was corrected in the revised version of the manuscript.
Comment 52: 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.
Response: We thank the reviewer for identifying these citation and reference-list inconsistencies. We have carefully reviewed and corrected the full bibliography and corresponding in-text citations throughout the manuscript. This included restoring proper capitalization for journal titles, datasets, geographic names, satellite/platform names, acronyms and technical terms; correcting author-name parsing errors such as “Colin Cameron and Miller (2015)” to “Cameron and Miller (2015)”; and replacing “Jpl, N. (2020)” with the correct institutional author “NASA JPL (2020)” in both the reference list and figure captions. We also checked the remaining entries for similar formatting and capitalization issues and corrected them throughout.
Comment 53: 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.
Response: Thanks for the suggestion. In the revised Figure 1 the typo was corrected and stations were labelled.
Comment 54: 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.
Response: We agree that the original Figure 3, based on two illustrative glaciers, was insufficient as the primary evaluation of the snow-classification procedure. In the revised manuscript, classification performance is now evaluated quantitatively against independently interpreted Sentinel-2 imagery across 21 glacier-date validation cases, comprising 3,599 reference points and spanning contrasting glacier and snow-cover conditions. The new Figure 3 presents four representative validation cases together with their confusion matrices and the pooled validation results. These examples include challenging conditions and an explicit failure case in which surface water adjacent to the glacier is misclassified as snow. Because validation was performed over an area extending beyond the glacier outlines, these cases also capture commission errors associated with non-glacier surfaces that are not included in the subsequent glacier-scale proxy calculations. Detailed performance statistics for all validation cases are provided in the Supplementary Material. The original Figure 3 has been retained as the new Figure 4, but its role has been revised: it is now presented only as an illustration of the temporal behaviour of the classification across different Landsat sensors and acquisition conditions, rather than as evidence of classification accuracy.
Comment 55: Figure 4. Caption typo: "dhashed" should be "dashed".
Response: Thanks for noticing the typo. It was corrected in the revision.
Comment 56: Figure 5. Legend entry "Rocky and small glaciers" is undefined in the methods. Font sizes are too small.
Response: We agree. The previous legend entry “Rocky and small glaciers” was imprecise because glacier size was not used as an exclusion criterion. We have replaced this category with “Other glacier classes (excluded)” and clarified in the Methods that the analysis retained valley and mountain glacier units, while rock glaciers, glaciarets, and other glacier classes were excluded. No minimum glacier-area threshold was applied. Figure 5 has also been redesigned by combining the regional panels and associated histograms, allowing the maps, symbols, legends, and text to be substantially enlarged and improving overall readability.
Comment 57: 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.
Response: We agree. The CCF insets in the original Figures 10 and 11 were too small for effective interpretation and displayed only a subset of stations without an explicit selection criterion. We have therefore removed these insets from the main figures and added supplementary panel-grid figures showing the complete CCFs for all 17 gauged catchments included in the lag analysis, separately for streamflow and vegetation responses. All panels use consistent formatting and display the analysed lag range, approximate significance bounds, and the lag corresponding to the maximum absolute correlation retained in the analysis. The revised main figures now focus on the spatial distribution of the resulting lag metrics.
Comment 58: Table 1. Please add scenes per year (or a time-series plot of observation density in the supplement).
Response: We agree and have added a new supplementary time-series figure showing annual observation density throughout the study period, separated by Landsat sensor. During this assessment, we also identified that the original Table 1 inadvertently omitted two WRS-2 path–row combinations represented in the processed dataset. Table 1 has therefore been corrected to report the complete set of Landsat scenes and path–row combinations used in the analysis. The new supplementary figure explicitly documents the temporal increase in observation density associated with the progressive overlap of Landsat missions.
Comment 59: New table required. A station/catchment table (see comment on L117–121). Could be added to supplement.
Response: We agree. As also noted in our response to Comment 23, we have added a new supplementary station/catchment table (Table S1) for all 17 gauging stations retained in the lagged analysis. The table reports station and catchment characteristics
Citation: https://doi.org/10.5194/egusphere-2026-2386-AC2
- Static glacier outlines
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 200 | 83 | 18 | 301 | 37 | 14 | 13 |
- HTML: 200
- PDF: 83
- XML: 18
- Total: 301
- Supplement: 37
- 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.