Glacier mass balance and its climatic and nonclimatic drivers in the Ladakh region during 2000–2021 from remote sensing data

Abstract

This study investigates the geodetic mass balance of nearly all glaciers in the Ladakh region, which are crucial for local water security. Utilizing multiple digital elevation models from 2000 and 2021, we estimate glacier mass balances. Climatic drivers of glacier mass balances are explored using ERA5-Land reanalysis data, evaluated by in situ climate data. The study also examines the role of nonclimatic (morphological) variables on glacier mass balances. Results indicate Ladakh glaciers experienced negative mass balances during 2000–2021, with significant spatial variability. Western Ladakh glaciers lost slightly higher mass (−0.35 ± 0.07 to −0.37 ± 0.07 m w.e. a⁻¹) than eastern Ladakh glaciers (−0.21 ± 0.07 to −0.33 ± 0.05 m w.e. a⁻¹). While warming is the main driver of widespread mass loss in Ladakh, the spatial variability in mass loss is attributed to changes in regional precipitation and glacier morphological settings. Eastern Ladakh glaciers, being smaller and at higher elevations, experience lower mass loss, whereas western Ladakh glaciers, larger and at lower elevations, are more susceptible to the impact of temperature, resulting in higher mass loss. The study underscores the potentially greater vulnerability of western Ladakh glaciers to a warming climate compared to their eastern counterparts.

1. Introduction

The impact of global climate change on high mountain regions has prompted many scientists to investigate glaciers and their associated changes (You and others, 2020; Azam and others, 2021). The Himalaya-Karakoram (HK) region hosts more than 40 000 glaciers that exert an important control on regional water resources (Lutz and others, 2014; Vishwakarma and others, 2022). The HK region has large variability in topographic settings across its east-to-west stretch that gives diverse climates (Bookhagen and Burbank, 2006; Maussion and others, 2014) and hydrological regimes (Thayyen and Gergan, 2010).

The Ladakh and surrounding areas in the western Himalaya stand out as the coldest arid landscapes in the entire HK region. Also, they are home to many glaciers of varying sizes (Schmidt and Nüsser, 2017; Soheb and others, 2022). Although human settlements are sparsely distributed across Ladakh, they heavily rely on snow and glacier meltwater for their socio-economic needs (Nüsser and others, 2012, 2019; Mukherji and others, 2019). The Ladakh region receives very little precipitation compared to other regions (Archer and Fowler, 2004). Therefore, glacier meltwater is an essential lifeline for the region, especially when water availability or streamflow is very low (early summer or spring) and during drought years (Thayyen and Gergan, 2010; Pritchard, 2019; Balasubramanian and others, 2022). Any changes in regional climate and glacier meltwater will directly impact the region's livelihood (Schmidt and Nüsser, 2017). Despite the paramount importance of glaciers in the region, a comprehensive understanding of spatiotemporal changes in their mass balance and the factors regulating them remains elusive.

Glacier area changes and retreats in Ladakh and neighbouring areas have been relatively well documented. Schmidt and Nüsser (2017) noted a 10–20% decrease in glacierized areas across different mountain ranges in the Ladakh region over the last four decades (1969–2016). Shukla and others (2020) reported a 6% shrinkage in the glacier area of the Suru basin from 1971 to 2017, with individual glaciers experiencing up to ~45%. While glacier area change is a delayed signal of climate warming (Bhambri and others, 2011; Pandey and Venkataraman, 2013), glacier mass balance measurements are usually preferred for understanding the linkages between climate change and glacier response (Armstrong, 2010). Moreover, accurate measurement of glacier area or length changes is challenging for debris-covered glaciers, which limits the confidence in applying automated methods across large spatial scales (Bhardwaj and others, 2014; Racoviteanu and others, 2022). Glacier mass balance is considered an undelayed and contemporary indicator of atmospheric change and is classified as one of the seven headline climate monitoring indicators by the World Meteorological Organization (WMO; Trewin and others, 2021). So far, in situ mass balance observations are available for only three small to medium-sized glaciers across the entire Ladakh region, spanning a maximum temporal coverage of only five years.

Soheb and others (2020) reported a glaciological mass balance of −0.39 ± 0.38 m w.e. a⁻¹ for the Stok Glacier in the Stok Range during 2015–2019. Mehta and others (2021) reported a glaciological mass balance of −0.36 ± 0.04 m w.e. a⁻¹ for the Pensilungpa Glacier in the Suru Basin during 2016–2019. Srivastava and others (1999) reported a slight mass loss of −0.11 m w.e. a⁻¹ for the Rulung Glacier in the eastern Zanskar region during the period from 1979 to 1981. Together, these three glaciers cover only ~17 km² out of ~3500 km² glacierized area in the whole of Ladakh, accounting for less than 1% of the total glacier cover. Hence, existing field-based measurements are insufficient to represent the glacier mass change rates in the region.

In this context, geodetic mass balances computed using multi-temporal satellite-based digital elevation models (DEMs) provide a unique opportunity to assess regional and sub-regional changes. Since the accuracy of the geodetic mass balance approach improves as the length of observation data increases (~10 years or more), it currently stands as one of the most popular methods for estimating glacier mass balance using satellite data (Zemp and others, 2009; Berthier and others, 2023). Additionally, it is recommended to compare the existing glaciological mass balance series with the geodetic estimates, as the glaciological mass balances carry large uncertainties (Zemp and others, 2013).

Existing geodetic mass balance studies conducted in or around the Ladakh region suffer from one or more of the following drawbacks: (i) limited to small areas or selected glaciers, (ii) based on coarse resolution satellite data, for example, 90 m DEMs, (iii) limited temporal coverage. Table 1 presents the key information from existing geodetic studies conducted in and around the Ladakh region. Most studies, except the work by Abdullah and others (2020), focused on small to medium-sized areas or selected glaciers, lacking a comprehensive representation of glacier response from the region. The complete spatial representation is essential, especially in Ladakh's unique cold-desert climate, for comprehending the overall response of glaciers to climate change. The Ladakh-wide glacier mass change rates from Abdullah and others (2020) were based on 90 m resolution DEMs, which pose challenges in detecting local-scale variations and topographic variability within glaciers, especially for small-sized glaciers (<1 km²) constituting over 80% of the regional glacier cover (Schmidt and Nüsser, 2017). While very high-resolution DEMs would be ideal, 30 m resolution DEMs still provide a better resource for adequately capturing the mass balance characteristics of small-sized glaciers. Furthermore, existing studies cover temporal periods until 2012, and thus do not account for the recent heat waves and warmest years on record globally (WMO, 2022, 2023).

Table 1. Compilation of the existing remote sensing DEM-based glacier geodetic elevation change studies around the Ladakh region

WL/EL refers to the western/eastern Ladakh, and their spatial extents are shown in Figure 1. The spatial domains of previous studies are shown in Fig. 11 (in the Appendix).

^a Geodetic mass balance (m w.e. yr⁻¹).

