Quantitative prioritization of potentially critical glacial Lakes in the Indus River basin using satellite derived parameters

Abstract Glacial Lake Outburst Floods (GLOFs) in a deglaciating environment, representing significant threat to downstream people and infrastructure, increasing the necessity of prioritization and assessment of potentially critical glacial lakes (PCGLs). In this current study, a comprehensive updated inventory of glacial lakes for the Indus River basin has been prepared, resulting in 5335 lakes (≥ 0.0025 ± 0.0006 km2), covering a total area of 173.95 ± 10.13 km2. Parameters have been derived using satellite data for two purposes, ‘preliminary screening’ and ‘quantitative assessment’. Using 4 parameters, 367 lakes were preliminary screened-in for further analysis, and 6 quantitative parameters were analysed using Equal Weights (EW) and Unequal Weights methods (UEW). These both methods were applied for two scenarios, scenario-1 which considers 73 lakes (≥ 0.1 km2), and scenario-2 which considers all 367 preliminary screened-in lakes (≥ 0.02 km2). UEW outperformed EW analysis and identified 20 and 48 high potential lakes for scenario-1 and scenario-2 respectively, where all 20 lakes of scenario-1 have been found in common with those of scenario-2, with a difference in their rank. But only 10 of them have storage volume ≥ 10 MCM and has been flagged as high-risk lakes, of which 6 are end-moraine dammed lakes and have been considered as PCGLs, due to their high bursting potential and self-destructive nature. Rest 4 high risk lakes along with remaining 38 lakes with volume < 10 MCM, are considered as critical lakes in case of dynamic mode of failure, caused due to mass movement of avalanche, landslide, heavy precipitation, etc. Highlights 152 historic GLOF events has been reported from 28 sites that occurred in Indus River basin Updated glacial lake inventory of Indus River basin using high resolution IRS satellite data Total of 5335 GL’s were mapped, and 48 of them are identified as high potential lakes UEW outperformed EW method to evaluate weighted integrated index for potentiality assessment 6 lakes are found to be as potentially critical glacial lakes (PCGLs) in the Indus River basin Graphical Abstract


Introduction
Third Pole is geomorphologically the largest and highest mountain region on Earth which includes the Tibetan Plateau, Himalayas, Karakoram, Hindu Kush, etc. , which also contains the world's largest number of glaciers and snow outside the Polar Regions (Zheng et al. 2021). Many studies undertaken globally showed that, with few exceptions, glaciers around the world have been thinning and retreating at unprecedented rates since the industrial revolution, which began around eighteenth century (Williams and Ferrigno 2012;NSIDC 2020), and also due to climate change, accelerating glacier melt which results in expansion of existing lakes and formation of new lakes in mountainous regions (Ageta et al. 2000;Benn et al. 2000;Reynolds 2000), are generally known as glacial lakes. Sudden discharge of water from a glacial lake triggered due to ice or rock avalanche falling/sliding, extreme climatic events, emergence of upstream lake, etc. (Worni et al. 2013;Rounce et al. 2017), can lead to breaching of the unstable debris dam resulting in flash floods with an exceptional erosion/transport potential in downstream areas (Cenderelli and Wohl 2001;Emmer and Vil ımek 2013), known as Glacial Lake Outburst Floods (GLOF). Given the risk posed to downstream communities, industry, and infrastructure in deglaciating mountains, there has been an increase in GLOFs research (Emmer and Vil ımek 2013), with many studies identifying dangerous glacial lakes, while many others estimating GLOF hazard and risk (Worni et al. 2013;Prakash and Nagarajan 2017;Dubey and Goyal 2020;Mal et al. 2021;Zheng et al. 2021).
But the key differences between these studies are their preliminary screening choice of large glacial lakes (> 0.05 km 2 or > 0.1 km 2 ) and subjective choice of parameters and corresponding weights assigned to them. According to Rounce et al. (2016) and Emmer et al. (2016), confusion among the stakeholders or decision/policy makers may get created due to such inconsistent hazard and risk classification done using various parameters and their weights. Hence, there are numerous systematic glacial lake inventories available and frequent studies on GLOFs have been carried out in different regions of the Himalayas, but very minimal or few relevant studies are available for the Western Himalayas (Prakash and Nagarajan 2017). Therefore, the objective of this study is to develop a novel framework to assess non-subjective quantitative prioritization of glacial lakes in part of Western Himalayas, accounting for primarily downstream colonised mechanisms to facilitate uniform conclusion. With the two main objectives, a) to derive remote sensing-based parameters of glacial lakes using high resolution satellite data in Indus River basin; b) to quantitatively prioritize potentially critical glacial lakes (PCGLs) using integrated index method for two scenarios.
We believe, this present investigation is unique and first to prioritize glacial lakes susceptible to GLOF in a quantitative manner accounting for two-scenarios which includes both small and large glacial lakes, accompanied with the determination of exposed and nearest downstream areas like bridges, buildings, croplands, and hydropower structure over the Indus River basin. Perhaps, limited and only few earlier studies in the Indus River basin were based only on large lakes; nonetheless, to characterise the multifarious impact of GLOF with respect to external triggering mechanisms, a more comprehensive framework that accounts for smaller lakes is required.

