Spatially heterogeneous glacier elevation change in the Jankar Chhu Watershed, Lahaul Himalaya, India derived using ASTER DEMs

Abstract This study investigates elevation change (dh) and geodetic mass budget of glaciers in the Jankar Chhu Watershed (JCW), Lahaul Himalaya, India, based on the difference in ASTER DEMs during 2002–2018. Glacier-wide spatially heterogeneous dh patterns were evaluated in the context of morphological and topographical settings. During 2002–2018, glaciers show a mean annual elevation change rate (dh/dt) of −0.38 ± 0.10 m a−1, resulting in a specific mass budget of −0.32 ± 0.09 m w.e.a−1, close to the previously reported estimates in western Himalaya. Nearly stagnant thick debris-covered tongue (>10% debris cover) characterized by melt hotspots exhibits maximum dh at up-glacier instead of the terminus. Debris-free glaciers (<10% debris cover) show maximum dh near the terminus. Spatially heterogeneous dh under varying debris cover is interpreted as an insulating effect of debris thickness as validated by field measurements. We suggest that elevation change of debris-covered glaciers cannot be generalized and glacier-wide spatially detailed mapping of dh is necessary to better understand the control of different surface morphology under warming climatic conditions in the western Himalayas. Highlights We present a spatially heterogeneous elevation change rate of 33 glaciers (150.9 km2) in the Jankar Chhu Watershed (JCW), Chandrabhaga basin, Western Himalaya based on the differences of ASTER digital elevation models (DEMs). We estimate the altitude-dependent elevation change of glaciers under varying debris thickness and surface morphology. Glaciers in the JCW show mean dh/dt of –0.38 ± 0.10 m a−1, resulting in a specific mass budget of –0.32 ± 0.09 m w.e.a−1, close to the previously reported estimates in western Himalaya. Debris-covered ice exhibits higher elevation change than debris-free ones. Thick debris-covered tongues show maximum surface lowering at up glaciers while debris-free tongue exhibits maximum lowering near the terminus. Debris thickness controls the altitude-wise spatial elevation change pattern in the JCW.


