Thinning-induced glacier deceleration in the Zanskar Himalayas

Trends in surface velocity provide insight into glacier response to climate change as well as local drivers of ice dynamics. The Zanskar Himalayas are heavily glacierised, and retreating glaciers pose a threat to local and regional water security. Satellite image correlation provides a tool for observing surface velocity over multiple glaciers in a remote and challenging area for field work, providing key observations for tracking changes in this important region. This study provides a comprehensive analysis of long-term (1992–2023) interannual glacier surface velocity and elevation change for 12 selected glaciers with varied morphological characteristics in the Zanskar Basin of the Ladakh Himalayas. We show a statistically significant, region-wide deceleration in glacier flow, with a mean rate of 2.43 m yr⁻¹ per decade. This slowdown is accompanied by progressively increasing surface thinning, from ∼ 0.22 m yr⁻¹ between 2000 and 2005 to ∼ 0.57 m yr⁻¹ by 2015 and 2020, with the strongest thinning and velocity reductions occurring at lower elevations. The spatial and temporal agreement between thinning and velocity decline, together with the expected reduction in driving stress associated with thinning, suggests that ice mass loss is the dominant control on the observed deceleration. Glacier-specific characteristics such as geometry, topography, debris cover, and terminus type further influenced glacier response. For example, lake-terminating glaciers exhibited local increases in ice velocity near their termini. Overall, our results demonstrate a sustained multi-decadal slowdown in glacier flow that is closely linked to ice mass loss, providing new insight into the mechanisms governing glacier dynamics in the Zanskar Himalaya under ongoing climatic change.

1 Introduction

Climate change has severely impacted glaciers across the planet (Bolch et al., 2012a; Hugonnet et al., 2021a; Immerzeel et al., 2010; The GlaMBIE Team et al., 2025) as they are highly sensitive to climate forcings and thus serve as a significant indicator of climate change (Mackintosh et al., 2017; Oerlemans, 1989). Globally, glacier mass outside of Greenland and Antarctic Ice Sheets is projected to decline by 26 ± 6 % to 41 ± 11 % by the end of the 21st century (relative to 2015), under emission scenarios that correspond to +1.5 and 4 °C warming, respectively (Rounce et al., 2023). The Himalayan-Karakoram region, which is often called the Third Pole or the Water Tower of Asia (Immerzeel et al., 2010; Viviroli et al., 2011), hosts one of the largest volumes of glacierized regions outside of Greenland and Antarctica (Wester et al., 2019). As with glaciers worldwide, Himalayan glaciers have also experienced accelerated glacier mass loss over the last few decades (Brun et al., 2017; Maurer et al., 2019; Shean et al., 2020). They serve as a source of fresh water and play an important role in the global water cycle. For example, the meltwater generated from Himalayan glaciers and snow modulates the discharge of major river systems such as the Indus, Ganges, and Brahmaputra (Pritchard, 2019). While the glacial contribution to total runoff varies across the catchments, it underpins the water security of upstream populations and the broader hydrological systems upon which over a billion people downstream depend, recharges river-fed aquifers, and contributes to global sea level rise (Azam et al., 2021; Barnett et al., 2005; Bolch, 2017; Bolch et al., 2012b; Immerzeel et al., 2010).

As the glaciers retreat and thin, ice flow velocities can reduce (Dehecq et al., 2019). Recent studies indicate that many mountain glaciers are experiencing significant deceleration (Dehecq et al., 2019; Wu et al., 2020; Zhou et al., 2021). As glaciers flow due to their weight and gravity, by the processes of internal deformation and sliding at the bed, basal hydrology can substantially complicate this relationship (Bindschadler, 1983; Minchew and Joughin, 2020). Variations in meltwater supply to the bed modify the basal water pressure, which in turn regulates basal sliding (Bindschadler, 1983). As a result, many mountain glaciers exhibit pronounced seasonal velocity cycles, with speed-ups during periods of enhanced surface meltwater input to the bed (Greene and Gardner, 2025). These short-term hydrologically driven fluctuations can obscure longer-term dynamic trends if not carefully accounted for. Such effects can be mitigated using long time series spanning multiple decades, which average out seasonal variability and highlight persistent interannual to decadal trends. Furthermore, glacier-specific, process-based understanding can also help better interpret these results. These flow processes are also influenced by other factors, such as ice temperature, glacier geometry, bed characteristics, air temperature and precipitation (Cuffey and Paterson, 2010; Iken and Bindschadler, 1986). It follows that glacier velocity estimates are a key proxy for understanding the mass balance of a glacier, where few or no direct measurements are available (Millan et al., 2022; World Glacier Monitoring Service (WGMS), 2023).

Historically, glacier velocity has been monitored through the simultaneous development of on-field and satellite remote sensing-based measurements, both of which are established standards in glaciology (Hooke et al., 1989; Stevens et al., 2023; Vincent et al., 2022). For example, Differential Global Navigation Satellite Systems (DGNSS) can be used to track the position of ground stakes on a glacier over time, from which ice surface velocity can be estimated (Azam et al., 2012; Sugiyama et al., 2013). This method is reliable, but logistically expensive and time-consuming, and is generally limited to accessible areas. Most importantly, this method returns only point-based measurements, which limits both the spatial extent and temporal coverage that is needed to characterise the evolution of a glacier system, particularly mountain glaciers, which are remote and difficult to access (Azam et al., 2014; Dematteis et al., 2021; Patel et al., 2022; Sugiyama et al., 2013; Wagnon et al., 2007). By contrast, satellite-based remote sensing methods, such as feature tracking, provide a wide spatial coverage and an easy-to-access dataset. Such methods can be used to estimate glacier-wide surface velocity at a variety of scales from regional to global, and are efficient and robust (Li et al., 1998; Scherler et al., 2008; Scherler and Strecker, 2012; Satyabala, 2016; Bhushan et al., 2017; Dehecq et al., 2019; Millan et al., 2022).

