The authors are thankful to the reviewer for constructive comments and suggestions. The authors’ responses are in italics.

Review

The paper tries to evaluate snow cover change over 23 years in the Italian Western Alps using MODIS images. The topic, in general, can be a good analytical paper, but not a methodological one. In the results, showing and explanation, there are significant discussions as follows:

Title

I have a question: why did you title your paper “GEE,” though GEE was just used for downloading the data? The platform is not for downloading data, and all the data you used could be downloaded from other resources.

We appreciate the reviewer’s insightful comment. We agree that Google Earth Engine (GEE) is not designed merely as a data downloading platform, but as a cloud-based environment for large-scale geospatial processing and analysis. In our study, GEE was not used only to access MODIS products, but it represented the core computational environment for the preprocessing and analysis of all satellite data.

Specifically, GEE was employed to apply the MOD10A1 cloud mask, to correct the Snow-Covered Area (SCA) for terrain slope effects, and to compute the annual aggregated products such as integral Snow Cover Area (iSCA) and Snow Persistence (SP). Moreover, the Mean Daily Snowed Area (MDSA) was derived directly in GEE as the ratio between iSCA and SP. These steps were crucial for generating the consistent and spatially explicit snow metrics used in the subsequent trend analysis. Only the final estimation of trend slopes (gain) and offsets, together with their statistical significance (t-tests on regression coefficients), was carried out outside the GEE environment, since these statistical functionalities are not currently implemented within the platform. For these reasons, GEE was instrumental in all phases of the work, from data acquisition and filtering to correction, aggregation, and index computation, serving as the analytical backbone of our methodology.  For this reason, we consider it appropriate to include “GEE” in the paper title, as it reflects the methodological and operational relevance of the platform in the study workflow.

Beyond what has been reported, we have decided to change the title of the paper to:

“A Remote Sensing Approach for Measuring Climatic Change Effects on Snow Cover Dynamics”

Introduction

In the introduction, you mention your goal is “to reconstruct the evolution of snow cover, from 2000 to 2023, using remote sensing technologies.” Please look at https://doi.org/10.3390/rs13152945  (which you also cite) as an example; there should be more as well. What makes your study different from theirs? You state that other studies do not use in situ data for validation, but in fact they do, for example, see https://doi.org/10.5194/hess-10-679-2006  and https://doi.org/10.1016/j.rse.2023.113877

Moreover, MODIS pixels are not directly comparable with in situ measurements because of the large difference in spatial resolution. For MODIS validation, researchers usually compare with higher-resolution satellite data such as Landsat.

Moreover, you mention your goal “using remote sensing” without giving a literature review on how RS is used for snow cover. Your literature review does not specify which data each study uses.

We thank the reviewer for their insightful comments regarding the novelty, validation approach, and literature review. We recognize that these aspects require further clarification and refinement in the manuscript.

The reviewer correctly identifies Fugazza et al. (2021) as we cited as a key study in Alpine snow cover change. While that research is pivotal for understanding snow seasonality (Start, End, and Length of Season) across the entire Alpine region over the 2000–2019 period, our study offers a distinct and complementary perspective focused on the cumulative persistence of snow cover at the local/regional scale within the Piedmont and Aosta Valley region, extending the temporal analysis to 2023.

Crucially, our novelty is centered on the derived metric, denoted as nT​, which quantifies the change of a snow covered area in mean annual snow event over the study period. This metric provides a concise and spatially explicit representation of the accumulated change, which differs from analyses focused on phenological metrics (SOS/EOS). The use of nT​ allows for a deeper understanding of how the intensity of snow cover of a mean snow event has evolved annually. Importantly, because this product reflects the evolution of the trend within each individual pixel, it is inherently comparable across the entire study area. Consequently, the interpretation of the nT​ map is not directly affected by the confounding influence of elevation or other environmental drivers when assessing the rate of change within that specific location. We believe that parameterizing this behavior is of paramount importance for regional impact assessment.

We have revised the Introduction to clearly establish the need for such a parameterized approach, and the Discussion section now emphasizes the innovative nature of the nT​ index as the central core of our work.

We appreciate the reference to the cited validation studies and acknowledge that our initial statement regarding the lack of in situ validation in "local studies" was poorly formulated. We will revise the text to reflect the existence of robust validation literature and will incorporate the suggested references into the discussion.

We are well aware of the scale mismatch between the ≈500m MODIS pixel and the point-scale in situ weather station data. However, we deliberately chose to include these in situ ground observations for the following reason:

