Research ArticleGEOPHYSICS

Seismic observations, numerical modeling, and geomorphic analysis of a glacier lake outburst flood in the Himalayas

See allHide authors and affiliations

Science Advances  16 Sep 2020:
Vol. 6, no. 38, eaba3645
DOI: 10.1126/sciadv.aba3645


Glacial lake outburst floods (GLOFs) are a substantial hazard for downstream communities in vulnerable regions, yet unpredictable triggers and remote source locations make GLOF dynamics difficult to measure and quantify. Here, we revisit a destructive GLOF that occurred in Bhutan in 1994 and apply cross-correlation–based seismic analyses to track the evolution of the GLOF remotely (~100 kilometers from the source region). We use the seismic observations along with eyewitness reports and a downstream gauge station to constrain a numerical flood model and then assess geomorphic change and current state of the unstable lakes via satellite imagery. Coherent seismic energy is evident from 1 to 5 hertz beginning approximately 5 hours before the flood impacted Punakha village, which originated at the source lake and advanced down the valley during the GLOF duration. Our analysis highlights potential benefits of using real-time seismic monitoring to improve early warning systems.


Glacial lakes in the Himalayas are rapidly growing due to climate change and acceleration of glacier melt in recent decades (1, 2). Many of these lakes are dammed by unstable, often ice-cored moraines (3), surrounded by steep topography prone to landslides, and frequently subjected to seismic events and intense monsoonal precipitation. Risks associated with glacial lake outburst floods (GLOFs) are substantial for downstream inhabitants in these regions, and thousands of fatalities have already occurred as a direct result of sudden catastrophic releases of water (4). Societies and economies of GLOF-prone regions are severely impacted, including destruction of infrastructure, disruption to communities, and loss of life (25). The nation of Bhutan, in particular, has been classified as having some of the greatest national-level economic consequences of glacier flood impacts, as hydropower dominates the nation’s gross domestic product and socioeconomic development potential (4, 6). The electric power demands of Himalayan nations are on a steep rise with rapid economic growth, and hydropower development continues to expand into higher sites closer to glaciers (2).

Numerous GLOFs have been documented in the Himalayas by firsthand observation and satellite imagery after their occurrence (79). The most comprehensive analysis to date of the Landsat imagery archive finds an average GLOF frequency in High-Mountain Asia (HMA) of 1.3 GLOFs per year since the 1980s (9). Despite this steady rate, quantitative in situ observations of GLOFs are scarce (3) and limited to rare situations where preinstalled instruments are coincidentally located in the same valley as the GLOF (10). Efficient strategies to mitigate outbursts, design reliable early warning systems, and minimize destructive impacts would benefit from continuous observation of flood evolution through time. While satellite observations can help identify regions of high GLOF risk and quantify geomorphic impacts after occurrence, they cannot capture GLOF events in real time. Numerical flood models can be used to simulate dam outbursts and flood waves (7, 1113) yet require many physical parameters as input, which are often unknown or poorly constrained (12). For example, the shape of a breach hydrograph can strongly influence flood characteristics downstream (14). These essential upstream boundary conditions in flood routing simulations are usually unknown and must be estimated using physical or parametric dam breach models with large uncertainties (15). Here, we demonstrate how flood models can be tested using continuous time-stamped seismic observations of GLOF mass movement, thus helping to determine whether estimates of downstream flood arrival times are realistic.

The Pho Chhu (river) flows from the pristine mountain region of Lunana [4500 meters above sea level (masl)], southward through the Bhutan foothills where it eventually joins the Brahmaputra River in northern India (Fig. 1). The valley has a typical step-like elevation profile, alternating between sections of very steep and relatively flat terrain separated by river knickpoints. Fed by seasonal snow melt, glacier melt, and summer monsoon rains, the river plays a vital role in the welfare and livelihood of the people of Bhutan. In the Pho Chhu valley, residents of small villages live and farm along the banks of the river, with rice paddies being a primary seasonal crop. In the river source region, several large proglacial lakes are continually expanding in size as a result of accelerated glacier melting in recent decades (1, 2). Approximately 90 km downstream from the proglacial lakes is the village of Punakha (1200 masl), where a large 17th century Buddhist temple (the Punakha Dzong) is situated along the riverbank at the confluence of the Mo Chhu and Pho Chhu. This temple and surrounding region are historically and culturally important, and there is a great deal of concern about catastrophic outburst floods from the unstable lakes above.

Fig. 1 Region of study.

Lugge Tsho was the source of the 1994 GLOF event, while Raphstreng Tsho and Thorthormi Tsho are also considered high risk for future GLOFs. Locations of the five seismometers used in the study are denoted by the yellow circles in (A). (B) Location of the Punakha and Wangdue villages downstream of the GLOF source area. (C) Unstable glacial lakes in the Lunana region. (D) Mid-Holocene glacier moraines that were left largely intact by the GLOF. (E) Lugge Tsho and moraine breach zone.

