Evaluating the Scale and Potential of GLOF in the Bhutan Himalayas Using a Satellite-Based Integral Glacier – Glacial Lake Inventory

A comprehensive glacier–glacial lake inventory was developed for the Bhutan Himalayas based on satellite observations between 1987–1990 and 2006–2011. In total, 733 lakes (covering 82.6 km2) were delineated between 4000 and 6000 m a.s.l. and their relationships to associated glaciers were documented. Using this new inventory, the scale and potential for glacial lake outburst flooding (GLOF) based on multiple criteria was examined. This included a history of connectivity characteristics of glacial lakes to mother glaciers, potential flood volumes, and debris-cover of mother glaciers in addition to the conventional criteria of expansion rate and lake size. The majority of the lakes with high expansion rates (more than double in size) and large areas (>0.1 km2) met the conditions of being continuously in contact with a mother debris-covered glacier for nearly 20 years. Based on these multiple criteria, two lakes were identified as having potential for large-scale GLOF. Potentially dangerous glacial lakes listed in the International Centre for Integrated Mountain Development (ICIMOD) study were re-visited, and some overlaps with the glacier–glacial lake inventory were found.


Introduction
Nearly 10,000 glacial lakes are located in the Hindu Kush Himalayan (HKH) region and on the Tibetan Plateau [1,2].Previous studies have reported changes in the surface areas of these lakes in recent decades.Compared to lakes in the western HKH region, glacial lakes in the eastern HKH region have exhibited larger increases in total and individual areas; notably, rapidly expanding supraand pro-glacial lakes in contact with glaciers [3].A previous study compiled a glacial lake inventory containing 30 lakes in the Bhutan Himalayas in the easternmost HKH region, and recorded a history of growth and integration of small lakes into fewer but larger lakes since the 1950s [4].Lakes on the northern side of the Himalayan Range have shown slow expansion since the 1950s.In contrast, those on the southern side formed after the 1950s have rapidly expanded [5].
Glacial lake outburst flooding (GLOF) occurs when glacial lakes collapse, drain large volumes of water, and cause severe downstream damage (e.g., [6][7][8][9]).To understand the recent status of glacial lakes and their GLOF potential in the Bhutan Himalayas, the International Centre for Integrated Mountain Development (ICIMOD) compiled a glacial lake inventory by means of topographic maps and satellite images.They identified 24 potentially dangerous glacial lakes (PDGLs) in the Bhutan Himalayas [10], which were selected using the criteria of: (1) expansion since 1968; (2) threshold lake area (>0.1 km 2 ) and distance to glacier (<0.5 km); (3) water depth from in situ survey; and (4) various moraine-dam conditions and characteristics, such as materials, shapes, drainage, topography, presence/absence of internal ice, slope stability, signs of past breaching and GLOF, and seepage flow.
Since the ICIMOD study, there have been rapid improvements in observations and continuous progress in the understanding of GLOF processes [11].For example, potential GLOF was assessed in the Cordillera Blanca, Peru, and some lakes were identified as having higher GLOF potential [12].In general, there are both external/short-term and internal/long-term factors responsible for the collapse of moraine-dammed lakes.Short-term factors can be the presence/absence of slope movements into lakes, earthquakes, high-waves associated with flooding in upstream lakes, blocking of underground outflow channels, and intensive rainfall and/or snowmelt.Long-term factors can be the destruction of a dam by melting of buried ice, hydrostatic pressure effects, and gradual degradation of dams [13].Nonetheless, it is extremely difficult to collect accurate information about many of these factors because of a region's remoteness.Therefore, prior to any assessment of GLOF, it is necessary to build a geospatial database to gather basic statistical information and to evaluate possible GLOF scales and potential by means of satellite remote-sensing techniques.Following this strategy, an up-to-date comprehensive glacial lake inventory of the Bhutan Himalayas was generated from high-resolution satellite images acquired by the Advanced Land Observing Satellite (ALOS) operated by the Japan Aerospace Exploration Agency (JAXA).Another glacier data set was generated to describe the connectivity of the glacial lakes identified in the glacial lake inventory with respect to their associated glaciers, and the concept of potential flood volume (PFV) was introduced, which enabled the estimation of drainable water volumes from damming topography [14].The current study re-assesses GLOF scale and potential in the Bhutan Himalayas using multiple criteria, including connectedness of glacial lakes to associated glaciers and PFVs.

