Glacier change at the ISM-MLW transition: multi-decadal evidence of spatial heterogeneity from the Sara Umga watershed, western Himalaya

1 Introduction

Glaciers are among the largest reservoirs of freshwater on the planet and play a crucial role in sustaining downstream hydrology. Outside the polar regions, the Hindu Kush Himalaya (HKH) region contains the greatest concentration of glaciers (Bolch et al., 2012), serving as vital sources of water for densely populated regions of south Asia (Azam et al., 2021). Given their sensitivity to climatic fluctuations, glaciers are widely recognized as one of the most reliable indicators of climate change (Haeberli et al., 2007). Limited ground-based surveys and extensive remote sensing observations have confirmed that glaciers are thinning at a rapid pace in the western Himalaya (Hugonnet et al., 2021; Das and Sharma, 2021; Bhambri et al., 2023), altering runoff regimes and enhancing cryospheric hazards (Singh et al., 2021; Bolch et al., 2019; Kumar et al., 2024; Das et al., 2024). Long-term monitoring of glacier fluctuations is therefore indispensable not only for anticipating regional water and energy security but also for understanding their contribution to global sea-level change (Scherler et al., 2011). Over the past century, average air temperature in the Himalaya has increased by ∼0.1 °C per decade (1901–2014), accelerating to ∼0.2 °C per decade between 1951 and 2014 (Sabin et al., 2020). This sustained warming has driven widespread glacier recession and permafrost thaw since the mid-19th century across western Himalaya (Rowan, 2017; Nüsser and Schmidt, 2021; Deswal et al., 2025; Burn et al., 2025).

Unlike glaciers in the Alps and Andes, Himalayan glaciers commonly occur in steep, tectonically active terrain and are extensively debris-covered, leading to complex surface morphology and fundamentally different ablation regimes and dynamic behaviour. Active tectonics and geology typically contribute copious quantities of debris to Himalayan glaciers (Scherler et al., 2011). If debris input from hillslope processes exceeds the glacier’s capacity to transport or remove it, the glacier surface becomes increasingly debris-covered (Anderson et al., 2018). This is especially evident in the ablation zone, where progressive vertical thinning exposes englacial debris to the surface (Benn et al., 2012; Rowan et al., 2015). In western Himalaya, studies have reported that, alongside accelerated ice loss, there has also been a notable increase in debris cover (Bhambri et al., 2023). Debris cover exerts a strong control on ablation as thin layers tend to accelerate surface melt (Anderson and Anderson, 2016), whereas debris thicker than a few centimetres insulates the ice beneath and suppresses melt rates (Reznichenko et al., 2010). In theory, the expansion of thick debris should lower thinning rates compared to debris-free ice (Anderson et al., 2021).

Influence of debris cover extends beyond glacial melt, as it also affects glacier dynamics in many ways, such as its impact on ice thickness, flow velocity, and surface heterogeneities, i.e., ice cliffs and supraglacial lakes (Racoviteanu et al., 2022). Debris-covered glaciers of the western Himalaya respond differently to climate warming compared to debris-free ones (Chand et al., 2020; Mandal and Sharma, 2020). As these glaciers are sensitive to changes in thickness and velocity, estimates based solely on length or areal fluctuations often fail to capture the true extent of ice loss (Das and Sharma, 2019; Chand et al., 2020). The impact of enhanced mass loss, reflected by numerous supraglacial lakes and ice cliffs accelerates localized ablation and alters glacier dynamics (Bhambri et al., 2023). This study stems from the recognition that debris-covered valley glaciers respond to a combination of topographic, climatic, and intra-glacial forcing, creating a chain of feedback mechanisms (Das et al., 2022; Rowan et al., 2015). Unlike clean-ice glaciers, where terminus retreat or areal change may indicate mass loss, debris-covered glaciers require multi-parameter assessments. Studies that integrate areal and terminal changes with mass loss, and velocity attenuation provide a far more comprehensive understanding of their dynamics (Shukla and Garg, 2019; Bhambri et al., 2023).

Studies have investigated glacier changes in surrounding basins such as the Chandra (particularly Bara Shigri), Parvati, Bhaga and the Miyar catchment (Chand et al., 2017; Chand et al., 2020; Deswal et al., 2017; Mandal and Sharma, 2020; Das and Sharma, 2019; Das et al., 2023). Comparable basin-scale studies are inadequate for the Sara-Umga Watershed (SUW). This is a critical gap as the SUW lies at the periphery of the ISM-influenced Parvati Basin and the MLW-dominated Chandra Basin, making it sensitive to shifts in regional climate. Existing reports on the SUW are contrasting in nature. Dobhal and Kumar (1996) provided the first length estimate of the GL01(∼17.5 km), derived from 30 m resolution Landsat Thematic Mapper imagery acquired between 1988 and 1989. Bahuguna and Kulkarni (2004) estimated a frontal retreat of ∼304 m between 1963 and 2003, whereas the Geological Survey of India’s (GSI) ‘Inventory of Himalayan Glaciers’ documented an average retreat rate of ∼41 m yr⁻¹ over 1963–2004 (Sangewar and Shukla, 2009). This inconsistency highlights a significant knowledge gap, particularly given the valley’s geographic proximity to some of the most extensively studied glaciers in the western Himalayas and its hydrological importance as a sub-basin where multiple hydropower projects are already established.

The present study provides the first comprehensive reconstruction of glacier change in the SUW. The objectives of this study are: (a) reconstructing glacier recession in terms of length, area, and debris cover; (b) quantifying geodetic thinning and surface velocities, with particular attention to velocity slowdown on debris-covered glacier tongues; (c) attributing patterns of change to local topographic factors (i.e., hypsometry, slope, aspect, and debris cover) and climatic factors. Additionally, we discuss the hydrological and hazard implications, particularly with regard to downstream water availability, the expansion of proglacial and supraglacial lakes, and associated cryospheric risks. By integrating multi-temporal satellite datasets, DEM differencing, and glacier velocity analyses, this research not only situates SUW within the broader context of western Himalayan glacier change but also offers insights into the drivers of recession and attenuation in a marginal climatic zone (i.e., ISM and MLW).

2 Study area

The SUW is situated in Kullu District, Himachal Pradesh, within the western Himalaya (Figure 1). It forms a sub-catchment of the Parvati basin, which drains into the Beas River near Bhuntar, and lies in the sub-humid middle Himalaya (Chand et al., 2020). The SUW hosts a total of 54 glaciers (>0.05 km²) in 2024 (Supplementary Figure S1), the largest being the Sara Umga Glacier (GL01; 17.4 km in length), Badgan-I Glacier (GL39; >11.6 km in length), and Tichu Glacier (GL28; >11.9 km in length), respectively (Figure 1). The SUW is drained by the Tosh Nala, originating from the GL01(∼4022 m a.s.l.) and joined by several tributaries downstream.

In 2024, the SUW had a total glacierized area of 128.38 ± 2.98 km², amongst which 24.84 ± 1.15 km² is classified as debris-covered glacier surface, concentrated mainly in the glacier ablation zones (Table 1). A majority of the glaciers (n = 14) are attenuated towards the southwestern aspect, covering a cumulative glacierized area of 75.59 ± 2.91 km² (Supplementary Figure S2). The location of the basin at the interface of two major climatic regimes, ISM and the MLW, makes it highly sensitive to seasonal variability, with summer precipitation dominated by the ISM and winter snowfall driven by MLW system (Bookhagen and Burbank, 2010).