Such an event occurred on the night of 7 October 1994, when the moraine dam of Lugge Tsho breached and debris-laden flood waters surged down the Pho Chhu valley. Twenty-one lives were lost; and the flood destroyed an estimated 12 houses, 5 water mills, 816 acres of crops, 965 acres of pasture land, 16 yaks, 6 tons of stored food grains, 4 bridges, 2 stupas, and damaged part of the Punakha Dzong (16). While the exact trigger of the Lugge Tsho breach is unknown, various causes of moraine failure have been hypothesized such as the melting of its ice core (8), a gradual increase in hydrostatic pressure as the lake depth increased due to melting (5), or collapse of part of the right lateral hillside into the lake, causing a sudden increase in hydrostatic pressure (17). After the tragic incident, much effort was put forth by the government of Bhutan to assess the risk of future GLOFs in the region, establish an emergency warning system, artificially lower lake water levels, and study Lugge Tsho in more detail. Post-GLOF field investigations found a mean lake depth of approximately 50 m, a lake volume of 58.3 million m3, and a typical discharge at the lake outlet varying from 2.5 to 5 m3 s−1 during September and October 2002. The total volume of water released during the GLOF was estimated by Yamada et al. (18) as 17.2 ± 5.3 million m3 based on a differential Global Positioning System (GPS) survey of the lowering of the lake level by 16.9 ± 3.2 m, while an integration of the Wangdue station hydrograph (fig. S1) during the GLOF yields approximately 25 million m3 (19). In addition, dam breach, flood propagation, and debris flow models (sediment-water mixtures) have also been used to simulate GLOF scenarios in this region, including sequential and simultaneous (as may happen during a large earthquake) breaches occurring below high-risk lakes in Bhutan. These models suggest that downstream villages including Punakha and the major portion of Wangdue Phodrang are at risk for severe inundation if another large GLOF occurs (1113, 19, 20).

Across the Himalayan region as a whole, the destructive potential of GLOFs is increasing. Yet the scarcity of observational measurements deters robust validation of numerical simulations, hinders quantification of GLOF dynamics, limits real-time warning, and increases uncertainties regarding societal impacts. Continuous observational coverage offered by seismic monitoring is one potential avenue for addressing this problem. Displacement of mass at Earth’s surface generates elastic seismic waves, which carry information about the source and can be recorded by seismometers at high temporal resolution across large spatial scales (21). Proof-of-concept studies have already shown the potential of seismic monitoring for diverse types of surface activities including river bedload transport and debris flows (2124) and further demonstrate the ability of seismic records to specifically provide insight into flood mechanics (10, 23). Here, we extend the application of seismic data to a Himalayan GLOF using data from the International Deep Profiling of Tibet and the Himalaya (INDEPTH) II experiment (25, 26). These data were collected by a passive broadband seismic array situated on the Tibetan Plateau, which was coincidentally recording when the GLOF occurred in Bhutan in 1994. Insofar as the authors are aware, INDEPTH II is the only seismic data available from 1994 in the region (within a 150-km radius). We perform a time-frequency analysis of the seismic signal produced by the GLOF and use cross-correlation functions (CCFs) between seismic stations to locate and track the source of coherent seismic energy through time. With the seismic data, firsthand accounts, and gauge station measurements, we constrain the progression of the flood from initial outburst to arrival in populated villages using a numerical flood model. To further quantify geomorphic impacts of the GLOF, we apply historical spy satellite images, Landsat, and modern high-resolution imagery to analyze lake area changes over time, the extent of flood-deposited sediments, the rate of vegetation regrowth post-GLOF, and pre- and post-flood morphology of the source area dated by 10Be.


Analysis of seismic data

Seismic energy generated during the GLOF was recorded by five seismometers (with locations ranging from approximately 75- to 130-km distance from the breach point) as a clear high-frequency (1 to 4 Hz) signal lasting several hours (Fig. 2 and fig. S2). Seismic energy at these frequencies was most likely excited by energy from turbulent flow and bedload transport processes, as observed in previous studies of seismic signals generated during high-flow regimes (10, 2123, 27). The GLOF signal strength was 5 to 15 dB above typical background noise levels (fig. S10) and occurred during the night and early morning when local anthropogenic noise was at a minimum. We also note several earthquakes (with epicenters in the Kuril Islands and Chile) and a nuclear explosion (Southern Xinjiang, China) that were recorded by the seismic array during these hours (Fig. 2A). Spectrograms show a similar pattern on vertical and horizontal components across all five stations (fig. S2). Some interstation variability in peak frequency content is observed during the two GLOF phases, likely due to lateral heterogeneity in attenuation structure and differences in distance from the source. The first detectable high-frequency signal arrived at approximately 1:45 a.m., beginning with relatively weak amplitude and limited primarily to frequencies between 1 and 3 Hz. Over the next several hours, the seismic energy varied somewhat through time, as the flood wave passed through sections of the valley with different slope and river channel characteristics. Around 25 min after this first arrival, an increase in spectral amplitude occurred across a wider range of frequencies (1 to 4 Hz) (fig. S2) and then subsequently tapered off over the next 1.5 hours. At approximately 3:50 a.m., the signal power began increasing again and reached a maximum at 6:00 a.m. with frequency content ranging from approximately 1 to 5 Hz. During this interval, the flood wave passed through the main branch of the Pho Chhu and impacted Punakha at approximately 7:00 a.m. based on eyewitness accounts.

Fig. 2 Key seismic observations on 7 October 1994.

