Retrospective Analysis of Glacial Lake Outburst Flood (GLOF) Using AI Earth InSAR and Optical Images: A Case Study of South Lhonak Lake, Sikkim

Abstract

On 4 October 2023, a glacier lake outburst flood (GLOF) occurred at South Lhonak Lake in the northwest of Sikkim, India, posing a severe threat to downstream lives and property. Given the serious consequences of GLOFs, understanding their triggering factors is urgent. This paper conducts a comprehensive analysis of optical imagery and InSAR deformation results to study changes in the surrounding surface of the glacial lake before and after the GLOF event. To expedite the processing of massive InSAR data, an InSAR processing system based on the SBAS-InSAR data processing flow and the AI Earth cloud platform was developed. Sentinel-1 SAR images spanning from January 2021 to March 2024 were used to calculate surface deformation velocity. The evolution of the lake area and surface variations in the landslide area were observed using optical images. The results reveal a significant deformation area within the moraine encircling the lake before the GLOF, aligning with the area where the landslide ultimately occurred. Further research suggests a certain correlation between InSAR deformation results and multiple factors, such as rainfall, lake area, and slope. We speculate that heavy rainfall triggering landslides in the moraine may have contributed to breaching the moraine dam and causing the GLOF. Although the landslide region is relatively stable overall, the presence of a crack in the toparea of landslide raises concerns about potential secondary landslides. Our study may improve GLOF risk assessment and management, thereby mitigating or preventing their hazards.

1. Introduction

On 4 October 2023, a catastrophic glacial lake dam failure occurred at South Lhonak Lake in Sikkim, triggering a flash flood in the Testa Valley and causing extensive damage in northern Sikkim. The flood resulted in significant casualties, destruction of downstream infrastructure including bridges and roads, and widespread damage to residential properties. This disaster had a profound impact on various towns and infrastructure across Sikkim, with Chungtang Town being the worst affected, where over 80% of the area suffered damage.

Glacial lakes, formed by depressions carved out by glaciers (i.e., cirques) or moraines blocking water in ice valleys, are typically surrounded by unstable glacial moraines consisting of loose glacial till. The moraine dams of these lakes are inherently unstable, and dam break events occur under specific conditions, resulting in the sudden discharge of water and debris known as glacial lake outburst floods (GLOFs). Glacial lakes serve as sensitive indicators of climate change in alpine and cold regions and represent a potential source of disasters triggering mountain outburst floods or debris flows. Most glacial lakes are formed during glacier retreat [1,2]. Due to the serious consequences of GLOFs, there is an urgent need to gain a deeper understanding of the spatiotemporal evolution characteristics and triggering factors of GLOFs.

In recent years, analysis of triggering factors and discussion of the impact of moraine collapse on GLOF disasters have been on the rise [3,4]. Interferometric Synthetic Aperture Radar (InSAR) possesses all-day and all-weather monitoring capabilities, enabling the rapid acquisition of regional deformation information and facilitating regular and continuous observations. Consequently, it has emerged as a vital tool for monitoring the stability of surrounding slopes of glacial lakes. Numerous scholars have utilized D-InSAR, Permanent Scatter Interferometry (PSI) [5], Small Baseline Subsets InSAR (SBAS-InSAR) [6,7], offset tracking methods, and pixel offset (PO) tracking for small baseline subsets (PO-SBAS) to acquire surface deformation results around glacial lakes and assess their stability [8,9,10,11,12]. Yang et al. (2022) provides a systematic analysis of the contributing factors that lead to GLOF disasters through the utilization of meteorological data, SAR, and optical images [10]. Currently, the rapid advancement of SAR platforms has resulted in explosive data growth [13], and the field of InSAR image processing is shifting towards big data analysis. The influx of remote sensing big data presents both opportunities and challenges in data processing, management, and analysis [14]. Processing massive InSAR data within a limited time constraint poses a significant challenge [15].