TABLE 1. Glaciers parameters according to different size classes for all glaciers.

The basin is primarily underlain by the Vaikrita Group of the Higher Himalayan Crystalline (HHC) Unit, with the Kharo Formation exposed in most of the area, except for a small south-eastern intrusion of the Morang Formation from the same group (Bhukosh, GSI). The Kharo Formation is composed mainly of Schist, Gneiss, Migmatite, Quartzite, and Marble, while the Morang Formation is characterized by Kyanite-bearing Schists and Quartzite. Within the valley, foliated-Biotite-rich Quartzite forms the predominant lithology (Khan et al., 2021). The region is tectonically active, as evidenced by numerous lineaments, predominantly with a NE-SW orientation, along with the presence of the Vaikrita Thrust in the southern sector (Khan et al., 2021).

3 Data and methods

3.1 Image collection and rectification

To delineate glacier extents in the SUW, we employed multi-temporal and multi-spectral satellite data spanning over 3 decades (Supplementary Table S1). Landsat 5 Thematic Mapper (TM) image from 1994 was pan-sharpened using the Brovey transformation fusion method and subsequently co-registered to a Sentinel-2 base image through projective transformation (Mandal and Sharma, 2020). We have also included Landsat 7 ETM+ (2000), Landsat 5 Thematic Mapper (TM; 2010), and Sentinel-2B (2024). All datasets were aligned to a common reference using the Auto Sync Workstation in ERDAS Imagine (Yuan et al., 2013). The most recent Sentinel-2B image (2024) provided the baseline for mapping current glacier outlines. For glacier mapping, shortwave infrared bands 11 and 12 (20 m) were resampled to 10 m resolution to enhance accuracy (Das et al., 2023).

3.2 Glacier area and terminus change

Automatic and semi-automatic delineation techniques are effective for clean-ice zones but remain unreliable for the debris-covered areas (Bhambri et al., 2011). Since the majority of the glaciers in the SUW are debris-mantled, outlines were manually digitized by on-screen vectorization (Mandal and Sharma, 2020). Visual indicators such as emergent meltwater streams, variations in spectral tone and texture, ice cliffs, and exposed ice faces were used to refine glacier boundaries in the debris-covered ablation zones (cf. Bhambri et al., 2011; Chand and Sharma, 2015). Outlines from 2024 served as a reference for adjustments across earlier periods. Areal changes were calculated by assessing multi-year outlines, while terminus positions were mapped using two methods: (a) using the central flow line and (b) parallel transects at 50 m intervals along former (1994) and contemporary (2024) snout boundaries (Figure 2). The second method has been adopted to account for the irregularities of the glacial snout undergoing melt. GL01 showed considerable inconsistencies in length change, and to understand this spatial variability, we have used two methods, and we reported both mean retreat distances as average frontal change (Chand and Sharma, 2015).

A glacier inventory was produced from Sentinel-2B (2024), adopting a minimum mapping threshold of >0.05 km² following RGI standards (Frey et al., 2012). Contiguous ice bodies were separated into individual units using hydrological watersheds derived from the Shuttle Radar Topographic Mission (SRTM) 1arc second within the GIS platform. Polygons were validated against 3D visualization in Google Earth and attributed with topographic parameters (i.e., elevation, slope, aspect). To assess consistency, our 2000 outlines were compared with the GAMDAM (Sakai, 2019; Nuimura et al., 2015), RGI v7 (RGI 7.0 Consortium, 2023) and ICIMOD (Bajracharya and Shrestha, 2011) inventories (Supplementary Figure S5).

To ensure reliability, uncertainties for each parameter were estimated separately. Mapping error was calculated using a buffer half the sensor’s pixel size (Das et al., 2023): 15 m for Landsat TM, 7.5 m for ETM+, and 5 m for Sentinel-2. This yielded overall mapping uncertainties <5%, consistent with earlier studies (Chand and Sharma, 2015; Mandal and Sharma, 2020). Terminus change uncertainty (U T e r m i n u s) was computed using the formula from Hall et al. (2003):where a and b are the pixel resolutions of images 1 and 2, and R E G E R R O R is the registration error. The estimated uncertainties were 45.80 m for 1994, 38.67 m for 2000 and 23.13 m for 2010 with respect to the base image (2024).

U T e r m i n u s = a 2 + b 2 + R E G E r r o r (1)

3.3 Surface velocity extraction

Glacier surface velocities (v) were derived using the Co-registration of Optically Sensed Images and Correlation (COSI-Corr) tool in ENVI (Leprince et al., 2007). This pixel-tracking algorithm has been successfully applied in the western Himalaya (Das et al., 2022). We used a frequency correlation algorithm for image correlation. The frequency correlation approach in COSI-Corr requires defining several input parameters (Das and Sharma, 2021). For this analysis, we used an initial window size of 64 pixels and a final size of 32 pixels, with a step size of 4 × 4 pixels, yielding an effective ground resolution of ∼60 m. Ideally, the selected window size should be at least twice the magnitude of the expected displacement, but in the case of the SUW no prior field constraints were available. To minimize high-frequency noise and improve the robustness of the correlation results, a frequency mask threshold of 0.9 and two robustness iterations were adopted, consistent with values successfully tested in the adjacent valleys in western Himalaya (Chand et al., 2020; Das and Sharma, 2021; Das et al., 2022). Using these settings, we processed ten annual image pairs from 1993 to 2022 (Supplementary Table S1), covering the late-ablation period (August-October) (Supplementary Figure S3).

COSI-Corr generates three primary outputs: East-West Displacement (EWD), North-South Displacement (NSD), and a Signal-to-Noise Ratio (SNR) image (Chand et al., 2020; Das and Sharma, 2021). Total ground displacement D was derived by combining the two orthogonal displacement components through the root-sum-square method (Equation 2). Surface velocity v was then calculated by dividing displacement by the time interval t between image acquisitions (Equation 3). The SNR product was further used to assess the quality and reliability of the displacement fields.

D = E W D 2 + N S D 2 (2)

v = D / t (3)

Corrections were applied for decorrelation noise, orbital errors, jitter distortions, and topography-dependent artifacts (Das and Sharma, 2021). Low SNR values (<0.9), cloud/shadow zones, and miscorrelated vectors were masked or removed (Millan et al., 2019). Polynomial detrending and directional filters were applied, followed by 3 × 3-pixel median filtering (Copland et al., 2009; Ding et al., 2016).

To evaluate image correlation errors, we adopted the procedure of Scherler et al. (2008), which is based on two diagnostic checks: (i) analysis of displacements over stable, off-glacier terrain assumed to remain motionless during the observation period, and (ii) assessment of displacement vectors to ensure alignment with glacier flow streamlines. The uncertainty in measured velocities over non-glaciated areas was then quantified using the formulation described by Berthier et al. (2003):

σ o f f = S E 2 + M E A N 2 (4)

