Mapping seasonal snow melting in Karakoram using SAR and topographic data

Mapping seasonal snow melting is crucial for assessing its impacts on water resources, natural hazards, and regional climate in Karakoram. However, complex terrain in the high-mountain region poses great challenges to remote-sensing-based wet snow mapping methods. In this study, we developed a novel framework that incorporates synthetic aperture radar (SAR) and topographic data for robust and accurate mapping of wet snow over Karakoram. Our method adopts the Gaussian mixture model (GMM) to adaptively determine a wet snow index (WSI) and a computed topographic snow index (TSI) considering the impact of terrain on wet snow distribution to improve the accuracy of mapping. We validated the mapping results against Sentinel-2 snow cover maps, which demonstrated significantly improved accuracy using the proposed method. Applied across three major water basins in Karakoram, our method generated large-scale wet snow maps and provided valuable insights into the temporal dynamics of regional snow melting extent and duration. This study offers a practical and robust method for snow melting monitoring over challenging terrains. It can contribute to a significant step forward in better managing water resources under climate change in vulnerable regions.

1 Introduction

Monitoring seasonal snow melting is of global importance within cryosphere studies, given the profound and far-reaching impacts of snow on climate, hydrology, and ecosystems. Snow cover plays a crucial role in modulating the Earth's energy balance by altering surface albedo, thereby exerting cooling effects on the terrestrial surface and influencing climate patterns at local and regional scales. Notably, in high-altitude areas, the accumulation of snow serves as a primary water source for downstream flows and governs the runoff dynamics in many mountainous basins (Barnett et al., 2005).

The Karakoram region, characterized by its elevated topography and unique climatic conditions, is of exceptional significance in snow cover monitoring. Situated at the center of the Hindukush–Karakoram–Himalaya (HKH) mountain system, Karakoram is renowned for hosting some of the world's highest peaks and harboring the largest alpine glacier system outside the polar regions (Nie et al., 2021). Across its expansive landscape, snow and ice reserves are substantial, encompassing an area exceeding 20 000 km², with seasonal snow covering nearly 90 % of this expanse (Lund et al., 2020; Hasson et al., 2014; Xie et al., 2023).

Figure 1 Geolocation of the Karakoram region overlaid on the COP-30 DEM. The study region includes three major water basins: Hunza, Shigar, and Shyok, which are delineated in red. Elevation histograms of the tree basins are shown in the right panel. Median elevations of basins are indicated with vertical red lines in histograms. Footprints of S1 images used in the study are shown with black and blue boxes (black: relative orbit 27; blue: relative orbit 129).

In the expansive and challenging terrain of the Karakoram region, ground-based observation methods struggle to effectively cover the vast and rugged terrain, and hence remote sensing techniques have become the primary approach for snow cover mapping. While spaceborne optical and multi-spectral sensors like MODIS, Landsat-7/8, and Sentinel-2 (S2) have been employed in numerous studies, their reliability is often compromised by cloud cover and illumination conditions (Immerzeel et al., 2009; Hasson et al., 2014; Tahir et al., 2011, 2015; Khan et al., 2015; Wu et al., 2021). To overcome this limitation, synthetic aperture radar (SAR) presents a valuable alternative for monitoring snow regardless of weather and daylight conditions.

The foundation of snow classification in SAR imagery was established based on the unique backscattering responses generated from the distinctive interactions of snow with SAR signals. For dry snow, radar signals can penetrate through the snowpack down to a specific depth depending on the signal wavelength and thus generate a strong backscattering signal at the snow–ground boundary (Baghdadi et al., 1997; Millan et al., 2015). As the snowpack undergoes melting, its liquid water content increases in the wet snow pack and causes high dielectric losses, resulting in significant reductions in backscattering intensities (Baghdadi et al., 1997; Nagler et al., 2016). Based on the change in backscattering intensities, early wet snow detection methods were developed using the ratio of SAR backscattering coefficients (Baghdadi et al., 1997; Nagler and Rott, 2000). The ratio was derived from SAR images acquired during the snow melting season and reference images obtained under snow-free or dry snow conditions. An empirical threshold of −2 dB on C-band Sentinel-1 (S1) ratio images was proposed to classify snow, distinguishing wet from dry snow, and has proven to be effective (Nagler et al., 2016).

Subsequent refinements were proposed to enhance the algorithm for robust application on various ground surface types. For example, sigmoid functions were introduced as a soft threshold to replace −2 dB, and it was shown to be effective in resolving the uncertainties arising from mixed pixels of wet snow and other constituents (Malnes and Guneriussen, 2002; Longepe et al., 2009; Rondeau-Genesse et al., 2016). Various strategies for selecting bias-free reference images were devised, such as choosing a specific reference scene during winter (Koskinen et al., 1997), averaging multiple images over a defined period (Luojus et al., 2006), or employing linear interpolation between images at the beginning and end of the melting period (Pettinato et al., 2014). In practice, auxiliary data were usually combined to improve the accuracy of snow detection, such as digital elevation models (DEMs), land category maps, meteorological model, and snow cover maps derived from optical multi-spectral sensors (Nagler and Rott, 2000; Liu et al., 2022b, a; Nagler et al., 2016; Tsai et al., 2019). Recent developments in machine learning also brought opportunities to further improve SAR snow mapping. Supervised classification algorithms, such as support vector machine and random forest, were applied to different SAR products and have shown promising results (Tsai et al., 2019; Huang et al., 2011). Deep learning algorithms were also exploited and have exhibited great potential in achieving more robust and accurate wet snow classification (Liu et al., 2023).