Study area
Indus River basin from its origin up to the Indian administrative boundary only is the present study area, covering a total of 3,42,841 km 2 area, extends geographically from 30.32 to 37.09 N and longitudes 72.50 to 82.45 E (as shown in Figure 1). Study area encompasses majorly within Western Himalayas and Karakoram, and some portion of Hindu Kush and Inner Tibet Mountain range (GTN-G 2017), primarily located in India and transboundary region. The Indus River originates from Bokhar Chu (glacier) in the northern slopes of Mt. Kailash (6714 m) in the Tibet Autonomous Region of China (Husain 2012;CWC 2019), and travels for a total of 2880 km before it merges into the Arabian sea, of which only 709 km lies within India. The river follows a long and straight course in Ladakh region (flowing in northwest direction in India), running between the Ladakh and the Zaskar ranges, and then turns southwest near Chillar in the Dardistan region. Major rivers flowing in the basin are Satluj, Beas, Ravi, Chenab, and Jhelum joining the main river on the left, whereas rivers like Shyok and Gilgit joining on the right. Topography in the study area varies from 222 m to 8564 m amsl, with mean elevation of 4114 m amsl. Whereas, slope varies from 0 to 87.76 , while the mean slope is 22.23 . Annual precipitation ranges between 100 and 500 mm in the lowlands to a maximum of 2000 mm on mountain slopes. Whereas, the average annual rainfall in the Indus plains is about 230 mm, and snowfall at higher altitudes accounts for most of the river runoff (Ojeh 2006).

Historic GLOF events
Many researchers worldwide had reported and documented historic GLOF events occurred across the Third Pole over the century or beyond. Whereas, in a past decade, those reported historic events were restudied and verified again using local geomorphology and satellite data by various researchers, and several inventories of such events with detailed documentation has been prepared (Hewitt and Liu 2010;Emmer 2018;Nie et al. 2018;Bhambri et al. 2019;Veh et al. 2019;Bazai et al. 2021;Zheng et al. 2021). To get an overview of the occurrence and impact of GLOF in this present study, a list of historic events has been compiled and prepared using up-to-date inventories available on GLOF events mentioned above, (which only enlists and depicts total nos. of events occurred in the Indus River basin), including one most recent GLOF event occurred on 21 st August 2021 in the Rumbak valley (near Leh, India) from a lake of area 0.03 km 2 (located approx. 17 kms upstream from its confluence with the Indus River), has been shown in Table S1 (supplementary material). On the basis of which, it has been found that, in the Indus River basin, GLOF events are occurring since 1533, with a total of 153 events from 29 sites till present, spatial distribution of which has been represented in Figure S1 (supplementary material).

