Bhutan is a landlocked country bordered by China in the north and India in the south. The landscape is mainly dominated by high mountains and rugged terrain with an elevation ranging from 200 to 7,000 m a.s.l. (Mool et al., 2001). It has a total area of 38,394 km² and, as of 2018, a population of 734,374 (NSB, 2019). It consists of three broad climatic zones; subtropical in the southern foothills, temperate in the central valleys and inner hills, and alpine in the northern parts. The south-western monsoon originating in the Bay of Bengal governs the summer monsoon (Mool et al., 2001), which starts in late May and lasts until August. In 2017, the total annual rainfall was 1,816.74 mm, while the average yearly maximum and minimum temperatures are 22.3 and 12.3°C, respectively (NCHM, 2018). Figure 1 shows how the monthly mean, minimum and maximum temperatures and precipitation all increase from January to reach their peak values in July or August before falling in the remainder of the year.

Numerous glacier-fed rivers originate from the northern Himalayan mountains, providing essential ecological and livelihood services to settlements downstream. The rapid flow of these rivers coupled with the rugged terrain also presents an excellent opportunity to generate hydropower electricity (Tshering and Tamang, 2004; Farinotti et al., 2019). The river system in Bhutan is commonly divided into eight basins: Amochu, Wangchu, Phochu, Mochu, Chamkharchu, Mangdechu, Kurichu, and Drangmechu. Three of these river systems, namely Amochu, Kurichu, and Drangmechu, are transboundary. Amochu, which flows along the most western corner of Bhutan, originates in the TAR. In the east, the Kurichu flows from TAR, while the Drangmechu flows partly from TAR and partly from the state of Arunachal Pradesh in India (Katel et al., 2015).

The Bhutan Himalaya (Figure 1) contains 1,583 glaciers (1,487 ± 235 km²), of which 219 are debris-covered glaciers (951 ± 193 km²), and 1,364 are clean-ice glaciers (536 ± 42 km²) (Nagai et al., 2016). The debris-covered glaciers account for about 64% of the total. The total glacier cover accounts for about 4% of the total land cover in Bhutan (Nagai et al., 2016). In the Bhutan Himalaya, the glaciers have suffered an area loss of 23.3 ± 0.9% from 1980 to 2010 (Bajracharya et al., 2014), and, from 2000 to 2010, a mass loss ranging from 0.14 to 0.25 ± 0.13 m w.e.a⁻¹, with a thinning rate of 0.50 ± 0.14 m a⁻¹ (Gardelle et al., 2011). In-situ measurements have revealed that the Bhutanese glaciers have experienced comparatively higher mass loss than glaciers in the neighboring eastern Himalayan and southeastern Tibetan Plateau regions (Tshering and Fujita, 2016). Unfortunately, it is projected that the Bhutanese glaciers will undergo continuous and more rapid melting in the future, in response to current climate warming (Rupper et al., 2012). The recent inventory of High Mountain Asia lists 1,576 glacial lakes with a total area of 127.97 km².

Sentinel-2 MSI imagery with 10 m ground resolution provides superior earth surface information to Landsat imagery or any other freely available optical imagery (Copernicus, 2021a). For example, the band combinations and indices available in Sentinel-2 imagery, such as NDWI, false and natural composite, are nine times more detailed than corresponding band combinations from Landsat-8 (Nagy and Andreassen, 2019). The combination of the Sentinel-2A and Sentinel-2B constellation is also superior in terms of temporal resolution, with its 5-day revisit cycle at the equator and once in 2–3 days in mid-latitudes (Copernicus, 2021a ). Here, Sentinel-2 multispectral image, Level-1c, top-of-atmosphere (TOA) reflectance corrected data was used for mapping glacial lakes from 2016 to 2020. The data were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/). We chose images from November and October when the glacial lakes reach their post-monsoon maximum extent. This period also has the advantage of having low perennial snow and cloud coverage (Zhang et al., 2015), making the image selection process easier.

The Sentinel-1 mission consists of a constellation (A and B) of two polar-orbiting satellites, each with a 12 days revisit cycle, thus enabling collective temporal coverage every 6 days. It operates day and night and performs C-band synthetic aperture radar imaging, allowing image acquisition regardless of weather conditions (Copernicus, 2021b). Here, Sentinel-1 data was used to compliment the Sentinel-2 data to enable the mapping of glacial lakes from 2016 to 2020 which were missed by Sentinel-2 due to cloud coverage and other associated complications. The interferometric wide (IW) swath mode GRD-1 products with a spatial resolution of 20 m × 22 m in the range and azimuth directions were used. We chose VV polarization, as the water surface is more sensitive to the co-polarized wave (Wangchuk et al., 2019).