(A) Example of a seismic trace (station BB20 BHZ) band-pass filtered between 1 and 5 Hz and spectrogram from 1 to 5 Hz during the GLOF duration. Black arrows mark seismic events that correlate in time with events in the USGS catalog, and gray arrows mark potential local earthquakes that appear on all five stations (fig. S2) but are not in the catalog. mb, body wave magnitude; Mw, moment magnitude. (B) Twelve-hour time series of signal coherence for the day containing the GLOF event (red; 1994-10-07) and the day before the GLOF event (gray; 1994-10-06). The top and bottom bounds of each shaded region represent the 95th and 75th percentile of coherence values across all lag times, respectively. The blue shaded regions mark 90 min of most significant GLOF signal relative to background levels. (C) Migrations of the CCFs during the two intervals using a velocity of 3.0 km s−1 (likely indicating regional short-period Rayleigh waves; fig. S3) for several station pairs is summed together to form these final coherence maps (see Materials and Methods and fig. S3). These two images illustrate the downstream progression of seismic energy generated by the GLOF. The coherence map on the left corresponds to seismic energy generated during the initial GLOF breach, while the map on the right corresponds to seismic energy generated approximately 4 hours later and ~70 km downstream. Times are Asia/Thimphu local time [universal time coordinated (UTC) +6].

We further examined the signal correlation across stations to constrain the origin of the seismic energy in space and time (21, 28, 29). Computing CCFs for every station pair across a series of frequencies ranging from 1 to 5 Hz, we find strong correlation of waveforms during two distinct intervals (Fig. 2B). The first spans from approximately 1:45 to 3:15 a.m., during which a strong peak in CCF amplitude is apparent (fig. S3). Migrating the CCFs during this interval and subsequently summing them together, a region of high coherence emerges, focused directly on the GLOF breach location at Lugge Tsho (Fig. 2), indicating that during this time, the outburst event was the dominant source of seismic energy at these frequencies. Approximately 4 hours later, a second prominent interval of high coherence spans from around 5:45 to 7:15 a.m. Migration of the CCFs during this later interval indicates that the GLOF-induced seismic energy originated from a lower (downstream) section of river, indicating that the flood wave had reached this point in the valley (approximately 70 km from the breach and 20 km above Punakha) by around 5:45 a.m.

Seismic signal generation

Generation of seismic waves from fluvial processes is understood to occur by two main processes: (i) transport of sedimentary grains that stochastically impact the riverbed (30) and (ii) turbulent fluid flow that interacts with the riverbed (31). Previous observations of seismic energy associated with bedload transport and turbulent flow processes made at much smaller distances from the seismic source (typically only hundreds of meters to a few kilometers away) show peak frequencies ranging from 1 to 100 Hz depending on distance from the source and local seismic attenuation structure (10, 2123, 27). Here, we demonstrate that coherent 1- to 5-Hz seismic energy is generated and can propagate to distances of ~100 km from the source region during large GLOF events.

Numerical models that predict seismic energy excitation due to turbulent flow and bedload transport processes suggest that >1-Hz seismic energy is difficult to produce at such large distances because of seismic attenuation of the high frequencies (3032). However, the high river flow rates (~2500 m3 s−1) and thick water flow depths associated with the GLOF represent extreme conditions that have not been explored in detail using numerical models and may well violate their underlying physical assumptions. A GLOF scenario of this magnitude may contain physics not captured by bedload transport and turbulent flow processes alone, such as flexure of the riverbed due to the initial flood wave. Additional work is needed to fully understand the physical processes responsible for the observed seismic signal generation during the GLOF event and how these observations compare with predictions from recent numerical models. In this study, we attribute the primary seismic signal to the initial flood wave as it proceeded down the valley.

Constraining a flood model

Hydropower viability, disaster preparedness, and paleoseismic investigations have previously simulated flood events in this region using numerical models (1113, 19, 20). Here, we build on results from these earlier studies by calibrating a flood model using the new set of independent observations: (i) an estimated start time and location based on the beginning of detectable seismic energy and migrated CCFs, (ii) the second interval of correlated seismic energy that occurred several hours later and ~70 km downstream of the initial breach, (iii) the approximate arrival time of the flood at Punakha from firsthand observations, and (iv) the flood hydrograph from Wangdue station (around 110 km downstream from the breach). Together, these provide key constraints on the location of the flood wave through the duration of the GLOF and allow us to parameterize a model with a higher degree of confidence than previously possible. The start time of the GLOF, in particular, is a key aspect that was previously unknown. Because of the sensitivity of flood models to input parameters (fig. S1), these independent observations are especially useful for validating and selecting model runs that are most realistic (see Materials and Methods).We use the U.S. Army Corps of Engineers Hydrologic Engineering Center’s River Analysis System (HEC-RAS) software to perform a series of two-dimensional (2D) unsteady flow simulations; select model runs that agree with the independent observations within a ±30-min threshold (see Materials and Methods and Fig. 3); and report ranges of simulated breach-to-arrival times for the main populated villages along the river valley for Thanza (0.4 to 0.6 hours), Tshojo (1.0 to 1.3 hours), Lhedi (1.4 to 1.8 hours), Samdingkha (3.9 to 4.6 hours), Punakha (4.4 to 5.2 hours), and Wangdue (5.7 to 6.5 hours) (Table 1). We note that around 6:30 a.m., the peak flows in the best-fit model runs precede (in time) the peak in seismic coherence, approximately 70 km downstream from the breach. While this may reflect the actual order of events, the model may also be slightly overestimating the flood wave velocity within this part of the valley. In all model runs, the duration between arrival and peak flow gradually decreases as the flood travels down the section of valley above Punakha. This consistent aspect of the flood illustrates the manner in which topography can shape a flood wave, in this case, causing the wave to become steeper and more prominent while traveling down a steep valley gorge. Furthermore, a reduction in hillslope angle of the arriving flood wave correlates in time with the second peak in coherent seismic energy at approximately 5:45 to 7:15 (fig. S13), perhaps indicating that changes in topographic slope influence the degree to which flood wave energy is transferred into the solid earth as seismic waves. Results from model runs with different breach hydrographs all converge to a similar shape before reaching the larger villages downstream. In the village of Samdingkha, for example (~8 km above Punakha), the duration between first arrival and peak flow ranges from 10 to 30 min. Without an external early warning system, this sudden rise in water level permits only a short time for inhabitants to move to safe ground, particularly if the early stages of the rise go unnoticed for several minutes.

