Mapping and monitoring of glacier lake outburst floods using geospatial modelling approach for Darkut valley, Pakistan

Climate change and human activities have resulted in the receding of glaciers throughout the world including Pakistan. Glacier lake outburst floods (GLOFs) are amongst the most common climate‐change‐induced hazards in northern Pakistan. In the present study, GLOF mapping and modelling was carried out using remote sensing and geographical information system techniques coupled with ground‐truthing. Change detection techniques such as the normalized difference water index were applied on Landsat imagery for the identification of the temporal behaviour of Darkut glacial lake for the last 25 years. The depth of the lake was estimated to be 81 m and the volume of the lake was calculated using a digital terrain model and extracted as 9.79 × 106 m3. The glacial lake extent has increased from 0.045 to 0.154 km2 in the last two decades. Two GLOF scenarios (peak and extreme flood) were developed on an existing volume of water in the study. There are 14 households exposed to medium flood and 10 to low flood risk while one helipad and one school are also in the low flood zone in the first scenario (i.e. peak flood) based on 87.84 m3·s−1 of water. The second scenario (i.e. extreme flood) was executed on 3,128 m3·s−1 of water, in which 14 households are at high flood risk, eight at medium and 35 in a low flood zone, as well as one school, a helipad and a community stockpile which are exposed to low flood. The outcomes of the study will help in the development of risk management plans, preparedness strategies and risk reduction from GLOF hazard.


| INTRODUCTION
The change in temperature has affected global glaciers in general but it has been noticed that, among all, the glaciers of the Hindu Kush Himalayan (HKH) regions are majorly impacted. They are receding with a high pace and rate of snowfall due to which less snow is accumulated in the glaciers. These gradual variations in the glaciers lead to the development of new water bodies and the extension of old lakes that can pose a great threat to downstream vulnerable communities in the case of glacier lake outburst floods (GLOFs) (Xu et al., 2009;Ashraf et al., 2017).
The glacial lakes are classified on the basis of their position relative to the glacier and the damming mechanism, e.g. pro glacier, en glacier, sub glacier, supra glacier and moraine-dammed lakes (Salerno et al., 2012;Zhang et al., 2015). In 2005, the International Centre for Integrated Mountain Development (ICIMOD) published an inventory of glaciers and glacial lakes which reported that the trend of rising temperature in the HKH regions is higher than the global average. The inventory listed a total of 2,420 glacial lakes that have been identified in 10 river basins; the maximum number of glacial lakes is identified in the Gilgit River basin (614) followed by the Indus basin (574) (Campbell and Pradesh, 2005). In the Gilgit River basin, out of the 614 glacial lakes 380 lakes have been identified as major lakes, which are about 62% of the total lakes (Ashraf et al., 2012). Around 93% of the lake area of the basin is contributed by these major lakes. There were 52 potentially dangerous lakes in the HKH mountain regions of Pakistan; among these potentially dangerous lakes 16 lakes were identified in the Karakoram ranges (Zaidi et al., 2013). Himalaya is known as one of the most vulnerable regions throughout the world due to rapid climate change (Liu et al., 2014a;2014b;Nie et al., 2017). In 2007, the Intergovernmental Panel on Climate Change (IPCC) reported that there will be an increase in temperature in the Himalayan region of between 1 C and 6 C by 2100 (Parry et al., 2007). This will result in the high melting of glacier ice and the snow cover will reduce by 43% to 81% leaving behind more glacial lakes in the future (Bajracharya et al., 2015).
The HKH ranges are the biggest source of glacier ice and snow, which contribute to feeding the largest river systems. Consistent mapping and monitoring of these dynamic resources is pivotal to keeping track of the pattern and trend of changes. The temporal mapping of glacial lakes provides an indicator for monitoring changing climate patterns (Jilani et al., 2008). Gilgit-Baltistan has witnessed more than 35 GLOF events in the past 200 years (Watkins, 2007). In 1994 a huge GLOF occurred in Sosot village, Gupis valley, and a lake was formed due to the accumulation of river water which blocked the Ghizar river. This event took five lives and caused heavy destruction to cropland and other property. In the same valley, another GLOF event occurred in 1999 at Gupis Khalti village, which also blocked the Ghizar river and created a lake 1.5 km long known as the "Khalti Lake" which was breached after a few months and caused huge destruction in downstream settlements. In 2007In , 2008 two GLOF events from Ghulkin glacier caused destruction of property and blocked Karakorum Highway (PMD, 2016).
The analysis of the GLOF disasters shows that these events are linked with extreme weather conditions such as a rise in temperature and intensive rainfall (Khalsa et al., 2004;Bajracharya et al., 2007). Remote sensing (RS) and geographical information systems (GISs) are marvellous resources available due to the development in the field of space technologies; they provide an insight into the different aspects and characteristics of highaltitude remote glaciers and glacial lakes that are physically very difficult to access to assess the situation (Banerjee, 2003;Zhang et al., 2019).
RS and GIS tools and techniques make it possible to study, analyse and monitor consistently the changes and developments taking place at such sites. Many researchers have used GIS and RS techniques for the assessment of GLOF events, e.g. Somos-Valenzuela et al. (2015) and Zaidi et al. (2013). Different characteristics of glaciers and lakes such as identification of the snow line, glacier moraines, lake boundaries and area can easily be demarcated using satellite images which can be validated by ground survey. However, for the development of the GLOF model, several other parameters are taken into account that need ground assessments. Therefore, the best way is to apply both approaches to achieve good results (Banerjee, 2003;Mal et al., 2016). For unsteady GLOF simulation and modelling, a digital terrain model and field survey data are crucial. The Hydrologic Engineering Center's River Analysis System (HEC-RAS and HEC-GeoRAS) is one of the emerging advanced tools for flood modelling and identification of downstream flood zones (Matkan et al., 2009).
The study aimed to assess and analyse the dynamics of glacial lakes and to develop a model of possible GLOFs using RS and GIS technology coupled with HEC-RAS and field assessments and surveys. The objectives of the study were (a) to analyse temporal changes in the glacial lake using time series  RS data (Landsat) and (b) to develop GLOF scenarios and identify downstream inundation zones.