Declassified satellite data, Corona KH-4, was used to map historical glacial lakes from the 1960s. The Corona KH-4 system carried two panchromatic cameras and produced high-resolution images (∼2.7 m) between 1962 and 1972 (U.S. Geological Survey, 2008). After declassification of the Corona satellite data in 1995, KH-4 has been used in various environmental change studies (Song et al., 2015; Goerlich et al., 2017). Because of its high ground resolution, it has also been used in mapping historical glacial lakes such as in the northern Tien Shan (Bolch et al., 2011) and Poiqu River basin (Central Himalaya) (Zhang et al., 2019). KH-4 data covering the entire Bhutan Himalaya was downloaded from the USGS online portal (https://earthexplorer.usgs.gov/). Cloud cover was the major challenge while selecting the KH-4 data, but as our study area was relatively small, we were able to find multiple numbers of tiles covering the same area, which meant that it was not too difficult to obtain appropriate images. A few lakes could not be identified in the KH-4 data due to extreme cloud coverage in the available imagery. In such instances, we used hexagon KH-9 data to fill the gaps. The KH-4 data were from 1962 to 1972, while the KH-9 data were from 1973 to 1974, so there was not a substantial time difference between these two image sets.

We used the SRTM-30 m digital elevation model (DEM) as the primary remote sensing data to determine the geomorphometric characteristics of the lake surroundings for GLOF hazard assessment (Huggel et al., 2002; Huggel et al., 2003; Romstad et al., 2008; Wang et al., 2015). Specifically, we employed the DEM to quantify the lake-surrounding topographic potential for landslide or ice/snow avalanche into the lake and to determine the steep lakefront area (SLA). SRTM-30 m was also used as ancillary data for glacial lake mapping from 2016 to 2020. The SRTM-30 m data were available from the USGS online portal system.

We used glacier polygons to determine the lake-glacier distance threshold and avalanche-prone areas. The Randolph Glacier Inventory version 6.0 (RGI Consortium, 2017) and GAMDAM (Sakai, 2018) provides a glacier data set covering the high mountains throughout the world. However, its accuracy is relatively low when considering highly localized areas such as the Bhutan Himalaya, and so, for avalance mapping, we used the glacier inventory of Bhutan by Nagai et al. (2016) which is purely focused on the Bhutan Himalaya and has better local accuracy. However, owing to their multi-temporal coverage, the RGI v6.0 glacier data sets (RGI Consortium, 2017) were used to determine the threshold lake-glacier distance. Table 1 presents a summary of all the data used in the study.

We used the GlakeMap python package (Wangchuk and Bolch, 2020) to map glacial lakes from 2016 to 2020. GlakeMap is an automated machine learning-based method developed using the sci-kit-learn, snappy, and arcpy python libraries. The package enables the complementary combination of multi-source source data such as normalized difference water indices (blue and green) and radar backscatter and reaps the benefits of the strength of each data source. Such complementary combinations are principally chosen through rules-based image segmentation and the machine learning algorithm “random forest classifier.” The proposed method has been tested in several high mountain areas worldwide, including the Bhutan Himalaya, and has achieved average detection and delineation accuracies of 97.87 and 98.95%, respectively. Firstly, NDWI_(green) and NDWI_(blue) from Sentinel-2 were calculated using Eqs 1 and 2, and then the Sentinel-1 data were processed to calculate radar backscatter. We performed rules-based image segmentation of NDWI indices and radar backscatter to derive the rough polygons for all the glacial lakes. Six predictor data sets were prepared: NDWI (green), NDWI (blue), NIR infrared band, radar backscatter, slope, and compactness ratio (calculated using Eq. 3). The glacial lakes were then delineated based on the initial lake polygons derived from image segmentation and attributes in the predictor data sets by the random forest classifier model. All these procedures were executed automatically using the GlakeMap package. Wangchuk and Bolch (2020) provide more methodological details. The glacial lake polygons were then laid over the false-color composite of Sentinel-2 and manually verified. Further validation of the lake polygons was done using high-resolution imagery available in Google Earth. N D W I ( b l u e ) = ( B l u e   −   N I R ) ( B l u e   +   N I R ) (1) N D W I ( g r e e n ) = ( G r e e n   −   N I R ) ( G r e e n + N I R ) (2) C R = 4 π A p 2 (3) where CR is compactness ratio, A is the area of a lake, and p is the perimeter of a lake.

We used the Corona KH-4 data to map glacial lakes in the 1960s. These images are not geometrically corrected and suffer from a significant geometric distortion. The images were co-registered by collecting 90–100 ground control points from the Sentinel-2 data using a geo-referencing toolset in the ArcGIS software. The high (10 m) resolution of the Sentinel-2 data makes it more suitable to use for this purpose than other freely available remote sensing imageries. The KH-4 images are single band, panchromatic images (Dashora et al., 2007), so the glacial lakes were manually digitized based on expert judgment. However, a single band image does provide substantial contrast between the lake and surrounding areas due to its high resolution. When the scenes were complicated by cloud cover and shadows cast by mountains, alternative scenes were used from other months or years. For the process of finding alternative images, the multiple numbers of images available between 1962 and 1970 covering the Bhutan Himalaya was advantageous, but when no alternative images were found, we used Hexagon KH-9 data to fill the data gap. A total of three hexagon KH-9 were used, and 78 lakes, missing in the KH-4 imagery, were mapped from the KH-9 data. The glacial lake mapping procedure is summarized in Figure 2.

The glacial lakes are mostly the meltwater from current or past glaciers and ice fields/caps accumulated on a favorable topography (Dimri et al., 2016; Otto, 2019). Here, glacial lakes were identified as any lake located within the 10 km buffer distance from the glacier terminus following the earlier studies (Wang et al., 2013; Zhang et al., 2015). The RGI v6.0 glacier inventory data (RGI Consortium, 2017) was used to determine this threshold distance. Based on the hydrological relationship between the lake and the nearest glacier terminus identified from high-resolution Google Earth images, the glacial lakes were classified into four categories: non-glacier-fed, glacier-unconnected, glacier-connected, and supraglacial lakes (Ageta et al., 2000). The minimum lake area delineated was chosen to be 0.003 km², as the main purpose of this inventory was to produce baseline glacial-lake data for the Bhutan Himalaya.

Since the lake polygons were manually corrected using the Sentinel-2 data after the automatic mapping process, we estimated the error associated with manual mapping from Sentinel-2. Based on the quality of satellite imagery, the error area of the glacial lakes is about 0.5 pixels (Fujita et al., 2009). Because the pure lake water pixels are usually surrounded by a mixture of water and non-water pixels, the maximum area error of glacial boundary extraction is estimated roughly to be half the area of the edge pixels (Wang et al., 2020). The manually delineated lakes are assumed to have a regular/Gaussian distributed error area. So, the lake error within one standard deviation (1σ) was calculated using the equation (Hanshaw and Bookhagen, 2014). E r r o r ( 1 σ ) = P G × G 2 2 × 0.6872 (4) E = E r r o r ( 1 σ ) A × 100   % (5) where p is the perimeter of the glacial lake (m), G is the spatial resolution of a satellite image, 0.6872 is the coefficient under 1σ, E is the relative error of the glacial lake, and A is the total area of the glacial lake (Wang et al., 2020).

GLOF hazard assessment is mainly concerned with quantifying the likelihood and magnitude of occurrence of GLOF (Allen et al., 2016; Wang et al., 2015). Indicators such as lake area and volume, mass movement entering into glacial lakes, the moraine width to height ratio, the height of freeboard, lake geometry, and distance between the lake and the glacier terminus are commonly used (Bolch et al., 2011; Emmer and Vilímek, 2013; Worni et al., 2013; Aggarwal et al., 2017; Nagai et al., 2017; Frey et al., 2018). However, for first-order hazard assessment, attributes such as causes of historical GLOF, availability and applicability of remote sensing data, and suitability of automation are commonly employed (Bolch et al., 2011) (Supplementary Table S2). Accordingly, here we selected the following hazard assessment factors: 1) ice/snow avalanche, 2) landslide/rockfall, 3) steep lakefront area (SLA), 4) lake area expansion 5) upstream GLOF, 6), glacier-lake distance 7) lake area. Assuming that the differing intensity of the particular factor poses a varying degree of GLOF hazard to the lake, the quantified value of each element was categorized into three or two alternative classes: high, medium, and low (Table 2). Based on the previous literature, we assigned an index value of 1 to high, 0.5 to medium, and 0.25 to the low alternative classes (Bolch et al., 2011). A summary of this GLOF hazard assessment process is provided in Figure 2.

We calculated the GLOF hazard factors 1) ice/snow avalanche and 2) landslide/rockfall using the topographic tsunami generating potential of the lake surrounding terrain. This method identifies any potential release cells within the vicinity of a lake, considering slope as the single most important criterion. It is based on the assumption that the probability of detachment of a mass of tsunami generating potential is higher on the steep slopes directly above the lake than on more distal slopes. Initially developed by Romstad et al. (2008), the model has been used to map topographic potential for glacial lakes in the Swiss Alps (Schaub, 2015), northeast India (Allen et al., 2016), Tibetan Plateau (Allen et al., 2019) and the Third Pole (Zheng et al., 2021). Here, following Allen et al. (2019), the topographic potential area was defined as a 30° slope within the lake’s watershed area and with a threshold reach angle of 15°. The slope of the area surrounding the lake was derived from the SRTM-30 m data. The pixels within the potential topographic area which are covered by glacier were then identified as the avalanche-prone area, while those without glacier covering were considered to form the landslide/rockfall prone area. An avalanche area of more than 0.1 km² was assigned a high index value, 0.05–0.1 km² a medium value, and less than 0.05 km² a low value. Here, the lower limit avalanche of 0.05 km² was considered at par with the minimum size of the lake selected for this hazard assessment assuming that an avalanche of any volume may cause GLOF magnitude of the same volume. Because large landslide/rockfall areas were widespread across the sampled lakes, with for example, 78.4% of the sampled lakes found to have a landslide area of 0.1 km², while only 21.9% of the sampled lakes were found to have a simillarly sized avalanche area, we followed a slightly different weighting categorization for landslide-zone quantification. Areas larger than 0.5 km² were assigned a high value, those in the range 0.1–0.5 km² a medium value and those less than 0.1 km² a low value.