Here, σ o f f represents the error of off-glacier velocity, MEAN signifies the mean stack velocity of random points, SE denotes the standard error of the mean velocity, calculated as follows:

S E = S T D N o f f (5)

In this context, STD refers to the standard deviation of off-glacier area velocity, while N o f f represents the number of measurements in the off-glacier area. Following Das and Sharma (2021), we derived ∼10,000 random points from the off-glacier area and evaluated off-glacier displacements based on the frequency distribution curve. Mean and standard deviation values of off-glacier displacement were computed for each annual and seasonal image pair.

3.4 Glacier elevation change

Glacier surface elevation change (Δh) was quantified using globally consistent satellite-derived elevation change products developed by Hugonnet et al. (2021). These datasets are available at a spatial resolution of 100 m and provide multi-temporal estimates of glacier surface elevation change (Hugonnet et al., 2021). In this study, Δh rate maps corresponding to two decadal periods (2000–2010 and 2010–2020) were used (Supplementary Figure S7). The associated standard error layers were also incorporated to assess uncertainty. As this study directly utilizes the final, quality-controlled Δh rate grids provided by Hugonnet et al. (2021), no additional DEM processing steps (generation, co-registration, bias correction) were required.

To derive glacier-wide thickness change, each glacier was subdivided into 25 m elevation bands (e.g., 4500–4525 m a.s.l., 4525–4550 m a.s.l.). Mean Δh values were calculated for each band to represent elevation-dependent thinning or thickening; and the glacier-wide mean Δh was obtained by averaging these band-wise values across the full elevation range. Data voids were present in some upper accumulation zones; therefore, both the total number of pixels and the number of valid pixels contributing to the mean Δh were quantified for each elevation band to assess data robustness. Spatial patterns of thinning were further analysed by overlaying Δh maps with surface velocity, ice-cliff and supraglacial lake inventories, as well as slope and aspect layers.

Estimating uncertainty in glacier surface elevation change is challenging due to spatial autocorrelation, sensor noise, and natural variability in glacier surfaces. No universally accepted approach exists (Vijay and Braun, 2016; Bhambri et al., 2023). While some studies use standard deviation or root mean square error (RMSE) (Bolch et al., 2011; Gardelle et al., 2013), others employ normalized median absolute deviation to reduce sensitivity to outliers (Höhle and Höhle, 2009). In this study, uncertainty estimation follows the conservative approach proposed by Gardelle et al. (2013), which is widely used in glacier mass balance studies.

The uncertainty in pixel-wise elevation change U_Δh was assumed to be equal to the standard deviation σ_Δh of Δh values over stable, non-glaciated terrain within each elevation band:

U Δ h i = σ Δ h i (6)

Depending on altitude, σ_Δh reached values of up to ±40 m, particularly in higher accumulation zones. This approach incorporates both measurement uncertainty and natural spatial variability, thereby providing a conservative estimate.

Glacier-wide uncertainty was calculated as:

U Δ h = 1 N ∑ i = 1 n U Δ h i 2 (7)

Where U Δ h i is the uncertainty for the i t h elevation band, n is the total number of elevation bands, N is the effective number of independent observations. The effective number of independent measurements N was calculated as:where, N t o t is the total number of pixels within the elevation band, r is the spatial resolution of the Δh grid (100 m), d is the spatial autocorrelation distance. An autocorrelation distance of 510 m was adopted, consistent with previous studies in the Himalayan region (Bhattacharya et al., 2021; Bhambri et al., 2023). This formulation accounts for spatial dependency among neighboring pixels and avoids overestimation of the effective sample size.

N = N t o t · r 2 d 2 (8)

3.5 Hypsometric index

The Hypsometric Index (HI) provides an effective morphometric measure to evaluate the distribution of glacierized areas with elevation, offering insights into glacier stability and sensitivity to climate forcing. It is calculated as the ratio of mean glacier elevation relative to its maximum and minimum elevations (Jiskoot et al., 2009; King et al., 2017; Robson et al., 2018).

HI = H max ‐ H med H med ‐ H min and if 0 < HI < 1 then HI = ‐ 1 HI (9)

The glaciers were then classified into five groups: very top-heavy (HI < −1.5), top-heavy (−1.5 < HI < −1.2), equidimensional (−1.2 < HI < 1.2), bottom-heavy (1.2 < HI < 1.5), and very bottom-heavy (HI > 1.5). These categories distinguish glaciers with extensive accumulation zones from those dominated by elongated, ablation-heavy tongues (Jiskoot et al., 2009; King et al., 2017; Robson et al., 2018; Deswal et al., 2023).

3.6 Ice cliff density and supraglacial lake mapping

To examine whether surface features such as ice cliffs and supraglacial lakes correspond with zones of enhanced ice flow and thinning, we prepared an inventory of these features within the frontal 5 km of the GL01 (measured from the snout). Mapping was carried out manually using Sentinel-2B imagery (2024) and cross-checked against high-resolution Google Earth images for validation. For mapping the lakes, semi-automated processes (e.g., normalized difference water index), cannot be applied for the lakes as they were very small (mean area <0.003 km²) in size and were mapped by manual digitization. All identified features were converted into vector format for subsequent analysis. Ice cliffs were further processed by applying a line-density operation, which was then converted into a raster layer for direct comparison with glacier surface velocity from 2021 to 2022. For spatial analysis, we also plotted the mapped ice cliffs with every 100 m elevation bands, for comparison with velocity and thickness data.

3.7 Climate analysis

Long-term climatic trends were analysed using Climate Research Unit TS (v4.09) gridded monthly precipitation and temperature (0.5 ° × 0.5 °) (Harris et al., 2020) data spanning from 1901 to 2024 for the grid 32.25° N and 77.75° E. Monthly records were aggregated into seasonal indices: Pre-monsoon (April-June: AMJ), monsoon (July-September: JAS), post-monsoon (October-November: ON) and winter (December-March: DJFM) seasonal periods were considered, following previous such studies (Jain et al., 2009; Rani and Sreekesh, 2018). Statistical significance of trends was tested using the non-parametric Mann-Kendall test (Mann, 1945; Kendall, 1970). To see the magnitude/slope in the trend, Sen’s slope estimator (Sen, 1968) has been employed.

4 Results

4.1 Area change

The glacierized area of the SUW shows a consistent decline during the observation period (Table 2), accompanied by progressive fragmentation of larger valley glaciers. The total glacier area decreased from 132.92 ± 9.28 km² in 1994 to 128.38 ± 2.98 km² in 2024, representing a net loss of ∼4.54 ± 0.04 km² (∼3.4%). Area loss accelerated after 2010, with ∼2.94 km² lost between 2010 and 2024, compared to ∼0.70 km² (1994–2000) and ∼0.91 km² (2000–2010). Among the large glaciers, GL01 exhibited the greatest absolute area reduction, decreasing from 35.30 ± 1.72 km² in 1994 to 33.97 ± 0.55 km² in 2024, accounting for ∼30% of the total basin-wide loss (Supplementary Figure S4-a5). Comparable but smaller reductions were observed for GL28, GL39 (Supplementary Figure S4-d5), and GL44 (Supplementary Figure S4-c5).

TABLE 2. Absolute change in glacial area and debris cover between the observation periods (1994–2024).

