Response of the Thick and Thin Debris-Covered Glaciers between 1971 and 2019 in Ladakh Himalaya, India—A Case Study from Pensilungpa and Durung-Drung Glaciers

: This paper aims to broadly understand the response of glaciers to thick and thin debris cover from one of the less explored regions (Zanskar) of the Himalaya. The present study is based on ground-based measurements (from 2015 to 2019), satellite data (since 1971), and available topographic maps (at a 1:50,000 scale). The study includes snout retreat, changes in equilibrium line altitude (ELA), surface elevation, and modeled mass balance of thick and thin debris-covered Pen-silungpa (Suru River basin) and Durung-Drung (Doda River basin) glaciers in the western Indian Himalaya, Ladakh, for the past five decades. The Durung-Drung Glacier (DDG) receded ~−624 ± 547 m with an average rate of −12 ± 11 m a − 1 between 1971 and 2019. The frontal part of the DDG is broad (~2 km wide), which shows wide discrepancies in its retreat. Compared to DDG, the small and narrow snout of the Pensilungpa Glacier (PG) retreated −270.5 ± 27.5 m (1971 to 2019), with an average rate of −5.6 ± 0.57 m a − 1 . Similarly, the four years (2015–2019) of field observations suggest that the retreat rate of PG and DDG is −6.7 ± 3 and −18 ± 15 m a − 1 , and the rate of modeled glacier mass loss is −0.29 ± 0.3 and −0.3 ± 0.3 m w.e. a − 1 , respectively. Furthermore, the ELA of the DDG and PG be-tween 1971 and 2019 increased by ~59 ± 38 and ~23 ± 19 m, respectively. The change in the longitudinal profile of the glaciers along the centerline between 2000 and 2017 shows the DDG and PG lost ~17 and 15 m surface ice thickness. The change in debris cover plays a critical role in the glacier surface lowering, shrinkage, retreat, and mass balance. Hence, we quantitatively evaluated the influence of the debris cover on summer ablation and terminus recession on two different characteristic glaciers (DDG and PG) with its potential effect on the mass balance process (area-volume loss).