Fig. 3 Results from the HEC-RAS 2D unsteady flood model.

(A) Elevation profile of the river valley, with hourly locations of flood arrival time from the best-fit model run. (B) Distance from the moraine breach versus time. The gray square symbols and brackets are independent observations from seismic, eyewitness, and gauge station sources. The color-shaded regions represent the range of model outputs that match observations within ±30 min, and the colored curves represent the single best-fit model run. The orange curve is the simulated arrival of the flood wave, the blue curve is the peak flow, and the yellow curve indicates when the flow has subsided and reached 1/e (~37%) of the peak. The horizontal separation of the orange and blue curves indicates the duration between flood arrival and peak flow for a given location. The gray (dashed) boxes indicate the intervals in Fig. 2 during which the peaks in coherent seismic energy were detected. (C) Map view of region also with modeled flood arrival times. Times are Asia/Thimphu local time (UTC+6).

Table 1 Range of modeled flood arrival and peak flow times at populated villages.

View this table:

Geomorphic change caused by the GLOF

Satellite imagery reveals prominent changes in the Lunana region both before and after the GLOF. The most consistent change is a steady increase in the area of Lugge Tsho over the past 45 years, starting at approximately 0.42 km2 in January 1976. The lake area increased at a rate of 0.038 km2 year−1 as the glacier receded and melted, reached 1.1 km2 in September 1994, and then underwent a sudden decrease to 0.87 km2 during the GLOF in October 1994. Subsequently, the lake area continued increasing at a rate of 0.025 km2 year−1 and has since reached a size of 1.42 km2 as of September 2018 (fig. S4). We note that the lake expansion is primarily due to the receding glacier rather than rising lake level (see Discussion). The location of the 1994 breach is evident in both visible imagery and digital elevation model (DEM) on the lower left lateral moraine, where a cross section along the moraine crest shows a channel approximately 180 m wide and 40 m deep. Declassified satellite imagery from 1976 clearly shows that this location was a preexisting outlet for Lugge Tsho (fig. S5). Below the breach, changes in spectral reflectance are visible in post-GLOF Landsat imagery where the flood deposited debris and sediment along large regions of the valley floor. Upon first breaching the lake-fringing moraine, the flood waters flowed into another small lake approximately 500 m downstream. This seasonal lake was full at the time the GLOF occurred due to accumulated snowmelt and monsoonal precipitation from the prior summer months. The flood washed out the natural dam of this small lake basin; thus, now it no longer accumulates water as it did in years prior.

Around 10 km from the breach, a prominent set of glacier moraines is situated below Thanza (Fig. 4 and table S2). Here, we analyzed three boulders on the well-preserved outermost moraine ridge using 10Be surface exposure dating, following procedures in Schaefer et al. (33). Results yielded three consistent ages ranging from 4.4 to 4.7 thousand years (ka), indicating that the moraines were deposited during the mid-Holocene and left largely intact by the 1994 flood (fig. S5). After passing through the moraines, the flood wave spilled over the Tshojo plain, where it subsequently deposited a 2.2-km2 swath of sediment. In total, we estimate approximately 4.8 km2 of valley floor that was covered by sediment as a result of the GLOF. This occurred primarily along the upper 25 km of the river (fig. S6), as below this upper region, the channel is highly confined, and sediment deposition was minimal. In the decades following, the impacted vegetation has slowly recovered (fig. S7). In 1994, the mean enhanced vegetation index (EVI) of the Tshojo plain dropped to around 25% of the pre-GLOF value as a result of the flood. From 1994 to 2005, the EVI steadily increased as the vegetation began reclaiming the area, attaining around 75% of the pre-GLOF value in 2005 and remaining steady for several years afterward. Another increase occurred from 2011 to 2013, during which time an EVI approaching that of the pre-GLOF conditions was attained.

Fig. 4 Geomorphic map of the Lunana area.

The lower (red) moraines were deposited by the glacier more than 4000 years ago during the mid-Holocene, as indicated by the 10Be ages of three sampled boulders given with 1-s analytical error. The younger lake-fringing (purple) moraines are late-Holocene in age.


Arrival time of flood in Punakha