Increase in glacier number results from fragmentation of larger glaciers.

Uncertainties derive from a buffer of half the sensor pixel size and therefore decrease with improving image resolution.

Comparison with earlier inventories supports our calculations for 2000 (Supplementary Figure S5): GAMDAM (134.55 km², 2001) and RGI v7 (133.73 km², 2002) exceed our calculated glacierized area for 2000 (132.23 km²) by 1.8% and 1.1%, respectively. This deviation is well within the calculated mapping uncertainty of ±6.9% (Equation 1). These small positive offsets from the GAMDAM and RGI v7 sources arise from the placement of ice divides (Supplementary Figure S5a) in the upper accumulation areas, as well as the inclusion of seasonal snow-patches (Supplementary Figure S5b) as glacier ice. In comparison, ICIMOD (2005) reports 108.90 km² of glacierized area in the SUW, i.e., 23.33 km² less than our calculated area for 2000. The debris-covered glacierized area in the present study for 2000 was calculated at 14.51 km², which appears to have been omitted from the ICIMOD (2005) estimate of the total glacierized area in the catchment (Supplementary Figure S5a). Omission of debris cover from the total glacierized area calculation is a major contributor alongside other factors, including differences in ice-divide delineation, coarseness of vectorization, and minimum area thresholding parameters.

Glacial fragmentation increased after 2000, particularly in the ablation zone of large debris-covered glaciers. Many glaciers split into smaller units, including GL01A, GL22A, and GL24A (Supplementary Figure S1), which have been recorded as independent glaciers in later (post- 2010) inventories. These fragmented glaciers are small and exhibit minimal area change after detachment. The emergence of these fragmented units shows structural disintegration, indicating physical separation of the stagnant tongues from the active trunk sections.

Debris cover increased significantly over the past 3 decades, reflecting a fundamental change in the ablation zone (Table 2). In 1994, debris covered 13.67 ± 2.26 km², increasing slightly to 14.51 ± 2.80 km² by 2000 and 15.04 ± 1.54 km² by 2010. During this early period (1994–2010), debris-cover expansion was relatively modest, with rates of 0.14 km² yr⁻¹ (1994–2000) and 0.05 km² yr⁻¹ (2000–2010), indicating gradual exposure of englacial debris driven by surface lowering. In contrast, the period 2010–2024 shows a pronounced acceleration, with debris-covered area increasing sharply by 9.80 ± 0.011 km² at a rate of 0.70 km² yr⁻¹. By 2024, debris cover reached 24.84 ± 1.15 km², representing a net increase of 11.17 ± 0.014 km² since 1994 and an average long-term expansion rate of 0.37 km² yr⁻¹. This nonlinear growth highlights the increasing dominance of debris-covered ablation zones as glacier thinning intensifies.

4.2 Length changes and terminus retreat

Glacier length change was quantified using two independent approaches; measurement along the central flowline, and averaging across parallel strips spaced 50 m apart along the snout boundary (Section 3.2). For the nine glaciers analysed, the two methods agree closely (r = 0.966), with no statistically significant systematic bias between them. Agreement of these two methods improves with increasing glacier size. For example, for the largest glacier (GL01), results differ by only 10 m (1.3%), whereas the smallest glaciers (GL24, GL25) differ by 36%–49%, which reflects a greater sensitivity of short, irregular termini to the choice of measurement transect. Retreat is therefore reported as the mean of both methods, with uncertainty combining the co-registration error (±45.80 m for 1994–2024) with half the method spread for each glacier (Supplementary Table S3). Variability in glacial retreat rate indicates uneven frontal response, consistent with debris cover heterogeneity and local topographic controls.

Between 1994 and 2024, all glaciers analysed show a pattern of steady frontal retreat along the entire snout. GL22 recorded the largest cumulative retreat (1469.5 m; 48.98 m yr⁻¹), equivalent to approximately 32.3% of its 1994 length, followed by GL01 (780.4 m; 26.0 m yr⁻¹), GL24 (703.0 m; 23.43 m yr⁻¹), GL28 (436.0 m; 14.53 m yr⁻¹), GL39 (472.3 m; 15.74 m yr⁻¹), GL44 (421 m; 14.05 m yr⁻¹) and GL07 (301.5 m; 10.05 m yr⁻¹). The smallest retreats were recorded by GL25 (280.8 m; 9.36 m yr⁻¹) and GL09 (145.3 m; 4.84 m yr⁻¹). The disproportionate retreat of GL22 relative to its size is consistent with its minimal debris cover (<1% of its 2024 area) and steep surface gradient (∼26.4°), which together promote rapid frontal recession.

4.3 Surface velocity

Surface velocities of the major glaciers in the SUW were quantified for the period 1993–2022, revealing pronounced temporal variability and glacier-specific flow regimes. During 1993–94, basin-wide mean velocity was relatively high (∼23 m yr⁻¹), with GL01 (29.5 m yr⁻¹) and GL07 (27.0 m yr⁻¹) exhibiting the fastest flow, indicating dynamically active glacier trunks. Longitudinal velocity profiles of GL01 (Figures 3b, 4A) show that peak velocities were concentrated in the mid-ablation zone, particularly near the transition between clean ice and debris-covered ice. Here, the steep surface gradients and greater ice thickness possibly enhance ice flux. In contrast, the debris-covered lower tongue exhibited reduced velocities, signalling early dynamic decoupling (Figure 4A).

The cross-profile (Figure 3a) velocity measurements along GL01 reveal pronounced lateral and longitudinal variability in surface ice velocities. Profiles L1-L4 (Figure 3a; Figure 4B) show relatively high velocities (≈25–37 m yr⁻¹) in the central parts of the glacier, reflecting focused ice flux along the main flowline. A systematic decline of surface velocity towards the margins can be observed where lateral drag and debris cover increase. In contrast, profile L5 exhibits consistently lower velocities (≈3–20 m yr⁻¹), indicating a more stagnant, debris-dominated section of the glacier tongue. The asymmetric velocity patterns across profiles suggest spatially heterogeneous responses to basal, englacial, and surface debris conditions, especially in the ice-marginal regions (Benn et al., 2012; Rowan et al., 2015), alongside ice thickness.

By 2000–01, the mean basin velocity declined to 20.6 m yr⁻¹, driven by a marked reduction in surface velocities of GL28 (17.4 m yr⁻¹) and GL44 (12.4 m yr⁻¹), particularly in their lower reaches. GL01 and GL07 maintained near-constant velocities. This slowdown intensified during 2009–2011, when basin-wide velocities reached a minimum of ∼8.8 m yr⁻¹, and every glacier recorded its lowest mean velocity of the observation period. In GL01, velocity reduced from 30.8 m yr⁻¹ (2000–01) to ∼17 m yr⁻¹ (2009–11), with near-stagnant conditions developing toward the snout, where the terminal 2 km averaged only 4.6 m yr⁻¹. GL07 recorded the most pronounced reduction, from 26.9 to 5.3 m yr⁻¹, whereas GL28 (7.0 m yr⁻¹) and GL39 (5.6 m yr⁻¹) sustained comparatively low velocities.