Despite the efforts to improve the method for robust SAR-based snow mapping, challenges remain in mountainous regions such as Karakoram, where complex topography may strongly distort SAR signals and thus lead to considerable uncertainties when applying a single-value threshold for snow classification (Snapir et al., 2019; Lund et al., 2020; Karbou et al., 2021b, a). Furthermore, large-scale application of wet snow mapping in Karakoram requires a method to be adaptively responsive to basin-specific variations, posing practical challenges to efficient method development (Rondeau-Genesse et al., 2016). ML techniques may offer versatile solutions, but their application in this region is limited by the availability of training data (Tsai et al., 2019; Liu et al., 2023).

To address these challenges, this study proposes a novel framework that integrates SAR and topographic data for versatile and robust wet snow mapping covering three major water basins of Karakoram. In the first step, we employed an unsupervised learning algorithm, namely the Gaussian mixture model (GMM), to adaptively determine the wet snow index (WSI). Secondly, a topographic snow index (TSI) was calculated to account for the influence of topography on snow distribution. We applied the proposed method to a time series of SAR imagery acquired by Sentinel-1 (S1) between 2017–2021 and generated region-wide wet snow maps of high spatial and temporal resolution. The validation using S2 images showed considerable improvements compared to conventional threshold-based methods. Crucial snow dynamic variables including the wet snow extent (WSE) and snow melting duration (SMD) were derived from the snow maps, demonstrating the significance of closely monitoring wet snow in watershed management.

The paper is organized as follows. Section 2 introduces the study site and data. Section 3 explains the proposed method in detail. Section 4 presents the result of the study, including the validation and the snow dynamic variables. Section 5 discusses the method and its implications for snow mapping. Section 6 concludes the study and provides an outlook for future development.

Table 1 Information on Sentinel-1 (S1) images used for the generation of wet snow maps. The ascending pass crossed over the study region in the late afternoon (around 18:00 local time).

2 Study area and data

2.1 Study area

The Karakoram region, spanning extensively across parts of Pakistan, India, and China, is bordered by some of the highest mountain systems on Earth, including the Himalayas, the Pamir Mountains, and the Hindu Kush mountain ranges (Fig. 1). The study area encompasses the majority of Karakoram, covering an expansive geographical domain of approximately 450 000 km². This region extends from approximately 35 to 38° N latitude and 76 to 78° E longitude, characterized by a wide range of altitudes from around 2000 m a.s.l. (above sea level) in the valleys to well over 8000 m a.s.l. at the highest summits. The extreme topographic variation gives rise to rugged terrain, including steep valleys and towering mountain peaks.

Karakoram is situated upstream of both the Upper Indus Basin (UIB) and the Tarim River Basin. It covers three significant watersheds: Hunza, Shigar, and Shyok, as shown in Fig. 1. These basins serve as the upstream sources of the UIB and contribute 65 %–85 % of annual flows to the Indus River with the melting water from snow and glaciers, sustaining the livelihoods of millions of residents residing within these basins (Archer, 2003; Immerzeel et al., 2009; Forsythe et al., 2012). The hydrological importance of Karakoram emphasizes the important role of mapping snow melting in the region.

2.2 Data

2.2.1 Sentinel-1 SAR imagery

The SAR imagery used in our research was acquired by the spaceborne C-band S1 SAR sensor. The S1 SAR satellite provides high geolocation accuracy and a short orbit repeat cycle of 12 d, facilitating precise and frequent monitoring of snow melting (Nagler et al., 2016).

In this study, we use single-look-complex (SLC) data from S1 acquired in interferometric wide (IW) swath mode. These data cover scenes with a swath width of 250 km at a spatial resolution of $\sim \mathrm{5}\times \mathrm{20}$ m in the range and azimuth direction. Both VV and VH dual-polarization data were employed for the following analysis. The details of the S1 images utilized in this study, including orbit numbers and acquisition dates, are provided in Table 1.

2.2.2 Copernicus digital elevation model

Digital elevation models (DEMs) provide essential topographic information crucial for both SAR image pre-processing and snow distribution mapping, particularly in the rugged landscapes of mountainous regions. In this study, we employed the Copernicus Global 1 arcsec (COP-30) DEM, recently released by the European Space Agency (ESA) in 2020, to facilitate the SAR image processing and snow mapping. The COP-30 DEM product was derived from the TanDEM-X SAR data acquired between 2010 and 2015, providing global coverage at a resolution of 30 m and a vertical root mean square error as low as 1.68 m (Li et al., 2022; Guth and Geoffroy, 2021).

The COP-30 DEM data covering the study regions were downloaded through the Copernicus Space Component Data Access PANDA Catalogue (European Space Agency and Airbus, 2022), as shown in the base map in Fig. 1. The DEM products are referenced in geographic coordinates using the World Geodetic System 1984 (WGS84). The vertical heights are referenced to the EGM2008 geoid model.

2.2.3 Sentinel-2 Level-2A imagery

To validate the snow maps generated using the proposed method, we also derived snow cover maps on selected dates using multi-spectral S2 images. The S2 sensor operates in a sun-synchronous orbit with a revisit time of 5 d (Drusch et al., 2012). Equipped with 13 spectral bands ranging from visible and near-infrared to shortwave infrared, S2 images offer valuable spectral information for land cover characterization and have been widely used in snow cover mapping. In this study, we used the S2 Level-2A (L2A) products to generate snow cover maps for validation purposes. The L2A product is orthorectified bottom-of-atmosphere surface reflectance data that are derived through the atmospheric correction of the Level-1C products using the Sen2Cor method (Main-Knorn et al., 2017). Practically useful supplementary data are also included in the L2A product, including cloud and snow confidence maps and a scene classification map that identifies elements like clouds, cloud shadows, and snow. The spatial resolution of images under different spectral bands varies from 10 to 60 m.

We selected S2 L2A images taken during the summer months with minimal cloud cover as potential candidates. From these, only images that could be matched with an S1 SAR image within a window of ±3 d were used to generate reference snow maps with the Let-it-snow (LIS) algorithm (Sect. 3.5; Gascoin et al., 2019). The LIS algorithm employed RGB spectral bands B2 (blue), B3 (green), and B4 (red), as well as infrared bands B8 (NIR) and B11 (SWIR), for snow–cloud–ice classification. We accessed the S2 data through the Copernicus Open Access Hub.