The type of damming material and its geometry determine a dam’s stability (Wang et al., 2011). The dam width-to-height ratio, dam crest width, presence of permafrost, buried ice, and freeboard height (Worni et al., 2013; Worni et al., 2014) are widely used to determine dam stability. However, identifying these dam parameters is difficult with existing freely accessible satellite data. Fortunately, the steep lakefront area (SLA), which is based on the depression angle between the flat lake surface and the surrounding terrain, is an excellent alternative measurement of the width-to-height ratio. The steepest elevation angle from a particular point towards the lake surface is considered as a value for SLA (Fujita et al., 2013). The SLA has been widely considered for hazard assessment, for example in the Himalayas in India (Dubey and Goyal, 2020) and Nepal (Rounce et al., 2016). Here, we identified the SLA by following the model developed by Fujita et al. (2013), which uses the minimum threshold angle of 10° within the buffer distance of 1 km from the lake shoreline. We gave a high index value to lakes with an SLA and low values to those without. However, bedrock-dammed lakes were considered to be stable (Huggel et al., 2002) and were, therefore, assigned a low hazard value even if they had an SLA.

Lake expansion is one of the crucial criteria for hazard assessment (Bolch et al., 2011; Anacona et al., 2014; Rounce et al., 2016; Allen et al., 2019) as it influences the volume of water available for GLOF (Richardson and Reynolds 2000). Glacial lake expansion also intensifies the lake’s vulnerability to potential hazards such as mass movement (Rounce et al., 2016). The threshold expansion percentage used by previous studies was based on the length of the study periods. For example, a threshold of more than 100% expansion was employed for a 40°years study period in the Sikkim Himalayas (Aggarwal et al., 2017), while a threshold of 20–50% expansion over a 20°year period was used in Nepal Himalayas (Khadka et al., 2018). We did not use glacial-lake area change from the 1960s to 2020, despite the data being available, because some glacial lakes have expanded significantly since the 1960s but have now attained their full extent and stopped expanding. Instead, when the expansion within the last 5 years (2016–2020) was greater than the uncertainty range, we assigned a high index value. When the growth was less than the uncertainty range or the lake had undergone recession a low index value was assigned.

A GLOF from upstream lakes is one form of mass movement that can cause a breach of the main lake downstream. Because it has high erosion and transport potential (Breien et al., 2008), it can produce displacement or tsunami waves equivalent to those caused by slope moments of the same magnitude (Emmer and Cochachine, 2013). In Bhutan, the most recent impactful GLOF of Lemthang Tsho was caused by the breaching of upstream supraglacial lakes (Gurung et al., 2017). We considered any upstream lake of size ≥0.01 km² as a potential danger to downstream glacial lakes. According to Emmer and Cochachine (2013), all glacial lakes are hazardous as they are susceptible to at least one form of GLOF scenario. Therefore, we considered all upstream glacial lakes as a potential GLOF hazard source to the downstream lake without performing separate hazard assessments. We assigned a high index to lakes with ≥3 upstream lakes, a medium index to lakes with 1–2 upstream lakes, and a low index to lakes with no upstream lakes.

Glacial lakes expand due to the addition of melt-water from ice or thinning of the moraine dam and by a cumulative discharge of water from the ablation and calving of the glacier at the proximal end of the lake (Richardson and Reynolds, 2000). The distance and slope between the lake and the mother glacier determine the lake-glacier interaction (Wang, 2011). GLOFs in the Himalayas (Khadka et al., 2021) and on the Tibetan Plateau (Wang et al., 2011) have been within 0–700 and 0–800 m lake-glacier distance, respectively. In the Bhutan Himalaya, out of 13 identifiable lakes with a GLOF history, at least two lakes have a lake-glacier distance of 1,000 m. Therefore, here, we used a maximum threshold lake-glacier distance of 1,000 m. This is also the value used by Allen et al. (2019) for the region-wide hazard assessment of lakes in the Tibetan Plateau. We measured glacier-lake distance using the RGI v6.0 (RGI Consortium, 2017) glacier dataset. High, medium, and low indices were assigned for glacier-lake distances of 0 m (contact with glacier), 0–500 m, and 500–1,000 m, respectively.

Lake area is a widely used factor for assessing GLOF susceptibility, as larger lakes present a greater GLOF hazard than smaller ones. Larger glacial lakes are typically deeper with a larger volume, which means higher potential flood volume and greater hydrostatic pressure on the damming wall (Richardson and Reynolds 2000). Here, lake area was calculated from the lake polygon extracted from the Sentinel-1 and Sentinel-2 data. A glacial lake size of 0.01 km² is considered large enough to present a significant hazard to downstream settlements (Nagai et al., 2017). When considering regional scale hazard assessment, an area threshold of 0.1 km² is also reasonably used (Allen et al., 2019). Based on 51 past GLOF events in the Himalayan region, Nie et al. (2018) have recommended a minimum area threshold of 0.05 km² for hazard assessment. In the Bhutan Himalaya, out of 13 existing lakes with a GLOF history, 10 have an area ranging from 0.01 to 0.06 km². Therefore, we used a minimum area threshold of 0.05 km². We assumed that GLOF magnitude will increase with increasing area, and so we assigned indices as follows: low for lakes <0.05 km², medium for lakes in the range 0.05–0.1 km², and high for those with area >0.1 km².

The different hazard factors have varying degrees of contribution to the GLOF hazard potential of lakes. Thus, assigning a weight to each variable is deemed necessary (Wang et al., 2011) to realize the objective hazard classification of the lakes. The seven factors determined in this study were weighted using AHP (Table 3). AHP is a multicriteria decision-making technique that quantifies subjective expert judgment by assigning values through pairwise comparison. The pairwise comparison enables prioritizing the available alternative options by estimating their relative significance (Saaty, 2008). The AHP technique also ensures consistency of the subjective judgment (Saaty, 2008) which is challenging to maintain in other semi-quantitative methods of GLOF hazard assessment (Emmer and Vilímek, 2013). AHP has been widely used in GLOF hazard assessment worldwide and particularly in the Himalayas (Anacona et al., 2014; Aggarwal et al., 2017; Prakash and Nagarajan, 2017b; Khadka et al., 2021).