The emergence and evolution of high-performance computing (HPC) presents a novel approach to tackling the challenge of large-scale InSAR processing. As a quintessential technology in HPC, graphics processing units (GPUs) have attracted increasing attention from researchers owing to their potent parallel processing capabilities, high memory bandwidth, low power consumption, and excellent performance in handling computationally intensive problems [13]. Cloud computing can more effectively tackle these challenges by offering on-demand access to a shared pool of computing resources, allowing users to leverage the capabilities of distributed computing and storage without requiring substantial upfront investments in hardware and software infrastructure [16,17,18]. We employed the AI Earth cloud platform to develop an InSAR processing system, which integrated the platform’s massive SAR remote sensing data and the advantages of cloud computing with our automated, high-performance InSAR processing algorithms [19], aiming at measuring surface deformation around glacial lakes efficiently.

This paper is organized as follows: Section 2 introduces the regional overview and disaster situation of the glacial lake studied in this paper. Section 3 elaborates on the data and methodologies employed in this study. Section 4 primarily focuses on observing and analyzing the InSAR deformation results alongside optical images of the landslide around the studied glacial lake. Section 5 discusses the correlation between InSAR deformation results and variables such as rainfall, lake area, and slope, while also speculating on the possible causes of this GLOF and the risk of secondary landslides.

2. Study Setting

The South Lhonak Lake (27°54′41.27″N, 88°11′29.38″E) is a moraine lake situated in the northwest of Sikkim, India. It stretches from east to west and is positioned on the southern Lhonak ice tongue in the Lhonak Valley, with an elevation of 5200 m (Figure 1). It stands out as one of the largest and fastest-growing glacial lakes in the region. In recent decades, global warming has led to an observable trend of retreating and thinning glaciers in high-altitude regions [20,21,22,23,24,25]. If this trend persists, it is projected that by 2100, approximately 80% of the current Himalayan glaciers will have vanished [26]. The rapid glacier retreat provides enormous space and sufficient water sources for the expansion of glacial lakes [27]. Sikkim, which has abundant hydraulic resources, experiences substantial rainfall from June to September, leading to frequent natural disasters such as landslides and floods [28]. Landslides triggered by earthquakes and heavy rainfall are common in this area [29].

The GLOF event studied in this paper occurred in the Teesta River basin area of northern Sikkim. A study by Raj and Remya et al. (2013) [30] suggested that heavy rainfall might induce GLOFs on moraine lakes. Over the 46 years from 1962 to 2008, the South Lhonak Glacier, the parent glacier of South Lhonak Lake, retreated by 1.9 km, leading to the formation of a moraine-dammed lake at the glacier snout. Subsequently, from 2008 to 2019, it further retreated by approximately 400 m. Concurrently with glacier retreat, the lake has exhibited significant expansion over the years, with its area escalating from 0.42 km² in 1990 to 1.35 km² in 2019 [28]. Sharma et al. conducted field investigations of South Lhonak Lake in 2014 and 2016, and based on on-site water depth, the research indicates that the storage capacity of South Lhonak Lake is 65.81 ± 2.5 million m³, with maximum and average depths of 131 ± 2.5 and 67.05 ± 2.5 m, respectively [31]. As lakes expand, the risk of catastrophic GLOFs escalates. In the Sikkim Himalaya, there are a total of 320 glacial lakes, of which 14 are potentially vulnerable to GLOFs [30]. A substantial inventory of glacial lakes and statewide hazard assessments [30,32,33,34] consistently identified South Lhonak Lake as a potential hazard with a high probability of eruption.

3. Data and Method

3.1. Data

3.1.1. SAR Images

To monitor the surface deformation of the slope surrounding South Lhonak Lake in Sikkim using InSAR technology, we utilized the Sentinel-1 SAR satellite with C-band, launched by the European Space Agency. A total of 83 images from ascending Track 12 and 79 images from descending Track 48 were acquired from Sentinel-1 in Interferometric Wide (IW) swath mode between January 2021 and September 2023, covering the study area before the disaster. Subsequently, 15 Sentinel-1 images from descending Track 48, collected between October 2023 and March 2024, were applied to discover the surface deformation after the disaster. We employed the Shuttle Radar Topography Mission (SRTM) DEM with a 30 m resolution as the reference DEM to remove the topographic phase.

3.1.2. Optical Images

To observe changes in the glacial lake area and its surrounding surface conditions, we utilized 35 Landsat 8/9 images without cloud cover, spanning from January 2021 to September 2023. Band combinations primarily involved bands 3, 5, and 7, with a spatial resolution of 30 m. These images were then employed to delineate the outline of the studied glacial lake and analyze alterations in its area. Additionally, to further scrutinize surface changes, we selected two GF-2 images with a spatial resolution of 4 m captured before and after the GLOF, specifically on 26 December 2020 and 27 October 2023.