3 Methodology

This section describes the proposed method for integrating SAR and topographic data to map melting wet snow. The key steps of the method are summarized in Fig. 2.

The raw S1 SLC data were first pre-processed to generate backscatter images. Same-orbit SLC time series were co-registered to a common reference scene to accurately align the geolocation of pixels. Then we multi-looked each SLC with a window size of 12×1 (rg × az) to obtain backscattering intensity images with a squared pixel spacing of approximately 14 × 14 m. The intensity images were further converted to γ⁰ images with terrain-based radiometric correction (Frey et al., 2013). All pre-processing steps were conducted using the GAMMA software (Werner et al., 2000).

Figure 2 The developed method for mapping snow melting. Critical processing steps are shown in the gray box. The COP-30 DEM was used as the reference DEM. GMM stands for Gaussian mixture model.

3.1 SAR backscattering ratio

With the pre-processed γ⁰ images, we derived the composite backscatter ratio (R_c) following the method proposed by Nagler et al. (2016). The R_c metric combines the backscatter ratios from both the VV and VH channels to comprehensively assess surface condition changes associated with snow melting. This approach incorporates a weighting factor W, which is determined by the local incidence angle (LIA, θ), to account for variations in backscatter due to differing incidence angles. This adjustment enables more robust snowmelt detection across varying terrain geometries.

To compute R_c, we first calculated the SAR backscatter ratio R_i for each polarization $i\in \mathrm{vv},\mathrm{vh}$ using the following equation:

$$\begin{array}{}\text{(1)}& R_i={\mathit{\gamma }}_i^{\mathrm{0}}/{\mathit{\gamma }}_{i,\mathrm{ref}}^{\mathrm{0}},\end{array}$$

where ${\mathit{\gamma }}_{i,\mathrm{ref}}^{\mathrm{0}}$ represents a reference image constructed using the average of multiyear winter images. Note that this differs from other alpine regions, such as the Alps, where summer months are often used as the reference due to snow-free conditions. In Karakoram, due to the all-year-long snow cover at the higher elevations, we used winter images as a reference to leverage the dry snow conditions in winter for highlighting the contrast in backscattering intensity between summer wet snow (lower intensity due to water absorption) and winter dry snow.

With R_i for each polarization, the composite ratio R_c was calculated as a weighted sum of R_vv and R_vh:

$$\begin{array}{}\text{(2)}& R_{\mathrm{c}}=W{R}_{\mathrm{vv}}+(\mathrm{1}-W)R_{\mathrm{vh}},\end{array}$$

where the weighting factor W is defined based on the LIA as

$$\begin{array}{}\text{(3)}& {\displaystyle }W{\displaystyle }=\left\{\begin{array}{ll}\mathrm{0}& (\mathit{\theta }<{\mathit{\theta }}_{\mathrm{1}})\\ \mathrm{0.5}(\mathrm{1}+\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_{\mathrm{1}}})& ({\mathit{\theta }}_{\mathrm{1}}\le \mathit{\theta }\le {\mathit{\theta }}_{\mathrm{2}})\\ \mathrm{0.5}& (\mathit{\theta }>{\mathit{\theta }}_{\mathrm{2}}),\end{array}\right.\end{array}$$

with θ₁=20° and θ₂=45°.

3.2 Wet snow index

While R_c alone effectively indicates surface condition changes, it can be sensitive to local variations and does not inherently incorporate adaptive boundaries for wet snow classification. Therefore, instead of directly applying a threshold to R_c for wet snow classification, we propose using a GMM to convert R_c into a wet snow index (WSI) to have a probabilistic measure that better captures the varying conditions of wet snow across different terrains. By leveraging the density distribution of R_c values, the WSI enables a dynamic scaling of the classification based on the underlying distribution of R_c values.

The GMM is an unsupervised probabilistic model for clustering and density estimation. We configured the GMM to identify two clusters corresponding to the wet snow pixels and the no-snow (or dry snow) pixels as

$$\begin{array}{}\text{(4)}& P\left(R_{\mathrm{c}}\right)=\sum _{i=\mathrm{1}}^{\mathrm{2}}{\mathit{\pi }}_i\cdot N\left(R_{\mathrm{c}}\right|{\mathit{\mu }}_i,{\mathit{\sigma }}_i),\end{array}$$

where P(⋅) the probability density function (PDF) of R_c, π_i the ith Gaussian component's mixing coefficient for wet snow (i=1) and no-snow (or dry snow) (i=2) clusters, and $N\left(R_{\mathrm{c}}\right|{\mathit{\mu }}_i,{\mathit{\sigma }}_i)$ the Gaussian distribution with mean μ_i and standard deviation σ_i. We used the full covariance structure in the GMM, i.e., each component has its own general covariance matrix, after testing different types of covariance structures.

To train the GMM and determine its parameters (π_i, μ_i, and σ_i), we randomly sampled 10⁶ unlabeled pixels from the summer R_c images of each basin and employed the expectation-maximization (EM) algorithm to fit the model for each basin with the sampled pixels (Dempster et al., 1977). During the training, the maximum number of EM iterations was set to 100, and the convergence threshold was set to 10⁻³. With the estimated GMM parameters, the WSI can be determined using a modified logistic function:

$$\begin{array}{}\text{(5)}& \mathrm{WSI}={\displaystyle \frac{L}{\mathrm{1}+e^{k(x-x_{\mathrm{0}})}}},\end{array}$$

with

$$\left\{\begin{array}{rl}{\displaystyle }k& {\displaystyle }={\displaystyle \frac{|{\mathit{\mu }}_{\mathrm{1}}-{\mathit{\mu }}_{\mathrm{2}}|}{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}}}\\ {\displaystyle }{x}_{\mathrm{0}}& {\displaystyle }={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{1}}+{\mathit{\mu }}_{\mathrm{2}}}{\mathrm{2}}}\end{array}\right.,$$