Bhushan and others (2018) conducted a statistical analysis of the relationship between morphological variables and glacier mass balance in the Zanskar region, Ladakh but only considered 43 glaciers larger than 1 km². A survey of existing studies (Table 1) yields a general mass balance pattern of Ladakh glaciers, but the role of climatic and nonclimatic drivers controlling mass balance variability remains unknown. Moreover, larger and smaller glaciers (<1 km²) exhibit distinct climate responses and mass balance characteristics (e.g. DeBeer and Sharp, 2009; Shean and others, 2020). Initial studies suggest a higher topographic control on the mass balance characteristics of smaller glaciers compared to the general climate influence (DeBeer and Sharp, 2009; Hussain and others, 2022). Given the predominance of smaller glaciers in Ladakh, studying the response of small and large glaciers separately is necessary to better understand the relative importance of climatic and nonclimatic drivers controlling mass balance variability.

Additionally, previous studies have spatial data gaps due to limited spatial coverage of DEMs acquired in a single-year (elevation data acquired in a single year). To overcome this limitation, we incorporate alternative data sources such as ICESat-2 laser altimetry data, into our study (e.g. Berthier and others, 2023).

Recent large-scale and global geodetic glacier mass balance studies, predominantly using ASTER DEMs (e.g. Brun and others, 2017; Shean and others, 2020; Hugonnet and others, 2021), have effectively addressed spatial gaps and provided broad insights into glacier mass balances. These studies also play a critical role in estimating sea-level rise contributions and projecting future scenarios under various climate change conditions (Kraaijenbrink and others, 2017; Rounce and others, 2023). The latest study by Hugonnet and others (2021) provides mass balance estimates up to 2019. Consequently, a comprehensive study, backed by current regional knowledge and employing meticulous methodologies for DEM generation and elevation change computation across all Ladakh glaciers, is imperative to comprehend the up-to-date characteristics of glacier mass balance and its response to climate change.

In this study, we adopted a manual approach, carefully selecting stereo DEMs with high accuracy from a single-year acquisition during the ablation period to estimate the mass balance of nearly all Ladakh glaciers. A crucial subsequent step involves a comparative analysis between mass balance values derived from the manual approach and those generated by automated time series methodologies (e.g. Shean and others, 2020; Hugonnet and others, 2021). This comparison is particularly relevant in the context of democratizing the approach, as it allows for the utilization of single-year DEMs in similar studies conducted for different time periods, not aligning with available mass balance datasets (e.g. Hugonnet and others, 2021). Moreover, automated methods demand the processing of vast amounts of data and substantial computing resources.

Given the limited spatiotemporal understanding of Ladakh's glacier mass balances, a restricted comprehension of the utility of global-scale elevation change datasets for sub-regional scale applications, and an absence of a comprehensive understanding of factors controlling mass balance (both climatic and nonclimatic), this study aims to:

• • Estimate the surface elevation change and mass balance of nearly all glaciers in the Ladakh region, covering more than 3000 km² of glacierized area, by differencing 30 m resolution DEMs between 2000 and 2021,

• • Compare the SRTM-ASTER-based mass balances (2000–2021) with SRTM-ICESat-2 altimetry-based results (2000–2021) to verify the recent single-year regional mass balance estimates,

• • Investigate the governing factors, both climatic and nonclimatic, that control the spatial variability of glacier mass balances.

2. Study area

The Ladakh region, situated in the western Himalaya (Fig. 1), is characterized by a high-altitude cold desert landscape, with an average elevation of over 3000 m a.s.l. and featuring prominent mountain peaks, such as Nun, Kun, Stok Kangri, Kang Yatze, all soaring between 6000 and 7000 m a.s.l. The region is primarily influenced by the western disturbances during winter months, which contribute to more than two-thirds of the annual precipitation (Lee and others, 2014). Indian summer monsoon penetration is significantly weaker in these areas due to the high-mountain reliefs which act as formidable orographic barriers (Dimri, 2021; Phartiyal and Nag, 2022).

Figure 1. Ladakh region, encompassing both the WL and EL sub-regions, showing the glaciers assessed in this study. The maroon dashed polygon represents the ASTER DEM extent, while the pink and orange dashed polygons correspond to the footprints of the HMA DEM and Pléiades DEM, respectively. The arrow indicates the Stok Glacier, and glacier outlines are based on RGI v6.0 (blue). The background is the hillshade view, derived from the SRTM DEM. The inset shows the study area within the HMA domain, with the study region (red) and major rivers (cyan) marked.

Given the extensive coverage of the Ladakh region, we divided it into two sub-regions, namely eastern Ladakh (EL) and western Ladakh (WL), to analyse the spatial variability of glacier mass balances. Basin-wise, EL encompasses the majority of the Leh, Tsokar and Tsomoriri basins, along with smaller parts of the Pangong and Shayok basins. Whereas WL covers most of the Zanskar and Suru basins, with minor contributions from the Drass and Shingo basins. All the rivers in this region ultimately flow into the Indus. For a detailed sub-division of river basins around the Ladakh region, refer to Soheb and others (2022). EL exhibits a notably drier climate and a greater range of temperatures than WL (Table 2).

Table 2. Physiographic and climate characteristics of WL and EL

Glacier statistics are based on RGI v6.0 and Herreid and Pellicciotti (2020). Debris cover details are taken from Scherler and others (2018) and Herreid and Pellicciotti (2020). Climate statistics are based on Drass and Leh stations from WL and EL, respectively and extracted from previous studies (e.g. Osmaston and others, 1994; Raina and Koul, 2011; Lee and others, 2014; Mehta and others, 2021).

According to the Randolph Glacier Inventory (RGI v6.0; RGI Consortium, 2017) and the inventory from recent satellite images (Soheb and others, 2022), the glacierized area in EL and WL is over 700 km² and 2500 km², respectively (Fig. 1 and Table 2). In terms of morphological setting, individual glacier sizes in EL are significantly smaller, with over 80% of them measuring less than 0.75 km² and located at very high elevations, terminating at an average of ~5500 m a.s.l. (based on RGI v6.0). On the other hand, in WL, individual glacier sizes vary widely, with more than 30 glaciers covering an area of 10 to 69 km², extending to lower elevations and terminating at an average elevation of ~4800 m a.s.l. (RGI v6.0). Glaciers in WL are often covered by debris in their ablation zones, consisting of about 24% of the total glacierized area, whereas glaciers in EL are less debris covered, consisting of 16% (Scherler and others, 2018; Herreid and Pellicciotti, 2020). Only two glaciers in the study area, Stok (EL) and Pensilungpa (WL), have glaciological mass balance records for 5 years or less and are currently being measured (Soheb and others, 2020; Mehta and others, 2021).

3. Datasets and method

3.1. DEM (elevation) data

Glacier surface elevation changes between 2000 and ~2017/2020/2021 were analysed using the DEMs generated in this study and those available for the region. For the year 2000, the SRTM DEM (SRTM GL1; source: United States Geological Survey) was used, which was acquired through interferometric synthetic aperture radar (InSAR) data in X- and C-band radar frequencies between 11 and 22 February 2000 (Farr and others, 2007). The void-filled SRTM GL1 data has a spatial resolution of 30 m with a reported vertical accuracy of 6.9 ± 0.5 m (Mukul and others, 2015). This SRTM DEM was subtracted from the recent DEMs (2020/2021), generated from the ASTER Level-1A V003 (L1A) stereo pair images. ASTER stereo image coverage over WL was incomplete in 2020/2021, and no suitable images were available in 2018–2019, so 2017 stereopair imagery was used to complete the DEM coverage for WL.