Study Site
The current study focused on glacial lakes in the Bhutan Himalayas (27.496 • -28.166 • N; 88.791 • -92.161 • E) mostly located in the Kingdom of Bhutan (Druk Yul is the local name) (Figure 1).Bhutan is a landlocked country in southern Asia surrounded by India and China with a land area of approximately 40,000 km 2 and a population of approximately 700,000.Glaciers and glacial lakes exist along the highly glacierised area of the main Himalayan Range and its branching ridges (within 4000 to 7500 m a.s.l.).This west-to-east high mountain range separates climatic conditions from north to south.The northern side has an arid climate with an annual precipitation of 500 to 1000 mm, while the southern side is influenced by the Indian Monsoon, especially in its high season between June and September, with a high level of precipitation (1000 to 2500 mm) [10].On the southern side, a river system, typically divided into several sub-basins, flows north to south.Tectonically formed gentle hillslopes are dominant on the northern side (i.e., a southern rim of the Tibetan Plateau) (e.g., [15,16]).All of the river basins eventually connect to the Brahmaputra River system.Glacier surface areas compose approximately 642 km 2 of the Bhutanese domain, and have exhibited rapid shrinkage (>20%) between 1990 and 2010 [18].The glaciers in this region are summeraccumulation-type, whereby accumulation and ablation occur simultaneously during the monsoon season [19], and thus are highly sensitive to changing climate.In situ measurements by [20] revealed that glacier mass loss in the Bhutan Himalayas (−1.12 to −2.04 m w.e. a −1 ) is greater than in neighboring regions (−0.22 to −0.51 m w.e. a −1 in Khumbu, Nepal, for example).

Subsection Completion of Glacier and Glacial Lake Inventory
A former glacial lake inventory was generated by the ICIMOD covering 2674 lakes in the Bhutan Himalayas [10] based on topographic maps published in the 1950s-1970s by the Survey of India [21].However, this inventory had a problem of significant horizontal offset [21].The non-systematic horizontal shift and the 20-year temporal range were not ideal for recording areal evolution, present location, and other geographic parameters.One of the main outcomes from the Bhutan/Japan joint research project carried out between 2009 and 2012 was the production of an up-to-date glacial lake inventory generated from high-resolution satellite images acquired by ALOS as operated by JAXA [22].Glacier surface areas compose approximately 642 km 2 of the Bhutanese domain, and have exhibited rapid shrinkage (>20%) between 1990 and 2010 [18].The glaciers in this region are summer-accumulation-type, whereby accumulation and ablation occur simultaneously during the monsoon season [19], and thus are highly sensitive to changing climate.In situ measurements by [20] revealed that glacier mass loss in the Bhutan Himalayas (−1.12 to −2.04 m w.e. a −1 ) is greater than in neighboring regions (−0.22 to −0.51 m w.e. a −1 in Khumbu, Nepal, for example).

Subsection Completion of Glacier and Glacial Lake Inventory
A former glacial lake inventory was generated by the ICIMOD covering 2674 lakes in the Bhutan Himalayas [10] based on topographic maps published in the 1950s-1970s by the Survey of India [21].However, this inventory had a problem of significant horizontal offset [21].The non-systematic horizontal shift and the 20-year temporal range were not ideal for recording areal evolution, present location, and other geographic parameters.One of the main outcomes from the Bhutan/Japan joint research project carried out between 2009 and 2012 was the production of an up-to-date glacial lake inventory generated from high-resolution satellite images acquired by ALOS as operated by JAXA [22].
The definition of a glacial lake and the delineation method for this study are explained in detail in Ukita et al. [21], which compiled the lakes in four basins in our study site.The visual identification and manual delineation of lake outlines were performed on ortho-rectified satellite images that had been pan-sharpened from ALOS PRISM and AVNIR-2 images acquired between 2006 and 2011 (Figure 2a,b).A PRISM-derived digital surface model (DSM) was used to distinguish flat lake surfaces from the surrounding topography.The specific dates of the ALOS images used are listed in the attached attribute table.Because our previous study was strongly motivated by the disaster management of GLOF, small lakes of <0.01 km 2 were excluded from the analysis [21].Moraine-dammed and bedrock dammed lakes are included, whereas those situated farther downstream from the Little Ice Age moraines (>2 km) were also excluded as they are less related to glaciation.On the other hand, lakes both in contact with glaciers and not in contact were included.
The definition of a glacial lake and the delineation method for this study are explained in detail in Ukita et al. [21], which compiled the lakes in four basins in our study site.The visual identification and manual delineation of lake outlines were performed on ortho-rectified satellite images that had been pan-sharpened from ALOS PRISM and AVNIR-2 images acquired between 2006 and 2011 (Figure 2a,b).A PRISM-derived digital surface model (DSM) was used to distinguish flat lake surfaces from the surrounding topography.The specific dates of the ALOS images used are listed in the attached attribute table.Because our previous study was strongly motivated by the disaster management of GLOF, small lakes of <0.01 km 2 were excluded from the analysis [21].Morainedammed and bedrock dammed lakes are included, whereas those situated farther downstream from the Little Ice Age moraines (>2 km) were also excluded as they are less related to glaciation.On the other hand, lakes both in contact with glaciers and not in contact were included.The details of ALOS-derived image processing, DSMs, and in situ validation of delineated lake outlines are described in [21,23].They reported that the PRISM DSM used for the ALOS image ortho- The details of ALOS-derived image processing, DSMs, and in situ validation of delineated lake outlines are described in [21,23].They reported that the PRISM DSM used for the ALOS image ortho-rectification had a vertical accuracy of 2.3 ± 7.8 m against ground control points (N = 3268) and a delineated lake outline had a mean distance of 9.5 m (11.7 m in root mean square) from a global positioning system (GPS)-measured outline.It was noted, for the purposes of the current study, that these error values were sufficiently small for delineating lakes >0.01 km 2 .
The glacier inventory, used in the current study, was also developed from the same ALOS images from the period of 2006 to 2011 [16] (Figure 2a,b).The outlines of the glaciers and debris-covered areas were manually delineated.In total, 1364 debris-free glaciers and 219 debris-covered glaciers were identified.There were some difficulties in delineation associated with steep adjacent slopes.However, judging from a comparison and validation with other glacier inventories, these difficulties were not considered to pose serious problems for the current study.