3.2. InSAR Calculation in Cloud Platform

The time-series analysis algorithm of the InSAR processing system in this study is built upon the SBAS algorithm, with the main workflow illustrated in Figure 2. This system integrates InSAR GPU-assisted processing technology and a fast InSAR time-series analysis method with full resolution [19]. The massive InSAR data processing system operates on a cloud platform that provides open and free online services available at https://aiearth.aliyun.com (accessed on 3 April 2024) [14]. Users can directly utilize the system in a web browser without the need to download, install any software, or configure complex software environments. Furthermore, the utilization of cloud computing and distributed processing technologies ensures the efficient and reliability of large-scale InSAR processing, with low hardware requirements. The user interface is easy to use without requiring programming skills or specialized knowledge of InSAR. This allows researchers from various disciplines to conveniently utilize the system for further analysis. Through a web-based Graphical User Interface (GUI), users can select the area of interest (AOI), configure start and end times for data processing, and adjust basic interferometry processing parameters easily and customarily. These interferometric parameters include maximum and minimum temporal baselines, valid value ratios of corresponding points in the interferograms, and filtering parameters. The system optimizes the InSAR processing workflow by simplifying processing parameters, which reduces human intervention and increases processing efficiency. InSAR processing and deformation result displays from massive Sentinel-1 data can be promptly completed on the cloud, with exportable computed results for further analysis. This approach addresses the challenges in traditional InSAR technology, such as data downloading, storage, parameter setting, and limitations in computing resources, thereby meeting the demands of diverse natural disaster surveys and emergency monitoring efforts.

3.2.1. GPU-Assisted InSAR Processing Module

The InSAR data processing involves two primary modules: interferometry processing and InSAR time-series analysis. This study utilizes GPU parallel technology developed under the CUDA framework to optimize the key time-consuming steps in interferometric modules. These steps include geometric coregistration, resampling, and enhanced spectral diversity (ESD) [35]. The Sentinel-1 SAR satellite constellation employs terrain observation using the progressive scans (TOPS) mode to acquire IW SLC products. In this mode, the varied squint angle in the azimuth direction introduces an azimuth-varied Doppler centroid frequency, imposing stricter coregistration requirements. An inaccurate coregistration leads to phase jumps between bursts. Therefore, a very high azimuth registration accuracy is needed [36]. However, precise coregistration methods have high computational complexity. To attain this level of accuracy, distinct coregistration methods for initial and fine coregistration are necessary, typically utilizing the precision orbit-based geometric positioning method and the ESD method, respectively. Geometric coregistration entails determining the azimuth/range offset between an image pair by back-projecting DEM into master and slave image coordinates [35]. Subsequently, the slave image is resampled with the master image, followed by the utilization of ESD to refine the azimuth shifts derived from geometric coregistration, thereby enhancing registration accuracy further. ESD leverages the phase difference in the overlapping area of master and resampled slave images between adjacent bursts within a subswath to estimate misregistration residuals [36,37].

The CUDA implementation of geometric coregistration is structured into three stages [35]: data initialization, offset estimation, and warp function estimation. During the data initialization, both the DEM and orbit-fitted polynomials are transferred to GPU global memory. The GPU constant memory is utilized to store the variables, such as the bounding box and intervals of the DEM product. In the offset estimation stage, DEM pixels are radar coded into master and slave image coordinates on the GPU. Finally, the least squares (LS) fit is performed on the GPU. Resampling is achieved by two kernels: DerampDemod and Interpolation. Before kernel execution, the slave image and deramping parameters are transferred into the GPU global memory. After calling the DerampDemod kernel, the estimated phases and the deramped slave image are stored in the global memory and subsequently used by the Interpolation kernel. During 2D interpolation, the interpolation convolution kernels and warp functions are transferred into the GPU global memory. In the ESD approach, overlaps, along with the image azimuth sampling frequency and Doppler centroid frequency difference, are transferred into GPU global memory.

3.2.2. Automated Full-Resolution Fast InSAR Time-Series Analysis Method