With the availability of more satellite datasets, many ice velocity studies have been carried out in the Himalayas, providing critical insights into glacier spatio-temporal patterns and evolution. Several studies have revealed heterogeneous patterns in the velocity of the Himalayan glaciers (Bhambri et al., 2011; Dehecq et al., 2019; Garg et al., 2025; Tripathi et al., 2023). Findings from these studies show that glacier velocity varies spatially and temporally, region-wise and within the same glacier, depending on factors such as elevation, slope, size, debris cover fraction, land vs lake terminating, mass budget and other local conditions such as ice thickness and basal conditions (Agarwal et al., 2023; Bhushan et al., 2018; Dehecq et al., 2019). Recent studies integrating glacier surface velocity and elevation changes in the Himalayas found substantial heterogeneity in their pattern and trends (Bhambri et al., 2023; Garg et al., 2025). Glaciers in the Garhwal Himalayas exhibited significant surface lowering, associated with a reduction in glacier surface velocity, with Gangotri glacier being one exception with an active terminus (Bhambri et al., 2023). Debris-covered glacier in the region was found to have a heterogeneous effect on ice melt. In another study, focused on glaciers in the Chenab basin in the Western Himalaya, researchers reported a significant slowdown by 54 % and 20 % in Bhut and Warwan sub-basins, respectively (Garg et al., 2025).

Ladakh is located between the Himalayas and the Karakoram, known as the Trans-Himalayan region. Due to its high elevation and position in the lee of several mountain ranges, it has a semi-arid climate, with lower precipitation than many other Himalayan regions (Archer and Fowler, 2004). Most glaciers are relatively small glaciers (< 0.75 km²) and high in elevation (Schmidt and Nüsser, 2017; Soheb et al., 2022), but due to low precipitation, they play a crucial role in water and food security (Nüsser et al., 2012; Schmidt and Nüsser, 2017). For example, in drought years, glaciers and snowmelt may become the only source of water for the region (Thayyen and Gergan, 2010). Various in-situ mass balance measurements suggest these glaciers are losing mass (Azam et al., 2025; Mehta et al., 2021; Shrivastava et al., 1999; Soheb et al., 2020). The implications for regional water security from this ice loss are governed by the timing of 'peak water'- the tipping point where annual meltwater discharge reaches its maximum before declining due to reduced glacier volume (Huss and Hock, 2018). With peak water for the region projected to occur around ∼ 2028, as reported by Huss and Hock, (2018), the current acceleration in mass loss represents a critical transition from a short-term increase in water availability to a long-term reduction of streamflow, making this region extremely vulnerable to water shortages.

Even though many glacier velocity studies are focused on the Western Himalaya (Bhambri et al., 2011; Garg et al., 2025; Satyabala, 2016; Tripathi et al., 2023), the Ladakh region remains understudied, especially regarding long-term dynamic changes in glaciers. Previous studies are restricted to short-term velocity assessment or comparison of velocity with data from two different periods with a significant gap (Bhushan et al., 2018; Rana et al., 2023; Singh Jasrotia et al., 2024). Despite the insights these studies provide into discrete periods of glacier surface velocity, gaps remain in our understanding of continuous, long-term velocity evolution and thinning. Thus, a comprehensive analysis of the glacier dynamics in the Ladakh region is still missing.

The overall aim of this study is to understand the evolution of glacier velocity and its controlling factors. In particular, we consider the following specific objectives:

1. Characterise the long-term (1992–2023) glacier surface velocity trend for selected glaciers in the Zanskar basin, Western Ladakh and evaluate the glacier surface elevation change using the existing dataset by Hugonnet et al. (2021a).

2. Investigate the factors influencing glacier velocity evolution, including changes in surface elevation, glacier morphology and extent, and assess whether glacier velocity is linked to a single climatic factor or a combination of interacting factors.

2 Study Area

The Zanskar Basin is a high-altitude, cold desert in the Ladakh region of the Western Himalayas (Fig. 1). The region can be generally classified as a cold-arid climate, as it is in the dry trans-Himalayan region, where the penetration of the Indian Summer Monsoon is very weak. The majority of the precipitation occurs during winter in the form of snowfall (Kamp et al., 2011; Lee et al., 2014).

Figure 1 (a) Study area map of Zanskar Basin in the Ladakh Himalaya (black dashed lines showing the boundary), highlighting the glaciers selected for this study (in green) labelled as DDG (Drang Drung Glacier), HG (Hagshu Glacier) and other unnamed glaciers as G3–G12. Neighbouring glaciers around the region are shown in orange and are taken from Randolph Glacier Inventory (RGI) version 7. The inset map shows the sub-regions of High Mountain Asia (HMA) in black outline. (b) shows the field photo of Hagshu Glacier, and (c) shows the field photo of Drang Drung Glacier taken during October 2023.