Multiple Criteria for Scale and Potential
In the Bhutan Himalayas, pro-glacial lakes in contact with glaciers are expanding at faster rates than supra-glacial lakes and those without contact are expanding more slowly [4].Many of the new lakes, formed after 1990, are supra-glacial and located on debris-covered glaciers [24].This implies that the connection type and the presence of supra-glacial debris are strongly related to continuous lake expansion.In previous studies [2][3][4]25], lakes still in contact with associated glaciers were classified as supra-glacial (Cs) and pro-glacial (Cp) lakes.In the current study, four categories were used, including the Cs and Cp types.Because no satellite remote-sensing technique currently provides information on sub-lake ice bodies, lakes surrounded by a glacier surface were categorized as Cs types and those with a boundary to non-glaciated terrain as Cp types.Lakes with distant glaciers were classified as glacier-fed (Df) type and those without upstream glaciers as non-glacier-fed (Dn) type.Glacier outlines and information were obtained from the ALOS-derived glacier inventory described above (see [16]).For example, Figure 2a,b presents spatial distributions of glaciers and associated debris-covered areas.In addition, moraine-damming or bedrock-damming are distinguished by geomorphological visual interpretation in Google Earth (Google, Menlo Park, CA, USA).
A history of lake expansion is essential for GLOF assessment.To obtain past lake outlines, the ALOS-derived lake outlines were back-traced to ortho-rectified panchromatic (10-m spatial resolution) images acquired between 1987 and 1990 by Satellite Pour l'Observation de la Terre-1 (SPOT-1) (Table 1, Figure 2c,d).To describe the expansion history over an interval of around 20 years, an expansion factor was calculated as a ratio of the ALOS-derived lake area (2006-2010) to the SPOT-1-derived lake area (1987)(1988)(1989)(1990), if a lake had been formed in the SPOT-1 image.The connection types during the early period were also identified in reference to the SPOT-1 imagery (without glacier outline delineation) using the same definition as the later ALOS period.PFVs were also calculated according to the definition given in [14].It was assumed that lake water would be discharged until the lake surface dropped to a critical level.This critical level was linked to an elevation gradient defined with respect to the angle from downstream of the moraine (or the damming bedrock) to the shoreline of the lake after hypothetical drainage.Note that this is not the same as the slope gradient of the front side of the moraine (or the damming bedrock), which is significantly steeper.This method was validated from GLOF of four lakes in the Himalayas.Based on these observations, 10 • was adopted as the critical value of the elevation gradient for outburst termination [14].