The traditional SBAS approach may reduce the coverage of coherent points in some cases due to the multi-look procedure [7,38]. Moreover, it is susceptible to losing information when the images have rich details or the spatial deformation signal becomes complex [39]. To address these limitations, the full-resolution InSAR method utilizes single-look data to circumvent the multi-look procedure. It performs time-series analysis directly in the SAR coordinate system, ensuring both high processing efficiency and precise monitoring results. Additionally, the statistically homogeneous pixel selection (SHPS) algorithm is used to filter the differential interferograms. The workflow of the full-resolution fast InSAR time-series analysis method for processing Sentinel-1 data is illustrated in Figure 3.

• Employ the small baseline principle to select interferometric pairs and generate the optimal interferometry network [40].
• Calculate burst offsets between each image and the reference image, generating a burst offset file and determining the burst offsets of each slave image based on the AOI of the reference image.
• Automatically download the corresponding orbit auxiliary files and external DEM files. SRTM DEM with a resolution of 30 m was utilized to subsequently mitigate terrain phase effects.
• Utilize GPU to accelerate the generation of differential interferograms; details of GPU-accelerated InSAR processing are available in Section 3.2.1. Subsequently, all generated differential interferograms are resampled based on the registration parameters to ensure consistency with the SAR coordinate system of the reference image.
• Image cutting. Interferograms are cropped according to the specified range of the AOI.
• SHPS phase filtering and phase unwrapping. Utilize the SHPS algorithm to reduce noise in the interferograms while preserving the spatial resolution of SAR images. Coherent points surrounding each reference pixel are selected, aiming to retain interferogram details while eliminating phase noise from incoherent and low-coherence areas. Then, phase unwrapping of interferograms was achieved using minimum cost flow (MCF) networks [41].
• Corrections for orbital error and terrain-related atmospheric delay errors.
• Time-series analysis in SAR coordinate system. With high-pass and low-pass filters, the average deformation rate is calculated using the linear least squares (LS) method. Subsequently, a time-series analysis is performed. The InSAR time-series analysis module follows the traditional method, employing the Small Baseline Subset method to derive deformation time series through the singular value decomposition (SVD) algorithm [6].

4. Results and Analysis

4.1. Analysis of InSAR Deformation Results

Figure 4 shows the spatial-temporal baseline distribution, which was generated according to the small baseline principle described in Section 3.2. This paper limited the temporal baseline to 5–36 days to keep high coherence of interferograms. The Sentinel-1 satellite exhibits high orbit control accuracy, with its orbit tube diameter being significantly smaller than the critical spatial baseline. Therefore, we did not constrain its spatial baseline. Before the disaster, 240 interferograms were created from ascending track 12 and 220 from descending track 48; after the disaster, 39 interferograms came from descending track 48. The mean deformation rates and time-series deformation of South Lhonak Lake were measured via line-of-sight (LOS), as depicted in Figure 5a,b. The positive values indicate displacement towards the satellite, while negative values indicate movement away; blue signifies positive displacements while red shows negative. The underlying image is a GF-2 optical image taken on 26 December 2020.

The InSAR analysis has revealed the presence of two distinct zones of deformation within the moraine encircling South Lhonak Lake. A significant deformation zone, which is positioned in the northern sector of the moraine, displayed persistent deformation patterns well in advance of the catastrophic event. This area exhibited a gradual and sustained subsidence along with localized ground displacements, indicating a protracted period of instability. Figure 5a illustrates that the maximum deformation rate derived from the ascending images is about −80 mm/yr, whereas the overall deformation rate from the descending images is slightly lower (Figure 5b), with a maximum deformation rate of around −65 mm/yr. The significant deformation area results from both ascending and descending images generally align with the landslide area observed through optical images. In the deformation results from the ascending images, the area with the highest mean deformation rate is within the landslide area, while a portion of the area with the highest mean deformation rate from the descending images lies to the east of the landslide area. The significant deformation area from the ascending images is more congruent with the landslide area.