The glacierised area within the basin is approximately ∼ 1700 km², comprising around ∼ 1755 glaciers (Randolph Glacier Inventory V7; RGI v7.0; (RGI Consortium, 2023). Most of the glaciers in this basin are small cirque glaciers (< 0.75 km²), although several valley-type glaciers with larger areas are also present (RGI Consortium, 2023). As a significant sub-basin of the Indus River system, the Zanskar Basin plays an important role in regional water resources and is therefore a key site for understanding the impacts of climate change on future Himalayan water availability. For example, one of the major rivers in this basin is the Doda River, which originates from the Drang Drung Glacier and flows into the Zanskar River, a major tributary of the Indus River.

In this study, we have selected 12 glaciers to study the evolution of inter-annual glacier surface velocity and surface elevation change. We selected the glaciers to capture morphological and glaciological diversity in the region based on their size, elevation, degree of debris cover, terminus type, orientation, and slope. From small mountain glaciers (G12; 3 km²) to large valley glaciers (DDG; 68 km²) and a mean elevation ranging from 4800 to 5400 m (Table 1). The glaciers are given unique IDs based on their known name (DDG- Drang Drung Glacier, HG- Hagshu Glacier); the remainder are denoted G3 to G12 (shown in Fig. 1). Refer to Table 1 for further details.

Table 1 Characteristics of the 12 selected glaciers in our region. We detail the: mean elevation (m), area (km²), length of the glacier (km), mean aspect (°), mean slope (°) and whether the glacier is debris covered or not. The debris cover percentage presented here is estimated from the Supraglacial debris cover dataset v1.0 (Scherler et al., 2018a), based on linear spectral unmixing-derived fractional debris cover (FDC). All the information is from RGI v7.

3 Data and Methods

3.1 Satellite Data

We selected satellite images from the Landsat series of sensors to generate annual velocity fields from 1992 to 2023, downloaded from the United States Geological Survey (USGS) EarthExplorer: https://earthexplorer.usgs.gov/ (last access: 25 February 2025). In this study, we used the panchromatic band (B8) of the Landsat series 7 Enhanced Thematic Mapper Plus (ETM+), Landsat 8 and 9 Operational Land Imager (OLI) and the Green (B2) band from Landsat 5 Thematic Mapper (TM), which is closest to the wavelength range of the panchromatic band. Previously, various studies have used Landsat images for velocity field estimation (Altena et al., 2019; Bhambri et al., 2023; Dehecq et al., 2015; Garg et al., 2025; Nanni et al., 2023). 24 pairs of satellite images for inter-annual velocity (Supplement Table S1) from Landsat 5, 7, 8 and 9 were manually selected based on the quality of the image available (minimum cloud cover with a 20 % threshold, minimum snow cover over glaciers) to reduce noise in the results. Satellite images were selected mostly in the months of September and October, as the end of the melting season (the melt season is generally from May–September) is typically followed by minimal snow cover and less cloud cover fraction and is a consistent season that satisfies the image quality conditions mentioned above. The presence of heterogeneous snow cover over two images may reduce or enhance the textural features and induce a mismatch in the correlation algorithm, leading to erroneous results. A previous study reported the horizontal accuracy of the Landsat sensors to be < 6 m (Tucker et al., 2004), which is within the expected range of observation uncertainty; i.e., uncertainty associated with glacier surface features introduces errors at or above the observational horizontal accuracy (Bhattacharya et al., 2016).

3.2 Glacier Velocity Estimation

3.2.1 Image Correlation

We used a pairwise image correlation algorithm to generate all the annual velocity fields from the satellite images. Horizontal displacement between a pair of images was estimated using ENVI's plugin software package COSI-Corr (freely available at http://www.tectonics.caltech.edu/slip_history/spot_coseis/, last access: 24 October 2024) (Leprince et al., 2007b). This software utilises a cross-correlation feature tracking algorithm in the frequency domain. The algorithm performs sub-pixel correlation of remote sensing datasets to capture all kinds of displacement on the Earth's surface (Leprince et al., 2007a). The phase correlation in the frequency domain was used to estimate the relative shift between the time gaps in the image pairs. Selecting an appropriate window size for feature tracking is often challenging, as it involves balancing noise suppression with spatial details. After multiple iterations and careful visual inspection of the results- particularly assessing the presence of noise, we determined that using an initial window size of 64×64 pixels and a final window of 32×32 pixels provided the most reliable outcomes, which was also found to be used by previous studies (Kraaijenbrink et al., 2016; Nanni et al., 2023; Bhushan et al., 2018). A larger initial search window is used to find a rough match between the pixel displacement with less noise and finally refined with a smaller window size. A step size (refers to the spatial resolution or the distance between the centres of the image patches) of 4×4 pixels for all except for Landsat 5 (30 m spatial resolution), where a 2×2 pixel step size was set, resulting in a resultant output of 60 m spatial resolution velocity fields. Step sizes smaller than the search window size may introduce redundancy in measurement (Leprince et al., 2007a, b), thereby improving the accuracy of the velocity estimates compared to steps that are of a similar magnitude to the search window (Fahnestock et al., 2016). The outputs from the COSI-Corr algorithm consist of 3 images: N–S (North–South) displacement component, E–W (East–West) displacement component and SNR (Signal to Noise ratio), which assesses the quality of correlation for each pixel (Leprince et al., 2007a).

3.2.2 Filtering Process

To minimize the effect of noise, which impacts the accuracy of the correlation outputs, we applied basic filtering to the N–S and E–W surface displacement components. To remove unwanted correlation results, a threshold of 0.9 was set for the Signal-to-Noise (SNR) ratio (Scherler et al., 2008). The Robustness iteration, which is a quality control loop that ensures the derived displacement or velocity vectors are consistent and robust against errors or noise, was set to 4. First, an initial displacement is calculated from image matching, followed by outlier removal. This process is repeated for the number of iterations defined. 2–4 iterations were found to be working well in most cases of glacier velocity estimation (Leprince et al., 2007a). As the three outputs are generated (N–S, E–W displacement and SNR) we used a replace/discard tool in the COSI-Corr algorithm, which removes unwanted pixel values if the velocities defined at these pixels exceeds a threshold. Here, a value of 200 m was set based on existing literature and the maximum glacier velocity in the Western Himalayan region, which is estimated to be approximately ranging from ∼ 50–90 m yr⁻¹, with a few exceptions in our case (Bhushan et al., 2018).

Finally, the resultant displacement was calculated using Eq. (1), based on the N–S and E–W velocity components for each image pair, and the final velocity was evaluated using Eq. (2).

$$\begin{array}{}\text{(1)}& {\displaystyle }{D}_{xy}=\sqrt{(D_{\mathrm{NS}}{)}^{\mathrm{2}}+(D_{\mathrm{EW}}{)}^{\mathrm{2}},}\text{(2)}& {\displaystyle }V={\displaystyle \frac{D_{xy}}{t}},\end{array}$$

Here, D_xy is the horizontal resultant displacement (meters), V is the final surface velocity (m yr⁻¹), and t is the time (years) between the two image acquisitions.

The final velocity maps are produced at 60 m resolution. The same parameter values were maintained throughout all the image pair processing to ensure the velocity outputs are consistent. Finally, we used a 3-pixel × 3-pixels median filter on the velocity maps to remove any outliers (noise) within the maps and to smooth the velocity outputs with an acceptable loss of details. Due to poor correlation in some of the pair-wise correlation processes, data gaps exist in the velocity outputs from 2003–2007 and 2011–2012.

3.2.3 Uncertainty Estimation

The uncertainty in the velocity field obtained from remote sensing could be influenced by several factors, including satellite scene characteristics (e.g., cloud cover over the region of interest, shadow, etc.), surface heterogeneity (difference in snow cover), image co-registration error, and performance of the feature-tracking algorithm. To minimise the error due to the first two factors, we selected images with a minimal cloud cover fraction < 20 %, and snow cover over the glacier surface by manually checking the scenes (September–October). Due to a lack of in situ data, it was not possible to validate our results with GNSS measurements.

Ideally, the ice-free ground (off-glacier) should have a zero vertical or horizontal displacement, which is an indicator of no errors produced in the processing. However, we observed that due to limitations in the algorithm, these ice-free areas exhibited certain displacements. For example, it may arise due to a change in surface characteristics between the two images, such as snow cover change, varying sun angles and shadows, and also due to a horizontal shift in the two images. Hence, we used the displacement in ice-free areas as a measure of uncertainty in our processing. The uncertainty was estimated based on the following method (Koblet et al., 2010; Zhang et al., 2020).

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{off}=\sqrt{{\mathrm{SE}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{X}}^{\mathrm{2}}}},\text{(4)}& {\displaystyle }\mathrm{SE}={\displaystyle \frac{\mathrm{SD}}{\sqrt{N_{\mathrm{off}}}}},\end{array}$$

