Dynamics of Lower-Barun Glacier and Glacial Lake and its GLOF Susceptibility Using Geospatial Analysis and Modelling

Shrinkage of some of the glaciers has direct impacts on the formation and expansion of glacial lakes. Sudden glacial lake outburst floods (GLOFs) are a major threat to lives and livelihoods downstream as they can cause catastrophic damage. In this study, we present the dynamics of the Lower-Barun glacier and glacial lakes and their GLOF susceptibility. We used multitemporal Landsat and Sentinel satellite imagery and extracted the lake outlines using the Normalized Difference Water Index (NDWI) with manual post-correction while the glacier outline was digitized manually. Multi-criteria decision-based method was used to assess the GLOF susceptibility. For the estimation of peak discharge and failure time, an empirical model developed by Froelich (1995) was used. The surface area of the Lower-Barun glacial lake was increased by 86% in the last 40 yrs (from 1979 to 2018), with a mean increase of 0.0432 km2/yr. The shrinkage in the glacier area Journal home page: https://cdhmtu.edu.np/journal/index.php/jalawaayu hnjfo' JALAWAAYU Research Article (Interdisciplinary Journal of Atmospheric and Hydrospheric Sciences)


Introduction
Himalaya is the hotspot for the effects of climate change and slight change in climate possess various effects on the Himalaya and its ecosystem. High Himalaya glaciers are retreating rapidly, often resulting in formation and expansion of glacial lakes, which later can cause glacial lake outburst floods (GLOFs) (Bolch et al. 2011). Over the past 30-year period, the glaciers had lost almost a quarter of their total area in Hindu-Kush Himalaya (HKH)region (Bajracharya et al. 2015). In Nepal Himalaya, the total glacial area had decreased by 24% between 1977 and 2010, while the number of glacial lakes increased by 11% (Bajracharya et al. 2014b). Before the 1950s, most presentday moraine-dammed lakes did not exist in Nepal Himalaya (Chalise et al. 2006). During the mid-1950s to 1960s, small supraglacial lakes started to coalesced and began to grow rapidly in the 1970s to form a present-day large moraine-dammed glacial lake (Watanabe et al. 2009;Bajracharya & Mool, 2009;ICIMOD, 2011). In the Himalayan region, glacial lakes are increasing in number and volume due to thinning and recession of the glaciers (Khadka et al. 2018;Mool et al. 2001).
The glacial mass loss rate in High Mountain Asia (HMA) and global average mass loss are comparable and suggest that glacial loss is prominent all around the world in this century (Cogley, 2012). A study suggested that HMA will lose around 65% of its ice mass this century as per the climate model projections driven by high emissions scenarios (Kraaijenbrink et al. 2017) that could result in the formations of new glacial lakes.
Recent studies report that the Koshi, Gandaki, and Karnali basin of Nepal, the Tibet Autonomous Region (TAR) of China, and India has 3,624 glacial lakes in total, of which 2,070 lakes are in Nepal, 1,509 lakes in the TAR, China, and 45 lakes in India. A total of 47 glacial lakes are identified as potentially dangerous glacial lakes (PDGLs). 21 out of 47 PDGLs are situated in Nepal, which is the same number as the previous number identified in 2011(ICIMOD, 2011Bajracharya et al. 2020). However, only 13 of them are the same as previous, and 8 of the previous PDGLs have been removed from the list, and new 8 PDGLs are added (Bajracharya et al. 2020). Likewise, a recent study on Mahalangur Himalaya displays 345 glacial lakes (>0.001 km 2 ) in 2018. Furthermore, the lakes with an area of more than 0.045 km 2 were assessed; seven of those lakes were identified with very high GLOF susceptibility (Khadka et al. 2020).
GLOFs refers to the sudden release of stored lake water, which can have devastating socioeconomic consequences, including loss of life, buildings, bridges, transportation routes, arable land, and hydropower systems. Countries and people in the Himalayan region face GLOFs as one of the crucial problems (Ageta et al. 2000). Himalayan GLOFs develop at high altitudes, outspread for long distances, and cause monstrous damage to downstream infrastructure and livelihood (Chen et al. 2007;Liu et al. 2014). Nepal had experienced 26 recorded GLOFs events that caused noticeable damage and with a reported loss of life. The GLOFs that occurred at Nare Lake in 1977, Dig Tsho in 1985, Tam Pokhari in 1988, and Chubung lake in 1991 are examples of the devastating socioeconomic impacts that can be caused by GLOFs, which include loss of life in addition to the damage to livelihoods, hydropower systems, and infrastructure (Bajracharya & Mool 2009;Bajracharya et al. 2008;Dwivedi et al. 2000;Vuichard & Zimmermann, 1987).
In the Himalaya, the main trigger of GLOFs is an avalanche entering a lake, which may cause a tsunami-like wave that can overtop and erode the terminal moraine and generate the ensuing flood (Richardson & Reynolds, 2000). Furthermore, the factors like glacier microclimate, glacier dynamics, topography, glacier debris cover, proximity to downstream infrastructure, potential flood volume, and peak discharge intend to variability to risks and lake growth rates (Benn et al. 2012). The high-relief areas around Mount Everest and Makalu Range in Nepal have a high risk of exposure to glacial lakes. These regions have many lakes and experience frequent slope failure and flash floods due to monsoon (Haritashya et al. 2018). The Himalaya have been reported as the hotspot of GLOFs. More than 50 recorded GLOF events in Himalaya, and only a few have detailed reports (Nie et al. 2018).
The Barun valley has experienced at least two known GLOFs (USAID, 2014;Byers et al. 2019;Sattar et al. 2021). The first GLOF event was in 1964, not recorded but witnessed by local people (Yamada and Sharma, 1993). The second GLOF was triggered on April 17, 2017, due to the rockfall and avalanche in Langmale lake (Byers et al. 2019). The majority of the GLOF events reported in the Himalaya were caused by snow/ice/ rock avalanches (Nie et al. 2018). There are several GLOF susceptibility assessment methods, and they vary based on the parameters considered input data and the level of subjectivity (Prakash and Nagarajan, 2017). For instance, Methods for evaluating the GLOF susceptibility for the Asian high mountains are developed by Bolch et al. (2011), for the Himalayan region (Wang et al. 2011). Several empirical models for the estimation of peak discharge and failure time are available. Some of the models are discussed in Singh and Snorrason (1984), MacDonald and Langridge-Monopolis (1984), US Bureau of Reclamation (1988), Van Thun andGillette (1990), andFroehlich (1995).
A recent study found that Nepal is one of two countries with the most significant socioeconomic consequences from glacier outburst floods in the world (Carrivick & Tweed, 2016). Engineering interventions have been applied in the glacial lakes of Peru to reduce the potential risk of GLOF, but such efforts are limited in HMA (Haritashya et al. 2018;Cuellar & McKinney, 2017). Imja (ICIMOD, 2011;Rana et al. 2000) and Tsho Rolpa glacial lakes have been in the spotlight. The GLOF mitigation work has been successfully installed in Nepal (Somos-Valenzuela et al. 2015;Cuellar & McKinney, 2017). Furthermore, an early warning system has been implemented to warn communities downstream about GLOF. This study focuses on assessing the dynamic of Lower-Barun glacier and glacial lake and GLOF susceptibility using geospatial tools and modelling.