• Reliable Daily Presence: The meteorological stations provide a highly reliable, ground-truth daily record of snow presence on the ground, which is less susceptible to cloud-obscuration issues compared to optical satellite data.
• Temporal Resolution Constraint: Higher-resolution data (like Landsat or Sentinel-2) have longer revisit times. Given the high probability of cloud cover and the short duration of many snow events, relying solely on these sensors would compromise our primary objective of analyzing daily snow cover persistence over a long-term time series (2000–2023). The in situ data, despite the scale difference, provided a necessary temporal fidelity benchmark for our daily product analysis.

The reviewer is correct that the literature review lacks sufficient detail regarding the specific Remote Sensing (RS) methodologies and data used in the cited works. We will extensively revise the Introduction to provide a more structured overview of how RS technologies are applied to snow cover studies. Specifically, the revised introduction will: specify the sensor and product used in each referenced study; outline the different RS methodologies (e.g., phenological metrics, frequency/duration analysis, change detection) applicable to snow cover research, thereby clearly contextualizing our approach.

Area of Interest

In Fig. 1 you can see some “+” on the image. Please remove them, as outer latitude and longitude locations are enough.

Ok we will correct it

Data

Why do you use MOD10A1, as we have gap-filled MODIS snow data MOD10A1F / MYD10A1F?

We acknowledge the existence of the gap-filled products, MOD10A1F and MYD10A1F, however, our primary computational environment was the Google Earth Engine (GEE) platform, chosen for its cloud-computing capabilities crucial for the large-scale, long-term processing (as detailed in our response to Question 1).

At the time of data acquisition and initial processing, the MODIS Collection 6 (C6) and later C6.1 daily Snow Cover product, MOD10A1, was the only high-frequency MODIS snow product readily available in the GEE data catalog that covered the entire 2000–2023 study period

2.2.2 Ground Data

Please show the location of in situ data. As I understood from the abstract, there is only 7 station data. Do you think this is enough? Moreover, show them on the map so we can see the elevation of the stations. If the stations are mostly in low altitudes, it doesn’t give an accurate evaluation. Besides, 500 m vs in situ data is not a precise evaluation. How do you find which pixel represents the station?

The selection of only seven in situ stations in the initial manuscript was based on strict methodological criteria to ensure the highest reliability and suitability for our long-term study:

• Temporal Coverage: Only stations providing complete and continuous data across the entire 2000–2023 study period were retained.
• Location and Spatial Sampling: Stations were selected to achieve the best possible spatial distribution across the study area (Piedmont and Aosta Valley) and were chosen to ensure they fell entirely within a single MODIS pixel and not near the pixel boundaries, mitigating geometric uncertainty.
• Data Reliability (Manual Stations): Crucially, only non-automatic (manual) stations were chosen. These stations, were favored because they provided a more reliable ground truth, especially during the summer, where automatic sensors often exhibit false positives due to dense vegetation or other land cover features misinterpreted as snow.

The initial seven stations already spanned a significant elevational range, from approximately 200m to 2000m above sea level (a.s.l.), ensuring that validation was conducted across different climatic zones of the study area.

We recognize, as highlighted by multiple reviewers, that the number of stations is limited. We will perform the revision expanding the validation to include a larger number of in situ stations (9 stations) by relaxing the stringent requirement for the entire 2000–2023 period.

The reviewer is correct that the spatial distribution and elevation of the stations must be visually presented and justified. The revised manuscript will include a new figure showing the location and elevation of all in situ stations. We already produced it but we removed it due to paper length limitations

Data Processing

Linear interpolation for filling the gap for a variable like snow is not suitable, as snow changes rapidly and also varies from one place to another. Use MOD10A1F / MYD10A1F, which is already gap-filled.

We thank the reviewer for this relevant comment. We are aware of the existence of the MOD10A1F and MYD10A1F gap-filled MODIS products. However, these datasets are not currently available within the GEE public catalogue for the entire 2000–2023 period and for the spatial coverage of our study area (the Western Italian Alps). For this reason, their use would require external preprocessing pipelines and large data transfers, which would compromise the efficiency, reproducibility, and consistency of our GEE-based workflow.

