GLOF in the South Lhonak Lake, India: photogrammetric analysis and estimation

Glaciers in the Himalayan region are considered the freshwater towers of South and East Asia and are strongly affected by climate change (Vohra 2007). In recent decades, most of the world's mountain glaciers have experienced adverse mass balance and terminal recessions (Haeberli et al. 1999). Swift declines have been observed across the Greater Himalayas, causing the formation of new lakes (Ageta et al. 2001). Lakes at the glacier's terminus are primarily confined by lateral or end moraines, posing a heightened risk of breaching. Due to their potential for holding significant water volumes, these lakes can be hazardous. Breaching and the sudden release of water from these lakes can generate flash floods, leading to substantial damage in downstream areas and loss of lives and infrastructure (Allen et al. 2016; Aggarwal et al. 2017). The 2013 Kedarnath disaster (central Himalaya), which was a glacier lake outburst flood (GLOF) from the Chorabari Lake associated with flash-flooding and landslides triggered by intense precipitation in the region, led to over 6,000 fatalities and washed away a major part of a settlement located just over 1 km downstream of the lake (Allen et al. 2019). The risk that a potentially dangerous lake may present in low-lying areas can only be recognized if a detailed hazard assessment is undertaken (Nie et al. 2018). The probability of a GLOF event is, however, difficult to estimate (Wang et al. 2015; Maskey et al. 2020) due to rapid changes in the glacial systems, the low frequency of such events, and the high complexity of the involved processes. Therefore, evaluating the characteristics of a glacial lake, its surroundings, and local hydrodynamical flood modeling allow us to understand how a valley will behave in case of a potential GLOF event (Khanal et al. 2015; Wang et al. 2015). Glacier lakes are inaccessible for evaluating the cause and quantification of geomorphologic changes, emphasizing the importance of advanced remote sensing techniques for hazard quantification (Budhathoki et al. 2010). The surface area of a lake can be determined using remotely sensed data. However, when estimating the potential peak discharge due to a possible outburst, it is essential to know the volume of the glacier lake rather than just its area (Benny & Dawson 1983).

The reliability of these remote sensing-based bathymetric studies depends on finding a significant correlation between water depth and reflected energy (Baban 1993). Collecting glacier lake bathymetry data in the Himalayas is challenging, time-consuming, and costly, with no studies using remote sensing for this purpose. Cao et al. (2019) and Anees et al. (2022) reviewed the various photogrammetric methods for obtaining bathymetry of lakes that are inaccessible and shallow in depth. The widespread availability of affordable remote sensing platforms with high spatial and temporal resolution facilitates swift, semi-automated, and cost-effective evaluations of glacier parameter changes over large areas (Racoviteanu et al. 2008). Halder & Bose (2024) have demonstrated the integration of Google Earth Engine (GEE) with Sentinel-1 SAR data and Landsat 9 imagery for a comprehensive, near-real-time flood mapping and analysis of GLOFs.

This remote sensing methodology involves extracting topographic information from satellite images (Kumar & Bhardwaj 2021) and volume from the satellite stereo pairs, which enables the generation of precise digital elevation models (DEMs) (Janet et al. 2006; Tsutsui 2007; Giribabu et al. 2013). However, pre- and post-stereo pairs are to be available to estimate the volume of water or debris displaced in the event, which is only sometimes possible. Willneff et al. (2005) demonstrated the possibility of deriving the depth from single-image photogrammetry using precise refined rational polynomial coefficients (RPCs), image coordinates, and iterative image coordinates tracing methods. Claudio Bozzini (2014) and Bhushan et al. (2021) studied the accuracy of the DEMs extracted using stereo-matching techniques from the stereo images. Also, Atefi & Miura (2021) and Robson et al. (2022) have demonstrated the computation of volumes of landslides using satellite stereo data of before and after events. Hence, geomorphology changes and land patterns can be assessed by comparing these DEMs. In Sikkim, glaciers terminating in lakes have increased growth compared to those without lakes (Marti et al. 2016; Aggarwal et al. 2017). Therefore, evaluating the characteristics of the South Lhonak Lake, nestled in the Himalayan region, has faced environmental challenges, with glacial and landslide dynamics being key contributors to its vulnerability. Sharma et al. (2018) and Sattara et al. (2021) presented the susceptibility of damage due to GLOF and future growth of the South Lhonak Lake using satellite data, DEM, and other ancillary data.

