Dulung Proglacial Lake, Suru Sub-Basin, Western Himalaya: Evolution, Controls and Impacts on Glacier Stability

Proglacial lakes are continually developing and expanding across the Himalayan glaciered terrain in response to climate change. These lakes are known to destabilize the glaciers by enhancing their frontal ablation, causing higher than average glacier area and mass losses. Thus, to comprehend the dynamics of proglacial lakes and their influence on the overall glacier health, we study the lake-terminating Dulung Glacier located in the Suru sub-basin, Ladakh, western Himalaya and compare it with the adjacent land-terminating Chilung Glacier. The pronounced melting of the Dulung Glacier, supported by glacier topography (surface gradient between accumulation and ablation zone) and valley morphology (wider near the snout and narrower downwards), seems to be the prime reason for the formation, accommodation and sustenance of the proglacial lake. The expansion in proglacial lake (.008 km2a−1) during 1977–2018 is accompanied by an enhanced degeneration of the Dulung Glacier (mass balance: −.47 ± .06 m w.e.a−1, shrinkage rate: .3 ± .001% a−1; retreat rate: 32 ± .7 ma−1, surface ice velocity reduction: 16%), which has accelerated post-1993. In comparison, land-terminating Chilung Glacier shows lower degeneration rates (mass balance: −.28 ± .02 m w.e.a−1; shrinkage rate: .2 ± .001% a−1; retreat rate: 17 ± 0.7 ma−1, surface ice velocity reduction: 8%) during 1971–2018. This suggests a substantial impact of the proglacial lake in enhancing the Dulung Glacier’s sensitivity towards climate change compared to the Chilung Glacier. If the current rate of lake expansion continues, it would further enhance the Dulung Glaciers’ degeneration rates, thus impacting its stability.


INTRODUCTION
The Himalayan glaciers' response has markedly changed due to climate fluctuations in the last few decades (Kumar et al., 2005). This alteration is evident from the enhanced mass wastage, shrinkage and retreat of most Himalayan glaciers, however, with few exceptions Kulkarni and Karyakarte, 2014;Azam et al., 2018;Shukla et al., 2020b). Such glacier degeneration poses a serious problem for regional-scale water budget and resource management (Barnett et al., 2005;Immerzeel et al., 2020). Another significant threat to the Himalayan cryosphere as a consequence of glacier retreat and shrinkage is the formation of glacial lakes (Nie et al., 2017;Shugar et al., 2020). Glacial lakes have become quite frequent in the Himalayan glaciated environment, with majority of them formed in the last 50 years (Raj, 2010). According to a report by the International Centre for Integrated Mountain Development, there are 15,000 glaciers and 9,000 glacial lakes in the Himalayan region (Mool, 2005). Amongst all these lakes, the proglacial lakes (PGLs), specifically the moraine-dammed lakes, have been recognized as the most disastrous ones. The PGLs typically develop in over-deepened glacier beds behind the lateral and terminal moraines that were formed during the Little Ice Age (Benn et al., 2012;Westoby et al., 2014). Largely, the PGLs receive input from the meltwater generated supra-, en-or subglacially. Their origin can be directly linked to the climateinduced thinning/negative mass balance of the associated laketerminating glaciers (glaciers that reach or terminate in a lake) (Reynolds, 2000;Singh et al., 2011;Kääb et al., 2012;Gardelle et al., 2013;Neckel et al., 2014;King et al., 2019). Once formed, these lakes grow over time due to the differential thermal regimes of ice and warm water that instigate calving, leading to increased terminus retreat of the mother (lake-terminating) glacier (Kirkbride and Warren, 1997;Sakai et al., 2000;Nie et al., 2017;Shukla et al., 2018). In contrast, the dominant ablation process on land-terminating glaciers (glaciers that reach or terminate on dry land) is through surface melting (Singh et al., 2011). The lakes also have the ability to drain out a huge volume of water stored within them to the downstream areas resulting in glacial lake outburst floods (Richardson and Reynolds, 2000;Quincey et al., 2007;Liu et al., 2008). Therefore, to understand the contribution of the PGLs in altering the glaciers' response and future prediction of any hazard, it is imperative to regularly monitor and study the dynamics of PGLs.
Glacial lakes, being typically situated at remote locations, are difficult to monitor through field investigations owing to obvious challenges of harsh weather, inadequate workforce and logistics. Under such conditions, remote sensing provides the suitable platform to assess lake dynamics and hazard potential on a regular basis (Huggel et al., 2002;Schneider, 2004;Quincey et al., 2005;Bolch et al., 2008;Zheng et al., 2021). Remote sensing techniques have been widely employed in the Himalayan region to monitor the glacial lakes (Quincey and Richardson, 2007;Allen et al., 2009;Gardelle et al., 2011;Raj et al., 2012;Worni et al., 2013;Zhang et al., 2015;Nie et al., 2017;Shukla et al., 2018). These studies have either inventoried the glacial lakes (Worni et al., 2013;Zhang et al., 2015;Aggarwal et al., 2017;Shukla et al., 2018;Pandey and AliChamati, 2021), or investigated their impact on two-dimensional (area, retreat) (Basnett et al., 2013;Mir et al., 2018;Kumar et al., 2021;Pandey and AliChamati, 2021), or three-dimensional (ice thickness, velocity) glacier parameters (King et al., 2018;Liu et al., 2020;Pronk et al., 2021). Though there are many studies on mapping and monitoring Himalayan glacial lakes (Richardson and Reynolds, 2000;Salerno et al., 2012;Majeed et al., 2020), questions on their evolution and role in modifying the overall glacier response still exist. With the second-highest count of glacial lakes after the Brahmaputra, the Indus basin shows an increase in their number during 1990-2010 (Zhang et al., 2015). Meanwhile, though less frequent, glacial lake outburst floods are well reported from the latter, with a recent event (2014) at the Gya village of Ladakh region (Schmidt et al., 2020). Therefore, studies focusing on lake dynamics and their impact on glacier stability need earnest attention to not only predict any hazardous situation in the future but also comprehending the associated glaciers' state.
The present study focuses on the Dulung Glacier and associated PGL at its snout, located in the Suru sub-basin, Ladakh, western Himalaya. Recently, Mir et al. (2018) also evaluated the Dulung Glacier lake and its impact on glacier response (i.e., shrinkage and retreat rates). Though the study could partly explain the glaciers' behavior, the comprehensive glacier response remains unexplored. Also, the causative factors influencing the lake formation were only partly assessed. Therefore, to have an in-depth understanding of the glacial lake evolution, apart from Dulung, we also investigated a neighboring glacier named Chilung. Even though Dulung and the adjacent Chilung Glacier share similar geographical and climatic settings, the latter does not have any glacial lake. Therefore, it was interesting to assess and compare the respective valley and other topographic conditions possibly contributing to the lake evolution on one glacier and not on the other. From this study, we aim to: 1) assess the conditions linked with the PGL formation and 2) understand the impact of the PGL expansion on the response of the lake-terminating (Dulung) in contrast to the land-terminating (Chilung) glacier in the adjoining valley.