Current Area and Other Basic Statistics
The current study's new Bhutan Himalayas glacial lake inventory contains 733 lakes with a total lake area of 82.6 km 2 within the range of 4000 to 6000 m a.s.l.(Figure 1a).Individual lake areas range from 0.01 to 5.8 km 2 with a mean of 0.1 km 2 (Figure 3).The smallest lake-size category (0.01 to 0.05 km 2 ) accounted for 58% (n = 427) of the total lakes; smaller lakes are more abundant than larger lakes.The mode was found in the size range of 0.1 to 0.5 km 2 , which occupied 31.7 km 2 or 38% of the total lake area (Figure 3b).Hydrological boundaries were defined by nine river basins as was done by [21] (Figure 1a).The Kuri Chu Basin had the highest number of lakes (n = 146), followed by the Pho Chu Basin (n = 142), which has the highest density of lake number divided by basin area (6.1 × 10 −2 lakes per km 2 ) (Table 2).Fewer lakes were distributed in the Amo Chu (n = 21) and Wang Chu Basins (n = 13).In terms of lake area, the Kuri Chu Basin contained the greatest total lake area (18.3 km 2 , 22.1%) followed by the Tibetan Basin (17.1 km 2 ).The Wang Chu and Amo Chu Basins had smaller lake-to-basin densities of 2.6 × 10 −3 and 3.4 × 10 −3 lakes per km 2 as well as smaller lake-coverages of 0.5 km 2 and 1.4 km 2 , respectively.
The vertical distribution of the lakes ranged between 4134 and 5885 m a.s.l. with a mean elevation of 5008 m a.s.l.(Table 2, Figure 4).The elevation category between 5000 and 5100 m a.s.l. had the largest number of lakes (n = 146, 20%).The elevation distribution of glacial lakes for each river basin showed that the largest lake was in the Tibetan Basin at 5129 m a.s.l.The mean elevation was highest in the Tibetan Basin (5137 m a.s.l.) and lowest in the Wang Chu Basin (4571 m a.s.l.) (Figure 4j).Glacial lakes in the Kuri Chu Basin had the largest vertical range of 1700 m followed by those in the Pho Chu Basin of 1233 m.Lakes in the Mangde Chu Basin had the smallest range of 708 m.PFVs were also calculated according to the definition given in [14].It was assumed that lake water would be discharged until the lake surface dropped to a critical level.This critical level was linked to an elevation gradient defined with respect to the angle from downstream of the moraine (or the damming bedrock) to the shoreline of the lake after hypothetical drainage.Note that this is not the same as the slope gradient of the front side of the moraine (or the damming bedrock), which is significantly steeper.This method was validated from GLOF of four lakes in the Himalayas.Based on these observations, 10° was adopted as the critical value of the elevation gradient for outburst termination [14].

Current Area and Other Basic Statistics
The current study's new Bhutan Himalayas glacial lake inventory contains 733 lakes with a total lake area of 82.6 km 2 within the range of 4000 to 6000 m a.s.l.(Figure 1a).Individual lake areas range from 0.01 to 5.8 km 2 with a mean of 0.1 km 2 (Figure 3).The smallest lake-size category (0.01 to 0.05 km 2 ) accounted for 58% (n = 427) of the total lakes; smaller lakes are more abundant than larger lakes.The mode was found in the size range of 0.1 to 0.5 km 2 , which occupied 31.7 km 2 or 38% of the total lake area (Figure 3b).Hydrological boundaries were defined by nine river basins as was done by [21] (Figure 1a).The Kuri Chu Basin had the highest number of lakes (n = 146), followed by the Pho Chu Basin (n = 142), which has the highest density of lake number divided by basin area (6.1 × 10 −2 lakes per km 2 ) (Table 2).Fewer lakes were distributed in the Amo Chu (n = 21) and Wang Chu Basins (n = 13).In terms of lake area, the Kuri Chu Basin contained the greatest total lake area (18.3 km 2 , 22.1%) followed by the Tibetan Basin (17.1 km 2 ).The Wang Chu and Amo Chu Basins had smaller lake-to-basin densities of 2.6 × 10 −3 and 3.4 × 10 −3 lakes per km 2 as well as smaller lake-coverages of 0.5 km 2 and 1.4 km 2 , respectively.The vertical distribution of the lakes ranged between 4134 and 5885 m a.s.l. with a mean elevation of 5008 m a.s.l.(Table 2, Figure 4).The elevation category between 5000 and 5100 m a.s.l. had the largest number of lakes (n = 146, 20%).The elevation distribution of glacial lakes for each river basin showed that the largest lake was in the Tibetan Basin at 5129 m a.s.l.The mean elevation was highest in the Tibetan Basin (5137 m a.s.l.) and lowest in the Wang Chu Basin (4571 m a.s.l.) (Figure 4j).Glacial lakes in the Kuri Chu Basin had the largest vertical range of 1700 m followed by those in the Pho Chu Basin of 1233 m.Lakes in the Mangde Chu Basin had the smallest range of 708 m.