Here, σ_off is the uncertainty in the velocity (m yr⁻¹), $\stackrel{\mathrm{‾}}{X}$ which is the mean of the velocity of pixels extracted from the off-glacier area, where N_off= 100 is the number of off-glacier pixel measurements. SE is the Standard Error of the mean velocity taken across the off-glacier points (m yr⁻¹), and SD is the standard deviation (m yr⁻¹). This method was applied to all the velocity datasets generated. The uncertainties range from 0.85 to 5.49 m yr⁻¹, which are of similar magnitude to previous studies (Bhushan et al., 2018; Das and Sharma, 2021; Shukla and Garg, 2020). All the pair-wise uncertainty is listed under Table S2. To show the reliability of the data generated, we compared it with the ITS_LIVE velocity dataset, as shown in Fig. S2 in the Supplement.

3.3 Glacier Surface Elevation Change

We utilised the open-access available global glacier elevation change datasets from Hugonnet et al. (2021b; accessed from https://doi.org/10.6096/13). The dataset is available for the period of 2000 and 2019, with the mean glacier surface elevation change rate data generated using ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) stereo pairs calculated across four different timeframes: 2000–2005, 2005–2010, 2010–2015 and 2015–2019, at a spatial resolution of 100 m.

We analyzed the glacier surface elevation change rate dataset using a widely used method of binning the data into multiple elevation bins (dh = 50 m). Bin-wise elevation change analysis is an effective approach for capturing the spatial heterogeneity of glacier responses to climate forcing, revealing how processes operating at different elevations shape overall surface elevation change and dynamics. Aggregating the data into elevation bins acts as a form of suppressing spatial noise while preserving the dominant elevation-dependent signals for better interpretations. Bin-wise analysis reveals where the most pronounced thinning or thickening occurs, often highlighting elevation-dependent feedback in melt, accumulation, and ice flow. For example, as surface elevation lowers and ice thins, the gravitational driving force decreases, leading to a progressive slowdown of ice flow (Cuffey and Paterson, 2010). However, there could be some exceptions, for example, differential surface thinning rate, leading to increased surface slope, may increase the driving stress or some complex feedback mechanisms from proglacial lake-ice interaction. Mathematically, the driving stress τ(x) is expressed as,

$$\begin{array}{}\text{(5)}& \mathit{\tau }\left(x\right)=\mathit{\rho }gH\left(x\right){\displaystyle \frac{\partial S}{\partial x}}\left(x\right),\end{array}$$

In this equation, τ(x) is the driving stress, ρ is the density of ice, g is the acceleration due to gravity, H(x) is ice thickness, and S(x) refers to the ice surface at position x and $\frac{\partial S}{\partial x}\left(x\right)$ represents the surface slope of the glacier at position x along a given flow line.

By contrast, sliding at the bed of the glacier can also play a significant role in ice flow, but sliding is poorly constrained due to a lack of observations. It depends on bed roughness, thermal regime of the glacier and subglacial hydrology (Bindschadler, 1983; Weertman, 1957).

3.4 Data analysis and supporting datasets

The central flowlines of the glaciers were manually delineated from high-resolution satellite data for the analysis of velocity profiles. For all other analyses, including glacier elevation changes vs glacier velocity, the datasets were resampled to the spatial resolution of elevation change data (100 m), and a reference DEM (Digital Elevation Model) was used to make all the elevation-wise analyses, such as binning the data into different elevation zones. Additionally, we used the Hypsometric Index to calculate and classify the glaciers as very bottom-heavy (HI > 1.5), bottom-heavy geometries (1.2 < HI ≤ 1.5), very top-heavy (HI < −1.5), top-heavy (−1.5 < HI < −1.2), and equidimensional glaciers (−1.2 < HI < 1.2) following Jiskoot et al. (2009). While top-heavy glaciers concentrate their mass at higher elevations within expansive accumulation zones, bottom-heavy glaciers feature a disproportionately large area at lower elevations. The glacier boundary from RGI Version 7 (Randolph Glacier Inventory V7) was used for the glacier extents and for calculating all the geospatial statistics. A T-statistic test that compares the means of two groups while considering how much the data “spreads” (variance) was used to calculate the p-value in our data analysis (Kim, 2015). The p-value, derived from this statistic, represents the probability of observing results as extreme as the sample, assuming the null hypothesis. A p-value below 0.05 generally indicates statistically significant results (rejecting the null hypothesis).

For analysing the climate trend (temperature and precipitation) of the region, we utilised the ERA5-Land Reanalysis climate dataset. ERA5-Land has a spatial resolution of ∼ 9 km, and temporally it provides an hourly dataset from 1950 to the present (Copernicus Climate Change Service, 2019). Data from the nearest grid to the study location were extracted for further analysis. The uncertainty was estimated based on the 1σ (standard deviation) of the data.

4 Results

4.1 Glacier Velocity Trends

We analysed interannual glacier surface velocity changes across 12 glaciers using the satellite-derived datasets described in section 3, with some data gaps. A spatial map of glacier velocity for the year 2020–2021 is shown in Fig. 2. Overall, the glacier flow speeds were generally greater in larger glaciers than in smaller ones. We observed a typical velocity distribution, where ice flow increased from the margins toward the centreline and decreased from the accumulation area towards the terminus. In addition, all glaciers demonstrated a very similar pattern of low glacier velocity near the terminus (∼ < 15 m yr⁻¹). In general, the velocity gradually increased with the surface slope, at times by as much as 50 % around sudden changes in slope, also coincident with crevassing.

Figure 2 Glacier velocity in m yr⁻¹ from 2020–2021 for the selected glaciers (with extents outlined in black), subdivided into three different groups. The left vertical panel (a) shows all the glaciers investigated with insets on the right outlined by yellow boxes; (b) velocity for glaciers DDG, HG, G3, G4 and G5; (c) velocity for glaciers G6, G7 and G8; (d) velocity for glaciers G10 and G11.

The results reveal temporal and spatial heterogeneities in ice flow. That is, while some glaciers exhibited a marked deceleration over time, others maintained a relatively constant velocity or showed inter-annual accelerations during specific periods. Despite interannual velocity variations, the mean flow speed over the central flowline across all the glaciers varied from 31.1±8.5 m yr⁻¹ in 1992 to 26.2±1.0 m yr⁻¹ in 2023 (∼ −16 %), with a minimum mean velocity of 24.5±5.7 m yr⁻¹ during 2010, and a maximum of 46±5.1 m yr⁻¹ from 2000 and 2002. The overall region-wide analysis for all the selected glaciers showed a statistically significant decrease in mean glacier flow speed over the study period (p < 0.01, estimated using the t-statistic test; Kim, 2015), at a rate of −2.43 m yr⁻¹ per decade (Fig. 3). The fastest flowing glacier is the Drang Drung Glacier (DDG) with a mean velocity of ∼ 57.9 m yr⁻¹ along the flowline during 1992 and 2023, followed by G10 with a mean velocity of 38.5 m yr⁻¹. Interestingly, both glaciers are lake-terminating (discussed more in Sect. 5.3). By contrast, the slowest-flowing glacier is G12 with a mean velocity of ∼ 9.6 m yr⁻¹, which could be related to the glacier geometry and ice thickness.

Figure 3 a) Mean velocity (m yr⁻¹) for all the glaciers. (b–m) Mean velocities (m yr⁻¹) for individual glaciers from DDG to G12. Each panel shows the interquartile range (shaded orange and blue box) with the median denoted by the black horizontal line, lines and whiskers denote the maximum value (upper) and minimum value (lower), excluding the outliers. The trend is reported in the upper right corner of each panel (which considers the data gap).