Regarding the interpolation method, we recognize that snow cover can change rapidly in time and space. Nevertheless, we adopted a linear interpolation approach because it provides a robust, reproducible, and well-controlled solution to fill data gaps caused mainly by cloud masking. Linear interpolation uses only the immediately preceding and following valid observations to estimate the missing values, avoiding the introduction of artificial patterns or smoothing artifacts. In contrast, higher-order or non-linear approaches such as LOESS, splines, or Savitzky–Golay filters often generate unrealistic oscillations, overshooting, or temporal distortions, especially when applied to large, heterogeneous datasets like MODIS time series over complex mountainous terrain.

The use of more complex interpolation methods would require a dedicated evaluation to verify whether they actually outperform the linear approach under the specific conditions of our dataset, which is beyond the scope of the present study. 

Within our objectives, linear interpolation offers the most balanced trade-off between accuracy, simplicity, and reproducibility, ensuring that the resulting snow cover metrics (iSCA, SP, and MDSA) remain physically consistent and comparable across the entire study domain.

For these reasons, we consider linear interpolation the most appropriate and robust method for the purpose and scale of this study.

Snow Cover Trends Quantification

In general, you cannot consider a trend by evaluating some years and comparing them to one year (2000 in your case). As one year can be dry, wet, etc., you can do a Cumulative Frequency Distribution for example.

We thank the reviewer for this constructive comment. We fully agree that directly comparing individual years to a single reference year (e.g., 2000) could be misleading, as any given year may be characterized by specific climatic anomalies (e.g., being unusually dry or wet). However, this is not the rationale behind our normalization procedure.

In our analysis, the intercept normalization is applied as a methodological step to remove the influence of local baseline conditions rather than to compare actual annual values. Different sites or regions can present markedly different environmental or climatic starting points — for instance, areas that are inherently drier, wetter, or subject to different baseline vegetation or soil moisture levels. Such local effects can bias the magnitude of the estimated gain (trend), making it difficult to compare temporal patterns across sites.

By normalizing the series with respect to the model intercept (which corresponds to the modeled level at "t=0" , i.e., the year 2000 in our parametrization), we express the gain as a relative change rather than as an absolute value. This operation effectively removes site-specific baseline effects, allowing a consistent and dimensionless comparison of trends across spatially heterogeneous locations.

Therefore, the use of the intercept in the normalization does not imply that the year 2000 is taken as a physical or climatic reference year. It is simply the mathematical origin of the time variable in our model, and serves as a standardized reference point for comparing the relative magnitude of temporal changes across sites.

Regarding the suggestion to use a Cumulative Frequency Distribution (CFD), we believe that this approach would not be appropriate in our specific context. A CFD provides information about the statistical distribution of observed values, typically aggregated over time or across samples, and is not designed to capture or represent temporal trends at the pixel level. In our case, the analysis is explicitly spatiotemporal, and each pixel is modeled independently in time to estimate its individual rate of change. Therefore, computing and mapping a CFD for each pixel would not yield a meaningful representation of the underlying temporal dynamics, as the cumulative distribution does not convey directional or temporal information. For our purposes, the linear modeling approach combined with intercept normalization provides a more direct and interpretable framework for quantifying and comparing temporal trends across heterogeneous spatial domains.

In Fig. 2, what does “m a.s.l” stand for?

The abbreviation 'm a.s.l.' stands for meters above sea level, a common geographic unit used to indicate elevation or altitude. We will ensure this is clearly defined either in its first appearance. In this specific case it appears in the legend of the scatter plot (Figure 2). It was intentionally used to represent the elevation of the data points (different shades of blue). The specific tones were chosen to illustrate the variability and correlation of the data as a function of their elevation. This visualization is critical for addressing the previously raised concerns regarding the uncertainties between the point-scale weather station data and the coarse-resolution MODIS pixels. 

In Fig. 4, I also have a question: in high elevations, climate change affects the duration and amount of snow. However, in 23 years, climate change will mostly affect low or mid-elevation. In high elevation, as we do not consider SWE and our goal is just snow cover, how are all the nT values negative, even in high elevation? Because even with climate change, we still have below-zero temperatures to have snow cover, even with lower SWE, for example.

The same question also applies to Fig. 5.

We thank the reviewer for this insightful observation concerning the trends in high-elevation areas and the physical mechanisms driving snow cover changes, particularly over a 23-year period.

We first clarify that the values shown in Figures 4 and 5 represent the spatial and altitudinal  variability of nT​, only for those areas where the calculated trend was determined to be statistically significant. The reviewer is correct that the most dramatic effects of climate change are often concentrated at low and mid-elevations. However, our results demonstrate that, even at high elevations, a discernible negative trend in snow cover is evident, leading to the overall negative nT​ values observed in these areas (albeit the magnitude is much smaller compared to lower altitudes).