Glacier-Glacial Lake Connections and Recent Expansion
As described above, the integrated glacier-glacial lake inventory provided information on the connectivity of lakes and glaciers into four types (Cs, Cp, Df, and Dn).By comparing the images from the early SPOT-1 period and the later ALOS period, there were 16 different combinations (i.e., 4 × 4 from the set of Cs, Cp, Df, and Dn) to describe changes in the lake-glacier relationship.Table 3 summarizes these combinations, where 662 lake-glacier relationships were identified.The others were difficult to confirm due to exclusion from SPOT-1 coverage or cloud cover, or were thought to be formed prior to the lake formation stage.The lake relationships comprise: 516 Df lakes, of which 403 lakes were already in the Df type before the 1990s; 100 Cp lakes, of which 67 lakes were in the Cp type before the 1990s; 58 Dn lakes, of which 43 lakes were in the Dn type; and 14 Cs lakes, of which 10 lakes were not yet formed by the 1990s.One lake was changed from Cs to Cp, 69 lakes showed a history of changing from Cp to Df, and nine lakes changed from Df to Dn.There were no lakes showing a history of moving from Dn to Df or Dn, to Df, to Cp, or to Cs. Figure 5 plots the relationship between lake area ~1990 and expansion factors over the last 20 years with colors indicating changes in glacier-glacial lake connectivity.The expansion factors for the Cs lakes were the highest at 11.1 ± 12.3 times (n = 4) followed by those of the Cp-from-Cp lakes at 2.8 ± 3.5 times (n = 67) (Figure 5a).The Df-from-Df lakes had smaller mean expansion factors and a smaller standard deviation of 1.3 ± 0.7 times compared to those of the Df-from-Cp lakes (1.9 ± 5.2 times) (Figure 5b,c).The Dn-from-Df lakes and the Dn lakes showed only minor variation in their expansion factors concentrated around 1.0 (Figure 5d).

Glacier-Glacial Lake Connections and Recent Expansion
As described above, the integrated glacier-glacial lake inventory provided information on the connectivity of lakes and glaciers into four types (Cs, Cp, Df, and Dn).By comparing the images from the early SPOT-1 period and the later ALOS period, there were 16 different combinations (i.e., 4 × 4 from the set of Cs, Cp, Df, and Dn) to describe changes in the lake-glacier relationship.Table 3 summarizes these combinations, where 662 lake-glacier relationships were identified.The others were difficult to confirm due to exclusion from SPOT-1 coverage or cloud cover, or were thought to be formed prior to the lake formation stage.The lake relationships comprise: 516 Df lakes, of which 403 lakes were already in the Df type before the 1990s; 100 Cp lakes, of which 67 lakes were in the Cp type before the 1990s; 58 Dn lakes, of which 43 lakes were in the Dn type; and 14 Cs lakes, of which 10 lakes were not yet formed by the 1990s.One lake was changed from Cs to Cp, 69 lakes showed a history of changing from Cp to Df, and nine lakes changed from Df to Dn.There were no lakes showing a history of moving from Dn to Df or Dn, to Df, to Cp, or to Cs. Figure 5 plots the relationship between lake area ~1990 and expansion factors over the last 20 years with colors indicating changes in glacier-glacial lake connectivity.The expansion factors for the Cs lakes were the highest at 11.1 ± 12.3 times (n = 4) followed by those of the Cp-from-Cp lakes at 2.8 ± 3.5 times (n = 67) (Figure 5a).The Df-from-Df lakes had smaller mean expansion factors and a smaller standard deviation of 1.3 ± 0.7 times compared to those of the Df-from-Cp lakes (1.9 ± 5.2 times) (Figure 5b,c).The Dn-from-Df lakes and the Dn lakes showed only minor variation in their expansion factors concentrated around 1.0 (Figure 5d).

Potential Flood Volumes
The analysis revealed that 169 lakes have zero PFV and 226 lakes have positive PFV ranging from 0.04 to 33.7 million m 3 (Figure 1b).Other lakes were not analyzed by [14].Lakes with larger PFVs (>5 million m 3 ) were found in all nine basins: one in Amo Chu; three in Chamkhar Chu; three in Dangme Chu; eight in Kuri Chu; five in Mangde Chu; three in Mo Chu; three in Pho Chu; six in Tibet; and one in Wang Chu.They were sorted from the largest to the smallest PFVs (Table 4).Expansion rates, lake areas, ID, PDGL (or not), connectivity, moraine-damming or bedrockdamming, presence or absence of debris-cover, and glacial features are also shown.Table 4 also shows the rest of the ICIMOD-suggested PDGLs in a descending order of PFVs.These PDGL comprise 11 lakes with positive PFVs and 11 lakes with zero PFV.

Potential Flood Volumes
The analysis revealed that 169 lakes have zero PFV and 226 lakes have positive PFV ranging from 0.04 to 33.7 million m 3 (Figure 1b).Other lakes were not analyzed by [14].Lakes with larger PFVs (>5 million m 3 ) were found in all nine basins: one in Amo Chu; three in Chamkhar Chu; three in Dangme Chu; eight in Kuri Chu; five in Mangde Chu; three in Mo Chu; three in Pho Chu; six in Tibet; and one in Wang Chu.They were sorted from the largest to the smallest PFVs (Table 4).Expansion rates, lake areas, ID, PDGL (or not), connectivity, moraine-damming or bedrock-damming, presence or absence of debris-cover, and glacial features are also shown.Table 4 also shows the rest of the ICIMOD-suggested PDGLs in a descending order of PFVs.These PDGL comprise 11 lakes with positive PFVs and 11 lakes with zero PFV.