Accurate estimates of GLOF arrival times are vital for disaster preparedness. Existing model estimates for Punakha include 4.75 hours (19), 5.75 hours (20), and 7 hours or more (5). Given the onset of seismic energy at 1:45 a.m. at Lugge Tsho and eyewitness accounts of rising waters around 7:00 a.m. in the village, we estimate that the flood took approximately 5 hours to reach Punakha. Assuming that the delay between the initial breach and generation of detectable seismic energy was short (see Materials and Methods), our results suggest that some previously published simulations give reasonably accurate arrival times (within ±30 min) and that the shorter (4.75 to 5.75 hours, 17 km hour−1 average) times are more accurate than the 7-hour (12 km hour−1 average) estimate (5). These results show how a simple estimate of the GLOF start time based on the onset of seismic energy can be highly useful for validation of numerical models, by providing an estimate of the average velocity of the flood wave. However, a future breach may have different hydrograph characteristics depending on the trigger mechanism and nature of the moraine dam failure. The 1994 GLOF also cleared out a substantial amount of blocking debris, which would allow a future GLOF to travel more quickly down the valley. Further research toward constraining probable breach hydrograph characteristics and valley roughness parameters will be crucial for refinement of GLOF models. With the vast amount of existing seismic data [in databases such as the IRIS DMC (Incorporated Research Institutions for Seismology Data Management Center), for example] and increasing number of seismic networks, it is likely that other GLOF events have been recorded by seismic instruments but have yet to be investigated. A comprehensive search and deeper analysis of any available seismic records in GLOF-prone locations may reveal new insights into quantifying and modeling GLOF trigger events and flood waves in these regions.

Sediment deposition and vegetation recovery

The area covered by new sediment from the GLOF (4.8 km2) combined with existing sediment-covered areas (3.0 km2) amount to 7.8 km2. This approximately agrees with previous damage assessments, which estimated a total of 7.2 km2 of crops and pastures affected by the flood. The 2.2-km2 region of the most prominent GLOF sediment deposition is located in a scrub alpine vegetation zone, which receives around 1 m year−1 of precipitation primarily during the summer monsoon months (34). The flora in this area are composed of sedges, mosses, accessorial herbs, and some patches of woody vegetation (rhododendrons, junipers, and spireas) (35). The nonlinear vegetation recovery rates (fig. S7) are likely due to factors such as soil moisture, nutrient availability, competition between species, and seasonal precipitation in this high-elevation ecosystem. The observed recovery patterns may inform future studies regarding resilience of aromatic medicinal plants to changing climate, as these flora play key roles in the lives of local inhabitants (36).

Comparison with other GLOFs in the Himalayas

Compared to other GLOF occurrences in HMA during the past three decades, this 1994 event is a top contender for the largest volume of water released. Yet in terms of the percentage of lake volume released, it is on the lower end (fig. S8). Yamada et al. (18) surveyed the lake bathymetry and estimated the total volume of Lugge Tsho to be 58.3 million m3 in September 2002. Accounting for the lake growth between 1994 and 2002 (3 million m3) and volume of water released during the GLOF (17 to 25 million m3), we estimate that the lake volume was 73 to 80 million m3 in 1994. On the basis of this approximation, around 24 to 30% of the volume of Lugge Tsho was released. Because of the large size of the lake, substantial rate of outflow during the breach, and considerable vertical relief, this smaller-than-average percentage of draining resulted in a major destructive GLOF downstream. This is consistent with the flood simulations that suggest very low dampening and low deceleration of the flood wave peak due to high relief energy and the gorge character of the Pho River (20).

Regarding the seismic signal, we note that the GLOF event was exceptionally large, resulting in a high signal-to-noise ratio. Smaller GLOF events that produce less seismic energy would require seismometers to be located closer to the GLOF source to observe the weaker signal. This trade-off between geographical coverage of a seismic array versus the capability of detecting smaller events is an optimization problem, which future studies could address. There is no detectable coherent seismic energy generated from the river valley on the day before the GLOF event, suggesting that background-level river transport processes are not strong enough to be detected at these distances (fig. S12). A search for smaller GLOF events, which have occurred in other locations, along with any associated seismic signals (observed by nearby stations) may help further constrain signal strength versus distance from the source. Such observations could also be used to constrain numerical models that simulate seismic energy generated by various flood magnitudes.

Current state of high-risk lakes in Lunana

We find that Lugge Tsho surpassed the pre-GLOF size of 1.1 km2 in 2005 and reached an area of 1.4 km2 in 2018 due to retreat of the glacier terminus (fig. S4). If a future event were to cause the same 16 to 23 m of lake surface lowering as occurred during the 1994 GLOF, then this would translate to 22 to 32 million m3 of water released or approximately 30% more than the 1994 flood volume. While satellite imagery indicates that the lake boundaries (excluding the retreating glacier terminus) are nearly unchanged since the GLOF, the ongoing glacier retreat has also exposed unstable steep valley walls and lateral moraines above the lake (37, 38). A large mass movement into the water could result in a sudden increase in hydrostatic pressure and subsequent overtopping or structural failure of the Lugge Tsho moraine dam. The adjacent Thorthormi and Raphstreng lakes are also vulnerable to the same problem. Thorthormi sits topographically above Raphstreng by approximately 80 m, and a breach in its moraine could result in a cascading combined GLOF from both lakes (3). Efforts to artificially lower the level of Raphstreng by enlarging the outlet channel were undertaken during 2009–2012 in an attempt to reduce the risk, but this dangerous manual work by the local people was extremely difficult with uncertain effectiveness. On 20 June 2019, a minor breach occurred below Thorthormi lake, during which residents of Lunana were evacuated and no deaths or serious injuries occurred. While the increased flow was relatively minimal (water level increased by approximately 1 m during ~6:00 to 7:20 p.m.), it was large enough to wash away two bridges in Thanza and Tenchey. A field investigation after the event found that enhanced melting and basal sliding of the Thorthormi glacier had caused it to surge—this displaced water from Thorthormi lake, which overtopped the primary moraine, spilled into the subsidiary lakes, and breached the lower subsidiary lake that drained completely. Satellite imagery before and after the event confirmed the draining of the subsidiary lake, but the stability of the primary Thorthormi moraine dam is uncertain (39).