In addition to the ASTER DEMs, we used a Pléiades 2-m resolution DEM from 2021 (source: National Centre for Space Studies, France; Berthier and others, 2014) and an HMA 8-m DEM from 2013 (source: National Snow and Ice Data Center, NSIDC; Shean and others, 2020) to compare the geodetic and glaciological mass balances on the Stok Glacier. Pléiades and HMA DEMs were used to closely cover the glaciological mass balance measurement periods. No suitable ASTER DEM was available for the same period.

As the ASTER DEMs over Ladakh covers multiple years, ICESat-2 laser altimetry data were used to compare the sub-regional geodetic mass balance from SRTM-ASTER and SRTM-ICESat-2 laser altimetry. ICESat-2, the successor mission to ICESat (refer to Kääb and others, 2012), has been acquiring data since 15 September 2018. We used ICESat-2 L3A ATL06 (land ice height; version 5) product, which comes with a window size ranging from 40 to 80 m to segment photons from global geo-located photon data (ATL03) to estimate the geolocated land ice surface height above the WGS 84 ellipsoid. ATL06 data has a spatial resolution of 20 m and an elevation retrieval accuracy of 0.1 m (Smith and others, 2019; Fan and others, 2022; Zhao and others, 2022).

3.2. ASTER DEM generation and co-registration

We selected 16 stereo pair images from the end of the hydrological year (September/October; Table 4 in Appendix) with less than 20% cloud cover, obtained from NASA's Earthdata portal, to derive DEMs employing the method outlined by Shean and others (2020) with the Ames Stereo Pipeline (ASP) v3.0.0 (Beyer and others, 2018; Shean and others, 2016). The generated DEMs were posted at a spatial resolution of 30 m, with elevations relative to the WGS84 ellipsoid and Universal Transverse Mercator (UTM) zone 43N projection (EPSG: 32 643). Subsequently, we co-registered each ASTER DEM to the SRTM DEM by iteratively minimizing slope and aspect-dependent elevation differences over nonglacierized stable terrain (static) areas, as defined by the RGI v6.0 glacier polygons. This process followed the Nuth and Kaab method (Nuth and Kääb, 2011) and was implemented using the demcoreg package (Shean and others, 2019). The same approach was employed to co-register the Pleiades DEM to the HMA 8 m DEM. Table 4 provides the horizontal (x and y) and vertical (z) offsets, the median elevation difference before and after co-registration and normalized median absolute deviation (NMAD) over static surfaces of all DEMs.

3.3. Glacier elevation and volume change, and mass balance calculation

We differenced each pair of co-registered DEMs to get the elevation change map for the study area over the period 2000 to ~2017/2020/2021. Before averaging elevation changes, we discarded all unexpected or spurious elevation change pixels exceeding ± 100 m, both for glacierized and nonglacierized terrains, following previous studies (Brun and others, 2017; Lv and others, 2020). Additionally, pixels with surface slopes > 45° (calculated from the SRTM DEM) were excluded, as DEM errors tend to increase rapidly with slope (Berthier and Brun, 2019). Glaciers truncated at the boundary of a DEM are also excluded from elevation change calculations to prevent potential bias, especially if only their accumulation or ablation area is sampled (Gardelle and others, 2013). Processed elevation change raster tiles were then mosaicked to produce a single elevation change raster for the EL and WL sub-regions using ASP's dem_mosaic utility (Shean and others, 2020). The mosaicked elevation difference raster was subsequently used for volume change and mass balance calculations.

We computed glacier volume change and mass balance at two spatial scales: sub-regions (EL and WL) and individual glacier polygons. The total volume change (Δv, in m³) for the entire period was calculated by multiplying the mean glacier elevation change (Δh, in m) for all individual glaciers or total glacierized areas (A, in m²):

(1)$${\Delta v = \Delta h \cdot A} $$

The volume change was then converted to specific glacier geodetic mass balance rate (Δm, in m w.e. a⁻¹):

(2)$$\Delta m = \Delta h\cdot \rho \cdot \Delta t$$

where, ρ is the density, and Δt is the time interval of DEM difference in years. For all the calculations (both individual glacier and spatial units), we considered a mean density of 850 kg m⁻³ following Huss (2013).

We did not impose a minimum glacier area threshold in regional mass balance or elevation change calculation and consider all the glacier polygons in RGI v6.0 for the entire Ladakh region. The outlines of glaciers with surface area > 1 km² over the WL were obtained from Herreid and Pellicciotti (2020) for better accuracy. To understand the glacier size/area sensitivity in regional mass balance aggregation, we provided estimates categorized by size groups in Sect. 4.2. This was done to account for the different responses of glaciers of varying sizes.

3.4. ICESat-2 laser altimetry data processing

We used Python's icepyx library (Scheick and others, 2019, 2023) to download ALT06 land ice height data from NSIDC for the latest available period over the ablation season of 2021 between 1 August and 30 September, to calculate the surface elevation change. The ablation season of 2021 was chosen to compare mass balance estimates from different datasets (SRTM-ASTER and SRTM-ICESat-2) over roughly the same cumulative time intervals.

ALT06 data were processed by separating them into each beam (ground tracks). Elevation data were extracted from SRTM GL1 DEM corresponding to ICESat-2 data points to calculate elevation change between 2000 and 2021. Before the final elevation change calculation, we discarded elevation change values larger than (i) ± 100 m, assuming that this value was the maximum possible elevation change, and (ii) five standard deviations (SD, e.g. values > 5 × SD) to filter out outliers (Hugonnet and others, 2022). Processed elevation change values were segregated into on-ice or off-ice elevation changes based on glacier outlines. On-glacier data points were used for mass balance calculation, whereas off-glacier points were used to assess the associated uncertainty based on the standard approach (see Sect. 3.5.3).

After data processing and cleaning, a total of 37 206 data points (2000–2021; 21.5 years) were used for elevation change calculation in the WL, of which 9992 (27%) were on ice. For the EL, 5412 data points were used for elevation change calculation, covering the Kang Yatze and Stok ranges, of which 836 (~16%) were on ice. Uncertainty estimation involved utilizing a total of 27 214 and 4576 off-glacier (static area) data points for WL and EL, respectively, to calculate std dev. statistics.

3.5. SRTM C-band penetration bias correction, seasonality correction, and overall uncertainty estimation

3.5.1. Penetration bias correction

The C-band radar wave penetrates to varying depths in different glacier surfaces (snow, ice, and supraglacial debris), introducing biases in the geodetic elevation change measurements using SRTM-C DEMs (Gardelle and others, 2013; Vijay and Braun, 2016; Berthier and others, 2023). It is essential to address the penetration bias before calculating and interpreting glacier elevation change measurements. Following Kääb and others (2012), we corrected for the penetration bias over various glacier surfaces. The correction values were 2.3 ± 0.6 m, 1.7 ± 0.6 m and 0.4 ± 0.8 m for snow, clean ice, and debris-covered ice, respectively, in both WL and EL sub-regions. This penetration bias estimate was originally calculated by Kääb and others (2012) for the Jammu and Kashmir region, which encompasses the Ladakh region.