| STUDY AREA
Darkut glacial lake (Figure 1) was selected for this research due to the risk it poses downstream. It is located in Yasin valley, District Ghizar, in northern Pakistan. It lies between the latitudes 36 37˝0 0 -36 44 00 N and longitudes 73 24 00 -73 27 00 E. Gilgit city (the capital city of Gilgit-Baltistan) is located at a ground distance of about 140 km southeast of Darkut village. The village is inhabited by 260 households.
Darkut Glacier which is also known as Ghamu Bar Glacier is located at an altitude of about 2,751.12 m above sea level. The glacier is about 3,900 m long and the average width is 644 m. The total area covered by glacier ice is approximately 2.51 km 2 . The snout of the glacier is partially covered in debris and saturated crevasses; however, the debris is saturated towards lateral moraines. Seasonal snow on the mountains contributes a lot to water runoff. There is no supra glacier drainage on the surface but glacial meltwater flows through inglacier paths. The height of the glacier ice at the snout is 39 m and the distance of glacier ice from new terminal moraine is 151 m (Focus, 2013).
Darkut Lake is formed at the terminus of the glacier ice; ice melted water is accumulated by the old terminal lateral moraines. Due to the western lying position of the glacier, it goes through high melting and continuously supplies water to the lake. The study conducted by the Pakistan Meteorological Department (PMD) in 2016 reveals that the lake is spread over 0.48 km in length and 0.30 km wide with an approximate depth of 81 m.