Potential for early warning systems

In Bhutan and elsewhere across the Himalayas, numerous glacial lakes pose immediate GLOF threats (2). Existing warning systems usually consist of automatic water level (AWL) stations installed in priority locations to monitor lake levels and river flows and transmit data in real time using Global System for Mobile Communications (GSM) or Iridium satellite technologies. If lake levels are detected to suddenly drop and/or stream levels rise, emergency warnings are issued (via mobile text messages), and a network of warning sirens is sounded. However, AWL sensors are known to be somewhat unreliable and susceptible to false alarms (40). Our results demonstrate the feasibility of seismic monitoring as another way to remotely detect GLOFs, which could potentially improve the next generation of early warning systems. The CCF methodology would need to be automated and tested more robustly to ensure reliable distinction between GLOF seismic signatures and other tectonic, meteorological, anthropogenic, or geomorphic sources (21), but with further refinement, a network of seismometers strategically deployed across a region could hypothetically monitor for sustained signals originating from probable GLOF source locations. To detect smaller GLOF events, a priori information should be used to focus detection efforts on known hot spots. Implementation of probability density functions in coherence map calculations may help remedy artifacts due to asymmetric distribution of seismic stations. During the CCF peak at the onset of the GLOF (around 1:45 a.m.), we find that only a few minutes of seismic data are sufficient to detect the anomalous high-frequency signal originating from the upper Lunana valley (fig. S9). Unlike AWL sensors, seismometers can be installed in safer and more accessible sites with the capability to monitor large regions across multiple valleys, although further research is needed to determine the optimal trade-off between array density and event detection capabilities. The use of seismometers and AWL sensors jointly could substantially improve existing early warning systems, with cross-validation of the two independent detection methods helping to minimize occurrence of false alarms and maximize warning time.


The risk of larger and potentially more destructive floods is rapidly increasing across the Himalayas, because of growing glacial lakes and ongoing construction of hydroelectric dams and other infrastructure in vulnerable regions (2). We have demonstrated that a GLOF can be detected remotely from seismometers located many kilometers away from the source, which has the potential to improve efficiency and maximize warning time. Our robust spatiotemporal constraints indicate at least 5 hours between source and main impact areas in Punakha, which affords sufficient warning time for many downstream inhabitants if a GLOF is detected in its earliest phase. We find that the 1994 flood event eroded, transported, and deposited large quantities of sediment in the Lunana valley but left the local moraine ridges largely intact near the source region. We also estimate a post-GLOF recovery of two to three decades for the affected alpine scrub vegetation, which may help to quantify resilience of the local ecosystem in future studies. Given the current situation and ongoing GLOF risks in the Himalayas, future research could focus on (i) deeper analysis and characterization of tectonic, meteorological, anthropogenic, and geomorphic seismic signatures to ensure clear distinction of natural hazard signals and prevention of false alarms; (ii) continued development of efficient algorithms for automated real-time processing of environmental seismic data; and (iii) deployment of optimized seismic arrays in vulnerable regions to detect GLOF events and provide efficient early warning systems.


Seismic data

Recent proliferation of high-quality broadband seismic data in addition to developments in the analysis of the ambient seismic wavefield and other seismic signals have forged new avenues in studying characteristics of seismic energy generated by environmental processes (10, 2124, 27). For example, time-frequency analyses of passive broadband seismic data have been used to quantify increases in high-frequency energy associated with high-flow regimes in rivers and cross-correlation between multiple stations used to isolate coherent seismic phases and provide estimates of their origin (28, 41). Here, we use seismic data from the INDEPTH II experiment as a tool to investigate the 1994 GLOF in Bhutan, which coincidentally occurred while this temporary seismic network was actively recording. INDEPTH II was a collaborative geoscience project between the Chinese Academy of Geological Sciences and investigators from U.S., German, and Canadian Geoscience institutions to investigate the deep structure and mechanics of the Himalaya-Tibet region (25). In 1994, the second phase of the project acquired passive seismic data in southern Tibet and the Himalayas, continuously recording three-component broadband and short-period 24-bit data along a ~350-km linear array at a sample rate of 50 Hz (42). For our analysis, we used data from a total of five broadband and short-period INDEPTH stations ranging from 75 to 135 km in distance (99 km average) to the northwest of the GLOF source area (Lugge Tsho). We downloaded the data from the IRIS DMC using a window for the approximate time of the GLOF occurrence (7 October 1994, from 00:00 to 12:00 hours, Asia/Thimphu time zone) for stations BB18, BB20, BB23, SP25, and SP27 (network code XR). The corresponding seismic traces were detrended, and instrument responses were removed to obtain units of velocity (m s−1) using the opensource Python framework ObsPy (43) and band-pass filtered between 1 and 5 Hz. This frequency range corresponded to the coherent high-frequency signal observed across all five stations (fig. S2) during the GLOF duration and also excluded lower frequency bands associated with noise sources such as ocean-generated microseisms (44, 45) as well as higher-frequency anthropogenic noise. Previous studies also observed a similar increase in seismic energy in these same frequency bands originating from turbulence and sediment transport by rivers and flood events (10, 22, 23, 27). For the purposes of this study, we assume that the beginning of detectable seismic energy marked the initiation of the GLOF event. We note that the actual outflow may have begun slightly earlier, but the seismic energy was below the threshold of detection at first due to a gradual increase in outflow through the moraine breach. This remains difficult to constrain as the exact shape of the breach hydrograph is unknown; thus, our evaluation of the time between the breach and downstream arrival of the flood wave is a minimum estimate.