Study area
The Lower-Barun Glacial Lake falls under the Mahalangur Himalaya. It is located in eastern Nepal at 27˚47'51" N latitude and 87˚05'26" E longitude at 4450 m above sea level (Haritashya et al. 2018;Maskey et al. 2020). It is located southeast of Imja Glacial Lake in the upper reach of the Barun River and adjacent to the Makalu massif. Lower-Barun Glacier is approximately 2.7 km long from the calving front to the base of an icefall. The lake has consistent calving of the terminus glacier. The discharge from the lake feeds the Barun river and finally to Arun River.

Spatial data
Multi-temporal (extending from the 1979s to recent time) and multi-spectral (blue, green, red, near-infrared, and shortwave channels) optical satellite imagery was used for the lake evolution and glacier change analysis for the entire MBNP. The opensource data, such as Landsat TM, ETM+, and recent Landsat-8 OLI; Sentinel 2A/2B, 5 m ALOS digital elevation model (DEM) was used to obtain the basin characteristics, analyzing glacial lakes and glacier changes.

Glacial Lake Mapping
The multi-stage process with a semi-automatic classification method was used to delineate glacier and lake boundaries (Bajracharya et al. 2014). Many studies have used manual digitization (Salerno et al. 2012;Wang et al. 2015;Zhang et al. 2015) and the Normalized Difference Water Index (NDWI) based fully automated or semiautomated methods the delineation of glacial lake outlines. The NDWI is the most widely used method to delineate the water bodies from satellite imagery. We used semiautomated NDWI (Eq. 1) to delineate glacial lake boundary with manual post-correction by visually checking (Huggel et al. 2002).
Where, ρ G, and ρ NIR represent the top-of-atmosphere (TOA) reflectance for the green, red, and Near Infrared (NIR) bands, respectively.