where k is the slope factor, x₀ the logistic curve's midpoint, and L the carrying capacity representing the supremum of the function. Here the carry capacity L was set to 10 to amplify the differences between pixels of different R_c values.

In the WSI logistic function, the slope factor k is determined by the separation between the two clusters, providing a flexible and adaptive control over the sensitivity of WSI to the difference between snow conditions. In the case when the two clusters are perfectly distinct from each other (i.e., $|{\mathit{\mu }}_{\mathrm{1}}-{\mathit{\mu }}_{\mathrm{2}}|\gg {\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}$), the WSI is transformed into a step function and thus effectively acts as a single-value hard threshold at x₀. Therefore, Nagler's method can be taken as a special case under this assumption with $x_{\mathrm{0}}=-\mathrm{2}$ dB. Contrarily, the mixed clusters (i.e., $|{\mathit{\mu }}_{\mathrm{1}}-{\mathit{\mu }}_{\mathrm{2}}|\ll {\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}$) would lead to progressively flattened WSI and soft segmentation boundaries. A similar form of logistic function was proposed by Rondeau-Genesse et al. (2016), which was determined with empirical parameters and used as a soft threshold to replace the hard threshold of −2 dB on R_c. In our approach, the GMM allows for an adaptive choice of the parameter k based on the distribution density of R_c, thereby enabling flexible and robust applications at large scale.

Figure 3 Parameters of the GMM and WSI used for (a) Hunza, (b) Shigar, and (c) Shyok. The gray shaded histogram shows the density of R_c of the sampled points. Dashed red and blue curves represent the PDF (scaled on the left y axis) of the two clusters in the GMM. Solid black lines (scaled on the right y axis) show the WSI. The vertical dashed black line marks the center of WSI curves (where R_c=x₀). The mean (μ) and standard deviation (σ) of the GMM are reported above the panel.

3.3 Topographic snow index

Given the strong impact of terrain properties on snow distribution, we introduce the TSI as a component of the proposed method to account for the terrain influence on snow presence.

Terrain properties, including the elevation, slope, and aspect, collectively influence the spatial and temporal distribution of snow. Compared to lower altitudes, regions of higher elevation experience lower temperatures and are conducive to more snow accumulation. The steepness of slopes and the orientation of aspects, on the other hand, impact the snow distribution through solar illumination and wind exposure. To take these factors into account, we calculated the TSI in two steps. First, we derived the discrete topographic bin maps using the COP-30 DEM by partitioning the terrain with its elevation, slope, and aspect. The partition was based on slope below or above 20°, elevation in every 100 m, and aspect in every 15°. The topographic bin map can effectively capture the localized terrain attributes that influence the occurrence of snow. The median WSI value within each topographic bin was then calculated to obtain the TSI so that the regional propensity of wet snow occurrence under the specific topographic conditions could be quantified.

An example TSI distribution for Hunza at different elevation, slope, and aspect is presented in Fig. 4. In this example, TSI values show different patterns across seasons and topographic conditions. In spring, strong TSI signals are found around 4000 m a.s.l. for east-facing slopes (aspect 0∼190°) over flat terrain (slope θ<20°), while no obvious TSI signal is observed for steep terrain (slope θ>20°). This can be explained by the limited snowmelt during the spring season of Karakoram. In summer, strong TSI signals are observed above ∼ 5000 m a.s.l. for all slope and aspect conditions, indicating the presence of wet snow. However, steep slopes showed an unevenly distributed TSI across slope aspects, which can be attributed to the shadowing effect of the surrounding terrain. In autumn, TSI signals generally decrease due to the absence of wet snow, and the snow line retreats to higher elevations. This example indicates that the dynamic influence of topography on snowmelt can be effectively captured by the designed TSI signal.

Figure 4 TSI values at different elevation (y axis), aspect (x axis), and slope classes (a–c: flat slope with θ<20°, d–f: step slope with θ<20°) for the Hunza basin. Three observation dates in spring (a, d), summer (b, e), and late autumn (c, f) are displayed. Aspects of 0, 90, 180, and 270° align with the north, east, south, and west direction, respectively.

3.4 Integrated snow index

As illustrated in Fig. 2, the final step generated an integrated snow index (SI) map through pixel-wise multiplication of WSI and TSI. This multiplication scales the WSI by incorporating terrain characteristics, thereby linking the observed SAR backscattering ratio directly to terrain properties.

In order to classify the integrated SI into binary snow maps, it is crucial to apply an adaptive threshold that accounts for the variation in topographic features across different basins. The variation in SAR backscatter response within a basin is inherently handled by the GMM when deriving the WSI. In contrast, the TSI is time-varying and basin-specific, requiring an optimal coefficient to condition the SI for classification. This coefficient should moderate the influence of TSI in the SI formulation and compensate for the concentration of TSI values within a limited range so that the impact of TSI can be normalized to align with the traditional −2 dB reference. To determine this coefficient, we performed a sensitivity analysis, evaluating F₁ score, precision, and recall across different values using the S2 validation snow map. The results (Fig. A1) demonstrate that Hunza and Shyok exhibit similar responses, with optimal coefficients close to 3.5, while Shigar reaches its optimum at approximately 2.5. However, to avoid overfitting to specific basins or validation dates, we selected 3.5 as an overall coefficient to balance classification performance across all basins. This coefficient also reflects a moderate threshold applied to the TSI to determine the overall SI threshold for each basin.

The threshold was then calculated using the following equation:

$$\begin{array}{}\text{(6)}& \text{SI Threshold}=\mathrm{3.5}\times \mathrm{WSI}{|}_{R_{\mathrm{c}}=-\mathrm{2}},\end{array}$$