Time-frequency analysis and CCFs

To explore the spectral characteristics of the event and quantify the temporal variation of the seismic signal generated by the GLOF, we estimated the power spectral density of the time series for each station using Welch’s averaging method. We first divided each seismic trace into 2-min segments each with 50% overlap. Then, for each 2-min segment, we used Welch’s method to average modified periodograms computed using 10-s windows, also overlapping by 50%. To approximate source locations of the GLOF energy, we followed an approach similar to those outlined by previous studies for locating coherent seismic noise sources (28, 41, 46). We first applied a 1-bit normalization to reduce the influence of punctual sources of seismic energy, such as earthquakes, anthropogenic noise, or instrument issues. This simply means keeping the sign of the time series (−1 if less than 0 and +1 if greater than 0) and discarding the magnitude (46). We then calculated the normalized cross-correlation of 20-min segments (overlapping by 50%) in the time series and computed their envelopes (hereafter referred to as the CCFs) for every station pair along the seismic array for time lags ranging from −40 to +40 s for a series of frequencies ranging from 1 to 5 Hz, using a window size of 0.5 Hz. Time information from each CCF envelope was migrated to positions in space as follows: We defined a regular grid of potential source locations in the region, and for each station pair, we calculated the theoretical time delays between the two stations for every grid point. The CCF amplitude at each corresponding lag time was then mapped to positions in space. The resulting coherence map Aij for stations i and j is given byAij(x,y)=CCFij (didjv)(1)where d is the distance between the hypothetical source and corresponding station and v is the assumed seismic velocity (see below). Thus, each CCF delay time (along the vertical axes in fig. S3) maps to a hyperbola, where the amplitude of the hyperbola is simply the CCF amplitude. This was repeated for each 20-min segment and station pair. The resulting maps for hand-selected station pairs and frequencies (those with distinct correlation peaks; see fig. S3 and table S1) were summed together to form a final coherence map for a given interval, where high coherence values indicate the most probable source locations (Fig. 2). Similar results are also obtained if all station pairs are included in the stack (fig. S11). To determine an appropriate velocity during migration, we calculated coherence maps for a range of velocities between 1 and 5 km s−1. We found that a velocity of 3.0 km s−1 resulted in the highest coherence (fig. S3) for these frequencies, similar to previous studies (28), and likely indicates short-period Rayleigh wave energy. We found that two distinct peaks in coherence occurred at approximately 2:00 and 6:30 a.m. and chose to focus on these peaks for further analysis. We defined time windows of 90-min duration spanning each respective peak, which we found to be a sufficient length for capturing the intervals of strongest GLOF signal relative to background levels. We also note that quasi-linear placement of the seismometers parallel to the valley northwest of the GLOF means that the source location was well constrained along the valley but poorly constrained perpendicular to the valley, leading to blurring in that direction. This artifact may be remedied by invoking a probability density function centered on the river channel to better localize the signal, although the GLOF event studied here was large enough such that this was not necessary. While this and other sources of error such as lateral velocity heterogeneities and varying surface topography caused some blurring of the coherence maps, we found this basic methodology precise enough to clearly track the start and a subsequent down-valley shift in the location of seismic energy during the GLOF.

Flood model

We implemented the U.S. Army Corps of Engineers HEC-RAS software to perform a 2D unsteady flow simulation using the diffusion-wave equation (47). This equation is a simplified version of a full dynamic wave model (neglects inertial force and advective accelerations) and has been found to be a satisfactory approximation in many situations (48). Simulations were run at a nominal mesh resolution of 30 m using a time step of 1 s, solved using an implicit finite-volume approach (47). We used the 30-m Advanced Land Observing Satellite (ALOS) DEM as terrain input, preprocessed using standard carve and fill operations to remove any spurious elevation artifacts that may cause unrealistic damming and pooling in localized sections where the river channel may be narrower than the DEM resolution (12, 49). To approximate the normal (pre-GLOF) river flow conditions, we specified inflows for 15 major tributaries between Lugge Tsho and Wangdue station and allowed the model to come to equilibrium. The (relative) contributions of each tributary were estimated by performing a flow accumulation analysis for the region based on upstream watershed areas (49). We then multiplied all the relative inflow values with a single scale factor to estimate (absolute) contributions, so that the river flow at Wangdue station (located downstream of all tributaries) matched the observed pre-GLOF conditions of approximately 290 m3 s−1 (fig. S1). In HEC-RAS, we expressed the tributaries as inflow boundary conditions and allowed the model to run for 24 hours to establish initial conditions before simulating the flood. We then performed multiple model runs using a range of Manning roughness coefficient values (n) and various breach hydrograph shapes (fig. S1). We tested values of n spanning from 0.05 to 0.07 (in increments of 0.01), which is the typical range for mountain streams with cobbles and large boulders (47). For the breach hydrographs, we used a simple triangular approximation scaled to have ramp-up times (tru) ranging from 15 to 120 min (in increments of 15 min). On the basis of previously published differential GPS survey of the lowering of the lake level (18) and the Wangdue station GLOF hydrograph (19), we assumed 17 to 25 million m3 as a probable range for the total volume of water released during the GLOF and constrained the simulated breach hydrographs to 25 million m3. As a conservative threshold, we discarded any model runs, which did not agree with all independent constraints within ±30 min, and report the corresponding range of input parameters and model outputs of the remaining ones. We note that the breach may have initiated before the seismic signal was detectable (see “Seismic data” section); thus, we also include simulated breach times of 1:15, 1:30, and 1:45 a.m. in our analysis. Out of a total of 72 model runs (3 values of n, 8 values of tru, and 3 breach times), 14 runs produced output that satisfied the specified ±30-min threshold. Last, we quantified the model sensitivity of arrival time estimates in Punakha (fig. S1), and found that a more gradual release of water and greater channel roughness both resulted in a slower-moving flood wave. Increasing tru from 45 to 60 min delayed arrival in Punakha by 20 to 25 min (fig. S1E), and increasing n from 0.05 to 0.06 delayed arrival time in Punakha by 25 to 35 min (fig. S1F). In general, we found that the HEC-RAS model performed well in satisfying the independent observations within the specified range of model parameters. However, we stress that the overall sensitivity to breach hydrograph characteristics requires caution if external constraints are not available.