Bathymetry survey and field data collection
A detailed field survey was conducted in 2020. Lower-Barun Glacial Lake volume, water depth, dam height, cross-sections of the lakes, and moraine characteristics was collected from the field. The topographic and morphological characteristics of lakeshore (including mountain slope, inlets, outlets, land use) were collected.
For the bathymetry survey Humminbird Depth Sounder HDR 650 and Rafting boat were used. The survey was conducted along transverse and longitudinal transect lines as per the designed protocol. However, because of bad weather conditions and wind, the prior lines were deviated in real-time surveying. Bathymetric models were created within ArcGIS 10.3 using natural neighbor interpolation (Haritashya et al. 2018;Thompson et al. 2016), and lake volume was calculated. Bathymetric modeling requires an outline with zero depth. The outline zero-depth data were generated lake extents using multispectral satellite imagery.

GLOF Susceptibility
Multi-criteria decision-based method is used to assess the GLOF susceptibility of the Lower-Barun Glacial Lake. Various potential factors that govern the mechanism of potential GLOF were identified and compiled. Six factors, lake area, expansion rate, glacier-lake distance, dam front slope gradient, ice/snow avalanche potential, and the potential of rockfall with upstream GLOF were selected to assess the GLOF susceptibility of the Lower-Barun Lake (Prakash and Nagarajan, 2017;Khadka et al. 2020).
The lake area, expansion rate, and glacier, and lake proximity are calculated using the multi-temporal and satellite image and the DEM analysis. The dam front slope gradient (DFSGS) is calculated by using the difference between the elevation of a dam, elevation of the base, and horizontal length (Eq. 2) (Fujita et al. 2013).
The avalanche is the main trigger of the GLOFs in the Himalaya (Nie et al. 2018). So, this factor was ranked as the highest priority by expert judgment, and maximum weight is given to this factor in Table.

GLOF Susceptibility Index
The weights of each factor are assigned using the Analytical Hierarchical Process and developed by Khadka et al. (2021). Three classes are provided with index values of 1, 0.5, and 0.25 (Khadka et al. 2020;Prakash and Nagarajan, 2017). The final weights for each factor were calculated by multiplying the factor weight (W i ) with the class index values (C i ) based on the measurements of lake factors. The final weights of each factor were summed to evaluate the GLOF susceptibility of the lake (Eq. 3).

Moraine Dam Breach Modelling
Lower-Barun Glacial Lake is one of the six potentially critical glacial lakes of Nepal Himalaya, identified by the Bajracharya et al. (2020). The field survey data of the lake will support to assess the potential GLOF risks. The potential moraine dam breach hnjfo' JALAWAAYU Volume 1, Issue 2, 2021 modeling performed for Lower-Barun Glacial Lake (Aggarwal et al. 2016;Kougkoulos et al. 2018) using different dam break parameters collected from the field, satellite data, and digital elevation model. The morphometric characteristics of the lake and failure time of the breach are estimated from empirical equations and field observation. Several empirical dam breach models are available in the literature (Westoby et al. 2014;Wahl, 2014). In this study, the equations developed by Froehlich (1995) was used since this method has the lowest uncertainty (Wahl, 2004;Basheer et al. 2017).
The peak discharge (Q p , m 3 /s) and dam failure time (T p , h) were estimated as follows (Froehlich, 1995):  Lake surface area km 2 2 Volume of the lake m 3 3 Breach height m 4 Breach width m