where $\mathrm{WSI}{|}_{R_{\mathrm{c}}=-\mathrm{2}}$ represents the WSI at a backscatter ratio of $R_{\mathrm{c}}=-\mathrm{2}$ dB for each basin. This value is basin-specific, allowing the threshold to adapt based on each basin’s distinct characteristics. Together, these conditions form an integrated, basin-adaptive thresholding mechanism, combining SAR backscatter and topographic information into a single index to determine the SI threshold.

Note that while the SI threshold is basin-specific, it is time-independent. The WSI is derived from a GMM trained on samples collected from multiple summer scenes over several years, ensuring that it represents an aggregated measure for each basin and is not tied to individual scenes or seasons. This design ensures robustness to seasonal variations in liquid water content and enables consistent application across different validation dates.

3.5 Sentinel-2 snow cover maps

The proposed method was validated using S2 snow cover maps generated following the LIS algorithm proposed by Gascoin et al. (2019). Before running the LIS algorithm, the input S2 multi-spectral bands were resampled to a pixel size of 20 m×20 m to match the resolution of different bands. The COP-30 DEM was also resampled to the same pixel size as the S2 images.

The LIS algorithm started by generating provisional snow masks by applying thresholds on the normalized difference snow index (NDSI) and the red band reflectance (ρ_red) with the condition

$$\begin{array}{}\text{(7)}& (\mathrm{NDSI}>n_i)\mathrm{and}({\mathit{\rho }}_{\mathrm{red}}>r_i),\end{array}$$

where n_i=0.4 and r_i=0.2 (Gascoin et al., 2019). This step was designed to identify snow-covered areas while excluding non-snow surfaces such as vegetation and bare ground. However, this approach could sometimes exclude snow-covered pixels due to errors in cloud masking. To correct the errors, a refinement step was introduced to reassign dark cloud pixels that were initially misclassified. Following Gascoin et al. (2019), dark clouds were identified by applying a threshold of 0.3 on the bilinearly downsampled red band, which reduced the resolution of the red band from 20 to 240 m by a factor of 12. This process helped to exclude cloud shadows and high-altitude cirrus clouds from the snow classification. Afterwards, the provisional snow masks were further refined using the basin snow line calculated from the COP-30 DEM. We calculated the total snow cover fraction (SCF) within every 100 m elevation band using the provisional snow mask and defined the snow line using the lowest elevation band where the SCF exceeded 30 %. For pixels identified as dark clouds above the determined snow line, the conditions used in Eq. (7) were reapplied with adjusted thresholds to account for the unique conditions at high altitudes. Specifically, the relaxed thresholds of n_i=0.15 and r_i=0.04 were used to classify snow pixels, and dark cloud pixels with ρ_red>0.1 were reassigned as cloud, while other pixels were categorized as “no snow”. These adjusted thresholds help to differentiate snow from dark clouds in challenging high-altitude environments, ensuring a more accurate classification. Following the adjustment of the snow mask, we extended the LIS algorithm by further applying a threshold on the NIR band with ρ_NIR>0.4 to classify glacier ice and water bodies from snow (Paul et al., 2016).

4 Results

4.1 Validation of snow classification maps

The validation of snow classification results was carried out for three basins on specifically suited summer dates (Table 2). These dates were chosen based on conditions of minimal cloud cover and the shortest possible intervals between acquisitions of S1 and S2 images. As S2 snow cover maps classify snow (both dry and wet) rather than only wet snow, we have chosen only the summer dates as listed in Table 2 to ensure that the S2 snow cover maps were mostly covered by wet snow. For Hunza and Shyok, two S2 images with an acquisition interval of 2 d were combined to achieve complete coverage of the basin, whereas same-day acquisitions of S1 and S2 were found for Shigar in 2019.

In the three basins, adaptive SI thresholds were used to generate the snow classification maps. The threshold values for each basin are also reported in Table 2. These thresholds provided basin-adaptive and time-independent classification boundaries to distinguish wet snow and dry snow or snow-free pixels.

Table 2 Sentinel-1 (S1) and Sentinel-2 (S2) data as well as SI thresholds used for the validation. Different acquisition dates and adaptive SI thresholds were used for basins.

Figure 5 (a) R_c, (b) integrated SI, (c) SI classification map, and (d) S2 snow cover map (as a reference) for all three basins (Hunza, Shigar, and Shyok). The S2 snow map was generated using the Let-it-snow (LIT) algorithm, as described in Sect. 3.5. Zoomed insets provide a closer view of selected locations in each basin, highlighting performance differences on glacier surfaces.

Figure 5 shows maps of R_c and SI, as well as the SI snow map and the S2 snow map over the three basins. Compared to the R_c map, the SI map shows a much clearer boundary that separates wet snow pixels from no-snow or dry snow pixels. Over glacier surfaces and valley regions (shown with the zoomed-in insets), R_c falls in the value range of −4 to 0 dB, making it sensitive to the choice of threshold values. The R_c map over these regions presents noisy patterns, likely due to the complex scattering mechanisms on glacier surfaces. Over glacier surface, especially in the ablation zone, radar backscatter responses are highly variable due to refreezing, supraglacial features (e.g., crevasses, sun cups, debris cover), and the presence of wet debris or bare ice (Scher et al., 2021). These features contribute to significant spatial variability in backscatter within a single pixel, making it challenging to distinguish dry snow, wet snow, and ice with R_c alone. Additionally, glacier movement between winter and summer scenes introduces further variability, compounding the uncertainty in detection.

In Fig. 5b, it is observed that the variation of SI values is within different ranges across the three basins. This is due to the different density distributions of R_c, WSI, TSI, and SI in the three basins as shown in Fig. 7, underscoring the necessity of applying dynamically adaptive thresholds for each basin to generate classification results.