During 2020–22, the basin-wide mean velocity partially increased to 14.2 m yr⁻¹, although spatially uneven. GL01 again recorded the highest velocity (19.2 m yr⁻¹), followed by GL07 (16.1 m yr⁻¹), GL44 (13.9 m yr⁻¹), GL39 (10.8 m yr⁻¹), and GL28 (10.6 m yr⁻¹). In GL01, higher velocities were confined to the clean-ice and transition zones near the approximate equilibrium-line altitude (ELA), where velocities peaked at 31.3 m yr⁻¹ (6–10 km from the snout), while the debris-covered lower tongue remained comparatively sluggish (11.0 m yr⁻¹ within 2 km of the terminus; Figure 4A).

Over the full observation period (1993–2022), GL01 maintained the highest long-term mean surface velocity (23.4 m yr⁻¹), followed by GL07 (19.0 m yr⁻¹) and GL28 (14.4 m yr⁻¹), while GL44 (11.9 m yr⁻¹) and GL39 (11.4 m yr⁻¹) recorded the lowest long-term means. Although velocities increased partially after the 2009–11 minimum, the basin-wide mean remained well below early-1990s values (22.9 m yr⁻¹ in 1993–94 versus 14.2 m yr⁻¹ in 2020–22), indicating a net deceleration of ∼38% over 3 decades. These velocity patterns reflect a transition from a dynamically connected valley glacier system, increasingly falling prey to debris-mediated stagnation. It is presumed to reduce ice flux at the lower elevations, causing spatial decoupling between active upper trunks and near-stagnant glacier termini (Bhambri et al., 2023).

4.4 Elevation change

Glacier-wide surface elevation change (Δh; reported as rates, m yr⁻¹) for GL01 during 2000–2020 reveals a strong altitudinal control on thinning and thickening patterns. The lower ablation zone (3900–4400 m a.s.l.) experienced the most consistent and pronounced thinning, with mean elevation-change rates ranging from −1.22 ± 0.95 to −1.67 ± 0.39 m yr⁻¹ (between L5(a)-L5(b) of Figure 3a). Maximum thinning within these bands locally exceeded −3.0 m yr⁻¹, indicating sustained mass loss near the glacier tongue (Figure 3C). Thinning remains substantial up to ∼4300–4400 m a.s.l. (−1.45 ± 0.40 m yr⁻¹), highlighting the sensitivity of the debris-covered lower tongue to enhanced melt processes (Figure 4D; Supplementary Figure S6).

Between 4400 and 4700 m a.s.l., the glacier shows a transition zone characterized by high variability in Δh. Mean values in this zone are close to zero, but the standard deviations of Δh are 2–8 times larger than the means (e.g., 0.36 ± 2.25 m yr⁻¹), indicating that this transition zone is dominated by scatter rather than a robust thickening or thinning signal. Sixteen pixels with |dh/dt| > 5 m yr⁻¹ in this elevation range are unlikely for surface-elevation-change rates and were treated as MMASTER processing artifacts; band means were recomputed with these pixels masked. The residual scatter is consistent with heterogeneous debris cover, ice cliffs, supraglacial ponds, and dynamic ice inputs from upstream areas. Apparent mean thickening is observed between 4700 and 4800 m a.s.l. (+1.70 ± 1.22 m yr⁻¹ after artifact masking), but 36% of the pixels in this band were rejected during quality control, and this value should be treated as tentative. Above ∼4800 m a.s.l. (between L3(a)-L3(b) of Figure 3A), elevation changes progressively diminish, with mean Δh values oscillating around zero (−0.07 ± 0.54 m to 0.06 ± 0.33 m yr⁻¹) up to ∼5600 m a.s.l., indicating near-equilibrium conditions. The upper accumulation zone (>5700 m a.s.l.) exhibits minimal net change, with mean values close to zero and relatively low variability (≤±0.92 m yr⁻¹), reflecting limited melt and reduced dynamic thinning at higher elevations. Above ∼6000 m a.s.l., Δh values are effectively negligible, though the small pixel count in these elevation bands limits robust interpretation.

4.5 Hypsometric index

Retreat distances reported in this section are derived from the central flowline method (Section 3.2). Among the studied glaciers, GL01 exhibited a bottom-heavy distribution (HI = 1.25), and retreated 780.37 m along the central flow line between 1994 and 2024 (at a rate of ∼26.01 m yr⁻¹) (Figure 5A). The very bottom-heavy glacier GL24 (HI = 1.64) with the steepest mean slope of ∼26.9° displayed a frontal retreat of 702.99 m (at ∼23.43 m yr⁻¹). The equidimensional GL39 (HI = −1.07) displayed a retreat of 472.31 m in the same time frame, with a retreat rate of 15.74 m yr⁻¹ (Figure 5F). However, GL22 (Figure 5C), with HI of 1.09, showed the highest retreat among the sampled glaciers (1469.46 m at a rate of 48.98 m yr⁻¹). The rapid retreat of GL22, despite its classification as an equidimensional glacier, can be attributed to its minimal debris cover (<1% of its 2024 area) and one of the steepest surface gradients (∼26.4°) among all the glaciers analysed (Figure 5c). Local factors such as debris cover, supraglacial features (ice cliffs and lakes), slope, and glacier size likely modulate responses, sometimes overriding hypsometric controls (Deswal et al., 2023). The very top-heavy GL44 (HI = −1.86) and top-heavy GL07 (HI = −1.48) show snout retreats of 421 m and 301.55 m at rates of 14.05 m yr⁻¹ and 10.05 m yr⁻¹ respectively (Figure 5e; Figure 5b).

HI provides an important first-order indicator of glacier health by describing the distribution of glacierized area with respect to elevation (Jiskoot et al., 2009; King et al., 2017). In theory, glaciers with highly negative HI values (top-heavy or very top-heavy) retain much of their area in the accumulation zone (Equation 9), therefore are expected to be more resilient to climatic warming, while equidimensional or bottom-heavy glaciers are expected to respond more rapidly to rising temperatures (Davies et al., 2012) (Table 3). However, results from the SUW suggest a more complex relationship, with coupling of multiple factors governing glacier behaviour. The findings suggest that while HI captures the overall susceptibility of glaciers, it does not fully account for the variability in retreat rates (Figure 5). Considering the small sample size, hypsometric analysis should be regarded only as indicators of tendency rather than a direct controlling relationship governing glacial retreat (Robson et al., 2018).

TABLE 3. Topographical details of selected glaciers.

Z_min, Z_max, Z_mean and Z_med are the minimum, maximum, mean and median elevation of the glaciers, respectively.

HI, is the hypsometric index calculated after Jiskoot et al., 2009; Robson et al., 2018.

LCCF: Length change in m along central-flow line (1994–2024) ± 45.80.

ALCS: Average length change in m along strips at glacial snout boundary (1994–2024) ± 1.05.

4.6 Supraglacial features and ice cliff density for the Sara Umga Glacier