Introduction
It is widely recognized that glaciers are highly sensitive to climatic change, and a lot of research has been focused on the response to increased global temperature (since 1950 onwards) and on evaluating future changes under the present climate change scenarios [1].Changes in glacier mass balance and its components are the most reliable indicators of climatic change [2].These are determined by winter accumulation (snowfall) and summer ablation (melting), which are, in turn, driven by anomalies in global and regional atmospheric circulation, as well as changes in local weather patterns and specific glacier characteristics (e.g., altitude, aspect, and supraglacial debris) [3][4][5].The Himalaya contains thousands of glaciers of widely varying properties, which are spread over nearly 37,000 km 2 with an east-west range of >2000 km [6].One of the significant characteristics of the Himalayan glaciers is that the glaciers are mainly debris-covered (70-80%) and have been receding since the end of the Little Ice Age [5,[7][8][9][10][11][12][13].The supraglacial debris on the surface of glaciers is commonly found to have significant control over the rate of ice ablation and recession [4,9,[14][15][16][17][18][19][20].It has been observed that the thickness of supraglacial debris significantly alters the glacier response to climate forcing [11,21,22].Using remote sensing and field-based data in different parts of the Himalaya, several studies, such as mapping debris cover and its correlation with glacier melting and recession, have been carried out [10,11,16,19,[21][22][23][24].As a consequence of the complex climate system, glacial geometry, glacier surface properties, and geology, the recession rates of the glaciers are variable [21].The studies suggest that most of the Himalayan glaciers are retreating between 5 and 20 m a −1 with the negative mass balance varying between 0.2 to 1.2 m w. e. a −1 and ELA fluctuations from 4900 and 5300 m asl [5,10].Glaciers are key markers of understanding climate change, so monitoring of glaciers in terms of estimating the overall changes of mass, volume, area, and length over space and time is essential.Despite the importance of the Himalayan glaciation, the knowledge of the glacial dynamics and forcing factors is scanty [10,25].Recent studies of Himalayan glaciers indicate wide variability in terminus retreat rate and mass balance in different sectors of the mountain range, primarily linked to the topography and climate of the region [10,[26][27][28].However, variable retreat rates of glacier termini and inadequate supporting field data (e.g., mass balance, ice thickness, velocity, etc.,) of the Himalayan glaciers make it challenging to develop a coherent picture of climate change impact.
In this study, we selected two glaciers with different characteristics, i.e., the Pensilungpa Glacier (PG) in Suru River and Durung-Drung Glacier (DDG) in Doda River basins, for a comparative study of glacier fluctuations between 1971 and 2019.The PG is characterized by a thick debris cover, while the DDG has a thin debris cover.The purpose of this study is to (i) measure the frontal retreat and area loss by the glaciers, (ii) reconstruct the present and past ELA followed by validation with field-based ELA, and (iii) estimate the centerline thickness lost by the glacier and mass balance.

Study Site
The DDG (30° 46′11.93″N-76°18′45.14″E)and PG (33° 49′5.37″N-76°17′06.44″E)are located in the SE and SW direction of Pensi-La pass (4454 m asl), respectively (Figure 1).DDG is the source of the Doda River, which originates from 4121 m asl and flows in the SE direction.The Doda River is the largest tributary of the Zanskar River (a tributary of the Indus River) and receives meltwater from 697 glaciers covering ~1080 km 2 glacierized area [29].The NE-facing, 23 km long DDG is the largest glacier in the valley occupying an area of ~72 km 2 .The glacier is almost clean, only 5.5% of the ablation area is covered by thin and scattered debris, and the average glacier slope is ~6°.The DDG is a compound valley-type glacier.The accumulation zone is broad with a gentle slope, formed by >6 cirques (Figure 1).The ablation zone of the glacier is narrower than the accumulation zone.It is bounded by well-preserved series of lateral and recessional moraines, indicating the glacier's dynamic behavior and past extent.Similarly, the Suru River (a tributary of the Indus River) originates from the PG at an altitude of 4670 m asl.The Suru Basin is supported by 252 glaciers, covering an area of 481 km 2 , 11% of the total area [20].The entire length of PG is ~8 km, covering an area of ~16 km 2 (including tributary glaciers), and ~17% of its overall glaciated area is covered with debris of varying thickness, which constitutes 30% of the ablation area.The PG is northflowing and originates from the valley's southwest slope with an average gradient of 8°, and terminates at an elevation of ~4670 m asl.Several longitudinal and transverse crevasses, small supraglacial ponds, and ice cliffs are present on the PG's surface (in the ablation zone).
The distribution of debris cover over the glacier follows a common pattern: the debris thickness decreases with increasing distance from the snout and elevation.Field observations (2015-2019) suggested that in the lower ablation zone (at the snout ~4700 m asl), the PG has a thick debris cover, and the thickness ranges between 15 and 40 cm.However, the upper ablation zone of PG (>4800 m asl) is thin debris-covered, and the debris thickness varies from 2 to 7 cm.In comparison, the ablation zone of DDG is partially debriscovered.The peculiar feature of DDG is that the glacier's center flowline is almost clean, and the left and right margins of the glacier are debris-covered.

Climate of the Study Area
The climate of the study area is controlled by the Indian Summer Monsoon and westerly circulations [5,30], and is temperate during summer and extremely dry and cold in winter.Due to the absence of instrumental data in the study area, we have used the Climate Research Unit Time Series (CRU-TS 4.03) data [5,31] to understand the climatic conditions.Here we present the temperature and precipitation of the interpolated gridded data (30.755N and 79.250 E) derived from CRU-TS 4 for 117 years (1901-2018) for peak months of accumulation period (January, February, November, and December) and melting period (June-September).The result shows that the February, November, and December temperature gradually increases, and the precipitation in January and December decreases.However, the temperature of June, July, and August shows decreasing trend, while the September temperature shows an increasing trend.
Similarly, the precipitation in June shows an increase, whereas no changes have been observed for the rest of the summer months.The increased temperature of early winter and late summer demonstrated the prolonging of the summer period and the shortening of winters.This increase in temperature (shift from solid to liquid precipitation) has been identified as one of the major factors responsible for the decrease in snowfall and enhanced melting of the glaciers over the Indian Himalayan region [5,10,20].The mean monthly minimum temperature (−10 °C) was observed in January, and the average minimum monthly precipitation (11 mm) occurred in November.The maximum monthly mean temperature (18 °C) and average precipitation (71 mm) were recorded in July.Overall, the long-term  climate data shows an increased winter temperature trend and a decrease in winter precipitation after 2000 (Figure 2).The data are taken from the peak months of accumulation period (January (A), February (B), November (C), and December (D)) and melting period (June (F), July (G), August (H), and September (I)).(E,J) showing the winter and summer variability, respectively.The result shows that the temperature of February (B), November (C), and December (D) is gradually increasing and the precipitation of January (A) and December (D) is decreasing.However, the temperature of June (F), July (G), and August (H) shows decreasing trend, while the September (I) temperature shows an increasing trend.Similarly, the precipitation of June shows an increase, whereas for the rest of the summer month no change has been observed.Overall, the long-term  climate data show an increased winter (E) temperature trend and a decrease in winter precipitation after 2000.

Datasets and Methods
In this study, the snout retreat, changes in the ELA, surface elevation, and modeled mass balance of the thick and thin debris-covered Pensilungpa and Durung-Drung glaciers, western Himalaya, Ladakh, India, are estimated using multiple sources (satellite images and field observations) (Table 1).The satellite images were acquired from the US Geological Survey (USGS, website) with minimal cloud cover.We used medium-resolution declassified Corona KH-4B (spatial resolution of 5 m), Hexagon-KH-9 (7 m), low-resolution Landsat series sensors TM/ETM + (30 m), Sentinel multispectral instrument (MSI) scenes from 2017 and 2019 (10 m) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)-L1A imageries.The satellite images selected for this study are based on their suitability (minimum cloud (<10%) and seasonal snow cover) for the most accurate glacier mapping.
Primarily, the acquired satellite images were co-registered by the projective transformation method using the Landsat ETM+ for 1999.However, the CORONA KH-4B and Hexagon KH-9 images were co-registered using a two-step approach (i) Projective transformation using 30 GCPs, and (ii) Spline adjustment of the subset image using the methodology used by many authors [20,[32][33][34].The identification of glaciers boundary (historic and recent) and mapping principles was followed by Shukla et al. [20,34] and Mehta et al. [5].Glacier boundaries were mapped manually for different points in time (1971, 1977, 1994, 2000, 2015, 2016, 2017, 2018, 2019, and 2021).Furthermore, high-resolution imagery from Google Earth™ was also used for the accurate demarcation of the boundary of the glacier.

Field Survey
From 2015 onwards, the snout position, frontal area loss, and glaciological mass balance measurements are ongoing on the glaciers in the Suru and Doda River basins [5,28,[34][35][36].More than ten permanent reference points were fixed near the frontal area of the DDG and PG snout on stable ground (Figure 3).The length from these reference points was then measured using a measuring scale up to the base of glacier ice from the left, right, and center of the snout.After that, the average length change was calculated based on the average recession from the glacier snout's left, right, and center parts for every studied year.The geographical position of each reference point was determined by Differential Global Positioning System (DGPS) and the real-time kinematic (RTK) points were taken with an accuracy of less than 0.1 m, with the base station set up at a distance of around 300 m from the glaciers snout.The obtained values were transferred on the map.We collected ten DGPS points, every year (2015 to 2019) at the margin of the snout of the DDG and PG and demarcated the frontal boundary of the glaciers.The glaciers front outlines of different years (e.g., 1971, 1977, 1994, 2000, and 2019) were plotted on the map to quantify the cumulative glacier frontal changes (Figure 4).Moreover, the field-based DGPS data of the glacier front were plotted on the map to calculate the frontal recession of the glacier between 2015 and 2019.In addition, we measured the recession from 2004 to 2019 using Google Earth imagery (Figure 4).

Estimation of Equilibrium Line Altitude (ELA)
The mass balance of the Himalayan glaciers is controlled by terrain, aspect, geographical position, and climatic condition of the area.The fluctuation of the ELA depends on the mass balance gradients of the glaciers.The modern ELA of the Himalayan glaciers ranges between 4400 to 5800 m asl [5,10,11,25,37,38].The valley morphology, moraine location, and snout position help to understand the glacier's past and present ELA.To determine the valley morphology, moraine locations, and snout position of the glaciers we used satellite images, referenced topographic map (Survey of India 1:50,000) along with the field surveys (GPS/DGPS), and hypsometry was calculated from SRTM (90 m) data.
Area Altitude Balance Ratio (AABR) Spreadsheet [40,42,47,48] We used an Excel spreadsheet for the calculation the ELAs for the study glaciers using the method of AABR given by Osmaston (2005).

Glacier Volume Estimation
The statistical regression relation of different approaches for estimating the ice volume of the glaciers was used in the present study.We used the following formula for ice volume calculation: V = 0.04 S 1.35   (1) V is glacier volume in km 3 , and S is glacier area in km 2 .This equation was derived by Liu et al. [51] based on thickness measurements by an ice-penetrating radar on 22 glaciers in the Tien Shan and Qilian Shan, western China.
V representing the glacier volume, A the glacier area, and c and γ are the two scaling parameters.
Here, we used three sets of scaling parameters, which have been applied by Cogley [52] and Frey et al. [53] for the Himalaya-Karakoram region.This equation is also used for all ICIMOD glacier inventories [54,55].This equation was derived by Chen and Ohmura [56], who used measurements from 63 glaciers to determine the related scaling parameters.Bahr et al. [57] derived the parameters in a theoretical study and LIGG et al. [58] established a thickness-area relation based on ice-thickness measurements on 15 glaciers [59].According to Cogley [52], the scaling parameters are given in Table 3. Table 3. Scaling parameters for thickness estimation are taken from Cogley [52].

Factor c
Exponent Haeberli and Hoelzle [60] presented an approach for estimating glacier volume based on the average surface slope and vertical glacier extent.The following formula is used to calculate the volume based on mean ice thickness along the central flowline: where τ is the average basal shear stress along the central flow line, ƒ is a shape factor, g is the gravitational acceleration (9.81 m/s 2 ),  is the density of ice (917 kg m 3 ), and α is the mean surface slope along the central flowline.ƒ is chosen constant as 0.8, considered the typical value for valley glaciers [2].The estimation of the mean ice-thickness of the entire glacier (hF) came out through the extrapolation from the mean ice-thickness along the central flowline (hƒ), following the assumptions of a semi-elliptic cross-sectional geometry and a non-branched glacier, a multiplication with (π/4) is added to Equation (3):

Glacier Mass Balance and Surface Elevation Change
Several researchers have used AAR and ELA methods to estimate the glacier mass balance in different parts of the world [11,[61][62][63][64][65][66].These studies have suggested that the specific mass balance can be calculated from the ELA and AAR [11,61,62,[66][67][68].Kulkarni et al. [67] developed a regression equation between AAR and specific mass balance to estimate the mass balance in the Indian Himalaya (Equation ( 5)).This regression equation between AAR and specific mass balance was derived from the in-situ measurements of two glaciers (Shaune Garang and Gor Garang) in the western Himalaya.This study suggested 0.5 AAR for zero mass balance.The relationship is Y = 2.4301 × X−1.20187 (5) where Y = mass balance of the glacier (m w.e.a −1 ) and X is the glacier's Area Accumulation Ratio (AAR).The AAR-based empirical relationship was derived by Pratap et al. [11] from the mass balance data of 11 glaciers of the Indian Himalayas from NE to NW.They suggested that the steady-state AAR of the glaciers is ~0.55 and the balanced budget ELA of the glaciers is ~5000 m asl.Hence, the relationship between specific mass balance (in m w.e.) and AAR takes the form (Figure 5A).Y = 1.9433 × X−1.3149 (6) where Y = Specific mass balance (m w. e. a −1 ) and X is the AAR.Similarly, we developed the empirical relationship between ELA and specific mass balance derived from the in situ measurements of the 11 glaciers in the Indian Himalayas to estimate the glacier mass balance.The empirical relationship is (Figure 5B) Y = −0.0011× Z + 5.0055 (7) where Y = Specific mass balance (m w.e.a −1 ) and Z is the ELA of the glacier.
This study uses these three empirical equations for specific mass balance estimation of the PG and DDG.
The longitudinal profile to estimate the surface elevation change is constructed by using SRTM version -3 (v3) of 2000 and ASTER DEM of 2017 (Figure 6).The surface elevation change was measured using the longitudinal profile of two periods (2000 and 2017).The longitudinal profiles along the centerline of the glacier were drawn and compared.After comparing these two longitudinal profiles, surface lowering was measured in an Excel spreadsheet that helped to compare these two longitudinal profiles; surface lowering was measured in an Excel spreadsheet that helped estimate the thickness loss by the glaciers between 2000 and 2017.For a generation of ASTER DEMs 2017, we adopted the method of Garg et al. [28] using the stereo pair bands: 3N (nadir-looking) and 3B (backward-looking).This process involved inputting ground control points, generation of tie points, creation of epipolar images, image matching, and DEM geocoding.With SRTM DEM as vertical and Sentinel MSI as horizontal references, nearly ten ground control points were collected mainly from stream confluence, glacier peaks, mountain ridges, permanent boulders, or any other prominent feature for enhancing the accuracy of the DEMs.Nearly 97 tie points were used to generate the epipolar images (right and left), obtaining a y-parallax of 1.8 in 2017.Finally, the output DEMs were extracted at 30 m spatial resolution, possessing high terrain details (level 4).The DEMs were then edited using the spline method and median filter with kernel size 5 × 5 to rectify the unnatural peaks and void and smoothen the DEMs.The generated DEMs covered the entire glacier; however, they comprised certain voids (data gaps) primarily due to cloud or snow cover and shadowed regions.The DEM co-registration (horizontal and vertical adjustments) and surface elevation change were estimated by Garg et al. [28].

Terminus Retreat and Area Changes
The terminus retreat and area change of PG between 1971 and 2019 is divided into three categories: (i) remote sensing-based total terminus retreat from 1971 to 2019 was ~270 ± 27 m at an average rate of 5.6 ± 0.57 m a −1 , (ii) Google Earth-based total terminus retreat from 2004 to 2017 was ~110.5 ± 18.25 m at the average rate of 7.4 ± 1.23 m a −1 and (iii) field observations from 2015 to 2019 with glacier retreat of ~27 ± 11.5 m at an average rate of 6.7 ± 3 m a −1 .Overall, the glacier lost about 3% of its total length (Figure 4; Table 4).The area loss by the PG in its pro-glacial region for the period between 1971 and 2019 was ~(−) 1.5 ± 0.02 km 2 (8% of total area), an average rate of (−) 0.03 ± 0.004 km 2 a −1 , out of which ~(−) 0.09 ± 0.07 km 2 was frontal area loss.The measured field data suggested that the glacier lost ~1043 ± 179 m 2 (0.0010 km 2 ) in the proglacial area between 2015 and 2019.
The frontal geometry of the DDG is very peculiar, and it appears like a colossal lobe.The glacier's snout is very broad, approximately 2.0 km wide (Figure 4).Due to the large extension of the terminus, the retreat of the glacier snout is reflected as oscillation or quasistagnant.The broad snout of the DDG creates wide discrepancies for the frontal retreat.The left part (NW) of the glacier retreat was ~576 m, and the right part (NE) of the glacier retreat was ~104 m, while, at the same time, the center part of the glacier retreat was ~1195 m between 1971 and 2019.The estimation from Google Earth shows that the DDG retreated by ~227 ± 131 m with an average rate of 17.46 ± 10 m a −1 between 2004 and 2017.Further, the field observation suggests that the DDG retreated ~61.5 ± 33 m with an average rate of 15.4 ± 8 m a −1 from 2015 to 2019.Overall, the DDG reduced ~(−) 624 ± 547 m length with an average rate of ~(−) 12 ± 11 ma −1 which is ~3% of the total length in the last 48 years (Figure 4).
During the study period, the DDG lost ~7.78 ± 0.05 km 2 area (10% of total area) with a rate of ~0.16 ± 0.001 km 2 a −1 , out of which 0.91 ± 0.07 km 2 was the frontal area (which is 12% of the total area lost).

Equilibrium Line Altitude (ELA) and Accumulation Area Ratio (AAR) changes
The results of ELA and AAR estimated for DDG and PG, based on AAR, THAR, AABR, and MELM methods, are presented in Table 5.It shows that the ELA of DDG and PG in 1971 was at ~5192 ± 29 and 5198 ± 23 m asl, respectively, while the AAR of the glacier was ~0.54 for DDG and ~0.51 for PG.The terminus of the DDG and PG was located at the altitude of 4090 and 4610 m asl in 1971, respectively.Between 1971 and 2019 (48 years), the ELA of DDG and PG ascended by ~59 ± 38 and ~23 ± 19 m, respectively.The recent ELA of the glaciers is located at the ~5251 ± 52 (DDG) and 5221 ± 31 m asl (PG), respectively (Table 5).Moreover, the estimated AAR is ~0.46 for the DDG and ~0.43 for the PG in 2019.An in situ study of PG suggested that the ELA of the glacier fluctuated between 5215 and 5223 m asl between 2016 and 2019, and the mean ELA of the glacier was at 5221 m asl [5].The 2019 estimated ELA of PG is well correlated with the in situ (calculated) ELA.

Glacier Mass Balance and Thickness Change
We used the AAR/ELA method to estimate the glacier mass balance rate (specific balance) of DDG and PG for the years 1971, 2002, and 2019.The result shows the continuous mass loss since 1971 (Table 7).The modeled average specific balance of the DDG was −0.

Longitudinal Profiles and Elevation Change
For surface elevation change, a longitudinal profile along the centerline was drawn between 2000 and 2017 (Figure 6).The longitudinal profile shift between 2000 and 2017 shows that PG lost ~15 m surface ice thickness, while DDG lost ~17 m ice thickness (Table 8).The PG lost maximum thickness between 5100 and 4900 m asl, where the debris thickness is <5 cm, while the minimum thickness was lost between 4800 and 4700 m where the debris thickness was >10 cm (Table 8).However, between 4750 and 4700 m asl, no loss in ice thickness was observed.Compared to PG (thick debris-covered) and DDG (thin debriscovered), the ice thickness loss on the centerline shows linear trends except for some places (Figure 6).The upper ablation zone (4950-5150 m asl) of DDG lost higher ice thickness (average −23 m), while the lower ablation zone (4700-4900 m asl) lost less ice thickness (average −10 m), whereas DDG showed linearity with ice thickness loss, except in some places between 4950 and 5000 m asl (Figure 6).Generally, melting at lower elevations is higher than in the upper areas; at PG, that is not the case, probably due to the thick debris cover.The elevation-wise thickness lost by the glaciers is presented in Table 8 and Figure 6.

Discussion
The glaciers' response to ongoing climate change across the Himalayas is complicated to understand due to the complex topography and lack of systematic or long-term field measurements.Within each of the individual mountain ranges, the large amount of scattering in the pattern of relative changes suggests that the local factors are essential in controlling the behavior of the glaciers.These may include differences in local climatic history and/or differences in the intrinsic sensitivity of the glaciers to climatic change as a result of differences in slope [69][70][71], aspect, elevation, hypsometry, and presence of debris cover [72] or rate of mass turnover [73,74].Similar patterns have been observed in the Himalayan region [18,[75][76][77][78][79][80]; however, these studies have reported that the smaller glaciers' loss is much faster than the larger ones.The DDG and PG are located in the same climatic regime and flow in the same direction (NNE).The geomorphological evidence suggests that these two glaciers were connected in the past, but now they are dispersed in different directions [30].Geometrically and morphologically, the DDG and PG are different, DDG has a thin debris cover (only 5.5% area) and is the biggest glacier in the basins (Doda and Suru), having a broader snout (~2 km).DDG (72 km 2 ) is five times bigger than the PG (16 km 2 ), and both glaciers reduced in length ~624 ± 547 and 270 ± 27 m with the rate of 12 ± 11 and 5.6 ± 0.5 m a −1 between 1971 and 2019, respectively.The PG has lost its length at a two times slower rate than the DDG, while both glaciers lost a similar percentage of the length, which is ~3% of the total length.This comparative study shows that the thick debris-covered ablation front of PG retreats more slowly, whereas the thin debriscovered DDG front retreats at a faster rate.Scherler et al. [21] suggested that the glacier retreat rate varies with >20% debris cover zero, whereas it is high for the debris-free glacier.This can be correlated with the influence of debris (>17%) over PG with a slower retreat rate compared to DDG with less (5.5%) debris cover.A debris cover influences the terminus dynamics and modifies a glacier's response to climate change [21].Surface ablation rates are generally increased in the presence of a thin (<5 cm) debris cover but are significantly reduced when a thick (>5 cm) debris cover is present [5,10,15,21].Thin and patchy debris cover reduces the albedo and elevates shortwave radiation absorption, whereas ablation rates are strongly reduced further down the glacier due to the insulating effect of thicker debris [10,20].Additionally, we observed that the broad front of the DDG indicates the oscillatory nature.The left and central parts of the glacier lost higher length, but the right part of the snout is quasi-stagnant (Figure 3).Due to the oscillatory nature, the center part of the snout of the DDG advanced 33 ± 12 m between 2013 and 2017, whereas the glacier shows high uncertainty and a total retreat of ~36 ± 60 m in the same periods (2013)(2014)(2015)(2016)(2017).
An earlier study by Shukla and Qadir [81] suggested that the retreat of DDG was ~691 and PG was ~298 m between 1977 and 2013 with an average rate of ~19 and ~8 m a −1 , respectively.Similarly, Kamp et al. [24] suggested that the retreat of DDG was ~311 m with an average rate of ~9 m a −1 between 1975 and 2008.However, the retreat of PG was ~1887 m, with the average rate of the retreat being ~61 m a −1 between 1975 and 2006 [24].Similarly, a recent study by Rashid and Majeed [82] shows that the DDG retreated 925 m at the rate of 20 m a −1 between 1971 and 2017; it seems that the authors took the recession from the center of the snout.In comparison, our study suggested that the snout of DDG retreated ~576 m from the left side (NW), ~104 m from the right part (NE), and ~1195 m from the center part of the glacier between 1971 and 2019.The present study suggests that the center part of the snout of the DDG advanced between 2013 and 2017, whereas the PG was in continuous retreat.
Lee et al. [13] studied >14,000 glaciers in the Himalayas and suggested that the Himalayan glaciers lost 40% of their area between LIA (400 to 700 years ago) and recent years.Mehta et al. [5] studied the PG between the same period of LIA (400 and 700 years ago and recent years) and suggested that glaciers lost 16% of their total area and 27% of their total length based on the field data.The comparative studies showed that PG lost a lower percentage of the area.Lee et al. [13] also indicated that the ranges of mass loss rate of the Himalayan glaciers between LIA and the recent year were −0.011 and −0.020 m w.e.a −1 , which was a lower rate than current trends.Maurer et al. [12] estimated the mass balance of the Himalayan glaciers and suggested that the average mass balance rate of the Himalayan glaciers was −0.33 ± 0.13 m w.e.a −1 between 1975 and 2016.The present study based on ELA and AAR model suggested that the modeled mass balance rate of DDG is −0.36 ± 0.3 m w.e. a −1 and PG is −0.39 ± 0.26 m w.e. a −1 between 1971 and 2019.DDG and PG lost 12% and 17% of ice volume during the same periods, respectively.The in situ mass balance study of PG shows a rate of −0.36 m w.e. a −1 between 2016 and 2019, which supports the result of the estimated mass balance rate (−0.39 ± 0.26 m w.e. a −1 ) of the PG glacier by Mehta et al. [5].The ELA and AAR-based mass balance of the DDG and PG is well correlated with the modeled and in situ mass balance of the Himalayan glaciers.The in situ mass balance rate of the Dokriani (Uttarakhand), Chhota Shigri (Himachal Pradesh), Pokalade (Nepal), and Pensilungpa (Zanskar, Ladakh) is −0.69, −0.32, −0.56, and −0.36 m w.e. a −1 , respectively [5,10,83,84], whereas the remote sensing and geodetic-based studies suggested that the mass balance rate across the Himalaya was −0.33 ± 0.13 m w.e.a −1 between 1975 and 2016 [12], the Himalayan sub-region was −0.35 m w.e. a −1 between 1999 and 2008 [85], Hindu Kush-Karakoram-Himalaya region (HKKH) was −0.21 ± 0.05 m w.e. a −1 between 2003 and 2008 [86], High Mountain Asia lost ice with the rate of −0.34 ± 0.06 m w.e. a −1 between 2000 and 2016 [26].However, Garg et al. [87] suggested that the average mass balance of the glaciers in the Suru River basin was −0.69 m w.e. a −1 , which is higher than the present study.The ELA and AAR-based mass balance of DDG and PG are well correlated with above mentioned in situ (glaciological) and remote sensing (geodetic) based mass balance studies.
The combination of continuous retreat in length, reduction in volume, and upward shift of the ELA can be attributed to the rising temperature and changes in precipitation trends in the Himalaya.More definitive connectivity between glacier recession, ELA change, and recent climate change requires a regional-scale study.Nevertheless, as far as the DDG and PG are concerned, if the present climate trend continues, it is expected that the glaciers will continue to decrease in ice volume.
The morphological parameters of these two glaciers (DDG and PG) are different, while the topographical characteristics are almost similar.The rate of area and length loss by the DDG is about eight times higher than the PG.According to CRU-TS 4 meteorological data (1901-2018), both glaciers are located in the same grid and showing that the winter temperature and summer precipitation have increased since 1980 (Figure 2).The difference in glacier recession cannot be well explained by climate records individually.Thus, we calculated the topographical characteristic (slope, aspect, altitude, and size) of both glaciers to track their influences on glacier changes.The average slope of the PG is 8°, and DDG is 6° and both glaciers face the NE direction.The ELA of both glaciers is located at 5221 m asl for PG and 5251 m asl for DDG, while the mean median altitude of the PG is 5320 m asl and DDG is 5250 m asl.Aizen et al. [100] suggested that the melting of the glaciers is driven by solar radiation (up to 86%), based on the meteorological station records.However, Shukla et al. [20] and Kumar et al. [3,18] suggested that the aspects and slope are important factors affecting glacier changes.There is no significant difference in topographic and climate factors between PG and DDG.Thus, we extrapolate that the incoming solar radiation affected by the aspects and slopes cannot cause the difference in glacier changes in Suru and Doda River basins.We find that the DDG (72 km 2 ) is bigger than PG (16 km 2 ); consequently, we may speculate that different retreat rates of the glaciers are probably caused by the difference in glacier size.The debris thickness could also affect the ablation of glaciers [101][102][103].The presence of supraglacial debris covers influences glacier processes.Depending on thickness, debris cover may either enhance or retard the ablation process [10,20,21].PG covered 17% (1.8 km 2 ) of their area with debris, while the DDG covered 5.5% (3.8 km 2 ) of the area with debris.We observed during the period between 2015 and 2019 that the lower ablation zone of the PG is covered by thick debris and reduces the frontal retreat of the glacier, whereas the lower ablation zone of the DDG is covered by thin or scattered debris and enhances the frontal retreat of the glaciers.The supraglacial ponds/lakes and proglacial lakes act as ablation hotspots and increase the ablation rate of the glacier [81,104,105].Two large proglacial lakes format the front of the DDG, covering an area of about 2 km 2 .Due to the proglacial lake, the snout of the DDG retreats faster than PG, and the ice of the margin of the snout breaks continuously.Although we cannot find a strong correlation between glacier retreat, meteorology, and topography factors, both glaciers are located in the same climatic zone, aspect, and slope.Nevertheless, we expect that the glacier size, debris cover, and glacier lake probably influence the rate of glacier retreat in the Suru and Doda River basins.

Conclusions
In this study, we present the result of length, area, ELA, volume, and mass balance changes of two glaciers (DDG and PG) located in the Doda and Suru River basins.The DDG and PG are located in the same climatic zone, aspect, and slope.Nevertheless, there is a lot of heterogeneity in the recession patterns of both glaciers.The DDG lost an area of 7.78 km 2 (10% of the total area) which was five times higher than PG (1.5 km 2 ; 8% of the total glacier area).The retreat rate of the DDG is 13 m a −1 between 1971 and 2019, which is two times higher than the PG (5.6 m a −1 ).
Moreover, the ELA of DDG shifted 59 m upward between 1971 and 2019, while the ELA of PG shifted 23 m in the same period.However, the ELA of both glaciers was located in the same altitude zone (5200-5300 m asl) in 2019.We expect that the glacier size, debris cover, and glacier lake probably influenced the rate of glacier retreat in the Suru and Doda River basins.This trend in glacier retreat rates over time is similar to that for other glaciers in the western, central, and eastern Himalayas.The increase in temperature and reduction in precipitation mainly causes the retreat of the glaciers.The increasing amount of precipitation could not compensate for the glacier mass loss due to the temperature increase in Zanskar Himalaya (Suru and Doda basins).This suggested that the warming of the climate is also responsible for the glacier retreat in the study area.The rate of glacier retreat is not only controlled by climate change but also by the topographic setting and morphology of the glacier.The comparative study also confirms the possible influence of factors such as snout geometry, glacier size, elevation range, slope, aspect, debris cover, and presence of supra and proglacial lakes other than the climate in the observed heterogeneous response of the studied glaciers.Therefore, these factors need to be accounted for in more detail in future studies for a complete understanding of the observed glacier changes and responses.

Figure 1 .
Figure 1.Location and contour map of the study area, (A) showing glacial distribution in Suru and Doda basins and black dash line indicates the boundary of Pensilungpa Glacier (PG) and Durung-Drung Glacier (DDG) catchments.The inset India map shows the location of the study area.(B) Digital Elevation Model (DEM) of DDG and PG showing the glaciers and snow cover in the study area.

Figure 2 .
Figure 2. Annual and seasonal variability in the climatic conditions of the study area derived from CRU-TS 4.03 data (30.755N and 79.250 E) for the period 1901-2018 (Modified after Mehta et al. [5]).

Figure 3 .
Figure 3. Field photographs showing the snout position of the glacier in 2015 and 2019 (Modified after Mehta et al. [5]), ((A,B), respectively) PG and ((E,F), respectively) DDG.Red circles on photograph (B,E) is indicating the area vacated by the glaciers and red arrows (in photographs (A,B) and (E,F)) are showing the reference points.(C,D,G,H) the close-up view of the reference points used for measuring the snout retreat with the help of chain tape survey.

Figure 4 .
Figure 4. Temporal changes in the snout position of the PG and DDG during the period between 1971 and 2019 (Modified after Mehta et al. [5]).(A) The respective glacier boundaries for 28, September, 1971 (Corona KH-4B), (B) 6, October 1977 (Hexagon KH-9), (C) 27 August 1994 (Landsat TM), (D) 4 September 2002 (Landsat ETM + ), (E) 15 September 2019 (Sentinel Multispectral Imagery).(F) The glacier boundaries are dropped on the Corona KH-4B showing the glaciers recession between 1971 and 2019.(G) Snout retreat of PG measured using Google Earth Image during the periods between 2004 and 2017, (H) field photograph showing the snout of the PG and red arrow indicating reference point (big boulder).(I) Snout and proglacial kidney-shaped lake (blue arrow) of DDG, and (J) snout retreat of DDG during the periods between 2004 and 2017.

Figure 5 .
Figure 5. (A) The relationship between specific balance and Accumulation Area Ratio (AAR).(B) Relationship between specific balance and Equilibrium Line Attitude (ELA) of Himalayan glaciers.The equilibrium state ELA and AAR is indicated by the intersection of the best-fit line and the vertical axis is specific balance.The ELA and AAR for equilibrium state of the glacier are obtained to be 5000 m asl and 0.55, respectively.Note that 12 years of positive balance occurred in 1975/1976 and 1982/1983 (Gara Glacier), 1982/1983 (Gor Gorang Glacier), 1983 and 1989 (Shaune Gorang Glacier), and 2005, 2009, 2010, and 2011 (Chota Shigri).The average ELA and AAR of the Himalayan glaciers is about 5030 and 0.45, respectively between 1975 and 2011.

Figure 6 .
Figure 6.Surface elevation change of the PG and DDG along the centerline (AA' and BB' for PG and CC' and DD' for DDG) at different times.(A) SRTM DEM (2000) and (B) ASTER DEM (2017).The methodology of the surface elevation change is discussed in the text.The longitudinal profile along the centerline of the (C) PG and (D) DDG, showing the elevation change between 2000 and 2017.

27 A
recent study by Mehta et al.[5] based on glaciological mass balance suggested that the average specific balance of PG was −0.36 m w.e. a −1 during the periods 2016/2017-2018/2019.The mass wastage rate in the ablation zone is about 0.82 m w.e. a −1 , with the mass gain in the accumulation area, being about 0.27 m w.e. a −1 between 2016/2017 and 2018/2019.The modeled value of the specific mass balance of PG is near to in-situ measurement.

Table 1 .
Data used for monitoring of glaciers in this study.

Table 2 .
Morphometric indices and their equation for estimate the ELAs of PG and DDG between 1966 and 2017.

Table 4 .
Frontal retreat of DDG and PG glaciers between 1971 and 2019.

Table 5 .
Estimated ELAs, ELAs changes, and AAR of the DDG and PG glaciers between 1971 and 2019.

Table 6 .
Estimated volume of DDG and PG glaciers between 1971 and 2019.

Table 7 .
Modeled mass balance of DDG and PG glaciers between 1971 and 2019.

Table 8 .
The elevation-wise thickness lost by the glaciers.