Evaluation of Glacial Lake Outburst Flood Susceptibility Using Multi-Criteria Assessment Framework in Mahalangur Himalaya

Ongoing recession of glaciers in the Himalaya in response to global climate change has far-reaching impacts on the formation and expansion of glacial lakes. The subsequent glacial lake outburst floods (GLOFs) are a significant threat to lives and livelihoods as they can cause catastrophic damage up to hundreds of kilometres downstream. Previous studies have reported the rapid expansion of glacial lakes and several notable destructive past GLOF events in the Mahalangur Himalaya, suggesting a necessity of timely and updated GLOF susceptibility assessment. Here, an updated inventory of glacial lakes across the Mahalangur Himalaya is developed based on 10-m Sentinel-2 satellite data from 2018. Additionally, the GLOF susceptibilities of glacial lakes (≥0.045 km2) are evaluated using a multi-criteria-based assessment framework where six key factors are selected and analyzed. Weight for each factor was assigned from the analytical hierarchy process. Glacial lakes are classified into very low, low, medium, high, and very high GLOF susceptibility classes depending upon their susceptibility index based on analysis of three historical GLOF events in the study area. The result shows the existence of 345 glacial lakes (>0.001 km2) with a total area of 18.80 ± 1.35 km2 across the region in 2018. Furthermore, out of the 64 glacial lakes (≥0.045 km2) assessed, seven were identified with very high GLOF susceptibility. We underline that pronounced glacier-lake interaction will likely increase the GLOF susceptibility. Regular monitoring and more detailed fieldworks for these highly susceptible glacial lakes are necessary. This will benefit in early warning and disaster risk reduction of downstream communities.


INTRODUCTION
The Himalaya is characterized by a dense distribution of glacial lakes (Zhang et al., 2015;Maharjan et al., 2018), and they are rapidly expanding due to glaciers wasting in response to global climate change (Zhang et al., 2015;Nie et al., 2017). Such glacial lakes are located on the surface of glaciers (supraglacial lakes or ponds), behind the terminal or lateral moraines (proglacial lakes), and other lakes which are disconnected, but lie in the periphery to the glacier and are fed by the glacier ice and snowmelt (Salerno et al., 2012;Otto, 2019). These glacial lakes display characteristic differences in terms of formation, dam structure, lifespan, expansion, emergence, disappearance and impact upon burst (Richardson and Reynolds, 2000;Mertes et al., 2017;Nie et al., 2017). Glacial lake outburst flood (GLOF) occur when the sudden release of large water volumes from the glacial lakes appear either by dam failure or overtopping waves due to external triggering or self-destruction, such as ice/rock avalanche and extreme precipitation into the lake (Somos-Valenzuela et al., 2015;Harrison et al., 2018). GLOFs are often associated with sediments and debris flows, having the potential for catastrophic damage to downstream areas (Osti and Egashira, 2009;Liu et al., 2014). The Himalayan region has been documented as a hot spot of GLOFs primarily from morainedammed glacial lakes with several historical GLOF records (Carrivick and Tweed, 2016;Nie et al., 2018;Veh et al., 2019).
The formation of new glacial lakes and the rapid expansion of the existing ones increases the risk of GLOFs in the Himalaya Nie et al., 2017;Khadka et al., 2018). Mahalangur Himalaya, a section of central Himalaya, is a highly glacierized area (Bajracharya and Mool, 2009). Glaciers in this region are responding to rising air temperature and reduced annual precipitation with moderate negative mass balance (Salerno et al., 2015). Consequently, an accelerated areal expansion of glacial lakes has been observed between 1987 and 2017 (Khadka et al., 2018) and lakes such as Imja Tsho and Barun Tsho have the possibility of future expansion (Watson et al., 2020). This implies that the potential threat posed by glacial lakes in this region may also be increasing, and thus a timely and updated assessment of GLOF susceptibility, the tendency of future occurrence of GLOF hazards, is urgent. In addition, about six historical GLOF events have been wellrecorded in the Mahalangur region (Byers et al., 2018;Nie et al., 2018;Veh et al., 2019); for instance, notable GLOF events from Nare Lake (1977), Dig Tsho (1985), Tam Pokhari (1998 and Langmale Lake (2017) (ICIMOD, 2011;Byers et al., 2018). The avalanche triggered displacement waves in Dig Tsho resulted in an estimated peak discharge of >2,000 m 3 s −1 , destroying a small hydropower project, a trail, and causing the loss of life and property (Vuichard and Zimmermann, 1987;Richardson and Reynolds, 2000). Almost more than a decade later, the moraine dam of Tam Pokhari collapsed in 1998 due to avalanche/landslides (Byers et al., 2013;Lamsal et al., 2016b) with an estimated 10,000 m 3 s −1 peak discharge and approximately 18 million m 3 water being released (Osti and Egashira, 2009). This GLOF caused nearly US$ 2 million in economic damage and the loss of lives (Dwivedi et al., 2000;ICIMOD, 2011). A recent GLOF in 2017 from the Langmale Lake in the Barun Valley was triggered by rockfall and had a peak discharge of 4,400 ± 1,800 m 3 s −1 , damaging buildings, bridges, and farmlands (Byers et al., 2018). The triggering mechanisms such as ice, snow and rock avalanching, glacier calving, wave overtopping, atmospheric triggers (extreme precipitation), ice-cored moraine degradation and earthquakes are the causes of dam failure and thus GLOFs, which are considered as the components in the GLOF modeling cascade (Westoby et al., 2014).
Different methods are used to appraise the potentially dangerous glacial lakes (PDGL) across the Himalaya (e.g., Rounce et al. (2016); Worni et al. (2013); Prakash and Nagarajan (2017)), which were adopted and modified from the Swiss Alps (Huggel et al., 2004), British Columbia (McKillop andClague, 2007), and other parts of the world (Wang et al., 2011;Emmer and Vilímek, 2013). These methods are generally different among others, with the selection of contrasting parameters for hazard and risk assessment. Rounce et al. (2016) comprehensively reviewed the existing methods and developed a robust and unified evidence-based GLOF hazard assessment framework. Based on this method, GLOF hazard assessment of glacial lakes (≥0.1 km 2 ) in Nepalese Himalaya has been conducted and classified into four levels, where glacial lakes with steep moraine dam and susceptibility to avalanche were considered to be very high hazard level (Rounce et al., 2017b). Allen et al. (2019) assessed the GLOF risk of all glacial lakes (≥0.1 km 2 ) in the Tibetan Plateau as a function of hazard and exposure using an automated, large-scale assessment approach that integrated lake area, upstream watershed area of each lake, the topographic potential for ice/rock avalanches and mean slope of the dam to model and ranks the hazard level of each lake. They identified the exposure of buildings to the modelled GLOF flow path for each lake to represent its subsequent downstream impacts. Wang et al. (2015) identified the PDGLs in the Chinese Himalaya by selecting four factors and using the analytical hierarchy process (AHP) to determine their weights. Moreover, eleven factors were considered by Prakash and Nagarajan (2017) for identifying the PDGLs over the western Himalaya whose precise remote measurements are difficult. Hence, the selection and use of an appropriate method and factors for GLOF susceptibility are necessary to identify the PDGLs and provide the first-order ranking of glacial lakes. Furthermore, the recent destructive GLOF events in the Nepalese (Byers et al., 2018) and Indian Himalaya (Das et al., 2015;Rafiq et al., 2019) are reported from glacial lakes <0.1 km 2 . Thus, we aim at lowering the size threshold to 0.045 km 2 from 0.1 km 2 , which was mostly used by previous studies, to include more small glacial lakes in this study. The lowering of the size threshold has also been emphasized by previous studies (Nie et al., 2018;Veh et al., 2019).
In this study, we aim to 1) map the glacial lakes across the Mahalangur Himalaya using the medium to high-resolution Landsat, and Sentinel-2 satellite images, and 2) evaluate their GLOF susceptibility selecting influencing factors and assigning their weights using the AHP method. More importantly, Nepal is ranked as one of the countries in the Himalaya with the greatest national-level socio-economic impacts from the GLOFs (Carrivick and Tweed, 2016). This study will evaluate the GLOF susceptibility of the lakes with the aim of identifying the potentially dangerous glacial lakes for detail investigations and for safeguarding life and properties downstream through the informed-decision making.