The inventory of supraglacial features in the frontal zone of GL01 highlights their strong role in shaping localized ablation and ice dynamics. A total of 201 ice cliffs were mapped, with lengths ranging from ∼10 m to ∼500 m (Figure 6b; Figure 6c) and were accounted within 100 m elevation bands (Figure 3d). This wide variability suggests a heterogeneous melt environment, where numerous small to medium cliffs contribute to localized surface lowering, while a few exceptionally large cliffs dominate frontal ablation. Similarly, 33 supraglacial lakes were identified, ranging from ∼62 m² to ∼2088 m². The findings indicate that while most lakes are small (<500 m²), a handful of large, persistent water bodies account for a substantial proportion of supraglacial storage. When compared with velocity data (2021–2022), some zones of higher ice cliff density and large supraglacial lakes coincide with locally pre-eminent surface velocities. This spatial association is tentative pending verification, as the ice-cliff inventory is confined to the frontal 5 km of GL01. The velocity data reflect that this is the slowest-flowing section of the glacier. Enhanced meltwater input from these ablation hot spots may plausibly influence local ice flow through englacial meltwater pathways, which could contribute to localized acceleration, even in otherwise stagnant, debris-covered ice.

To quantify the influence of debris cover on sub-debris melt, elevation-band (100 m) dependent debris thickness and melt factor data for GL01 were extracted from Rounce et al. (2021). Debris thickness decreases systematically from ∼2.95 m near the snout (3900–4100 m a.s.l.) to ∼0.01 m above 4800 m a.s.l., while the corresponding sub-debris melt factor increases inversely from ∼0.04 to ∼1.36 (r = −0.92) (Supplementary Figure S6). Thick debris at the lower tongue therefore suppresses ablation to a fraction of the debris-free rate, consistent with the theoretical expectation that supraglacial debris exceeding a few centimetres in thickness insulates the underlying ice (Østrem, 1959; Reznichenko et al., 2010). The transition zone between 4200 and 4800 m a.s.l., where debris thins from ∼1.28 to ∼0.08 m, is characterised by high spatial variability in both thickness and melt factor. The standard deviation in this section exceeds the mean values in five out of six elevation bands, reflecting the patchy and heterogeneous debris distribution documented against the velocity and thinning analyses above (Supplementary Table S2). Interestingly, the zone of thickest debris and most strongly suppressed melt (3900–4200 m a.s.l.) coincides with the zone of maximum surface lowering (Section 4.4). Sub-debris melt here is reduced to a small fraction of the clean-ice rates; the observed thinning cannot be attributed primarily to surface ablation alone, and instead reflects dynamic thinning arising from the reduced ice flux documented in Section 4.3.

Despite this debris-controlled melt suppression, localised velocity peaks coincide with increased ice cliff density. Our analysis reveals (Figure 6) a moderate positive correlation (r = 0.50) between ice cliff density and GL01 glacier surface velocity. Zones with higher ice cliff densities (>2 km) generally exhibit enhanced velocities, in some cases exceeding 25 m yr⁻¹, though the ice-cliff inventory zone (frontal 5 km) overlaps with the slowest section of the GL01 and this association should be interpreted with caution. Ice cliffs act as localized ablation hot spots (Racoviteanu et al., 2022), intensifying melt and generating meltwater that may contribute to accelerated flow through subsurface pathways (Zoet et al., 2013), which was not directly observed in this study.

4.7 Climate

Long-term climate trend analysis using the Mann-Kendall test and Sen’s slope estimator for the period 1901–2024 reveals a pronounced seasonal reorganization of temperature and precipitation in the region, despite relatively stable annual precipitation totals. Annual mean precipitation shows a weak positive but statistically non-significant trend, indicating no robust long-term increase in total moisture input (Figure 7e). In contrast, seasonal precipitation trends are strongly differentiated: pre-monsoon (AMJ) precipitation exhibits a significant increasing trend (Z = 3.21; 99% confidence) at a rate of 0.354 mm yr⁻¹ (Figure 7a), while post-monsoon (ON) precipitation also increases significantly (Z = 2.99; 99% confidence), though with a smaller magnitude (0.070 mm yr⁻¹) (Figure 7c). Conversely, winter precipitation (DJFM) shows a statistically significant declining trend (Z = −2.75; 99% confidence) at −0.160 mm yr⁻¹, suggesting a long-term reduction in snowfall. Monsoon-season (JAS) precipitation displays a weak and non-significant negative trend. Temperature trends show a clearer and more consistent warming signal, with statistically significant increases in annual mean temperature (0.006 °C yr⁻¹) as well as in post-monsoon and winter seasons (both ∼0.010 °C yr⁻¹), with a 99% confidence level (Figure 8). In contrast, pre-monsoon and monsoon temperatures show positive but non-significant trends. The analysis shows a shift toward warmer winters, reduced snowfall-dominated precipitation, and increased rainfall during shoulder seasons (AMJ and ON), fundamentally altering the regional cryo-hydrological regime and conditions consistent with sustained negative mass balance.

5 Discussion

Between 1994 and 2024, the total glacierized area in the SUW reduced from 132.92 ± 9.28 km² to 128.38 ± 2.98 km², corresponding to a net loss of 4.54 ± 0.04 km² (−3.42%) (Table 2). Although this magnitude of area loss is modest compared to some Himalayan basins, its temporal distribution is highly uneven. Nearly 65% of the total basin-wide loss (∼2.94 km²) occurred after 2010, while earlier periods were characterized by relatively small reductions (−0.70 ± 0.026 km² during 1994–2000 and −0.91 ± 0.019 km² during 2000–2010). Comparable rates and post-2010 acceleration have been reported from nearby basins such as Chandra (Patel et al., 2021), Bhaga (Das et al., 2022), and Parvati (Chand et al., 2020; Tak and Keshari, 2020), indicating a regionally coherent glacier response to recent climatic forcing across the western Himalaya.

Beyond net area loss, glacier evolution in SUW is strongly characterized by progressive fragmentation of formerly continuous valley glacier units, particularly after 2010. The detachment of glacier units such as GL01A, GL22A, and GL24A (Supplementary Figure S1) reflects sustained thinning, stagnation, and dynamic decoupling of lower glacier tongues from their accumulation zones. These fragmented units persist as small, low-gradient ice bodies with near-zero velocities, effectively functioning as relict (dead) ice patches rather than dynamically active glaciers. Area trajectories indicate that once fragmentation occurs, subsequent area loss within detached units is minimal, suggesting that fragmentation represents a threshold response to prolonged downwasting, rather than a gradual continuation of areal shrinkage. Consequently, basin-scale area change alone underestimates the extent of cryospheric degradation, as dynamic connectivity and ice flux are disproportionately reduced (Bhambri et al., 2023).

Fragmentation disproportionately affects large debris-covered glaciers, where thick supraglacial debris locally suppresses surface ablation but promotes long-term thinning and stagnation (Kulkarni et al., 2007; Bhambri et al., 2011). Our analysis reveals that the main trunk of the GL01 continues to lose area, while its fragmented counterparts remain largely stagnant, highlighting a growing spatial decoupling of ice dynamics within a single glacier system. From a process perspective, this transition marks a shift towards reduced ice flux to lower elevations, enhanced surface heterogeneity, and diminished dynamic connectivity across glacier systems (Scherler et al., 2011). These changes carry important hydrological implications, as fragmented and stagnant ice bodies contribute little to sustained meltwater buffering while enhancing short-term meltwater pulses and local instability during warm seasons (Mishra et al., 2024).

