Glacier Change , Supraglacial Debris Expansion and Glacial Lake Evolution in the Gyirong River Basin , Central Himalayas , between 1988 and 2015

Himalayan glacier changes in the context of global climate change have attracted worldwide attention due to their profound cryo-hydrological ramifications. However, an integrated understanding of the debris-free and debris-covered glacier evolution and its interaction with glacial lake is still lacking. Using one case study in the Gyirong River Basin located in the central Himalayas, this paper applied archival Landsat imagery and an automated mapping method to understand how glaciers and glacial lakes interactively evolved between 1988 and 2015. Our analyses identified 467 glaciers in 1988, containing 435 debris-free and 32 debris-covered glaciers, with a total area of 614.09 ± 36.69 km2. These glaciers decreased by 16.45% in area from 1988 to 2015, with an accelerated retreat rate after 1994. Debris-free glaciers retreated faster than debris-covered glaciers. As a result of glacial downwasting, supraglacial debris coverage expanded upward by 17.79 km2 (24.44%). Concurrent with glacial retreat, glacial lakes increased in both number (+41) and area (+54.11%). Glacier-connected lakes likely accelerated the glacial retreat via thermal energy transmission and contributed to over 15% of the area loss in their connected glaciers. On the other hand, significant glacial retreats led to disconnections from their proglacial lakes, which appeared to stabilize the lake areas. Continuous expansions in the lakes connected with debris-covered glaciers, therefore, need additional attention due to their potential outbursts. In comparison with precipitation variation, temperature increase was the primary driver of such glacier and glacial lake changes. In addition, debris coverage, size, altitude, and connectivity with glacial lakes also affected the degree of glacial changes and resulted in the spatial heterogeneity of glacial wastage across the Gyirong River Basin.