STUDY AREA
The present study involved two mountain glaciers, namely Dulung (size: 14 km 2 ) and Chilung (size: 10 km 2 ) (Figure 1), located in the Suru sub-basin, Ladakh, western Himalaya. The study area extends between the latitude 33°52′ to 33°57′ N and longitude 76°09′ to 76°13′ E. The tributary glaciers contributing to the main trunk are considered single glacier entities (GLIMS, 2015). The Dulung Glacier comprises four glacierets (D1, D2, D3, and D4) along with the main trunk, while the Chilung Glacier comprises two main trunks (C1 and C2) and a glacieret (C3) (Figure 1). In the past, these smaller glacierets would have existed as a single glacier entity. The Dulung and Chilung Glacier's snouts are located at an elevation of~4,450 m above sea level (m asl) and~4,550 m asl, respectively, with the presence of a PGL at the formers' snout ( Figure 1). The mean elevation of the glacial lake is 4,407 m asl. The glaciers occupy the maximum altitude of 5,700 m asl (Dulung) and~5,600 m asl (Chilung). The total length of the Dulung Glacier is approximately 6 km, while that of the Chilung Glacier is approximately 3 km in 2018.
Both the glaciers drain their meltwater into the Suru River, a tributary of River Indus. Suru River originates from the Pensilungpa glacier, which lies near the Pensi La pass. The Suru River is the main hydrological source in the Kargil district of the Ladakh union territory. The villages situated in close proximity (downstream) to the Dulung Glacier are Rangdum, Parkachik, Tangole and Panikhar, with the respective aerial distance of about 21, 58, 67, and 74 km. The regional climatic conditions are harsh, with the occurrence of heavy snowfall during October to June.