This study utilizes advanced photogrammetric techniques and high-resolution satellite imagery to analyze the landslide debris and depletion of water storage in the South Lhonak Lake, Sikkim, India. A combination of both stereo and mono imaging is done, covering the image collection dates during pre- and post-GLOF events. In addition, the SAR images from the RISAT-1A satellite are also used to interpret the glacier lake area after the GLOF. The specific importance of this study is that it offers a comprehensive analysis of changes in the glacier lake, which is crucial for understanding the catastrophic event that occurred on 3 October 2023. The methodology involves data satellite acquisition, photogrammetric image processing, generation of pre- and post-DEMs, change detection, and estimation of the volume of the landslide debris and amount of water depletion from the lake due to GLOF.

The South Lhonak Lake is located at about 5,200 m above sea level in the southeast of the Indian state of Sikkim. The study area encompasses the South Lhonak glacial lake, situated at coordinates 27.913° N and 88.199° E, along with the downstream drainage area of the Chungthang Dam site (which is the confluence of the Lachen and Lachung rivers). The lake mainly receives melt water from the South Lhonak Glacier and is supplemented by the outflow from the North of Lhonak glacial lake. The glacier has a total area of 12.5 km², as mapped in 2019. In line with the glacier retreat, the lake has been exhibiting significant growth over the years as it grew from 0.42 km² in 1990 to 1.35 km² in 2019 (Sattara et al. 2021). According to measurements of the lake's depth, the lake's volume is 65.8 × 10⁶ m³, with the deepest part measuring 131 m (Sharma et al. 2018).

High-resolution stereo satellite images were used to analyze the terrain dynamics surrounding the South Lhonak Lake before and after a GLOF event. Data used for analysis are tabulated in Table 1.

Data selection and remarks of the satellite data used are listed in the Table 2.

Resolution variability: Different resolutions (60 cm–2.5 m) can lead to inconsistencies in DEM accuracy and terrain analysis, with lower-resolution data potentially missing finer details.

Sensor differences: Variations in sensor capabilities across different satellites can affect data quality, base-to-height ratios for stereo fusion and compatibility, requiring consistent processing to minimize discrepancies.

The selection of high-resolution satellite data from CartoSat-3 and CartoSat-2S ensures detailed and accurate analysis of terrain changes around the South Lhonak Lake. However, resolution variability, temporal gaps, and sensor differences pose challenges that need careful management to ensure reliable results. Consistent processing and careful data integration are essential for the reliability of the findings.

A DEM is a representation of the topography or terrain of a surface in digital format. It is typically a raster grid where each cell contains an elevation value representing the height of the surface at that location. DEMs are generated using multiple methods such as Satellite Stereo Imagery, Satellite Radar Interferometry (InSAR), airborne stereo photogrammetry, or Lidar (Light Detection and Ranging) methods. However, the present study uses satellite stereo imagery as there was an opportunity to capture pre- and post-event scenario of the South Lhonak Lake from the archived satellite data and current sceanrio through specific planning for satellite data acqusition (Evans et al. 2008).

Photogrammetry is a sophisticated technique used to derive three-dimensional (3D) information from two-dimensional (2D) images. This method involves capturing multiple images of an object or terrain from different angles and processing them to generate an accurate 3D model. Detailed breakdown of the process is explained in the following steps.

This initial step involves acquiring images of the target area or object. For the study area, satellite images were acquired from different satellites (CartoSat-1 on 9 October 2015; CartoSat-3 on 13 May 2023; and CartoSat-2S on 11 and 12 November 2023).

It is the process of re-establishing the precise relationship between all images and real-world coordinates as they were at the time of image acquisition. This is done by refining rational polynomial equations coefficiets, which map the relationship between the image plane, the ground plane, and the camera by principle of triangulation and it is achieved by following steps.

Image orientation: This involves aligning the images with each other by measuring tie points across overlapping images and RPCs. In this process, over 350 semi-automatic tie points were used to refine the RPCs across all the images. The tie points and root mean square error (RMSE) accuracy, which measure the alignment accuracy, resulted in a value of 0.537 pixels. This accuracy complies with the requirement of being less than one pixel, which is the target after triangulation and adjustment.

Resolution variability: The study utilizes satellite data with varying resolutions ranging from 60 cm to 2.5 m. This variability in resolution leads to inconsistencies in the DEMs generated from these images, as higher-resolution images can capture finer details compared to lower-resolution ones. These differences can affect the accuracy of terrain change detection, as smaller features might be visible in high-resolution images but missed in lower-resolution ones, potentially leading to inaccurate analyses of changes over time.

Temporal gaps in data acquisition: Significant temporal gaps between data acquisition dates, spanning from 2015 to 2023, present a challenge in capturing continuous changes in the terrain. These gaps can result in missed intermediate changes, which are crucial for accurate analysis. Additionally, the loss of stereo-image pairs necessary for generating accurate DEMs is a concern. In such cases, the study resorts to single-image photogrammetry, an iterative process for deriving DEMs from single images. However, this method is less accurate than using stereo pairs.