Although the overall regional slowdown in glacier velocity is statistically significant, not all individual glaciers exhibit significant trends. In particular, DDG, G5, and G10 show statistically non-significant velocity trends (Fig. 3).

We considered three decades from our study period (1992–2000, 2001–2010, 2011–2023) to understand decadal flow pattern changes. Despite the data gaps between 2001 and 2010, we can discern a broad understanding of changes in the flow trend over this 31-year period. From 1992 to 2000, all glaciers exhibited an increase in velocity. In contrast, from 2001 to 2010 and 2011 to 2023, the glaciers generally show reduced flow speeds. To analyze the relationship between ice velocity trend and elevation, we categorized the glaciers into two elevational zones based on the regional mean elevation: a lower-elevation zone (< 5000 m) and an upper-elevation zone (> 5000 m). It is important to note that this 5000 m threshold was chosen to represent the average elevational distribution of the glaciers in our study area and is not intended to represent the Equilibrium Line Altitude (ELA). The mean velocities within both elevation groups showed a declining trend at rates of −4.6 and −1.7 m yr⁻¹ per decade, respectively. This suggests the glacier velocity declined faster near the lower elevation zone (< 5000 m) compared to the higher elevation zone (Fig. 4), which is consistent with other studies (Wu et al., 2020).

Figure 4 Box plots for the mean velocity (m yr⁻¹) for the region of glaciers with elevation (a) below 5000 m; and (b) above 5000 m. In each panel, the median velocity trend is denoted by the black horizontal line, the interquartile range is shown in the coloured box, and the lines and whiskers denote the maximum value (upper) and minimum value (lower), excluding the outliers. Outliers are not included in the plot.

Velocity profiles along the central flowline are shown in Fig. 5. The patterns correspond with glacier hypsometry. Glaciers like G4, G8, and G11 exhibited a bottom-heavy geometry, while HG and G3 are very bottom-heavy (HI > 1.5), with a larger share of area at lower elevations. Such geometries can enhance driving stress in the ablation zone due to thinning-induced surface steepening and increased ice mass concentrated at lower elevations (Bhushan et al., 2018; Sam et al., 2018).