Current Status of Glacial Lakes
The current study analyzed 733 lakes in the Bhutan Himalayas, which included 278 lakes already delineated by [21].The largest number of glacial lakes (58%) were in the size class of 0.01-0.05km 2 (Figure 3a).Nie et al. [2] also showed that the highest number of lakes were in the size category of 0.01-0.05km 2 .This size category had the greatest total area of glacial lakes in the Bhutan Himalayas and the Himalayan Range.This indicates that glacial lakes in the Bhutan Himalayas are not different from the rest of the Himalayas in terms of size distribution.Veettil et al. [24] delineated 1065 glacial lakes from a much smaller domain than the current study's domain.This was a larger number than the current study because smaller lakes (<0.01 km 2 ) and lakes >2 km away from the Little Ice Age (LIA) moraine were included.
The vertical distribution of the glacial lakes varies between the nine basins (Table 2, Figure 4).The southern basins have lower lake elevations (i.e., the lowest mean elevation of 4571 m a.s.l. was found for the Wang Chu Basin followed by the Mo Chu and Pho Chu Basins with mean elevations between 4100 and 4200 m a.s.l.).In contrast, the Tibetan Basin had the highest mean elevation of 5137 m a.s.l. with individual lakes above 5800 m a.s.l.This south-to-north gradient reflects elevation differences in the glacier termini, as shown by [16].The elevation of the glacier terminus in the Bhutan Himalayas is higher than in the Tibetan Basin to the north and lower in the six southern basins, with the transverse basins (Kuri Chu and Dangme Chu) in between (Figure 10c in [16]).The Kuri Chu Basin had the largest standard deviation of the elevations (348 m), possibly because of a mixture of contrasting cold and arid climates in the north and warm and humid conditions in the south.The Pho Chu Basin had the second largest standard deviation of 279 m.There were two major elevation distribution groups, one around 5000-5100 m a.s.l., and the other around 4300-4400 m a.s.l.(Figure 4e).The higher group consisted of lakes from geographically neighboring basins.The lakes in the lower group were located in the Lunana region (Figure 2a; around (A) in Figure 1b), which includes the rapidly expanding Thorthomi and Lugge Lakes that caused GLOF in 1994.Relatively large glacial lakes connected to debris-covered glaciers are located in this region.The thermal insulation effect of debris mantles contributes to the debris-covered termini being located at a lower elevation than the debris-free termini [16,26,27].

Interaction and Integration of GLOF Scale Multiple Indices
Figure 6 shows the normalized frequency distributions for the expansion factors of the Cp, Df, and Dn classes, where an expansion factor of unity indicates no change over the course of approximately 20 years.Cp (in red) has a long tail toward large expansion factors followed by Df (in blue).In contrast, Dn is highly concentrated around 1.These findings strongly suggest that lake expansion is related to the retreat of associated glaciers.Thus, the connectivity between glaciers and glacial lakes, whose change is a proxy for glacial retreat, is critical for the expansion of glacial lakes.
As discussed in Section 5.1, debris-cover is another critical factor in lake expansion.Figure 5 illustrates the distributions of debris-covered glaciers for different connectivity groups.Most of the lakes which have expansion factors exceeding two (more than double in size) and lake sizes exceeding 0.1 km 2 belong to the group of Cp-from-Cp (Figure 5a).Furthermore, many of these lakes are associated with or connected to debris-covered glaciers.The lakes currently classified as Df, which comprise a large percentage of all the lakes (473 out of 733), showed large expansion, up to a 10-fold.However, these lakes are not larger than 0.1 km 2 nor debris-covered.This is strong evidence for the importance of debris-covered mother glaciers for rapid expansion of large glacial lakes, which is consistent with previous studies [1][2][3]28].
Finally, PFVs in terms of the four volumetric classes were added to Figure 5.The lakes in the largest PFV class with volumes exceeding 10 × 10 6 m 3 are also large lakes (>0.1 m 2 ).However, many of these lakes are not rapidly expanding and only two of them have expansion factors exceeding 2. ) where n is an integer.