STUDY AREA
Mahalangur Himalaya (Mahalangur Himal) is the section of Himalaya covering northeastern Nepal and south-central Tibet (China) bounded on the east by the Popti La and the Arun River and on the west by the Gyabrag Glacier, the Nangpa La, Nangpa Tsangpo and Dudh Koshi (Carter, 1985). This study is focused on the southern part of the main Himalaya range encompassing the Dudh Koshi and Barun River basins, covering a total area of 3,701 km 2 ( Figure 1). This region includes four of the six highest mountain peaks of the world: Mt. Everest (8,849 m a.s.l., Sagarmatha in Nepali), Mt. Lhotse (8,516 m a.s.l.), Mt. Makalu (8,485 m a.s.l.), and Mt. Cho You (8,188 m a.s.l.) (Figure 1). The climate of the study area ranges from the temperate in the lower elevation regions in the south to polar tundra climate in the higher elevation regions in the north (Karki et al., 2016). The climate is influenced by South Asian monsoons in summer and westerlies in winter (Hamal et al., 2020;Sharma et al., 2020b). The climatology of this region obtained from the Pyramid weather station (5,050 m a.s.l., Figure 1A) reveals the maximum precipitation and air temperature during the summer (June-September) and minimum during the winter (December-February) ( Figure 1C). The temporal analysis of maximum, minimum and mean monthly air temperature shows +1.75, -5.71, and -2.48°C month −1 during 1994-2013, while cumulated mean precipitation of ∼436 mm year −1 during 1994-2016 ( Figure 1C).
The study area covers a total glacierized area of ∼740 km 2 derived from an updated GAMDAM glacier inventory (Sakai, 2019). Most of the glaciers are valley glaciers with debris cover in the ablation zone that extends in broad elevation range from 4,817 to 6,732 m a.s.l. (Bajracharya and Mool, 2009;Rounce and McKinney, 2014;Salerno et al., 2015). The mean mass balance for nine large glaciers in the Dudh Koshi basin was found to be −0.58 ± 0.19 m w.e. a −1 between 2005 and 2015 (King et al., 2017). The glaciers in the southern slope of the Everest region were found to be shrinking with an overall area loss of 13 ± 3% during 1962-2011 (Thakuri et al., 2014). This region hosts the highest number of glacial lakes per square kilometer compared to other parts of Nepal (Khadka et al., 2018) and has previously experienced six GLOF events in 1977, 1985, 1998(Byers et al., 2018Nie et al., 2018).

