Research Methodology

Overall Research Framework

The study proposes an integrated approach used for the detection, identification, and analysis of temporal tendencies of snow cover in Hohe Tauern National Park, Austria, in 2018–2025. It includes methods used for multi-sensor and multi-source data analysis: multi-spectral optical data from Sentinel-2, SAR backscatter from Sentinel-1, topographical data as a Digital Elevation Model (DEM), and gridded climate data from the SPARTACUS reanalysis product. The general process consisted of four stages: (i) pre-processing of the input data and calculation of relevant spectral indices; (ii) analysis of the combined Sentinel-2 and Sentinel-1 data using Principal Component Analysis (PCA); (iii) creation of a rule-based algorithm for the classification of Non-Snow, Dry Snow, and Wet Snow; (iv) analysis of the elevation-based, seasonal, and climatic trends in snow cover, and accuracy assessment using AUC–ROC metrics.

Data Collection

Sentinel-2 Level-2A multispectral surface reflectance and Sentinel-1 GRD SAR data were downloaded from the Copernicus Open Access Hub and processed in Google Earth Engine. For each year of the study period, two seasonal composites were generated for each sensor: the snow accumulation period (December–February) and the snow depletion period (June–September). The Sentinel-2 data provided the spectral reflectance values in the visible (Blue B2, Green B3, Red B4), near-infrared (NIR B8), and shortwave infrared (SWIR1 B11) domains. For the Sentinel-1 SAR data, the C-band backscatter coefficients in VV and VH polarizations were used in Interferometric Wide (IW) swath mode. The 12.5 m ALOS PALSAR DEM data provided slope and aspect layers for altitude stratification. SPARTACUS reanalysis data from Geosphere Austria provided daily temperature and precipitation at 1 km resolution. MODIS (500 m) and Landsat 8/9 (30 m) NDSI products were used as independent references for accuracy assessment.

Data Preprocessing

Sentinel-2 Level-2A imagery was accessed as atmospherically corrected surface reflectance products. Cloud-contaminated pixels were removed using the QA60 quality assurance band. Cloud-free seasonal composite images were generated using the median reducer within GEE for each seasonal window. Sentinel-1 GRD products were filtered to include only IW swath mode acquisitions in ascending orbit configuration. Raw backscatter intensity values were converted to the logarithmic decibel (dB) scale, and seasonal median composites of VV and VH polarization backscatter were computed to reduce the influence of speckle noise. All optical and SAR datasets were resampled to a uniform spatial resolution of 10 m and co-registered to WGS84 / UTM Zone 33N. Prior to PCA, all input bands were subjected to mean-centering and z-score standardization independently to eliminate scale-related biases.

Spectral Indices

The Normalized Difference Snow Index was calculated as NDSI = (B3 − B11) / (B3 + B11), exploiting the high reflectance of snow in the visible spectrum and its low reflectance in the SWIR region. The Normalized Difference Vegetation Index was computed as NDVI = (B8 − B4) / (B8 + B4) and incorporated into the multi-sensor stack to minimize spectral confusion between snow and vegetated surfaces.

Principal Component Analysis

To evaluate the integrated information content of the combined optical and SAR dataset, an eight-band input stack was assembled: five Sentinel-2 spectral bands (B2, B3, B4, B8, B11), NDVI, and two Sentinel-1 backscatter bands (VV, VH). NDSI was deliberately excluded from the PCA input stack to avoid methodological circularity in the subsequent validation step. Following z-score normalization, PCA was performed using the getPrincipalComponents function within GEE at 10 m spatial resolution. Pearson's correlation coefficient was computed between each principal component and the independently derived Sentinel-2 NDSI layer to identify the component most effectively capturing snow-related spectral variability.

Rule-Based Snow Classification

A hierarchical rule-based classification algorithm was developed to categorize each pixel into one of three classes: Non-Snow, Dry Snow, or Wet Snow. In the first stage, pixels with NDSI ≥ 0.4 were classified as snow-covered; pixels below the threshold were assigned to Non-Snow. In the second stage, snow-covered pixels were differentiated into Dry Snow and Wet Snow using two SAR criteria: (i) an absolute backscatter threshold VV < −11 dB, and (ii) a relative change criterion ΔVV < −2.0 dB relative to the dry-season baseline. Pixels satisfying either criterion were labeled Wet Snow; all remaining snow pixels were Dry Snow. A minimum mapping unit filter removed clusters of fewer than ten connected pixels.

LST, Change Detection & Elevation Analysis

Land Surface Temperature (LST) was derived from Landsat 8/9 thermal infrared imagery (Band 10, TIRS). Seasonal mean LST composites were generated for each accumulation and depletion period from 2018 to 2025. Inter-annual snow cover change was quantified using pixel-wise post-classification comparison between consecutive annual pairs. Elevation-stratified time series (Zone I <1500 m, Zone II 1500–2000 m, Zone III 2000–3000 m, Zone IV >3000 m) were computed for each snow class, and slope and aspect layers were used to analyse terrain-orientation effects.

Climatic Correlation & Accuracy Assessment

SPARTACUS gridded daily temperature and precipitation were spatially aggregated for each seasonal window. Pearson's r was computed between climatic variables and the fractional extent of each snow class within each elevation zone, with statistical significance assessed via two-tailed t-tests at the 0.05 level. Accuracy was evaluated using AUC–ROC analysis comparing NDSI products derived from Sentinel-2 (10 m), Landsat 8/9 OLI (30 m), and MODIS (500 m); the optimal NDSI threshold was identified via the Youden index (TPR − FPR).

⤢ Click to zoom