Synthesis of Multiple Criteria
Previous studies have assessed the amount of water released by Sabai Tsho in Nepal in September, 1998, as 17.7 × 10 6 m 3 [29]; 17.2 × 10 6 m 3 by Luggye Tsho in Bhutan in November, 1994 [30], and 5.0 × 10 6 m 3 by Dig Tsho in Nepal in August, 1985 [7].Based on these observations, Table 4a shows integrated information of the GLOF-related factors discussed so far in a descending order of PFV, with a minimum value of 5.0 × 10 6 m 3 .In total, there were 33 lakes, with a maximum PFV of 33.7 × 10 6 m 3 by Lake Ph-016.Twelve lakes had PFVs exceeding 1.0 × 10 7 m 3 .The expansion factors ranged from 0.9 to 4.2, with only two lakes greater than 1.5.In terms of connectivity: 17 lakes were in the Dffrom-Df group; seven in Cp-from-Cp; six in Df-from-Cp; one for Dn-from-Dn; and two were unidentified.Of these 33 lakes, nine were associated with glaciers with debris-cover.In fact, some of them had a large debris-cover ratio, up to 24% of the glacier surface.Visual geomorphological interpretation with Google Earth identifies eight lakes dammed by bedrock and the others dammed by moraines.No remarkable relationship is found between this categorization and other indices (Table 4).
Having complied this information, it was tempting to extract lakes with high GLOF scales and potential.Using the four criteria: large PFV, high rate of expansion factor, currently either Cp or Cf and thus hydrologically feasible for expansion, and connected with a debris-covered glacier, two lakes emerged; Ph-051 and Da-021 had large PFVs (>15 million m 3 ), large expansion factors (>3 times), and had been in contact with debris-covered glaciers (Cp-from-Cp).
The former (Ph-051) is Thorthomi Lake located in the Lunana region in the Pho Chu Basin (Figures 1a,c and 2b).Numerous small supra-glacial ponds on the debris-covered terminus were

Synthesis of Multiple Criteria
Previous studies have assessed the amount of water released by Sabai Tsho in Nepal in September 1998, as 17.7 × 10 6 m 3 [29]; 17.2 × 10 6 m 3 by Luggye Tsho in Bhutan in November, 1994 [30], and 5.0 × 10 6 m 3 by Dig Tsho in Nepal in August, 1985 [7].Based on these observations, Table 4a shows integrated information of the GLOF-related factors discussed so far in a descending order of PFV, with a minimum value of 5.0 × 10 6 m 3 .In total, there were 33 lakes, with a maximum PFV of 33.7 × 10 6 m 3 by Lake Ph-016.Twelve lakes had PFVs exceeding 1.0 × 10 7 m 3 .The expansion factors ranged from 0.9 to 4.2, with only two lakes greater than 1.5.In terms of connectivity: 17 lakes were in the Df-from-Df group; seven in Cp-from-Cp; six in Df-from-Cp; one for Dn-from-Dn; and two were unidentified.Of these 33 lakes, nine were associated with glaciers with debris-cover.In fact, some of them had a large debris-cover ratio, up to 24% of the glacier surface.Visual geomorphological interpretation with Google Earth identifies eight lakes dammed by bedrock and the others dammed by moraines.No remarkable relationship is found between this categorization and other indices (Table 4).
Having complied this information, it was tempting to extract lakes with high GLOF scales and potential.Using the four criteria: large PFV, high rate of expansion factor, currently either Cp or Cf and thus hydrologically feasible for expansion, and connected with a debris-covered glacier, two lakes emerged; Ph-051 and Da-021 had large PFVs (>15 million m 3 ), large expansion factors (>3 times), and had been in contact with debris-covered glaciers (Cp-from-Cp).
The former (Ph-051) is Thorthomi Lake located in the Lunana region in the Pho Chu Basin (Figure 1a,c and Figure 2b).Numerous small supra-glacial ponds on the debris-covered terminus were integrated into several large lakes during the second half of the 20th century [5,6].Because the area is accessible and there is concern over such large supra-glacial lakes, the United Nations Development Programme supported the Government of Bhutan to install a water drainage outlet for this lake [31][32][33].Since 2009, about 350 local workers, every summer, have attempted to control the water level.Thorthomi Lake had higher multiple criteria related to a large-scale GLOF compared to other lakes in the Bhutan Himalayas.Attempting to control water levels is a reasonable step toward the prevention of GLOF hazards.
The latter (Da-021) is an unnamed lake connected to a north-facing debris-covered glacier in the Dangme Chu Basin, China (Figure 1b,d and Figure 3b).Further remote sensing and if necessary in situ observations are needed to assess the possibility of a large-scale GLOF event at this location.Flood simulation using numerical models with accurate downstream topographic data would be ideal to quantify possible damage as demonstrated in [34][35][36].
These two lakes were selected in the current study based on four criteria.However, other glacier-related factors are important for lake expansion.For example, calving of a glacial terminus, wind effects, englacial conduits, and glacial flow velocities could affect lake expansion [37].Further in situ and remote-sensing studies focused on these factors would help the understanding of why these two lakes have been rapidly expanding compared to other lakes.When the study region is extended from the Bhutan Himalaya to others, a variety of local topography should be considered in terms of topographic influence on lake expansion, because the north-to-south asymmetric feature in the Bhutan Himalayas illustrated in [38,39] has an influence on glacier advance/retreat.
Finally, it was noted that five of the 33 lakes in Table 4a were identified as PDGLs by the ICIMOD study [10].Table 4b lists another 19 ICIMOD PDGLs, of which six lakes have positive PFVs with a maximum of 2.7 × 10 6 m 3 (Ph-027).Of these lakes, Mo-015 has a relatively large expansion factor (3.5-fold) but the present lake area is small (0.06 km 2 ) and has been separated from the mother glacier (Df-from-Df).The other 11 PDGLs have zero PFVs, which means adequately low water levels relative to the foot of the moraine and/or gentle hillslope.