Lake Dynamics
The evolution of the Lower-Barun Glacial Lake from a small supraglacial lake to a substantial proglacial lake was assessed from 1979 to 2018 using Landsat and Sentinel imagery ( Figure 2). Our analysis shows the significant growth in the area of the lake over the last four decades. The lake area increased from 0.305 km 2 in 1979 to 2.019 km 2 in 2018 (Table 4). The average percentage change in the lake area per decade is 38%, and overall, 85% change is observed in 40 years' time period. An average increase of 0.0432 km 2 area per year is observed.

Glacier Dynamics
The decadal change in the Lower-Barun Glacier is more significant towards the calving fronts (Figure 3). The glacier area decreased from 24.61 km 2 to 22.63 km 2 over the past four decades (1997 to 2018) ( Table 5). The shrinkage in the glacier area is around 0.49 km 2 per decade in average for the last four decades. The glacier has shrunk 8.02% in last 40 years' time period. The retreat of Lower-Barun Glacier is 0.20% per year in the last four decades.  The Lower-Barun Glacier has retreated significantly. The terminus of the glacier has already retreated 2550 m during last four decades ( Figure 5). Around 830, 530, 640, and 720 m of the glacier length retreated in 1979-1990, 1990-1998, 1998-2008, and 2008-2018 Figure 6 exhibits the extensive field-based results showing the Lower-Barun Lake and its surrounding features responsible for GLOF susceptibility. The red area in the figure signifies the possible avalanche area. This avalanche area is justified by Figure  7a, where we observed the live avalanche during the field survey. This part has a high slope and snow cover. The same area is shown in figure 7b. It also shows the area with possible rockfall areas.

GLOF Susceptibility
For the GLOF susceptibility index, the rank of each factor was found to be high. The lake area is more than 0.5 km 2 , the expansion rate is around 0.85 over the last four decades, which is greater than 0.5, the lake is in contact with glaciers, the DFSG is calculated as around 14 degrees, which is more than 10 degrees. Likewise, the lake is susceptible to avalanche and rockfall. The susceptibility index was found to be 0.94, suggesting that the lake is highly susceptible to GLOF.  Figure 7, (a) shows the live avalanche in the Lower-Barun lake, (b) highlights the possible rockfall areas which has the high slope, (c) points out the onsite bathymetric survey, (d) with highlighted portion shows the steep glacier and (e) photo shows the Lower-Barun Lake with its terminus (Figure 7).

GLOF Modelling
The bathymetric survey was conducted in November 2020 at Lower-Barun Glacial Lake. The water volume of the lake is found to be ˳103.6 × 106 m 3 . Due to the weather variations survey near the glacier terminus was not possible, so in such a case, points from Haritashya et al. (2017) were used for verification and authentication. The deepest part of the lake is near the calving front of glacier. The water volume displayed by our result is found to be less than that of volume reported by a study in 2017 (112 * 106 m 3 ) (Haritashya et al. 2017).

Figure 8. Bathymetric map of the Lower-Barun Lake
The GLOF modelling of Lower-Barun Lake was conducted based on three potential moraine breach scenarios. The peak discharge and failure time is calculated using the empirical equation given by Froehlich (1995). The peak discharge of 5768 m 3 /s is produced when the breach depth is 20 m and the entire water volume is released. Likewise, in the case of 15 m breach depth, the peak discharge of 4038 m 3 /s is formed. Breach depth scenario of 10 m, peak discharge of 2442 m 3 /s is produced and in case of breach depth of 5 m produces the peak discharge of 1034 m 3 /s. The peak flow is achieved at 3 hrs in after breach initiation in 20 m breach.