Glacial Lake Delineation and Classification
The glacial lake extent was delineated by consistent manual digitization from the geometrically corrected false-colour composite Landsat (30 m) and Sentinel-2 images (10 m). The alternative images of the same season were used whenever the shading and cloud cover disturbed the lake delineation. Moreover, the scene classification algorithm in Sentinel Application Platform (SNAP) (Chand and Watanabe, 2019) and Google Earth were used to remove the effect of shading in lake delineations. All lakes within 10 km buffer zone from the nearest glacier were included and classified into four types based on their positions in relation to the glaciers: supraglacial lakes on the surface of glaciers; proglacial lakes at the lateral or terminal position of the glacier; and unconnected lakes far from glaciers that were further classified into unconnected glacial-fed and nonglacial fed lakes depending upon whether the lake gets fed by glacier melt Khadka et al., 2018). For glacial lake classification, each lake was assigned a unique code in ArcGIS and the lake, glacier, and lake watershed shapefiles were converted into KML files. The KML files were loaded into Google Earth and lakes were classified into different types (supraglacial, proglacial and unconnected lakes) according to visual interpretation using their respective codes (Khadka et al., 2018). In Google Earth, lakes usually far from glaciers with no glaciers in their watershed were classified as unconnected non-glacial fed (Zhang et al., 2015). The minimum size threshold of ∼0.001 km 2 (10 pixels of Sentinel-2 image) was used to map glacial lakes in 2018.