Figure 5 (a–l) Scatter plots showing the glacier velocity profile along the central flowline for all the glaciers from 1992–2023. The x-axis represents the distance along the flowline (km) from the start of the terminus (marked by the origin). The different colours represent data points from different years. The black thick line refers to the mean velocity across the time period.

Elevation bin-wise velocity analysis also revealed inter-glacier variability associated with glacier geometry changes and surface elevation lowering. We evaluated elevation bin-wise (bin size = 50 m) velocity evolution across different periods (1999–2000, 2009–2010, 2014–2015 and 2019–2020), finding velocity trends that coincided with specific elevation bands, which sheds light on the underlying dynamics influencing the trends. We found that most glaciers in our study exhibited similar patterns of glacier slowdown across the glacier, with the greatest reductions occurring near the ablation zone (Fig. S3). For instance, in DDG (panel b), the peak velocity during the 2000 and 2005 period exceeds 40 m yr⁻¹, while in the most recent period between 2015 and 2020), velocities have significantly reduced across the entire elevation range. In other cases (e.g., panels c, g, l), the velocity profiles exhibited a marked shift in both the peak magnitude and elevation position, suggesting changes in flow dynamics, which could be potentially driven by mass loss and evolving ice geometry. Some glaciers, such as G9 (panel f) and HG (panel h), exhibited more complex or less consistent trends, possibly reflecting the interplay of local topography, debris cover, and variations in surface elevation change. Several glaciers, such as HG, G4, and G5, showed velocity increases at their lower elevations in 2019–2020.

Table 2 Overall mean velocity trend (m yr⁻¹) for individual glaciers over the period 1992–2023. Uncertainty is calculated as one standard deviation (m yr⁻¹) and the p-value denotes a statistically significant trend, as calculated using the t statistic-test.

4.2 Glacier Surface Elevation Change

Glacier surface elevation change revealed a clear altitudinal dependency, with consistent thinning patterns within each elevation bin over the period 2000–2020. Across the 12 glaciers analyzed, thinning generally increased in the lower elevation bins, particularly below ∼ 4700 m, indicating stronger ablation and dynamic thinning towards glacier termini (Fig. 6). This pattern persisted across the period (2000–2020), with the lowest elevation bands consistently exhibiting the greatest thinning rate (dh $/$ dt) of ∼ −1.3 m yr⁻¹. While high-elevation regions (> 5000 m) also experienced surface lowering, the rates were relatively lower, likely reflecting reduced melt and less dynamic thinning at those elevations.

Figure 6 (a–l) Glacier hypsometric curve of surface elevation change. The left y-axis shows glacier area (km²), represented as grey histograms and the right y-axis shows the glacier surface elevation change rate, dh $/$ dt (m yr⁻¹) from Hugonnet et al. (2021a), data for four different periods, 2000–2005, 2005–2010, and 2015–2019 represented in red, blue, green and purple lines, respectively. The colored bands represent the 1σ (standard deviation).

Figure 7 Glacier-wide mean elevation change (dh $/$ dt) in m yr⁻¹ for different periods. Blue denotes thickening and red denotes thinning.

The overall mean glacier surface elevation change rates for all the glaciers in the study region show a clear acceleration of thinning, from −0.22 m yr⁻¹ (2000–2005) to −0.57 m yr⁻¹ (2015–2020). This pattern aligns with glacier hypsometry, as bottom-heavy glaciers (e.g., HG, G3, G4, G8, G11) experienced greater thinning, with HG showing especially strong losses in the most recent period (−0.72 m yr⁻¹; Fig. 7). The largest glacier – DDG, an equidimensional glacier (where glacier area is more evenly distributed across its elevation range), experienced a mean glacier-wide thinning rate of −0.06 m yr⁻¹ during 2000 and 2005, to −0.63 m yr⁻¹ during 2015 and 2020. Few glaciers (HG, G3 and G6) showed a mean positive rate during 2000 and 2005, which could be linked to the higher accumulation rate in higher elevations (Fig. 8).

4.3 Non-climatic factors

We used a correlation matrix to understand the interplay between glacier surface velocity and different non-climatic factors such as mean elevation, area, mean slope, mean aspect, and length (Fig. 9). The matrix indicates that glacier length and area have the strongest positive correlations with mean velocity, with a correlation coefficient of 0.75 and 0.67, respectively, which could be related to how large and thick the ice is. In contrast, slope ($r=-\mathrm{0.37}$) and aspect ($r=-\mathrm{0.36}$) shows a moderate negative correlation with velocity. While steeper slopes generally increase gravitational driving stress, the negative correlation observed here likely reflects the inverse relationship between slope and ice thickness; in this study area, steeper glaciers could be thinner, and the reduction in thickness exerts a stronger limiting influence on velocity than the slope gradient provides. The mean elevation shows almost no correlation (r=0.01) with surface velocity.

Figure 8 Glacier surface elevation change (Hugonnet et al., 2021a). (a–b) show the elevation change map for different periods (2000–2005 and 2015–2020). The colour bar represents the 2nd–98th percentile of elevation change. The inset in panels (a) and (b) show the tongue of the Drang Drung glacier and its neighbouring glaciers.

Figure 9 Correlation matrix between mean surface velocity and non-climatic variables such as Mean elevation (m), Area (km²), Slope (in °), Aspect (in °) and length (km).

4.4 Climatic trends

We analyzed climate data over the region from the ERA5-Land dataset, calculating trends over the period 1990–2025. The overall mean annual temperature increased by 0.02±0.01 °C yr⁻¹ over the period, with increases in both summer (May–September; 0.01±0.01 °C yr⁻¹) and winter (October–April; 0.03±0.01 °C yr⁻¹) temperatures (Fig. 10a). Summer precipitation (MJJAS) remains consistent with no significant change. In contrast, the winter precipitation shows a declining trend (Fig. 10b).