Lake Evolution and Dynamics
The glacial lakes in the Himalaya regions have expanded significantly (Nie et al. 2017). Our results also suggest that the growth of the Lower-Barun Glacial Lake is vibrantly noticeable in the last four decades, from 1990 to 2018. Like most of the glacial lakes in Nepal Himalaya, the evolution of this lake also started from the formation of the small supraglacial lakes, and later such small lakes coalesced to form large proglacial lakes (Haritashya et al. 2017;Maskey et al. 2020). The surface area of Lower-Barun Glacial Lake has been increasing at an alarming pace from 0.04 to 1.8 km 2 , with an average increase of 0.04 km 2 /yr (Haritasjya et al. 2017). Our analysis documents a similar average increase of more than 0.0432 km 2 /yr. A similar result was presented in a recent paper with an average growth of 0.054 ± 0.006 km 2 /yr in Lower-Barun Glacial Lake (Maskey et al. 2020). The lake expansion rate of Lower-Barun Lake was found to be 46.8% which most than that of lakes, like Imja Tsho and Lumding Tsho with expansion rate 42.1% and 32.9%, respectively during 2006 to 2016 (Khadka et al. 2019).
Lower-Barun Lake initially expanded slowly, but later with increased timeline area of the lake also enlarged. Lake expansion is pretty apparent in the future. Figure 6 displays the possible lake extent for the future. A study shows that the more excellent solar absorption by the more giant lakes in the summertime caused thermal erosion with enhanced ice calving and swift retreat of the glacier . Other factors like free convection and subaqueous melting allied with tidal action and upwelling also drive the melt of ice at calving front and calving rate (Rohl, 2008;Benn et al. 2007). The thermal erosion at the ice-contact water body has a vital influence on calving (Benn et al. 2007), and it suggests that at Lower-Barun Lake, it may have caused higher rates of retreat and augmented calving.
Furthermore, the lakes with a more calving ice front are assumed to have a relatively higher likelihood of failure and pose potential hazards downstream (Mergili & Schneider, 2011;Wilson et al. 2018). Our result reports the significant rate of lake expansion over the last four decades. Lake area expansion rate is considered one of the essential factors for potential flood hazards (Khadka et al. 2020).

Lower-Barun Glacier Dynamics
The Lower-Barun glacier shrunk by 8.02% in the last four decades. The glacier area has decreased and has supported the expansion of the Lower-Barun glacial lake. The Lake developed into a sizeable proglacial lake from a small supraglacial lake in the last four decades (Khadka et al. 2019). There is a significant amount of evidence of glacier shrinkage and retreat over past decades in the Nepal Himalayas (Kadota et al. 1997;Fujita et al. 2001;Bolch et al. 2008;Bajracharya et al. 2011). Studies also show that the overall glacier area decreases from the 1980s to 2010 (Bajracharya et al. 2014b). A study on the Everest region has shown that the glaciers are retreating at an average of 10-59 m/yr (Bajracharya & Mool, 2009). Temperature is one of the vital climatic variable responsible for the retreat and meltdown of the glacier (Thakuri et a. 2019) The average warming of 0.06 degrees Celsius per year over 1971-1994 has been reported (Shrestha et al. 1999). Over the last four decades from 1976-2015, on average, maximum air temperature increased +0.045 °C/yr over Nepal's entire country (Thakuri et al. 2019).
The mass budget of the Himalayan glacier has decreased over the last five decades as reported by the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC). This also justifies the fact that the Lower-Barun Glacier has retreated over the last four decades . With the increased retreat of the glacier, it has significantly affected the Lower-Barun glacial lake. The lake area has increased by around 85% over the last four decades, and the lake is more susceptible to potential GLOF hazards.