Geometric and radiometric corrections: Geometric corrections ensure that satellite images accurately represent the spatial relationships of terrain features, while radiometric corrections normalize the brightness and contrast of images. Differences in these corrections applied to images from various satellites can introduce errors and inconsistencies in the DEMs. Ensuring uniform geometric and radiometric corrections across all images is crucial to minimize discrepancies and maintain accuracy in the data processing.

Multiple sensor data: The lack of availability of a single stereo pair for the period just before the GLOF event necessitated forming stereo pairs from multiple dates and sensors. Data from different sensors acquired at various altitudes and angles lead to variations in base-to-height (B/H) ratios, which can cause inconsistencies in the accuracy of the DEMs. To address these issues, the study considers the lowest accuracy for both pre- and post-GLOF DEMs, acknowledging that the inherent differences in sensor data quality impact the reliability of terrain change detection.

Once the photogrammetric processing is complete, stereo pairs are prepared from oriented overlapping images. Corresponding matching points are then measured semi-automatically. A point cloud (elevation points) is generated by measuring these corresponding matching points in the overlapping stereo pairs, representing the 3D coordinates of points on the surface as shown in Figure 3(d). This step is essential for creating a detailed 3D model of the terrain. Using the point cloud data, the surface of the terrain is reconstructed in the form of a DEM at a 5-m spatial interval for both pre- and post-GLOF events, as explained in the following sections.

Part A DEM: A DEM of this region was created using a stereo pair of CartoSat-3 images from 13 May, along with CartoSat-1 data at a resolution of 2.5 m in the area known as Part A. This combination of data suggests that the terrain remained relatively unchanged between the dates of acquisition, providing consistent elevation data.

While the stereo pair yielded highly accurate elevation information for 13 May in Part A, there may be slight variations in terrain in other sections (B and C) where stereo imagery was not feasible, making this pair unsuitable for generating DEMs in Parts B and C.

To produce a precise DEM for Part A, mass points and break lines were manually extracted, and accuracy was meticulously verified using stereoscopic methods. Regarding DEM accuracy, the stereo pair had a convergence angle of approximately 27°, with a B/H of 0.52 and a stereo fusion accuracy of 2 m (as both images were re-sampled to 2 m). Consequently, the vertical accuracy of the DEM was determined to be 3.8 m.

Part B DEM: The stereo pair of images to create a DEM in this area is acquired on two dates. The first image was obtained on 13 May 2023, using CartoSat-3, while the second image was taken on 11 October 2023, with CartoSat-2S. Despite the excellent stereo configuration, we encountered difficulties working with the lake outlet and landslide zone. DEM in this area is generated by manually adding mass points and break lines with vertical accuracy of 2.0 m, having a B/H ratio of 0.3, and a stereo fusion of 60 cm.

Stereo images were acquired from CartoSat-2S on 11 and 12 November 2023, from CartoSat-2S, with a converging angle of 43.54°. These images were utilized in LPS Leica photogrammetric Software (LPS) and manual mass points and break lines were drawn on stereo pair to create elevation terrain data (Krishna et al. 2008; Kugler & Wendt 2021; Kumar & Bhardwaj 2021; Siva Subramanian et al. 2023; Barbarella et al. 2017). The resulting DEM, with a 5-m spacing, effectively depicts the height information of the South Lhonak Lake area and its surroundings following the lake breach, as illustrated in the accompanying image (Figure 4(j)).

To assess the accuracy of the pre- and post-event DEMs for the GLOF scenario, a detailed evaluation was conducted using manual stereoscopic analysis on satellite stereo models (Table 3).

This process involved superimposing the DEMs onto their respective stereo pairs and examining errors at specific points. For the pre-event DEM, errors were evaluated on 120 points, yielding a maximum residual of 2.6 m and a minimum residual of −1.5 m. The RMSE for the pre-event DEM was found to be 1.27 m, with a linear error at 90% confidence level (LE90) of 2.52 m. The pre-event DEM was created using a combination of two stereo pairs and single-image photogrammetry, resulting in an absolute vertical accuracy of about 2.5 m and a relative vertical accuracy of 2 m. Similarly, the post-event DEM, assessed on 100 points, showed a maximum residual of 1.8 m and a minimum residual of −0.75 m. The RMSE for the post-event DEM was calculated as 1.09 m, with an LE90 of 2.1 m. This DEM was generated using a single stereo pair of satellite sensor data, with a B/H ratio of 0.32 and a stereo fusing accuracy of 60 cm, leading to a vertical accuracy of approximately 2 m. Overall, the accuracy of volumetric calculations derived from these Digital Elevation Models (DEMs) was estimated to have vertical uncertainty of approximately ±4 meters. This estimate applies to the surface differences captured in the Difference of DEMs (DoDs) analysis. The uncertainty accounts for potential errors in vertical elevation measurements as well as limitations associated with stereoscopic analysis techniques.