GLOF Susceptibility
We designed and applied the multi-criteria decision-based method to assess the GLOF susceptibility across the Mahalangur Himalaya. The overall processes followed in the GLOF susceptibility assessment is shown in Figure 2. Here, we identified and compiled several factors governing the mechanism of GLOF potential based on previous literature (Supplementary Table  S1). Similarly, we complied and studied the details characteristics of three GLOFs in the study region ( Table 2). Various triggering factors can initiate the GLOFs through either by dam overtopping or dam failure. The GLOF initiating from these two cases are determined by dam characteristics (dam type, dam freeboard, dam crest width, dam width to height ratio, slope of the dam, piping), lake and glacier characteristics (lake area, lake perimeter, area expansion rate, glacier and lake distance), lake surrounding characteristics (potential for snow/rockfall avalanches), and other characteristics (earthquake, extreme precipitation events) (Emmer and Vilímek, 2014;Prakash and Nagarajan, 2017). Here, six factors ((lake area (F1), expansion rate (F2), glacier-lake distance (F3), dam front slope gradient (F4), ice/snow avalanche potential (F5) and potential of rockfall with upstream GLOF (F6)) were selected to examine the GLOF susceptibility as they represent lake and glacier characteristics, dam characteristics and surrounding terrain characteristics (Table 3 and Supplementary Table S1).
Eminently, factors F1-F6 are used to characterize the GLOF resulting from dam failure, whereas F1, F2, F5, and F6 are used for the GLOF resulting from dam overtopping (Prakash and Nagarajan, 2017). Furthermore, these factors well represent the characteristics of three well-documented past GLOF events within the study area ( Table 2) and can be remotely measured from freely available datasets. Each factor was classified into alternatives based on previous literature, the potentiality for triggering a GLOF and their potential role on GLOF impact depending upon detail characteristics of three GLOF events in the region (Table 3), nevertheless, partly influenced by subjective expert judgement (refer Factor selection and GLOF susceptibility).
In this study, glacial lakes with an area of ≥0.045 km 2 were considered in the GLOF susceptibility assessment, as suggested by Nie et al. (2018). Moreover, GLOF from Dig Tsho and Tam Pokhari glacial lakes (∼0.5 km 2 ) had caused catastrophic damage in the region; meanwhile, large lakes have the potential to generate large flood volume than smaller lakes (ICIMOD, 2011;Fujita et al., 2013). Furthermore, lakes with the large surface area are exposed to higher chances of impact from the potential mass movement from the surrounding terrain. Hence, the size threshold of ≥0.5 km 2 was set as the upper limit, and the F1 class was adjusted accordingly ( Table 3).
Previous studies had used the threshold of >100% expansion to determine the lake growth in 40 years period for GLOF susceptibility assessment (Bolch et al., 2011;Aggarwal et al., 2017). Following the previous study, here, the threshold value of >50% was chosen for determining the expansion rate of glacial lakes in the period of 20 years . Notably, glacial lakes were considered as rapidly expanding that had an expansion rate of 20% in 20 year period (Nie et al., 2013;Nie et al., 2017). Thus, the expansion rate >50% was considered as the upper limit and <20% as the lower limit to adjust the F2 class (Table 3). Distinctly, out of three, two glacial lakes had high expansion rates before their burst in the study area ( Figure 3). Seven glacial lakes evolved after 1998, so their expansion rate was considered as the highest alternative.
As all glacial lakes were chosen regardless of types (lakes connected or unconnected to glaciers), the glacier and lake proximity was considered as one of the factors in the GLOF susceptibility evaluation, which was determined from delineated lake polygon and glacier outline. Out of 12 Himalayan GLOF cases, nine lakes were in contact with glacier before burst, while others were far within 700 m ( Supplementary Table S2). Similarly, the previously drained glacial lakes in Tibetan Plateau were reported to be at the distance of 0-800 m from the glacier, and Wang et al. (2011) used <500 m threshold determining it to be more favourable value for a lake to outburst. Moreover, those glacial lakes that have experienced GLOFs in the study area were connected to their parent glaciers ( Table 2). Thus, glacial lakes in contact with glaciers were set as the upper limit and >500 m as the lower limit to adjust the F3 class (Table 3).
Dam characteristics, such as dam geometry (dam width to height ratio, width of crest, dam distal face slope), dam material properties, ice-cored moraine and freeboard conditions govern the stability of the dam (Huggel et al., 2004;Wang et al., 2011;Prakash and Nagarajan, 2017). However, owing to the difficulty for precise remote measurements of these parameters, we used a single feature, i.e., dam front slope gradient (DFSG) to model the stability of the moraine dam. DFSG was defined as a tangent angle determined by the difference in height between the lake dam and the dam base, overall considering a horizontal distance downstream of ∼0.6-1 km as the dam base differently to all lakes ( Figure 4). The DFSG was calculated along the terminal moraine and/or flow path of the glacial lake. The DFSG is reasonably similar to the average steep lakefront area slope (Fujita et al., 2013) and is mostly relevant to moraine-dammed glacial lakes, indicating possibilities of erosion, incision and failure. Thus, glacial lakes were classified according to their   Table S7) from visual interpretation in Google Earth following the classification scheme described by Otto (2019). Here, the glacial lakes with a steep dam (DFSG >10°) were classified as the highest alternative, and F4 has been adjusted accordingly ( Table 3). The value of DFSG for bedrock dammed glacial lakes was assigned as zero as they are relatively stable when compared to moraine dams. Moreover, the GLOF in the Himalaya are most often reported as breaches from moraine dam (ICIMOD, 2011;Nie et al., 2018). The impact of the mass movement has caused more than 70% of the known GLOF events (dynamic failure) in the Himalaya, either from ice/snow avalanche, rockfall/landslide or upstream GLOF (Supplementary Figure S1). The avalanche is the major trigger of GLOFs in the Himalaya (Nie et al., 2018). Thus, this factor was ranked as the highest priority by expert judgement, and the maximum weight assigned to this factor from AHP is justifiable (Table 3). Glacial lakes susceptible to avalanche and rockfall were classified as the highest alternative (F5 and F6), while those susceptible to upstream GLOF (lakes >0.05 km 2 located upstream) were classified as medium ( Table 3). Within the upstream watershed of each glacial lake, ice and rock avalanche-prone areas were calculated as the sum of the pixels fulfilling the criteria of slope >30°, where the overall trajectory slope is >17°(or angle of reach; tanα 0.25) . Further, glacier outlines were used to discriminate the probable ice and rock/landslide areas (Figure 4). The overall process was modelled in ArcGIS 10.5 software. Moreover, those lakes with large watersheds where the prone areas were located at horizontal distance >2 km 2 were discarded based on flow direction algorithm in ArcGIS 10.5, assuming that the avalanche/rockfall does not hit the lake directly. However, in exceptional cases, mass movements from more gentle slopes and obtaining larger run-out distances are possible.

Determination of Weights and GLOF Susceptibility Index
The AHP method was used to assign the weights of each factor. The AHP is a semi-quantitative decision-making approach using weights through pairwise comparison between different factors without inconsistencies (Saaty, 1990). The expert evaluation was conducted to find the relative importance of each factor against each other by establishing a comparison matrix ( Supplementary  Table S3). Finally, the comparison matrix was normalized to determine the rank and weight of each factor through AHP (Tables 3 and Supplementary Table S4). The consistency ratio of the computed AHP was ∼9% within an acceptable range for good consistency in the judgements and pairwise comparisons (Supplementary Table S5) (Saaty and Vargas, 2012). Three alternatives are provided with index values of 1, 0.5 and 0.25 (Prakash and Nagarajan, 2017). The final weights for each factor were computed by multiplying the factor weight (W i ) with the class index values (C i ) based on the measurement of lake factors. The final weights of each factor were summed to evaluate the outburst susceptibility of the lake (Eq. (1)).

GLOF susceptibility index
The value of GLOF susceptibility index falls between 0.25 and 1. To display the results of GLOF susceptibility, the values obtained from GLOF susceptibility index (0.25-1) are categorized into five classes using the equal interval classification method (Figure 2). This method is representative as revealed by previous three GLOFs in the study area, which had the susceptibility index exceeding 0.85 ( Table 2).