GLOF Susceptibility
Recent studies also report that the Lower-Barun Glacial Lake is very highly susceptible to GLOF and Rank I (Bajracharya et al. 2020;Haritashya et al. 2017;Khadka et al. 2020). The multi-criteria-based assessment was used to classify the glacial lake GLOF susceptibility. This process is systematic, easy to apply and repeat, and straightforward. However, this method also has certain limitations, such as subjective expert judgment and evaluation.
The live experience of an avalanche in the lake during the field expedition ( Figure 6a) and the high slope rock formations in the mountains exerts high risk to the lake and makes it more susceptible to GLOF. Since the majority of recorded GLOFs (+50) in the Himalaya were triggered by snow/ice/rock avalanches (Nie et al. 2018). However, the chances of an avalanche or rockfall approaching the lake and producing the GLOF depends upon their volume, velocity, trajectory, and positions concerning the lake (Bolch et al. 2011;Rounce et al. 2016). The hanging glaciers can also raise the GLOF susceptibility. In addition, large surge waves produced either by frontal calving or avalanche, if generated, could easily cause over-topping failure.
The Lower-Barun Lake is very highly susceptible to GLOF; it does not necessarily indicate that the lake is currently in a state of outburst; instead, this lake should be emphasized and require regular monitoring and extensive In-Situ investigation. The susceptibility of GLOF should be minimized to ensure the lower impacts downstream in case of GLOF outburst. Nepal Government has already practiced the engineering method, lowering the water level and stabilization of drainage outlet in Imja lake (UNDP Imja Report, 2013; http://cfgorrp.dhm.gov.np/). Installation of the early warning system downstream could help in the minimization of impacts due to GLOF hazards.

GLOF Modelling
The decrease in the water volume is due to the gap in the bathymetric measurements and due to the limitation of the depth sounder. The lake is deepest near the calving area and is slightly shallow near the lake terminus; this feature is signified from the previous study.
This study modeled four potential scenarios, but GLOF could occur in several other ways. The study carried out by Maskey et al. (2020), has used NWS-BREACH model to calculate the peak discharge using different scenarios of 5, 10, 15 and 20m breach. The peak discharge of Lower-Barun Lake ranged from 2,619 to 7,936 m 3 /s. The peak discharge calculated by this study is slightly less than that stated by Maskey et al. (2020) as the water volume of the Lower-Barun glacial lake is less as per this study. Likewise, the failure time obtained by this study for 20 m breach is 3 hrs which significantly matches with the failure time displayed by Maskey et al. (2020). A recent study has used the eight possible scenarios to model the GLOF in Lower-Barun Lake (Sattar et al. 2021) and this study reported that the peak discharge from Lower-Barun GLOFs at current condition could be high as 20,810 m3/s, with flow depths and velocities reaching up to 29 m and 11 m/s, respectively, and can exert disastrous risk downstream (Sattar et al. 2021). The places like Yangla Kharka, Pemathang, Barun Bazar, Gola, Arun III hydropower Project will be at great risk and could be fatal if GLOF at Lower-Barun Lake occurs.
The 100% release of the lake water is quite rare but can happen only in specific topographical conditions. For instance, 95% of the water volume was released during the GLOF event in Peruvian Lake Palcacocha in 1941 (Emmer, 2017). GLOF nature varies with the various reasons like triggering factors and lakes characteristics (Westoby et al. 2014). Lower-Barun Lake has an expansion rate of around 85% in the last four decades; this implies more GLOF hazards as it can expand more in the future due to the current glacier retreat trend. Nepal lies in a seismically active zone. So, Nepal Himalaya has experienced several mass movements events like avalanches, rockfall, landslides (Kargel et al. 2016).

Conclusions
The lake area increased from 0.305 km 2 in 1979 to 2.019 km 2 in 2018. The average percentage change in the lake area per decade is found to be 37.84%, and overall, 85.58% change is observed in 40 years' time period. An average increase of 0.043 km 2 area per year is observed. The glacier area decreased from 24.61 km 2 to 22.63 km 2 over the past four decades. The shrinkage in the glacier area is around -0.49 km 2 in average and has shrunk 8.02% in the last four decades. The glacier retreat in Lower-Barun is found to be 0.20% in a year in the last four decades. The GLOF susceptibility index found to be 0.94, which suggests that the lake is very high susceptible to GLOF. Furthermore, the water volume of the lake is found to be ～103.6 × 10 6 m 3 . The peak discharge of 5768 m 3 /s is produced when the breach depth is 20 m and the entire water volume is released. Likewise, in the case of 15 m breach depth, the peak discharge of 4038 m 3 /s is formed.