To further analyze the deformation area of the landslide, four points (P1 to P4) were selected and shown in Figure 5. The time-series deformation curves for these points, obtained from the ascending and descending images, are presented in Figure 6a,b, respectively. It is evident that the overall deformation trend of the selected points from both ascending and descending orbits indicates relatively small deformation signals before January 2022, with relatively flat curves. However, starting from 2022, the deformation rate markedly increased. The inflection point indicating acceleration in the ascending images is around January 2022. Deformation results derived from the ascending images reveal that the deformation on the eastern side of the landslide area is most pronounced, with a deformation rate exceeding −60 mm/yr. Among the four selected points mentioned, P1 exhibits the highest cumulative deformation, approximately −170 mm, with an average deformation rate of −70 mm/yr (Figure 6a). Deformation results from the descending images indicate two significant acceleration periods at points P1 and P2 from January 2022 (Figure 6b). Overall, deformation at P3 is minimal, showing an indistinct trend in deformation rate change.

4.2. Optical Image Analysis

From the GF-2 images in Figure 7, South Lhonak Lake exhibits a spindle-shaped form, gradually tapering from west to east. Positioned to the west of the lake is its parent glacier, the South Lhonak Glacier. The dashed red line delineates the landslide area, while the dashed blue line outlines the lake boundary of South Lhonak Lake before the GLOF. Significant changes occurred in the glacial lake area during this period due to the substantial time gap between the two GF-2 images. Consequently, Landsat 8 images from 25 June 2023 were utilized to delineate the lake outline before the disaster. By comparing the lake extent in Figure 7a with the boundary marked by the dashed blue line, we can visually discern the changes in South Lhonak Lake’s area from the end of 2020 to the middle of 2023. In Figure 7b, post-GLOF observations indicate a decrease in lake area and a significant drop in water level. The yellow circle highlights the breach in the glacial moraine dam, with evident erosional marks visible after the breach. A comparison of optical images before and after the GLOF reveals conspicuous signs of moraine landslide at the junction of the South Lhonak glacier’s terminus and the glacial lake. Measurements indicate that the axial length of South Lhonak Lake is approximately 2800 m, with a maximum width of about 700 m. We further enlarged the GF-2 images of the landslide area to observe more subtle surface changes, as depicted in Figure 8. According to measurements, the landslide area spans approximately 1000 m in length and 360 m in width, covering an area of approximately 0.29 km². In Figure 8b, a clear and nearly parallel tensile crack is visible at the trailing edge of the landslide. To the east of the landslide area, there is a gully formed by melting water from alpine ice and snow.

5. Discussion

5.1. Correlation between InSAR Deformation Results and Multiple Factors

5.1.1. Rainfall Factor

Figure 9 illustrates the monthly average rainfall in North Sikkim from 2018 to 2022 and the five-year average rainfall. Notably, the rainy season typically begins around May each year, with May to September experiencing the highest monthly rainfall, often exceeding 300 mm. The InSAR analysis in Section 4 indicates that during the monitoring period, deformation patterns obtained from descending images show a significant acceleration at point P1 in May 2023, coinciding with the onset of the rainy season. In April 2022, rainfall in North Sikkim significantly surpassed previous Aprils, reaching levels comparable to those usually seen in May, exceeding 400 mm/month. This suggests an earlier onset of the rainy season in 2022 compared to previous years. Concurrently, as depicted in Figure 6b, there was a notable acceleration in deformation at points P1 and P2 in April 2022, occurring relatively earlier compared to the acceleration observed at point P1 in 2023. Moreover, the average monthly rainfall in September and October 2022 was high, both exceeding 400 mm/month, thereby extending the rainy season compared to previous years. The accumulated deformation and duration of accelerated deformation from April to October 2022 at points P1–P4 over three years, spanning from 2021 to 2023, are the highest, likely due to the prolonged rainy season in 2022. Hence, we believe there is a correlation between the landslide area’s deformation rate and rainfall patterns.

5.1.2. Lake Area Factor