Introduction
Glacier elevation changes and mass budget data for the Himalayan region are available either from annual field measurements on a limited number of glaciers using the direct glaciological method (Dobhal et al. 2013;Azam et al. 2014) or from the comparison of the glacier surface topography of different years from digital elevation models (DEMs) followed by a density assumption for converting volume to mass change (K€ a€ ab et al. 2012;Gardelle et al. 2013;Brun et al. 2017;Shean et al. 2020;Hugonnet et al. 2021). Together with available DEMs, and the fact that inaccessible mountainous areas and entire glacier systems can be measured, the so-called geodetic method has become a popular approach to derive elevation and mass changes for a large number of glaciers at regional and global scales (Gardelle et al. 2013;Fischer et al. 2015;Hugonnet et al. 2021). In recent years, several DEMs from different sensors and satellite platforms have been deployed to determine glacier surface elevation and mass change in the Hindukush-Karakorum-Himalaya (HKH) region (K€ a€ ab et al. 2012;Gardelle et al. 2013;Shean et al. 2016;Hugonnet et al. 2021).
Based on the geodetic method, region-wide heterogeneous elevation change and mass loss have been reported throughout the HKH region during the last two decades (K€ a€ ab et al. 2012;Gardelle et al. 2013;Brun et al. 2017;Shean et al. 2020;Hugonnet et al. 2021). Karakorum region shows a slight mass loss or balanced condition (Brun et al. 2017; although the mass gain is more concentrated in the Kunlun Shan (Brun et al. 2017). Region-wide mass loss is higher in the Nyainqentanglha range (À0.62 ± 0.23 m w.e.a À1 ) in the eastern Himalayan glaciers compared to Lahaul Spiti region (À0.37 ± 0.09 m w.e.a À1 ) in the western part (Brun et al. 2017). Using ice, cloud and land elevation satellite (ICESat) laser altimetry and SRTM DEM, K€ a€ ab et al. (2012) reported almost balanced mass budgets for the Karakoram region during 2003-2008, while for surrounding regions, including the Lahaul-Spiti region, higher surface lowering was observed (Vijay and Braun 2016;Mukherjee et al. 2018;Vijay and Braun 2018;Abdullah et al. 2020). Heterogenous glacier mass loss has been reported in the Chandrabhaga basin of the western Himalayas (Berthier et al. 2007;Vijay and Braun 2016;Mukherjee et al. 2018;Vijay and Braun 2018). Differential dh/dt of debris-covered and debris-free glaciers was also reported for the Himalayan region (K€ a€ ab et al. 2012;Maurer et al. 2019). Previous studies noted much more spatially heterogeneous surface elevation change rates, even within a single basin, which is controlled by glacier morphological and topographic settings (Brun et al. 2019), supraglacial debris thickness and distribution Anderson and Anderson 2016;Salerno et al. 2017;Tsutaki et al. 2019;Bisset et al. 2020), supraglacial water bodies and ice cliffs (Benn et al. 2012;Ragettli et al. 2016;King et al. 2017;Buri et al. 2021) and emergence surface velocity (Anderson and Anderson 2016;Brun et al. 2018). In comparison with debris-free (clean) ice, the melt may enhance when the surface is covered with a thin layer of debris (1-2 cm) as a result of increased absorption of solar radiation and related heat transfer Fyffe et al. 2020). On the other hand, thick debris cover reduces ice melt rates as less surface heat will be conducted through the debris layer and transferred to the ice (Brock et al. 2010;Lambrecht et al. 2011;Vincent et al. 2016). Although  reported that thin debris layers do not enhance the surface melting of Karakorum glaciers in the lower altitudes. Studies have reported higher surface lowering on thick debris-covered ice compared to debris-free ice in the Himalayas (K€ a€ ab et al. 2012;Gardelle et al. 2013;Maurer et al. 2019). This characteristic of debris-covered ice in the Himalayan region is referred to as a 'debris-covered anomaly' (Salerno et al. 2017;Bisset et al. 2020). Similarly, supraglacial water bodies and ice cliffs act as melt hotspots enhancing ice melting (Sakai et al. 1998(Sakai et al. , 2000(Sakai et al. , 2002Benn et al. 2012;Reid and Brock 2014;Steiner et al. 2015;Buri et al. 2016;Miles et al. 2016Miles et al. , 2018Anderson et al. 2021). So, ascertaining the possible role of these factors on glacier mass loss at the basin scale will help better understand the heterogeneity of glacier dynamics, calibration of hydrological and glaciological models, and assessment of meltwater contribution to river discharge.
The outline above demonstrates that detailed studies on the influence of different surface features (supraglacial ponds, ice cliffs and debris thickness) on glacier elevation change and glacier-wide updated and improved estimates are required to evaluate the possible reasons for the heterogeneous elevation change pattern. Here, an attempt has been made to investigate the basin-wide and glacier-wide elevation change and mass budget of glaciers in the JCW of the Chenab basin in the Lahaul Himalayan region. Earlier studies on elevation change and mass budget in the Lahaul Himalayan region are mainly confined to the upper Chandrabhaga basin (Vijay and Braun 2016;Mukherjee et al. 2018).
To the best of our knowledge, no studies have investigated the elevation change of glaciers in the JCW based on remote sensing data and field measurements associated with several other topographic and morphological parameters. Hence, the aims of the study are: (1) to access basin-wide and glacier-wide elevation change and geodetic mass budget; (2) to investigate the glacier-wide spatial elevation change pattern under varying debris thickness and surface morphology; and (3) to examine the possible reasons of heterogeneous elevation change pattern.

Study area
The JCW (32 38'43.00"N-32 59'55.25"N and 76 56'41.17"E-77 14'35.48"E) lies in the western Himalayan region of the Lahaul Spiti district of the Himachal Pradesh state in India (Figure 1(a)). The Jankar Chhu (small river) is a tributary of the Bhaga River ( Figure 1b). Bhaga River confluence at Chandra forms Chandrabhaga, which ultimately contributes to the Indus River system. The trunk stream Jankar Chhu originates after the confluence of Dali and Bagrari Chhu. This watershed has an altitude ranging from $3300 to 6100 m a.s.l (above sea level). The JCW is sandwiched between the Pir-Panjal range in the south and the Zanskar range in the north Sharma 2019, 2021). This part of the western Himalayas falls under the monsoon-arid transition zone. The climate of the study area is dominated by the mid-latitude westerlies in winter and the south Asian monsoon in the summer (Das and Sharma 2019). However, mid-latitude westerlies contribute more than $80% of annual precipitation in this region (Vijay and Braun 2016;Das and Sharma 2019).
The JCW has 153 glaciers (>0.02 km 2 ) covering an area of 185.6 ± 3.8 km 2 ($27% of the total watershed area) as of 2016 (Das and Sharma 2019). Dali, Mayar and Bagrari are the largest glaciers in the JCW (Figure 1(b)). In this study, a total of 33 glaciers covering an area of 150.9 km 2 (81% of the glacierized area) were selected based on heterogeneous characteristics (morphological and topographical) and satisfying the following criteria: (1) minimum size of >1 km, 2 that lies entirely within the delineated catchment boundary; (2) Figure 1. (a) Divisions and major ranges of the Hindukush-Karakorum-Himalaya (HKH) region. White dot shows the location of the Jankar Chhu Watershed (JCW) in the western Himalayan region. Blue lines show the major river system in the HKH region. Background is a shaded relief map. (b) Distribution of studied glaciers with the respective ID in the JCW has overlaid SRTM DEM with 60% opacity. Stars denote the field-visited glaciers. Red and blue colours show the distribution of debris-covered and debris-free ice. Glacier outlines (2016) were obtained from the study of Das and Sharma (2019).
valley type glaciers with well-defined accumulation and ablation zone with varying debris thickness, presence of supra and proglacial lakes; (3) higher elevation range; and (4) availability of void-free DEMs.

Data sources
ASTER scenes of 28 October 2002 (late ablation season) and 01 May 2018 (early ablation season) were acquired from www.earthdata.nasa.gov. No cloud-free scenes are available for the late ablation season during 2018-present day. Several field visits were carried out in the JCW between 2015 and 2020. Elevation information was recorded using the Garmin Etrex GPS instrument for elevation change accuracy assessment on ice. Four glaciers (Dali, Mayar I, Mayar II and Bagrari glacier) were monitored during the field visits based on their heterogeneous surface morphology. Snout characteristics (position, shifting and presence of lakes) and surface morphology (ponds, ice cliffs, debris thickness) were mapped and measured using handheld GPS (Garmin Etrex; ±2 m accuracy). Ice temperature was measured using a handheld laser Infrared Thermometer under varying debris thickness at Bagrari glacier ( Figure 2). To analyse the spatially heterogeneous thinning patterns of the glacier, several other datasets were used (Supplementary Table 1). Surface area change rates were taken from the study of Das and Sharma (2019). Glacier velocity datasets were taken from the study of Das and Sharma (2021). Mean surface velocity between 1992 and 2020 was used. Ice thickness datasets were obtained from Farinotti et al. (2019). Debris thickness and sub-debris melt factors data were taken from Rounce et al. (2021).

ASTER DEM processing
The available global ASTER DEM is generated by NASA with the SilcAst software (http:// www.silc.co.jp/en/products.html) and is packaged in the AST14DMO product (Fujisada   Abrams et al. 2020). Recent studies have reported that ASTER14DMO DEMs are plagued with a high-frequency satellite jitter that induces altitude perturbations in both across-track and along-track directions (Girod et al. 2017;Hugonnet et al. 2021) resulting in lower vertical error/accuracy in a mountainous region (Fujisada et al. 2005). Previous studies have reported a vertical accuracy of ±20 m for AST14DMO products (Hirano et al. 2003;Fujisada et al. 2005;Berthier et al. 2007). However, sufficient geometric quality in DEMs is necessary for glacier volume change estimation where the expected elevation change is smaller than the accuracy of the product.
ASTER DEMs were generated from the stereo bands (3N and 3B) using the MicMac ASTER tool with a spatial resolution of 30 m. Girod et al. (2017) developed a DEM processing tool implemented in MicMac photogrammetric processing library for bias-corrected ASTER DEM generation. The MicMac tool helps to overcome the jitter effects (across and along track) in the ASTER stereo pairs, which significantly improves the altitude perturbations. Across-track jitter (roll) produces a large offset in the epipolar line and cannot be corrected using the correlation algorithm (Girod et al. 2017;. Whereas, along-track jitter (pitch) produces systematic bias (horizontal shift and altitudinal biases) and can be corrected using 3-D co-registration of DEMs over stable terrain . In this study, two across-track jitters corrected ASTER DEMs of 2002 and 2018, generated using the MicMac tool were obtained from the study of Hugonnet et al. (2021) at 30 m spatial resolution (Supplementary Table 1). Detailed methodology of ASTER DEM generation using the MicMac tool is presented in Girod et al. (2017). MicMac-generated DEMs were further processed for additional bias (horizontal shift, altitude and curvature dependent) removal and discussed in the next section.

DEM postprocessing
At first, off-glaciers (stable terrain) areas were obtained based on the criteria suggested by the previous studies (Berthier et al. 2007;Nuth and K€ a€ ab 2011;Gardelle et al. 2013;Berthier et al. 2014Berthier et al. , 2016. Glaciated areas were erased using GLIMS glacier outlines. Terrain characterized by slope <4 and >45 were excluded because DEM derived from the optical stereo images (ASTER) are less precise over steeper terrain (i.e. >45 ) (Berthier et al. 2007. On the other hand, nearly flat terrain (<4 slope) does not constrain the retrieval of mean horizontal shift between reference and target DEMs . In this study, a robust analytical 3-D co-registration method proposed by Nuth and K€ a€ ab (2011) was used to remove horizontal and vertical shifts based on the elevation difference in stable terrain (Supplementary Figure S1). The elevation differences along and across the satellite track were computed on stable areas (Supplementary Figure S2), and when necessary, the biases were corrected using a fifthorder polynomial fitting (Nuth and K€ a€ ab 2011). A systematic bias of ASTER DEM as a function of altitude has been reported in the European Alps and Himalayan region (Gardelle et al., 2013), where ASTER elevations were shown to be underestimated on ridges/summit and overestimated in the valleys. To detect the occurrence of similar biases in the JCW, we compute the mean elevation difference of two ASTER DEMs in the offglaciers region for each curvature and altitude range and then plot these differences as a function of altitude (Supplementary Figures S3-S4). Extreme biases occur in the lower (<3800 m) and higher elevations (>5800 m). These biases may be linked to the roughness of the topography. Thus, these apparent elevation biases are corrected based on the relationship between elevation difference and terrain curvature over stable terrains. We performed fifth-order polynomial fitting to minimize these biases. After the biases between the DEMs have been identified and corrected (Supplementary Figure S5), glacial elevation changes can be computed and analysed.
In the high mountainous regions, data gap in the images mainly occur due to cloud cover, cast shadow and the zone of fresh snow in the accumulation region (Pieczonka and Bolch 2015), hence needs to be excluded before proceeding for mass balance (MB) calculations. Outliers for the stable terrain were defined by 1.5-fold of the interquartile range (Pieczonka and Bolch, 2015). Outliers in the glacierized regions were excluded in three different steps: (1) zones with large outliers (DEM pixels for which elevation difference exceeded greater than 100 m) were excluded. This threshold value was chosen based on the previously reported elevation change measurements in this part of the Himalayas (Berthier et al. 2007;K€ a€ ab et al. 2012;Gardelle et al. 2013;Vijay and Braun 2016;Mukherjee et al. 2018); (2) altitude bins where the number of pixels is <60% were excluded for change analysis. Where no elevation change was available for an altitude bin, a mean elevation change value of its neighbouring altitude bands was assigned to assess the MB over the whole glacier area; (3) for each altitude bin where mean elevation change exceeds ±3 standard deviations were excluded. The resulting voids due to outlier removals were filled using bilinear interpolation (K€ a€ ab 2008; Zheng et al. 2018).

Altitude-dependent elevation change and mass budget
After the biases between master and reference DEM were identified and corrected, glacier elevation change was computed by subtracting the newer DEM (2018) from the older one (2002). Elevation changes were analysed for 25 m altitude bins. Within each altitude band, the percentage of valid pixels available for elevation change was computed. Mean dh and standard deviation for each altitude bin were calculated. The average dh of all altitude bands was considered for glacier-wide elevation change measurements.
To calculate the geodetic mass budget of the glaciers, the surface elevation changes of valid pixels within the glaciated area were averaged. First, the volume change was calculated by integrating the mean elevation change values with the corresponding glacier area for individually analysed glaciers. Later, for the volume-to-mass conversion, a constant ice density (q ice Þ of 850 ± 60 kg m À3 was used (Huss 2013). The geodetic mass budget was calculated using glacier outlines of 2000 (Das and Sharma 2019). The final geodetic MB was calculated using the following formula: where A surf is the total glacier area for a given glacier and q w is the water density (999.972 kg m À3 ).

Debris categorization
Studied glaciers were categorized into two groups (debris-covered glaciers and debris-free glaciers) based on their percentage of debris cover area to total ice area and surface morphological characteristics (melt hotspot, elevation and slope). Melt hotspots (ponds, ice cliffs and supraglacial channels) on the glacier surface were mapped for four glaciers using 3D images in Google Earth and further corrected using field photographs (Anderson et al. 2021 Ali et al. (2017) proposed three fractions: clean glaciers with a debris-cover fraction <25%; sparsely debris-covered glaciers with a debris-cover fraction between !25% and 50% and debris-covered glaciers with debris-cover fraction !50%. It is noteworthy that this distinction is arbitrary, as there is a continuum between debris-covered and debris-free glaciers. Most of the glaciers in the JCW are debris free and only 11% of total glacierized are covered by debris as of 2016 (Das and Sharma 2019). For debris categorization, we used the following two criteria (1) an optimal threshold fraction of 10% of debris cover, and (2) spatial distribution of debris in the ablation zone of the glacier. The homogeneous distribution of debris in the ablation zone was given priority rather than the distribution of debris along the valley edges. In addition, we analyse modelled debris thickness pattern on the glacier from the study of Rounce et al. (2021) and presented in supplementary Figure S6. We assign a debris thickness of <5 cm as a thin debris layer and >5 cm as a thick debris layer ).

Uncertainty estimations
4.5.1. Accuracy assessment of dh based on stable terrain Rigorous statistical analysis was performed based on elevation differences in the off-glacier areas (stable terrain) (Table 1). More than 20,000 random points from off-glacier regions were selected from the elevation difference map. The performance analysis indices are retained as follows: 1. The mean and standard deviation of off-glacier regions (Berthier et al. 2014); 2. The median and normalized median absolute deviation (NMAD) were employed as they are not sensitive to mean and standard deviation (Berthier et al. 2007; H€ ohle and H€ ohle 2009; Berthier et al. 2014;Vijay and Braun 2016). The NMAD can be expressed using the following equation: where Dh j is denotes the individual errors and median Dh is the median of the errors. Statistics of stable terrain elevation differences are shown in Table 1.

Uncertainty in elevation change and mass budget
In this study, surface elevation change uncertainty mainly arises due to the uncertainty of DEM difference (U dh ) and the uncertainty of the glacier area change (U Dsurf ). The uncertainty measurement of U dh was calculated using Equations (3-5) according to (Gardelle et al. 2013): and, where U dh i ð Þ is the uncertainty corresponding to the ith altitude band, d dh ðiÞ is the standard deviation of the mean elevation change of the stable terrain of the ith altitude band, N totðiÞ is the total number of pixels in the ith altitude band, dh PR is the pixel resolution of the elevation difference image (30 m in this study) and d is the spatial autocorrelation distance (in meters) of elevation change in stable terrain. We used a spatial autocorrelation distance of 400 ± 50 m. Glacier mapping uncertainty was calculated based on the buffer method (Granshaw and Fountain 2006). A buffer size of 5 m was used for individual glaciers. The overall uncertainty of elevation changes (U EC ) was calculated using the following formula of standard error propagation: The uncertainty of density of ice (U q ice ) was considered as ±60 kg m À3 for the mass budget measurement (Huss 2013). Finally, the mass budget uncertainty (U M ) was calculated using the following equation.
where D h is the elevation difference, and t is the time period in years.

Off-glacier elevation difference
Off-glacier statistics show the accuracy of the final elevation difference map in the JCW (Figure 3). Any apparent motion at stable terrain is due to image misregistration and distortion produced in the orthorectification process. The apparent thickness change in the stable terrain shows approximately normally distributed (Figure 3(c,d)). The spread of these normal distributions, as measured by their standard deviations (r), provides an estimate of random error in the DEM difference map. Off-glacier mean and standard deviation in the final elevation difference map is À0.02 and 0.93 m, respectively (Figure 3(c,d)).

On-glacier elevation difference
A total of 1247 GPS points were collected from the terminus zone of three glaciers (Mayar I, Mayar II and Bagrari). On-glacier average elevation change rate between GPS elevation (2018) and ASTER DEM (2002) is À2.16 ± 0.32 m a À1 (Figure 4). Terminus zone of Mayar I glacier experienced dh/dt of À2.35 ± 0.40 m a À1 (Figure 4(a)). Debris-covered ice experienced slightly higher (-2.26 ± 0.26 m a À1 ) dh/dt compared to the debrisfree part (-1.94 ± 0.37 m a À1 ) during the observation period. Mayar II glacier experienced a mean dh/dt of À1.85 ± 0.19 m a À1 (Figure 4(b)). Thick debris-covered tongue of the Bagrari glacier showed a mean dh/dt of À2.16 ± 0.38 m a À1 based on the difference between GPS (2018) and ASTER 2002 DEM difference (Figure 4(c)). Based on the same GPS points, dh/dt was extracted from the elevation difference map. Maryar I, Mayar II and Bagrari glaciers show a mean dh/dt of À2.37 ± 0.36, À1.84 ± 0.27 and À1.66 ± 0.24 m a À1 , respectively, during 2002-2018. Spatial distribution of the point-based elevation change map shows that ice covered by a thin debris layer (<5 cm) near the terminus experienced maximum dh/dt compared to ice under a thick debris layer (>5 cm), although debris morphology (boulder, sediments) needs to be accounted separately (Figure 4(c)).

Mayar I glacier.
Mayar I is the longest glacier (9.53 km), covering an area of 15.15 km 2 in 2000. The altitude of Mayar I glacier ranges from 4580 to 5950 m a.s.l. A total of 1.25 km 2 (7.93% of total ice area) of the surface area is covered by debris. Debris- covered areas are mainly confined to the glacier tongue edges (Supplementary Figure S9). Mean debris thickness ranges between 0.01 and 0.73 m. Mayar I glacier experiences a maximum dh/dt of À1.97 ± 0.45 m a À1 near the terminus, followed by a continuous decrease in dh/dt (Figures 7(d) and 8(e)). Near the line of equilibrium ($5250 m a.s.l), dh/dt reduces to À0.08 ± 0.1 m a À1 (Figure 8(d)). The altitude zone of 4580-4600 m a.s.l shows a mean debris thickness of 0.51 ± 0.24 m and dh/dt of À1.16 ± 0.23 m a À1 . The following altitudinal zone (4600-4650 m a.s.l) shows a higher dh/dt of À2.28 ± 0.45 m a À1 and a debris thickness of 0.05 ± 0.01 m. Mayar I glacier is nearly stagnant (4.82 m a À1 ) in the lower ablation zone (4580-4900 m a.s.l) with a mean ice thickness of 55.83 ± 8.56 m. Ice thickness is maximum (154.57 ± 32.54 m) between 5000 and 5200 m a.s.l ( Figures  7(d-f) and 8(e-h)).

Mayar II glacier.
Mayar II glacier covered an area of 12.97 km 2 in 2000. Only 3.64% of the total ice area is covered by debris. Mean debris thickness of Mayar II glacier ranges between 0.12 and 3.14 m. The mean slope of the glacier is 14 with a majority of the ice area oriented in south-southeast directions. Part of the Mayar II glacier terminates at a proglacial lake (Supplementary Figure S10). Mayar II glacier experienced a dh/dt of À1.07 ± 0.21 m a À1 near the terminus (4790-4850 m a.s.l) where the mean debris thickness is 2.64 m (Figure 7(g-i)). The zone of comparatively thin debris layer (4900-5000 m a.s.l; 0.53 m) shows a much higher dh/dt of À1.88 ± 0.41 m a À1 . The mean velocity of lower ablation is 3.63 m a À1 with a mean ice thickness of $80.6 m (Figures 7(g-i) and 8(i-l)).

Bagrari glacier.
Bagrari glacier covers a surface area of 13.73 km 2 (2000). Out of the total ice area, 16.75% of the area is covered by debris in the ablation zone. The altitude of Bagrari glacier ranges from 4590 to 6050 m a.s.l, with a mean surface slope of 15.35 . The nearly flat (<5 slope) lower debris-covered ablation zone of the Bagrari glacier is characterized by melt hotspots (Supplementary Figure S11). Near terminus (between 4590 and 4700 m a.s.l), Bagrari glacier experienced a comparatively lower dh/dt of À1.12 ± 0.32 m a À1 compared to the up-glacier zone (4700-4900 m a.s.l; À1.60 ± 0.36 m a À1 ) (Figure 7(j)). The 4590-4700 m a.sl zone shows the highest mean debris thickness of 2.82 m and the lowest velocity of 2.96 m a À1 with a mean ice thickness of 25.65 m (Figures 7(j-l) and 8(m-p)). In contrast, the 4700-4900 m a.s.l zone shows a debris thickness of 0.46 m with a mean velocity of 4.53 m a À1 . This shows that a zone of higher debris thickness is associated with lower elevation change and vice versa.

Reported elevation change and mass budget in Western Himalaya
Several studies have reported the geodetic MB of glaciers in high-mountain Asia (K€ a€ ab et al. 2012;Gardelle et al. 2013;Brun et al. 2017;Zhou et al. 2018;Maurer et al. 2019;Shean et al. 2020;Hugonnet et al. 2021) as well as in the western Himalaya (Berthier et al. 2007;Vijay and Braun 2016;Bhushan et al. 2018;Mukherjee et al. 2018;Vijay and Braun 2018) at various spatial and temporal scale. The detailed assessment of the previously reported elevation change rate and associated mass budget concentrating in the western Himalayan region are given in Supplementary Table S3. Berthier et al. (2007) reported a geodetic MB of À0.7 to À0.85 m w.e.a À1 (ice density of 900 kg m À3 for the glacier material below the firn line and 600 kg m À3 above the firn line) in the accumulation zone based on SRTM and SPOT5 datasets (1999)(2000)(2001)(2002)(2003)(2004) in the Lahaul-Spiti region. Higher mass loss of À1.31 and À1.12 m w.e.a À1 have been reported for Bara Shigri and Chhota Shigri glaciers (Berthier et al. 2007). During 2003-2009, K€ a€ ab et al. (2012 measured geodetic mass budget of À0.32 ± 0.06 m w.e.a À1 for the western Himalayan region based on the ICESat and SRTM DEM (Supplementary Table S3). The present mass budget estimate (-0.32 ± 0.08 m w.e.a À1 ) is similar to their findings. Gardelle et al. (2013) estimated a geodetic MB of À0.45 ± 0.14 m w.e.a À1 (an ice density of 850 kg m À3 ) for the glaciers in the Lahaul-Spiti region during 1999-2011. The present study shows a lower mass budget of À0.32 ± 0.08 m w.e.a À1 . Based on the difference in ASTER DEM during 2000-2016, Brun et al. (2017) reported a mass budget of À0.38 ± 0.08 m w.e.a À1 , close to our estimation (-0.32 ± 0.10 m w.e.a À1 ). Vijay and Braun (2016) reported a mass budget of À0.53 ± 0.37 m w.e.a À1 (ice density of 850 kg m À3 ) between 2000 (SRTM) and 2012 (TanDEM-X) for $1700 km 2 glacierized area in the Lahaul Spiti region (mainly confined to the upper Chandra basin). The present analysis shows a much lower mass budget (-0.32 ± 0.10 m w.e.a À1 ) based on the differences of ASTER DEM (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) for the 155 km 2 of the glacierized area in the JCW. Lower estimates of mass budget in the present study could be attributed to the differences in the studied area and heterogeneous topographical and morphological settings. Another study by Vijay and Braun (2018) Table S3). The above findings demonstrate that glaciers in the JCW have experienced similar mass loss as reported for the different basins in the western Himalayas during the last two decades, although at the regional scale mass budget varies significantly.
6.2. Heterogeneous elevation change and possible reasons 6.2.1. Differential lowering of debris-covered and debris-free ice The dh/dt of debris-covered ice is higher (-0.79 ± 0.23 m a À1 ) than clean ice (-0.57 ± 0.12 m a À1 ) for the studied glaciers in the JCW. K€ a€ ab et al. (2012) reported a similar dh/dt for the debris-covered and debris-free ice in the Himalayas during 2003-2008. Gardelle et al. (2013) also showed a similar dh/dt pattern in the Lahaul region during 1999-2011. Vijay and Braun (2016) also reported a higher down-wasting rate for debris-covered than debris-free ice in the Lahul Spiti region. Maurer et al. (2019) reported similar observations for glaciers in the Bhutan Himalayas. This characteristic of debriscovered ice is referred to as a 'debris-covered anomaly' in the Himalayan region (Salerno et al. 2017;Bisset et al. 2020). Thick debris-covered glacier tongues consist of melt hotspots (supraglacial ponds, ice cliffs) in the JCW and show maximum lowering (Supplementary Figures S8-11). Higher dh/dt on debris-covered ice is attributed to the mechanisms of heat transfer through melt hotspots reported by previous studies (Vijay and Braun 2016;Anderson and Anderson 2018;Bhushan et al. 2018;Maurer et al. 2019;Anderson et al. 2021). In classical theory, thick debris cover reduces melt (Anderson and Anderson 2016), while the presence of ice cliffs and supraglacial ponds can foster ice melts on debris-covered glaciers (Sakai et al. 1998(Sakai et al. , 2000K€ a€ ab et al. 2012;Pellicciotti et al. 2015;Steiner et al. 2015;Bhushan et al. 2018;Tsutaki et al. 2019;Anderson et al. 2021;Buri et al. 2021). Previous studies have reported that the supraglacial ponds act as a temperature storehouse and enhance the melting of underneath ice (Tsutaki et al. 2019;Anderson et al. 2021). Ice cliffs significantly increase the melt rate, even though it occupies only a tiny portion of the total debris-covered ablation zone (Reid and Brock 2014;Steiner et al. 2015;Buri et al. 2016Buri et al. , 2021. Further, orientation plays a significant role in surviving ice cliffs on debris-covered ice . Based on the numerical modelling, Buri and Pellicciotti (2018) calculated that south-facing ice cliffs disappear within a few months while north-facing ones sustain a longer time and act as a stable contributor to melting debris-covered glaciers. This study quantifies ice cliffs and melts ponds based on the high-resolution 3-dimensional Google Earth images and field measurements for Bagrari, Dali and Mayar II glaciers (Supplementary Figures S8-S11). The result shows that the spatial distribution of ice cliffs and supraglacial ponds corresponds with the high dh/dt for the studied glaciers, although morphology and dynamics of ice cliffs can vary from seasonal to annual scale Kneib, Miles, Jola, et al. 2021). Ice cliffs on the Bagrari glacier are covered by a very thin layer of dust particles (Figure 9). Field measurements on Dali and Bagrari glaciers in Lahaul Himalaya revealed the presence of very thick debris cover exceeding >2 m at some places ( Figure  9). Apart from that, different debris types (giant boulders, fine-grain sediments) should not be generalized as heat absorption capacity is different for varying materials (Figure 9). Nearly flat down glacier tongues with very low surface velocity are more prone to have thicker debris than the steep up-glacier where steep rock wall contributes debris and transport to the lower elevation due to gravity . Earlier studies have reported that nearly stagnant debris-covered tongues with gentle slopes are more prone to the growth of melt ponds in the Himalayas (Thompson et al. 2016;Irvine-Fynn et al. 2017). Further, some glaciers are characterized by a large proglacial lake (Mayar II) that further supports the melting and enhances surface lowering at the terminus (Supplementary Figure S9). For example, part of Mayar II glacier is terminating at large proglacial lakes (Supplementary Figure S8) and experienced a higher dh/dt of À2.4 ± 0.35 m a À1 between snout ($4250 m) and 4700 m altitude, while above 4700 m thinning rate is considerably low (-0.32 ± 0.18 m a À1 ). Vijay and Braun (2016) reported a similar surface lowing trend for lake-terminating glaciers in the Lahaul Spiti region. Besides, Maurer et al. (2019) reported higher thinning for lake-terminating glaciers across the Himalayan region.

Spatial elevation change patterns
We observed two separate altitude-wise elevation change patterns for debris-covered and debris-free glacier tongues in the JCW (Figure 10). Thick debris-covered glacier tongue shows maximum surface lowering up-glacier instead the terminus (referred to as thinning pattern I). Bagrari glacier shows this kind of thinning pattern (Figure 10(a)). The insulating effect of a thick debris cover is considered to be the probable reason for the reduced melt near the terminus. Down-glacier areas are prone to have thicker debris compared to up-glacier because of avalanche-driven debris accumulation from the potential debris supply slope (Nagai et al. 2013;Das and Sharma 2019) and the emergence of englacial debris (Boulton (f-h) Different shapes, types, and morphology of ice cliffs. Exposed bare ice is covered with thick dust particles which act as an insulator of heat. Figure 10. Elevation change rate (dh/dt) and debris thickness patterns of debris-covered and debris-free glaciers in the JCW. (a) Debris thickness distribution in the ablation zone of Bagrari glacier overlaid on Google Earth images and altitude-wise elevation change rate. (b) Debris thickness distribution and altitude-wise elevation change of Mayar I glacier. Green bar shows the zone of maximum thinning. Note: Thick debris-covered tongue shows maximum thing upglacier instead of the terminus (referred to as pattern I). In contrast, debris-free tongue shows maximum thinning near the terminus (referred to as pattern II).
1970). Consequently, maximum surface lowering takes place at higher altitudes, where either comparably thinner debris cover or supraglacial ponds and lakes, and ice cliffs trigger lowering (Ragettli et al. 2016). Dobhal et al. (2013) also observed nearly four times higher melt Figure 11. Ice temperature under varying debris thickness at the terminus region of the Bagrari glacier. Debris thickness and underneath ice temperature were measured using a measuring tape and infrared laser thermometer, respectively. rates at 4400 m compared to the terminus ($3900 m) of Chorabari glacier in the Indian Himalayas and linked this to the insulating effects of thick debris cover at lower altitudes. Some of the debris-free glaciers show maximum surface lowering near the terminus (referred to as thinning pattern II). This type of surface lowering is explained by the higher air temperatures at lower elevations and lower albedo of the bare ice compared to the snow and firn-covered areas. Mayar I glacier shows this kind of thinning pattern (Figure 10(b)). Visual appearance in the high-resolution Google Earth image and field photographs of Mayar I glacier shows that the lower ablation zone (4475-4775 m a.s.l) is characterized by a bare ice surface with numerous supraglacial streams and very thin debris layers (Supplementary Figure S9). Field investigations show the ice surface is dispersedly covered by dust particles and large boulders in several places (Supplementary Figure S9). Elevation change starts to increase near the Equilibrium Line Altitude (ELA) at $5350 m for Mayar I glacier. Latent heat transfer on the thin sediment layer is more effective, which starts appearing after the line of equilibrium. This supports the general theory that the thin debris layer absorbs more heat and transfers to the underneath ice, increasing the melting and vice versa. This observation is supported by the field measurements of ice temperature under varying debris thickness for the Bagrari glacier using an Infrared laser thermometer and measuring tape ( Figure 11). We observed that the ice temperature decreases as debris thickness increases ( Figure 11). Altitude-dependent dh/dt also reveals that the altitude band having higher debris thickness experiences a lower thinning rate and vice versa. We observed a negative relationship (À0.4) between debris thickness and elevation change, revealing that as debris thickness increases surface lowering decreases and vice versa. Studies in the Himalayas have reported a higher elevation change under a thin debris layer based on field measurements (Dobhal et al. 2013;Pratap et al. 2015).

Conclusions
This study presents basin-wide and glacier-wide elevation change and geodetic mass budget in the JCW, Lahaul Himalaya at spatial and temporal scales. The geodetic mass budget was estimated for 33 glaciers (>1 km 2 in size), covering $155 km 2 of glacierized area. During 2002-2018, the glaciers underwent a mean annual surface elevation change of À0.38 ± 0.10 m a À1 , resulting in a specific mass budget of À0.32 ± 0.08 m w.e. a À1 , close to the previously reported rates in Lahaul Spiti region. The rate of surface lowering of debris-covered ice was higher (-0.78 ± 0.35 m a À1 ) than clean ice (-0.58 ± 0.23 m a À1 ). Thick debris-covered tongue characterized by melt hotspots (supraglacial ponds and ice cliffs) shows maximum elevation change at up glacier instead of the terminus. Debris-free glaciers show maximum elevation change near the terminus, possibly due to the higher air temperature in the lower altitudes. This heterogeneity in glacier surface lowering is controlled by the interplay of several local factors such as slope, thickness of the debris, elevation, and topographic settings. The higher elevation change rate of debris-covered ice is attributed to nearly stagnant surface movement coupled with the presence of melt hotspots, as reported for the other part of the Himalayas. Future works should be concentrated on field-based precise debris thickness measurements, mapping, and monitoring of ice cliffs and supraglacial ponds.
New Delhi (3090/(NET-DEC.2014) for financial support during field visits in the Himalayas. We thank Prof. Bradley C. Rundquist and three anonymous reviewers for their insightful comments and suggestions which considerably improved the earlier draft. First author acknowledges Marin Kneib for his valuable suggestions which significantly improved the figure quality and readability of the draft.