For the penetration bias adjustment over various surfaces, debris-covered areas for the WL were classified using the debris outlines from Herreid and Pellicciotti (2020) for glaciers > 1 km², and for the remaining glaciers, from Scherler and others (2018). Data on debris cover for small glaciers (<1 km²) were unavailable in the former source. Snow and ice surfaces were classified using accumulation and ablation area polygons separately provided by Herreid and Pellicciotti (2020) for glaciers > 1 km², while for the remaining glaciers, the median elevation was used. The debris, ablation and accumulation area polygons from Herreid and Pellicciotti (2020) are considered more reliable due to manual inspection compared to automatic delineation. For the EL, debris-covered area polygons were obtained from Scherler and others (2018) and snow/ice surfaces were classified using glacier median elevations.

We adjusted the penetration depth bias for elevation change values obtained by differencing SRTM and ICESat-2 data points, using the same penetration bias values for different surfaces applied in the SRTM and ASTER differencing (Kääb and others, 2012).

3.5.2. Seasonality correction

We did not apply any seasonality correction since the majority of the DEMs were acquired around the end of the ablation season (refer to Table 4 for the acquisition dates). Specifically, a lone ASTER DEM acquired on 19 December 2021 was exclusively utilized for the Stok Range areas (Table 4). Given the notably low annual precipitation in the Leh area, situated ~20 km from the Stok Range, during October-December (<10 mm month⁻¹; Lee and others, 2014), we did not consider a seasonality correction for this particular DEM as well.

3.5.3. Glacier elevation change and mass balance uncertainty

The uncertainty in glacier elevation, volume and mass balance changes are primarily attributed to the uncertainties related to the DEMs used (SRTM GL1 and ASTER). To account for the uncertainty in glacier mass balance, we followed approaches similar to recent geodetic mass balance studies using similar DEMs (e.g. Fischer and others, 2015; Brun and others, 2017; Lv and others, 2020). We employed a simplified implementation of Rolstad and others (2009) method to estimate the spatial uncertainty in glacier elevation change (σΔh, in m).

In this approach, considering the influence of autocorrelation between DEM pixels, a spatially correlated area (A_corr, in m²) is required, defined as:

(3)$${A_{corr}} = \pi \cdot L^{2} $$

where L is the decorrelation length between pixels in the glacier elevation change map. We considered the value of L as 500 m following Brun and others (2017). This value is typically determined by analysing the semivariance (variogram) of randomly selected data pairs using DEMs (Rolstad and others, 2009; Menounos and others, 2019). The value of A_corr was ~0.79 km².

For glaciers with a surface area (A) > A_corr, we scaled the uncertainty of elevation change (σΔh) following Rolstad and others (2009):

(4)$$\sigma \Delta h = SD\sqrt {A_{corr}/5A} $$

where SD represents the std dev. of elevation changes over stable surfaces. For glaciers with a surface area smaller than A_corr, we assume that σΔh is equal to the SD (Shean and others, 2020). This involved considering stable terrain (off-glacier; shown in Fig. 12) statistics (SD) from a buffer region of 500 m surrounding each glacier polygon. These statistics were then applied to Eqn (4) for uncertainty computation.

For the region-wide analysis of elevation change uncertainty, we considered stable terrain statistics from all pixels within 50 m hypsometric bins across the entire study area and scale SD using Eqns (5) and (6) (Gardelle and others, 2013; Bolch and others, 2017):

(5)$$\sigma \Delta h = \displaystyle{{SD} \over {\sqrt {N_{eff}} }}$$

(6)$$N_{eff} = \displaystyle{{N_{tot} \times PS} \over {2d}}$$

where N_eff and N_tot are the effective and total number of observations, PS is the pixel size (30 m) and d is the distance of spatial autocorrelation, which was considered to be 500 m. We calculated the weighted averages of regional σΔh with respect to glacier hypsometry of 50 m elevation bands to quantify the uncertainty in elevation change over glacierized areas, ensuring that each glacierized elevation bin has a different uncertainty depending on the hypsometric distribution.

The total volume change uncertainty was then obtained as:

(7)$$\sigma \Delta v = \sigma \Delta h\cdot A$$

The overall uncertainty of the annual geodetic mass balance (σΔm, in m w.e. a⁻¹) was then calculated using the uncertainties discussed above and defined as:

(8)$$\sigma \Delta m = \displaystyle{{\sqrt {{\left({\Delta v\cdot \displaystyle{{\sigma_\rho } \over A}} \right)}^2 + {\left({\sigma \Delta v\cdot \displaystyle{\rho \over A}} \right)}^2 + {\left({\Delta v\cdot \rho \cdot \displaystyle{{\sigma_A} \over A}} \right)}^2} } \over {\Delta t}}$$

where ρ is the density (0.85) with an uncertainty (σ_ρ) of 0.06 (7% of ρ; Huss, 2013), and σ_A is the area uncertainty, representing a one-pixel buffer (30 m) following Dussaillant and others (2018). The level of area uncertainty is considered reasonable, except for debris-covered glaciers (Paul and others, 2015). We assumed that these three primary error sources (σΔh, σ_ρ and σ_A) are independent and uncorrelated (Fischer and others, 2015; Brun and others, 2017; Zhao and others, 2022). Throughout the text, uncertainties are reported at the 68% confidence interval (1-sigma) level.

The uncertainty associated with the penetration bias correction was considered to be ± 1.5 m, following Gardelle and others (2013). Subsequently, this penetration error was incorporated into the overall mass balance uncertainty, adhering to standard principles of error propagation (Eqn (8)).

3.6. Climate data

Long-term ERA5-Land reanalysis air temperature and precipitation datasets of 9-km resolution from 1950 to 2020 were used, to analyse the climate influence on the glacier mass balance over the study area. ERA5-Land datasets were downloaded freely from the Copernicus Climate Data Store and processed using Python's xarray package (Hoyer and others, 2022). The reliability of the ERA5-Land data was assessed using the nearest available gauged data from the Global Historical Climatology Network – Monthly (GHCN-M) version 4. We used GHCN-M's air temperature and precipitation from four available stations, for example, Srinagar (1587 m a.s.l.), Shimla (2202 m a.s.l.), Leh (3514 m a.s.l.) and Shiquanhe (4280 m a.s.l.). While three stations are situated within the western Himalaya, Shiquanhe is located on the western edge of the Tibetan Plateau, ~80 km east of the study area boundary.

To assess the reliability of the ERA5-Land datasets, monthly air temperature and precipitation values at the nearest gridpoint to each station location were compared with the station-based values from Srinagar, Shimla, Leh and Shiquanhe sites from GHCN-M. We found a significant strong to moderate coefficient of determination (r ²) between ERA5-Land and GHCN-M air temperature and precipitation data (temperature r ² = 0.93–0.97; precipitation r ² = 0.32–0.67), with a lower root mean squared error (RMSE; temperature RMSE = 2.3–14.9°C; precipitation RMSE = 17.2–96.7 mm). Figure 2 shows the comparison and reliability statistics.

Figure 2. Comparison of monthly air temperature (n = 600) and precipitation (n = 858/Srinagar, 456/Shimla, 236/Leh, 324/Shiquanhe) data between ERA5-Land (9-km resolution) and GHCN-M (gauged) data at Srinagar, Shimla, Leh and Shiquanhe sites over the period 1950–2021. The station locations are shown in Fig. 7a. Note the difference in the x- and y-axis across subplots.