Figure 10 presents Landsat optical images of the South Lhonak Lake region spanning from January 2021 to October 2023, with the final subgraph obtained on 15 October 2023, after the GLOF. These images visually and clearly depict the changes in the lake over time, illustrating the evolution of the glacial lake. Figure 11 depicts the correlation between the time-series deformation of point P1 derived from the descending images before the disaster and the variation in the area of South Lhonak Lake from January 2021 to July 2023. Through Figure 11, it is evident that the variation in glacial lake area exhibits periodicity. Around October each year, the lake area begins to decrease rapidly, reaching its minimum during winter. Subsequently, from approximately January, the area starts to expand rapidly, while during the months of April to September, the lake area remains relatively stable. Overall, the lake area shows a yearly increasing trend. At the same time, point P1 exhibits periods of rapid deformation from January to May 2022 and from February to May 2023 (indicated by the gray area in Figure 11), with deformation rates higher than those observed during other periods. It was observed that during 2022 and 2023, the lake area experienced rapid growth from January to March, corresponding to a significant increase in the deformation rate of P1 from January to May. However, the onset of deformation acceleration slightly lagged behind the rapid expansion of the lake area. Therefore, we believe there is a correlation between the deformation rate of P1 and the variation in the area of the glacial lake. Furthermore, the relatively high rainfall in April and May may have contributed to the accelerated deformation observed during these months.

5.1.3. Slope Factor

Figure 12 illustrates the slope map of the study area computed using NASADEM HGT v001. The dashed blue line represents the pre-disaster lake contour outlined using Landsat 8 images acquired on 25 June 2023, while the solid red line delineates the extent of the landslide. Analysis of deformation results derived from the ascending images reveals a significant acceleration at point P1 from May to July 2021, surpassing the velocities observed at other selected points. Notably, P1 exhibits the highest deformation rate among the four selected points (P1–P4). Furthermore, the slope map depicted in Figure 12 highlights that within the landslide area, the region housing point P1 exhibits the steepest slope.

From Figure 6, it is evident that the accumulated deformation in the vicinity of P1 and P2 exceeds that of P3 and P4, possibly due to the presence of a gully near the area of P1 and P2. Previous research has indicated that river erosion plays a pivotal role in triggering landslides [42,43], leading us to speculate that water flow erosion contributes to the substantial deformation observed in this area.

5.2. Possible Causes of Landslide and GLOF

The occurrence of GLOF may be affected by many factors, and the triggering mechanisms are intricate. Emmer et al. (2013) summarized the causes and mechanisms of moraine-dammed lake failures, identifying five dynamic factors and three long-term factors. The former includes earthquakes, slope movements into the lake, floods from upstream lakes, blockage of underground outflow channels, and intense rainfall or snowmelt. The latter can be considered as spontaneous “dam self-destruction” resulting from melting buried ice, hydrostatic pressure, and temporal effects [44,45,46].

The relationship between the circular area of landslides affected by earthquakes and the magnitude (M > 4) can be expressed as [47]:

log₁₀ A = M − 3.46 (±0.47)

(1)

Among them, A represents the affected area (km²), and M denotes the earthquake magnitude. We utilize the following formula to calculate the maximum critical distance R, in kilometers, of landslides affected by earthquakes:

$$R=\sqrt{\frac{{10}^{(M-2.99)}}{\pi }}$$

(2)