This observed snow loss at high altitudes is primarily attributable to two mechanisms:

• Increasing air temperatures are frequently driving the zero-degree isotherm to elevations well above the highest peaks in the study area (≈4500m a.s.l.). This results in non-snow precipitation and extreme melting of the snow during these events.
• The significant reduction in the extent of perennial snowfields within the study region is also captured by the MODIS snow cover product. Areas that historically maintained permanent snow (e.g., above 3500m a.s.l.) have experienced substantial ablation and retreat over the 23-year period. This reduction in historically persistent snow contributes to the overall negative nT​ values, even if the effect is small over a two-decade trend analysis.

Therefore, while temperatures remain below zero for much of the year at high altitudes, the increasing frequency of warm-air intrusions and the substantial reduction of permanent snow explain the overall negative, yet often minimal, trends observed in Figures 4 and 5.

Compare with Landsat/Sentinel snow cover data.

We thank the reviewer for the suggestion to compare our results with snow cover data derived from Landsat or Sentinel missions. While we recognize the scientific value of these higher-resolution datasets, such a comparison is not feasible within the framework and objectives of our study, primarily due to the substantial differences in temporal resolution and data continuity.

MODIS provides daily global snow cover observations which are fundamental for capturing the full temporal evolution of the snow season and for deriving time-integrated metrics. In contrast, Landsat (16-day revisit) and Sentinel-2 (5–10-day revisit) have much coarser temporal sampling, further limited by frequent cloud cover in mountain regions. As a result, these datasets cannot provide a continuous and temporally complete record suitable for computing daily or cumulative indicators such as the number of snow days or long-term persistence trends.

For this reason, MODIS represents the only consistent long-term dataset that allows for the reconstruction of daily snow dynamics over a 23-year period across the entire Western Italian Alps. To compensate for the absence of comparable high-temporal-resolution satellite products, we performed a validation of the snow persistence (SP) and snow-cover duration metrics using ground-based snow depth data from seven snowmeter stations distributed across the study area. This ground validation provides a robust and independent assessment of the MODIS-derived snow indicators, effectively bridging the validation gap that would otherwise remain if only coarser-temporal-resolution satellite datasets (e.g., Landsat or Sentinel) were used. Therefore, the combination of MODIS time series and ground observations ensures both temporal completeness and empirical reliability of our snow cover trend analysis.

Conclusion

The conclusion is too short. Discussion and conclusion can be combined.

To provide a more complete and exhaustive explanation of our findings and their implications, we have revised the manuscript by merging the original Discussion and Conclusion sections into a single, comprehensive section titled: “Discussion and Conclusions”.

All the changes will be update in the revised manuscript

The authors are thankful to the reviewer for constructive comments and suggestions. The authors’ responses are in italics.

The manuscript presents an analysis of snow cover changes in the Western Italian Alps (2000–2023) using MODIS data processed through Google Earth Engine (GEE). The topic is relevant, and the effort to quantify long-term snow dynamics at regional scale is commendable. However, the paper requires substantial revision to reach the scientific depth.

Conceptual framework and novelty



Although the introduction provides context on snow monitoring, it lacks a clear articulation of the research gap and the specific contribution of this study. Similar MODIS-based analyses have been published for Alpine regions, yet the manuscript does not sufficiently differentiate its approach or findings from prior work. The use of the nT index could represent a methodological contribution, but its conceptual basis and added value over standard trend analyses (e.g., Mann–Kendall or Sen’s slope) are not adequately explained.

We thank the reviewer for their valuable feedback regarding the clarity of the Introduction, the articulation of the research gap, and the justification for the use of the nT​ index. We agree that these sections require significant enhancement to clearly establish the novelty and specific contribution of our study.

We acknowledge that the initial manuscript did not sufficiently differentiate our work from prior studies in the Alpine region. While similar long-term analyses exist, our study offers a distinct regional and thematic focus that addresses a specific gap. Our goal is not to track snow phenology but to quantify the accumulated loss of intensity of snow cover for a single mean snow event.

The reviewer correctly identifies the nT​ index as a potential methodological contribution. We recognize that its conceptual basis and added value over standard trend analyses were poorly explained. The nT​ index is not intended to replace statistical trend methods (such as Mann–Kendall or Sen’s slope), but rather to offer a complementary, easily interpretable, and spatially explicit measure.