Subsequently, we computed the linear regression trend for mean annual air temperature and precipitation in each data grid around the study area. The statistical significance of these trends was assessed using a Student's t-test, considering trends significant when confidence levels exceeded 95% (p = 0.05).

3.7. Morphological (nonclimatic) data

We investigate the influence of morphological variables that have been identified as influential in glacier mass balance in previous studies in the HK region (e.g. Salerno and others, 2017; Bhushan and others, 2018; Brun and others, 2019). These variables include area, median slope, median aspect (converted to northness or cos(aspect)), median and minimum elevation of the glaciers, elevation range (the vertical difference from the glacier's top to bottom), and percentage of supraglacial debris cover. The median values of slope, aspect and area were computed using the SRTM DEM and then masked to individual glacier polygons using the rasterstats module in Python. The remaining morphological variables were obtained from the inventories corresponding to the respective glaciers. To establish the statistical relationships, the Pearson correlation analysis was conducted and represented through a heatmap.

4. Results

4.1. Surface elevation change and mass balance over WL/EL

The results show widespread glacier surface thinning across the entire Ladakh region over the past two decades, with noticeable spatial variability (Figs 3 and 4; Table 3). In the WL, we observe a cumulative glacier surface elevation change of −7.75 ± 1.53 m (−0.44 ± 0.09 m a⁻¹) for the period 2000–2017, for a total glacierized area of 2106 km² (Fig. 3). This corresponds to a mean mass loss of −0.37 ± 0.09 m w.e. a⁻¹. Over the same period, the total glacier volume loss was −16.32 km³, corresponding to a mean loss of −0.93 km³ a⁻¹. For a glacierized area of 1085 km², which was not covered in 2017 DEMs (Fig. 3), the cumulative glacier elevation change was found to be −8.98 ± 1.56 m (−0.41 ± 0.07 m a⁻¹) for the period 2000–2021. This translates to a mean mass loss of −0.35 ± 0.07 m w.e. a⁻¹ for WL glaciers during 2000–2021. During the same period, the total glacier volume loss in this area was −9.74 km³, corresponding to a mean loss of −0.45 km³ a⁻¹.

Figure 3. Map illustrating glacier elevation changes between February 2000 and October 2017 for the WL (covering 2106 km² of glacierized area). A part of the glacier elevation change map is for February 2000–October 2021 (dashed red outline; 1085 km²). Glacier outlines are based on Herreid and Pellicciotti (2020) for > 1 km² glaciers and from RGI v6.0 for the rest of the glaciers. The lower left inset shows an enlarged view of the Prul Glacier, the only surging glacier in the study area. The top-right inset shows an enlarged view of the Durung Drung and Pensilungpa glaciers area.

Figure 4. Map illustrating glacier elevation changes between February 2000 and September 2020 for the EL (covering 445 km² of glacierized area). A part of the glacier elevation change map is for February 2000–October 2021 (dashed red outline; 272 km²). Glacier outlines are based on RGI v6.0. Insets show an enlarged view of the Stok region where glaciological mass balance records are available. The bottom inset includes the Lato/Kang Yatze area, considering its extensive glacierized area.

Table 3. Glacier surface elevation changes (Δh) for the WL and EL between 2000 and 2017/2020/2021

The total glacierized area and error statistics (median of off-glacier stable terrain and NMAD) are also given. NMADs values pertain to stable areas (Fig. 12) after slope filtration, with pixels > 45° slope discarded.

The magnitude of surface thinning in the EL was about 30–50% lower compared to the WL (Table 3). For the period 2000–2020, the estimated cumulative glacier elevation change was −5.00 ± 1.51 m (−0.25 ± 0.07 m a⁻¹) for 445 km² glacierized area (Fig. 4). This corresponds to a mean mass loss of −0.21 ± 0.07 m w.e. a⁻¹. During the same period, the total glacier volume loss was −2.22 km³, corresponding to a mean loss of −0.11 km³ a⁻¹. For the remaining glacierized area of 163 km², which was not covered in 2020 DEMs (Fig. 4), the cumulative glacier elevation change was −8.55 ± 1.19 m (−0.39 ± 0.05 m a⁻¹) over 2000–2021. This translates to a mean mass loss of −0.33 ± 0.05 m w.e. a⁻¹ for EL glaciers during 2000–2021. The total glacier volume loss in this area during the same period is −1.40 km³, corresponding to a mean loss of −0.06 km³ a⁻¹.

According to the surging glacier inventory of HMA (Guillet and others, 2022), the only surging glacier in the study area is the Prul Glacier (RGI60–14.18918; shown in Fig. 3; with a total glacierized area of about 50 km²), located in the WL. The mass balance for this glacier was −0.23 ± 0.10 m w.e. a⁻¹ during 2000–2017.

To understand the altitudinal dependence of the glacier elevation change, we analysed the area-weighted elevation change rates within 50 m altitudinal bins (hypsometry) of the glacierized areas. We observed an overall stronger magnitude of elevation change rates (mass balance gradient) in the WL compared to that of the EL (Fig. 5). Most of the glaciers in the WL terminate at lower elevations than those in the EL, exposing them to warmer summer temperatures, and resulting in greater melt and mass loss magnitude. The debris-covered areas of WL glaciers, constituting about 24% of the total area, experienced a strong negative elevation change rate of −0.99 m yr⁻¹, although slightly reduced elevation change rates were observed at the glacier's terminus areas ~3500–3700 m a.s.l. (lowest elevation bins; Fig. 5a). Elevation changes over debris-covered areas in the EL glaciers (16%) were significantly lower at a rate of −0.25 m a⁻¹. Note here that the EL glaciers are located at higher elevations.

Figure 5. Hypsometric elevation change rates of the glacierized areas in the WL (Panel a; 2106 km²) and EL (Panel b; 445 km²) using on 50 m altitudinal bins. The blue vertical lines represent the median elevation of the glaciers in the corresponding sub-regions.

We compared both SRTM C-band penetration-corrected and uncorrected elevation change rates for both sub-regions to quantify the penetration bias (shown in Fig. 5). Our observations revealed that the elevation change rates for the entire glacierized area are underestimated by about 10–20% if the penetration bias is uncorrected, consistent with the previous finding of ~20% underestimation by Vijay and Braun (2016) for the neighbouring Lahaul-Spiti glaciers.

4.2. Contributions to region-wide mass balances by glaciers of varied size distributions

To examine the influence of glacier size, categorized by area, on regional mean mass balances, we categorized glaciers into different area groups (Fig. 6). Smaller glaciers in both sub-regions exhibited scattered (less consistent) mass balance patterns with higher uncertainties (see Fig. 13 for uncertainties). For instance, small glaciers (<1 km²) in the WL displayed a mean mass balance rate of −0.18 m w.e. a⁻¹, with an uncertainty range of ± 0.40 m w.e. a⁻¹. Similarly, small glaciers (<1 km²) in the EL showed a mean mass balance rate of −0.17 m w.e. a⁻¹, with an uncertainty range of ± 0.28 m w.e. a⁻¹. Here, it is important to note that the higher uncertainty associated with small-sized glaciers is mainly attributed to the uncertainty estimation method employed in this study. In Eqn (4), the SD, representing elevation change uncertainty, is scaled for glaciers with an area equal to or larger than the spatially correlated area (A_corr ≈ 0.79 km²). However, for glaciers with a smaller area, it remains equal to the SD. That means, as the area of the glaciers in question increases, the denominator in Eqn (4) becomes larger, leading to a decrease in the elevation change uncertainty (σΔh). Consequently, aggregating with larger glaciers significantly reduces mass balance uncertainty in both sub-regions (Fig. 6).

Figure 6. Mass balance rates for the WL (top rows) and EL (bottom rows) glaciers plotted against their surface area. Four area categories are shown for each sub-region. In each panel, the number of glaciers and mean mass balance rates are displayed. The colour code indicates the glacierized area sampled for surface elevation change calculations.

For the largest area group in the WL (>5 km²), the uncertainty range was considerably lower at ± 0.14 m w.e. a⁻¹ for a mean mass balance rate of −0.46 m w.e. a⁻¹. Similarly, in the EL's largest area category, the mass balance uncertainty range was ± 0.12 m w.e. a⁻¹ for a mean mass balance rate of −0.43 m w.e. a⁻¹. Overall, we note that the small-sized glaciers exhibit larger uncertainties in mass balances, while larger glaciers show lower uncertainties in both sub-regions, as discussed previously. The mean mass balance of smaller glaciers (<1 km²) in both sub-regions is roughly half of that of larger glaciers (glacier area > 3 or 5 km²; Fig. 6), consistent with large-scale geodetic results from HMA (e.g. Shean and others, 2020). Another crucial observation is that the aggregated mean mass balances for the smallest and largest area categories in the WL (−0.18 ± 0.40 and −0.46 ± 0.14 m w.e. a⁻¹, respectively) and EL (−0.17 ± 0.28 and −0.43 ± 0.12 m w.e. a⁻¹, respectively) are identical. However, it is noteworthy that the largest area category in the EL sub-region comprises only four glaciers, which may not be considered a true representation of such an area category. Overall, the magnitude of negative mass balance for individual WL glaciers is higher (Fig. 6a–d).

4.3. Comparison between ICESat-2-derived and DEM differencing-derived mass balances

To assess the consistency across different datasets during a comparable time frame, we compared glacier mass balances derived from SRTM-ASTER DEMs with those obtained from SRTM-ICESat-2 data points. Given the incomplete coverage of ASTER DEMs for the recent year 2021, especially in the WL glacierized area (as mostly covered by 2017 DEMs; Sect. 3.1), cross-checking the results using ICESat-2 data presented a valuable opportunity. Elevation change data points obtained by subtracting SRTM from ICESat-2 datasets are shown in Fig. 14.

For the period 2000–2021, the regional glacier mass balance rates derived from subtracting SRTM from ICESat-2 data points were −0.40 ± 0.53 and −0.25 ± 0.25 m w.e. a⁻¹ for the WL and EL, respectively. In comparison, the regional glacier mass balance rates obtained by subtracting SRTM from ASTER DEMs were −0.35 ± 0.07 and −0.27 ± 0.06 m w.e. a⁻¹ for the WL and EL, respectively. Both methods yielded comparable results within their uncertainty limits, affirming that the mass balance from ICESat-2 data serves as a reliable representation of regional glacier mass balance for the recent time period. These findings support the region-wide mass balance results obtained from SRTM-ASTER DEMs over the last two decades.

4.4. Trends in air temperature and precipitation across the study area

Over the past 71 years (1950–2021), mean air temperatures in the study area have increased by 0.10°C decade⁻¹ (Fig. 7a). While the warmest year on record was observed in the past two decades, temperatures between 2000 and 2020 have remained relatively stable (Fig. 7c). The absolute changes in mean decadal air temperature align closely with the long-term upward trend (Fig. 7e).

Figure 7. Air temperature and precipitation trends are shown for the study area (dashed rectangle in panels a and c) based on ERA5-Land datasets. Leh, Srinagar, Shimla and Shiquanhe GHCN-M gauging locations are marked with L, SR, SH and SQ, respectively, in panel a. The mean annual air temperature and precipitation from spatially averaged ERA5-Land data for the EL and WL sub-regions during 1950–2020 (long) and 2000–2020 (short) are also shown (b and d). The statistically significant trend values at the 95% confidence interval are marked with hatching in panels a and b. Linear trendlines are overlain for the entire long-term period (1950–2020; 70 years) and recent (2000–2020; 20 years) periods to highlight recent changes. Decadal changes in air temperature and precipitation for the EL and WL (e and f).

Mean annual precipitation between 1950 and 2021 decreased by about 10 mm decade⁻¹, with statistically significant values observed only over the WL (Fig. 7b). Intriguingly, the WL exhibited a slight positive precipitation trend over the last two decades, particularly in the 2010s (Fig. 7d), contrary to its long-term negative trend. In contrast, long-term precipitation in the EL exhibited minimal change. While most precipitation trends lacked statistical significance, a slight increase in precipitation was observed, especially in the north-eastern part of the EL (Fig. 7b).

Overall, the study area has witnessed widespread warming over the past decades, accompanied by a long-term decrease in precipitation primarily over the WL areas, while the EL areas have experienced minimal changes in precipitation.

4.5. Analysis of morphological (nonclimatic) variables and their influence on glacier mass balance

We observed varying influences of morphological variables on glacier mass balances in both sub-regions (Fig. 8, Figs 15 and 16). A strong relationship between glacier slope and mass balance was evident in both sub-regions, indicating that steeper glaciers in both areas exhibited a lower negative glacier mass balance.

Figure 8. A heatmap representation of the Pearson correlation matrix of glacier mass balance (Δm) and various morphological variables for the WL (n = 2203) and EL (n = 1264) sub-regions. Z _med, Z _min, Z _range and deb_percent represent glacier median elevation, minimum elevation, elevation range, and percentage of debris-covered area, respectively.

Elevation variables, such as median and minimum elevations, did not show a strong relationship with glacier mass balance in the WL sub-region. However, stronger relationships between median and minimum elevations and mass balances were observed in the EL sub-region. In the WL, debris cover is negatively correlated with mass balance, which suggests that glaciers with higher debris cover experience greater rates of mass loss.

Glacier surface area does not exhibit a strong relationship with glacier mass balance, especially in WL glaciers, but the relationship is slightly stronger in EL glaciers. Glacier elevation range and aspect do not show any relationship with glacier mass balances in both sub-regions.

In summary, the mass balance of EL glaciers demonstrates a stronger and more distinct relationship with morphological variables (slope, median and minimum elevations, and area) compared to WL glaciers. The mass balance of WL glaciers is distinctly related to mainly two variables: slope and debris cover.

5. Discussions

5.1. Results in the context of existing mass balance estimates

Glacier mass balance estimates in the Ladakh region have historically been limited by incomplete coverage, often focusing on specific glaciers or basins (refer to Table 1 for details about existing studies). This limitation arises mainly from the incomplete coverage of DEMs. In contrast, this study's mass balance estimates are based on 30 m resolution DEMs, providing nearly complete regional coverage and spanning a period of about 20 years or more between 2000 and 2021.

The results for both sub-regions are generally align with other studies conducted in and around the study area, as well as the western Himalaya-wide glacier mass balance results (Fig. 9). In the WL, where most existing estimates indicate a similar mean mass loss pattern (Fig. 9a), the results presented by Abdullah and others (2020) and Garg and others (2021), −0.99 ± 0.43 and −0.69 ± 0.28 m w.e. a⁻¹, respectively, appear relatively higher than our WL-wide geodetic mass balance rate of −0.35 ± 0.07 m w.e. a⁻¹, although they fall within the overlapping uncertainty range. For the EL, only Abdullah and others (2020) reported the geodetic glacier mass balance covering the entire sub-region (−0.39 ± 0.24 m w.e. a⁻¹ for 2000–2012), which is comparable to our estimate of −0.27 ± 0.06 m w.e. a⁻¹ for 2000–2020/2021 (Fig. 9b). Direct comparison of individual glacier mass balances from Hugonnet and others (2021) and Shean and others (2020) for the overlapping glaciers also demonstrates good agreement with our estimate (this study: −0.22 ± 0.36 m w.e. a⁻¹, Hugonnet and others (2021): −0.19 ± 0.25 m w.e. a⁻¹, and Shean and others (2020): −0.21 ± 0.44 m w.e. a⁻¹ for the WL, and −0.18 ± 0.27, −0.19 ± 0.22, and −0.16 ± 0.25 m w.e. a⁻¹ for the EL). Moreover, the direct comparison of individual glacier mass balances from Hugonnet and others (2021) and Shean and others (2020) with our estimates provides valuable insights into the robustness of calculating individual glacier mass balances in their large-scale automated workflow, cross-validated by our local-scale study. In summary, the comparison with existing studies suggests a general agreement with our estimate and indicates that glacier mass balance rates in the WL are slightly higher than those in the EL.

Figure 9. Glacier mass balance rates based on geodetic methods in the WL (a) and the EL (b) sub-regions from various studies. The legend provides information about the source, DEMs used, and region names. The spatial domains of previous studies conducted in and around the study area, as well as over the large-scale Himalayan sub-divisions, are shown in Fig. 11. Mass balance estimates shown here for Hugonnet and others (2021) and Shean and others (2020) correspond to the same glaciers as those for which we have provided estimates in our current study.

We attempted to compare the geodetic mass balance against the in situ glaciological series of the Stok Glacier (<1 km²), located in the EL, since the in situ mass balance record is available for five continuous years during 2014–2019 (Fig. 10). The Stok Glacier's five year mean glaciological mass balance was −0.39 ± 0.38 m w.e. a⁻¹ for the period 2014–2019 (Soheb and others, 2020). The mean geodetic mass balance over the period 2000–2021 (21 years) is calculated to be −0.44 ± 0.25 m w.e. a⁻¹. Direct comparison of results from both methods is, however, not straightforward because the in situ records cover only a quarter of the period for which our geodetic balance was calculated. To more closely represent the temporal period, we calculated the geodetic mass balance for an overlapping period of 2013–2021 using HMA 8-m and Pléiades DEMs. The mean geodetic mass balance for the period 2013–2021 (6.8 years) is found to be −0.27 ± 0.08 m w.e. a⁻¹. Both time-period geodetic estimates, 21 years and 6.8 years, provide similar results compared to the in situ glaciological record. Geodetic mass balance estimated using HMA and Pléiades DEMs exhibit a lower uncertainty range of ± 0.08 m w.e. a⁻¹ due to the high-resolution data (Fig. 10a). The comparison also highlights the suitability of the current glaciological measurement network and robust mass balance series of the Stok Glacier. This underscores the importance of examining and inter-comparing different methods, as recommended by Zemp and others (2013).

Figure 10. The mean (a) and cumulative (b) mass balance comparison from geodetic and glaciological methods for the Stok Glacier. DEM source, acquisition date, geodetic and glaciological mass balance values and differences are shown in the respective panels.

5.2. Complementary approaches for elevation change measurements

The single-year DEM differencing method, utilized in our study, and the temporal stacking of DEMs, employed in large-scale studies (i.e. Brun and others, 2017; Shean and others, 2020; Hugonnet and others, 2021), each has its own set of advantages and disadvantages. Temporal stacking requires substantial computing resources, which are not easily accessible to the broader community. In contrast, our approach, using nearly spatially complete single-year DEM differencing, yields an identical mass balance estimate for the Ladakh region. The inter-comparison with large-scale DEM stacking-based approach and reliable mass balance estimate provided in this study are particularly important for democratizing the single-year DEM-based method. This approach has the potential to be applied in estimating mass balances for other local-to-regional glacierized areas for different time periods not covered by existing large-scale datasets, such as those from Hugonnet and others (2021) and Shean and others (2020). Nevertheless, the single-year DEM differencing approach is not feasible for a chosen area or basin where DEM coverage in a single year is not complete, as was the case in the WL sub-region in this study.

As an alternative elevation data source, we utilized ICESat-2 altimetric datasets in our study to cross-verify the geodetic mass balance derived from the SRTM-ASTER DEM with SRTM-ICESat-2 based results. This approach demonstrated good agreement in estimating the region-wide glacier mass balance for a similar time period (Sect. 4.3). However, while ICESat-2 elevation measurements are valuable for region-wide assessments, they may not be suitable for obtaining individual glacier-wide mass balances for narrow, small-sized mountain glaciers due to the sparse footprints of ICESat-2 altimetric data (Treichler and Kääb, 2016; Fan and others, 2022), as opposed to the spatially complete DEMs.

5.3. Mass balance controlling factors: climatic and nonclimatic

Glacier mass balance fluctuation in the HK region is generally attributed to regional changes in air temperature and precipitation (e.g. Maurer and others, 2019; Bhattacharya and others, 2021). However, establishing a direct relationship between mass balance and climate variables in the Himalayan region is challenging due to the scarcity of long-term in situ meteorological data from high-altitude areas. To address this, we attempted to relate glacier mass balances to climatic changes using ERA5-Land reanalysis datasets, evaluated against four gauge-based long-term climate datasets (Sect. 4.4).

In general, glacier mass loss across the study area corresponds to widespread air temperature rise. We observed higher glacier mass loss in the WL (−0.35 to −0.37 m w.e. a⁻¹) than in the EL (−0.21 to −0.33 m w.e. a⁻¹), which is roughly aligned to a decrease in precipitation over the WL compared to the stable precipitation in the EL. Although long-term observational mass balance records are lacking in the region, historical geodetic mass balance records from Maurer and others (2019) for about 30–40 glaciers in the WL region (no records for the EL glaciers) showed a nearly doubled mass loss, reaching about −0.40 m w.e. a⁻¹ during 2000–2016 (recent period) in comparison to 1975–2000 (~−0.20 m w.e. a⁻¹; historical period). Similar observations were reported by Bhattacharya and others (2021) for selected glaciers in the central-western Himalayan region (Gurla Mandhata; ~300 km southeast of Ladakh), indicating higher mass loss during recent years, particularly over 2013–2019, compared to 1966–2000. This increased mass loss in the last two decades in and around the western Himalaya (Maurer and others, 2019; Bhattacharya and others, 2021), as also found in our study for the WL during 2000–2021, may be linked to decreased precipitation in the WL areas. This association is supported by our long-term climate data analysis (Fig. 7). Although reanalysis data has inherent uncertainties, in consistency with our analysis, several gauge-based studies (e.g. Bhutiyani and others, 2007, 2010; Shekhar and others, 2010; Dimri and Dash, 2012; Ren and others, 2017) reported a similar climatic fluctuation of strong and widespread temperature increase and heterogeneous precipitation/snowfall trends in and around the western Himalaya over the last few decades.

In addition to precipitation decrease (no-change) as an important factor for higher (lower) glacier mass loss in the WL (EL), the EL glaciers, being significantly smaller in size and situated at very high elevations (mean terminus elevation at ~5500 m a.s.l.), experience mostly negative temperatures, or the zero-degree isotherm reaches over glaciers for a very short period in summer. This likely contributes to the lower glacier mass loss in the EL compared to the WL. The EL glaciers have probably receded, and their overall geometry and local morphometric setting are now optimal to be close to an equilibrium state. This observation is consistent with DeBeer and Sharp (2009), who noted that small glaciers do not necessarily shrink in response to climate warming due to their favourable topographical locations in shadowed cirques and niches. Hussain and others (2022) found positive mass balance for some high-altitude fragmented glaciers in the Gangotri region of the central Himalaya, attributed to nonclimatic factors such as high elevation. Nonetheless, due to continued projected global warming, particularly in mountain areas, the EL glaciers, like other glaciers in the region, will experience increased mass loss in the future.

In addition to climatic influences, we found that nonclimatic (morphological) variables play a crucial role in explaining glacier mass balances in both sub-regions. Surface slope is a notable morphological factor explaining 30–42% of the variability in mass balances (Fig. 8, Figs 15 and 16), indicating that glaciers with flatter/gentler slopes experience more negative mass balances compared to steeper ones. This is because flatter glaciers' tongues are generally heavily debris-covered (due to favourable topography for debris accumulation and lower surface ice velocity due to gentler elevation gradient), leading to lower mass flux transport/divergence from the accumulation zone (e.g. Anderson and others, 2021; Miles and others, 2021; Rounce and others, 2021), which is reflected in the more negative net elevation change, especially evident in the WL glaciers, where the glaciers are large enough to experience such a phenomenon. Steeper glaciers, on the other hand, tend to have a larger portion of their surface area situated above the high-summer-temperature elevation line (zero-degree isotherm), resulting in lower mass loss. Brun and others (2019) also found such a relationship between glacier tongue slope and mass balance, where glacier tongue slope explains up to 49% of the mass balance variability across 12 HMA regions encompassing over 6000 glaciers. Additionally, Salerno and others, (2017) reported that a lower downstream surface gradient (flatter slope) leads to greater thinning of glaciers in the Everest region of the Nepalese Himalaya.

Further, in terms of nonclimatic variables, the WL glaciers generally feature a larger elevation range due to their extensive size, which ideally should contribute to less mass loss due to a greater area available for accumulation compensation through large mass turnover, provided that a large enough area is above the zero-degree isotherm. However, at the same time, they are more sensitive to air temperature increases due to more surface area (Oerlemans and Reichert, 2000). This sensitivity likely contributes to the WL's higher mass loss, where the elevation range of the glaciers is much larger than the EL glaciers. Notably, the presence of debris-covered tongues on the WL glaciers appears to play a significant role in the increased mass loss. Although the WL glaciers are relatively more debris-covered (24%) compared to the EL glaciers (16%), the larger debris-covered area does not provide the expected shielding effect against higher mass loss (Østrem, 1959; Nicholson and Benn, 2006; Vincent and others, 2016) (Fig. 17). One potential reason is that heavily debris-covered glacier tongues often feature numerous ice cliffs and supraglacial ponds, leading to higher (differential/localized) surface melt due to exposed ice surfaces and high energy absorbance capacity of supraglacial ponds, as indicated by previous studies (e.g. Steiner and others, 2015; Miles and others, 2018; Anderson and others, 2021). Additionally, the WL glaciers, especially the larger ones, terminate at lower elevations with a higher debris cover in their tongues and flatter tongue slopes (see Fig. 17), making them more sensitive to summer temperature.

Overall, nonclimatic factors, particularly the higher elevation of the EL glaciers, play a significant role in preventing accelerated mass loss, along with a relatively stable precipitation pattern in the EL sub-region over the past decades. On the other hand, debris cover and glacier slope stand out as nonclimatic factors that affect mass balance in the WL, with greater mass loss observed for low-angle and higher debris cover glaciers, alongside the significant precipitation decrease in the WL areas.

6. Summary and conclusion

In this study, we conducted a comprehensive analysis of glacier elevation change and mass balance in the Ladakh region, covering over 3000 km² of glacierized area. Leveraging 30 m resolution DEMs from SRTM and ASTER spanning the last two decades, we derived spatially complete mass balance estimates for nearly all glaciers in Ladakh. The SRTM-ASTER DEM differenced mass balance estimates were cross-verified against SRTM-ICESat-2-based mass balance estimates. Additionally, we systematically assessed both climatic and nonclimatic variables to elucidate the drivers of regional variability in glacier mass balances across Ladakh.

The key findings indicate that WL glaciers exhibited higher mass loss (0.35 ± 0.07 to −0.37 ± 0.09 m w.e. a⁻¹) compared to EL glaciers (−0.21 ± 0.07 to −0.33 ± 0.05 m w.e. a⁻¹) over the last two decades. Region-wide mass balances for both sub-regions were consistent with estimates derived from SRTM and ICESat-2 elevation datasets, affirming the reliability of our single-year DEM-based approach for both region-wide and individual glacier mass balances in Ladakh. While a general air temperature increase is an apparent factor contributing to widespread mass loss across the study area, and a decrease in precipitation in the WL areas is associated with higher mass loss, our findings highlight the significant role of nonclimatic factors. Specifically, the morphological setting of the EL glaciers, characterized by their existence at very high elevations and overall morphometric settings (i.e. smaller size), plays a crucial role in defying higher mass loss compared to the WL glaciers.

Despite the presence of a high debris cover on the WL glaciers, typically known to reduce melt, their low-angle debris-covered tongues experienced a greater impact from warming. Additionally, surface slope emerged as a critical factor influencing mass balance variability in both sub-regions, explaining 30–42% of the observed variance. This study underscores the dominance of nonclimatic factors, specifically topographic control, in shaping overall glacier mass balance, particularly in the EL.

The presented mass balance estimates offer the most recent insights into glacier mass changes in the region, serving as inputs for calibration and validation of individual glacier or region-wide glacio-hydrological modelling experiments. Future research should delve into understanding the complex response of debris-covered glaciers, especially given their prevalence in the WL areas and the broader Himalayan region. Utilizing high-resolution DEMs, such as HMA 8-m and Pléiades, will be crucial to discern relative mass loss contributions from debris-covered parts, ice cliffs and supraglacial pond areas. This approach will provide a more comprehensive understanding of how clean-ice and debris-covered glaciers distinctly respond to climate change.