Data preparation
There are many earth observation satellites capturing data repeatedly in various spectral ranges and at different spatial and radiometric resolution. For inventorying glacial lakes, high to medium resolution datasets are proved to be useful, and mainly satellite data of Landsat series has been used widely for mapping glacial lakes due to its wide coverage and free accessibility (Bolch et al. 2010;Mergili et al. 2013;Wang et al. 2013b;Zhang et al. 2015). Indian Remote Sensing (IRS) satellite data of sensors viz., AWiFS, LISS-III, and LISS-IV provides medium to high resolution imageries, and has also been used for making an inventory of glaciers and glacial lakes (NRSC 2011;Panigrahy et al. 2012;Gupta et al. 2019;Guru et al. 2019. In this present work, 157 high resolution orthorectified Resourcesat-2 LISS-IV (RS2-L4) scenes having spatial resolution of 5.8 m were procured via NRSC user order processing system web portal (https://uops.nrsc.gov.in/). Entire RS2-L4 coverage over the study area has been represented in Figure S2 (supplementary material). The RS2-L4 data selected for mapping lakes were from the months September to December of 2015-17 year. During this period in each year, RS2-L4 images featured less perennial snow and cloud coverage. Data were extended to previous years when cloud obscuration and snow coverage were high during 2015-17 and the maximum time span of the selected RS2-L4 data has been only two preceding years. Some RS2-L4 images for winter months were also included when data for previous years were unavailable. Additional details of data used in the study have been represented in Figure 2. The data predominantly spanned the September-December months for the period 2015-17 (78%), with 9%, 32%, 25%, and 14% of the scenes imaged during September, October, November, and December respectively. Digital Elevation Model (DEM) of Shuttle Radar Topography Mission (SRTM) with 30 m spatial resolution has also been used to generate topographic information and subbasin boundaries. Additionally, Randolph Glacier Inventory v6.0 has been used to identify what could be considered as glacial lakes (RGI Consortium 2017). Other information like name of lake and river has been gathered from open source digital toposheets available from Texas University Toposheet Library at 1:250,000 scale and Tibet Map Institute at 1:100,000 scale (U.S. Army Map Service 1955; Andre 2017). Google Earth and Google Map has also been used extensively for identifying lakes under shadow and/or snow-covered areas along with their types, and to generate attribute information respectively.

Glacial lake inventory
Presence of cloud, snow, and mountain shadows in the study region has complicated the automated satellite image classification of lake boundaries, and hence, glacial lakes were carefully identified and manually digitized, hence guaranteed high quality control. To demarcate lake boundaries clearly, an image processing technique that redistributes pixel values in order to adjust the contrast i.e. histogram equalization (nonlinear contrast stretching technique) has been used to enhance the satellite image. Using RS2-L4 high resolution data, lakes of size up to 0.001 km 2 are easily possible to identify due to more than 5 Â 5 pixel coverage, but incur more time due to labour-intensive work. Thus, in this present study, a threshold of ! 0.0025 km 2 is applied to identify glacial lakes so as to map lake boundaries clearly. The presence of glacial lake has been identified preliminary by using Normalized Difference Water Index (McFeeters 1996) and only those were considered which are within 10 km buffer distance from RGI glacier; as determined to be used in previous studies also Wang et al. 2020). However, glacial lakes in areas where glaciers are missing have been considered and mapped to ensure the veracity and consistency of the inventory. Uncertainty estimation of glacial lake area is crucial as the present study has been conducted using remote sensing data, and hence, estimation of error is important to determine the accuracy and significance of results (Bhambri et al. 2011;Mir et al. 2017;Jain and Mir 2019). To estimate uncertainty, developed due to manual digitization over a sensor resolution, Linear Resolution Error (LRE) is adopted, by assuming that the average lake margin passes through the centre of the pixel along its perimeter (Salerno et al. 2012;Thakuri et al. 2016). Gorman (1996) states the uncertainty of half a pixel (± 0.5 pixel) can be assigned for manual digitization. Consequently, in this study, we consider uncertainty of half pixel, thus adopting the LRE of ± 2.9 m and has been using the worst LRE of þ 2.9 m as buffer for uncertainty assessment of glacial lake extents.
Various researchers have proposed glacial lakes classification schemes based on dam type, process of formation, topographic feature, and geographical position (Liu et al. 1988;Mool et al. 2001aMool et al. , 2001bHuggel et al. 2002;Yao et al. 2018;Wang et al. 2020). Lakes located on the glacier surface can be mapped using satellite data, but there are englacial and subglacial lakes that may also exist, but cannot be mapped from aerial/optical satellite images, it requires a ground-based instrument (Yao et al. 2018). Following ICIMOD (2011) classification of glacial lakes, this study inventoried and classify glacial lakes broadly in 4 classes and 10 subclasses, i.e. a) moraine-dammed lake (end-moraine, lateral moraine, lateral moraine with ice, and other moraine dammed lake) which typically develop in over deepened glacier bed behind the terminal moraine, and also between the lateral moraine and the valley side; b) ice-dammed lake (glacier ice-dammed and supraglacial lake) which are dammed by glacier ice, either directly by another glacier or on surface of a glacier; c) glacier erosion lake also known as bed-rock dammed lake (cirque, glacier trough valley, and other glacial erosion lake) refers to the water body in the depression formed by erosion and abrasion of glacier; d) other glacial lakes which are formed in a glaciated valley, and fed by glacial melt, but damming material not directly part of the glacial process. Once mapping is completed manually, a total of 22 attributes has been given to all mapped features, broadly consisting of hydrological information (basin, subbasin, river name, and lake name), geometrical information (maximum length, mean width, surface area), geographical information (latitude, longitude, region, state/UT, district, and 250 K & 50 K toposheet no.), topographical information (elevation and aspect), lake information (feature types, glacial lake type, and lake ID), and data source information (source of database, source of elevation, and RS2-L4 date of pass).