Figure 10 (a) Trend of summer temperature (MJJAS), winter temperature (ONDJFMA), and annual mean temperature from ERA5 Land; (b) Summer precipitation and winter precipitation (1990–2025) from ERA5 Land. Each data point for temperature is the mean, whereas for precipitation, it is cumulative.

5 Discussion

5.1 Trends in Surface Velocity and Elevation Change

Our analysis shows an overall consistent slowdown in glacier velocity over recent decades. This trend aligns well with other regional and glacier-specific observations of glacier deceleration, which are linked to ongoing thinning and mass loss. For instance, Dehecq et al. (2019) reported a significant decadal slowdown (2000–2017) in glacier flow across High Mountain Asia. Especially, the Lahaul-Spiti region, within our study area, exhibits an average velocity trend of −4.6 ± 0.6 m yr⁻¹ per decade. Similar patterns have been reported for Parkachik Glacier in the Suru Basin, adjacent to our study area, where velocity decreased from 45.2 ± 1.8 m yr⁻¹ in 1999 and 2000 to 32.28 ± 0.80 m yr⁻¹ in 2020 and 2021, primarily linked to mass wastage, debris-cover increase, and reduced accumulation (Rana et al., 2023). Another study calculated the long-term change (1993–2018) of Parkachik Glacier, showing evidence of velocity slowdown by ∼ 35 % (10.5±5.7 m yr⁻¹) due to increased debris cover near the terminus and mass loss (Garg et al., 2022). A recent study focusing on the Chenab Basin showed glaciers slowing down by 54 % and 20 % between 1993 and 2021 in the Bhut and Warwan subbasins, respectively, which was primarily linked to increased debris-cover over glaciers (Garg et al., 2025). Another recent study by Bhattacharjee et al. (2025) reported a velocity of 0.26 ± 0.02 m d${}^{-\mathrm{1},}$ which corresponds to 94.97 ± 7.31 m yr⁻¹ for DDG close to the ELA (Equilibrium Line Altitude) and 0.04 ± 0.003 m d⁻¹, equivalent to 14.6 ± 1.1 m yr⁻¹ for G3 near the terminus, which aligns very well with our results. Although previous studies have reported a slowdown of glaciers in the Ladakh region, a contrasting increase in glacier velocity was found in the Drang Drung Glacier. For example, the results from the study by Singh Jasrotia et al. (2024) show velocity increased from 71±6.1 m yr⁻¹ in 1999 and 2000 to 140±7.4 m yr⁻¹ in 2019 and 2020 (by ∼ 50 %), which contrasts with our findings, and is likely an artefact of their consideration of two timeframes only. Overall, our declining glacier velocity trends align with patterns observed for nearby glaciers in other studies. However, our results contrast with the recorded flow acceleration of some Karakoram glaciers during 2000/2002 and 2009/2010, as reported by Heid and Kääb (2012), which was due to positive mass balance, and also glacier surging in this region.

The slowdown we observed could be primarily explained by sustained glacier thinning, resulting in reduced driving stress, especially at lower elevations in the ablation zones, often exceeding several meters per year in some cases. Detailed investigations of surface elevation change revealed strong glacier thinning across all the glaciers, with higher rates observed from 2010 to 2020. HG and G8 exhibited the highest negative elevation change (0.72 m yr⁻¹) between 2015 and 2020, despite them being debris-covered glaciers. This pattern matches the results reported by Bhushan et al. (2018), where debris-covered glaciers exhibited similar elevation change compared to clean ice glaciers. This could be likely due to melt-enhancing features on debris-covered ice, such as ice cliffs and supraglacial ponds that conduct a significant amount of surrounding heat into the ice, resulting in ablation (Brun et al., 2016; Pellicciotti et al., 2015; Reznichenko et al., 2010). A geodetic mass balance estimated by Mandal et al. (2024) also reports an elevation change rate of $-\mathrm{0.44}\pm \mathrm{0.09}$ m yr⁻¹ (2000–2017) for Western Ladakh glaciers that covers Zanskar Basin as well, which aligns well with Hugonnet et al. (2021a) trends.

Several glaciers in our domain (DDG, HG, G4 and G5) showed higher velocities in lower elevations in recent years, accompanied by lower thinning rates near those elevations (Fig. 6). For G4 and G5 (top-heavy glaciers) and DDG (equidimensional hypsometry), sustained ice flux from upstream likely contributes to maintaining flow toward lower elevations. However, the dynamic behaviour of these lower ablation zones cannot be explained by geometry alone. Enhanced surface melt in recent years likely increases meltwater inputs to the glacier bed through moulins or crevasses, promoting periods of elevated basal water pressure and reducing effective pressure at the ice-bed interface, which can lead to higher basal sliding. HG, a heavily debris-covered glacier with a bottom-heavy hypsometry is likely more sensitive to such hydrological controls near the ablation zone. While debris can modulate surface melt spatially, meltwater routed through crevasses, ice cliffs, and supraglacial ponds can still reach the bed and locally enhance sliding. Bhushan et al. (2018) also reported a dynamically active trunk of DDG (> 20 m yr⁻¹ in 2013–2014), consistent with our mean flowline velocities. This highlights the complex interplay between glacier geometry, debris, subglacial hydrology evolution, and dynamics.

Overall, our findings reinforce the broader picture of sustained glacier slowdown across the Zanskar region, which is closely tied to glacier mass loss due to thinning, and is consistent with the patterns observed in other studies and existing theory. The climate data analysis further supports this link, showing increasing summer temperatures and declining winter precipitation (snowfall) over recent decades (Fig. 10), which is likely to intensify in the future under continued warming.

Figure 11 Velocity near the terminus of the lake terminating glaciers. (a) DDG; and (b) G10 for the year 2020–2021. Highlighting the high velocity near the terminus for G10 and DDG. “DC” represents the Debris Cover area along the glacier tongue. The basemap is a Sentinel 2 satellite image from September 2021.

5.2 Velocity trend of lake-terminating glaciers