Conclusions
The current study developed a comprehensive glacier-glacial lake inventory for the Bhutan Himalayas based on satellite observations between 1987 and 1990 by SPOT-1 and ~2006-2011 by ALOS.This new integrated inventory provides much-needed information for a better understanding of glacier-related processes located in the eastern part of the Hindu Kush Himalayas.As a highly sought application of the inventory, GLOF scale and potential was analyzed using multiple criteria.In addition to the conventional criteria of expansion rate and lake size, the connectivity characteristic of the glacial lake to a mother glacier, potential flood volume, and debris-cover of the mother glacier, if any, for each of the glacial lakes in the inventory were examined.
With the availability of the two observational windows, not only the current connectivity condition but also its evolution was evaluated.The details of connectivity history using four different states: supra-glacial (Cs), pro-glacial (Cp), glacier-fed (Df), and no-glacier-fed (Dn) conditions were evaluated.The analysis of the connectivity evolution revealed that the evolution of the types Cp-from-Cs and Cf-from-Cp were more common, indicative of hydrological changes associated with the retreat of glaciers.
It was found that the majority of the lakes with high expansion rates (>2) and large size (>0.1 km 2 ) satisfied the following two conditions: (i) Cp-from-Cp, thus continuously in contact with mother glaciers; and (ii) being associated with debris-covered glaciers.A close examination based on multiple criteria including connectivity and debris-cover revealed that there were two lakes with high GLOF scale and potential.One of them was the much-publicized Thorthomi Lake, which justifies ongoing drainage efforts.The other is an unnamed lake in China, which may require further investigation.
This study is a step toward a better understanding of glacial processes, especially those related to GLOF scale and potential.The current study demonstrated the usefulness of the new integrated glacier-glacial lake inventory as a valuable tool for objective GLOF scale comparison.A potential for further studies is that possible trigger mechanisms for GLOF events, such as moraine-dam failure, mass movement entering a lake, and the wave-erosion or self-destruction of a terminal moraine, were not explicitly discussed [29,30].Future work needs to establish links between GLOF scale and GLOF potential and trigger mechanisms using observations and well-designed simulations.A comparison of local topography among different regions might explain regional difference of lake expansion trends.

Figure 3 . 19 Figure 4 .
Figure 3. Histograms of inventoried glacial lakes in the Bhutan Himalayas.(a) Number of lakes per basin.(b) Total area for each area class.Geosciences 2017, 7, 77 8 of 19

Figure 4 .
Figure 4. Vertical distribution of glacial lakes.(a) Total, and (b-j) the nine river basins.

Figure 5 .
Figure 5. Area of glacial lakes ~2010 versus expansion factors from ~1990 to ~2010.Dot colors correspond to connectivity changes initially separated in (a) present Cs and Cp, (b) Df from Df, (c) Df from Cp, and (d) present Dn.PFV, connection to debris-covered glacier, and PDGLs are noted by symbols.Two lakes with proposed multiple high criteria for large-scale glacial lake outburst flood (GLOF) annotated as (A) and (B) in (a) correspond to those in Figure 1.

Figure 5 .
Figure 5. Area of glacial lakes ~2010 versus expansion factors from ~1990 to ~2010.Dot colors correspond to connectivity changes initially separated in (a) present Cs and Cp, (b) Df from Df, (c) Df from Cp, and (d) present Dn.PFV, connection to debris-covered glacier, and PDGLs are noted by symbols.Two lakes with proposed multiple high criteria for large-scale glacial lake outburst flood (GLOF) annotated as (A) and (B) in (a) correspond to those in Figure 1.

Figure 6 .
Figure 6.Normalized frequency distribution of expansion factors between 2006 and 2010 and ~1990 smoothed using a moving window (w 10) where n is an integer.

Figure 6 .
Figure 6.Normalized frequency distribution of expansion factors between 2006 and 2010 and ~1990 smoothed using a moving window (w = 10 n 10 ) where n is an integer.

Table 2 .
Statistical and vertical distribution of glacial lakes.

Table 3 .
Summary of glacier-glacial lake connectivity.

Table 3 .
Summary of glacier-glacial lake connectivity.