Figure 6 Detailed comparison between the snow classification map obtained using SI (left column), R_c (middle column), and the S2 snow cover map (right column, reference) for Hunza (top row), Shigar (middle row), and Shyok (bottom row) basins. Class labels are indicated in the color bar. Note that the selected regions were located in the mid-elevation range of the area, and hence the snow class in both S1 and S2 results refers to wet snow.

Figure 7 Density distribution of R_c, WSI, TSI, and SI in three basins (left: Hunza, middle: Shyok, right: Shigar). The x axis of R_c is inverted as negative R_c corresponds to a higher probability of wet snow.

The snow maps in Fig. 5c and d offer a visual comparison between the snow maps produced by the proposed method and the S2 LIS algorithm. A detailed comparison is further illustrated in Fig. 6. Overall, the two methods show good agreement, demonstrating the effectiveness of the proposed approach. Notably, the SI-based classification exhibits reduced false positives and cleaner snow-free classifications in valleys and over glacier surfaces, as evident in the detailed comparison. This improvement is primarily attributed to the incorporation of topographic information through TSI, which smooths SI values within each topographic bin. However, a consistent mismatch in the ice/water category can be observed between the SAR (both the SI and ratio methods) and S2 results, particularly over glacier surfaces. This discrepancy arises from the differing detection principles of the two approaches: the S2 results classify glacier ice using thresholds on the NIR band, while the SAR-based methods do not explicitly resolve glacier ice. On the observed date, glacier ice in the ablation zone may have partially melted, reducing the SAR backscatter ratio and leading to its misclassification as wet snow in the SAR results. As discussed earlier, glacier surfaces present unique challenges for SAR-based methods due to their complex scattering mechanisms. While the inclusion of TSI improves the robustness of our method by integrating topographic controls, it does not explicitly account for the heterogeneity of glacier surfaces. Consequently, glacier-specific conditions, such as localized melting or scattering from mixed ice–snow surfaces, may lead to underestimation or misclassification.

We further quantified the comparison using a confusion matrix and the F₁ score, as listed in Table 3. In the confusion matrix, S2 snow-free (S2-SF) and S2 snow (S2-S) pixels are used as negative and positive labels, respectively, while S1 snow classification maps designate no-snow or dry snow (S1-N/D) and wet snow (S1-WS) pixels as their counterparts. To minimize errors caused by the presence of dry snow at high altitudes, areas above 5500 m a.s.l. were excluded from the calculation. This elevation threshold, representing 11.08 % of pixels in Hunza, 12.95 % in Shigar, and 30.43 % in Shyok (19.60 % in total), was chosen because such areas are consistently classified as non-melting zones due to persistently low air temperatures that inhibit snowmelt. Across all three basins, the proposed method shows an improvement in classification performance. In Hunza and Shyok, both true negative and true positive rates increased, while in Shigar, the true negative rate improved by 0.11, though the true positive rate decreased slightly by 0.04. Despite this, the overall F₁ scores improved for all three basins, highlighting the method’s enhanced precision and recall.

Figure 8 Snow coverage profile along elevation bands. Elevation bands were binned every 100 m, and the snow coverage within a band was calculated based on the aggregated snow pixel percentage within every band. Numbers reported in legends measured the mean absolute error (MAE) between the profile curve of SAR wet snow classification and the S2 snow cover map.

Table 3 Confusion matrix and F₁ score between S1 snow classification maps and S2 snow cover maps. Snow-free (S2-SF) and snow (S2-S) pixels in S2 snow cover maps correspond to the negative (N) and positive (P) labels, respectively. The associated labels in S1 snow classification maps are no-snow or dry snow (S1-N/D) pixels and wet snow (S1-WS) pixels. Results of the proposed method are highlighted in bold, whereas the R_c-threshold-based results are reported in normal font.

We also evaluated the accuracy of the classification map using the elevation profile of snow distributions. As illustrated in Fig. 8, the profile of snow coverage was analyzed along 100 m elevation bands. Below 4500 m a.s.l., snow coverage was overestimated by approximately 7 % when using the R_c thresholds, leading to a misrepresentation of the actual snow line, especially over the challenging conditions on glacier surfaces and valleys (Fig. 5). In contrast, the SI method provided a snow classification closer to the S2 profiles in the low-elevation regions. Between 4500 and 5500 m a.s.l., the SI curve exhibited a noticeably steeper slope than the S2 curve, with an underestimation of snow coverage between 4000 and 5000 m and an overestimation from 5000 to 5500 m. This pattern suggests that the SI uncertainties in mixed snow conditions within transition zones might have led to a nonlinear exaggeration of the TSI response to snow cover. A more precisely calibrated TSI model could further align snow coverage profiles between the SI method and S2 results. Above 5500 m a.s.l., where expansive dry snow cover predominates, both R_c and SI maps showed a greater reduction in wet snow coverage compared to S2, highlighting the differing sensitivities of SAR signals to dry snow conditions in these methods compared to S2's multi-spectral data.

Figure 9 Temporal and elevation dependence of WSE in (a) Hunza, (b) Shigar, and (c) Shyok, as well as (d) the total WSE of the three basins within the studied period. (e) The temperature and precipitation record obtained from the ERA5 dataset. In panels (a)–(c), the WSE is represented by the color scale. In panel (d), the blue line indicates the weekly average temperature of the three basins, and the shaded blue indicates the range between the weekly maximum and minimum temperature. Monthly averaged precipitation is shown with the green bar. The shaded gray zone shows the summer of the year. The WSE was calculated as the ratio of wet snow-covered area to the total area within each respective elevation band.

4.2 Temporal dynamics of snow melting

In this section, we further applied the proposed method to the collected S1 image time series between 2017–2021 to generate wet snow maps across all three basins. The generated wet snow maps were resampled from the original SAR image resolution of 14×14 to 30×30 m before the analysis to speed up the processing. The temporal interval of the time series was 12 d, the same as the acquisition interval of S1 images. We analyze two key properties of snow cover dynamics derived from the time series data: the wet snow extent (WSE) and the snow melting duration (SMD). It is worth noting that these properties represent just a subset of the potential insights that can be derived from this dataset.