Sensitivity Analysis
The sources of uncertainties in the GLOF susceptibility assessment can arise from the satellite data used for preparing the inventory and the DEM used for estimating the factors. The uncertainty in glacial lake delineation was calculated as a product of the perimeter of the lake and half of the resolution of data (Salerno et al., 2012). The uncertainty in lake delineation affects the input values for F1 (lake area) and F2 (lake expansion rate); however previous studies Rounce et al., 2017b) have successfully used the similar dataset for preparing the inventory and computing the factors for the GLOF hazard (Rounce et al., 2017b;Allen et al., 2019), susceptibility (Prakash and Nagarajan, 2017), and risk  assessments. Moreover, the use of 30 m SRTM DEM in the high mountainous region has limitation to reveal the topography precisely. In a multi-criteria decision-based GLOF assessment, sensitivity analysis (SA) plays a key role to determine how much uncertainty of the results of the assessment are influenced by the uncertainty of its input criteria (Saltelli et al., 1999). Here, we performed the SA in two ways: 1) shifting the class/alternative and 2) modifying the GLOF susceptibility class. We checked the sensitivity of GLOF assessment by shifting the low and medium class/alternative of F1 to medium and high classes, respectively keeping other factors constant. Similarly, SA was performed for the remaining four factors (F2, F3, F5, and F6) except F4. For F4, the SA was conducted by setting the high and low alternatives to >5°and <5°, respectively, keeping other factors constant. The second SA examines the extent to which GLOF susceptibility of lakes will change when the susceptibility class ( Figure 2) is modified by ±5%, i.e., ±0.05.

Inventory of Glacial Lakes
In 2018, a total of 345 high mountain lakes (>0.001 km 2 ), covering a total area of 18.80 ± 1.35 km 2 , are distributed across the Mahalangur Himalaya (Figure 1). Among the lakes, 114 (1.38 ± 0.27 km 2 ) are supraglacial, 113 (12.78 ± 0.68 km 2 ) are proglacial, 87 (3.74 ± 0.31 km 2 ) are unconnected glacier-fed, and 31 (0.89 ± 0.10 km 2 ) are unconnected non-glacial fed lakes ( Figure 5). The mean size of the glacial lakes is 0.05 km 2 , and the largest glacial lake in the region is the Lower Barun Tsho (1.92 ± 0.04 km 2 ). About 83% of the total glacial lakes are below the mean size (0.05 km 2 ); however, remaining lakes, i.e. those above the mean size constitute about 80% of the total lake area. The glacial lakes are distributed in elevation range between 3,772 and 5,592 m a.s.l. with mean elevation of 4,989 m a.s.l. The density of the glacial lakes in the study area is ∼10 glacial lake per 100 km 2 . The limnic ratio, the ratio of the total lake area to the total documented area where the inventory was performed, was 0.005 (Chou, 1995).

Expansion of Glacial Lakes
In this study, we selected 64 glacial lakes with an area ≥0.05 km 2 in 2018 for GLOF susceptibility assessment. Thus, their surface area changes between 1998 and 2018 were examined ( Figure 6). The results show that seven out of 64 lakes evolved after 1998, of which two were proglacial lakes, and the remaining five were supraglacial lakes. The total area of the glacial lakes in 1998 was 12.0 ± 1.6 km 2 , which expanded by ∼30%, reaching 15.7 ± 0.7 km 2 in 2018. The three typical end-moraine proglacial lakes (Imja Tsho, Lumding Tsho and Barun Tsho) connected with their parent glaciers expanded the most in the region, contributing to ∼65% of total expansion ( Figure 6, Supplementary Table S6). Besides, previous small supraglacial ponds on the surface of the glaciers grew and coalesced rapidly to form large glacial lakes that are usually dammed by the mix of ice and moraine ( Figure 6B). Such characteristics may involve in developing a large proglacial lake behind the terminal moraine (Watanabe et al., 1994). Notably, most of the supraglacial lakes are dynamic, transient and often drain through the expansion of englacial conduits during summer (Benn et al., 2001;Watson et al., 2016).

GLOF Susceptibility
The GLOF susceptibility assessment shows that out of 64 glacial lakes assessed, seven were classified with very high level of GLOF susceptibility (Table 4). Similarly, 11 glacial lakes were classified with high, 18 with medium, 24 with low and 4 with very low GLOF susceptibility (Figure 7 and Supplementary Table S7).
Glacial lakes within small size class have high level of GLOF susceptibility ( Figure 7B), implying that small size glacial lakes must be included in a comprehensive risk evaluation. Our method ranked notable glacial lakes like Imja Tsho as medium level and others (Lumding Tsho, Barun Tsho and South Chamlang Tsho) as very high outburst susceptibility. The details on these glacial lakes with their susceptibility classes are presented in Table 4 and Supplementary Table S7.