We have performed substantial revisions to the manuscript to address these points.

2. Methods and reproducibility



The methodological description is too general. The computation of iSCA, SP, MDSA, and nT indices should be fully explained, including equations, temporal aggregation steps (in data gap-filling), and cloud-masking criteria. The validation using seven snow stations might be insufficiently developed: statistical metrics (correlation coefficients, RMSE, bias) should be provided to demonstrate reliability. Without these details, the reproducibility and robustness of the findings remain uncertain. Also, GEE script should be provided. Data validation snow stations map should be added. A workflow figure is recommended. 

We thank the reviewer for their critical and constructive feedback, which highlights key areas for improvement regarding the methodological rigor, reproducibility, and validation of our study. We fully agree that the original description of the Materials and Methods was too general. We have implemented a major revision of the Methodology section to address all points raised and ensure the study's transparency and robustness.

We have expanded the Methodology to provide a full, detailed, and explicit explanation of all processing steps. Specifically, the revised manuscript now includes: 

• a full explanation of the cloud-masking criteria used and the precise spatio-temporal gap-filling/interpolation steps implemented within the Google Earth Engine (GEE) environment.
• The conceptual basis and the mathematical equations for the calculation of iSCA, SP, MDSA, nT.

We recognize that the initial validation using only seven stations was insufficient. The choice of the stations is strictly correlated to the need of a manual check. In fact not many stations of the AOI have this possibility. In facts most of the stations have only automatic check, creating lots of false positives that we cannot manage. To significantly enhance the reliability of our findings, we are implementing two key changes:

As suggested by another reviewer, we have increased the number of in situ snow stations used (to 9 stations) for validation to ensure a more robust and spatially representative evaluation of the MODIS product; inclusion of a comprehensive statistical analysis comparing the MODIS-derived presence/absence of snow cover with the in situ daily observations. In general we already presented some statistical analysis in figure 2 like the bias of 20 days, the correlation R. These were important to observe the solidity of the observation. We know that there is an important Bias, observable also in the other statistical analysis. We decided to insert the MAE value to improve the robustness of the observations.

To dramatically improve clarity and reproducibility, we have added a Workflow Diagram into the Materials and Methods section. This figure visually represents the entire data processing chain, from raw MODIS data ingestion in GEE, through preprocessing and gap-filling, to the final derivation of the nT​ index. A new map is also included, clearly showing the location and elevation of the expanded network of in situ snow stations used for validation.

Results and interpretation



The results section mainly provides descriptive summaries. Trends are mentioned but not quantified in a way that allows meaningful comparison across elevation bands or subregions. The discussion of climatic controls (temperature, precipitation, topography) is superficial and largely speculative, lacking direct evidence from meteorological datasets. The authors should more explicitly link observed snow-cover trends to climate variability and recent warming patterns in the Western Alps.

We appreciate the reviewer’s feedback regarding the need for a more rigorous quantification of the observed trends and a more in-depth discussion linking our snow cover findings to regional climatic controls and warming patterns.

The Results section has been revised to move beyond mere descriptive summaries. We wish to clarify that the quantification of the trends was already implicitly present in the original manuscript through the use of boxplots and the nT​ maps (Figures 4 and 5), which visually and spatially represented the rate of change and variability across the study area and elevation bands. However, we recognize that the explicit textual description of these quantified results was insufficient.

Therefore, the primary revision consists of enhancing the text to explicitly and numerically detail the quantified changes visualized in the figures. We now clearly report the rate of change derived from the nT​ index for different elevation bands and sectors within the manuscript body. This ensures that the magnitude and spatial variability of the observed trends are better articulated and readily accessible for comparison between sub-regions.

We agree that the discussion on climatic controls was superficial. To transform the discussion from being speculative to being evidence-based, we have performed a comprehensive revision integrating findings from published regional climate studies focused on temperature and precipitation trends in the Western Alps. Especially we compare our results with https://doi.org/10.5194/tc-15-1343-2021 and https://doi.org/10.1002/joc.8597. Bozzoli et al. studied the relationship between temperature, precipitations and snowfall in the Alps, with many observation in our study area. They investigated our AOI as SW area. In order to this we extended the discussion on the climatic controls.

These revisions ensure that our conclusions are solidly anchored in both the statistical quantification of the remote sensing data and the established body of regional climate science.

All the changes will be update in the revised manuscript