By substituting M = 5 into Formula (2), the maximum critical distance is calculated to be approximately 5.7 km. According to seismic records from the United States Geological Survey (USGS) (http://earthquake.usgs.gov/earthquakes accessed on 13 April 2024), no earthquakes of magnitude 5.0 or higher occurred within a 100 km radius of South Lhonak Lake from January 2023 to October 2023. Thus, it is unlikely that the South Lhonak Lake outburst studied in this article was caused by an earthquake. Furthermore, Sattar et al. (2021) [28] identified the cliffs around South Lhonak Lake most likely to be prone to avalanches based on mass movement potential. No avalanches were found through optical imaging in this area before the GLOF, essentially ruling out the possibility of the moraine-dammed lake failure being due to avalanches.

Numerous studies have highlighted rainfall as one of the primary triggers of landslides in many regions worldwide [48], due to its significant impact on slope stability [49,50]. The analysis in Section 5.1.1 also indicates a rapid increase in the slope’s deformation rate around the moraine-dammed lake as the rainy season commences. Sikkim received 101 mm of rainfall in the first five days of October 2023. Between 3–4 October alone, the rainfall in the state escalated to five times the normal level. On 4 October, daily rainfall across Sikkim reached 115–205 mm, coinciding with the occurrence of the moraine-dammed lake failure. Observation of Landsat 8/9 images from just before and after the GLOF revealed no significant landslides in the moraines on 29 September. However, on 7 October, a notable landslide in the moraine was visible. Therefore, during the period of sustained heavy rainfall from 29 September to 7 October, a large-scale landslide occurred in the moraine area, characterized by significant slope and cumulative deformation. Consequently, we attribute heavy rainfall as the primary probable cause of the landslide in this study. Simultaneously, the GLOF disaster occurred in South Lhonak Lake on 4 October 2023. In summary, we speculate that the landslide in the moraine induced by heavy rainfall may be one of the possible reasons for the eventual breach of the terminal moraine dam. Under such climate conditions, glaciers will persist in retreating and melting. Consequently, over time, the glacial lake area may revert to its original volume after the outburst, posing a potential hazard as a glacial lake susceptible to external drivers, potentially leading to recurrent GLOF events.

5.3. Secondary Landslide Risk

The LOS average deformation rate (Figure 13) and time-series deformation (Figure 14) around South Lhonak Lake after the GLOF were computed using the Sentinel-1 SLC data from October 2023 to March 2024 from the descending track. Significant deformation was observed both in the landslide area and its eastern vicinity, with the deformation within the landslide area notably surpassing that in its eastern counterpart. Six points were selected within the zone of significant deformation, comprising P5–P7 within the landslide region and P8–P10 situated to the east of it. The average slope of the eastern region, as depicted in Figure 12, is considerably smaller than that of the landslide area. From Figure 13, it is apparent that the average deformation rate surrounding the crack at the trailing edge of the landslide is relatively small. In Figure 14, the deformation trend of P5–P7 appears relatively consistent, without discernible points of acceleration. This suggests that the larger deformation observed in the landslide area may be attributed to the movement of surface debris and floating soil, indicating that the landslide zone is relatively stable. The cumulative deformation of P8–P10 is smaller than that of P5–P7, and since the end of December 2023, the deformation rate has noticeably decreased, with the curve remaining relatively flat, indicating a current state of relative stability in this region as well. Upon observing the Landsat 8 optical images after the disaster, no prominent secondary landslides were detected in the landslide area, and there were no observable landslide occurrences in the eastern area, aligning with our initial inference.

However, Yin et al. (2023) [51] pointed out that the cracks in the top area of landslides serve as significant indicators of the instability of the landslide headscarp and have a potential role in landslide hazard warning. Hence, a risk of secondary landslides in the area is inferred. The volume of earthwork for potential secondary landslides was measured and estimated using GF-2 images, resulting in approximately $2.5\times {10}^6$ m³. Future monitoring of surface deformation in this area will continue.

6. Conclusions

In this study, a cloud-based InSAR processing system was utilized to analyze surface deformation around South Lhonak Lake in Sikkim, India, before and after the GLOF event. Our findings reveal that a moraine landslide occurred at the junction of the South Lhonak glacier’s terminus and the lake, with observable deformation characteristics preceding the landslide. There exists a certain correlation between InSAR deformation results and variables such as rainfall, lake area, and slope. Notably, the surface deformation rate in the landslide area accelerated at the onset of the rainy season, with both the accumulated deformation and the duration of accelerated deformation showing correlation with the duration of rainfall. Soon after the area of South Lhonak Lake began to increase rapidly each year, the deformation rate of the points within the landslide area subsequently increased significantly. Additionally, areas with steeper slopes exhibited heightened susceptibility to accelerated deformation induced by heavy rainfall, with water flow erosion potentially contributing to moraine landslides. Through comprehensive analysis, we inferred that heavy rainfall-induced landslides in the moraine could be one of the primary factors leading to a breach in the terminal moraine and triggering the GLOF. Analysis of post-disaster InSAR deformation results and optical images suggests that significant deformation in the landslide area was primarily caused by the movement of loose soil and debris on its surface, with a relatively stable overall condition of the landslide region. Although deformation is evident in the eastern vicinity of the landslide area, it remains relatively stable at present. However, the presence of cracks in the top area suggests a potential risk of secondary landslides. Future research may leverage higher-resolution optical images, SAR data from various bands (e.g., Lutan-1 SAR satellite data and NISAR satellite data), and meteorological data to further investigate GLOF occurrences. Such efforts can mitigate the impact of GLOF events and inform infrastructure development planning.