Here, our expert judgment for the priority ranking was influenced by the previous causes of GLOFs in the Bhutan Himalaya and its neighboring Himalayan regions and following the convention in the existing literature. We assigned the highest rank to ice/ice snow avalanche followed by landslide/rockfall. The majority past GLOF events in Bhutan (55.6%) (Komori et al., 2012; NCHM, 2019a) (Supplementary Figure S1) and Himalayas (53%) (Richardson and Reynolds, 2000) were caused by ice avalanches or glacier advances. Landslide/rockfall was assigned the second rank as it is the key component of mass movement. In some previous studies (Allen et al., 2016; Allen et al., 2019; Zheng et al., 2021), ice/snow avalanches and landslides are often treated as one integrated factor. This assumes that the landslide area may provide an alternative representation of the avalanche area in case of uncertainty resulting from the data and method. We also assigned a relatively high ranking to SLA, considering it a first-order indicator of dam stability (Richardson and Reynolds, 2000; Hegglin and Huggel, 2008; Emmer and Vilímek, 2013). Lake area was ranked lowest because our initial screening criteria already considered a minimum area of lake that may cause considerable damage to downstream settlements. Accordingly, the paired-wise comparison matrix was constructed, and weights for each factor (F _w ) were determined (Table 3) with an overall consistency ratio of 0.085 (8.85%). The final index of each factor (F _i ) was calculated as the product of class index (C _i ) and factor weight (F _w ) using Eq. 6. The hazard score of the lake was calculated as an average score of all factor index values using Eq. 7. Finally, the hazard score of each lake was normalized to a value ranging from 0 to 1. Here, based on the hazard score, the hazard potential of the lake was categorized into five groups, namely, very low (0–0.2), low (0.2–0.4), medium (0.4–0.6), high (0.6–0.8) and very high (0.8–1.0). F i = C i × F w (6) H a z a r d   s c o r e   o f   l a k e = S u m   v a l u e   o f   F i   N u m b e r   o f   f a c t o r s ( 7 ) (7)

Our GLOF hazard assessment relies on well-established methods and past GLOF experience in the context of the Bhutan Himalaya. However, validation of the hazard-assessment method was still deemed necessary as we relied solely on remote sensing data, which has inevitable uncertainties. GLOF is a high-magnitude, low-frequency natural disaster, thus validation is often tricky (Allen et al., 2019). However, the Bhutan Himalaya has experienced multiple GLOF events which were substantial enough to validate our method. The validation of hazard assessment methods through the use of past GLOF events has commonly been done in previous studies (Anacona et al., 2014; Emmer and Vilímek, 2014; Allen et al., 2019). A total of 23 GLOFs events have been documented in the Bhutan Himalaya (Komori et al., 2012; Gurung et al., 2017; NCHM, 2019a) (Supplementary Table S1). However, five of these events were located on the TAR side of the Himalayas and do not have any hydrological connection with the glacier-fed rivers that flow into the Bhutanese interior. Four of the events are from glaciers, and for another, the lake involved no longer exists due to the complete emptying of water during the GLOF, according to the verification carried out for this study. Out of the 13 remaining events, eight of them fulfilled our initial screening criteria (area ≥0.05 km² and glacier-lake distance ≤1,000 m) and were used for validation. Although eight lakes out of a total of 278 represents only 2.9%, they do include all the past GLOF events in the region that caused damage to downstream settlements.

Using Sentinel-2 and Sentinel-1 imagery, we mapped 2,553 (151.81 ± 7.76 km²), 2,566 (152.64 ± 7.83 km²), 2,572 (153.94 ± 7.83 km²), 2,569 (153.97 ± 7.79 km²) and 2,574 (156.63 ± 7.95 km²) glacial lakes from 2016, 2017, 2018, 2019 and 2020, respectively. A total of 2,187 glacial lakes corresponding to 130.19 ± 2.09 km² were mapped using historical KH-4 data from the 1960s (Supplementary Table S3). Of the total in 2020, 1,118 (82.92 km²) were glacier-fed lakes. In this inventory, the size of the glacial lakes ranges from 0.003 to 4.30 km², with a mean area of 0.061 km² and a standard deviation of 0.17 km². There were slightly more non glacier-fed lakes (56.6%) than glacier-fed lakes, but the glacier-fed lakes were slightly dominant in terms of areal coverage (52.9%). Thorthormi Tsho, located in the Phochu basin, is the largest glacial lake in this inventory.

We used the 2020 inventory to analyze the distribution of glacial lakes in the BTH (Figure 3). The glacier-fed lakes were chiefly present in the elevation range 4,250–5,500 m a.s.l and had a mean elevation value of 4,810 m a.s.l. Their distribution was most dense in the central part of the BTH, becoming more sparse towards the eastern and western ends. Basin-wise, the Pchochu basin has the highest proportion 28.18% (area) and 22.63% (number), followed by Kurichu 26.35% (area), 20.66% (number). Wangchu with 1.34% (area), 0.78% (number), and Amochu with 2.24% (area), 1.71% (number) basins had the lowest proportion of glacier-fed lakes. Although Drangmechu had the most glacial lakes of all types by both number and area, it only contained 14.7% (number) and 8.9% (area) of the total number of glacier-fed lakes.

A type-wise comparison of the glacial lakes shows that the glacier-unconnected lakes dominate both the number and area (76.9 and 56.7% of the total, respectively). Supraglacial lakes are uncommon, accounting for only 1% of the total area and 4.6% of the total number. Although glacier-connected lakes account for only 18.5% of the total numbers, they make up 42.3% of the total area. The mean area was the largest for glacier-connected lakes (0.17 km²), followed by glacier-unconnected lakes (0.055 km²) and supraglacial lakes (0.016 km²). We did not detect any supraglacial lakes in the Drangmechu basin (Figure 4).

Following Maharjan et al. (2018), we examined the glacial lake distribution for six different size classes (km²) (<0.01, 0.0–0.05, 0.05–0.1, 0.1–0.5, 0.5–1, and 1–5) as shown in Figure 4. Small glacial lakes (<0.1 km²) were most frequent in terms of number (84.34%) while larger glacial lakes (>0.1 km²) formed a major portion of the total area (70.7%). The most common size, in terms of numbers, was 0.01–0.05 (km²) (47.58% of the total), while the 0.1–0.5 range predominance in area (38.0%). The glacial lakes less than 0.01 km² had the lowest areal coverage (2.33%), and the size class 1–5 km² had the lowest count (0.62% of the total). There were no glacier-fed lakes in size categories 0.5–1 and 1–5 (km²) in Amochu, Wangchu, or Kurichu basins, while Drangmechu lacked any lakes in the 1–5 (km²) category.

Over the last 6 decades (1960s−2020), the total number of glacial lakes increased from 2,187 (130.19 ± 2.09 km²) to 2,574 (156.63 km²). The number of glacier-fed lakes increased by 190 (13.77 km²), while the number of non-glacier-fed lakes grew by 200 (12.67 km²). No significant change in the mean size of the glacial lakes was observed between the 1960s (0.059 km²) and 2020 (0.060 km²). The standard deviation of the area of glacial lakes increased slightly, from 0.144 km² in the 1960s to 0.17 km² in 2020. The maximum size of the glacial lakes increased from 3.4 to 4.3 km² over the same period. The number and area of glacier-fed lakes increased across all the basins, with the exception of Amochu and Drangmechu. The further analysis of glacial lake change behavior and the results presented in the following sections concerns only the glacier-fed lakes.

The number of glacial lakes remained stable in the Amochu basin, while in Drangmechu, the glacial lakes shrank by 6.81 km². The Phochu basin had the largest expansion, with 51 (10.43 km²) new lakes emerging over the period, followed by the Kurichu basin where 41 (4.54 km²) new lakes formed (Figure 5).