Trends in glacier length change reinforce this interpretation of dynamic weakening. Mean cumulative frontal retreat across the nine glaciers analysed amounted to 556.69 m between 1994 and 2024 (18.55 ± 1.05 m yr⁻¹). Frontal retreat alone does not capture the full extent of change in these systems. The debris-covered tongue of GL01 retreated 780.37 ± 45.80 m over 3 decades, yet its terminal velocity fell from 25.5 to 11.0 m yr⁻¹, and the lower ablation zone thinned at the rate of 1.67 m yr⁻¹ over 2000–2020. Surface lowering therefore dominates over frontal retreat, consistent with the widely documented behaviour of debris-covered Himalayan glaciers where thinning proceeds despite limited terminus retreat (Bhambri et al., 2011; Dobhal et al., 2013; Das et al., 2022; Bhambri et al., 2023).

Surface velocity measurements provide further insight into the evolving glacier dynamics of SUW. Longitudinal velocity profiles reveal strong spatial and temporal variability across glaciers of different sizes and morphometric characteristics (Das and Sharma, 2021). The GL01 exhibits the highest overall activity, with velocities ranging from ∼0.6 to 112 m yr⁻¹ with a mean of ∼23.4 m yr⁻¹, indicating the persistence of a dynamic core despite overall slowdown. GL07 shows comparably high but more constrained motion (∼0.9–93.6 m yr⁻¹, mean ∼19 m yr⁻¹) with pronounced down-glacier deceleration. GL28 displays a wide velocity range (∼0.9–117.4 m yr⁻¹) but a lower mean (∼14.4 m yr⁻¹), suggesting episodic acceleration despite long-term weakening. GL39 and GL44 exhibit markedly reduced dynamics (relative to GL01, GL07 and GL28), with mean velocities of ∼11.3 m yr⁻¹ and ∼11.9 m yr⁻¹, and minimum values < 1 m yr⁻¹, indicative of near-stagnant lower tongues. A coherent temporal signal emerges across all glaciers, with relatively high velocities during 1993–94 and 2000–01, a pronounced basin-wide minimum during 2009–11, partial recovery from 2014 to 15, and a tendency toward stabilization or decline during 2020–22, particularly at lower elevations. The basin-wide mean velocity (five major glaciers) decreased from ∼23 m yr⁻¹ (1993–94) to ∼20.6 m yr⁻¹ (2000–01), fell to a minimum of ∼8.8 m yr⁻¹ in 2009–11, and recovered partially to ∼13.6 m yr⁻¹ (2014–15) and ∼14.2 m yr⁻¹ (2020–22), a net deceleration of ∼38% across the observation period. This pattern reflects the combined influence of short-term dynamic thinning and long-term mass loss, with large glaciers retaining limited dynamic activity while smaller and debris-covered glaciers transition toward stagnation (Das et al., 2022).

Hypsometric characteristics provide an important geomorphic context for these contrasting dynamic responses. Very top-heavy and top-heavy glaciers, such as GL44 (HI = −1.86) and GL07 (HI = −1.48) retain a substantial portion of their area at higher elevations (>5000 m. a. s l), favouring sustained accumulation and explaining the persistence of dynamic cores (Deswal et al., 2023). In contrast, equidimensional glaciers such as GL09 (HI = −1.04) and GL39 (HI = −1.07) exhibit limited high-elevation ice storage and show lower, more uniform mean velocities, making them more susceptible to stagnation and fragmentation under sustained warming (Robson et al., 2018; Deswal et al., 2025). Bottom-heavy glaciers (e.g., GL01, HI = 1.25; GL28, HI = 1.30) show transitional behaviour. According to hypsometric calculations, GL01 is classified as a bottom-heavy glacier, but displays a wide variance in mean velocity values, vis-à-vis mass transfer. Interestingly, most mass-contributing tributary glaciers join GL01 at elevations below ∼4500 m. a.s.l. This highlights the role of complicated topography in modulating geometric and statistical calculations performed especially in the Himalayan region. The current analysis helps ascertain that glacier responses are shaped by hypsometry in idealized conditions, but are more actively driven by valley geometry, debris cover, climatic forcing, supra- and subglacial features, and positions of ice stream intersections along the trunk glacier, overriding hypsometric controls. However, HI remains a useful diagnostic tool, but must be interpreted in conjunction with geomorphic and climatic parameters for an accurate assessment of glacier health in the western Himalaya (Deswal et al., 2025).

Debris-covered area increased from 13.67 ± 2.26 km² in 1994 to 24.84 ± 1.15 km² in 2024, corresponding to a net rise of 11.17 ± 0.014 km² (Table 2). Debris expansion was modest prior to 2010 (+0.84 ± 0.026 km² during 1994–2000 and +0.53 ± 0.019 km² during 2000–2010), but intensified sharply thereafter, with ∼9.80 ± 0.011 km² (∼88%) of the total increase occurring during 2010–2024. This timing closely mirrors the increase in fragmentation, velocity decline, and surface heterogeneity across the basin. The rapid expansion of debris cover reflects progressive downwasting, reduced ice flux, and increased exposure of englacial debris generated from destabilized surrounding slopes. As surface velocities decline, debris becomes increasingly concentrated on low-gradient glacier tongues, reinforcing stagnation through feedback between reduced ice flow and debris accumulation (Das et al., 2022). These results demonstrate that debris cover does not confer long-term stability; instead, it masks internal ice loss while promoting surface heterogeneity. Formation of ice cliffs and supraglacial lakes, act as ablation hotspots and has the potential to form glacial lakes, which are susceptible to floods (Racoviteanu et al., 2022). Hence, studies need to move beyond simplistic assumptions of debris-related insulation when assessing glacier resilience under warming climates (Das et al., 2022). This limitation extends to glacier mapping itself: the substantial area deficit of the ICIMOD inventory relative to our outlines for approximately the same period of observation (2000/2005) (Supplementary Figure S5), is largely attributable to unmapped debris-covered ice, suggesting that inventories derived from spectral-index methods may under-represent glacier extent in heavily debris-mantled basins. In parallel, periglacial landforms such as rock glaciers, ice-rich tills, and debris cones store substantial interstitial ice and are also sensitive to sustained warming (Pandey et al., 2024). Progressive thawing can destabilize slopes, alter drainage pathways, and increase the likelihood of floods, rockfalls, and debris-flow events, particularly in valleys undergoing rapid transformation (Burn et al., 2025).

High-altitude lakes mapped between 2020 and 2024 (Figure 9) exhibit a clear expansion trend across multiple locations. Lakes were mapped according to their dam type: Bedrock dammed (BDL), Moraine dammed (MDL), and Ice dammed (IDL). Nine of the 34 mapped lakes show an increase in area between 2020 and 2024; the remaining 25 were unchanged. The most pronounced expansions were observed in L27 (BDL, +101.9%), L21 (MDL, +84.9%), L17 (IDL, +83.9%), and L16 (MDL, +83.2%), respectively. Expansion rates varied across the dam types, i.e., 4 of the 5 IDL expanded (80%), compared with 4 of 20 MDL (20%) and 1 of 9 BDL (11%), indicating that ice-dammed lakes respond most readily to current melt conditions. L16 (MDL) occupies a precarious setting (Figures 9c,d), characterized by steep surrounding slopes, a rapidly thinning frontal moraine, and enhanced meltwater supply from an actively thinning glacier snout. The rapid growth of ice-dammed lakes, together with increasing surface stagnation, indicates localized destabilisation of lower-ablation zones. Such conditions are conducive to the development of cryospheric hazards and may increase the likelihood of future glacial lake outburst flood (GLOF) events, which are increasingly documented across the region (Emmer et al., 2022; Das et al., 2024).