Satellite derived parameters
Parameters used for identifying PCGLs in the Indus River basin have been derived using remotely sensed data for two purposes, namely 'preliminary screening' and 'quantitative assessment'. Preliminary screening has been done sequentially using 4 parameter criteria viz., lake type, lake area, surrounding characteristics, and water bodies in the river course. These 4 parameters were examined and derived using satellite data, the former two while generating attributes for the inventory, while the latter two only for lakes screened-in after applying the first two criterions. First, under a background of criticality of different damtypes, a) end moraine-dammed lakes appear to be the most common type for GLOF in the Third Pole region (Ashraf et al. 2012, Nie et al. 2018Zheng et al. 2021); b) laterally dammed lakes outburst due to sudden emergence of upstream lake(s) or by intensive precipitation events (Rao et al. 2014;Champati Ray et al. 2016); c) glacier ice-dammed lakes recorded numerous GLOFs due to increase in water level but mostly concentrated in Karakoram (Bhambri et al. 2019; d) other moraine dammed and supra-glacial lakes which are closely-spaced with same or another category of lake on a receding valley glacier can be considered as dangerous lakes, in case of upstream flow mergence (Mool et al. 2001a;Bajracharya et al. 2007aBajracharya et al. , 2007b; e) cirque erosion lakes being stable in nature but if associated with steep hanging glacier within its vicinity can also be considered as dangerous lakes, in case of mass impact causing overtopping waves (ICIMOD 2011;Ashraf et al. 2012;Worni et al. 2013;Gupta et al. 2019). Hence, all moraine-dammed and ice-dammed lakes along with cirque erosion lakes have been considered for preliminary screening.
Second, only a glacial lake greater than a certain size could possibly cause damage (Wang et al. 2013a). Hence, setting up a threshold for lake area of more than equal to 0.02 km 2 to obtain limited nos. of lake for prioritization, as lakes of this size should also be considered as dangerous lakes for identifying PCGLs (Ashraf et al. 2012). Third, surrounding characteristics which invite different indicators for different dam types, we evaluate this qualitatively separately for moraine-dammed and non-moraine-dammed lakes. As moraine-dammed lakes could be associated with a glacier but all may not have potential rockfall/avalanche sites nearby. Similarly, non-moraine-dammed lakes like supra-glacial lake and cirque erosion lake may or may not have lake(s) upstream and hanging glacier nearby at steep slopes respectively. Therefore, those with potential sites upstream have been considered. Fourth, despite being large or small in size, if a lake drains into a big naturally impounded water body without having any habitats and infrastructure in its course before mergence, then the possible susceptibility of a glacial lake could be negligible. Hence, lakes which have exposed settlements or hydropower sites in their course of flow are only considered.
For quantitative assessment, next 4 key parameters viz., distance and slope between glacier and lake (L€ u et al. 1999;Wang et al. 2012bWang et al. , 2012c, and distance and slope between lake and nearest exposed site (Bajracharya et al. 2007a), along with two formerly used parameters viz., lake type and lake area, has been selected to quantitatively evaluate glacial lakes obtained after preliminary screening and prioritize them in order of their potentiality. These 6 selected parameters have been included because they met criteria for evaluating the severity of exposed sites in case the event occurs. First, dam-type and lake area have been used worldwide to access critical glacial lakes (ICIMOD 2011;Worni et al. 2013;Wang et al. 2013a;Zheng et al. 2021). Second, the smaller the distance and higher the slope between a) lake and its parent glacier snout; and b) lake and exposed site in its course, more associated their dynamics will be (L€ u et al. 1999;Bajracharya et al. 2007a;Wang et al. 2012cWang et al. , 2013a. Third, all these 6 parameters can be derived using satellite data and could be used quantitatively for prioritization of glacial lakes. However, lake type being the only qualitative parameters in the assessment, ranking them [M(e), M(l), M(lg), M(o), I(d), I(s), E(c)] on the basis of their criticality explained above, between 7 and 1 (where 7 being the highly critical and vice-versa), is a way possible to prioritize in a quantitative manner.

Prioritization of PCGLs
Prioritization of preliminary screened-in lakes has been quantitatively assessed for two scenarios. First, scenario-1 comprises lakes which may self-destructively breach due to their high storage capacity i.e. preliminary screened-in lakes with area ! 0.1 km 2 . Second, scenario-2 involves all preliminary screened-in lakes with area ! 0.02 km 2 due to the fact that several recent GLOF incidents occurred from small lakes in conjunction with external triggering events like extreme precipitation, mass movement into the lake such as ice/rock avalanche, landslides, earthquakes, etc. Prioritization for both scenarios has been done using two integrated index methods viz., Equal Weights (EW) and Unequal Weights (UEW). Quantitative assessment is usually done by constructing a vulnerability index ( Zurovec et al. 2017) which is based on several sets of indicators that result in vulnerability of a region. Each indicator used is measured in different scales and units. Therefore, all measured parameter values need to be normalized (certifying that they are comparable) based on two types of functional relationship, between the indicators and vulnerability (UNDP 2006). If vulnerability increases with the increase in the value of the indicator (positive correlation x ij ), normalisation was carried out by using Eq. 1 for parameters viz., lake type, lake area, slope b/w glacier and lake, and slope b/w lake and nearest exposed site. If vulnerability decreases with an increase in the value of the indicator (negative correlation y ij ), the Eq. 2 will be used for parameters viz., distance b/w glacier snout and lake, and distance b/w lake and nearest exposed site.
where, X ij is the normalized value of parameter (j) with respect to lake (i), X i is the actual value of the parameter with respect to lake (i), and Min X ij and Max X ij are the minimum and maximum values, respectively, of parameters (j) among all the lakes (i).

Equal weights method
Weights are assigned based on their degree of influence, where simple average of all normalized scores is used to construct the composite index (CI) using number of parameters (K), and subsequently rank them:

Unequal weights method
A statistically sound and well-suited method proposed by Iyengar and Sudarshan (1982) has been used for the development of composite index, by using given formula: where, y i is a linear sum of normalised score i.e. x ij , and w is the weight determined by where, c is a normalizing constant and determined by The advantage of this method is it would ensure that large variation in any one of the indicators would not unduly dominate the contribution of the other parameters. The computed composite index using both methods, lies between 0 and 1, with 1 indicating maximum potentiality and 0 represents no potentiality. Assessment of potentiality has been done using logarithmic beta distribution, which is a suitable classification for probability distribution to obtain meaningful characterization of the different stages of potentiality, as distribution will be skewed and takes values in the interval (0, 1). The probability density is given by: where, B a, b ð Þ ¼ Ð 1 0 x aÀ1 1 À x ð Þ bÀ1 dx, a and b are the lower and upper bounds of composite index and a, b are beta functions. Beta distribution values obtained using above formulae, are divided into three and five linear ranges, such that each interval has the varied probability weight depending upon the density (number of lakes), for scenario-1 and scenario-2 respectively. These intervals have been used to categorise three stages of potentiality (%) due to less nos. of lakes analysed in scenario-1 viz., low potential (0-33), moderate potential (33-66), and high potential (66-100) and five stages of potentiality (%) due to more than 100 lakes analysed in scenario-2 viz., low potential (0-20), moderately low potential (20-40), moderate potential (40-60), moderately high potential (60-80) and high potential (80-100).
Criticality of a potential lakes is defined not only by its dam type and geometry, but also on its storage volume. Lakes with a large water volume definitely possess high risk to the downstream and vice-versa. Kougkoulos et al. (2018) analyzed a list of several potential lakes identified by many researchers globally along with several lakes which had outburst history, for assessing GLOF risk level. The model applied by Kougkoulos et al. (2018), asserts that lakes with the storage volume < 1 MCM have been categorized as low risk lakes, between 1 and 10 MCM are medium risk lakes, and ! 10 MCM are considered as high risk lakes. These three levels of risk obtained by categorizing the lake volume has been used in this study to identify the risk level of lakes in the Indus River basin also. But estimating the exact capacity of water contained in a lake, remotely, is a challenging task, as field investigation is required for bathymetry data collection and it is indeed complicated due to the lake's inaccessible geographical location and inclement weather condition. Hence, it needs to be derived by using empirical relationships. But as all glaciers are meant to have complex morphometries, meaning the associated lake to it may also have complex bathymetries, and hence any relationship developed could also be more unpredictable (Cook and Swift 2012). As a result of that, a handful number of studies had been conducted by many researchers to calculate approximate volume of lake(s) using remote sensing imageries in conjunction with bathymetric data, resulting an empirical formula based on relationships between lake depth, area, and volume. Details of such empirical formulas has been given in Table 1. These depth-area-volume relationships allows rapid and simple calculation of lake capacity from broadly existing satellite images, whilst avoiding the requirement of often challenging fieldwork. One such relationship developed by Huggel et al. (2002) using various types of lakes situated globally has gained significant attention and has been used by many researchers in several studies to estimate lake volume (Huggel et al. 2004;Bolch et al. 2011;Mergili and Schneider 2011;Jain et al. 2012;Byers et al. 2013;Gruber and Mergili 2013;Wilcox et al. 2013;Che et al. 2014). However, this equation had also been modified for specific locations by Loriaux and Casassa (2013) and Yao et al. (2012c). However, there is still a question on applying these empirical relationships confidently across a range of locations and to various types of lakes (Cook and Quincey 2015). Various assessments have been done to test the correlation and associated error between the lake volume calculated by using the above empirical relationships and the interpolated bathymetric measurements which was compiled from various published data sets. Findings of Cook and Quincey (2015) suggest that, empirical relationship of Huggel et al. (2002) gave reasonable approximations of lake volume, followed by that of Evans (1986). Whereas, Patel et al. (2017) derives a new relationship to estimate volume of moraine-dammed lakes which shows minimum difference with that of the observed volume despite large differences in volume derived by all other empirical relationships. However, it has been observed that, relationship of Patel et al. (2017) overestimates the lake volume for other types of lakes than that of moraine dammed lakes. In this present study, in order of high correlation obtained with that of the observed data in their own individual studies, lake volume has been estimated by three above relationships developed by Patel et al. (2017), Huggel et al. (2002), Evans (1986), to identify suitable risk level associated with a lake.

Glacial lake inventory
A total of 5335 glacial lakes (! 0.0025 ± 0.0006 km 2 ) have been inventoried in the Indus River basin, covering a total area of 173.95 ± 10.13 km 2 . Out of all 10 types of glacial lakes, Other Glacial Erosion Lakes are found to be the maximum in the entire study area with 2516 nos. (47.16%) covering a total area of 47.42 ± 0.10%. Two other categories of lake, namely, Other Moraine Dammed and Supra-glacial lakes are 905 (16.96%) and 809 (15.16%), extend over an area of 12.12 ± 1.22 km 2 (6.97 ± 0.28%) and 6.09 ± 0.93 km 2 (3.50 ± 0.32%) respectively. Detailed information on various types of lakes, their count, and total area has been summarized in Table 2, along with its spatial distribution categorised in 4 major classes of lake type shown in Figure 3. "Glacial Lake Atlas of Indus River  Huggel et al. (2002) has been finally used for estimating lake storage volume, Ã Area (A) used for Volume (V, m 3 ) estimation is in m 2 except for Sakai (2012) relationship where Area (A) is in km 2 .
Basin" was published in December 2020, contains more detailed information and distribution of glacial lakes in the Indus River basin (Rao et al. 2020). Glacial lake database can be accessed for visualisation and research purposes from GLOF Geoportal via link https:// bhuvan.nrsc.gov.in/nhp.
The size of glacial lakes in the basin ranges from a minimum of 0.0025 ± 0.0006 km 2 to maximum 2.62 ± 0.02 km 2 . Distribution of glacial lakes in 6 classes of area ranges viz., 0.0025-0.005 km 2 , 0.005-0.01 km 2 , 0.01-0.05 km 2 , 0.05-0.1 km 2 , 0.1-0.5 km 2 , and ! 0.5 km 2 is shown in Table S2 (supplementary material). Majority of lakes i.e. 4633 (86.84%) are small in size i.e. less than 0.05 km 2 , covering a total area of 34.04 ± 1.63%, and remaining 65.96 ± 1.63% of area has been covered by lakes having area ! 0.05 km 2 .  All lakes are situated in elevation ranging from 2521 m to 6086 m. Distribution of lakes in 4 elevation ranges viz., low altitude (up to 3000 m), medium altitude (3001-4000 m), high altitude (4001-5000 m), and very high altitude (> 5000 m) is shown in Table S3 (supplementary material). It has been observed that, majority of glacial lakes are situated above 4000 m elevation range i.e. 4657 (87.29%). Whereas, remaining 12.71% glacial lakes are below 4000 m elevation. Glacial lake distribution in the Indus River basin has also been analysed as per region-wise, details of which are shown in Table S4 (supplementary material). A total of 4280 (i.e. 80.22%) glacial lakes lies within Indian region covering 79.26 ± 0.02% of the total lake area, whereas remaining 19.78% of lakes are located in transboundary region with a total of 20.74 ± 0.02% area. Within India, Ladakh (UT) shares 75.21% of lake, followed by 12.75% and 11.99% in Jammu & Kashmir (UT) and Himachal Pradesh respectively, with a total area of 72.28 ± 0.15%, 20.90 ± 0.28%, and 6.73 ± 0.14% respectively. Whereas, remaining 0.05% of glacial lakes are found in Uttarakhand state, covering only 0.09 ± 0.001% of lake area within India.

Preliminary screening
Initial step was to reduce the number of lakes from the total inventory for prioritization. Hence, for preliminary screening of total glacial lakes mapped i.e. 5335 has been done step-wise. Firstly, 2613 screened-in on the basis of lake type. Then, 614 glacial lakes have been screened-in considering lake area. Further, these lakes were checked for surrounding characteristics which depict only 391 lakes to have potential triggering site in their surroundings. Lastly, 367 lakes were found out to have been flowing with exposed sites nearby to its course and are not draining into a large water body. Therefore, 367 preliminary screened-in glacial lakes has been used for further analysis in prioritization of PCGLs.

Prioritization of PCGLs
From 367 lakes, quantitative analysis for scenario-1 has been done for 73 lakes only (with area ! 0.1 km 2 ), whereas all 367 lakes have been analysed for scenario-2. Lakes for both scenarios undergo assessment using 6 key quantitative parameters by applying EW and UEW methods. For both scenarios, in the EW method, weights has been distributed and assigned equally amongst each parameter i.e. 16.66%. Whereas, in case of UEW method, parameter weight will be different for both scenarios as it depends upon density and variance. In scenario-1, weights obtained for lake type, lake area, distance b/ w glacier and lake, slope b/w glacier and lake, distance b/w lake and nearest exposed site, and slope b/w lake and nearest exposed site are 19.03%, 17.92%, 20.58%, 17.17%, 11.88%, and 13.42% respectively. Whereas, in scenario-2, weights obtained are 21.51%, 22.33%, 18.45%, 13.68%, 11.89%, and 12.15% respectively. As a result, in scenario-1, EW and UEW methods categorise, a) 19 and 20 lakes under high potential category, b) 29 and 30 moderate potential lakes, c) 25 and 23 are categorised as less potential lakes respectively. Whereas, in case of scenario-2, EW and UEW methods results in, a) 53 and 48 lakes are under high potential category, b) 95 and 92 lakes as moderately high potential, c) 87 and 102 are categorised as moderate potential lakes, d) 74 and 73 lakes as moderately less potential, e) 58 and 52 as less potential lakes respectively, has also represented in Table 3.