Glacier velocities across the study region show clear differences between lake-terminating and land-terminating glaciers, with the former generally maintaining higher surface flow speeds – often up to twice that of their land-terminating counterpart. For instance, mean velocity along the central flowline for DDG and G10 corresponds to 57.9 and 38.5 m yr⁻¹, respectively, which is significantly higher than other glaciers (Fig. 11). While glacier size and thickness have controls on driving stress and could partly explain these higher velocities, the dynamic influence of proglacial lakes likely provides an additional control on them.

The presence of a pro-glacial lake has been previously linked with enhanced glacier mass loss (King et al., 2019; Minowa et al., 2023; Zhang et al., 2023). That is, a lake-terminating glacier loses mass at its terminus by surface melt and frontal ablation, which includes mechanical calving and subaqueous melt (Carrivick and Tweed, 2013), both of which drive dynamic thinning and retreat. This is evident from the elevation change analysis of both DDG and G10 (Fig. 6). For DDG, retreat rates of 21.1 m year⁻¹, and a total retreat of 925 m have been observed between 1971 and 2017, with the pro-glacial lake developing around 2014 (Rashid and Majeed, 2018). This thinning and retreat could be related to calving, which was evidenced in 2023 (e.g. see Fig. S4). The pro-glacial lake near the G10 terminus was formed well before 1975 (Govindha Raj, 2010). In contrast to DDG, calving at G10 is not apparent from satellite observations, and it is less clear what mechanisms underlie the thinning of this glacier. Such lake expansion could alter the local force balance near the termini, promoting frontal acceleration (Pronk et al., 2021). A recent study on Bhutan lake-terminating glaciers found, despite high thinning rates, some glaciers showed acceleration near the terminus, which is potentially due to increased meltwater flux (Hyde et al., 2025). A similar conclusion was drawn by Wu et al. (2020) for a lake-terminating glacier in Kangri Karpo Mountains in South-Eastern Tibet.

Numerical modelling experiment on an alpine glacier has also explained that the presence of a pro-glacial lake enhances the glacier retreat by 4 times and induces ice acceleration by 8 times higher as compared to land-terminating glaciers forced with the same climatic parameters (Sutherland et al., 2020). The contrasting behaviour of lake-terminating glaciers may therefore be explained in two ways: (a) less buttressing and altered local force balance at the ice-lake margin, and (b) enhanced surface lowering near the terminus. The latter, in particular, influences the flow characteristics by regulating the glacier's dynamic thinning phenomenon, as also discussed and found in the majority of lake-terminating Himalayan glacier studies (King et al., 2019; Pronk et al., 2021). Glacier geometry is equally important in this regard, as longer and thicker glaciers inherently sustain greater driving stresses, enabling higher ice flow velocities irrespective of terminus conditions. However, disentangling the relative influence of these geometric controls from lake-driven dynamic processes remains beyond the scope of the present dataset, as the observations do not directly constrain either factor. This limitation warrants further investigation.

5.3 Limitations

While the study provides some critical insights into the evolution and control of interannual glacier velocity for selected glaciers in the Zanskar region of Ladakh, a few limitations need to be acknowledged. First, issues in velocity retrieval led to data gaps between 2002 and 2008, primarily due to the limited availability or poor quality of optical satellite imagery with persistent cloud cover and low surface contrast hindering reliable feature tracking. These data gaps in time series analysis may miss specific events during the acceleration or slowdown phase, especially during periods of high inter-annual variability, such as high-melting years, which can lead to biases in the trend estimates.

Second, the analysis focuses on interannual to interdecadal variability and does not resolve short-term or seasonal velocity fluctuations driven by meltwater input and subglacial hydrological processes. In addition, erroneous velocity due to the presence of cloud cover, shadows, or snow may occur. For example, spatially, velocity retrieval was found to be less reliable in accumulation zones and small tributaries for some years, where low surface texture and seasonal snow cover lead to noisy or non-physical velocity vectors. In contrast, the main glacier trunks show internally coherent and consistently downstream flow, supporting confidence in the primary velocity patterns while highlighting regions of reduced reliability (Troilo et al., 2024). The estimated uncertainties in our results are well within limits and less than the magnitudes of velocity observed in the terminus of DDG and G10.

Finally, due to the lack of reliable in situ meteorological observations, ERA5-Land data were not bias-corrected. While the use of raw ERA5-Land data may influence the magnitude of derived climatic trends (e.g. Satyapragyan et al., 2026), it does not affect the primary objective of this study, which focuses on the temporal evolution of glacier velocity.

6 Conclusions

We set out to understand the long-term evolution of glacier velocity and its potential drivers in the Zanskar Basin of the Ladakh Himalayas. By analyzing interannual glacier velocity changes from 1992 to 2023 alongside glacier surface elevation change trends over four distinct periods between 2000 and 2019, of 12 selected glaciers, we provide a detailed picture of how glacier dynamics have evolved over time.

The results highlight an overall trend of glacier slowdown across all the glaciers, accompanied by prolonged glacier thinning, highlighting the dynamic response of the glaciers to sustained mass loss. Surface thinning has been more pronounced at lower elevations (< 5000 m), while patterns of velocity change vary across glaciers and time periods.

While climate-driven thinning is the dominant control on glacier deceleration, glacier geometry, slope, debris cover, subglacial conditions, and terminus type significantly modulate the spatial and temporal patterns of ice flow. Though lake-terminating glaciers exhibited persistently higher velocities compared to land-terminating ones, consistent with previous findings across the Himalayas, the specific mechanisms driving this difference cannot be conclusively resolved here, given current data constraints. Further, climate records indicate rising summer temperatures and declining winter snowfall, trends that are expected to intensify, implying continued reductions in accumulation, enhanced ablation, and accelerated mass loss in the future. Our findings contribute to a broader understanding of long-term spatial and temporal variability in glacier dynamics in the region influenced by the complex interplay between glacier thinning and other non-climatic factors.