Dataset Used
A total of 16 Landsat [Multispectral Scanner (MSS), Thematic mapper (TM), Enhanced thematic mapper plus (ETM + ) and Operational land imager (OLI)] scenes were used (Supplementary Table S1) in this study to map the PGL boundary and extract glacier parameters such as lake and glacier area, retreat, debris cover, snow line altitude (SLA) and surface ice velocity (SIV) for the period 1971/77-2018. Historical glacier parameters (glacier area, length) for the year 1971 were mapped from declassified Corona KH-4B imagery. All the acquired satellite images were preferentially of the peak ablation period, i.e., August-September, having minimal cloud and seasonal snow cover to ensure accurate mapping and assessment of glacier and lake parameters (Paul et al., 2013). The Landsat, declassified Corona KH-4B and Sentinel Multispectral Instrument (MSI) satellite images were downloaded from the Earth Explorer website (http:// earthexplorer.usgs.gov). Besides, high-resolution Linear Imaging Self Scanning Sensor (LISS)-4 (5.8 m) imagery, obtained from nrsc. gov.in, along with the Sentinel MSI FIGURE 1 | Location map of the study area. (A) The study investigates the proglacial lake (PGL) located in the Suru sub-basin, western Himalaya and its impact on the associated Dulung Glacier's health and compare the response with neighboring Chilung Glacier. (B) The glacier boundaries marked with yellow (Dulung) and red (Chilung) colors are overlaid on the Landsat OLI image acquired on 29 August 2018 with the band combination (false-color composite: R: short-wave infrared, G: nearinfrared, B: red). The Dulung Glacier comprises a main trunk along with four glacierets (D1, D2, D3, and D4), while the Chilung Glacier comprises two main trunks (C1 and C2) and a glacieret (C3). To understand the conditions favoring the lake formation, the transverse profiles of the valleys were derived along the equally spaced colored (

Methodology
The study involved the assessment of multidimensional glacier parameters (area, retreat, debris cover, SLA, SIV, surface elevation changes and mass balance), during the period 1971/77-2018. The time frames for assessing glacier and lake parameters were kept the same, except for SIV and mass balance, which are evaluated for 1993-2018 and 2000-2017, respectively.

Glacier and Proglacial Lake Mapping
Before mapping and extraction of glacier parameters, the raw image data was pre-processed. For this, the satellite images were coregistered with respect to the reference/base image, which in this study is Landsat ETM + image, acquired on 4 September 2000. For 1971For , 1994For , 2000, the glacier boundaries were obtained from Shukla et al. (2019). The glacier boundaries for these years were delineated using a hybrid approach, i.e., normalized difference snow index (NDSI) for snow and ice boundary and debris cover delineated manually. However, we have not evaluated the lake dimensions for 1971 due to its shadowing in the concerned satellite imagery. Glacier boundaries (mainly the debris cover boundary) for the remaining years (1993, 2008, and 2018) and lake boundaries for all the studied years were manually digitized using false-color composite with RGB band combinations of SWIR-NIR-red, NIR-green-blue and red-green-blue. Band combination SWIR-NIR-red was used for the identification of the glaciers (appearing blue-green). Visual interpretation like terminus shadow, convex/concave shape, emergence of stream and proglacial lake was used to identify the glacier terminus. High-resolution imageries from Google Earth were also referred to identify the actual position of glacier snout and for 3-dimensional visualization. The upper boundaries of the glaciers were kept the same, assuming that they remained unchanged over the studied period. A band combination of NIR-red-green was used to identify water surfaces (Paul et al., 2013), as water shows maximum absorption in the near-infrared and middle infrared bands with slight reflectance in the visible bands. This sharp spectral contrast helps distinguish between the glacial lake and the surrounding topography/feature (Pietroniro and Leconte, 2000).
3.2.2 Extraction of Glacier and Lake Parameters 3.2.2.1 Estimation of Lake and Glacier Area, Retreat, Debris Cover and Snow Line Altitude The glacier area and supraglacial debris cover estimates were directly extracted from the manually delineated glaciers' boundaries. The glacier length was measured along the central flow line of the glacier. The frontal retreat was measured using the parallel line method following Garg et al. (2017) and Shukla et al. (2020b). Theoretically, snow line demarcates the boundary separating the accumulation and the ablation zones of the glacier. SLA fluctuates both seasonally and annually and when estimated at the end of hydrological year, it is considered equivalent to equilibrium line altitude (ELA) (Guo et al., 2014). The snow and ice facies show maximum contrast in their spectral signature in the SWIR and NIR bands. The same principle was used for delineating the annual fluctuations in snowlines during 1977-2018. The altitude of snowlines was extracted using the SRTM DEM as per Shukla et al. (2020b).

Assessment of Glacier Dynamics
Surface Ice Velocity. The SIV was estimated as per Leprince et al. (2007) and Shukla and Garg (2020), using the feature tracking method on optical imageries (1993/94, 1999/2000, 2009/10, 2017/ 18). The basic principle involved is the correlation between two orthorectified and co-registered imageries acquired at two different periods. Correlation of the image pairs was performed in Co-registration of Optically Sensed Images and Correlation (COSI-Corr) using the initial and final window size of 64 down to 32 pixels for ETM+ and OLI (PAN band) and 64 down to 16 pixels in case of TM (NIR band) images. The images were resampled to 60 m by selecting the step size of 4 and 2 pixels for ETM+/OLI and TM images, respectively. The correlation produced three components: north-south, east-west displacements and Signal to Noise Ratio images. The signal-tonoise ratio measures the established correlation's accuracy and may be introduced due to the poor correlation of the images. Therefore, these poorly correlated pixels need to be eliminated so that the error does not propagate to the displacement measurements. For this, the data filtering algorithm was used, removing pixels with a signal-to-noise ratio less than .9 (Saraswat et al., 2013;Gantayat et al., 2014). After that, de-striping was performed separately on the east-west and north-south displacement images, followed by a median filter to smoothen the images.
Further, the directional filter was applied to remove miscorrelations from the displacement images. The offsets in north-south and east-west directions were then utilized to obtain the net surface displacement using the Euclidean norm (Leprince et al., 2007). Finally, the annual SIV was estimated by normalizing the displacement for 365 days interval. Thereafter, the SIV was extracted along the central flow line considering the fact that the glacier velocity is maximum at the center due to the minimum frictional resistance. After extracting the SIV values, the magnitude filter was applied manually, considering that the glacier flow rate should not be abrupt but gradual (Scherler et al., 2008;Garg et al., 2017).
February 2022 | Volume 10 | Article 788359 Surface Elevation Change and Mass Balance. To estimate surface elevation change/"dh", the remote sensing-based geodetic method was adopted for 2000-2017. Before deriving the "dh", certain horizontal and vertical adjustments were made in the secondary DEM (ASTER DMO) with respect to the master (SRTM) DEM. In this regard, the secondary DEM was initially corrected for X-Y/horizontal shift (Berthier et al., 2007(Berthier et al., , 2016 to ensure maximum congruence between the DEMs. Following this, the secondary DEM was corrected for vertical biases possibly influencing the "dh". Prior to the implementation, the following masks were built, on which the corrections were applied: 1) outlier mask-the elevation difference of ≥±100 were eliminated, 2) slope mask-gentle (<4°) to steep (>45°) sloping areas were excluded (Berthier et al., 2016;Garg et al., 2019), 3) cloud and glacier mask-clouds and glacier boundaries (digitized manually) were removed.
Thereafter, the secondary DEM was corrected for: 1) sensorspecific biases, i.e., along/cross-track, which might have aroused due to distortions in sensor geometry (pitch/yaw/roll). To account for these errors, the coordinate system was rotated from X & Y to cross & along-track directions, respectively. The respective elevation difference was computed and added to the ASTER DEMs (Nuth and Kääb, 2011). 2) elevationdependent biases were minimized by fitting a higher-order polynomial between elevation and elevation difference and subtracted from the secondary DEM (Berthier et al., 2006;Schiefer et al., 2007;Bolch et al., 2008). 3) slope-dependent biases might be frequent in steeper terrain and have been corrected by establishing a relation between slope and elevation difference and fitting it into the ASTER DEM. Finally, the "dh" of the range ±100 was extracted, as the values less or greater than 100 are considered outliers . Another possible error may be introduced due to the penetration of the C-band radar waves in the glacier surface, which may vary from 0 to 10 m depending on the local climatic conditions (Schiefer et al., 2007). Due to this, the SRTM DEM underestimates the surface elevation values at higher and overestimates at lower elevations (Berthier et al., 2006;Schiefer et al., 2007). Previous studies have widely recognized this potential source of error and need to be accounted for while estimating the glacier elevation changes Vijay and Braun, 2016;Bhushan et al., 2018). Also, the depth of penetration varies for the different glacier classes, namely snow/firn, ice and debris . In this study, penetration depth values of 2.3 ± .6, 1.7 ± .6, and .4 ± .8 m for snow, clean ice, and debris-covered ice, respectively, were used separately, as estimated by Kääb et al. (2012).
The "dh" was computed by subtracting the old DEM from the newer DEM (Eq. 1). The "dh" was estimated along the central flow line, which has maximum changes in the elevation.
Another aspect that needs due consideration is the difference in the acquisition months of the DEM being compared (Gardelle et al., 2013;Zhou et al., 2018). The SRTM DEM was acquired during early February and ASTER DEM during September.
Hence, the seasonality correction is required to account for the possible mass change during these 4-5 months. For this, a winter accumulation rate of ± .15 m w.e (equivalent to .17 m) per winter month was used (Gardelle et al., 2013).
The obtained "dh" is further used to estimate mass balance by multiplying the "dh" with average snow/ice density, i.e., 850 kg m −3 (Huss, 2013).