| MATERIALS AND METHODS
3.1 | Data acquisition 3.1.1 | Field data collection (river cross-section and discharge) A detailed field survey was conducted in August 2015 to collect all the data required for mapping and modelling of the GLOFs. The field survey consisted of river profiling, i.e. river cross-section data collection using GPS, survey of flood plains and mitigation structures along the river or flood path, identification of flood marks by visual F I G U R E 1 Map of the study area (Darkut glacial valley) observations, discussion and interviews from local people regarding the history of GLOFs, photographs of different sites along the flood path and mapping of flood inundation zones on base maps. The "float method" was used for flow discharge calculation (Trimmer, 1994;Tsubaki et al., 2011). The standard equation used for discharge calculation is as follows: where Q is the discharge, V is the velocity and A is the width × depth.
To ensure accuracy, this step was repeated three times and, for each flow, travel time was recorded and an average was used in further calculations.

| Satellite images
Landsat images of 30 × 30 m spatial resolution were used to assess the lake dynamics and monitor the spatial and temporal changes. Three snow and cloud-free scenes of Landsat-5 Thematic Mapper for 1991, 1996 and 2010, two scenes of Landsat-7 Enhanced Thematic Mapper Plus for 2000 and 2005 and one scene of Landsat-8 Operational Land Imager for 2015 from the month of August were acquired from the USGS website (https://earthexplorer. usgs.gov). Green (0.52-0.60), red (0.63-0.69), near infrared (0.76-0.90) and short wave infrared (1.55-1.75) bands were used for analysis. The month of August was chosen due to the fact that it is peak summer time which leads to heavy melting of glaciers that maximizes the lake volume. Also, in July and August the monsoon season provides heavy rains that add more water to the glacial lake that increases its size and volume. The higher the volume of the lake is, the higher is the chance of outburst due to the increased pressure on the sediment material of the outlet and the chance of moraine break, leading to a GLOF.

| Digital terrain model (DTM)
A digital terrain model (DTM) of 10 m resolution was acquired from Focus Humanitarian Assistance Pakistan, clipped to the area of interest and used in flood modelling to calculate the glacial lake volume to develop GLOF simulation to anticipate downstream risk. The 10 m DTM was created by the data providers through a combination of Shuttle Radar Topography Mission (SRTM) C-band (SRTM radar contains two bands, i.e. C-band and X-band; digital elevation models (DEM) are made from C-band radar data), Global 30 Arc-Second Elevation (GTOPO-30), Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), and Ice, Cloud and Land Elevation Satellite (ICESat) datasets.
The data were acquired between November 2, 2000 and January 5, 2008. These data were projected to the latitude and longitude coordinate system, WGS84 horizontal datum, EGM96 vertical datum, with horizontal measurement units in decimal degrees and vertical measurement unit in meters.

| Lake assessment
The analysis of temporal variation in the lake area was carried out using the satellite images. Due to the unavailability of free high-resolution imagery, Landsat images were used for classification and change detection using two indices, i.e. the normalized difference water index (NDWI) and the modified normalized difference water index (MNDWI). The NDWI (McFeeters, 1996) and MNDWI (Han-Qiu, 2005) were calculated as: where NIR is near infrared and SWIR is short wave infrared.
To ensure satisfactory results for NDWI, the MNDWI was proposed by Han-Qiu in 2005 in which the NIR band was replaced with the SWIR band. MNDWI can effectively depress information from built-up and other background features while expressing water feature information and making it possible to extract water body information accurately from the target area due to the higher reflectance in the green and NIR bands. The threshold value for MNDWI was set to zero like the NDWI but manual editing brought more refinement and accuracy in the delineation of water bodies (Aggarwal et al., 2013).
Both NDWI and MNDWI use the highest reflectance from the green with the NIR and SWIR bands for the extraction of water bodies (Kaushal and Kamal, 2016). The outcome ratio images of both NDWI and MNDWI were further classified into a series of specified intervals for delineating the glacial lake. The resulting images were visually inspected and correlated with Google Imagery and available field photographs.
In this study, a land-water threshold was manually applied to classify images into two classes, i.e. land and water. Suitable land-water thresholds for each index were determined through trial and error and observation of the histogram. The NDWI values for the lake surface, in this study, range from 46 to 75 for glacial lake mapping. Extensive manual editing and delineation of the lake were carried out using Google Earth images and were adjusted with reference data, i.e. photographs and ground observations, for mapping the final lake.
These indices indicate the presence of water. The NDWI values range from −1 to +1. Most of the water features are found close to the +1 value. McFeeters (1996) set zero as the threshold value for water bodies. The values towards −1 indicate vegetation and bare soil or land features.