Glacial Lake Inventory
Glacial lakes are indicators of climate change in the alpine environment formed due to climate-induced glacier dynamics (Otto, 2019). There are various classification schemes to classify the glacial lakes relating to glacial and non-glacial processes, their location, damming materials and structures (ICIMOD, 2011;Maharjan et al., 2018). The conventional methods for identifying glacial lakes around 10 km buffer zone from the glaciers include non-glacial fed lakes (those lakes that do not receive glacial supply); thus, studies usually use the term "high mountain or alpine lakes" and classify the types of lakes (Zheng et al., 2019). The classification of lakes helps to distinguish between glacier-fed and non-glacier fed lakes and their damming types (Bhambri et al., 2018). Non-glacier fed lakes lie in the periphery of glacial environment and are mostly fed by seasonal snowmelt and rainfall (Zhang et al., 2015;Khadka et al., 2018), while other studies classify them as unconnected glacial lakes (Salerno et al., 2016). The automatic lake delineations in the Himalaya requires substantial post corrections ; thus, glacial lakes were consistently digitized from a single expert to explicitly produce a high-resolution inventory of glacial lakes in the study area. Comparing the delineation of 31 lakes between 10 m Sentinel-2 images and 0.5 m Pleiades imagery revealed that precise lake delineation from fine-resolution imagery yielded a slightly greater area (Supplementary Table S8). Moreover, two lakes joined by a narrow channel (strait) are visible in the fine resolution imagery featuring as a single lake, which is often misclassified as two lakes from medium to high-resolution images. The lake mapping from 10 m Sentinel-2 images yielded better results (Supplementary Table S8) with lower errors than mapped from 30 m Landsat images because the standard way of detecting an error in lake delineation (i.e., ±1/ 2 pixel × lake perimeter) depends upon the spatial resolution of images used (Salerno et al., 2012). A total of 103 additional lakes were mapped from high resolution Sentinel-2 images in the study area in 2018 when compared with the inventory in 2017 from Landsat imagery (Khadka et al., 2018). Furthermore, small supraglacial lakes are highly detected from Sentinel-2 imagery (Chand and Watanabe, 2019). However, the mapping depends upon the quality of image, threshold and method used. Here, the lake expansion rate was calculated utilizing different resolution images between 1998 (Landsat TM, 30 m) and 2018 (Sentinel-2, 10 m) instead of using same resolution images as the nominal  differences in the expansion rates did not change the class it falls under as the current method applies three classes for lake expansion rates (Table 3).