Estimation of Uncertainties
The present study assessed multidimensional glacier parameters involving multisource data, hence, liable to introducing errors (Paul et al., 2013). Therefore, to have confidence in our results, estimation of uncertainties is necessary (Racoviteanu et al., 2008) and was carried out using well-established methods. While mapping the lake or glacier boundaries, the possible error sources are locational and interpretational (Bhambri et al., 2013). The locational errors were rectified through coregistration of satellite images with minimum possible root mean square error (Supplementary Table S1). Meanwhile, the interpretational errors were rectified/minimized using high-resolution satellite data, i.e., LISS-4 and Sentinel MSI imageries for 2017 and 2018, respectively. Comparing our results with these high-resolution satellite data suggest the length and area mapping uncertainties of .5 and 3.8% (Liss-4) and .9 and 4.6% (Sentinel), respectively. Further, area mapping uncertainty has also been assessed using the buffer method by Granshaw and Fountain (2006), mentioned in detail in Shukla et al. (2020b).
Further, the glacier area and length change uncertainties were also estimated ( Table 1) utilizing the formula proposed by Hall et al. (2003): where "a" and "b" are the pixel resolution of images 1 and 2, respectively and 'σ' is the registration error.
where U T is the terminus uncertainty and "x" is the spatial resolution of the sensor.
Since the debris extents were delineated within the respective glacier boundaries, the proportionate errors are likely to have propagated in debris cover estimates as per Garg et al. (2017). The SLA uncertainty is equal to half the spatial resolution of the satellite images used (for creating the buffer on either side of the snow line demarcating the snow and ice facies) in the horizontal direction and ±17 m (vertical accuracy of SRTM DEM) in the vertical direction.
Surface ice velocity estimations may be subjected to errors arising primarily due to 1) poor correlation of the images, 2) data quality 3) errors during orthorectification. To improve the correlation, the signal-to-noise ratio less than .9 were eliminated. Furthermore, peak ablation and cloud-free images were used to minimize the errors introduced due to poor quality images (cloud and seasonal snow-covered). Also, the Landsat scenes were used here for SIV assessment, which are Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 788359 orthorectified, eliminating the chances of error (Gantayat et al., 2014;Bhushan et al., 2018). Uncertainty in glacier dynamics is estimated over the stable ground, which, ideally, should not have any movement ( Table 2). Uncertainty in mass balance (Δ B) is obtained by combining all the potential sources of error, which may have accumulated during its computation. It involves the input variables of surface elevation change, the density of ice and water (Bhushan et al., 2018) and is estimated to be .06 and .02 m w.e.a −1 for Dulung and Chilung Glaciers, respectively, during 2000-17.
Stage 3 (2008-2018): represents the continued increment in the lake expansion at the rate of .02 ± .01 km 2 a −1 .

