Glacial Lake Outburst Flood Hazard and Risk Assessment of Gangabal Lake in the Upper Jhelum Basin of Kashmir Himalaya Using Geospatial Technology and Hydrodynamic Modeling

Climate warming-induced glacier recession has resulted in the development and rapid expansion of glacial lakes in the Himalayan region. The increased melting has enhanced the susceptibility for Glacial Lake Outburst Floods (GLOFs) in the region. The catastrophic failure of potentially dangerous glacial lakes could be detrimental to human life and infrastructure in the adjacent lowlying areas. This study attempts to assess the GLOF hazard of Gangabal lake, located in the Upper Jhelum basin of Kashmir Himalaya, using the combined approaches of remote sensing, GIS, and dam break modeling. The parameters, such as area change, ice thickness, mass balance, and surface velocity of the Harmukh glacier, which feeds Gangabal lake, were also assessed using multitemporal satellite data, GlabTop-2, and the Cosi–Corr model. In the worst-case scenario, 100% volume (73 × 106 m3) of water was considered to be released from the lake with a breach formation time (bf) of 40 min, breach width (bw) of 60 m, and producing peak discharge of 16,601.03 m3/s. Our results reveal that the lake area has increased from 1.42 km2 in 1972 to 1.46 km2 in 1981, 1.58 km2 in 1992, 1.61 km2 in 2001, 1.64 km2 in 2010, and 1.66 km2 in 2020. The lake area experienced 17 ± 2% growth from 1972 to 2020 at an annual rate of 0.005 km2. The feeding glacier (Harmukh) contrarily indicated a significant area loss of 0.7 ± 0.03 km2 from 1990 (3.36 km2) to 2020 (2.9 km2). The glacier has a maximum, minimum, and average depth of 85, 7.3, and 23.46 m, respectively. In contrast, the average velocity was estimated to be 3.2 m/yr with a maximum of 7 m/yr. The results obtained from DEM differencing show an average ice thickness loss of 11.04 ± 4.8 m for Harmukh glacier at the rate of 0.92 ± 0.40 m/yr between 2000 and 2012. Assessment of GLOF propagation in the worst-case scenario (scenario-1) revealed that the maximum flood depth varies between 3.87 and 68 m, the maximum flow velocity between 4 and 75 m/s, and the maximum water surface elevation varies between 1548 and 3536 m. The resultant flood wave in the worst-case scenario will reach the nearest location (Naranaag temple) within 90 min after breach initiation with a maximum discharge of 12,896.52 m3 s−1 and maximum flood depth and velocity of 10.54 m and 10.05 m/s, respectively. After evaluation of GLOF impacts on surrounding areas, the area under each inundated landuse class was estimated through the LULC map generated for both scenarios 1 and 2. In scenario 1, the total potentially inundated area was estimated as 5.3 km2, which is somewhat larger than 3.46 km2 in scenario 2. We suggest a location-specific comprehensive investigation of Gangbal lake and Harmukh glacier by applying the advanced hazard and risk assessment models/methods for better predicting a probable future GLOF event.