Introduction
Glaciers are an important component of the global water cycle and are highly sensitive to climate change [1,2].The Tibetan Plateau (TP), widely known as the "Third Pole", hosts the largest coverage of glaciers except for the Antarctic and Arctic regions [3][4][5].Under current global warming, glaciers on the TP are losing their areas and masses through extensive retreats and downwasting [3,6].Of them, the Himalayan glaciers have exhibited the greatest deficit in mass balance during the recent decades and have attracted worldwide attention from both scientific and public communities [1,3,[7][8][9].Changing glaciers in the Himalayas affect not only regional water management [6,10] and chemical compositions of surface fresh water [11][12][13] but also the occurrence of glacier-related disasters, such as glacial lake outburst floods (GLOFs).Historical GLOFs result in catastrophic destructions and fatalities [14][15][16][17] in the Himalayas, which is labeled as one of the major GLOF-vulnerable regions in the world [18][19][20][21][22][23][24][25][26].
Due to difficult access to frigid alpine environments, remote sensing is the most feasible method for monitoring large-scale dynamics in glaciers and glacial lakes [27][28][29][30][31][32].The freely available Landsat imagery since 1972, owing to its multi-decadal record and relatively high spatial resolution (30 m), has been widely used in inventorying and analyzing glacier and glacial lake changes [33][34][35][36][37][38][39].Landsat imagery is also one of the primary data sources for many global and regional glacier products, such as the global land ice measurements from space (GLIMS) glacier data [40], the Gamdam glacier inventory (GGI) [41], the Randolph glacier inventory (RGI) [42], and the second China glacier inventory (CGI) [43].These glacier products improve our understanding of land glacial distributions and have been widely applied in cryo-hydrological modeling.However, these existing products vary by data source, acquisition time, mapping method, and accuracy, which may lead to secondary errors in estimating glacier changes in a specific region.Case studies in the Xainza Xiegang Mountains and Koshi River Basin, for instance, show that the second CGI [43] outperforms other glacier datasets in terms of glacier outlines [5,44].Such discrepancies often require the generation of new glacier datasets with consistent data sources and methods in order to reduce the uncertainties in change analysis.
Debris covers contribute to the spatial heterogeneity of glacier changes across the Himalayas [45][46][47].However, existing glacier products, including GLIMS, GGI, RGI, and CGI, as well as results from other relevant studies [48][49][50], did not distinguish between the debris-free portion (thereafter C-part) and the debris-covered portion (D-part) of a debris-covered glacier (D-glacier).It is difficult to automatically extract the extent of D-part glaciers using optical satellite images alone because of the spectral mixture between supraglacial debris and other moraines [5,49,[51][52][53][54][55].Visual delineation of D-part glaciers is extremely time-consuming.Although a recent study revealed that the debris coverage over glaciers on the south slope of the Mount Everest increased by 17.6% between 1962 and 2011 [56], we know little about the development of glacier debris covers across the Himalayas.Recently, a glacier inventory in the eastern Himalayas was completed using manual digitization on the advanced land observing satellite (ALOS) images acquired between 2006 and 2011 [57].In this regional inventory, D-glaciers were divided into the C-part and D-part.Debris mantle has implications for the glacier response to climate change and development of glacial lakes [46,58], which are also influenced by the slope and flow velocity of ice in the glacier ablation zones [59].Inversely, the dynamics of supraglacial lakes on nine D-glaciers in the Everest region revealed the effect of glacial lakes on promoting ablation via ice cliff calving [60].Inter-conversion of glacial lake types caused by continuous glacier retreat, such as from supraglacial lakes to moraine-dammed lakes, was observed in the Himalayas in previous studies [22,48].The integrated assessments for dynamic evolutions of D-parts and the interplay between glacier and glacial lake changes at multiple spatio-temporal scales are urgently required.
The Gyirong River Basin (GRB), a total area of 4640 km 2 , is located between Gyirong County, China and Bagmati, Nepal (Figure 1) within the central Himalayas.The GRB extends from the lowest elevation of 590 m a.s.l.(above sea level) to the highest 7360 m a.s.l., and contained a total glacier area of 614 km 2 in 1988.The Gyirong Pass in this basin has been the only open trade port between China and Nepal since the Nepal earthquake in 2015, which caused severe destruction and the closure of the previous Nyalam Pass.Historically, GLOFs, such as those originating from Longda Tsho and Zanaco, also destroyed downstream roads and villages in the GRB [23,61].With the booming of Sino-Nepal economic trade, infrastructures such as transboundary highways and railways are complete, under construction, or upcoming.The issues of water resource and glacier-related hazards at the watershed scale have become increasingly important.However, the latest status, changes and, interplay between glaciers and glacial lakes in the GRB during the past decades are not fully understood.
To this end, we aim to (1) investigate the latest distribution of glaciers and glacial lakes using advanced mapping methods in the entire GRB; (2) reveal the changing characteristics for glaciers and glacial lakes using consistent Landsat images at six discrete episodes (1988, 1994, 2001, 2006, 2010, and 2015); and (3) improve the awareness of the interaction between glaciers and glacial lakes in the central Himalayas at a river basin scale.the closure of the previous Nyalam Pass.Historically, GLOFs, such as those originating from Longda Tsho and Zanaco, also destroyed downstream roads and villages in the GRB [23,61].With the booming of Sino-Nepal economic trade, infrastructures such as transboundary highways and railways are complete, under construction, or upcoming.The issues of water resource and glacier-related hazards at the watershed scale have become increasingly important.However, the latest status, changes and, interplay between glaciers and glacial lakes in the GRB during the past decades are not fully understood.
To this end, we aim to (1) investigate the latest distribution of glaciers and glacial lakes using advanced mapping methods in the entire GRB; (2) reveal the changing characteristics for glaciers and glacial lakes using consistent Landsat images at six discrete episodes (1988, 1994, 2001, 2006, 2010, and 2015); and (3) improve the awareness of the interaction between glaciers and glacial lakes in the central Himalayas at a river basin scale.Glacial retreat occurred in all size classes from 1988 to 2015 (Figure 4A).Of the size classes in Figure 4, 1-5 km 2 has the largest area and also suffered the largest absolute loss of area, entirely due to loss by C-glaciers.D-glaciers are considerably larger than C-glaciers, and all glaciers larger than 10 km 2 are D-glaciers (Figure 4B,C).Figure 4B reveals that most glaciers in the GRB are C-glaciers, of which over 99% in number have an area less than 5 km 2 .Although there are only 16 out of a total of 37 D-glaciers with areas greater than 5 km 2 , they contribute more than 86% of the total area in 2015 (Figure 4C).These patterns indicate that the glacial status has close relationship with its size and type.

Data
A total of 21 Landsat images were employed to map glacier and glacial lake extents in the studied GRB (Table 1).These images with a spatial resolution of 30 m were acquired from Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) in six discrete years (episodes) between 1988 and 2015.Because of frequent cloud contamination over the Himalayas, it was practically difficult to obtain ideal-quality imageries during the exact same month [49,62].Based on existing hydrological observations, glacier ablation in the central Himalayas mainly occurs during the summer (between June and August) [30,63], and the seasonal changes in glaciers and glacial lakes from September and December are relatively minor.Therefore, we here acquired the best possible images during the latter period (between September and December).In addition, another 9 images were also selected as auxiliary data in order to reduce the uncertainties caused by sporadic noises such as seasonal snow, terrain shadow, and regional cloud covers.Other datasets used in this study include the second CGI data set [43], high-resolution images from Google Earth, the ASTER Global Digital Elevation Model (ASTER GDEM, 30 m) version 2 [64][65][66], and the annual mean air temperature and precipitation data from the nearest two meteorological stations (Tingri County and Nyalam County) between 1988 and 2015 provided by the China Meteorological Administration (http://www.cma.gov.cn/).

Methods
We mapped the extents of glaciers and glacier lakes using an object-oriented automated algorithm followed by a rigorous visual quality assurance and manual inspection.Specifically, the normalized difference snow/ice index (NDSII) [67,68] was employed for extracting glaciers while the normalized difference water index (NDWI) [69,70] was used for mapping glacial lakes.The detailed procedure of the automated algorithm was described previously for both glacier mapping [5,49] and lake mapping [22,[71][72][73].Both the automatically mapping methods contain pre-processing of Landsat images, object segmentation, setting the threshold and converting from raster data into vector format.Outputs of the automated mapping in the format of shapefile were then meticulously inspected against their source Landsat imagery and, when necessary, modified with assistance of a graphical user interface (GUI) as developed by Wang et al. [71].The GUI-based mapping tool was originally developed for mapping open surface water only [71,72,74].Here we further extended this tool by adding the function of delineating and revising glacier extents using NDSII.As we described previously, mapped glaciers were classified into debris-free glaciers (C-glaciers) and D-glaciers, with the latter further classified into C-part and D-part areas using NDSII with aid of the interactive mapping tool.The C-glacier and C-part areas were well mapped by NDSII and separated individually by drainage divides while debris cover was arduous to be automatically extracted.The D-part area in this study was digitized manually with reference to the CGI data in China and GGI in Nepal, the source Landsat images, and high-resolution Google Earth images.The C-part and D-part layers were saved as single layers and were also merged to generate the new D-glacier data (approximate 1-month processing).
Glacial lakes were classified into five groups based on their spatial relationships with 'mother' glaciers at the same acquisition time, namely (1) supraglacial lake; (2) C-glacier connected lake; (3) D-glacier connected lake; (4) glacier-fed but unconnected lake; and (5) non-glacier-fed lake.Our minimum mapping unit is five Landsat pixels (~0.005 km 2 ) and polygons smaller than this threshold were removed.Quality assurance procedures includes three steps: (1) inspecting each lake and glacier at one time point with reference to Landsat source image and Google Earth high-resolution images and correcting the extent of lake and glacier if an error was found; (2) multi-temporal cross-checking for glacier and glacial lake data over the six episodes between 1988 and 2015 to remove unreasonable changes; and (3) validating the topology for glacier and glacial lake data, including deletion of small silver polygons, repeated polygon removal et al.The margins of D-part glaciers were visually identified and modified with the aid of multi-temporal high-resolution images from Google Earth and Landsat images.The misclassification of glaciers and glacial lakes because of seasonal snow or cloud cover, topographic shadows and cloud shadows was corrected by auxiliary images acquired in the same or adjacent years and cross-check from multi-temporal images.Although the chronological cross-validation is a time-consuming post-editing (approximate 2-month processing), this is essential to generate high quality data with removal of remnant noises and to improve the reliability of change detections.
The uncertainties in mapped glaciers and glacier lake areas were calculated using an error of ±0.5 pixels and the perimeter of each mapped glacier or glacial lake as the method described in previous studies [5,22,59,75].Hypsometric and topographic characteristics were estimated for individual glaciers and glacial lakes in association with DEM, DEM-derived slopes and aspect data.

Glacier Distribution and Changes
Table 2 shows the statistics of our mapped glaciers between 1988 and 2015 across the GRB.During the past three decades, glaciers in the GRB exhibited a declining trend in both number and area: from 467 and 614.09 ± 36.69 km 2 in 1988 to 434 and 513.07 ± 35.11 km 2 in 2015.Glaciers in the GRB were dominated by type C in number, but type D in area, implying that D-glaciers tend to have a greater average area than C-glaciers.Within D-glaciers, the C-part area is predominant in each of the 6 studied episodes.However, the total C-part area shows a decreasing trend while the D-part area increased.Of the total areas of D-glaciers in 1988 (308.58 km 2 ) and in 2015 (282.13 km 2 ), the C-part accounts for 235.80 km 2 (76.42%) and 191.56 km 2 (67.90%), respectively, whereas debris cover increased by 17.79 km 2 , from 23.58% to 32.10% of D-glaciers.
As shown in Table 3, glaciers in the GRB retreated rapidly during 1988-2015, and the retreat rate appears to have accelerated further since 1994.A total of 38 C-glaciers disappeared from 1988 to 2015.Five D-glaciers emerged from the previous C-glaciers between 2001 and 2010, as a result of surface downwasting and then debris formation at the termini of C-glaciers.In general, the areas of C-glaciers, D-glaciers, and C-part glaciers continued to decrease throughout the six episodes from 1988 to 2015, which is in sharp contrast with the consecutive increase in the area of D-part glaciers.We further summarized four major scenarios for glacial changes between 1988 and 2015, which are no change, decrease, conversion, and increase (Figure 2 and Table 4).The conversion among C-glacier, C-part glacier and D-part glacier was first calculated by spatial overlap analysis, then each changing type was summed and categorized into these four scenarios.The most evident decrease was identified in C-glaciers (−76.77km 2 ) and C-part glaciers (−21.88 km 2 ).The conversions among C-part, D-part, and C-glaciers are also obvious.For example, quite a few C-part glaciers were converted to D-part glaciers during the process of glacial downwasting, and some C-part glaciers were disaggregated into smaller sizes owing to glacial splitting at the upper levels.Although glacier changes in the GRB are dominated by the decrease and conversion scenarios, minor expansions are also observed in C-glaciers and C-part glaciers (approximate 0.2 km 2 in total), which were likely caused by ice motion or ice avalanche.In summary, the total area decreased by 74.57km 2 (24.41%) in C-glaciers and by 44.24 km 2 (18.76%) in C-part glaciers between 1988 and 2015, while D-part glaciers increased by 17.79 km 2 (24.44%) (Tables 3 and 4).Glacier changes in the GRB exhibit a strong spatial heterogeneity from 1988 to 2015.In comparison with C-glaciers, individual D-glaciers show greater area losses (Figure 3A).The largest D-glacier decreased by more than 3.0 km 2 but the maximum area loss in C-glaciers is less than 2.0 km 2 .In addition, larger glaciers tend to have lower retreat rates (Figure 3B), whereas some small glaciers have completely perished.Most D-glaciers have a retreat rate lower than 20%, compared with higher than 20% of area loss for most C-glaciers.
Glacial retreat occurred in all size classes from 1988 to 2015 (Figure 4A).Of the size classes in Figure 4, 1-5 km 2 has the largest area and also suffered the largest absolute loss of area, entirely due to loss by C-glaciers.D-glaciers are considerably larger than C-glaciers, and all glaciers larger than 10 km 2 are D-glaciers (Figure 4B,C).Figure 4B reveals that most glaciers in the GRB are C-glaciers, of which over 99% in number have an area less than 5 km 2 .Although there are only 16 out of a total of 37 D-glaciers with areas greater than 5 km 2 , they contribute more than 86% of the total area in 2015 (Figure 4C).These patterns indicate that the glacial status has close relationship with its size and type.Glacial retreat occurred in all size classes from 1988 to 2015 (Figure 4A).Of the size classes in Figure 4, 1-5 km 2 has the largest area and also suffered the largest absolute loss of area, entirely due to loss by C-glaciers.D-glaciers are considerably larger than C-glaciers, and all glaciers larger than 10 km 2 are D-glaciers (Figure 4B,C).Figure 4B reveals that most glaciers in the GRB are C-glaciers, of which over 99% in number have an area less than 5 km 2 .Although there are only 16 out of a total of 37 D-glaciers with areas greater than 5 km 2 , they contribute more than 86% of the total area in 2015 (Figure 4C).These patterns indicate that the glacial status has close relationship with its size and type.

Hypsometry and Topographic Characteristics of Glaciers
Glaciers and glacier changes show distinct altitudinal distributions between 1988 and 2015 (Figure 5).In 2015, the altitudes of all glaciers in GRB range from 3600 m a.s.l. to 7300 m a.s.l., with a mean altitude of 5507 m a.s.l.The majority glacier area is distributed at 5000-6000 m a.s.l., occupying 71.22% of the total glacier area (Figure 5A).C-glaciers generally have a higher altitude than that of D-glaciers (Figure 5B), within which C-parts are distributed higher than D-parts (Figure 5C).In total, glacial area loss between 1988 and 2015 mainly occurred at 4800-6000 m a.s.l.(Figure 5A).Specifically, most area losses of C-glaciers are observed at 4800-5900 m a.s.l. and for D-glaciers at 5200-5800 m a.s.l.(Figure 5B).In contrast, the increasing areas of D-part are more prominent in higher elevations (mainly from 4900 m a.s.l. to 5500 m a.s.l.) on D-glaciers that previously be covered

Hypsometry and Topographic Characteristics of Glaciers
Glaciers and glacier changes show distinct altitudinal distributions between 1988 and 2015 (Figure 5).In 2015, the altitudes of all glaciers in GRB range from 3600 m a.s.l. to 7300 m a.s.l., with a mean altitude of 5507 m a.s.l.The majority glacier area is distributed at 5000-6000 m a.s.l., occupying 71.22% of the total glacier area (Figure 5A).C-glaciers generally have a higher altitude than that of D-glaciers (Figure 5B), within which C-parts are distributed higher than D-parts (Figure 5C).In total, glacial area loss between 1988 and 2015 mainly occurred at 4800-6000 m a.s.l.(Figure 5A).Specifically, most area losses of C-glaciers are observed at 4800-5900 m a.s.l. and for D-glaciers at 5200-5800 m a.s.l.(Figure 5B).In contrast, the increasing areas of D-part are more prominent in higher elevations (mainly from 4900 m a.s.l. to 5500 m a.s.l.) on D-glaciers that previously be covered by C-part glaciers (Figure 5C).This indicates a trend of upward progression of debris covers induced by glacial downwasting.In 2015, the slope of glaciers in the GRB ranges from 10° to 50°, with a mean of 31° (Figure 6).Most C-glaciers have an average slope between 10° and 35°.Compared with D-glaciers (with a mean slope of 27°), C-glaciers (mean slope of 31°) are distributed on steeper terrains.Similarly, C-parts of the D-glaciers are more likely to occur on steeper surfaces, which are consistent with their higher altitudes than D-parts.In 2015, the slope of glaciers in the GRB ranges from 10 • to 50 • , with a mean of 31 • (Figure 6).Most C-glaciers have an average slope between 10 • and 35 • .Compared with D-glaciers (with a mean slope of 27 • ), C-glaciers (mean slope of 31 • ) are distributed on steeper terrains.Similarly, C-parts of the D-glaciers are more likely to occur on steeper surfaces, which are consistent with their higher altitudes than D-parts.In 2015, the slope of glaciers in the GRB ranges from 10° to 50°, with a mean of 31° (Figure 6).Most C-glaciers have an average slope between 10° and 35°.Compared with D-glaciers (with a mean slope of 27°), C-glaciers (mean slope of 31°) are distributed on steeper terrains.Similarly, C-parts of the D-glaciers are more likely to occur on steeper surfaces, which are consistent with their higher altitudes than D-parts.In GRB, most glaciers are small C-glaciers and oriented towards north (Figure 7).The number of north facing glaciers decreased from 90 in 1988 to 78 in 2015 due to the disappearance of some In GRB, most glaciers are small C-glaciers and oriented towards north (Figure 7).The number of north facing glaciers decreased from 90 in 1988 to 78 in 2015 due to the disappearance of some small C-glaciers.In case of area, however, D-glaciers facing west, southwest, and south aspects contribute more than half of the total area since several largest glaciers in GRB are D-glaciers.Between 1988 and 2015, glaciers facing north, northwest, west, southwest, and south decreased more in area than those in the other aspects.

Glacial Lake Inventory and Lake Changes
A total of 107 (4.62 ± 1.12 km 2 ) glacial lakes in 1988 and 148 (7.12 ± 1.50 km 2 ) in 2015 were identified from our mapping (Table 5).The altitude of glacial lakes ranges from 3800 m a.s.l. to 5600 m a.s.l. in 2015 (Figure 8A).Most of the glacial lakes (130) in 2015 are smaller than 0.1 km 2 .The 18

Glacial Lake Inventory and Lake Changes
A total of 107 (4.62 ± 1.12 km 2 ) glacial lakes in 1988 and 148 (7.12 ± 1.50 km 2 ) in 2015 were identified from our mapping (Table 5).The altitude of glacial lakes ranges from 3800 m a.s.l. to 5600 m a.s.l. in 2015 (Figure 8A).Most of the glacial lakes (130) in 2015 are smaller than 0.1 km 2 .The 18 lakes larger than 0.1 km 2 , however, account for over 50% to the total area (Figure 8B).The distribution of lake area in elevation is concentrated lower than that of lake count.

Impact of Debris Cover on Glacier Changes
The retreat rate of D-glaciers in the GRB is lower than that of C-glaciers between 1988 and 2015 because of the role of debris coverage.Most of the D-glacier termini remain nearly stagnant, which was also observed in other Himalayan regions [44,45,58].The stagnant behavior of D-glaciers was likely induced by thermal insulation between the glacier surface and the supraglacial debris [76].The retreat of D-glaciers was further constrained by the thickening of supraglacial debris when they expanded upward [77].Meanwhile, formations of new ice cliffs and supraglacial lakes in the debris-covered zones possibly promoted glacial ablation [45, 59,78].Downwasting of D-part glaciers is a complex process of mass losses for most of Himalayan D-glaciers.

Effect of Glacial Lake Changes on Glaciers
Glacial lakes may accelerate the melting of glaciers owing to the glacio-hydrological interplay.To elucidate this point, we further categorized C-glaciers and D-glaciers by their connections with or without glacial lakes.As shown in Table 6, C-glaciers and D-glaciers connected with glacial lakes experienced higher rates of area loss than those without connections to glacial lakes.The total areas for C-glaciers and D-glaciers decreased by 0.77 km 2 and 0.75 km 2 (summed through the related spatial conversion matrix), which was induced by the expansions of terminally connected glacial lakes between 1988 and 2015.This means that lake expansions contribute 15.88% to the total area loss in the connected C-glaciers and 16.41% to the total area loss of the connected D-glaciers.Clearly, glacial lake expansions result in faster glacial retreat through thermal energy transmission to glaciers from lake water [22,36,79].In general, glacial lakes in GRB expanded rapidly between 1988 and 2015, but the expanding status varies among lake types.The total number of glacial lakes increased by 41 (38.32%) and the total area by 2.50 km 2 (54.11%).Except for C-glacier connected lakes, glacial lakes of the other types grew in number and expanded in area (Table 5).Integrating all lake types, the number and area experienced a net increase from 1988 to 2015 in each of the size classes (Figure 8B).

Impact of Debris Cover on Glacier Changes
The retreat rate of D-glaciers in the GRB is lower than that of C-glaciers between 1988 and 2015 because of the role of debris coverage.Most of the D-glacier termini remain nearly stagnant, which was also observed in other Himalayan regions [44,45,58].The stagnant behavior of D-glaciers was likely induced by thermal insulation between the glacier surface and the supraglacial debris [76].The retreat of D-glaciers was further constrained by the thickening of supraglacial debris when they expanded upward [77].Meanwhile, formations of new ice cliffs and supraglacial lakes in the debris-covered zones possibly promoted glacial ablation [45,59,78].Downwasting of D-part glaciers is a complex process of mass losses for most of Himalayan D-glaciers.

Effect of Glacial Lake Changes on Glaciers
Glacial lakes may accelerate the melting of glaciers owing to the glacio-hydrological interplay.To elucidate this point, we further categorized C-glaciers and D-glaciers by their connections with or without glacial lakes.As shown in Table 6, C-glaciers and D-glaciers connected with glacial lakes experienced higher rates of area loss than those without connections to glacial lakes.The total areas for C-glaciers and D-glaciers decreased by 0.77 km 2 and 0.75 km 2 (summed through the related spatial conversion matrix), which was induced by the expansions of terminally connected glacial lakes between 1988 and 2015.This means that lake expansions contribute 15.88% to the total area loss in the connected C-glaciers and 16.41% to the total area loss of the connected D-glaciers.Clearly, glacial lake expansions result in faster glacial retreat through thermal energy transmission to glaciers from lake water [22,36,79].

Impact of Glacier Changes on Glacial Lakes
In the context of glacial retreat, the statuses of glacier changes significantly influence the formation, conversion, and current standing for glacial lakes.For supraglacial lakes on D-glaciers, they vary in location, size, and shape, as reported in previous publications [22,23,36].A good case is the rapid evolution and expansion of supraglacial lakes on Lalaga Glacier (Figure 9A), where the lake labeled by green dot reached an area of 0.31 km 2 in 2017 from 0.01 km 2 in 1988.A total of 9 C-glacier connected lakes became glacier-fed unconnected lakes between 1988 and 2015 (such as the lake in Figure 9B), caused by continuous glacier retreat.D-glacier connected lakes (Figure 9C), generally named as moraine-dammed lakes, expanded by 1.77 times in area from 1988 to 2015, and are prone to form large lakes that are potentially dangerous to the downstream regions.Historically, this type of lake usually resulted in GLOFs.For example, Longda Tsho and Zanaco in the GRB burst and caused flooding hazards on 25 August 1964 and 7 June 1995, respectively [23,61].This implies that close attention should be paid to the expansion of D-glacier connected lakes in the future.
Generally, glacial lakes unconnected from glaciers are relatively stable, such as the lakes in Figure 9D,E-H .C-glacier connected lakes are the most dynamic among all lake types.The process from connection to disconnection between a newly formed lake and its feeding C-glacier is distinctly observed across the GRB (Figure 9E-H ).Such conversion from new formation to disconnection makes a major contribution to the increase in number and area for glacier-fed unconnected lakes (Table 5).Similarly, the increase in number of non-glacier-fed lakes was caused by the disappearance of upstream glaciers (Figure 9H).It is obvious that glacier changes affect the hydrological connectivity between lakes and their feeding glaciers and therefore exert a profound influence on the evolution and status of glacial lakes.As discussed above, warming-induced glacial retreat controls the evolution of glacial lakes in the study area.Retreating glaciers cause a large volume of meltwater and space for connected lakes to expand [15,86].The future glacier melt runoff was projected to increase by a glacio-hydrological model in the Langtang watershed, located at the south of our GRB [10].As glacial lakes usually continue to expand in the case of accelerated glacier retreat and wastage, the risk of GLOFs is likely to increase in the Himalayas.We agree with others who predict GLOF frequency will increase

Climate as a Major Factor Driving Glacial Retreat and Lacustrine Evolution
Glacial retreat is known to be a result of rising temperature in the Himalayas [3,7,49].Warming is significantly observed in the central Himalayas (Figure 10A) based on in situ meteorological data from the nearest weather stations in Nyalam County and Tingri County, where annual mean air temperature increased by 0.35 • C per decade (adjusted R 2 value of 0.29 and p-value of 0.002) during 1988-2015.Meanwhile, annual precipitation has remained stable (Figure 10B) at a negligible increasing rate of 2.83 mm per decade (adjusted R 2 value of −0.04 and p-value of 0.875).The gridded GPCP precipitation data, however, show a decreasing trend across the high altitude region of the Himalayas from 1979 to 2014 [22].This confirms that continuous warming is likely the dominant contributor to the accelerated glacial retreat in the GRB between 1988 and 2015 rather than the variation in precipitation.Debris cover is also considered as an important factor of glacier changes [45].Upward expansion of D-part glaciers reflects the response of glaciers to climate change in the central Himalayas.Other properties, such as size, steepness, aspect, altitude, and connectivity with glacial lakes, also influence the heterogeneous wastage of glaciers [80][81][82].The changing composition of glacier surface, such as the content of black sooty particles [83][84][85], is also a contributor to glacial retreat.All these factors indicate the complex nature of alpine glacial changes.during the next decades and into the 22nd century [87].Therefore, a continuous monitoring of glaciers and glacial lakes will be crucially important for GLOF risk assessments, hazard management and mitigation.

Conclusions
This study provides a comprehensive monitoring of glaciers and glacial lakes across the entire GRB using archival Landsat imagery acquired in the past nearly three decades.Our results reveal dramatic retreats of C-glaciers and C-part glaciers in association with significant upward expansions of D-part glaciers in the studied six episodes between 1988 and 2015.We attribute the expansion of supraglacial debris primarily to the internal conversion from C-part to D-part.While C-glaciers retreated rapidly because of the lack of debris cover, D-glaciers appeared more stable due to thermal insulation under the thickening supraglacial debris.
We also demonstrate that the dynamics and conversions among glacial lake types were controlled strongly by glacier changes via hydrological connectivity.On the one hand, expanding glacial lakes likely accelerated the retreat of connected mother glaciers through thermal transmits.On the other hand, glacier lakes were hardly influenced by the upstream glaciers once they were disconnected from the lakes.The interplay between glaciers and glacial lakes implies the complexity of glacial downwasting.Under the current climate change, glaciers will possibly continue to retreat and downwaste.As a result, D-glacier-connected lakes may expand in the near future.Therefore, hazard assessment for the potentially dangerous glacial lakes is highly recommended in this important river basin, which contains the Gyirong Pass, the only trading port that currently functions between China and Nepal.
Author Contributions: Y.N. designed the research and collected the data; S.J. processed the remote sensing data and created the glacier and glacial lake inventory; J.W. developed the mapping tool; Y.N. and J.W. analyzed the data; Y.N., Q.L., and J.W. discussed and commented the paper; all authors contributed to the writing of the manuscript and improving the review.As discussed above, warming-induced glacial retreat controls the evolution of glacial lakes in the study area.Retreating glaciers cause a large volume of meltwater and space for connected lakes to expand [15,86].The future glacier melt runoff was projected to increase by a glacio-hydrological model in the Langtang watershed, located at the south of our GRB [10].As glacial lakes usually continue to expand in the case of accelerated glacier retreat and wastage, the risk of GLOFs is likely to increase in the Himalayas.We agree with others who predict GLOF frequency will increase during the next decades and into the 22nd century [87].Therefore, a continuous monitoring of glaciers and glacial lakes will be crucially important for GLOF risk assessments, hazard management and mitigation.

Conclusions
This study provides a comprehensive monitoring of glaciers and glacial lakes across the entire GRB using archival Landsat imagery acquired in the past nearly three decades.Our results reveal dramatic retreats of C-glaciers and C-part glaciers in association with significant upward expansions of D-part glaciers in the studied six episodes between 1988 and 2015.We attribute the expansion of supraglacial debris primarily to the internal conversion from C-part to D-part.While C-glaciers retreated rapidly because of the lack of debris cover, D-glaciers appeared more stable due to thermal insulation under the thickening supraglacial debris.
We also demonstrate that the dynamics and conversions among glacial lake types were controlled strongly by glacier changes via hydrological connectivity.On the one hand, expanding glacial lakes likely accelerated the retreat of connected mother glaciers through thermal transmits.On the other hand, glacier lakes were hardly influenced by the upstream glaciers once they were disconnected

Figure 1 .
Figure 1.Study area and distribution of glaciers and glacial lakes in 2015.The inset map is from ESRI's world basemap and the background image was acquired on 25 October 2015 from Landsat 8 OLI.Rectangles 1 and 2 are the areas covered in Figures 2 and 3.
Figure 1.Study area and distribution of glaciers and glacial lakes in 2015.The inset map is from ESRI's world basemap and the background image was acquired on 25 October 2015 from Landsat 8 OLI.Rectangles 1 and 2 are the areas covered in Figures 2 and 3.

Figure 1 .
Figure 1.Study area and distribution of glaciers and glacial lakes in 2015.The inset map is from ESRI's world basemap and the background image was acquired on 25 October 2015 from Landsat 8 OLI.Rectangles 1 and 2 are the areas covered in Figures 2 and 3.

Figure 2 .
Figure 2. Spatial patterns of glacier changes between 1988 and 2015.

Figure 3 .
Figure 3. Glacier area loss (A) and change rate (B) between 1988 and 2015.The glacier extents shown in both panels were mapped from 1988 imagery.

Figure 3 .
Figure 3. Glacier area loss (A) and change rate (B) between 1988 and 2015.The glacier extents shown in both panels were mapped from 1988 imagery.

Figure 3 .
Figure 3. Glacier area loss (A) and change rate (B) between 1988 and 2015.The glacier extents shown in both panels were mapped from 1988 imagery.

Figure 4 .
Figure 4. Area and frequency of glaciers in various size classes: (A) glacier total areas in the studied six years; (B) C-glaciers in 1988 and 2015; (C) D-glaciers in 1988 and 2015.Glacier numbers are given in green.

Figure 4 .
Figure 4. Area and frequency of glaciers in various size classes: (A) glacier total areas in the studied six years; (B) C-glaciers in 1988 and 2015; (C) D-glaciers in 1988 and 2015.Glacier numbers are given in green.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 19by C-part glaciers (Figure5C).This indicates a trend of upward progression of debris covers induced by glacial downwasting.

Figure 5 .
Figure 5. Altitudinal distribution of glaciers and their changes between 1988 and 2015: (A) distributions of glacier areas and the net area loss; (B) area distributions of C-glaciers and D-glaciers; (C) area distributions and changes of C-part and D-part glaciers.The shaded area in (C) is only the area increase in D-parts.

Figure 5 .
Figure 5. Altitudinal distribution of glaciers and their changes between 1988 and 2015: (A) distributions of glacier areas and the net area loss; (B) area distributions of C-glaciers and D-glaciers; (C) area distributions and changes of C-part and D-part glaciers.The shaded area in (C) is only the area increase in D-parts.

Figure 5 .
Figure 5. Altitudinal distribution of glaciers and their changes between 1988 and 2015: (A) distributions of glacier areas and the net area loss; (B) area distributions of C-glaciers and D-glaciers; (C) area distributions and changes of C-part and D-part glaciers.The shaded area in (C) is only the area increase in D-parts.

Figure 6 .
Figure 6.Area distribution of glaciers at various slope classes in 2015: (A) area distributions of C-glaciers and D-glaciers; (B) area distributions of C-part and D-part of D-glaciers.The average slope was calculated based on DEM with a 30 m spatial resolution for individual C-glacier, D-glacier, C-part and D-part glaciers, and then we summed the area for each class.

Figure 6 .
Figure 6.Area distribution of glaciers at various slope classes in 2015: (A) area distributions of C-glaciers and D-glaciers; (B) area distributions of C-part and D-part of D-glaciers.The average slope was calculated based on DEM with a 30 m spatial resolution for individual C-glacier, D-glacier, C-part and D-part glaciers, and then we summed the area for each class.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 19 small C-glaciers.In case of area, however, D-glaciers facing west, southwest, and south aspects contribute more than half of the total area since several largest glaciers in GRB are D-glaciers.Between 1988 and 2015, glaciers facing north, northwest, west, southwest, and south decreased more in area than those in the other aspects.

Figure 7 .
Figure 7. Aspect distributions of glaciers in GRB: (A) areas and counts for all glaciers; (B) areas for C-glaciers and D-glaciers in 1988 and 2015.

Figure 7 .
Figure 7. Aspect distributions of glaciers in GRB: (A) areas and counts for all glaciers; (B) areas for C-glaciers and D-glaciers in 1988 and 2015.

Figure 8 .
Figure 8. Characteristics of altitudinal distributions and size groups of glacial lakes: (A) altitudinal distributions in area and frequency for glacial lakes in 2015; (B) area and frequency distributions of glacial lakes among size classes in 1988 and 2015.

Figure 8 .
Figure 8. Characteristics of altitudinal distributions and size groups of glacial lakes: (A) altitudinal distributions in area and frequency for glacial lakes in 2015; (B) area and frequency distributions of glacial lakes among size classes in 1988 and 2015.

Figure 10 .
Figure 10.Air temperature (A) and annual precipitation (B) changes in the central Himalayas observed from the nearest meteorological stations between 1988 and 2015.The blue line represents the linear regression.The red curve is the 5-year moving average.

Funding:
This research was funded in part by the National Natural Science Foundation of China [Grant No. 41571104], CAS STS Program [KFJ-STS-ZDTP-015], CAS "Light of West China" Program, CPC Sichuan Provincial Committee Organization Department Talents Program, and CAS IMHE 135 Program [SDS-135-1708].

Figure 10 .
Figure 10.Air temperature (A) and annual precipitation (B) changes in the central Himalayas observed from the nearest meteorological stations between 1988 and 2015.The blue line represents the linear regression.The red curve is the 5-year moving average.

Table 4 .
Summary of various scenarios of glacier changes between 1988 and 2015.

Figure 2 .
Spatial patterns of glacier changes between 1988 and 2015.

Table 1 .
Landsat images applied for glacier and glacial lake mapping.

Table 2 .
Statistics of glacier inventories in six episodes between 1988 and 2015.

Table 4 .
Summary of various scenarios of glacier changes between 1988 and 2015.

Table 5 .
Inventory of glacial lakes in various types at six time points between 1988 and 2015 (units: count (area, km 2 )).

Table 6 .
Comparison of glacier changes by their connection with or without glacial lakes.

Table 6 .
Comparison of glacier changes by their connection with or without glacial lakes.Note: two proglacial lakes in 1988 were at the same terminus of one D-glacier. a