Glacier Morphological Changes
To understand the impact of lake dynamics on glacier response, multidimensional glacier parameters, i.e., area, length, SLA, SIV, surface elevation change and mass balance, were analyzed from 1971 to 2018.

Frontal Retreat and Area Changes
The snout position of the glaciers monitored for 47 years (Figure 3) reveals an overall retreat of the Dulung Glacier by 1,327 m (retreat rate: 32 ± .7 ma −1 ). Meanwhile, the Chilung Glacier in the adjacent valley has retreated by 690 m (retreat rate: 17 ± .7 ma −1 ). Decadal analysis reveals the highest retreat during Stage 3, followed by Stage 2 and minimum during Stage 1 ( Table 3).

Snow Line Altitude Changes
The end-of-summer SLA obtained during 1977-2018 shows an average altitude of 5,021 ± 70 and 5,073 ± 90 m asl for Chilung and Dulung Glaciers, respectively ( Figure 4). Apart from the higher average SLA, a slightly increased overall rise in SLA (absolute increase: 204 m) has also been observed for the Dulung Glacier (Figure 4), compared to the Chilung Glacier (absolute increase: 191 m). On the decadal scale, maximum upshift in the SLA has been observed during Stage 3, followed by Stage 2 and minimum during Stage 1 ( Table 3).