Glacial lake inventory
This study identified 5335 glacial lakes (! 0.0025 km 2 ) in the Indus River basin, which asserts the highest numbers of glacial lakes mapped so far. In comparison with previous studies, our inventory depicts more lakes than others, for ex: inventory of ICIMOD pre- whereas we mapped 702 lakes (! 0.05 km 2 ). As mapping is done manually using high resolution satellite data, this inventory of glacial lakes consists of features with clear boundaries.
Even lakes which are frozen and are fully or partly covered by mountain shadow, have been mapped in conjunction with Google Earth. Places where glacier mass is not visible but traces of it can be seen, lakes within the vicinity of such glaciers have also been mapped. Furthermore, we classify 10 glacial lake types using 3D perspective view on Google Earth, hence types given to all lakes in this updated inventory is more diversified since the ICIMOD 2018 inventory, and also in comparison to the most recent inventory by Chen et al. (2021) and Zheng et al. (2021). However, some degree of subjectivity in classification of lake types cannot be avoided due to the complexity amongst themselves. We identified 1582 moraine-dammed lakes in the Indus River basin, which are 60 more in nos. than ICIMOD had reported in their 2018 inventory, which is reliably consistent. Moraine-dammed lakes are majorly small in size (< 0.05 km 2 ) i.e. 1417 (89.57%) and only 165 are large in size (! 0.05 km 2 ). But only 104 moraine-dammed lakes are situated below elevation 4000 m, and predominantly found above 4000 m elevation.