Looking at the changes in terms of glacial lake types, we observed the largest areal expansion (10.8 km²) in glacier-connected lakes, accounting for 75.4% of the overall total lake area growth in the BTH. However, 182 new glacier-unconnected lakes appeared over the period, contributing 96.3% of the overall increase in lake numbers in the BTH. The area of supraglacial lakes decreased by 35.47%, although their total number actually increased by 10.86% over the period (Figure 6).

The glacial lake changes across the different size classes revealed that smaller glacial lake classes had the largest increases in numbers, while the larger glacial lake classes had the largest areal expansion. For example, a total of 182 new glacier-fed lakes of area <0.1 km² were formed over the period, although their areal change contribution to the overall change was only 3.8 km². On the other hand, the glacier-fed lakes >0.1 km² expanded their area by 9.9 km² despite their number increasing by only eight (Figure 6).

The glacier-fed lake changes across the four elevation zones, namely 4,000–4,500, 4,500–5,000, 5,000–5,500, and 5,500–6,000 m a.s.l are as shown in Figure 7. Over the last 6 decades, we observed the highest relative expansion (96.35%) in the elevation range 4,000–4,500 m a.s.l, followed by the range 5,500–6,000 m a.s.l. In contrast, lake area in the 4,500–5,000 m a.s.l range suffered recession (by 13.20%). Overall, the mean elevation of glacier-fed lakes increased from 4,986.46 to 5,031 m a.s.l.

The glacier-fed lake changes across the three lake-glacier distance ranges (0–500, 500–1,000, and >1,000 m) are also shown in Figure 7. Over the period, the glacial lakes with a glacier-lake distance within the range 0–500 m expanded, while glacial lakes with glacier-lake distance in other categories recessed. However, absolute net expansion (30.0 km²) for lakes in the 0–500 m category was more significant than net recession (−16.34 km²) in the other lake-glacier distance ranges. The mean glacier-lake distance decreased from 691.3 to 405 m over the 60 year period.

In 2016–2020, the glacier-fed lakes expanded (by 5.79% in area, and 1.91% in number) while the non-glacier-fed lakes remained reasonably stable (Figure 5). The annual changes between 2016 and 2020 show the area expanding at the rate of 0.96 km² per year. A total of 19 (4.82 km²) new lakes were formed, with annual increments occurring except for 2018, when the total number of lakes decreased by three (Figure 5).

The number and area of glacier-fed lakes expanded between 2016 and 2020 in all basins except the Wangchu basin (Figure 6). The largest expansion occurred in the Phochu basin where six new lakes (1.55 km²) were formed. The Wangchu basin exhibited the smallest areal growth (0.03 km²), with its number of lakes decreasing by 1. However, there was irregular interannual variability in the number and area of lakes across the different basins. For example, in the Phochu and Wangchuk basins, the glacial lake counts increased in 2016–2017, remained stable in 2017–2018, and decreased between 2019 and 2020. By contrast, the Kurichu basin had a continuous annual expansion in area and number between 2016 and 2020 with average rates of 1.5 lakes and 0.2 km², respectively (Figure 6).

Glacier-fed lakes of all kinds, except supraglacial lakes, increased both in area and number. The number of supraglacial lakes increased during 2016–2017 (19.56%) and 2017–2018 (9.80%) but remained stable in 2018–2019. However, the area of supraglacial lakes expanded continuously between 2016 and 2020 at an average rate of 0.06 km² per year (Figure 6). The lakes across the different size categories exhibited an overall expansion trend during 2016–2020, although there is heterogeneous interannual variation (Figure 6). For example, the class of glacier-fed lake of area <0.01 km² increased between 2016 and 2017 and continuously shrank between 2017 and 2020. On the other hand, the 0.01–0.5 km² category expanded with an irregular annual pattern. The number of large glacial lakes (0.5–1 and 1–5 km²) remained stable, while their area underwent a continuous yearly expansion. We determined that the heterogeneous behaviour across the different classes was due to the rapid evolution of glacial lakes from small to large stages. For example, in 2017–2018, the drop of 24 in the number in the <0.01 km² category was comparable to the increase of 30 in the number in the 0.01–0.05 km² category (Figure 6).

Figure 7 shows that within the last 5 years, glacial lakes in all elevation ranges increased with a maximum relative expansion in the 4,000–4,500 m a.s.l range (10.48%) and a minimum (3.06%) in the 4,500–5,000 m a.s.l, range. However, there was inter-annual variability across the different elevation zones. For example, the glacial lakes in all years expanded in the elevation band 4,000–4,500 m a.s.l. In 2016–2017, the glacial lakes suffered recession in the 5,500–6,000 m a.s.l. elevation range. The maximum relative change was observed at 0–500 m glacier-lake distance category (7.26%), while there was a slight decrease in area in the >1,000 m category (−0.367 km²) in the last 5-year period. There was some interannual variability with glacial lakes in the <1,000 m glacier-lake distance category expanding in all the years, while those in the >1,000 m glacier-lake distance range, decreased in some years (2017–2018 and 2018–2019).

Table 4 shows the GLOF hazard score value and hazard level for the eight lakes that were used for the validation. Out of these eight lakes which have failed in the past and caused a GLOF, three fell in the high GLOF hazard-level category and three in the very high category. There was only one lake in each of the medium and very-low GLOF hazard categories. This validation established that our hazard assessment method was robust enough to use on other lakes in the BTH.

Based on the initial screening criteria, we identified 278 glacier-fed lakes for GLOF hazard assessment, out of which 84 lakes were transboundary (Figure 8). As per the GLOF hazard-level classification defined in the method section, we found 31 very high, 34 high, 50 medium, 102 low, and 62 very low GLOF-hazard lakes in BTH. The list of glacial lakes with high to very-high GLOF hazard levels is provided in Supplementary Table S4. The glacial lakes with high to very-high GLOF hazards were primarily located in the Kurichu basin (35.4%), followed by Phochu (18.5%), Drangmechu (15.4%), and Mochu (12.3%) basins. The number of glacial lakes with high to very-high GLOF hazard levels was smallest in the Amochu and Wanghcu basins. However, all the lakes with high to very-high GLOF hazard levels in Kurichu, Amochu, and Dranmgemechu were transboundary, located on the TAR side of the Himalayas. There were no lakes with a very-high GLOF hazard level in the Mangdechu basin. We found a positive correlation between the total number of glacial lakes and high GLOF hazard lakes. For example, Kurichu and Phochu basins exhibited the highest number of glacier-fed lakes (231 and 253, respectively), and also the highest number of very-high GLOF hazard-level lakes. Similarly, the opposite pattern was true for the Wangchu basin, which had the lowest number of glacier-fed lakes and lowest number of very-high GLOF hazard-level lakes.

Overall, the GLOF hazard-score distribution was slightly skewed towards the lower hazard levels, as shown in Figure 9. The very-high to high GLOF hazard-level lakes were mainly distributed in the 5,000–5,500 m a.s.l elevation range (47.69%) followed by the 4,500–5,000 range (35.38%). There was only one lake with a very-high hazard level in the 5,500–6,000 m a.s.l band. More than half (53.8%) of the lakes with a high to very-high hazard level fall in the size category of 0.1–0.5 km² while the size category of 0.05–0.1 accounts for about 27.69%. Across the different types of glacial lakes, more than half (56.9%) of very-high to high GLOF hazard lakes were glacier-connected lakes followed by glacier-unconnected lakes (41.53%). There was at least one high hazard level supraglacial lake (Figure 9).