Factor Selection and GLOF Susceptibility
The Himalayas have been documented as a hotspot of GLOFs and more than 50 GLOFs are recorded, among them, few are reported with details (Nie et al., 2018). The factors used in the present study might not hold for all GLOFs in the Himalaya; instead, they are focused on the study area. Moreover, the details and exact triggering factors, lake area, expansion rate and glacier-lake activity before GLOF are still difficult to be obtained for all GLOF cases. Furthermore, Nie et al. (2018) heavily rely on old documents/reports/literature to report details about GLOF, and Veh et al. (2019) has reported some unreported GLOFs without details mechanism of lake failure. Thus, it is very difficult in the current study to analyze all the factors related to previous GLOFs in the Himalayan region to select key parameters and thresholds that govern the outburst. Notably, GLOF reported with details in the Himalaya had high susceptibility score (Supplementary Table S2). We recommend for improving the methodology in future studies in other regions as GLOF mechanism, impacts and risk differs from one region to other. For instance, destructive GLOF from lateral moraine dammed Chorabari lake (∼0.05 km 2 ) was triggered by extreme precipitation event (149 mm rainfall) in the Indian Himalaya which is very rare in our study area as precipitation analysis reveal only four extreme precipitation days (25-28 mm rainfall) during 1994-2016. In this study, we used a pre-defined lake area threshold of ≥0.045 km 2 to assess GLOF susceptibility. Thus, only 19% of the lakes were considered in the study due to the size threshold. We prioritized large lakes within the highest alternative because the mean area of the lakes that previously produced GLOF was 0.4 km 2 (Supplementary Table  S2). This does not imply that small lakes (<0.045 km 2 ) are not susceptible to burst. For example, 2016 GLOF from Gongbatongsha Tsho (0.017 km 2 ) was transboundary and Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 601288 destructive as it destroyed infrastructures and settlements along the SINO-Nepal highway (Liu et al., 2020). We used lake area rather than lake volume as a factor in determining the possible magnitude of a GLOF because most existing empirical formulas for volume calculation are areadependent (Fujita et al., 2013;Sakai, 2012). The expansion rate of the lake area is also essential because the growing lake area increases the potential flood volume. Besides, the growth of the proglacial lakes such as the Imja Tsho will bring the lakes closer to areas prone to ice/rock avalanche (Linsbauer et al., 2016;Watson et al., 2020), altering the GLOF susceptibility in the future. The expansion of South Chamlang Tsho is restricted by bedrock and lateral moraines, yet it is prone to avalanches (Figure 4) (Khadka et al., 2019;Lamsal et al., 2016a). Generally, the stability of a moraine dam depends upon its geometry, geomorphology and buried ice-core (Haritashya et al., 2018). The concept of using DFSG or steep lakefront (SLA) slope to analyze the stability of moraine dam is reasonable (Fujita et al., 2013) as piping failure due to seepages and crack on dam is difficult to be perceived through remote observations (Rounce et al., 2017a). Similarly, dam crest width and lake freeboard determine the magnitude of outburst resulting from dam overtopping; in this case, DFSG alone may not represent the dam conditions. These two factors were not included in the assessment since they are challenging to analyze from remote sensing, and most GLOFs in Nepal were due to dam failure. The atmospheric trigger (such as extreme rainfall) was not considered as it is inconvenient to determine for the lakes individually. Moreover, during 2014-2016, Sharma et al. (2020a) found that the Pyramid station only recorded 15 events >10 mm/ days (heavy precipitation days). Furthermore, it only recorded four R25 mm events (extreme precipitation days with >25 mm rainfall) during 1994-2016. However, atmospheric factors (heavy and extreme precipitation events) have the potentiality to trigger GLOF in two ways: directly by triggering a GLOF itself (Supplementary Figure S1) and indirectly by initiating sudden heavy snow and ice melt and triggering mass movements. This complexity is a limitation, even though the indirect trigger is captured by other factors, i.e. avalanche susceptibility used in the study. The size and topography of upstream watershed of each glacial lake determine the accumulation volume of precipitation or snowmelt that will alter the lake water level if it does not have outflow. Thus, the upstream watershed area was included in the assessment over Tibetan Plateau as a proxy to extreme precipitation events ; however, we did not consider this factor as extreme events are low in the study area. We excluded the seismic factor which was considered by a previous study (Prakash and Nagarajan, 2017) since it is difficult to predict accurately, although the earthquake can trigger a GLOF by initiating avalanches and weakening the moraine dam. It was reported that the 2015 Nepal earthquake and its aftershock had destabilized the glacial lakes by forming new cracks and triggering landslide activities within lateral and terminal moraines (Byers et al., 2017). Nevertheless, a subsequent study in the aftermath of a devastating earthquake did not find any evidence of GLOF from 491 glacial lakes (Kargel et al., 2016). Nepal lies in a high seismic zone with the same values of peak ground acceleration from the earthquake for the study area (Rahman and Bai, 2018); thus, it will not affect the final GLOF score for individual lakes.
We used a multi-criteria-based assessment framework to prioritize glacial lakes with high and very high GLOF susceptibility. This method is straightforward, systematic, and easy to implement and repeat. It is inevitable, of course, that it has certain limitations, such as subjective expert evaluation and judgement. Additionally, the weights of each factor obtained from AHP depends upon the two-way comparison matrix and is subjective (Supplementary Table S3). However, we have maintained its compatibility and consistency ratio (Supplementary Table S5) by adding a high sampling number (n 11) from the experts in the same field for the two-way comparison matrix (Supplementary Table S3). It should be emphasized that those glacial lakes classified as high and very high susceptibility do not necessarily indicate that they are currently in a state of outburst, but they require priority for regular monitoring and detailed ground investigation. The timing of a potential GLOF is challenging to estimate, however, some possible triggers such as snow/ice/rock avalanche and ice-cored moraine degradation will increase the risk of dam collapse (Westoby et al., 2014). The majority of recorded GLOFs in the Himalaya were produced by snow/ice/rock avalanches (Nie et al., 2018), and the probability of an avalanche or rockfall entering the lake and triggering the GLOF depends upon their volume, trajectory, velocity and positions in relation to the lake (Bolch et al., 2011;Rounce et al., 2016). For instance, the current risk of avalanches hitting into Imja Tsho is low due to the steep relief surrounding the lake is located far from the lake. In contrast, the steep hanging glaciers above Tam Pokhari glacial lake increases the susceptibility of the repeated outburst flood.
The GLOF susceptibility assessment is sensitive to change in thresholds of class/alternatives of factors. The results of the first sensitivity analysis revealed that 13-22% of the lakes (out of 64) changed their susceptibility class when the class/alternative level of F1-F3 are shifted. Similarly, 20% of lakes changed their class when high, and low class of F4 was altered to >5°and <5°, respectively keeping other factors constant. Moreover, 44% of the lakes changed their class when F5 was altered, revealing that ice/snow avalanche (F5) is the most sensitive parameter, among others. No change was observed in GLOF class of lakes when F6 was altered because almost 85% of lakes were prone to probable rockfall/landslides.
The results of the second sensitivity analysis show that if the susceptibility class ( Figure 2) is decreased by 0.05 then seven high GLOF class lakes (M30, M31, M34, M45, M61, M62, M63) will be classified into very high class. Similarly, if susceptibility class is increased by 0.05, four (M20, M47, M52, M53) out of seven very high lakes (Table 4) will fall into high class. This shows that the lakes are very sensitive to the small change in GLOF class.

