1. Introduction
As the largest mid-latitude repository of glaciers, High Mountain Asia (HMA) has experienced accelerated ice mass loss under climate warming, with a mass change rate of −20.1 Gt a−1 from 2000 and 2019 (Hugonnet and others, Reference Hugonnet2021). This melting has driven substantial meltwater runoff, promoting the formation and expansion of glacial lakes. Satellite-based inventories document an 11% increase in the number of glacial lake (n = 27205 to 30121) and 15% expansion in total area (1806.5 ± 2.1 and 2080.1 ± 2.3 km2) across HMA between 1990 and 2018 (Wang and others, Reference Wang2020). These proglacial water bodies now exert fundamental controls on glacier flow dynamics (Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021) and regional mass loss patterns (Tsutaki and others, Reference Tsutaki2019).
Proglacial lakes form behind the terminal moraine dams, originating from the coalescence of supraglacial ponds or lakes near glacier termini (Quincey and others, Reference Quincey2007; Carrivick and Tweed, Reference Carrivick and Tweed2013). These lakes amplify glacier mass loss rates and modify ice flow dynamics through two synergistic mechanisms: (1) subaqueous melting caused by thermal undercutting (Truffer and Motyka, Reference Truffer and Motyka2016) and the calving process (Benn and others, Reference Benn, Warren and Mottram2007a) amplify frontal retreat; and (2) basal lubrication from reduced effective pressure elevates ice flow (Sugiyama and others, Reference Sugiyama2011; Sutherland others, Reference Sutherland2020). Observations confirm these impacts: in the Everest region, lake-terminating glaciers lost mass 32% faster than land-terminating types during 2000–2015 (King and others, Reference King, Quincey, Carrivick and Rowan2017). Across HMA, the mass loss rates of lake-terminating glaciers exceeded regional averages by 18–97% (Brun and others, Reference Brun2019), with recent evidence suggesting substantial underestimation of subaqueous melting impacts (Zhang and others, Reference Zhang2023). Correspondingly, lake-terminating glaciers exhibit two- to threefold higher velocities than their land-terminating counterparts (Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021). Despite extensive research highlighting the profound impacts of proglacial lakes on glaciers, lake-induced effects propagate only through the lowermost ∼30% of glacier length, with significant behavioral changes in glaciers occurring during advanced stages of lake evolution (Scoffield and others, Reference Scoffield2024). This glacier–lake coupling establishes a positive feedback loop: lake expansion accelerates retreat, which in turn facilitates further lake growth. The evolving glacier–lake contact, driven by progressive terminus retreat, is anticipated to amplify spatial heterogeneity in glacial response behaviors. Crucially, understanding glacier behavior under changing glacier–lake interaction not only reveals the nonlinear characteristics of glacial responses to climate warming but also provides mechanistic support for water resource (Warren and Kirkbride, Reference Warren and Kirkbride2003) security assessments and disaster risk management.
The Himalaya and Southeastern Tibet host the highest concentration of lake-terminating glaciers within the HMA (Wang and others, Reference Wang2020; Chen and others, Reference Chen2021), with glaciers in these areas exhibiting the most negative mass balance in recent decades (Brun and others, Reference Brun, Berthier, Wagnon, Kaab and Treichler2017; Brun and others, Reference Brun2019; Wang and others, Reference Wang2020). With the continuing increase in both the number and size of proglacial lakes, investigating the distribution and changes of lake-terminating glaciers in the two regions is crucial for assessing and projecting glaciers’ responses to climate change. Several studies have examined the spatiotemporal distribution of glaciers and glacial lakes in these regions. For example, Dou and others, (Reference Dou2023) who updated the glacier–lake inventory in Southeastern Tibet and analyzed its response to climate change, and Kumar and others, (Reference Kumar, Sharma and Bahuguna2024) who monitored the spatiotemporal dynamics of glacial lakes in the Himachal Pradesh, Himalaya. However, systematic studies addressing the spatiotemporal distribution of the lake-terminating glacier–proglacial lake system and tracking changes in glacier–lake contact remain limited.
Therefore, this study focuses on lake-terminating glaciers with the following objectives: (1) to investigate their spatial distribution and spatiotemporal variability from 1990 to 2020 using multisource optical remote sensing datasets in the Himalaya and Southeastern Tibet; (2) to track the evolution of glacier-lake contact; and (3) to analyze the dynamics differentiation among lake-terminating glaciers.
2. The Himalaya and Southeastern Tibet regions
The Himalaya extends as an arcuate belt from Nanga Parvata in the west (about 8125 m a.s.l.) to Namcha Barwa in the east (about 7782 m a.s.l.). This cryospheric hotspot sustains 18874 glaciers, which are classified into temperate and continental cold glaciers. The region’s climate system is governed by the Indian Ocean monsoon and secondary contributions from the westerlies. The region experiences wet seasons during both winter and summer, with the summer monsoon affecting the Eastern Himalaya and the winter monsoon primarily impacting the Western Himalaya (Bookhagen and Burbank, Reference Bookhagen and Burbank2010). The extreme topography of the Himalaya blocks monsoonal air, resulting in relatively dry conditions on the northern slopes (Ageta and Higuchi, Reference Ageta and Higuchi1984).
The Southeastern Tibet (ST) encompasses the eastern regions of Nyainqentangla, Gangri Garpo, Goikarla Rigyu, and the western section of the Hengduan mountains. The region has an average elevation exceeding 3900 m a.s.l., with some areas rising above 7000 m a.s.l. Its climate system is governed by the Indian monsoon and the East Asian monsoon during the summer and spring, bringing abundant precipitation accompanied by warm conditions. The rainy season typically lasts from May to October each year. This region contains 7472 temperate glaciers, most of which are distributed above 3000 m a.s.l. and concentrated in the eastern Nyainqentanglha Mountains.
We divided the Himalaya and Southeastern Tibet into four subregions (Fig. 1a): Western Himalaya (WH), Central Himalaya (CH), and Eastern Himalaya (EH) and Southeastern Tibet (ST, including eastern of Nyainqentangla, Gangri Garpo, Goikarla Rigyu and western of Hengduan). In total, the study area contains 26346 glaciers with a combined area of 25900.5 km2 (RGI6.0).
Figure 1. (a) Lake-terminating glacier distribution overview; (b) and (c) Numbers and aspect distributions of three different type of lake-terminating glaciers in four subregions. Type 1 change: Glaciers that remain in contact with proglacial lakes. Type 2 change: Glaciers with newly formed proglacial lakes. Type 3 change: Glaciers detached from proglacial lakes. WH: Western Himalaya. CH: Central Himalaya. EH: Eastern Himalaya. ST: Southeastern Tibet.
3. Data and methods
3.1. Proglacial lakes inventory
The methodological framework for constructing a lake-terminating glacier inventory requires accurate delineation of proglacial lakes as a prerequisite step. Proglacial lakes were automatically mapped using a workflow implemented in Google Earth Engine (GEE).
Landsat (TM/ OLI) imagery (Table 1) was employed to map proglacial lakes due to its five-decade temporal coverage (1972–present), suitable spatial resolution (30 m), global availability and open-access policy. The Landsat products provided by GEE were preprocessed, including radiometric calibration, atmospheric correction and geometric correction. To minimize the influence of ice and seasonal snow cover, as well as annual changes of lake-terminating glaciers and proglacial lakes, we selected images acquired during summer and autumn (July to November). Two-time windows, 1990 ± 2 and 2020 ± 1, were chosen to reveal the present and past state of glacial lakes. This study focuses on identifying glacier–contact lakes, therefore the investigation was limited to a 2 km buffer zone around glaciers. In total, 841 Landsat Thematic Mapper (TM) scenes were used for the 1990 period (Figures S1 and S2) and 1133 Landsat Operational Land Imager (OLI) scenes for the 2020 period (Figures S3 and S4). To reduce the effect of cloud coverage, we applied the C Function of Mask (CFMask) algorithm (Foga and others, Reference Foga2017) to detect and remove the cloud before mosaicking the imagery.
Table 1. Overview of the dataset used in this study
An automated mapping algorithm based on hierarchical image segmentation and terrain analysis was employed to delineate glacial lake extents (Li and Sheng, Reference Li and Sheng2012; Zhang and others, Reference Zhang2017). The Normalized Difference Water Index (NDWI; [(GREEN − NIR)/(GREEN + NIR)]) was employed to extract proglacial lakes. To reduce disturbances from mountain shadows, a slope threshold of <20° and a shaded relief threshold of >0.25 were applied (Zheng and others, Reference Zheng2021). An example is given in Figure S5. Previous studies exhibit variability in minimum area threshold for glacial lake identification. For example, Wang and others, (Reference Wang2020) set the threshold to 0.0054 km2, Chen and others, (Reference Chen2021) set it to 0.0081 km2, Luo and others, (Reference Luo, Zhang, Chen and Xu2020) set it to 0.0036 km2 and Li and others, (Reference Li, Shangguan and Anjum2020) set it to 0.01 km2. Smaller lake areas generally lead to greater error identification (Salerno and others, Reference Salerno2012). To identify lake-terminating glaciers as accurately as possible, this study adopts a minimum lake area threshold of 0.0036 km2 (at least four pixels), following Luo and others, (Reference Luo, Zhang, Chen and Xu2020). Proglacial lake datasets were cross-referenced with RGI6.0 glacier inventory to identify lake-terminating glaciers. Finally, careful visual inspection and manual re-editing were performed to correct misclassified glacial lakes and identify lake-terminating glaciers, using Landsat images, Planet Labs (Table 1), online maps (Google Earth, Esri basemap, etc.) and other glacial lake datasets (Wang and others, Reference Wang2020; Chen and others, Reference Chen2021; Zheng and others, Reference Zheng2021).
The uncertainty (
$\lambda $) and relative error (
${E_g}$) in glacier area was estimated to using the equation (Bolch and others, Reference Bolch, Menounos and Wheate2010):
\begin{equation}\lambda = N \times \frac{{{G^2}}}{2} \end{equation}
\begin{equation}{E_g} = \frac{\lambda }{S} \times 100\% \end{equation} where 
$N$ is the total count of pixels along the outline of ice coverage, 
$G$ is the spatial resolution of the images used and 
$S$ is the glacier area. The uncertainty (
$\delta $) and relative error (
${E_l}$) of glacial lake area was estimated using the equation (Hanshaw and Bookhagen, Reference Hanshaw and Bookhagen2014):
\begin{equation}\delta = \frac{P}{G} \times \frac{{{G^2}}}{2} \times 0.6872 \end{equation}
\begin{equation}{E_l} = \frac{\delta }{A} \times 100\% \end{equation} where 
$P$ is the perimeter of the glacial lake and 
$A$ is the glacial lake area.
Then the accumulation of the study region can be calculated using the following equation based on error propagation theory (Wang and others, Reference Wang2020):
\begin{equation}{E_T} = \sqrt {\mathop \sum \limits_{i = 1}^n a_i^2} \end{equation} where 
${E_T}$ is the area error of the entire study region or subregions, 
$i$ is the glacier or lake of no.
${\text{ }}i$ in the entire study region or subregion and 
$a$ is the error area. The error for calculation of changes in glacier and glacial lake area (
$ds$) is obtained by:
\begin{equation}ds = \sqrt {{S_1}^2 + {S_2}^2} \end{equation} where 
${S_1}$ and 
${S_2}$ are the glacier or glacial lake areas at the beginning and end of the period.
3.2. Lake-terminating glacier identification and classification
In this study, lake-terminating glaciers are defined as glacier forming proglacial lakes along the direction of ice flow. Their classification was based on the glacier–lake contact change from 1990 to 2020. The identification of glacier–lake contact followed a two-stage procedure: (1) Preliminary screening was conducted using spatial intersection analysis between glacier boundaries and ice–marginal lakes, with a 500 m buffer tolerance. (2) Manual verification employed period-specific criteria. For the 2020 assessments, multi-source moderate- to high-resolution imagery (Planet Labs, Landsat, Google Earth, Esri basemaps, etc.) was used. Glacier–lake contact was confirmed when proglacial lakes overlapped with glacier terminus and exhibited diagnostic features such as terminal ice cliffs or transverse crevasses perpendicular to the flow direction. Evaluations for 1990 faced inherent uncertainties due to the limited spatial resolution of Landsat products (30 m), particularly for smaller glaciers where boundary delineation errors increased with decreasing glacier size. To address this, temporal cross-validation was applied: glaciers with ambiguous glacier–lake contacts were classified as interactive if sequential imagery from 1990 to 2020 showed lake expansion toward glacier termini. Lake-terminating glaciers were further categorized into three types: persistent glacier–lake connectivity (Type 1 change, Fig. 2a), transition from supraglacial lake to proglacial lake (Type 2 change, Fig. 2b) and detachment of proglacial lakes from parent glaciers (Type 3 change, Fig. 2c).
Figure 2. Examples of the three developing stages of lake-terminating glaciers from 1990 to 2020. Group a: Type 1 change (Jiongpu glacier as example), terminus keeping connected with proglacial lake. Group b: Type 2 change (Unnamed Glacier, RGI id: RGI60-15.03150, as example), terminus experienced transition from supraglacial lake to proglacial lake. Group c: Type 3 change (Jialong Co glacier as example), terminus experienced detach from proglacial lake.
3.3. Analysis of glacier changes
Glacier surface elevation change data from Hugonnet and others, (Reference Hugonnet2021) (Table 1) were extracted for glaciers larger than 2 km2. Data points were excluded if the 19-year elevation change exceeded five times the normalized median absolute deviation within each elevation bin (Agarwal and others, Reference Agarwal2023), following the methodology of Hugonnet and others, (Reference Hugonnet2021). Glacier velocity data were obtained from ITS_LIVE (Table 1) for glaciers exceeding 2 km2. Pixels with velocity errors greater than 5 m a−1 were removed in accordance with the method of Dehecq and others, (Reference Dehecq2018).
To facilitate the delineation of glacier ablation zones and the interpretation of surface velocity patterns, we followed the approach of Hugonnet and others, (Reference Hugonnet2021) and used the SRTM DEM to divide each glacier into 10 standardized elevation bins. These bins were then indexed from 0.1 to 1.0 according to their relative distance from the glacier terminus.
Due to higher measurement uncertainties in accumulation zones caused by low image contrast, the analysis of velocity anomaly was confined to the lower half of each glacier (0.1–0.5 elevation bins), which approximates the ablation zone (Dehecq and others, Reference Dehecq2018). Velocity anomaly was calculated as the difference between the annual glacier velocity and the multi-year mean velocity.
Center flowlines were generated using the Open Global Glacier Model (OGGM) (Maussion and others, Reference Maussion2019) with a slightly adapted algorithm from (Kienholz and others, Reference Kienholz, Rich, Arendt and Hock2014), and manually adjusted using Landsat-5/8 satellite imagery and ice flow velocity data from ITS_LIVE to ensure center flowlines terminated at the glacier termini in 1990 and 2020 and followed the main flow tributaries. Two sets of center flowline data, corresponding to the glacier terminus positions in 1990 and 2020, were then used to calculate the retreat of lake-terminating glaciers over this period.
4. Results
4.1. Distribution of lake-terminating glaciers
Our 2020 inventory identified 577 lake-terminating glaciers, representing approximately 2% of the glaciers in the RGI 6.0 database, with a cumulative area of 2561.5 ± 11.8 km2. Among these, 396 glaciers (∼69%) had areas below 2 km2 (Fig. 3a). The Central Himalaya hosted the highest number of lake-terminating glaciers, totaling 215 (∼37%, Fig. 1b), whereas Southeastern Tibet exhibited the largest mean glacier size (6.7 km2), excluding the Yanong Glacier in Gangri Garpo, which covers approximately 179.2 ± 2.1 km2 with a length of ∼19.35 km. From 1990 to 2020, the total glacier area decreased by ∼1.5 ± 3.0 km2, and the glacier terminus retreated by about 0.31 km (Fig. 3a). Lake-terminating glaciers are predominantly located at altitudes of 5000 to 6000 m a.s.l. (Fig. 3b). Regionally, those in the Central Himalaya exhibit the highest average elevation (5669 m a.s.l), in contrast to the Western Himalaya, which shows the lowest mean elevation (5295 m a.s.l, Fig. 3c). Most lake-terminating glaciers are north-facing (Fig. 1c), with 385 (∼66.61%) oriented towards north, northeast, or northwest.
Figure 3. (a) Sub-regional size distribution of lake-terminating glaciers and proglacial lakes in 2020. (b) Hypsometry of lake-terminating glaciers and proglacial lakes in 2020. (c) Elevation distribution of lake-terminating glaciers in four subregions in 2020. (d) The area distribution of different types of glaciers in 2020. Type 1 change: Glaciers that remain in contact with proglacial lakes. Type 2 change: Glaciers with newly formed proglacial lakes. Type 3 change: Glaciers detached from proglacial lakes. In figures (c) and (d), the boxplots represent the interquartile range (IQR), the bars indicate the median and the dots represent the mean. WH: Western Himalaya. CH: Central Himalaya. EH: Eastern Himalaya. ST: Southeastern Tibet.
Observations indicate that 583 proglacial lakes are connected to glaciers, with an average area of 0.28 ± 0.03 km2. Although only 217 proglacial lakes exceed 0.1 km2, they account for 92% of the total proglacial lake area (Fig. 3a). Most lakes (∼64%) are situated at elevations between 4800 and 5500 m a.s.l. (Fig. 3b). Notably, the number of lake-terminating glaciers does not correspond directly to proglacial lakes. Our results revealed that there are three types of glacier–lake contact (Fig. 4): 564 glaciers host a single proglacial lake, 7 glaciers develop multiple proglacial lakes and 2 proglacial lakes interact and receive inflow from multiple glaciers.
Figure 4. Examples showing multiple glaciers contributing to a single proglacial lake, and one glacier connecting with multiple proglacial lakes.
4.2. Change of lake-terminating glaciers between 1990 and 2020
Between 1990 and 2020, the number of lake-terminating glaciers increased by 158, representing a 38% growth. This change resulted from two concurrent processes: the formation of new proglacial lakes (quantitative increase) and glacier–lake detachment events (quantitative decrease). Among the three identified glacier change types, Type 2 change glaciers (transition from supraglacial to proglacial lakes) are the most prevalent (Fig. 1b), totaling 331. These glaciers also exhibit the smallest mean area (2.4 km2, Fig. 3d). In contrast, Type 3 change glaciers (proglacial lakes became detached from their parent glacier) have the largest mean area (7.7 km2). These glaciers area primarily face north-facing, including north, northeast and northwest orientation (Fig. 1c), and show similar elevation distributions, primarily clustered between 5000 and 6000 m a.s.l. (Fig. 3b). Over the study period, lake-terminating glaciers collectively lost 161.2 ± 18.1 km2 in terminus area. The greatest decrease was observed in Type 1 change glaciers (73.8 ± 13.1 km2), followed by Type 2 change glaciers (56.0 ± 10.9 km2) and Type 3 change glaciers (31.4 ± 6.0 km2).
We compared glacier dynamics—including length, surface elevation and flow velocity changes—across different glacier types within each subregion. To reduce uncertainties, glaciers smaller than 2 km2 were excluded from the analysis. The Western Himalaya subregion was omitted due to insufficient sample size (n < 5) meeting the 2 km2 threshold.
4.2.1 Glacier retreat
Across three subregions, Type 2 change glaciers exhibited greater median terminus retreat compared to Type 1 change and Type 3 change glaciers (Fig. 5). Statistically significant differences (Table S1 and S2) were observed in the Central and Eastern Himalaya. In the Central Himalaya, Type 2 change glaciers retreated 1.13 ± 0.24 km, significantly more than that Type 1 change (0.72 ± 0.06 km, p < 0.001) and Type 3 change (0.72 ± 0.09 km, p = 0.003). Similarly, in the Eastern Himalaya, Type 2 change glaciers retreated 1.09 ± 0.06 km, significantly exceeding the retreat of Type 1 change glaciers (0.66 ± 0.08 km, p = 0.003). Regarding relative retreat (Fig. 5), significant differences were observed in the Central Himalaya between Type 1 change glaciers (10.8 ± 0.9%) and both Type 2 change (17.4 ± 3.0%, p = 0.019) and Type 3 change (16.6 ± 1.6%, p = 0.047) glaciers. In the Eastern Himalaya, Type 2 change glacier (28.5 ± 2.8%, p < 0.001) and Type 3 change glaciers (16.8 ± 2.0%, p = 0.036) also exhibited significantly higher relative retreat than Type 1 change glaciers (8.9 ± 1.6%).
Figure 5. The boxplots for glacial terminus retreat (a) and percentage of retreat (b) illustrate the distribution of retreat distances and relative retreat percentages across different glacier types. The interquartile range (IQR) is represented by the box, the median by the central bar and the mean by a dot. Blue, red and green colors correspond to Type 1 change, Type 2 change and Type 3 change glacier, respectively. CH: Central Himalaya. EH: Eastern Himalaya. ST: Southeastern Tibet.
4.2.2 Glacier surface elevation change and ice flow velocities
Statistical analysis revealed no significant differences in glacier surface elevation changes among the three types of lake-terminating glaciers across the four subregions (Table S6). In the Eastern Himalaya, Type 2 change glaciers (−0.85 ± 0.10 m a−1) exhibited more pronounced thinning compared to Type 1 change glacier (−0.63 ± 0.07 m a−1) and Type 3 change glaciers (−0.77 ± 0.07 m a−1), as visually illustrated in Fig. 6.
Figure 6. Boxplot illustrates the distribution of median surface elevation changes within glacier extents for different glacier types across each subregion. The interquartile range (IQR) is represented by the box, the median by the central bar and the mean by a dot. Blue, red and green colors correspond to Type 1 change, Type 2 change and Type 3 change glacier, respectively. CH: Central Himalaya. EH: Eastern Himalaya. ST: Southeastern Tibet.
Glacier velocity analyses revealed distinct elevation distribution patterns during the study period (Fig. 7). In the Central Himalaya, the median peak velocity of Type 1 change glaciers (5.32 ± 1.02 m a−1, bin 0.1) occurs at a lower elevation than that of Type 2 change glacier (3.86 ± 1.40 m a−1, 0.2) and Type 3 change glaciers (3.62 ± 0.43 m a−1, bin 0.2). Similarly, in the Eastern Himalaya, Type 2 change glaciers (4.56 ± 2.84 m a−1, bin 0.1) reach their median peak velocity at a lower elevation than Type 1 change (5.45 ± 2.16 m a−1, bin 0.2) and Type 3 change glaciers (5.56 ± 0.69 m a−1, bin 0.3). In Southeastern Tibet, Type 1 change glaciers (5.73 ± 1.57 m a−1, bin 0.1) also exhibit median peak velocity at lower elevation compared to Type 2 change glacier (7.39 ± 1.93 m a−1, bin 0.2) and Type 3 change glaciers (4.18 ± 0.97 m a−1, bin 0.5).
Figure 7. Median velocity across normalized glacier elevation between 1990 and 2018 in (a) Central Himalaya, (b) Eastern Himalaya and (c) Southeastern Tibet. The interquartile range (IQR) is represented by the box, the median by the central bar. Blue, red and green colors correspond to Type 1 change, Type 2 change and Type 3 change glacier, respectively.
During 1990–2018, all glacier types exhibited negative trends in median velocity anomalies (Fig. 8), except for Type 1 change glaciers in the Eastern Himalaya, which showed a marginal acceleration (0.06 ± 0.01 m a−1 decade−1, p = 0.59). Type 2 change glaciers displayed the most pronounced negative velocity anomaly trends across subregions, with statistical significance. In the Central Himalaya, this rate was −0.34 ± 0.05 m a−1 decade−1 (p < 0.001), while in the Eastern Himalaya, it reached −0.45 ± 0.11 m a−1 decade−1 (p < 0.001). In Southeastern Tibet, Type 2 change glaciers exhibited a rate of −0.38 ± 0.08 m a−1 decade−1 (p < 0.001).
Figure 8. Annual glacier velocity anomalies for different glaciers (1990–2018). Black lines indicate the median anomaly, color bars represent the interquartile range and colored lines depict the linear trend. CH: Central Himalaya. EH: Eastern Himalaya. ST: Southeastern Tibet.
To further examine the impact of changes in glacier–lake contact on glacier terminus dynamics, we selected glaciers larger than 5 km2. This threshold ensured sufficient valid pixels for analysis despite glacier shrinking and the limited spatial resolution of the dataset (240 m). Velocity change trends were then tracked in the terminal 25% region of Type 1 change glaciers (n = 79, Table S7) and Type 3 change glaciers (n = 27, Table S8). Among Type 1 change glaciers, 56 exhibited deceleration trends, with 38 showing statistically significant trends, while 23 exhibited acceleration trends, of which 11 were statistically significant. In contrast, among Type 3 change glaciers, 23 showed deceleration and 4 exhibited acceleration trends, with 18 displaying statistically significant trends—all of which were associated with deceleration. Type 2 change glaciers were excluded from velocity trend analysis due to insufficient temporal resolution of lake formation data.
5. Discussion
5.1. Uncertainties in the lake-terminating glacier dataset and comparison with other lake datasets
The calculated error indicates that the total absolute area error in lake-terminating glaciers was approximately ± 13.1 km2 in 1990 and ± 12.3 km2 in 2020, with mean relative error was ± 12.5% and ± 13.8%, respectively. Individual glacier relative errors varied considerably (1–68%), and exhibited a strong power law relationship with glacier size (
${E_g} = 0.10988 \times {S^{ - 0.42384}}$, 
${R^2}{\text{ }}$= 0.88, Fig. 9a). Parallel analysis of proglacial lakes showed total absolute area errors of ± 0.57 km2 and ± 0.73 km2, with mean relative errors of ± 18.02% and ± 21.10% in 1990 and 2020, respectively. Relative errors for individual lakes ranged from 2% to 57% and similarly followed a significant power law dependency on lake area (
${E_l} = 0.05286 \times {A^{ - 0.45000}}$, 
${R^2}{\text{ }}$= 0.94, Fig. 9b). These analyses revealed an inverse relationship between feature size and measurement precision: both glacier and proglacial lake dimensions were significantly negative correlated with relative area errors (p < 0.01), indicating smaller features were associated with greater uncertainty, consistent with the observed power law dependencies.
Figure 9. Relationships of relative area error against size of glacier (a) and glacial lakes (b).
Among publicly available datasets covering the Himalaya and Southeastern Tibet, the inventory compiled by Zhang and others, (Reference Zhang2023) provides detailed classification of proglacial lakes with temporal coverage matching our study period. Therefore, we conducted a comparative analysis primarily against their dataset. According to their records, the study area contained 512 proglacial lakes (total area: 121.87 km2) in 1990 and 816 lakes (178.51 km2) in 2020. These figures exceed our inventory by 95 lakes (12.10 km2) and 233 lakes (15.65 km2) for the respective years. For lake-terminating glaciers, we selected candidates from the RGI 6.0 dataset using the glacial lake inventory provided by Zhang and others, (Reference Zhang2023) for comparison. Based on their dataset, the number of identified lake-terminating glaciers was 491 in 1990 and 758 in 2020, exceeding our corresponding figures by 72 and 181 glaciers, respectively. In 1990, 331 glaciers were consistently identified as lake-terminating in both datasets. Additionally, 88 glaciers were classified as lake-terminating only in our study, while 160 were identified exclusively Zhang and others, (Reference Zhang2023). Similarly, in 2020, 485 glaciers were jointly recognized as lake-terminating. Meanwhile, 92 glaciers were identified as lake-terminating only in our study and 273 only in Zhang and others, (Reference Zhang2023).
The observed discrepancies are likely attributed to the stringent criteria adopted for lake-terminating glacier identification (Section 3.2). Specifically, our identification strategy emphasizes three aspects: (1) the relative positioning of the glacier and its proglacial lake, requiring that the lake be located at the glacier terminus along the flow direction; (2) a comprehensive assessment of glacier–lake contact, based on the temporal evolution of the lake and glacier surface geomorphology; and (3) the exclusion of ambiguous cases to ensure the reliability of classification.
5.2. The changing glacier–lake interactions
Glacier–lake contact is a dynamic process, encompassing both the formation of proglacial lakes (Type 2 change glaciers) and the separation of the glacier from their lake (Type 3 change glaciers). In 2020, many Type 2 change glaciers were observed, accounting for nearly 60% of all lake-terminating glaciers. Notably, almost 70% of Type 2 change glaciers had an area smaller than 1 km2, with a median area of 0.51 ± 0.41 km2. The sensitivity of small glaciers to climate change varies considerably. In some regions, while larger glaciers undergo rapid melting and retreating, small glaciers may remain relatively stable (Granshaw and Fountain, Reference Granshaw and Fountain2006; DeBeer and Sharp, Reference DeBeer and Sharp2009). Huss and Fischer (Reference Huss and Fischer2016) reported in the Alps that small glaciers are particularly sensitive to climate change when situated on gentler slopes, covered by debris, or located at lower elevations. As proglacial lakes form on small glaciers, their impact on these glaciers cannot be ignored. This may lead to further differentiation in regional glacier changes.
The separation of proglacial lakes from glaciers is primarily driven by changes in glacier terminus slope and a decline in proglacial lake water level. However, only eight glaciers experienced separation from their proglacial lakes due to lake level decline (Figure S7), six of which were associated with documented GLOFs with clear historical flood records (Lützow and others, Reference Lützow, Veh and Korup2023; Shrestha and others, Reference Shrestha2023). Slope analysis along the normalized flowlines of Type 1 change and Type 3 change glaciers shows that Type 3 change glaciers exhibit steeper terminal slopes than Type 1 change glaciers (Fig. 10). Such steepening indicates a contraction of the lake basin, which restricts further lake expansion towards the glaciers. Consequently, as the glacier retreats, the proglacial lake becomes detached from the glacier terminus.
Figure 10. Glacier slope along the normalized glacier center flow line at the ablation zone. The lines represent the mean values, while the shaded areas indicate the interquartile range. Type 1 change: Glaciers that keep in contact with proglacial lakes (n=146). Type 3 change: Glaciers become detached from proglacial lakes (n=173).
5.3. Influence of glacier–lake contact change on glacier evolution
Although the differences in glacier surface elevation changes among glacier types within each subregion did not reach statistical significance, observational data reveal noteworthy variations. In the Eastern Himalaya, Type 2 change glaciers exhibit a more negative mean glacier surface elevation change rate (−0.85 ± 0.10 m a−1) compared to Type 1 change glaciers (−0.63 ± 0.07 m a−1), with a visually apparent distinction (Fig. 6). Additionally, the terminus retreat distance of Type 2 change glaciers (1.09 ± 0.06 km) is greater than that of Type 1 change glaciers (0.66 ± 0.08 km), with statistical significance (p = 0.003). This suggests a negative impact of proglacial lake expansion on glacier surface elevation change, a phenomenon also reported in the Alps (Carrivick and others, Reference Carrivick2022) and Patagonia (Minowa and others, Reference Minowa, Schaefer and Skvarca2023). In the Central Himalaya, Type 2 change glaciers exhibit statistically significant retreat compared to Type 1 (p < 0.001) and Type 3 (p = 0.003) change glaciers and, yet their surface elevation changes show no substantial divergence. This discrepancy likely reflects a temporal misalignment between the phase of rapid glacier retreat and the observational period for surface elevation change.
Analysis revealed diverse patterns in the terminal velocity change of lake-terminating glaciers. While glacier motion in HMA is generally decelerating due to enhanced surface melting (Neckel and others, Reference Neckel, Loibl and Rankl2017; Dehecq and others, Reference Dehecq2018), our results show that among Type 1 change glaciers (>5 km2), 49 exhibited significant terminal velocity trends, with 11 accelerating and 38 decelerating. It is commonly understood that when glacier termini approach flotation, basal stress decreases sharply, resulting in glacier acceleration (Benn and others, Reference Benn, Hulton and Mottram2007b). Conversely, glaciers respond dynamically to the stress imbalance at the terminus through dynamic thinning, leading to a subsequent slowdown of flow (Nick and others, Reference Nick, Vieli, Howat and Joughin2009; Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021). However, this process appears to represent a relatively short-term and cyclical behavior, as also observed in the periodic flow variations of the Longbasaba Glacier (Liu and others, Reference Liu2020). Such cyclicity typically induces glacier acceleration during phases of dynamic thinning, but deceleration emerges once thinning approaches completion. In addition, terminus deceleration may also be attributed to several factors, including enhanced lateral friction due to glacier narrowing (King and others, Reference King, Dehecq, Quincey and Carrivick2018), ice thickness decreases (Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021), or insufficient lake depth (Minowa and others, Reference Minowa, Schaefer and Skvarca2023). As proglacial lakes detach from glaciers, lake influence diminishes, leading to a pronounced reduction in terminal velocities. Among Type 3 change glaciers (>5 km2), all 18 with significant velocity trends decelerated. Additionally, the migration of peak flow velocities to higher elevation bins (Fig. 7) of all Type 3 change glaciers. Overall, glacier velocity changes reflect complex, stage-dependent interactions between retreat dynamics and lake characteristics. These trends exhibit considerable individual variability, governed by factors such as lake basin geometry, surrounding topography, glacier width, surface slope and ice thickness. And, variations in glacier–lake contact significantly influence glacier dynamics. To improve understanding of proglacial lake impacts on glacier dynamics, high-resolution, spatiotemporal analyses of glacier–lake co-evolution are required. Critical measurements include lake depth, basin morphology, water temperatures near the calving front and glacier thickness, which together will provide essential constraints on the processes driving the behavior of lake-terminating glaciers.
6. Conclusions
We investigated the spatiotemporal distribution and dynamics of lake-terminating glaciers in the Himalaya and Southeastern Tibet from 1990 to 2020. Over the past three decades, the number of lake-terminating glaciers increased by 37.7% (from 419 to 577), whereas their total area was decreased by 161.2 ± 18.1 km2 due to pronounced calving retreat. During this period, 173 glaciers detached from their proglacial lakes, while 331 developed new proglacial lakes. Terminus steepening was the primary driver of glacier–lake separation, with reductions or loss of lake water levels, though rare (n = 8), also contributing. Based on glacier–lake contact changes, glaciers were classified as Type 1 change (persistent contact, n = 246), Type 2 change (new lake formation, n = 331) and Type 3 change (detachment, n = 173). Type 1 change glaciers experienced the largest area loss (73.8 ± 13.1 km2), while Type 2 change glaciers were most numerous. Glacier–lake contact variations significantly modulate glacier dynamics. Among Type 1 change glaciers (>5 km2), 49 exhibited significant velocity changes, with 22.4% (n = 11) accelerating and 77.6% (n = 38) decelerating, whereas all 18 Type 3 change glaciers (>5 km2) with significant trends decelerated. These findings highlight the crucial role of glacier–lake interactions in modulating glacier behavior. However, the lack of in situ observational data limits our understanding of glacier–lake interactions and their variations in this region, particularly regarding lake basin morphology, water depth changes and high temporal resolution monitoring of lake expansion. This is essential for accurately assessing glacier dynamics, developing glacier dynamics and hydro-mathematical models, comprehending glacial chain-related hazards and managing the risks associated with GLOFs.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10088.
Acknowledgements
This work was funded by the National Key R&D Program of China (Grant Nos. 2024YFC3013400 and 2023YFE0102800), National Science Foundation of China (Grant Nos. 42361144874 and 42401101), Science and Technology Plan Projects of Tibet Autonomous Region (XZ202301YD0002C-03) and Institute of Mountain Hazards and Environment (Grant No. IMHE-CXTD-02). We would like to express our thanks to the editors and two anonymous reviewers for their constructive comments and suggestions which ultimately led to an improved manuscript.
Data availability statement
The data sets of the glacier and glacial lake for this study can be found and downloaded from https://doi.org/10.5281/zenodo.15300213
Open access