The observed glacier recession, surface lowering, expansion of supraglacial debris cover, and reduction in glacier flow velocities across the SUW are consistent with long-term climatic changes in the region. Climate records from 1901 to 2024 indicate a significant increase in annual mean temperature by approximately 0.006 °C yr⁻¹, with the strongest warming occurring during the winter (DJFM: +0.010 °C yr⁻¹) and post-monsoon (ON: +0.010 °C yr⁻¹) seasons. In contrast, pre-monsoon (AMJ: +0.004 °C yr⁻¹) and monsoon (JAS: +0.002 °C yr⁻¹) temperatures exhibit comparatively weak and statistically insignificant trends (Figure 8). Precipitation records reveal contrasting seasonal changes, with significant increases during the pre-monsoon (+0.354 mm yr⁻¹) and post-monsoon (+0.070 mm yr⁻¹) periods, whereas winter precipitation shows a significant decline (−0.160 mm yr⁻¹). Monsoon precipitation displays a weak and statistically insignificant decrease (−0.137 mm yr⁻¹), resulting in only a marginal increase in annual precipitation rates (+0.229 mm yr⁻¹) (Figure 7e). SUW is particularly sensitive to short- and medium-term hydroclimatic shifts in seasonality and volume of incoming precipitation. Although the basin is classified as monsoon-dominated (Figure 7), the SUW also receives a considerable amount of winter precipitation from the MLW activity. However, the combination of warmer winters and declining winter precipitation, as also observed in neighbouring basins (Chand et al., 2020; Patel et al., 2021), likely reduces effective snow accumulation. Subsequent initiation of melt with insignificant but rising pre-monsoon temperatures thereby lengthens the ablation season, enhancing glacier mass loss. Long-term climate analysis therefore provides the broader forcing context for the ongoing cryospheric changes in SUW. Continued glacier thinning, debris expansion, and dynamic slowdown signal a gradual reduction in the cryosphere’s capacity to function as a sustained freshwater reservoir, with potential consequences for seasonal runoff regulation, irrigation, hydropower generation, and water security downstream.

Our findings highlight the need for evidence-based mountain water and hazard management, which strongly aligns with SDG 6 (water) and SDG 13 (climate) and broader sustainability frameworks. Basin-scale glacier inventories integrating glacier extent, surface velocity, debris cover, thinning rates, and hypsometric characteristics provide robust datasets for planners, disaster-management authorities, and water-resource agencies. Identifying zones of stagnation, enhanced surface lowering, and supraglacial lake development can help prioritize long-term monitoring programmes and the implementation of early-warning systems.

Beyond direct cryospheric hazards, the observed changes in glacier parameters have important implications for future mountain water resources. Recent studies increasingly recognize groundwater as a critical but poorly quantified component of mountain hydrological systems, with glacial and periglacial meltwater contributing to groundwater recharge and sustaining streamflow during dry periods (van Tiel et al., 2024). As glaciers continue to lose mass and transition towards increasingly stagnant conditions, shifts in cryosphere-groundwater connectivity may alter groundwater recharge pathways, storage, and catchment-scale hydrological processes. While the relative importance of groundwater may increase as the cryosphere diminishes, the long-term decline in glacier meltwater inputs could ultimately reduce groundwater availability and modify seasonal runoff regimes. The Sara Umga case therefore demonstrates how process-based glacier research can bridge cryosphere science, mountain hydrology, sustainability science, and adaptive policy responses in climatically sensitive Himalayan zones.

6 Conclusion

This study presents the first comprehensive multi-decadal assessment of glacier change in the SUW, in the western Himalaya. By integrating multi-temporal satellite datasets, the analysis combines glacier geometry (area and length), surface velocity, elevation change (Δh), debris-cover evolution, and hypsometric characteristics to move beyond conventional area-based assessments. The strategic location of the SUW at the margin of contrasting accumulation regimes (ISM-dominated summer accumulation and MLW-dominated winter accumulation) enables direct observation of glacier responses in a region where relatively small shifts in atmospheric circulation may produce disproportionately large impacts on glacier mass balance and dynamics. This multi-parameter framework enables a process-based understanding of glacier response by explicitly linking surface morphology, ice dynamics, and climatic forcing across the region.

Between 1994 and 2024, glacierized area declined modestly from 132.92 ± 9.28 km² to 128.38 ± 2.98 km² (∼3.42%), but this limited area loss reveals substantial internal reorganization. Nearly 65% of total area loss occurred after 2010, accompanied by increasing fragmentation of large valley glaciers into smaller, disconnected units. In contrast, debris-covered areas expanded significantly from 13.67 ± 2.26 km² to 24.84 ± 1.15 km², with most of this increase occurring after 2010. Surface velocities show a clear long-term slowdown, with basin-wide mean velocities decreasing from ∼22.9 m yr⁻¹ in the early 1990s to ∼14.2 m yr⁻¹ by 2022, following a pronounced minimum (∼8.8 m yr⁻¹) in 2009–11, and near-stagnant conditions developing in debris-covered lower tongues. Elevation-change analysis reveals strongly heterogeneous and accelerated thinning, with maximum surface lowering in lower ablation zones (∼4000–4200 m a.s.l.; up to −1.67 ± 0.39 m yr⁻¹) and near-equilibrium conditions at higher elevations. Together, these trends indicate a transition from dynamically active glaciers toward systems characterized by fragmentation, reduced ice flux, and spatial decoupling between upper accumulation zones and stagnant lower tongues.

Our findings suggest that glacier change in the SUW is governed not by rapid area loss, but by progressive weakening and internal downwasting. The observed slowdown, thinning, and debris expansion indicate a declining capacity of glaciers to sustain long-term meltwater contribution, with implications for increasing seasonal variability in water availability and emerging cryospheric hazards. In climatically sensitive watersheds such as Sara-Umga, continued glacier degradation may amplify both water scarcity and hydrological extremes. Future research should prioritize ground-based long-term monitoring of glacier dynamics, debris evolution, and surface heterogeneity complemented by high-resolution remote sensing and modelling approaches to better understand glacier-climate interactions and support sustainable water-resource management and hazard mitigation in the Himalayan region.

Correction note

This article has been corrected with minor changes. These changes do not impact the scientific content of the article.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. During the preparation of the work, the authors used Quill Bot and Grammarly in order to refine the language of the manuscript. However, the content was created by authors, and authors reviewed and edited the content as required, and take full responsibility for the content of the published article.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2026.1896483/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Atanu Bhattacharya, JIS University, India

Reviewed by: Shruti Dutta, Sikkim University, India

Masoom Reza, University of Ladakh, India