Once the pre-event and post-event DEMs are created, they are compared by differencing to identify changes in topography. This process involves subtracting the pre-event DEM from the post-event DEM to detect alterations in the landscape. In this specific case, the analysis reveals a portion of land that has slid into the lake, triggering a wave surge. This surge subsequently caused a breach at the lake mouth, leading to significant depletion of water in lake. The topographical changes and their impacts, including the land slide, lake breach, and water depletion, are illustrated in Figure 5(a).

Observations and results of the analysis derived in this study are presented in Tables 4 and 5.

In this study, the debris slide volume was estimated by subtracting the pre-event DEM from the post-event DEM and integrating the elevation differences in the affected areas. This resulted in a calculated total volume of 11.5014 million m³. The total projected area of debris slides up to the water level after the event is 0.2238 km², and the maximum height of debris above the water level is 202 m. These factors contributed to the mass outflow of water from the lake, causing the GLOF. The total water level drop after the event is 32 m, and the total water area and perimeter of the lake after depletion are 1.491 km² and 6.679 km, respectively. The maximum length of the lake is 2,600 m (horizontal), the maximum width is 760 m (cross-section), and the minimum width is 18 m (cross-section). After the incident, erosion occurred at the mouth of the lake, measuring 30 m. The cross-sectional area through which water emptied at the mouth is 3,840 m². Through this mouth cut, the total volume of water emptied from the lake is 43.2172 million cubic meters (MCM), with an average discharge of 4,004.63 m³/s. The total duration for the complete lake depletion is 3 h which is taken from news papers.

The methodology used in this paper enhances understanding and ability to assess GLOF events in several significant ways, with the following implications for future use of this approach.

Enhanced DEM accuracy: The use of high-resolution satellite imagery and accurate DEMs with a vertical accuracy of around 2 m allows for precise monitoring of terrain changes. This capability enhances the ability to predict potential GLOF events by identifying subtle deformations and volume changes in glacier lakes and surrounding terrains.

Volumetric calculations: Accurate volumetric calculations of debris slides and water depletion (e.g., 11.5014 million m³ of debris, 43.2172 MCM of water) provide critical data for assessing the potential impact of GLOFs. This quantitative analysis is essential for risk assessment and designing mitigation measures for future GLOFs.

Detailed terrain and water body metrics: Metrics such as lake perimeter, maximum length, width, and the area of erosion at the lake mouth help in understanding the scale and impact of GLOF events. These measurements can inform the development of models to predict the behavior of future GLOFs.

Identification of high-risk areas: The study's findings highlight specific regions with notable changes in elevation and volume patterns, which are critical for targeted monitoring and mitigation efforts. By focusing on these high-risk areas, resources can be allocated more efficiently to reduce the potential impact of future GLOFs.

Infrastructure planning: Detailed knowledge of the lake's dimensions and the dynamics of water depletion (e.g., 32 m water level drop) can inform the design and placement of infrastructure such as dams, drainage channels, and early warning systems to manage and mitigate flood risks. The study underscores the importance of continuous monitoring and technological advancements in satellite imagery and data analysis. Continued research into improving the accuracy and timeliness of GLOF predictions can drive innovation in remote sensing and geospatial technologies.

In conclusion, the novelty of this study is that it demonstrates the significance of utilizing stereo photogrammetric techniques for integrating satellite images from multiple sensors for generating pre- and post-GLOF DEMs to obtain insights into glacier lake breach dynamics, which no other paper in this context has addressed. The study has established the potential of satellite-based methods to capture a comprehensive view of glacier terrains and monitor minor changes that could lead to potential hazards. The work has employed stereo and single-image 3D photogrammetry which is first of its kind to estimate the photogrammetric volume of landslide debris, which was found to be 11.5014 MCM. Additionally, it is estimated that the volume of water emptied from the South Lhonak Lake after the breach to be 43.5 MCM. Cross-sectional changes are also derived at the lake outlet and land slide locations along with computing the depth of erosion to understand the phenomenon better. The stereo-image photogrammetry allows for multi-perspective analysis, enabling a 3D reconstruction of glacier lake landscapes before and after the event. Although the single-image 3D photogrammetric method is time-consuming due to manual mono-plotting and an iterative 3D Ray tracing process from ground cube coordinates along with refined RPCs was employed to derive the precise elevation in Part C of pre-DEM in the absence of a stereo pair. Overall, the study highlights the critical role of satellite-driven analysis in advancing the understanding of glacier lake breach dynamics, which can provide a foundation for cause and effect of glacier hazard.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.