Satellite imagery

For analyzing GLOF-induced changes in land cover, we used the U.S. Geological Survey (USGS) Landsat 5 through 8 Thematic Mapper (TM) Collection 1 Tier 1 calibrated top-of-atmosphere (TOA) reflectance product in Google Earth Engine (50). We quantified Lugge Tsho area changes over time by manual delineation of the lake boundaries and glacier front between 1976 and 2018 using declassified spy satellite imagery (KH-9 Hexagon) and Landsat. The Hexagon images were downloaded from the USGS Earth Explorer website ( after digital scanning of the images by the USGS and have a ground resolution of approximately 15 m. To measure the extent of sediments deposited by the GLOF, we selected two Landsat 5 scenes acquired in August 1994 (pre-GLOF) and September 1995 (post-GLOF) and used supervised classification with manually defined training samples of vegetation, sediment, water, ice, and clouds from the pre-GLOF scene. We classified both images using the maximum likelihood algorithm and retained the pixels classified as sediment along the valley bottom before and after the GLOF. To analyze the post-GLOF vegetation recovery trend, we focused on the largest swath of sediment deposition over the Tshojo plain. For this region, we computed the EVI for all Landsat 5 and 7 scenes, excluding those acquired during monsoon season months (May to October). We evaluated the topography of the Lunana region using the High Mountain Asia 8-m DEM data (version 1) distributed by the National Snow and Ice Data Center (NSIDC) and obtained the valley profile from the 30-m ALOS global digital surface model dataset distributed by the Japan Aerospace Exploration Agency (JAXA).

Cosmogenic 10Be surface exposure dating

We applied cosmogenic 10Be dating (51, 52) to three boulders sampled in 2014 from the prominent and well-preserved lateral-terminal moraines near Thanza, about 10 km downstream from the GLOF breach location (Fig. 4). Geochemical processing was performed at the Cosmogenic Nuclide Laboratory at Lamont-Doherty Earth Observatory (LDEO) following standard protocols given in Schaefer et al. (33), and 10Be/9Be measurements were completed at the Center for Accelerator Mass Spectrometry at the Lawrence Livermore National Laboratory. The background correction for these measurements was below 1%. We used version 3 of the online cosmogenic nuclide calculator (53) with the default production rate and time-dependent Stone/Lal scaling scheme for exposure age calculations (52, 54). Geographic and analytical data are given in table S2, and the geomorphic map with the 10Be ages is shown in Fig. 4. This glacier chronology will be discussed in detail in forthcoming papers.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank V. Tsai and two anonymous reviewers for helping improve this manuscript. We would also like to acknowledge the National Center for Hydrology and Meteorology, Royal Government of Bhutan for providing field logistics and support for this study. Funding: J.M.M. and J.M.S. acknowledge support by a NASA Earth and Space Science Fellowship (no. NNX16AO59H). S.R. and J.M.S. acknowledge support by the NSF Geography and Spatial Sciences award 17-566. J.M.S. also acknowledges support by Global Change award EAR 10-574. S.R. acknowledges support by NASA HMA no. NNX16AQ61G and NSF no. 1853881. J.M.S. acknowledges funding by the G. Unger Vetlesen Foundation. Author contributions: J.M.M., J.M.S., and S.R. framed the research questions and designed the study. J.M.M. and J.B.R. performed the seismic analyses. J.M.M. performed the remote sensing and flood modeling analyses with guidance from J.M.S. and S.R. N.W. provided logistical support, hydrological data, and a compilation of news articles and technical reports relating to the 1994 GLOF. A.E.P., N.Y., and J.M.S. provided the 10Be surface exposure ages. J.M.M. led the writing of the manuscript with input and contributions from all coauthors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Seismic data used in the paper are available for download from the IRIS DMC website ( This is LDEO publication #8438.
View Abstract

Stay Connected to Science Advances

Navigate This Article