4.2.1 Wet snow extent

The temporal patterns and elevation dependencies of WSE across the Hunza, Shigar, and Shyok basins are depicted in Fig. 9a–c. WSE was calculated as the percentage of wet snow pixels within each 100 m elevation band, offering a granular view of snowmelt progression. Over the 5-year period analyzed, a consistent interannual pattern was observed in all three basins. Melting is typically initiated by the end of March to early April and is concluded by late September to early October. As temperatures rose from spring (i.e., April to May) into summer, the melting front (e.g., the upper and lower elevation boundary of the melting area) ascended along the altitude gradient. Specifically, the lower elevation boundary of wet snow extended upwards as snow at lower altitudes fully melted, while the upper boundary extended to higher altitudes as higher temperatures resulted in melting at greater elevations. In the peak melting months (i.e., July and August), the upper boundary of melting snow reached its maximum altitude before descending, whereas the lower boundary extended to its highest extent and stabilized by the end of the melt season.

Figure 9d–e illustrate the interactions between total WSE, temperature, and precipitation within the region using data from the ERA5 reanalysis dataset (Hersbach et al., 2020). The temperature data, averaged weekly, include mean, maximum, and minimum temperatures of air at 2 m above the surface of land across the Karakoram region, providing insights into the thermal conditions influencing snow melting. The precipitation data are the accumulated liquid and frozen water falling to the Earth's surface. They were compiled and averaged monthly to complement the temperature analysis by revealing precipitation trends and their impact on snowpack. While the three basins demonstrated similar interannual variability, annual discrepancies were pronounced within the time series. For instance, the peak WSE in 2018, around 40 %, was notably lower than the approximately 50 % observed in other years. This reduction in WSE is related to diminished winter precipitation during the 2017–2018 season, as indicated by the precipitation data. The onset of snow melting aligned well with the period when maximum temperatures rise above freezing, suggesting that peak temperatures were a more sensitive indicator for the onset of snow melting than mean temperatures. These complex interannual fluctuations underscore the snowpack's responsiveness to immediate weather conditions, such as temperature spikes and precipitation events.

4.2.2 Snow melting duration

The SMD reflects the temporal persistence of wet snow cover within a given year, allowing for consistent comparisons across years with varying numbers of observation days. To compute the SMD for each year, we first determined the ratio of days with wet snow cover (M) to the total number of observed days (N) for each pixel. Since the number of observation days (N) varied each year and was typically less than 365, we re-scaled this ratio to a 365 d basis using the formula $(M/N)\times \mathrm{365}$ to standardize the annual average of wet snow cover days and enable consistent comparisons between years.

Figure 10 presents the annual average SMD across the study region from 2017 to 2021. The average SMD displays a pronounced terrain dependency: valley areas at lower altitudes typically exhibit an SMD of fewer than 60 d, while higher altitudes, such as glacier accumulation zones, generally experience SMD exceeding 120 d. In certain high-altitude regions, the wet snow cover can persist for more than 180 d annually.

Figure 10 Annual average SMD of the study region, calculated based on the 5-year observation.

To assess the temporal and spatial dynamics of SMD, we evaluated the annual fraction of SMD for each basin, as depicted in Fig. 11. SMD was categorized into four ranges: 0–60 d (blue), 60–120 d (orange), 120–180 d (green), and 180–240 d (red). In Hunza, the area with an SMD of less than 60 d saw a decline from 2017 to 2019, followed by an increase from 2020 to 2021. This change was inversely related to the 60–120 d category, which expanded from 2017 to 2019 before contracting. The fraction of WSD exceeding 120 d initially decreased from 2017 to 2019 and then increased as the 60–120 d category diminished. A pronounced peak within the 180–240 d range occurred in 2017. The Shigar basin exhibited more pronounced annual oscillations in SMD. The area with an SMD below 60 d fluctuated around an average of 0.5, while the 60–180 d category mirrored the pattern in Hunza, increasing from 2017 to 2019 before a subsequent decline. The SMD range of 120–180 d remained relatively stable at about 0.2, with periods exceeding 180 d noted only in 2017 and 2020. The Shyok basin experienced the most substantial temporal variation in SMD. The 0–60 d category showed strong variance around the 0.5 level. The 60–120 d category peaked in 2019, representing a larger fraction (∼ 0.5) compared to those in Hunza (∼ 0.3) and Shigar (∼ 0.35). The 180–240 d range was present exclusively in 2017.

Figure 11 Fraction of snow melting duration in the three basins between 2017–2021. SMD was segmented into four categories as represented by different colors.

5 Discussion

5.1 Classification performance

Compared to the conventional methods where only a single-value threshold is used on the R_c map, the proposed method effectively improved the mapping accuracy in the validation. This is primarily attributed to the transformation of R_c into WSI and the incorporation of TSI.

The GMM enabled a data-driven approach that allowed adaptively transforming the R_c into WSI based on the local SAR signal responses. As shown in Fig. 3, the WSI function for each basin was characterized by distinct center (x₀) and slope (k) parameters determined from the GMM, indicating that varied local responses of SAR backscattering were raised by the diverse nature of wet snow distribution in different basins. This provides the flexibility required for large-scale application and thus offers an approach for robust wet snow mapping in complex and diverse landscapes.