Uncertainty in glacial lake mapping usually results from the data quality used for mapping and other factors such as cloud coverage and seasonal snow. The Sentinel-2 and Sentinel-1 data used here have the highest ground resolution of any freely available data. Although glacial lakes were initially mapped with an automatic method involving both the imageries, manual validation was done using a false-color composite of Sentinel-2. We calculated the uncertainty value based on the resolution of the Sentinel-2 images (10 m) following the method of Hanshaw and Bookhagen (2014). Our uncertainty range (±5.04% of total area) was lower than the values in previous studies [e.g., ±16.1% in Zhang et al. (2015), 13.47% in Wang et al. (2020), ±30% in Khadka et al. (2018), and 15.96% in Nie et al. (2017)] which used Landsat data. Salerno et al. (2012) projected that the error associated with lake area can increase by an order of magnitude when the imagery resolution increase from 10 to 30 m which is corroborated by the reasonably minimal uncertainty range in this study. This substantiates the fact that the glacial lake mapping using Sentinel-2 in this study provides a data set with low uncertainty over the BTH, which is essential in the context of glacial lake hazard assessment (Allen et al., 2019). Obstruction by cloud cover is one of the main challenges in glacial lake mapping in the Himalayas (Zhang et al., 2015), and cloud cover is almost perennial over the BTH. The average cloud cover for all Sentinel-2 scenes in the months of October and November was 56.50 and 61.27%, respectively. Here, the use of Sentinel-1, SAR data was helpful in mapping glacial lakes obstructed by cloud cover in the Sentinel-2 data. For example, 107 lakes missed in the 2019 Sentinel-2 imagery were mapped using Sentinel-1 data. Hence, we have reaffirmed the proposal of Wangchuk and Bolch (2020) that using Sentinel-1 as a complement to optical data is a valuable and novel approach to overcome the cloud-cover problem.

The number of glacial lakes mapped in this study was higher than in previous studies, as indicated in Table 5. We compared our 2018 glacial lake data set with that of Wang et al. (2020) in 2018 and that of Chen et al. (2021) in 2017, both of which form the latest available inventory covering the BTH. When we compare the number of glacial lakes with areas between 0.0085 and 0.05 km², Chen et al. (2021) reported 757 (16.59 km²), and Wang et al. (2020) mapped 970 (20.1 km²), while the current study found 1,475 (32.22 km²). The comparison indicates that the earlier inventories have missed glacial lakes, especially in the smaller size ranges, and implies that our study provides a more robust (both qualitative and quantitative) representation of the glacial lake census in the Bhutan Himalaya than previous studies at the regional scale.

Basin-wise, we found that glacier-fed lakes were predominant in the Phochu and Kurichu basins, in agreement with earlier distribution patterns (Ukita et al., 2011). Gardelle et al. (2011) suggested that highly glaciated regions are likely to have more glacier lakes. Here we found that Kurichu (709 km²) and Phochu (507 km²) have the largest areas of glacier and, correspondingly, the largest numbers of glacial lakes. An opposite trend was observed in the Wangchu (32.39 km²) and Amochu (12.8 km²) basins, which had the lowest numbers of glacial lakes. These trends indicate that glacial lake formation is primarily related to glacier dynamics, as claimed by Richardson and Reynolds (2000).

We observed the dominance of small glacial lakes (0.01–0.05 km²) in the BTH, which is in line with previous records (Ukita et al., 2011), and the general distribution pattern of glacial lakes in the entire Third Pole region (Zhang et al., 2015). However, glacial lakes in the larger size categories represent the major portion of the total area, in agreement with earlier findings in the Himalayas (Nie et al., 2017). Comparing the distribution patterns of glacial lake types, we found glacier-unconnected lakes were dominant both in area and number, followed by glacier-connected lakes. In contrast, there were very few supraglacial lakes. This distribution trend was also observed in the Nepal Himalaya (Khadka et al., 2018) and in the entire Himalayan region (Nie et al., 2017). Moreover, an earlier study across the Hindukush Himalayas also observed the sparse distribution of supraglacial lakes in the eastern Himalayas, which corresponds with our study site (Gardelle et al., 2011). This sparse distribution was shown to occur because the eastern Himalaya contain smaller-sized, less debris-covered glaciers than its western counterpart, which provide less favorable conditions for the formation and sustenance of supraglacial ponds (Gardelle et al., 2011; 2013). Using high resolution (3 m), Planet Labs satellite data for 2016–2018, 5,754 supraglacial ponds with an average size of 1,206 m² with high temporal and spatial variations were reported on three prominent glaciers in the Mangdechu and Chamkhar basins in the Bhutan Himalaya. This indicates that, because of the relatively smaller size and high variability of supraglacial lakes, they are usually unrepresented when the inventory is done using freely available satellite data with high to medium ground resolution.

Several researchers have studied glacial lake changes across the Himalayas. However, the previous studies are limited to the Landsat era, and typically restricted to 1990 or later (Zhang et al., 2015; Nie et al., 2017; Wang et al., 2020; Chen et al., 2021). The inventory for the period since the 1960s has been focussed on specific regions and has provided minimal information in terms of area and number at larger spatial scales. For example, only 73 lakes from the Poiqu basin, in the central Himalaya (Zhang et al., 2019), and 50 moraine-dammed lakes across the transboundary area between Bhutan and TAR have been studied (Komori, 2008). Glacial lake changes since the 1960s have, therefore, remained largely unknown due to the limitations imposed by the lack of high-resolution imageries. This study provides the first comprehensive assessment of glacial lake change since the 1960s for the entire Bhutan Himalaya.

We observed that the increase in non-glacier-fed lakes was comparable to the rise in glacier-fed lakes within the last 6 decades. This finding contrasts with regional-scale studies where a small change in the trend of non-glacier-fed lakes observed since 1990 (Zhang et al., 2015) or a relatively stable situation (Nie et al., 2017). This difference could be due to the substantial difference in the temporal range between the previous studies (since 1990) and the current study (since 1960s). Here, following Ageta et al. (2000), we defined non-glacier-fed lakes as those lakes which have lost hydrological connection with their mother glacier. Based on this definition, the increase in the number and area of non-glacier-fed lakes was primarily due to the transition of glacier-fed lakes to non-glacier-fed lakes rather than the expansion of existing lakes. For example, 157 glacier-fed-lakes evolved into non-glacier-fed lakes between the 1960s and 2016 as they lost hydrological connection with their glacier due to rapid shrinkage of the glacier or complete disappearance of the glaciers. Further, 192 new glacier-unconnected lakes were formed over the period, corroborating the idea of rapid delinking of the glacial lakes from their mother glacier. We determined that the decrease in glacier-fed lake area in the elevation range 4,500–5,000 m a.s.l was also caused due to transition of glacier-fed to non-glacier-fed lakes. The retreat of the glaciers between the 1960s and 1993 was found to be more prominent in the Bhutan Himalaya (Karma et al., 2003). Other studies have also reported unprecedented glacier melt and mass loss in the BTH after the 1990s (Gardelle et al., 2011; Bajracharya et al., 2014; Tshering and Fujita, 2016) in agreement with our results regarding the disconnection of lakes from their glaciers over time. Our study has revealed that extending the study period back to the 1960s, which marks the formation stage of most of the glacial lakes in the Bhutan Himalaya (Komori, 2008), provides new insights on the transition from glacier-fed to non-glacier-fed lakes.