| Calculation of lake volume
The volume of the lake is mostly calculated by measuring the area of the depth contours using area measurement devices on hard format maps. Bathymetric data are prerequisites for the generation of contours that can be used to measure the depth. The availability of new technology helped for calculating lake volume with high accuracy; for example, the digital elevation and DTM are widely used for such calculations and modelling purposes. In the absence of the technology and accurate bathymetric data, the volume of the lake can be estimated using depth as a function of distance based on the assumption of a linear increase in depth with distance from the lake shoreline, and lake depth can be measured at any given location by using the following simple equation: where D is the distance from the shoreline, Z max is the measured maximum depth at any given location and D max is the maximum distance from the shoreline. The F I G U R E 2 Methodology flow chart volume calculation using this method is basically the calculation of the area and depth (Hollister et al., 2010). The PMD conducted an assessment of Darkut Lake in 2016 using the boat and weighted cord method; therefore, the depth data were acquired from the PMD for further analysis. The GPS points containing the latitude, longitude and depth values were interpolated in the ArcGIS software environment using inverse distance weighted interpolation. The polygon of the shoreline boundary of the lake was created using survey points and the attribute was populated with elevation information (i.e. Z values). The lake boundary or shoreline depth was taken as zero. The output surface resulting from the interpolation was converted into integers and created the surface of a triangular irregular network. The volume of the lake and area were calculated using the surface volume tool in the 3D Analyst toolbar of ArcGIS.

| GLOF modelling using HEC-GeoRAS and HEC-RAS
The collected data were used to draw a sketch of the river profile/cross-section in HEC-GeoRAS (plug-in) in the ArcGIS environment to understand the river system. Geo-RAS required the DTM as a base layer that needs to be in a projected coordinate system. The flow path centreline was extracted from the DTM, and bank lines and cross-sections were delineated based on field GPS data. Cross-sections were digitized perpendicular to the flow path. Another necessary parameter to add is the Manning roughness value (n). The HEC-GeoRAS tool uses land-use features and n values against each class. Therefore, for mountainous streams and rivers a Manning value of 0.05 was used (Chow, 1959).
The geometric data and unsteady flow conditions were entered in HEC-RAS and simulation was carried out. As HEC-RAS functions on the basis of the Saint-Venant equation, it is one of the appropriate tools for generating flood scenarios caused by lake outburst. In HEC-RAS, after the analysis of unsteady flow, these data were used for the development of flood maps. The data were imported into the ArcGIS environment to produce maps. The flood polygons were overlaid with infrastructure data (collected in the field) and identified infrastructure at risk of possible flood zones.

| Land use and land cover (LULC) maps
The maps for LULC of Darkut village were produced using supervised classification techniques on Landsat satellite imagery for 1991 and 2015. The classified raster output was transformed into vector format, i.e. shapefile, to carryout analysis of LULC and to develop maps. The area of each class was calculated using the Field Calculator tool in ArcGIS. Graphs for the area of different land use classes were produced in MS Excel and LULC maps of Darkut were prepared to show the inundation area in the case of GLOF.

| Downstream flood inundation zone maps
The cross-sections and the surface elevation data were imported into ArcGIS HEC-GeoRAS and a water depth surface map, an intensity surface map and a time arrival zones map were created. A risk map was also developed in ArcGIS by overlapping flood zones and infrastructure GPS data. A detailed flow chart of the methodology is shown in Figure 2.