Response of Lake-Terminating Dulung Glacier and Synchronicity With Other Western Himalayan Glaciers
The overall response of the Dulung Glacier to climate change suggest its degeneration (mass loss: −.47 ± .06 m w.e.a −1 ) during 2000-2017. The glacier has also shown an upshift in the SLA (absolute increase: 204 m) during 1977-2018, suggesting its melting during the study period. This overall melting has resulted in the glacier shrinkage (.3 ± .001% a −1 ) during 1971-2018 and reduction in velocity (by 2 ma −1 Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 788359 equivalent to 16%) during 1993-2018. The overall glacier degeneration has also altered its frontal dynamics (retreat: 1,327 m; retreat rate:~32 ± .7 ma −1 ). These results are indicative of the overall negative health of the Dulung Glacier. Further, to comprehensively understand the Dulung Glaciers' state and its synchronicity with the regional trend, the glaciers' response is evaluated with respect to the other western Himalayan glaciers (discussed in this section). Compared to the Dulung, the neighboring land-terminating Chilung Glacier has experienced less mass wastage (−.28 ± .02 m w.e.a −1 ) during 2000-2017. The mass balance estimate of Dulung Glacier is comparable with the nearby Pensilungpa and Lalung Glaciers of the Suru sub-basin (sharing similar  climatic conditions as per (Shukla et al., 2020b)), showing an average mass loss of −.41 m w.e.a −1 , respectively during 1999-2016 (Bhushan et al., 2018). Comparing the Dulung Glaciers' mass loss with other lake terminating glaciers of the western Himalaya (Drung Drang and Samudra Tapu), a fair coherence is also observed (Mukherjee et al., 2018;Ramsankaran et al., 2019; Figure 7A; Supplementary  Table S3).
As a result of mass loss and shrinkage, the Himalayan glaciers accumulate debris cover (dominantly in the lower ablation region), causing their slowdown (Dehecq et al., 2019). The SIV estimates obtained from the previous investigations in the western Himalaya suggest an average velocity of 18 ± 13 ma −1 , with a slowdown rate of 2.1 ± 1.4% a −1 during 1992-2018 ( Figure 7B; Supplementary Table S3). Amongst these studies, the lake-terminating glaciers have shown higher average velocity (22 ± 19 ma −1 ) and slowdown (3 ± 2.4% a −1 ), compared to the land-terminating ones (average velocity: 14 ± 8 ma −1 , slowdown: 1.9 ± .6% a −1 ). The velocity change rates of these glaciers are also in accordance with our estimates of Dulung, showing average flow rates of 14.6 ma −1 (slowdown: 16%), almost double that of Chilung Glacier (average velocity: 9 ma −1 , slowdown: 8%) during 1993-2018. These studies support the findings, which suggest enhanced flow rates of the lake-terminating glaciers compared to their land-terminating counterparts (Pronk et al., 2021).
Other than mass, area and velocity changes, Dulung Glacier has also retreated during 1971-2018 (Figure 3), like the most other Himalayan glaciers (Kamp et al., 2011;Bhambri et al., 2013;Kulkarni and Karyakarte, 2014;Shukla and Qadir, 2016). On the regional scale, the western Himalayan glaciers have retreated at an average rate of 11.6 ± 8 ma −1 during 1962-2019 ( Figure 7B; Supplementary Table S3). Amongst these western Himalayan glaciers, the land-terminating ones have shown an average retreat rate of 10.4 ± 7 ma −1 . Whereas, the lake-terminating ones (Samudra Tapu, Gepang Gath, Panchi Nala and Drung Drang) have retreated significantly at an average rate of 18.6 ± 10 ma −1 . In agreement with these findings, a lower retreat rate was observed for Chilung (~17 ± .7 ma −1 ), the land-terminating glacier, compared to the Dulung Glacier (~32 ± .7 ma −1 ). Mir et al. (2018) and Kamp et al. (2011) suggested an average retreat rate of 23.3 and 26.6 ma −1 for the Dulung Glacier during 1962-2015 and 1975-2003, respectively. The reported difference in the estimates may be attributed to the different time durations considered by the respective studies. Furthermore, the retreat rate of Dulung has not only exceeded that of Chilung but with respect to all other western Himalayan glaciers except Gepang Gath ( Figure 7B; Supplementary Table S3).
These studies, thus, reveal the crucial role of PGL in influencing the calving rates (Figure 3), leading to the rapid retreat of the lake-terminating glaciers compared to their counterparts (R'ohl, 2006;Maurer et al., 2019;Shugar et al., 2020). Similarly, the enhanced degeneration of the Dulung Glacier would have led to the formation of a lake at its snout. However, this cannot be the sole reason, therefore, we investigated other possible factors that could have possibly contributed in the lake evolution.

Evolution of Proglacial Lake: Forming Factors and Dynamics
The Dulung PGL existed as a small pond in 1971, suggesting its formation during or before the 1970s. However, before understanding the impact of PGL on glacier response, certain queries regarding its formation and dynamics need to be answered.
Often the origin of the PGLs is closely linked to initial occurrences as supraglacial ponds, with the latter acting as the precursor for its (PGL) formation. Therefore, similar to supraglacial, the PGLs' evolution should also depend upon factors such as glacier topography and its (mother glacier hosting the lake) response to climate change. Besides, valley morphology should also have an important role in PGL formation, as it is the valley that ultimately accommodates the lake. Amongst these factors, the enhanced melting of the mother glacier is an obvious prerequisite for forming any type of glacial lake. This is because of the reason that the impoundment of water anywhere on or near the glacier suggests that the melt produced is in excess of the system's capacity to drain it. These factors are thus examined in detail to understand their role in the lakes' formation and subsequent expansion. Also, these factors are assessed for both, i.e., the Dulung and Chilung Glaciers, to comprehend the reason for lake formation only on Dulung.
Glacier surface gradient plays an important role in controlling the glacier dynamics, which, in turn, has a profound influence on glacial lake development (Salerno et al., 2012). The assessment of glacier topography suggests that the Dulung Glacier has attained a gentler mean slope of 18°compared to the Chilung Glacier (20°).
Meanwhile, the PGL has a mean slope of 6°. Besides, the difference between the accumulation and the ablation region gradients of Dulung is also more (21.1°/12.6°) compared to the Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 788359 11 Chilung (19.2°/16.8°). This implies that, unlike Chilung, the steeper slope in the accumulation and gentler in the ablation portions of the Dulung Glacier must have been one of the important factors in the formation of the PGL. Salerno et al. (2012) also suggested that the difference in the surface gradients between the two zones and not the overall slope is the key responsive factor contributing to the lake formation.
The longitudinal valley profile of the glaciers has been studied in three main sections (Figure 8). In the first section, i.e., 500 m from the snout downglacier, the Dulung Glacier exhibits a steeper slope of 5.4°than Chilung (4°). In the second section (500-1,000 m downglacier), the Dulung valley shows a flat profile, restricting the lake. This flat region ends into a mound (elevation of~4,357 m) at a distance of~1.8 km from the snout of Dulung. This mound is the moraine dam, which is~9 m high and~195 ± 5 m wide and seen as a loose and unconsolidated material in the Google Earth Imagery. On the other hand, Chilung Glacier valley is gently sloping for the rest of its length. It is important to note that the Chilung valley has not shown any mound in this section; rather, it follows a steep (less than the first section) profile downstream. Thereafter till 4,500 m, both the glaciers follow almost similar profiles. In the third section (4,500 m and further downglacier), Dulung valley differed from Chilung in having a lower elevation and a steeper slope (Figure 8).
The glacier valleys have a typical "U" shape, with a wide and flat floor and steep sides due to the intense erosion of the floor than sides. Investigation of the transverse profile ( Figure 1) reveals a wider (226 ± 5 m) and narrower (87 ± 5 m) valley for the initial distance of~4 km of the Dulung and Chilung FIGURE 8 | Longitudinal profiles of the Dulung (red) and Chilung (blue) glacier valleys. Dulung has an initial steeper section followed by a flat region that ends in a mound. The proglacial lake has formed behind this mound. Glaciers' snout downglacier, respectively. After this distance in downward direction, the trend is opposite, with a narrower (109 ± 5 m) and wider (104 ± 5 m) valley for Dulung and Chilung, respectively ( Figure 9). From these investigations of the glacier topography and valley morphology, it is evident that the relatively higher difference in the surface gradient between the accumulation and the ablation zones of the Dulung Glacier (compared to the Chilung Glacier) would have facilitated the formation of the PGL. Meanwhile, the Dulung Glacier valley, which is initially wider (for a distance of about 400 m from the glacier snout) and narrower down later, must have contributed to the accommodation and sustenance of the lake. Besides accommodating the PGL, valley morphology (shape, slope, erosion capacity, avalanche zones) also has a crucial role in determining its hazard potential (Schneider et al., 2014;Patel et al., 2017).
Once formed, these lakes grow over time due to the differential thermal regimes of ice and warm water at the glacier's terminus, thus, altering the glaciers' response to climate change. Thus, to understand the expansion of the PGL, we assessed the long-term end-of-summer SLA (or ELA), which gives a close approximation of the mass balance condition of the studied glaciers (specifically during 1977-2000, when the mass balance data is not available). During 1977-2000, a higher rise in the SLA was observed for the Dulung (increased by 13 ma −1 ) compared to the Chilung Glacier (increased by 11 ma −1 ). The relatively higher increase in the Dulung Glacier's SLA (compared to the Chilung Glacier) indicates its enhanced melting, which would have been a contributing factor in initiating the PGL formation. Further, this trend (of glacier melting) continued post-2000 and is reflected from the glaciers' mass balance condition depicting an increased mass loss of the Dulung (−.47 ± .06 m w.e.a −1 ) glacier compared to the Chilung (−.28 ± .02 m w.e.a −1 ) glacier. Also, the variations in SLA can very well be correlated with the lake area changes. Initially, the slow growth rate of the lake area (.008 km 2 a −1 ) during 1977-1993 may be attributed to the slight upshift in the SLA during this period (Section 4.2.2). Meanwhile, in the subsequent years, i.e., during 1993years, i.e., during -2008years, i.e., during and 2009years, i.e., during -2018 drastically by about 12 ma −1 and 25 ma −1 , respectively, for Dulung compared to the Chilung Glacier (9 ma −1 and 19 ma −1 , respectively). This upshift is reflected from the rate of increase in lake area during the respective time frames, i.e., .01 km 2 a −1 and .02 km 2 a −1 . Also, a good positive correlation (R 2 = .69) has been observed between the SLA changes and the rate of increase in the lake area. All these evidences reveal a synchronized relation between the Dulung Glaciers' melting and the glacial lake dynamics. All these inferences are in agreement that the initiation, formation and enlargement of the PGLs have a much closer and direct relationship with volumetric changes of the glacier, as they are ultimately fed from meltwater discharged from the mother glaciers (Yao et al., 2010).

Impact of Proglacial Lake on Glacier Stability
Glacier dynamics influence the PGL formation and vice versa. To understand the impact of the lake on its mother glacier, i.e., Dulung, its response was compared to the adjacent landterminating, Chilung Glacier. The investigation reveals that both the glaciers have degenerated during 1971-2018 (Table 3; Figures 5, 6). All these evidences show the sensitivity of both the Dulung and Chilung Glaciers towards climate change (Shukla et al., 2020b).
Besides the overall degenerative responses of these glaciers, their degeneration rates have also varied, with Dulung being more responsive towards climate change than the Chilung Glacier. Despite sharing similar geographic and climate settings (Shukla et al., 2020a), the Dulung Glacier has shown an increased mass loss, nearly twice that of the Chilung from 2000-2017. Meanwhile, the SLA change also reveals a higher increase for Dulung compared to the Chilung Glacier during 1977-2018 ( Figure 4). Other than the overall changes, the variability in glacier response is also observed on different time scales. The glaciers have degenerated more during Stage 2 and 3 compared to the Stage 1 (Section 4). This heterogeneity is also clearly visible in the lake area changes, showing an increase stage 2 onwards (Section 5.1). All these evidences reveal that the PGL has significantly modified the response of the Dulung Glacier over the last four decades. Besides the glacier-wide negative mass balance of the glaciers, the zonal analysis reveals maximum mass loss in the snout region of the Dulung in contrast to the Chilung Glacier (maximum mass loss further up in the ablation zone) (Figure 6). This may be attributed to the different mechanisms operating for lake and land-terminating glaciers. The accelerated melting at the snout portion of the former is owing to the transmission of thermal energy due to the increased absorption of solar energy at the ice-water interface, leading to enhanced terminal melting and, thus, loss in glacier area (King et al., 2018;Liu et al., 2020).
The shrinkage and mass loss of the glaciers have resulted in the increment of debris coverage of the glaciers (average 36%), predominantly in the snout portion. However, the Dulung Glacier has accumulated less debris content (33%) compared to the Chilung Glacier (39%) during the period 1971-2017 (Shukla et al., 2020b). Given the enhanced melt rates of the Dulung Glacier, it should have accumulated more debris cover at its snout (Ali et al., 2017); however, it is not so. This might have occurred due to the magnified mass loss at the snout-lake interface of the Dulung Glacier (Figure 6), thus transferring the debris cover to the moraine dam (surrounding the lake) rather than its accumulation in the lower ablation portion .
Consequent to the observed mass loss, both the glaciers have slowed down during the 1993-2018 period ( Figure 5), however, at variable rates (Dulung: 16%, Chilung: 8%). Despite the overall slowdown, a periodic acceleration was also noticed during 1993-2000 (Section 4.2.3). Relatively higher flow rates are observed for Dulung (14.6 ± 4 ma −1 ), almost double that of the Chilung (9 ± 4 ma −1 ) glacier. The velocity has also varied spatially, with higher flow rates observed at the Dulung Glacier's snout (14.4 ma −1 ), reducing further to 12.3 ± 3 ma −1 in the lower ablation zone ( Figure 5). Compared to Dulung, the Chilung Glacier showed nearly half the velocity rates in the snout region (5.9 ma −1 , Figure 5), increasing to 8.8 ± 2 ma −1 further up in the Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 788359 lower ablation zone. Meanwhile, a maximum velocity is observed near the ELA for both the glaciers (Dulung: 16.2 ± 1 ma− 1 ; Chilung: 14.2 ± 3 ma −1 ), further reducing in the accumulation zone (Dulung: 14.1 ± 4 ma −1 ; Chilung: 10.3 ± 5 ma −1 ). The heterogeneity in the velocity profiles is mainly noticed in the snout portion of the glaciers, which can be attributed to the dynamic thinning of the lake-terminating glaciers, also reported by Patel et al. (2021); Pronk et al. (2021). It occurs as a consequence of the reduction in the effective pressure, which may further reduce the compressive stress in the terminus region, leading to dynamic thinning of the glacier (Benn et al., 2012;Liu et al., 2020). Whereas the dynamics of Chilung Glacier differs from the Dulung Glacier, which is not influenced by any glacier lake and similar to other glaciers, flow according to the topography. The glaciers have also retreated significantly during 1971-2018 ( Figure 3). Compared to the Chilung (17 ± .7 ma −1 ), the Dulung Glacier has retreated at an accelerated rate of 32 ± .7 ma −1 during 1971-2018. Besides, the overall retreat of the glacier during 1971-2018, the retreat rates have also varied on the decadal scale. Decadal analysis reveals that Dulung's retreat rate has enhanced manifolds compared to the Chilung, showing only a slight increase during different periods (Section 4.2.1). Also, the regression coefficient presented a very good correlation (r 2 = .9) between the lake area and the length of Dulung, which indicates the role of the PGL in modifying the Dulung Glaciers' retreat.
All these evidences thus prove that the presence of PGL at the Dulung Glaciers' snout significantly influenced its stability as compared to the Chilung Glacier.

CONCLUSIONS
This study investigated the evolution of a proglacial lake (PGL) associated with the Dulung Glacier, Suru sub-basin, western Himalaya and its impact on the glacier's response during 1971-2018. Further, to better comprehend the impact of PGL dynamics on the glacier stability, the adjacent landterminating Chilung Glacier was also assessed. The PGL has dynamically evolved since 1970, favored by topography, valley morphology and glacier morphological changes. The expansion rates of the PGL almost doubled during the investigation periods, i.e., for stage 2 (1993-2008) and 3 (2008-2018) compared to stage 1 (1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993). The continuous PGL expansion during 1971-2018 has induced enhanced degeneration of the Dulung Glacier, which shows 1.5 times higher shrinkage and mass loss rates and twice higher retreat rate and velocity reduction than the land-terminating Chilung Glacier. A clear heterogeneity in the glaciers velocity profiles is also noticed in the snout portion (higher flow rates in case of Dulung Glacier), attributed to the dynamic thinning of the lake-terminating glaciers. The accelerated degeneration of the Dulung Glacier is observed stage 2 onwards, which is in sync with the increased lake expansion rate, implying the control of PGL in influencing the Dulung Glaciers' behavior. The degeneration rate estimated for the Dulung Glacier accords well with other lake-terminating glaciers of western Himalaya and is observed to be higher compared to the land-terminating glaciers, suggesting the highly responsive nature of the lake-terminating glaciers. With the continued melting of the Dulung Glacier, its PGL has a higher potential to grow further. Considering these changes, regular monitoring of the Dulung Glacier and its PGL is recommended. In the future, the GLOF risk assessment should also be performed to establish early warning systems in the susceptible areas in case of any future glacial hazard.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SG and AS designed the study and led the writing of the manuscript. SG executed the experiments and prepared the figures. AS supervised the study and improved the data analysis and interpretation. PG, BY, and US thoroughly read the paper and contributed to the interpretation of results.