Greater areal expansion was observed for glacier-connected lakes than for glacier-unconnected lakes, which resonates with earlier findings for the entire Himalayan region (Nie et al., 2017; Wang et al., 2020; Zhang et al., 2015). Such differences between connected and unconnected lakes could be explained by considering the sources of water for the lakes. When a glacial lake is no longer in contact with its glacier, the only water source is drainage of meltwater from the upstream glacier and precipitation. In contrast, a lake in contact with its glacier receives additional water from subaerial melting, water-line melting, and calving of the glacier front into the lake (Gardelle et al., 2011). Surface area expansion, through a positive feedback mechanism between the lake and ice, is favored for glacial lakes connected to glaciers (King et al., 2018; King et al., 2019). On the other hand, the decreasing trend of precipitation between the 1970s and 2010 across the Himalayas in general (Yao et al., 2012), and in Bhutan in particular (Supplementary Figure S2), was observed to be an essential contributing factor to the growth of non-glacier-connected lakes. Other possible factors contributing to the difference in growth rates are water loss through evaporation, infiltration, and discharge through drainage channels (Xu and Feng, 1994). The importance of such factors is supported by our findings that glacial lakes close to a glacier exhibit expansion, while lakes further away have exhibited recession. This finding supports the idea that glacier-melt water plays a major role in expanding lakes, as claimed in earlier studies (Zhang et al., 2015).

In contrast to earlier findings (Zhang et al., 2015; Nie et al., 2017), we observed that the area and number of supraglacial lakes has decreased over the last 6 decades. This contrasting behavior of supraglacial lakes could also be due to ∼60 years of evolution, as the highly transient nature of supraglacial lakes (Benn et al., 2001; Röhl, 2008; Gardelle et al., 2011; Taylor et al., 2021) limits their long-term expansion. Komori (2008) has suggested that most of the proglacial lakes in the Bhutan Himalaya would have started to form in the 1950s−1970s, and has also proposed three stages of typical proglacial lake formation: 1) appearance of the supraglacial lakes, 2) merging of supraglacial lakes to a single coalesced lake, and 3) stable expansion of coalesced lakes. If this is correct, it is possible that most of the present-day proglacial lakes were in the supraglacial lake stage in the 1960s and have evolved gradually since. For example, we have observed that the typical present-day proglacial lakes, Lugge Tsho, Rapstreng Tsho, and Chubda Tsho, were in the proglacial lake stage in the 1960s (Figure 10).

Across the different size classes, the number increase over the last 6 decades was found to be highest amongst the smaller glacial lakes, a similar trend to that found in earlier studies (Zhang et al., 2015; Khadka et al., 2018). We observed that total area expansion was almost the same across all the size classes, which is not in agreement with the trend in other regions such as the Nepal Himalaya (Khadka et al., 2018) and the Third Pole (Zhang et al., 2015). Most of the large lakes (0.5–1 and 1–5 km²) in the BTH were glacier-connected lakes, which constitute 32.8% of the total area despite their small numbers. The rapid expansion of those glacier-connected lakes has significantly influenced the overall expansion statistics across the different size categories.

Glacial lake expansion, in general, is in line with the increased atmospheric warming and eventual negative glacier mass balance. For example, the pronounced temperature rise in the higher elevations across the TAR and surrounding regions (Liu et al., 2009) corresponds well to a negative trend in mean glacier mass balance from the 1960s to 2000 (Bolch et al., 2012) and greater expansion of glacial lakes from 1990 to 2020 in the Himalayan regions (Zhang et al., 2015). Several past studies have revealed the rapid retreat of glaciers in the Bhutan Himalaya (Karma et al., 2003; Gardelle et al., 2013; Bajracharya et al., 2014; Tshering and Fujita, 2016), an accelerated glacier mass loss which is consistent with the current finding of glacial lake expansion since the 1960s. On the other hand, a decrease in precipitation (summer snow) accelerates glacier melting by lowering surface albedo (Fujita, 2008). Analysis of temperature and precipitation data (1997–2017) collected from eight weather stations located above 2,500 m. a.s.l (NCHM, 2018) indicates overall increasing temperature and decreasing precipitation trends, which fit with the current findings of glacial lake expansion (Supplementary Figure S2). However, summer and winter precipitation were not differentiated in this analysis.

The extent of glacial lake expansion varied from one basin to another. Such differential glacial lake evolution across the different basins is possibly due to local variability in glacier mass loss. Across the HKH, glacial lakes have grown in the region where glaciers have lost mass (the Himalayas), and shrunk in the part where glaciers have lost limited mass or surged (Karakorum) (Gardelle et al., 2011). The most notable observation was the drastic decrease in size of the glacier-fed area in the Drangmechu basin (−49.817%) which is due to 49 lakes in the Drangmechu basin evolving from glacier-fed to non-glacier-fed lakes. However, the expansion of glacial lakes is also influenced by the local climate regime, such as a change in precipitation (Khadka et al., 2018). Future studies on glacier-lake interaction against a backdrop of local climate fluctuation in the Bhutan Himalaya may be able to provide definite insights into whether such an interaction is occurring there.

Using Landsat satellite data, previous studies have examined glacial lake change over relatively long temporal ranges. For example, decadal change (1990–2000−2010) (Zhang et al., 2015), a 28 year period (1990–2018) (Wang et al., 2020), 25 years (1990–2015) (Nie et al., 2017), and 20 years (1990–2009) (Gardelle et al., 2011). Glacial lake change behavior at shorter temporal resolution is relatively unknown, but knowledge of the interannual variation of glacial lakes is essential for understanding glacial lake evolution against the background ever-increasing GLOF hazard (Zheng et al., 2021; Veh et al., 2020). Chen et al. (2021) have generated annual glacial lake data sets for High Mountain Asia from 2008 to 2017 from Landsat data. However, for the Bhutan Himalaya, their data sets have high levels of uncertainty for determining annual glacial lake behavior, as well as underrepresenting the overall statistics. For example, in their data, the total area increased from 102.38 to 111.5 km² between 2014 and 2015, while falling sharply back to 102.78 km² between 2015 and 2016. In contrast, we provide annual glacial lake data sets from 2016 to 2020 in the Bhutan Himalaya using high resolution (10 m) Sentinel-2 and Sentinel-1 using a highly automated method coupled with careful manual inspection.

Overall, the annual mapping of glacial lakes from 2016 to 2020 indicates that glacial lakes in the Bhutan Himalaya have expanded at a rate of 0.96 km² per year. Considering the uncertainty of ±5% in the total area, the annual expansion rate was within the uncertainty range. Comparing the interannual variation between glacier-fed and nonglacier-fed lakes, we found that only the former expanded while the latter has remained constant. Further, the rate of annual expansion was more prominent with the supraglacial lakes (8.94%). Typical proglacial lakes usually expand by progressing the upstream end of the shoreline (Komori, 2008). We found 102 glacier-connected lakes showing interannual upper-shoreline progression, a result which also indicate the interannual expansion from 2016 to 2020 (examples in Figure 10). Going by these observations, our results indicate that there was an interannual expansion of glacier-connected lakes although the increase rate for all lakes were within the uncertainty range.