Preliminary screening
Glacial lake inventory contains a large number of lakes of varied size and types, and all cannot be taken up for detailed investigation, thus making preliminary screening inevitable. Depending on the variety of parameters involved in analysis, criticality of a lake can be defined. Not only geomorphological and geo-technical characteristics makes a lake critical, other criterions like condition of lakes, damming material, associated glaciers, topographic features around the glaciers and lakes etc., along with long-term trend of temperature and precipitation, also defines the criticality (ICIMOD 2010). Analysis of trend and spatio-temporal variability of precipitation and temperature has been assessed using various techniques like Mann-Kendall test, innovative trend analysis, inverse distance weighted, ordinary and Empirical Bayesian kriging techniques, etc., using measured and modelled data in many plain regions of India (Gupta et al. 2017a(Gupta et al. , 2017bMachiwal et al. 2019;Gupta et al. 2020;Machiwal et al. 2020) for various applications, and also in hilly regions like Himalayas (Yao et al. 2012b;Mir et al. 2014, Yang et al. 2014 for relating glacier shrinkage to increase in glacial lakes area. Depending upon the availability and study extent, various parameters consisting of lake and glacier morphometric characteristics like lake surface area, nearby avalanche site, moraine dam height, mother glacier area, slope, aspect, etc. can be obtained from aerial observation and satellite remote sensing (Casassa et al. 2002;K€ a€ ab 2002;Singh et al. 2014), field survey for measuring lake bathymetry (Cook and Quincey 2015), derived modelled products like annual average precipitation and temperature (Fick and Hijmans 2017;Kougkoulos et al. 2018), and relationship based empirical formulas (e.g. Box and Ski 2007).
In this present study, parameters for preliminary screening and quantitative assessment have been used which can only be derived using satellite data and more importantly parameters comprise information of glacier, lake, and downstream characteristics. Lake type has been considered for both preliminary screening and quantitative assessment because it has been used worldwide to assess potential glacial lakes in glaciated regions (Worni et al. 2013;Wang et al. 2013a;Aggarwal et al. 2017;Prakash and Nagarajan 2017;Zheng et al. 2021). Studies conducted by Wang et al. (2012aWang et al. ( , 2013a to choose PCGLs emphasises to set a threshold of 0.1 km 2 . Whereas, Nie et al. (2018) suggested to use a threshold of 0.05 km 2 , given that considerable damages could also be caused by the outburst of smaller lakes, such as Chorabari Lake (Rao et al. 2014), Zanaco (Bajracharya et al. 2008;ICIMOD 2011;Yao et al. 2014), and Geiqu (Yao et al. 2014). Ashraf et al. (2012) emphasises on considering lakes greater than 0.02 km 2 for potentiality assessment as in conjunction with any hazardous event, such small lakes are reasonable enough to cause damage to its nearby downstream surrounding area. Whereas, lakes having area greater than 0.1 km 2 , are considered to be potentially dangerous and can cause high damage to the downstream settlements, infrastructures, habitats, cropped land etc., even if they are situated far away from the lake (Mool et al. 2001a;Bhagat et al. 2004;Huggel et al. 2004;Bajracharya et al. 2007b;Ashraf et al. 2012;Worni et al. 2013;Allen et al. 2016;Gupta et al. 2019;Zheng et al. 2021).
However, Ashraf et al. (2012) asserts that for identification of PCGLs, those associated with glaciers such as supraglacial, cirque and/or dammed by lateral-moraine or endmoraine with an area > 0.02 km 2 is usually considered as dangerous lakes. But this damage could be heightened if those lakes are associated with a hanging glacier or upstream lakes, and could be negligible if the lake drains into a big water body without having any habitats and infrastructure in its course. Distance and slope b/w glacier and lake along with distance and slope b/w lake and nearest exposed sites gives details about the proximity of a lake with respect to its parent glacier and the nearest downstream settlement or infrastructure, along with their respective gradient. Studies by Wang et al. (2011), Shijin et al. (2015, Cook et al. (2016), Alean (1985) asserts considering distance and slope b/w glacier and lake as this indicate the possibility of ice/snow avalanches and ice calving into the glacial lake which could also generate overtopping waves. However, studies by Haeberli (1983) Huggel et al. (2002), Hegglin and Huggel (2008) suggested using distance and slope b/w lake and nearest exposed sites as it implies the downstream impacts such as potential loss of human lives and damage to infrastructure if situated in a relatively close proximity to downstream of the lakes then the risk was considered as high.

Prioritization of PCGLs
The current study, in contrast to previous studies on assessing GLOF hazard and risk qualitative, has implemented the first quantitative assessment of glacial lakes to prioritize them in order of their impact to the downstream exposed sites. Integrated index-based approach has been considered as it is based on several sets of indicators, which produces a single index number, specifying good internal correlation between them. To obtain this, EW and UEW methods have been applied in this present study. In the case of EW, weights are assigned equally to all parameters, whereas in UEW, different weights to each parameter have been computed using variance and density, giving a more reliable estimate of weights. In a recent study by Zheng et al. (2021) asserts considering EW methods for assessing potential GLOF hazard and risk. But as the EW method gives equal weightage to all parameters which may not be a true case always for any scenario, it should not be considered for index estimation for such studies. Here we propose to consider the UEW method for such assessment, as it estimates weight on the basis of nos. of lake taken up for assessment in conjunction with the variability in data for each parameter used.
Therefore, in the present study, the UEW method is adopted for prioritization of lakes in various categories of potentiality, resulting in 20 and 48 high potential lakes for scenario-1 and scenario-2 respectively. All 20 lakes of scenario-1 have been found in common with those of scenario-2 but with a difference in their rank, hence considering all 48 high potential lakes as critical. From the empirical relationship given by Patel et al. (2017), Huggel et al. (2002), Evans (1986) to estimate volume, 10 out of 48 lakes are with storage volume ! 10 MCM. However, 18 lakes have estimated storage volume between 1 and 10 MCM, and 20 lakes are within < 1 MCM. It has been observed that, volume estimation by using relationship of Huggel et al. (2002) gives more reliable estimates in comparison to other relationship, as it performs well for all types of lakes globally, whereas Patel et al. (2017) overestimates volume for ice-dammed and cirque erosion lakes. However, estimates of Evans (1986) is close to that of Huggel et al. (2002), but is only lake specific, which is also not more in nos. in this study area. Therefore, for volume estimation, we have considered the empirical relationship of Huggel et al. (2002). According to Kougkoulos et al. (2018) categorization of risk level associated with lake volume, 10 lakes are flagged as high risk lakes, 18 as medium risk, and 20 as low risk lakes in the Indus River basin.
But due high bursting potential and self-destructive nature of end-moraine dammed category of lakes, only 6 lakes are being considered as PCGLs out of 10 high risk lakes. Whereas, remaining 4 lakes (2 cirque erosion, 1 glacier ice-dammed, and 1 other-moraine dammed lake) also have volume ! 10 MCM has also been flagged as high-risk lakes but not considered as PCGLs, as cirque erosion and other-moraine dammed lakes are stable in nature and probability of bursting is very less, whereas ice-dammed lake is dammed by glaciers on both sides, and there is no possibility of breaching and overtopping of this lake. However, rest 38 lakes with volume < 10 MCM, along with remaining 4 high risk lakes, are considered as critical lakes only in case of dynamic mode of failure which includes external triggering events, mass movement into the lake, etc. Spatial distribution of all 48 high potential lakes of scenario-2 which includes 20 high potential lakes of scenario-1 also has been represented in Figure 4 and their approximate lake storage volume estimated by three established relationships along with their categorised risk level is given in Table S5 (supplementary material).

Conclusions
In this present study, we demonstrated the efficacy of high resolution Indian remote sensing satellite imageries to map glacial lakes ! 0.0025 ± 0.0006 km 2 and effectiveness of satellite derived quantitative parameters for prioritization of PCGLs. This updated glacial lake inventory provides various information through 22 attributes, highest in nos. of all available inventory, making it an authentic database to be used by Central and State Disaster Management Authorities and other stakeholders, for regular monitoring of changes in spatial extent. The implemented methodology assesses quantitative parameters based on literature review which actually prosper potentiality of the lake due to its surrounding characteristics and criticality due to availability of exposed sites in the downstream valley. To ease out the analysis, preliminary screening helps to reduce the number of lakes for quantitative assessment. Furthermore, results indicated that a total of 48 lakes are of high potential in the Indus River basin, out of which 6 are considered as PCGLs due to their bursting potentiality. Risk level which is just a representative of storage volume in this study indicates 10 high, 18 medium, and 20 low risk lakes. For 6 PCGLs, it is very important to carry out GLOF modelling and associated risk to downstream people and property for planning suitable disaster mitigation measures. Whereas, other high potential lakes should be of high priority for detailed field investigation, moderate potential lakes less priority for field visit but requiring regular monitoring, and low potential lakes could be observed sporadically.
This potential behaviour of lakes has been evaluated using a quantitative approach, used for the first time for prioritization of PCGLs, where categorisation of integrated index defines the potentiality. In this analysis, unequal weight methods assert satisfactory results in assessing large and small sized lakes. We recommend that the use of surrounding characteristics should not be used only for outburst potentiality assessment, but also as a screening parameter, and proximity and gradient of downstream exposed sites from lake as an assessing parameter. Furthermore, we also recommend the use of unequal . Spatial distribution of 48 high potential lakes of scenario-2, which includes 20 high potential lakes of scenario-1 also, representing types of glacial lake along with their size distribution in the Indus River basin, subbasin and administrative boundaries (where HP is Himachal Pradesh, JK is Jammu & Kashmir (UT), LA is Ladakh (UT), UK is Uttarakhand, and TR is Transboundary Region).
weight method for prioritization of PCGLs, as it gives quantitative weights depending upon population density and variance, instead of qualitative and subjective weights. Additionally, as various high resolution digital elevation models become available, such assessment should be repeated periodically to significantly improve the identification and modeling of severe phenomena surrounding the lake and glacier, and then subsequent risk modeling for dynamic failures.