Comparison With the Previous Assessments
ICIMOD (2011) identified ten potentially dangerous glacial lakes, ranking four as first, three of each as second and third-order critical lakes highly depending upon expert judgment based on Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 601288 physical and socio-economic criteria. In contrast, our assessment ranked glacial lakes based on its susceptibility score. The first order Lower Barun and Lumding Tsho are classified as very high while Imja Tsho as medium GLOF susceptibility. An assessment of GLOF hazard has been conducted across the Nepalese Himalaya to classify the hazard level of glacial lakes solely depending upon one or two parameters or a combination (Rounce et al., 2017b). For instance, glacial lakes prone to snow avalanche and with steep moraine dam were considered as very high hazardous. Glacial lakes with medium to high hazard level were considered either susceptible to avalanche/rockfall or steep moraine dam or ice-cored dam or upstream GLOF (Rounce et al., 2017b). Out of 131 glacial lakes (>0.1 km 2 ) across the Nepalese Himalaya, 19 glacial lakes were categorized with a very high hazard (Rounce et al., 2017b). Our assessment identified high GLOF susceptible lakes of size <0.1 km 2 ( Figure 7B), which was excluded in previous assessment due to size threshold. Both studies have identified three common glacial lakes with a very high level of hazard (Table 4). Our study differs from their study in terms of methodology, as our method considered the collective response of different factors to classify lakes.
We used updated GAMDAM glaciers outline, considering areas >30°to know ice avalanche areas assuming that the breaking of hanging glaciers (Margreth et al., 2017) in the Himalayas can also initiate GLOFs. The uncertainty associated with DEM and GAMDAM outline affected the mass movement calculations as glacier outlines are prepared from the Landsat images of 1990-2010 (Sakai, 2019). However, since our methodology applied slope >30°to find avalanche or rockfall areas and classified into two alternatives (Table 3), the quality of the glacier outline will have no qualitative and quantitative effects on the results. Here, F5 (ice/snow avalanche) played a dominant role in the final susceptibility score because of its highest factor weight. The GLOF susceptibility classification largely depends on factors selections, quality of the data used, classification scheme, weights and scores (Bolch et al., 2011) and the alternation or adjustment of them would change the susceptibility ranking.

Glacier-Lake Interaction Likely Increases the GLOF Susceptibility
Any variation in the morphology of the glacial lakes (e.g., growth) can change the lake area (F1) and its expansion rate (F2), subsequently, resulting to the shift in their classes in which they fall at present. Due to such changes, GLOF susceptibility classes can vary in the future. Coupled glacier-lake interactions play a significant role in the expansion of the proglacial and the dynamics of the supraglacial lakes, while the expansion of glacial lakes promotes the retreat of glaciers (Song et al., 2017;Zhang et al., 2019). Proglacial lakes with wide calving front have high expansion rates (Haritashya et al., 2018); however, the expansion rate varies between glaciers from high to steady and diminishing rates (King et al., 2018). In the study region, the future lake expansion modeling shows that notable calving glacial lakes like Imja Tsho and Lower Barun Tsho, when expanded at current rates, would expand volumetrically to store 90 × 10 6 and 62 × 106 m 3 of additional water by 2040 and 2060 than their estimate in 2018 (Watson et al., 2020). This expansion would increase the GLOF susceptibility in two ways: the expanding lakes may reach up to the areas prone to an avalanche, and volumetric expansion (area × depth) could increase the hydrostatic pressure to the moraine dam. Consequently, growing lake volume will increase the potential flood volume, which can cause more damage to downstream areas. In the case of Lumding Tsho, ranked as a very high GLOF susceptible level (Table 4), has the capability for future expansion when demonstrated from modelled glacier thickness (Farinotti et al., 2019) and slope derived from DEM ( Figure 8). The lake, when expanded at its current rate (0.025 km 2 yr −1 between 1998 and 2018) will attend the extent shown in Figure 8A by 2040s, increasing the water volume. Although its expansion rate is lower than Imja Tsho and Lower Barun Tsho (Supplementary Table S6), it needs to be monitored with a focus as given to Imja Tsho by several researchers. It is recommended that GLOF hazard assessment is needed to be improved by including the modeling of possible future lake extents and their volumetric analysis.

CONCLUSION
This study was conducted in the Mahalangur section of the Himalaya to investigate the glacial lakes, and GLOF susceptibility as this region includes the largest concentration of glacial lakes, high elevation glaciers, and several previously recorded GLOF events. A high-resolution inventory of glacial lakes from Sentinel-2 images shows 345 lakes with a total area of 18.80 ± 1.35 km 2 in 2018 across the study area. The multicriteria decision-based analytical hierarchy process was used to evaluate the degree of GLOF susceptibility from glacial lakes in the region discussing the several factors that regulate the probability of an outburst. With the eminent GLOF events from small size glacial lakes (<0.1 km 2 ), the size threshold was reduced to 0.05 km 2 , and a GLOF susceptibility assessment was conducted with the basis of six factors and detail characteristics of three past GLOF events in the study area.
The results show that out of 64 glacial lakes (≥0.05 km 2 ) in the region, seven were assessed to have very high outburst susceptibility while the remaining ranged from high to very low, providing more detail insights into GLOF susceptibility in the region. The very high GLOF susceptible lakes must be prioritized for detail investigation and studies. The methodology is replicable for basin level or regional level studies and can be improved with the detailed characteristics of past GLOF and associated parameters. The results are significant to the decision-makers and scientists as it has implications for integrated disaster risk reduction. This study emphasizes that the GLOF susceptibility increases as glacierlake interaction becomes pronounced, adding to the risk of GLOFs. Furthermore, the application of high-resolution data, with robust methodologies combined with in-situ studies, will help in improving GLOF susceptibility assessment.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: Landsat images were obtained from the U.S. Geological Survey (USGS) data portal (http://earthexplorer.usgs. gov) and Sentinel-2 images were downloaded from the Copernicus open access hub of the European Space Agency (https://scihub.copernicus.eu/).