Identification of the hazard potential of glacial lakes is the first step towards GLOF mitigation and early warning preparation programs (Ives, Shrestha, and Mool, 2010). Here, we employed seven hazard factors and multiple-criteria decision-making tools (AHP) (Saaty, 2008) to produce a GLOF hazard-potential ranking for 278 glacial lakes with area ≥0.05 km² and with a glacier-lake distance of ≤1,000 m. Based on the glacial lake data sets of 1950–1999, ICIMOD has identified 24 potentially dangerous glacial lakes in the Bhutan Himalaya. However, its evaluation is based on qualitative judgment, and is highly subjective (Emmer and Vilímek, 2013). Moreover, it was based on old inventory data sets derived from topographic maps, which suffer from severe non-systematic shifts (Nagai et al., 2017). A regional study covering the entire Third Pole has identified 85 lakes in the Bhutan Himalaya as highly hazardous and 115 as very-highly hazardous, (Zheng et al., 2021). The numbers are higher than the findings of the current study. The discrepancy can be mainly put down to differences in the hazard-factor consideration. In contrast to the current study, they have integrated ice/snow avalanche and landslides/rockfall into a single “topographic potential” factor, which has led to an overestimation of the GLOF hazard level of the lakes. For example, we found only seven lakes out of 278 with zero topographic potential but 85 with zero avalanche area in the Bhutan Himalaya. While almost all glacial lakes pose flood hazards to downstream communities, it is essential that the hazard ranking prioritizes limited resources for flood mitigation activities (Ives et al., 2010). In this light, our result presents a robust illustration of the GLOF hazard potential of the glacial lakes in the Bhutan Himalaya. Supplementary Table S5 provides a comparison between current and previous GLOF hazard-assessment studies covering the Bhutan Himalayas.

The multi-criteria decision-based hazard assessment ranked Thorthormi Tsho as the most hazardous lake in the Bhutan Himalaya. Thorthomi Tsho is one of the most well-studied and monitored lakes due to its high flood potential. The lake has been observed since 2005 after being reported as highly dangerous by Häusler et al. (2000). Being aware of the grave danger it poses to downstream settlements and particularly to the Punakha Dzong and two large, nearly completed, hydropower power plants (Punatsangchu I and II), the Royal Government of Bhutan, in partnership with UNDP, started in 2008 to manually lower the lake level. Moreover, one of the subsidiary lakes (lake II) of the Thorthormi Tsho produced a GLOF in 2019, causing alarm in the downstream settlement, although no causalities or economic loss were reported. A recent field verification by the National Centre for Hydrology and Meteorology of Bhutan has further reemphasized that this lake is in a perilous state with unstable moraine and an actively receding mother glacier, and may drain into Raphstreng Tsho causing a cascading flood (NCHM, 2019b). Here, Raphstreng Tsho, with an area of 1.344 km^2, is also given a high flood potential with a hazard score of 0.72. The GLOF from Lugge Tsho was the most devastating flood recorded so far in Bhutan. After the GLOF of 1994, Lugge Tsho expanded rapidly to the current area of 1.39 km² (determined in this inventory) and a volume of 65.19 million m³ (NCHM, 2019c). It is also considered one of the most dangerous proglacial lakes in the BTH (NCHM, 2019c). Our hazard assessment gives Lugge Tsho a very high flood potential with a hazard score of 0.83. The Tarina Tsho I, Tarina II, and Chubdha Tsho, all previously recognized as the PDGLs (NCHM, 2019c), are all designated here as high to very high GLOF hazard lakes. These findings provide adequate evidence that the GLOF hazard assessment method adopted here was robust enough.

However, two previously known PDGLs in the Mangdechu basin, Metatshota, and Mang-gl-385 fell to medium and very-low GLOF hazard in this study. Metatshota has no SLA and was not actively expanding. However, it is a proglacial lake with a large avalanche zone, which is considered the most crucial GLOF hazard factor in the Himalayas. Moreover, a previous study has revealed that a flood from the Metatshota is likely to impact some farmlands in the downstream settlement of Bjizam and Tingtingbi (Koike and Takenaka, 2012). We, therefore, do not wholly rule out its GLOF hazard in the future. ICIMOD considered Mang-gl-385 to be a PDGL as it is an ice-dammed lake. However, our assessment revealed that Mang-gl-385 has no avalanche, landslide, or SLA areas, which are the factors with the most influence in determining a lake’s GLOF hazard.

Hydro-power plants have been the backbone of Bhutan’s economy, playing a vital role in enhancing the country’s GDP and overall living standard. Being dependent on the glacier-fed rivers (MOEA, 2021), the hydropower plants are inevitably vulnerable to any potential future GLOF. The GLOF may directly damage the plant as in Rishiganga, and Dhaulinganga, Uttarakhand, India (Shugar et al., 2021) or disturb the sustainable water flow in the long term (Xianbao et al., 2021; Nie et al., 2021). The National Commission for Hydrology and Meteorology, the agency that looks after the GLOF threat in Bhutan, conducts annual monitoring of the glacial lakes and, based on the ICIMOD hazard assessment report of 2001, it has identified 17 PDGLs for regular monitoring. As mentioned above, the ICIMOD assessment has drawbacks which do not render it ideal for considering the current GLOF hazard status of the glacial lakes. We have determined a GLOF hazard ranking of glacial lakes through well-established methods and using the latest updated inventory and our results should form the basis for reconsidering which PDGLs require monitoring in the future.

We used multiple factors and a well-established semi-quantitative decision-making approach, AHP, to generate the GLOF hazard ranking of the Bhutanese glacial lakes. Our study identified the previously well-known PDGLs as high GLOF hazard lakes, underlining the reliability of our method and results. While adding hazard factors does not necessarily improve the efficiency of the GLOF evaluation methods, we acknowledge that in future studies, integrating ground temperature and permafrost modelling, may provide even more robust results (Allen et al., 2019). Alpine regions elsewhere have exhibited increasing ground temperature and thawing of permafrost leading to slope failure (Haberkorn et al., 2021; Swanson, 2021), which is the key cause of dam overtopping of glacial lakes. However, our recommendations are made within the scope of freely available remote sensing imagery, although we acknowledge that novel field-based data can provide more appropriate information. One potential source of uncertainty in our results may stem from the glacier data sets used in this study. We used the glacier data set of Nagai et al. (2016), which is based purely on the Bhutan Himalaya and has a higher accuracy than other data sets such as RGI v6.0 (RGI Consortium, 2017) and GAMDAM (Sakai, 2018). However, the data set was produced using ALOS imagery from 2006 to 2011, which is almost a decade old, and could, therefore, have led to avalanche areas being slightly overestimated. Since we used multiple criteria, such a minor variation should not have affected our results significantly. We are also mindful of the uncertainty of the other remote sensing data used, such as SRTM-30. Thus, we reemphasize that our research is mainly concerned with first-order hazard assessment, and should be complemented with further detailed in-situ studies and high-resolution satellite observations.