Introduction
Glaciers are highly sensitive to climatic fluctuations [1][2][3][4][5] and are therefore considered the key indicators of climate change [3,[6][7][8]. The Himalayan region is one of the few regions where climate change effects are most prominent [9,10]. Most studies have reported that Himalayan glaciers are in a state of continuous recession [11][12][13][14][15], except the glaciers in the Karakoram region, which are stable or advancing [12]. Retreating glaciers at a faster rate are believed to have triggered many catastrophic events in recent decades; for example, GLOFs are becoming a serious and emerging threat in the high mountainous regions of the world in general and the Himalayan region in particular [4,13,14,16]. GLOF is a type of cryospheric natural hazard that can cause catastrophic fatalities and damage to neighboring areas, and the Himalayan region is one of the well-known GLOF-vulnerable regions in the world [17].
The Himalayan region, being geologically young and fragile, has witnessed significant warming in the present century compared to the past century. Warming has been reported throughout the Himalayan region at rates ranging from 0.15 • C to 0.60 • C/10 yr, which is very high when compared to the global mean warming rate of 0.74 • C/100 yr [18][19][20]. As a result of this warming trend, Himalayan glaciers are losing mass at a faster pace [21][22][23], resulting in the evolution and expansion of various glacial lake types [24][25][26]. Among these, potentially dangerous glacial lakes pose a major threat to the communities and infrastructure in the neighboring regions.
The study of Himalayan glaciers has substantially advanced in the last decade, and the spatiotemporal frequency and degree of changes in glacial lakes have been reported as well [27]. According to climate model projections based on high-emission scenarios, it has been reported that approximately 65% of the glaciers' ice mass in High Mountain Asia (HMA) will vanish by the end of the 21st century [28]. As a result, the size and number of proglacial lakes (glacier-fed) and supraglacial may increase, which coalesce to form large glacier lakes. Further, it is essential to mention that most of these lakes are formed behind unconsolidated (loose) and inherently unstable moraine dams.
Moraine dam failures induced by the high rate of lake area expansion, seepage, dead-ice melting, degradation of permafrost, seismic activity, cloudburst, a flood wave generated by a rock avalanche/fall or due to the snow/ice avalanche into the lake, and an abrupt rise in the water level may cause these lakes to burst [29][30][31]. Such outburst events often have very high destructive potential, as an enormous amount of water is discharged in a short span of time with a capacity to erode the unconsolidated/loose material, potentially resulting in a strong and violent flow that can travel a long distance in the downstream regions.
Earlier, glacial lake studies in the former state of Jammu and Kashmir were carried out at different levels, for example, at the international level (glacial lake inventory of HMA) by Ives et al. [32], Wang et al. [33], and Chen et al. [34]; at the national level by the Central Water Commission [35], National Remote Sensing Center Hyderabad [36], Ives et al. [32], Fujita et al. [37], and Worni et al. [38]; and at the basin level by Dar et al. [39] and Ahmed et al. [40,41]. The Potentially Dangerous Glacial Lakes (PDGLs) in the area have been identified by seven previous studies. According to Kumar and Goel [42], Gangabal lake and Konsar-Nag lake are susceptible lakes. In recent investigations on glacier lakes in the Indian Himalayas, Mal et al. [43] categorized 25 glacial lakes as PDGLs. Among these, four lakes are located in Jammu and Kashmir; however, only one of them, namely Sona-Sar lake, is found in the area under study. Sona-Sar lake was also confirmed as a high-risk Lake (HRL) by Dubey and Goel [44] and a high-hazard glacial lake (HHGL) by Ahmed et al. [41]. Dubey and Goel [44] have reported four glacial lakes (Hirubagwan, Krishna-Sar, Gad-Sar, and Har-Nag) as very high-risk lakes (VHRL) in the upper Jhelum basin. Dimri et al. [45] have identified 40 lakes as potentially dangerous and two as critical lakes in Jammu and Kashmir. GLOFs pose a severe threat to populations living in mountainous regions. To demarcate vulnerable areas for GLOF hazard management and mitigation, it is necessary to identify potentially dangerous lakes and simulate outburst floods.
Thus, keeping in view the hazard level of the dangerous lakes assessed by various regional studies, in-depth research evaluating the GLOF hazard and risk assessment of the individual lakes is quite rare. In the present study, we particularly focus on Gangabal lake (potentially dangerous), the largest lake in volume and area, located in the Sind catchment of the UJB of the northwestern Himalayas. The specific objectives of the current study include: (1) analyze the spatiotemporal dynamics of the lake (Gangabal) and its feeding glacier (Harmukh) by utilizing multitemporal satellite data between 1972 and 2020; (2) study the surface velocity, geodetic mass changes, and ice thickness of the Harmukh glacier; (3) use coupled HEC-RAS 1D and 2D unsteady flow analysis to model a worstcase GLOF scenario and locate exposed parts in the closest settlement, which is situated 9.2 km from the outlet of the lake; and (4) analyze the impact of the probable GLOF event on downstream landuse and landcover.
GLOFs are frequently associated with glacier mass loss and are thought to be one of the most essential geomorphological and hazardous consequences of melting glaciers [61]. Considering the ongoing glacier mass-loss induced by climate change, the study of past and probable future GLOFs is a very pertinent topic [62,63]. GLOFs and their related processes have major ramifications in a variety of domains, for example, risk management, disaster risk reduction, designing hydropower plants, and geomorphology and sediment yield [64,65]. GLOFs in the Himalayas have previously resulted in catastrophic deaths and destruction in neighboring regions [66]. With the continuous glacier mass loss in the region [67], GLOF risk has been exacerbated and needs much more attention owing to the potentially catastrophic damage [17]. It is necessary to construct a database of GLOF events from the existing literature in order to expose the process and mechanism of GLOFs for mitigation, hazard assessment, and risk management of glacial lakes. We compiled an inventory of historical GLOFs in the Himalayan region (see Figure 1). A total of 158 GLOF events were documented in the Indian Himalayan region, and the analysis revealed that most of the GLOF events had been seen in the Ladakh region of the northwestern Himalayas ( Figure 1). However, few of them are repeated events.
Many glacial lakes surrounded by unconsolidated moraines pose a severe threat to the people living in the world's mountainous regions. The dams could fail unexpectedly, causing devastating floods with peak discharges far beyond average. Glacier-dammed and moraine-dammed lakes are the two main sources of GLOFs in the Himalayan region [68]. The lake outbursts occur when glacier ice initiates to float on the dammed lake or when the stored water erodes the dam laterally, subglacially, or superficially [69]. Moraine-dammed lakes are susceptible to outbursts due to a variety of conditioning factors and triggering mechanisms. The majority of these glacial lakes are inaccessible due to their remote locations, tough terrain, and harsh weather conditions. Therefore, continuous monitoring and assessing these glacial lakes at regular intervals through field-based investigations and observations are impossible. Satellite remote sensing and GIS serve the purpose of monitoring and assessing these glacial lakes at different spatial and temporal scales. Furthermore, to prevent and monitor the GLOF hazard and to evaluate the damage and destruction that may occur in the future, the integration of various geospatial techniques and modeling approaches could play a key role in identifying hazardous/dangerous lakes and monitoring GLOF events in real-time. Such robust techniques have been effectively used worldwide to monitor and assess glaciers, glacial lakes, and associated hazards. The GIS and RS-based tools and techniques reasonably allow rapid investigation of huge glaciated terrains and are therefore considered suitable for glacial lake and GLOF studies [40].
Geospatial techniques coupled with hydraulic and hydrologic simulations can be used to obtain information about the status of glaciers and glacier lakes at high spatial and temporal resolutions and map flood extent [70]. Satellite images have already been utilized in several global studies to monitor glaciers and glacier lakes, as well as hydraulic and hydrologic models such as HEC-RAS, HEC-HMS, HEC-GeoRAS, HEC-Geo-HMS, Arc Hydro, etc., for GLOF mapping. The HEC-RAS model helps us analyze the possible GLOF impacts, especially in generating flood hydrographs and several other flood wave parameters such as flood water depth, water surface elevation, and flow velocity. The key objective of this study was to use RS, GIS, and hydrodynamic modeling (HEC-RAS 5.0.7) to analyze the hazard and risk associated with the potentially dangerous glacial lake (Gangabal) located in the Sind catchment of the Jhelum basin. GLOF modeling has the ability to provide crucial insights for predicting future events, which can aid in damage assessment, mitigation, and planning.

Study Area
This study focused on the largest glacial lake, Gangabal, located in the Sindh catchment of the Jhelum basin in the Western Himalayas, with an area of 1.65 km 2 , per recent estimates. The lake is positioned between 34 • 25 55.49"N longitude and 74 • 55 27.12"E lattitude at an altitude of 3576 m asl [71]. The lake is directly fed by the meltwater of the Harmukh glacier, one of the largest glaciers in the Jhelum basin [14]. Water from the lake is directly supplied to another lake in the downstream region, Nundkol, located at a distance of 0.9 km from the outlet. It then merges with the Sind tributary through Wangath Nalah ( Figure 2). The distance between the lake and the feeding glacier is 0.78 km, with a maximum slope of 32 degrees and an elevation gain of 393 m from the upper point of the lake to the glacier front. The lake has a maximum length of 2.5 km and maximum width of 1 km. It is located in the Ganderbal district near Ganga Gali. The lake has an average depth of 42.5 m, obtained through the improved area volume scale technique. The water level of the Gangabal lakes is rising at a negligible rate of 0.009 m/yr, as observed by Srivastava et al. [72]. Gangabal lake has glacial lake id no 0143J1501489 as per the Glacial Lake Atlas of Indus River Basin prepared by NRSC in 2020 [73]. Long-term climate data (1981-2020) of the meteorological observatory nearest to Gangabal lake, i.e., Shalimar, located in the Srinagar district, shows that T max has increased by 1.79 • C from 19.0 • C (1981) to 20.79 • C (2020), with a slope of 0.04 and R2 value of 0.19 ( Figure 3). T min , on the other hand, does not depict any remarkable impact, decreasing by 0.40 • C in the last 40 years. The mean temperature has shown a significant increased mean of 0.7 • C during the previous 40 years, with an R2 value of 0.13. A decreasing trend was also observed for the annual precipitation, as depicted in Figure 3.  Gangabal lake has been classified as a potentially dangerous lake by various researchers in the past, for example, a high-risk lake by Gupta et al. [36], a potentially dangerous lake by Guru et al. [74], medium-risk lake by Ahmed et al. [41], and important GLOF lakes by Kumar and Goel [42] and Kumar et al. [15]. The GLOF susceptibility of Gangabal lake is provided in Section 5.6. Considering the GLOF potential of the lake, the present study prioritized it for hazard and risk assessment using the integrated approaches of geospatial technology and hydrodynamic modeling.

Data
Area change analysis was performed for Gangabal lake and its feeding glacier using multitemporal satellite data. In this study, we used high-resolution Resourcesat-2 LISS-IV data scenes (5.8 m), Sentinel-2A (12 m), Landsat TM (30 m), Landsat ETM (30 m), Landsat ETM+ (30 m), and Landsat OLI (30 m) to map the glacial lake outlines from 1972 to 2020 and glacier outlines from 1990 to 2020. All the procured scenes of Resourcesat-2 were downloaded from the NRSC Hyderabad web portal https://www.nrsc.gov.in/, accessed on 25 October 2021. Sentinel-2A data scenes were acquired from the "Copernicus internet hub of European Space Agency" (ESA), and Landsat imageries were downloaded from the USGS Earth Explorer data portal https://earthexplorer.usgs.gov/, accessed on 12 July 2021. Cloud cover seems to be the major hindrance during the selection of scenes; therefore, the imagery from the adjacent month or year was preferred.
Advanced Land Observing Satellite (ALOS) Phased Array Synthetic-Aperture Radar (PALSAR) radiometrically terrain corrected (RTC) Digital Elevation Model (DEM) was obtained from the web portal, i.e., https://asf.alaska.edu/ (accessed on 12 July 2021) and was utilized to obtain the topographic information (elevation, slope, aspect, hill shade, etc.) for the GLOF modeling. ALOS PALSAR has a spatial resolution of 12.5 m. It has been widely used for hydrodynamic modeling to generate geometry data such as river center lines, bank lines, flow lines, cross sections, and breach invert levels. Google Earth images were also used in this study for validation and cross-checking. SRTM-DEM for 2002 and TAN-DEM X for 2012, with a spatial resolution of 90 m, were used for the mass balance calculation of the glacier. Recent Landsat-8 panchromatic (scene-8) scenes were utilized to estimate glacier velocity. RGI glacier extents that possess the required parameters for the depth estimation were used as input in the GlabTop-2 model to derive the depth values of the glacier. The details of satellite data and other ancillary data sources used in the study, including their source, number of bands, spatial and temporal resolution, and purpose, are given in Table 1. The freely available CubeSat data downloaded from the Planet web portal https://www.planet.com/ accessed on 2 August 2021 were used to calculate various glacier and glacial lake parameters. A step-by-step methodology was employed to understand and monitor the glacier and glacial lake-associated hazards in the study region [44,66,[75][76][77][78]. The methodology for assessing the status of the Harmukh glacier and the hazard assessment of Gangabal lake and analyzing the glaciological characteristics of its feeding glacier (Harmukh) are summarized in the following section. The overall methodology of this work is illustrated in Figure 4.

Characteristics of the Adjacent Glacier (Harmukh Glacier)
The glacier mapping and change detection analysis was carried out using multitemporal satellite data. In the past studies, glacier mapping has been done using different methods and techniques. In this study, we have used NDGI technique combined with visual interpretation to identify and delineate glaciers as suggested by [3,79], which is considered the most accurate method for glacier mapping [80].

Glacier Depth
Glacier ice thickness estimation is generally obtained through mathematical equations, for example, volume-area scaling [81][82][83][84]. In this study, the depth of the Harmukh glacier was calculated using the GlabTop-2 model. GlabTop-2 is a slope-dependent, gridded ice thickness estimation model that is based on the flow mechanics of the average basal shear stress (τ), shape factor (f ), and surface slope [85,86]. The formula is as follows (Haeberli and Hoelzle [87]): Ice thickness (h) = τ f ρg sin(α) (1) where in Equation (1) h represents the ice thickness, τ indicates the basal shear stress, f denotes the shape factor, and ρ is the ice density. Glacier polygons are transformed into a raster matching the DEM cells. Cells are differentiated as inner glacier cells (light yellow), glacier marginalized cells (orange), cells adjacent to the glacier (dark cyan), and nonglacier cells (white). Dark green cells depict randomly selected cells (r) for which local ice thickness is measured; the black square represents the buffer of variable size, which is enlarged (black-dashed square), until an elevation extent of hmin is attained inside the buffer. The value for parameter f is considered 0.8, which generally works well for Himalayan glaciers [86,88]. The hmin is considered as 20, and fractions of cells are taken as 30% for random cell pick-up, g is the acceleration due to gravity (9.8 m s −2 ), α is the slope, and ρ is the ice density, which is 900 kg m −3 [89]. A similar projection system is necessary for all the datasets used as an input file in the Python-based GlabTop-2 model to run the code smoothly. The working principle of the GalbTop-2 model is illustrated in Figure 5.

Glacier Surface Velocity
Glaciers move downslope under the influence of gravity and their underlying mass, which is highly dependent upon the topographical and geomorphological characteristics [90]. Therefore, glacier velocity (movement) plays a significant role in mass distribution received at the accumulation and ablation zone of the respective glaciers. The surface displacement/velocity of the Harmukh glacier was derived through the "Co-registration of the Optically Sensed Images and Correlation" (Cosi-Corr) model, a plug-in for the Environment for Visualizing Images (ENVI) interface. Cosi-Corr is freely available software from the California Institute of Technology and provides an important tool to accurately coregister and orthorectify and correlate optical remotely sensed scenes/images to obtain deformations on the ground surface from multitemporal satellite images. The Cosi-Corr model was initially developed and released by the tectonics observatory at Caltech in 2008 [91] to extract the crustal deformations caused by earthquake movements.
Local displacements among the satellite images (pre-and post-event) with different spatial resolutions can be obtained through the Cosi-Corr software [92] and have been extensively used for velocity estimation in various glaciated regions of the world [93][94][95][96][97]. The method coregisters optical satellite images and correlates them to determine the resulting displacement. The present study used two Landsat 8 scenes (pan bands) having 15 m spatial resolution for two consecutive years, i.e., 2020 and 2021, to obtain the glacier velocity. The surface displacement images were thoroughly checked and corrected for possible errors before the glacier surface velocity estimation. Three bands-east-west displacement, north-south displacement, and signal-to-noise ratio-are provided by the Cosi-Corr generated output. The correlation of two Landsat scenes (pre-event and post-event) was carried out using 32 × 32 initial and 16 × 16 final window size with 2 XY steps. The mask threshold was set at 0.90, and robust iteration 2 was used, as suggested by [98]. The results of the correlation yield a displacement image of east-west and north-south directions together with a signal-to-noise ratio image (SNR). The image with north-south and east-west directions was used to calculate the resultant surface movement with an SNR value greater than 0.9. The vector file showing direction was also generated from the discarded image with an SNR value of less than 0.9. The errors related to the calculated surface velocity of the glacier were assessed over the stable patches surrounding the glacier body by presuming less or no movement over the stable areas of a glacier [98,99]. The mean and standard deviation of displacements were summed across stable areas of a glacier and then divided by the interval of time between the two correlated satellite images. The final accuracy of the results generated from the model generally depends upon the exact coregistration, orthorectification, less cloud cover, and shadow of the satellite scenes.

Glacier Mass Balance
Understanding glacier mass balance is important to decode the response of glaciers to climate change and how glacier changes affect water supplies and sea level rise [99]. The glacier mass balance can be estimated through several techniques, such as Accumulation Area Ratio (AAR), Improved Accumulation Area Ratio (IAAR), Geodetic, energy balance, and glaciological methods. Among all these techniques, a glaciological method for mass balance estimation is the most recommended owing to its accuracy, reliability, and insight into the field conditions [100]. However, the glaciological method was not possible for the present glacier due to its remote location and financial and instrumental constraints. Therefore, in this study we used the geodetic balance assessment technique as it is considered a much simpler procedure and quick method to analyze glacier mass changes through time [101][102][103][104]. The geodetic technique quantifies the ice thickness changes by differencing the DEMs of two-time points. The derived thickness changes can be subsequently converted to mass loss by employing the glacier area and density assumptions [104,105]. The geodetic mass changes of Harmukh glacier were estimated using SRTM-DEM for the years 2000 and 2015. The inaccurate or no coregistration of the DEMs may result in enormous mass change values that lead to uncertainties in the final results. Therefore, coregistration of the DEMs was performed in the ERDAS imagine software before the DEM differencing to lessen the uncertainties. Furthermore, the interquartile range method recommended by Huber et al. [106] was used to find and eliminate the outliers in the off-glacier and on-glacier elevations. The widely applied deterministic interpolation technique, i.e., the Inverse Distance Weighted (IDW) method, was utilized for the voided areas derived from the DEM differencing results.
The errors in the glacier mass changes resulted due to the geodetic mass balance method (∇M) associated with the glacier area (A) and the density of ice (ρ), which is assumed to be 850 ± 60 kg m −3 , together with the surface elevation change h, as stated in [102]. Uncertainty associated with the elevation changes (σ ∆H ) was estimated by Equation (3).
where σ 2 dem in the equation represents the uncertainty in the DEM and σ 2 Void f ill denotes the uncertainty that resulted during the void-filling process. The uncertainty related to the date is represented in the equation by σ 2 TDXdate , which was calculated as ±two times the elevation change per year.

Glacial Lake Mapping
Several automated and semiautomated methods and techniques are available for extracting and mapping glacial lakes using geospatial technology, e.g., supervised and unsupervised classification, band rationing technique, NDWI, etc., as used by previous studies in different parts of the Himalayas [13,107]. NDWI, combined with visual interpretation, seems to be the most appropriate technique for identifying and mapping glacial lakes. We used the NDWI technique coupled with manual delineation, as suggested by Huggel et al. [108]. Shadow, cloud cover, and water-rich soils are the main challenges faced during the creation of NDWI maps. Meanwhile, DEM-derived hill shade, aspect, and slope values were removed from the final NDWI map to overcome such problems. High-resolution Google Earth imageries were used for smooth visual interpretation. The uncertainty in the present study was computed for the glacier and lake according to their areal extents following similar methods applied by [109,110]. The resolution of the satellite imagery determines the accuracy of the glacier and lake delineation. High-resolution images cause fewer uncertainties in comparison to low-resolution images.

Glacial Lake Depth, Volume, and Peak Discharge Estimation
The surface area of a glacier lake becomes a key indicator for measuring lake volume and hazard because the depth of a glacial lake is difficult to determine using remote sensing-based methods. In such cases, empirical equations developed in the various glacial lake-bearing river catchments are frequently used. Estimating parameters such as lake depth, volume, and potential peak discharge is essential before assessing and simulating a GLOF event. There are various equations available in the literature to estimate such parameters. Here, we used the widely employed empirical formula given by Huggel et al. [108] to calculate the lake depth and were subsequently compared with other formulae given by Fujita et al. [37] and Patel et al. [111]. In contrast, volume was calculated by several formulas generated for different regions of the Himalayas, such as for the Koshi basin of Nepal Himalayas by Liu et al. [112], the Sikkim Himalayas by Sharma et al. 2018, the Himalayan region by Fujita et al. [37], and the Chandra basin in the western Himalayas by Patel et al. [111]. Additionally, empirical equations developed for regions outside the Himalayas were also included, such as for the Swiss Alps by Huggel et al. [108] and the Cordillera of western Canada by Evans [113]. The results obtained by employing all these equations showed considerable variation. Thus, an average of all these equations was used in the present study, as suggested by Emmer and Vilmek [114] and Chowdhury et al. [115]. Similarly, average peak discharge was calculated using three empirical equations given by Evans [113], Popov [116], and Huggel et al. [108]. The details of the equations used to calculate the depth, volume, and peak discharge are given in Table 2.

S. No. Empirical Equation Source
Volume Hydraulic modeling is an essential part of developing a reliable flood forecasting system. The results of hydraulic model simulations may be used to create inundation maps that can be utilized by community officials or the general public to assess their flood risk. GLOF hazard and risk modeling has been carried out across the world using a variety of methods. For instance, few studies have employed geometric models such as modified single flow [118,119], Monte Carlo Least Cost Path [120,121], and Random Walk Process [122] to analyze and estimate probable GLOF flood inundation. Such types of models need minimal input data and may be applied to a wide range of lakes, but they can only provide a first-order evaluation and cannot provide accurate flood maps [123]. However, more complex modeling software such as HEC-RAS 1D, HEC-RAS-2D, FLO-2D, RAMS, SOBEK, DAM BREAK, and BASEMENT has also been used in different parts of the world [38,50,101,[124][125][126][127][128]. These models can be used to produce flood hydrographs and other important flood wave data, such as water depth, flow velocity, and water surface height, which can be used to develop risk management and mitigation plans.
The HEC-RAS software is an integrated and open-source software developed in a multitasking environment [54]. It has been widely used to study glacial lake hazards [129]. It allows for performing complex one-dimensional (steady flow) and two-dimensional (unsteady flow) calculations, water temperature, water quality modeling, and sediment transport computations. Graphical User Interface (GUI), data management and storage system, separate analysis components, graphics, mapping, and reporting tools are essential tools embedded within the system. It is user-friendly software designed to perform hydraulic calculations for natural and artificial channels, floodplain areas, over banks, and leveeprotected areas. The model is also utilized to simulate flood scenarios resulting from glacial lake breach outburst and is based on one-dimensional St. Venant equations. In the present study, we used the 2D hydrodynamic model, i.e., Hydrologic Engineering Center's River Analysis System (HEC-RAS) version 5.0.7, to simulate the GLOF event in the Sind catchment of the Jhelum basin. In this study, two hydrodynamic modeling GLOF scenarios, i.e., 100% and 75%, were performed for Gangabal lake based on the volume calculated through the empirical equations (see Section 4.1) with varied breach formation time and breach width. In the first stage, 1D HEC-RAS was used to create a breach hydrograph at the lake site. The essential breach parameters that are used as input to calculate the breach hydrograph include breach formation time, breach width, breach height, and volume of the lake. Breach width and breach formation time were derived through Froehlich's empirical equations: where B w and V w are the breach width (m) and volume of the lake (m 3 ), whereas h b and T f are the breach height (m) and breach formation time (h). Equation (1), mentioned in the above section, was used to compute the volume of the lake. B w is the function V w and h b . These equations have been widely used in past studies to approximate such parameters [104,130]. The volume of the lake was estimated through the improved volume area equation given by Qi et al. 2022. The breach heights (h b ) were estimated through the ALOS-DEM [104]. In the present study, the first potential GLOF scenario was considered a worst-case event wherein 100% volume (73 × 10 6 m 3 ) of water was considered to be released from the lake with breach formation time (bf) of 40 min and breach width (bw) 80 m. In the second scenario, 75% of the lake volume was considered for the potential GLOF routing with a breach formation time (bf) of 30 min and breach width (bw) of 40 m (Table 3). Subsequently, generated breach hydrographs were routed through a river starting at the lake site to the last selected location, namely Drag Tanga village, located at a distance of 29.4 km. Unsteady flow analysis was performed for Gangabal lake to observe the peak flow propagation, extent of flood inundation, water surface elevation, and flood water depth. Flow and geometry data are the two important datasets required to run the unsteady flow analysis in the HEC-RAS environment. The complete details of the software are given in the web link https://www.hec.usace.army.mil/. We used an implicit finite volume approach to simulate the flood as an unsteady clear water flow. In actuality, sediment may be entertained in the water flow, altering its rheology and reducing the traveled distance. As a result, assuming a clean water flow reduces the validity of our findings, which is a typical issue in GLOF research throughout the world. However, since modeling erosion, deposition, and transport requires comprehensive field-based data, which is not accessible, this simplification must be accepted [131]. The drainage basin and other elevation parameters were derived from the DEM. Geometry data such as centerline (river), bank lines, flow paths, storage area, and cross sections were created in the HEC-GeoRAS interface and were subsequently imported to the HEC-RAS model. The cross-section profile, XYZ perspective plot, and river bed profile downstream from Gangabal lake are depicted in Figure 6.
The unsteady flow analysis requires upstream and downstream boundary conditions to run the model. The inflow hydrograph is used as an upstream boundary condition, whereas normal depth derived from the river bed profile is taken as a downstream boundary condition. The downstream GLOF routing and the assessment of the flood hydraulics via the channel were done by using the outflow hydrograph (model output) as a boundary condition. Subsequently, 2D hydrodynamic modeling was carried out across the river channel. The 2D flood routing mainly requires the DEM as primary input data. The other data sets include cross sections (500 m on both sides of the channel), flow lines, meanings coefficient, normal depth, etc. Chow [132] suggested a Manning coefficient ranging between 0.03 and 0.07 for the hilly areas having steep slopes with less or no vegetation cover, gravels, boulders, cobbles, and bushes present along the river banks. Considering this, we used a Manning coefficient range of 0.04 to 0.06 for the GLOF routing in the current study.
The computational time was set as 24 h so that the flow wave (modeled) has enough time to propagate downstream to the last selected station. The model provides output in the form of flood inundations, peak flood, flood heights, maximum flood depths, and velocity and routed flood hydrographs throughout the flood channel. In the final stage, the potential impacts of the GLOF event were evaluated at various important locations in the downstream region.

Downstream Impacts of GLOFs
Landuse and landcover maps of the downstream area of the lake were generated to analyze the impact of the GLOF event on the various landuse classes, e.g., agriculture, forest, settlements, etc. Sentinel-2A data from the recent year with a fine spatial resolution of 10 m were utilized for this purpose. A maximum likelihood supervised classification algorithm was used to classify the satellite imagery. Supervised classification has been widely used for the classification of satellite imageries across the world [133,134]. Supervised classification is a three-step process that includes training samples, allocation of classes, and testing. ArcGIS 10.3 software was used to carry out this classification, and at the final stage, the results derived from the classification algorithm were tested for accuracy using the kappa coefficient.

Glacier and Lake Changes
Multitemporal satellite data for the years 1972, 1981, 1992, 2001, 2010, and 2020 were used to analyze the changes in the area of Gangabal lake. The data analysis shows that Gangabal lake has progressively increased in the area during the last 48 years. The lake area increased from 1.42 km 2 in 1972 to 1.46 km 2 in 1981, 1.58 km 2 in 1992, 1.61 km 2 in 2001, 1.64 km 2 in 2010, and 1.66 km 2 in 2020. Overall, the lake area experienced a 17% growth in its area between 1972 and 2020 at an annual rate of 0.005 km 2 per year. Previously, Dar et al. [39] have shown a 0.47 km 2 (36%) increase in the area of Gangabal lake between 1980 and 2016. The variation in the results may have been caused by digitization errors and the satellite imagery used. However, in this study, the post-correction (manual corrections) of the vector polygon layer derived from NDWI was performed by more than one interpreter to lessen the uncertainty in the lake area. Furthermore, the lake area for the years 2010 and 2020 was validated with the area of the lake polygon digitized on the high-resolution Google Earth imagery. The lake expansion is depicted in Figure 7 on different Landsat images. The spatiotemporal analysis of the glacier, namely the Harmukh glacier, associated with Gangabal lake, was also carried out for a period of 30 years during 1990-2020. The data analysis indicated that the area of the glacier in 1990 was 3.6 km 2 , which decreased to 2.9 in 2020. The glacier area for the years 1992, 2001, 2010, and 2020 was estimated at 3.6, 3.3, 3.1, and 2.9 km 2 . Deglaciation of 0.7 km 2 (19 ± 5%) since 1990 at the rate of 0.023 km 2 per year was observed ( Figure 8). The glacier also revealed an upward shift in ELA from 35 to 98 m during the last three decades. The glacier snout has also shown a retreating trend, with a 234 m retreat between 1990 and 2020 at 7.8 m per year. In a recent study, an area loss of 17% (0.6 km 2 ) between 1980 and 2018 was also reported by Mir et al. [135]. An area loss of 0.67 km 2 between 1969 and 2020 was also reported in past studies by Ahmad et al. [82], with a length change retreat of 5.12 m per year. Murtaza et al. [136] observed total deglaciation of 1.23 km 2 (16 ± 0.24%) between 1992 and 2018 in the Harmukh range of the northwestern Himalayas. The Harmukh glacier, the largest glacier in the entire range, lost an area of 0.34 km 2 (10 ± 0.88%), with a snout retreat of 187 m observed during the 28-year study period [136]. The other individual glacier studies across the Kashmir Himalayan basin have shown a retreating pattern; for example, Nehnar glacier has shown a retreat of 1.41 km 2 (50.35%) during 55 years [84], and the Kolahoi glacier lost an area of 23% since 1962 [137], Bodpathri glacier lost 0.5 km 2 (19.5%) during the past 56 years Ahmad et al. [82]. Additionally, glacier area loss, snout retreat, and length changes have been reported in various glaciers of the central and eastern Himalayas.

Geodetic Mass Balance
Glaciological (direct) and geodetic (indirect) techniques are two important methods to estimate the mass balance of a glacier. In the present study, the geodetic method was employed to determine the mass balance of the Harmukh glacier, as direct field-based measurements are not feasible owing to difficult and unapproachable terrain. The results derived from the DEM differencing show an average ice thickness loss of 11.04 ± 4.8 m for the Harmukh glacier at the rate of 0.92 ± 0.40 m/yr between 2000 and 2012. The ice thickness changes over the glacier surface range from 0.31 to 1.02 m/yr. The spatial distribution of elevation changes (ice thickness) is illustrated in Figure 9. The findings reveal negative mass changes over the ablation parts of the glacier and positive changes over the accumulation zone. The findings on average ice thickness loss of the Harmukh glacier are in close agreement with a study carried out by Abdullah et al. [138] over the Upper Indus basin, which reports an average ice thickness loss of 1.12 ± 0.40 m/yr for the Jhelum basin glaciers. However, the average ice thickness changes derived from the DEM differencing technique vary greatly across the different river catchments of the Himalayan region. For example, an average ice thickness loss of 0.46 ± 0.26 m/yr in the Ladakh range, 1.17 ± 0.4 m/yr in the Zanskar range, 1.12 ± 0.40 m/yr in the Greater Himalayan range, 1.69 ± 0.60 m/yr were reported by Abdullah et al. [138]. Another study by Majeed et al. [104] showed an average ice thickness loss of 0.38 m/yr for the Pangong group of glaciers in Trans-Himalayan Ladakh. The accuracy of the glacier thickness changes derived from the DEM differencing technique largely depends upon the accuracy (coregistration) of the utilized DEMs and their spatial resolution. We recommend that with the use of high-resolution DEMs, better results can be inferred through the geodetic mass balance estimation technique.

Glacier Depth
The glacier outlines of "Randolph Glacier Inventory (RGI) version 6.0" were extracted and used to calculate the ice thickness of the Harmukh glacier. The RGI is an open-source database that consists of important parameters (flowlines, elevation, slope, etc.) that are prerequisites for the GlabTop-2 model to compute the glacier ice thickness along with the DEM. The minimum and maximum ice thickness for the Harmukh glacier were estimated to be 7.3 and 85 m, respectively, whereas the mean thickness was calculated as 23.46 m. The spatial distribution of ice thickness was also analyzed by reclassifying the ice thickness map into four classes at equal intervals to estimate the area in each thickness class. As mentioned in Table 4, the maximum area of 47.18% (1.36 km 2 ) of the Harmukh glacier area comes under ice thickness class 1 (0-30 m), 41.85% area is covered with ice thickness class 31-60 m, whereas only 10.95% area of the glacier is covered with the ice thickness greater than 60 m. The maximum ice thickness was observed in the central parts (central flow lines) and near the ELA of the glacier. In contrast, the lowest thickness (pixel value) was observed at the glacier boundaries. The model has been widely used and tested for different glaciers worldwide to calculate the ice thickness and bed over-deepening sites [88,97,102,103,[139][140][141][142][143][144]. The model calibration and validation performed on the Chota Shigri and Bara Shigri glaciers located in Central Himalaya by Sattar et al. [97] revealed a reasonable agreement between GPR data and estimated ice thickness.

Glacier Velocity
A remote sensing-based technique was used in the present study to estimate glacier surface velocity/displacement due to the lack of field-based velocity measurement, which is laborious, time-consuming, and requires financial support. Among various tools and techniques (for example, Im-GRAFT, CIA, and IMCOR) that are used for the estimation of the surface velocity of the glacier, the Cosi-Corr model is appropriate and reliable for the estimation of glacier surface displacement (direction and magnitude) by employing optical images [145], and it has various processing tools that provide additional information [91]. The model yielded a mean velocity of 3.2 m/year using the pre-event (2019) and post-event (2020) Landsat-8 pan bands with a spatial resolution of 15 m. The minimum and maximum surface velocity vary between 0 and 7 m/yr. The velocity was observed to be higher in the upper parts of the ablation zones, varying between 4 and 6 m/yr and then gradually decreasing toward the glacier's snout, as depicted in Figure 10. The same phenomena have been observed in the Dhauliganga basin of central Himalaya by Sattar et al. [97]. The Harmukh glacier is a clean surface glacier; therefore, surface velocity may be greater than other debris-covered glaciers of the Jhelum basin. As illustrated in Figure 10, a significant variation in velocity and flow was observed across the glacier (accumulation zone, ablation zone, terminus). Several studies have carried the glacier surface velocity using satellite images at the catchment level [33,93,94,146] and on individual glaciers across the Himalayan region [147][148][149]. Garg et al. [150] have studied the glacier surface velocity of the Sakchum, Chota Shigri, and Bara Shigri glaciers, revealing an average velocity of 10.2, 20.9, and 25.3 m/yr, respectively, using ASTER images. The field-based measurements of the surface velocity of the Batal glacier, which is in a similar size class to the Harmukh glacier, showed an average velocity of 6.2 m/yr [151]. The surface flow of the glaciers mainly depends upon the ice thickness, surface gradient, and morphology of the ice surface. The surface velocity was higher near the crevasses (slope break), as illustrated in the red and yellow squares shown on the velocity map and validated on the Google Earth imagery. At the same time, the homogenous flow was observed on the gentle slopes.  The depth of the lake was calculated using three empirical equations given in Table 5. The results show a huge variation; therefore, an average of all these equations was taken as the mean depth of Gangabal lake, as suggested by Emmer and Vilmek [114] and Chowdhury et al. [115]. The results further suggest that the depth derived from equations given by Patel et al. [111] and Fujita et al. [37] shows higher values as compared to Huggel et al. [108].

Volume
The volume of the proglacial lake (Gangabal) was derived from five empirical equations developed in different regions of the world. The results show that the volume obtained from the equations given by Evans [113], Huggel et al. [108], and Liu et al. [112] are almost in a similar range. However, volume calculated through the equations of Fujita et al. [37] and Patel et al. [111] shows overestimated values of 102.4 × 10 6 and 117 × 10 6 m 3 . On the other hand, the volume obtained by Sharma et al. 2018 shows a value of 1.54 × 10 6 m 3 , which indicates underestimated values (Table 5). Thus, in this study, an average of all these equations was taken from the volume of the proglacial lake. The volume of Gangabal lake calculated through an average of all these equations was estimated as 73.76 × 10 6 m 3 , which is approximately equal to the volume estimates derived from the recently developed equation given by Qi et al. [117]. Therefore, it can be concluded that an empirical equation for volume estimation given by Qi et al. [117] is suitable for estimating the volume of those lakes that are not accessible for the bathymetry survey.

Peak Discharge
To ascertain the peak discharge of the glacial lake, three equations given in Table 5 were employed. Average volume, i.e., 73.54 × 10 6 m 3 calculated from the equations given in Table 5, was used as input data to obtain the peak discharges from these equations individually. However, like depth and volume, peak discharge values derived from the empirical equations show a huge variation. Therefore an average of all these equations was taken as a maximum peak discharge of the lake. Allen et al. [75] have suggested that the equation of Fujita et al. [37] is more reliable than the equations given by Huggel et al. [108] and Popov [116] because, in the case of the latter, the peak discharge is overestimated. Although we did not consider the equation of Fujita et al. [37] for the peak flow, instead, we used the average of all equations as suggested by Emmer and Vilmek [114] and Chowdhury et al. [115].

GLOF Susceptibility
As a result of the current glacier recession, proglacial lakes are forming and expanding all over the world. Between 1980 and 2018, the number and area of proglacial lakes in the world rose by almost 50% [152]. With the ongoing climate warming, hazards related to glaciers and glacial lakes are expected to increase in the future. GLOFs pose a severe threat to the people living in the mountain regions of the world. Therefore, it has received serious attention in recent years, primarily because of catastrophic damage and fatalities associated with them. The Himalayan region has witnessed several glacier-associated flood events [153] which have flooded the downstream valleys. The valley of Kashmir, where the present lake is located, is no exception. It has also witnessed many flood events dating back to July 1959, when a massive glacial flood hit the valley because of a persistent rain for four consecutive days [154], and a recent being the September 2014 flood event [155]. Because of previous GLOF episodes in Jammu and Kashmir (including Ladakh) and the flood events, it is necessary to monitor glacial lakes associated with glaciers in order to identify lakes that may cause downstream area destruction. According to Mal et al. [43], Jammu and Kashmir have the highest aggregated level of GLOF hazard, followed by Arunachal Pradesh (AP), Sikkim, Himachal Pradesh (HP), and Uttarakhand. A recent study by Dubey and Goyal [44] also categorized JK at the second position in terms of the GLOF hazard level, whereas Sikkim is on the top and followed by HP and AP. Gupta et al. 2021 reported four lakes as potentially dangerous in the Upper Jhelum basin. The present lake has been classified as a potentially dangerous lake by various researchers in the past, for example, a high-risk lake by Gupta et al. [36], potentially dangerous lake by Guru et al. [74], medium-risk lake by Ahmed et al. [41], and important GLOF lakes by Kumar and Goel [42] and Kumar et al. [15].
The lake has an area of 1.66 km 2 according to recent estimates (2020), which is greater than the threshold of 0.1 km 2 reported in past studies for categorizing the potentially dangerous glacial lake [31,[156][157][158][159]. The lake has expanded by 17% during the period 1972 to 2020. The lake is fed by the Harmukh glacier, which has retreated more than 15% during the past 48 years. With the increasing temperature and glacier recession, the area of the lake is anticipated to increase in the coming years. The water supply from Gangabal lake directly flows into Nund-Kol lake, located in the downstream region at a distance of 880 m from the lake outlet. In case of a lake breach outburst event or overtopping failure, the presence of a lake in the downstream course of Gangabal lake may witness a cascading effect, resulting in a catastrophic disaster, as the breaching of one lake can cause the bursting of another. The average slope between the lake and the glacier was estimated at 20.30 degrees (Figure 11a). A steep slope between a lake and its feeding glacier may cause mass movement and ice calving, leading to dam collapse [13,76,78]. The lake is also susceptible to potential rock/ice avalanches and other mass movements owning to its steep slope between the lake and glacier snout (Figure 11b). Rock and ice avalanche zones are visible on high-resolution Google Earth imagery. Glacial lakes with water volumes more than 10 × 10 6 m 3 are regarded as dangerous as they may cause a devastating flash flood downstream. The volume of Gangabal lake estimated using the improved volume estimation equation reveals 73 × 10 6 m 3 of water in the lake. Overall, the lake has been categorized as a high susceptibility lake. A detailed analysis of the GLOF susceptibility of Gangabal lake is available in Ahmed et al. [41].

Hazard Assessment of Gangabal Lake
The occurrence of the GLOF event is governed mainly by two factors (1) dam stability and (2) the possibility of a triggering factor. As far as the dam condition of Gangabal lake is concerned, we have considered it partly unstable in this study and, therefore, susceptible to dam breach/failure. The lake is also susceptible to potential ice/rock avalanches owing to the steep slope between the glacier and the glacial lake, as depicted in Figure 11a,b. Typically, triggering factors such as ice/snow avalanches, rock falls, landslides, or glacier-calving processes produce impulse waves that might start an overtopping failure, eventually resulting in a GLOF event [97]. The other factors that may result in the GLOF event are cloudbursts, a sudden expansion of the lake under warming temperatures, and incessant precipitation, which may cause the lake's water level to rise rapidly, leading to the dam overtopping or a piping failure. The mean water level of Gangabal lake using high-resolution GE imagery and the DEM was measured as 3571 m; the dam is almost flat with its elevation rising from 3576 to 3583.
Using Froehlich's equation, a series of breach hydrographs were calculated for the different breach scenarios with varied breach widths, viz, 20, 40, 60, and 80 m ( Figure 12). All these scenarios produced different peak flood discharge estimates mainly due to the different breach parameters. The peak flood at the time of the lake outburst varied from 10,620.03 m 3 /s (20 m) to 17,787.79 m 3 /s (80 m). In the first scenario (worst-case), 100% volume (73 × 10 6 m 3 ) of the lake was considered to be released for the GLOF routing with a breach formation time (bf) of 40 min and breach width (bw) of 80 m. The resultant breach hydrograph produced a peak flood of 17,787.79 m 3 /s that was achieved within 10 min of initiation of the breach event. The modeled breach hydrographs at various breach formation times are depicted in Figure 12. In the second potential GLOF scenario, 75% of the lake volume (52 × 10 6 m 3 ) was modeled with a 40 m breach width and 30 min of breach formation time (Figure 13). In scenario 2, the peak discharge of 9155.36 m 3 was attained after 12 min after the breach initiation.  Two-dimensional hydraulic routing of generated breach hydrographs using HEC-RAS 5.0.7 was performed for both scenarios (a) scenario 1 (bf = 40 min and bw = 80 m) and (b) scenario 2 (bf = 30 min and bw = 40 m) to assess the potential impacts of GLOFs from Gangabal lake on the downstream valley of the lake up to a distance of 29.3 km. The breach flood, i.e., GLOFs generated by the model, was routed through the river channel, including various locations (villages, bridges, roads, agricultural fields, etc.).
The assessment of the GLOF hydraulics in worst case scenario (scenario 1) revealed that maximum flood depth varies between 3.87 and 68 m, maximum flow velocity between 4 and 75 m/s, and maximum water surface elevation varies between 1548 and 3536 m (Figure 14a-c). The resultant flood potentially inundated areas at nine settlement locations, including 98-104 buildings. It is pertinent to mention that most buildings and other structures were inundated at the end of the channel, which is within the elevation zone of 1726 to 1869 m asl. However, depth, velocity, and peak flow were at their highest in the upstream regions due to the steep slope of the channel walls. The spatial distribution of flood inundation is shown in Figures 15-17.  The peak flood and other important breach flood wave parameters such as flood depth, flood velocity, and WSE derived from the model output at important cross-sections (locations) were found to be maximum at Naranaag temple (XCS-45500), situated at a distance of 11.67 km from the lake site and attained at 1:00 h after the breach initiation. The peak flood was 13,995.36 m 3 s −1 for the same cross-section (Naranaag temple) with a maximum flood depth of 14.43 m and maximum flow velocity of 8.83 m/s. The details of the derived parameters at each important location (cross-section) are provided in Table 6, whereas GLOF hydrographs are provided in Figure 17.  The analysis further revealed that flood depths at various cross-sections shows sudden fluctuations: the flood height increases at 13.37 km from the lake outlet, then suddenly decreases to 9.2 m at a distance of 14.67 km. Again, the flood heights gradually rise to 17.6 m near the Baba Nagre Khanam bridge and then gradually decrease to 3.8 m at the Sind river confluence. A similar phenomenon has been observed in the GLOF-routed channel of Thulagi lake [46,160] and Lower Barun lake [160]. This phenomenon resulted due to the narrow and shallow depth of the routed channel. At the narrow channel locations, the depth, velocity, and peak flood were observed as high, while it decreases at the locations where the width of the river channel is broader, which results in an increase in water spread and, thereby, a decrease in the magnitude of flood wave parameters. Further downstream, the peak discharge of 10,848 m 3 s −1 was estimated at bridge site 5 with 8.98 m flood depth and 6.86 m/s flood velocity, respectively. At the last selected location, i.e., Drug Tanga village, the flood discharge was estimated as 10,223.98 m 3 s −1 , whereas flood height and flood velocity were observed as 9.26 m and 5.37 m/s, respectively, shown in Figure 17.   In the second scenario, the maximum flow depth and velocity were modeled as 48 m and 52 m/s, respectively, whereas the maximum WSE rises to 3546 m at the lake site, which decreases to 1689 m at the Drag Tanga village. The peak flood of 5231 m 3 s −1 was estimated at the Naranaag temple cross-section with a peak depth and velocity of 9 m and 5.2 m/s, respectively. The total number of inundated buildings in the second scenario decreased to 75-80. For example, in the case of scenario 1, a total number of nine buildings were inundated near the unnamed settlement, as depicted in Figure 18a. However, in the second scenario, only four buildings were potentially inundated due to the resultant GLOF (Figure 18b). In the first scenario, the total inundated area was estimated as 5.31 km 2 , which also decreased to 3.46 km 2 in scenario 2. The inundation map derived from scenario 2 is depicted in Figure 18. Flow depth, velocity, and peak discharge decrease as the routed flow reaches the downstream areas having lower elevations. The analysis revealed that the breach hydrograph is the primary factor that governs the depth, velocity, WSE, and peak flow in the downstream channel. The maximum flow depth and velocity decrease as the Manning coefficient is increased. Flow velocity was observed higher at the places with steep slopes because it is dependent upon the slope of the channel. Depth was found highest in the upstream locations compared to the downstream areas; in the former case, the width of the channel is less, while it is more in the latter case. The channel widens as it proceeds toward the low-elevation regions, thus providing a sufficient water spread area.
The maximum area is inundated in the downstream region due to the gentle slope and greater channel width, but values of flood wave parameters remain minimum in such areas. The investigation showed that the flow hydraulics of a routed GLOF is determined by the initial breach hydrograph. The initial breach hydrograph significantly impacts the hydraulic characteristics of a GLOF at any point downstream. The flow hydraulics were also examined in relation to the channel's top width. The top width of the channel increases in a similar manner to the flow velocity, decreasing both inundation depth and flow velocity [161].

Downstream Impact Assessment
The threat of GLOFs is increasing not only in the hilly terrains but also in the downstream valley regions. Therefore, a thorough evaluation of the GLOF hazard, risk, and susceptibility using reliable scientific methodologies is the current need (Khanal et al. 2015). In this study, after evaluating the GLOF impacts on the downstream course of the river and its adjoining areas, the area under each inundated landuse class was estimated through the LULC map generated for both scenarios (scenario 1 and 2). The LULC map prepared for scenarios 1 and 2 is shown in Figure 19a,b. The analysis reveals that in the worst-case scenario, 5.3 km 2 of the geographical area would be potentially inundated due to the outburst flood from the lake outlet to the last selected location (Drag Tanga village).
The inundated landuse classes in scenario 1 in terms of area and percentage include water area 0.8 km 2 (15.09%), forestland 0.7 km 2 (13.2%), grassland 0.4 km 2 (7.55%), cropland 0.6 km 2 (11.32%), scrub/shrub 1.1 km 2 (20.75%), built-up area 0.3 km 2 (5.66%), and barren 1.4 km 2 (26.42%). In the second scenario, wherein 75% of the lake water volume was considered for the GLOF routing, the inundated area decreased by 1.84 km 2 which is 34% of the total inundated area in scenario 1 ( Table 7). The total area potentially inundated in scenario 2 was estimated as 3.46 km 2 , which includes a water area of 0.58 km 2 (16.76%), forestland of 0.51 km 2 (14.73%), grassland of 0.43 km 2 (12.42%), cropland 0.35 km 2 (10.11%), scrub/shrub 0.62 km 2 (17.91%), built-up area 0.19 km 2 (5.49%), and barren 0.78 km 2 (22.54%). The investigation shows that the exposure level was greater in the low-lying areas compared to the areas situated at higher elevations. Validation and calibration are crucial parameters in flood inundation models to evaluate and match inundation extents and depths. However, owing to the lack of field data and observed flow data, it is often difficult to validate the hydraulic properties derived from the flooding models. The derived hydraulic calculations cannot be validated in the present study because no past GLOF event has been witnessed or reported for Gangabal lake, although the model reliability has been verified with various studies that have documented GLOF events in various parts of the world, such as Peru, China, India, Nepal, and Tibet.  Glacial lakes are susceptible to breach outbursts under the conditions of rapid warming and incessant rainfall (extreme weather events). GLOFs may also occur due to the sudden ice/rock avalanche, landslides, cloud bursts, glacier surges, ice core melting, and glacier snout calving [33]. The temperatures in the entire Indian Himalayan region are significantly increasing, as evident from recent studies [45,162,163]. Averaged temperature mean, minimum, and maximum over India are showing a warming trend with an increase of 0.15, 0.15, and 0.13 • C per decade, respectively, between 1986 and 2015 [164]. The changing precipitation regimes and increasing temperatures over the Indian Himalayan region are major concerns for the region's snow cover and glacier health [10]. The temperatures in the Jhelum basin, where Gangabal lake and its feeding glacier are located, have shown a significant increasing trend and decreasing precipitation trend [84,165,166], while glaciers, on the other hand, are rapidly retreating [41,84,166], and glacial lakes are expanding, and their frequency is increasing [39][40][41]. The warming temperatures and, thereby, rapid glacier melting in the area may result in sudden glacial lake outburst events, enhancing the GLOF susceptibility of the region.

Limitations
The HEC-RAS model was successfully applied in this study and in various past studies to obtain the various hydraulic properties (flood depth, flood velocity, water surface profiles) of the outburst flood and aid in the development of inundation maps. Evaluating the accuracy of the outputs derived from the physically based hydraulic models is difficult since they are sensitive to numerous input parameters [125]. The topographic and geomorphological characteristics of the routed channel also affect the nature of GLOFs. The peak flood discharge and flood velocity through a particular channel are directly affected by Manning's roughness coefficient, which substantially affects GLOFs. The increase or decrease in Manning's "n" values plays a significant role in the peak discharge variations. GLOF modeling through hydrological models is nothing more than an approximation of an event. Despite that, these models can be utilized to analyze a physical phenomenon and its corresponding effects on the development of water resource structures and flood management. Most assumptions in GLOF modeling are based on breach parameters (breach depth and breach width), which impact peak flood discharges and the time of flood arrival. The results contain errors and uncertainties since there is a paucity of field-based bathymetry data and a number of other parameters related to the dam break and the routing of the breach hydrograph.
The use of a DEM with coarse spatial resolution also overestimates or underestimates the flood inundation and the magnitude of the hydraulic parameters. The DEM simply defines the elevation of a point in a given area at a particular spatial resolution. Based on the variation observed in the model output, we suggest considering uncertainties associated with the DEM when evaluating the GLOF flow routing and its potential socioeconomic impacts. The accuracy of the water depth is also influenced by the DEM data used to extract the cross-sections; if the DEM data are very accurate, the error will be minimal. The surface area of the glacial lake, considered a critical input parameter to generate a breach hydrograph, was obtained through the empirical equation. The empirical formulae used to calculate the volume and depth of Gangabal lake are based on the regression function of its area. The depth and volume estimates derived from the empirical equations for Gangabal lake may provide rough estimates and overestimates the values due to its larger area. However, when field-based data are not available or possible to obtain, empirical equations can be used to gather first-hand information. The model's main limitation lies in its taking into account the clear water while simulating unsteady flow. In actuality, sediment may be entrained in the flow, altering its rheology and reducing the distance traveled. As a result, assuming a clean water flow reduces the validity of our findings, which is a typical issue in GLOF research worldwide. The mechanism of glacial lake bursting and the genesis of glacial lake breaches are not well understood. Despite the various drawbacks mentioned above and assumptions, hydrodynamic modeling is useful because it provides some valuable estimates of the glacial lake breach outburst flood, enables an accurate estimation of the design flood, and offers some guidance for the prevention and mitigation of the potential flood.

Conclusions
Generally, glaciers across the Himalayas are retreating, resulting in the formation and growth of numerous glacial lakes under the direct influence of climate change. The glacial lakes are expanding at a remarkable pace. In this study, a detailed study of the potentially dangerous glacial lake, Gangabal, and its feeding glacier, Harmukh, was conducted using the integrated approaches of remote sensing, GIS, and hydraulic modeling using HEC-RAS. Our data analysis reveals that the glacial lake has gradually expanded from 1.42 km 2 in 1972 to 1.46 km 2 in 1981, 1.58 km 2 in 1992, 1.61 km 2 in 2001, 1.64 km 2 in 2010, and 1.66 km 2 in 2020. Overall, the lake area experienced a 17% growth in its area from 1972 to 2020 at the rate of 0.005 km 2 per year. On the other hand, the glacier has shown a substantial decrease in its area, from 3.6 km 2 in 1990 to 3.3 km 2 in 2000, 3.1 km 2 in 2010, and 2.9 km 2 in 2020. The region observed a deglaciation of 0.7 km 2 (19 ± 5%) since 1990 at the rate of 0.023 km 2 per year. The results derived from the DEM differencing show an average ice thickness loss of 11.04 ± 4.8 m for the Harmukh glacier at the rate of 0.92 ± 0.40 m/yr during the years 2000 and 2012. The minimum and maximum ice thickness for the Harmukh glacier derived from the GlabTop-2 model was estimated to be 7.3 and 85 m, respectively, whereas the mean thickness was calculated as 23.46 m. The Cosi-Corr model yielded a mean velocity of 3.2 m/yr using the pre-event (2019) and post-event (2020) pan bands with a spatial resolution of 15 m from Landsat-8. The minimum and maximum surface velocity varied between 0 and 7 m/yr. Assessment of GLOF propagation in the worst-case scenario (scenario 1) revealed that the maximum flood depth varies between 3.87 and 68 m, the maximum flow velocity ranges between 4 and 75 m/s, and the maximum water surface elevation varies between 1548 and 3536 m. The resultant flood wave in the worst-case scenario would reach the nearest location (Naranaag temple) within 1 h 10 min after the initiation of the breach with a maximum flood discharge of 12,896.52 m 3 s −1 and maximum flood depth and velocity of 10.54 m and 10.05 m/s, respectively. After the GLOF impact assessment on downstream exposed areas, the area under each inundated landuse class was estimated through the LULC map generated for both scenarios 1 and 2. In the first scenario, the total potentially inundated area was estimated as 5.3 km 2 , which was reduced to 3.46 km 2 in scenario 2. Owing to the tremendous hazard potential and incidence of GLOFs in the Himalayan region, it is crucial to consider the probable impacts of GLOFs while planning, designing, and building key infrastructure along the flow channels. This is especially important when constructing dams, as they are located on the flow path of a river channel. We attempted a preliminary assessment of both the Harmukh glacier and Gangabal glacial lake by applying an integrated research approach. We suggest utilizing high-resolution satellite data, climate data, field investigations, and robust hydrological models to comprehensively investigate benchmark glaciers and dangerous glacial lakes in the region. Data Availability Statement: Data will be available from the corresponding author on request.