| RESULTS AND DISCUSSION
The study was conducted to assess and analyse the impact of possible GLOFs on downstream settlements, using field survey data coupled with RS and GIS techniques. The field survey revealed that there was no supra  glacial lake on the glacier and no supra glacial drainage was observed. It is estimated from the watermarks that the water level in the lake is decreasing but in the summer season, due to the melting of the glacier, the inflow of water increases, which can lead to flood situations downstream in the case of outbursts. The slope on the north of the lake is highly inclined compared to the opposite side. The moraines barrier contains dominantly gravels, silt and trace boulders; the height of the moraine barrier is 15 m from the current water level in the lake. The current year (2015) spillway of the lake is 5 m wide and the stream bed is mature which prevents erosion of spillway during normal flow conditions. It can be eroded during intensive melting of the glacier, feeding a large amount of water into the lake which may cause flood conditions downstream. In short, the lake poses a potential threat of a GLOF to Darkut village.
Temporal NDWI maps were developed. A threshold was selected in the NDWI values of each year and values between 0.46 and 0.75 were found suitable to map the spatio-temporal variation in the glacial lake. The results from the NDWI for 1991 revealed that the lake area was 0.045 km 2 , which was taken as a base value to compare other values to infer the difference in lake growth in each year. Compared to 1991, 2.7% of the lake area was reduced in 2000. According to the historical data collected during participatory risk assessment, the glacier ablation zone has experienced a drastic melting which formed crevasses in the glacier debris. With the passage of time, small water bodies were developed inside the glacier debris which gradually melted and increased the lake area and volume. In 2015, the lake area increased almost triple and measured 0.154 km 2 . Figure 3 shows the temporal change in the lake area.
Time series LULC mapping shows that the glacial lake area and other features witnessed a drastic change in the last 24 years. The lake area was 0.04 km 2 in 1991, which was increased to 0.15 km 2 in 2015. Agricultural land and rangeland were also estimated to be increased in this time period. Glacier area was decreased to 26.32 from 33.15 km 2 ; similarly, barren land was reduced to 20.38 from 22.75 km 2 as agricultural land and rangeland were increased due to increased population and availability of water for irrigation from the melting glacier, as shown in Table 1. The temporal LULC map of the study area is shown in Figure 4.
A 2D unsteady flow model was used for simulation of the glacial lake. In this regard the depth of the lake was calculated using GPS surveying data acquired from the PMD; the interpolation of GPS points was performed in the ArcGIS environment. The depth was estimated to be 81 m. As the depth of the lake increases, the area decreases simultaneously. Most of the lake area is within a depth range of 1-16 m whereas a small area has a depth within the range 62-81 m. Due to the silting of soft material in summer, the depth of the lake is reducing which may trigger a flood due to high pressure on the lake outlet.
The estimation of Darkut glacial lake breach location, dimensions and time scenarios was developed using field data acquired through surveying of upstream and downstream areas of the glacial lake. Two scenarios, peak and extreme, based on GLOF potential risk to loss of life and property damage were developed. These two scenarios were developed using 2D unsteady flow routing (the full Saint-Venant equation, also called the diffusion wave equation) and were based on the amount of water and the location and length (size) of the breach of the outlet area. The results of both scenarios were in the form of four types of maps: a depth map, a velocity map, a flood time arrival map and a hazard/ risk map. Peak flood scenarios (i.e. the first scenario) of flow hydrography were constructed on 87.84 m 3 Ás −1 of water received every 2 min. The simulation timing was set at August 10 01:00 p.m. to 04:18 p.m. The computational interval was 15 s and 2 min was set for mapping the hydrograph interval, the hydrograph output interval and the detail output interval. All four maps of the first scenario are shown in Figure 5. Red in the flood time arrival map ( Figure 5) indicates a shorter time span, while yellow and orange show a longer time span taken by the peak flood to reach that area. The peak flood velocity map shows the risk of high intensity flood near the outlet, whereas downstream areas seem to have less intensity flood; this is due to the agricultural activities and mitigation measures. The risk map ( Figure 5) depicts the infrastructure at risk of flood with respect to flood water depths and the velocity of water in the case of a peak flood event.
During peak flood, the depth of the water may vary at different locations based on the land use. Table 2 shows different land use categories under various flood depths. The results show that most of the land types are affected by a low depth flood due to the low gradient downstream. A total area of 1.93 km 2 was estimated to be flooded to different depths. Built-up and agriculture land of about 0.48 km 2 and 0.63 km 2 respectively was calculated to be flooded at the depths between 2 -7 m. Similarly, 0.51 km 2 barren area was estimated to be flooded by between 2 and 15 m floods. Barren land was estimated to be flooded under different depths of flood water between 2 and 20 m.
According to the elements at risk table (Table 4), a total of 24 households, one school, and one helipad facility are at risk of flood if a GLOF event occurs during peak season.
An extreme scenario (second scenario) of flow hydrography was constructed on 3,128 m 3 Ás −1 water received every 2 min. The simulation timing was set as August 10 01:00 p.m. to 04:18 p.m. The computational interval was 15 s and a 2 min interval was set for mapping the hydrograph, the hydrograph output interval and the detail output interval. The second scenario results are also shown in four types of maps ( Figure 6). The second scenario of a GLOF event would affect a wider area and more facilities due to the wider spread and the high velocity and elevation. An extreme flood would impact barren land the most, as it may affect 8.94 km 2 of land, out of which 0.94 km 2 would be under 14 m of flood water. The velocity map depicts the situation based on intensity differentiation, the higher the velocity, the higher the impact downstream ( Figure 6). Table 3 shows that a total of 21.71 km 2 of land downstream would be flooded if extreme flood events occurred due to lake outburst. This includes a 2 km 2 built-up area which would be submerged with between 0 and 29 m of flood water and an area of 0.004 km 2 which would be affected with between 30 and 43 m. Agricultural land of 2.70 km 2 would be inundated under 0-43 m, out of which 1.21 km 2 would be under 14 m, 1.48 km 2 would be under between 15 and 29 m and 0.008 km 2 would be submerged with between 0 and 43 m. This flood would impact barren land the most, as it would flood 8.94 km 2 land, out of which 0.94 km 2 would be under 14 m, 7.97 km 2 would be under between 15 and 29 m and 0.024 km 2 would be under between 30 and 333 m of flood. Similarly, 1.14 km 2 of rangeland and 6.92 km 2 of freshwater sources would be affected with between 0 and 29 m of flood.
The infrastructure at risk due to an extreme GLOF event is shown in Table 4. According to the statistics 14 houses would be exposed to high flood, eight to medium and 35 houses to a low intensity flood. Besides that, one helipad and a school and community stockpile would be exposed to low flood.

| CONCLUSIONS
A glacial lake outburst flood (GLOF) is among the most devastating hazards in Gilgit-Baltistan. Due to changing climate and human-induced hazards, glaciers in northern Pakistan are receding and ultimately forming lakes, which pose a great threat in the form of a flash flood in the case of a GLOF for downstream settlements. In the present study, the findings are derived based on a detailed field survey, remote sensing data and secondary data acquired from the Pakistan Meteorological Department and Focus Humanitarian Assistance Pakistan. Downstream inundation areas in the case of a GLOF were identified and mapped through modelling techniques using HEC-RAS and HEC-GeoRAS software coupled with ArcGIS.
Based on the analysis of the data and simulation results, it is concluded that Darkut glacier is retreating and the lake is expanding at an unprecedented speed, which poses a great danger for Darkut settlement. Visual observation during a field survey also shows that the slopes at the north of the lake are unstable and any triggering event, e.g. an avalanche or earthquake, may lead to a GLOF. Therefore, it is recommended to conduct a detailed study to examine the capacity of lake outlet and a strategy needs to be developed to relocate households at risk.