Furthermore, the proposed method enriched the terrain analysis in snow mapping by incorporating multiple topographical factors (i.e., elevation, slope, and aspect) beyond the traditional use of snow line elevation. The add-on information from slope and aspect enables accurate capture of snow spatial distribution over the complex terrain, which is particularly important in the Karakorum region. As exemplified in Fig. 4 for the Hunza basin, the TSI values vary with elevation, aspect, and slope across different times of the year. At the onset of snowmelt, flat and east-facing slopes at altitudes between 3500 and 5000 m a.s.l. exhibited higher TSI values than the other regions, suggesting that snow in such areas was more susceptible to melting, likely due to solar exposure. As the season progressed into summer, altitude increasingly dictated the snowmelt on flat terrains, while on steeper slopes, aspect also played a significant role. Approaching late autumn, the snow line stabilized at similar altitudes across the two slope classes, yet variations were observed across the aspect with higher TSI values on south-facing slopes (approximate aspect of 75–255°), indicating conditions more conducive to melting. This level of detail in our analysis demonstrated the potential of our method to provide a more comprehensive understanding of wet snow dynamics than analysis using only snow line elevations.

While the proposed method demonstrates strong performance across the three basins in Karakoram, several limitations should be acknowledged. First, the choice of the SI threshold relies on the coefficient determined through sensitivity analysis for the selected study areas. While this approach balances classification performance across the basins, it is not fully dynamic and may require adaptation for larger, more topographically diverse regions, such as global-scale applications. Future work could explore supervised learning models, such as random forests or neural networks, to capture more complex, nonlinear relationships between SAR backscatter, topography, and snow conditions, enabling dynamic threshold adaptation.

Second, glacier surfaces present unique challenges for SAR-based snow classification. As shown in the results, glacier-specific scattering mechanisms, including contributions from wet debris, bare ice, and supraglacial features, introduce variability in radar backscatter that is not explicitly modeled in the current method. This limitation may lead to underestimation of wet snow on glacier surfaces and highlights the need for further refinement of the method to better handle glacier-specific conditions, potentially by incorporating land surface type information.

Finally, the use of TSI introduces potential bias in regions where topographic conditions deviate significantly from the assumptions underlying its calculation. As shown in Fig. 7, the long-tailed distribution of TSI values reflects the cumulative statistical nature of TSI, which relies on using a single median value to represent each topographic bin. This approach may be inadequate in bins with strong terrain variations, especially when the TSI distribution of pixels within a bin is highly skewed. In such cases, using the median can lead to systematic overestimation or underestimation – overestimating in left-skewed distributions and underestimating in right-skewed ones – ultimately affecting precision and recall in the classification results (see Sect. S2 of the Supplement for detailed visualizations). For instance, low TSI values in ablation zones may lead to underestimation of wet snow, as noted in the comparison with S2 results. While the inclusion of TSI improves overall robustness by integrating terrain characteristics, future studies should evaluate strategies to mitigate these biases, especially in regions with complex topography or unique land surface characteristics.

It is also important to note that while we followed Nagler's method (Nagler et al., 2016) to generate the R_c image, we did not apply the same post-processing steps, such as median filtering and land cover masking. These smoothing and filtering steps may influence accuracy, and future work could incorporate them to further evaluate their impact and potentially improve snow mapping performance.

5.2 Implications of wet snow maps

Large-scale wet snow maps, especially the ones with high spatial and temporal resolution, have significant implications for hydrological studies, water resource management, and climate impact assessments (Helmert et al., 2018). Snow data obtained from remote sensing and field site stations have proven to be fundamental for the development, calibration, and validation of snowpack, hydrology, and runoff prediction models (Schmugge et al., 2002; Andreadis and Lettenmaier, 2006; Dressler et al., 2006; Griessinger et al., 2019; Cluzet et al., 2024).

Using the wet snow maps generated from the proposed method, our study has extracted and analyzed two critical snow variables, i.e., WSE and SMD, which are crucial for understanding regional snow melting dynamics. The analysis of WSE uncovered detailed patterns of snowmelt changes over time and across elevations, which can provide valuable observations for calibrating snowpack models or forecasting runoff events (Cluzet et al., 2024). The interpretation of SMD highlighted the yearly differences in snow melting duration across the basins, with Hunza exhibiting relative stability and Shyok demonstrating the most variability. A long-term SMD observation record will provide key insights into the changes of the regional climate pattern.

6 Conclusions

In this study, we proposed a novel approach for mapping wet snow in complex mountainous regions, such as Karakoram, by effectively combining S1 SAR data with topographic information. We first adopted the GMM to adaptively transform the SAR backscattering ratio R_c into the WSI as a robust representation of wet snow under complex surface conditions. Then, we introduced the TSI to capture the likelihood of snow presence influenced by topographic conditions. Validation with S2 snow cover maps demonstrated a notable improvement in the accuracy of wet snow classification.

With the collected time series of S1 images over the three major water basins in Karakoram, we produced large-scale wet snow maps using the proposed method. The wet snow maps have enabled detailed analysis of crucial snow variables including the WSE and SMD. Analysis of the two variables revealed the dynamic pattern of the temporal–spatial distribution of wet snow in Karakoram, suggesting that the comprehensive dataset produced with this study can offer further enhancement for hydrological model calibrations and validation, thereby ensuring informed water resource management and climate modeling.

Future work involves integrating the approach with in situ observations and hydrological models to further improve the accuracy and utility of water resource planning tools. Continuing to advance this research would provide results that are greatly beneficial for fostering climate resilience and sustainability in Karakoram.

Appendix A: Sensitivity analysis of the SI threshold

We applied a sensitivity analysis to TSI coefficients to find the optimal SI threshold. A series of TSI coefficients were tested using validation images to examine how the TSI coefficient in the SI threshold affected the classification results. Three metrics, including the F₁ score, precision, and recall, were evaluated and used as the selection criteria. The result (Fig. A1) showed that the optimal coefficients for Hunza and Shyok were close to 3.5, while a lower coefficient of around 2.5 was found to be optimal for Shigar.

Figure A1 Sensitivity analysis of the choice of TSI coefficients for SI thresholds. Panel (a) shows the F₁ score, and panel (b) is for precision (blue) and recall (red). Dashed lines are results using $R_{\mathrm{c}}=-\mathrm{2}$ dB as the threshold. The dotted vertical line indicates where the coefficient is 3.5.