# Hazard from Himalayan glacier lake outburst floods

See allHide authors and affiliations

Edited by Andrea Rinaldo, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, and approved November 27, 2019 (received for review August 27, 2019)

## Significance

Glacier lake outburst floods (GLOFs) have become emblematic of a changing mountain cryosphere. The Himalayas suffered the highest losses from these sudden pulses of meltwater but lack a quantitative appraisal of GLOF hazard. We express the hazard from Himalayan glacier lakes by the peak discharge for a given return period. The 100-y GLOF has a mean discharge of ∼15,600 m^{3}⋅s^{−1}, comparable to monsoonal river discharges hundreds of kilometers downstream. The Eastern Himalayas are a hotspot of GLOF hazard that is 3 times higher than in any other Himalayan region. The size of growing glacier lakes and the frequency of lake outbursts determine GLOF hazard, which needs to be acknowledged better in flood hazard studies.

## Abstract

Sustained glacier melt in the Himalayas has gradually spawned more than 5,000 glacier lakes that are dammed by potentially unstable moraines. When such dams break, glacier lake outburst floods (GLOFs) can cause catastrophic societal and geomorphic impacts. We present a robust probabilistic estimate of average GLOFs return periods in the Himalayan region, drawing on 5.4 billion simulations. We find that the 100-y outburst flood has an average volume of 33.5^{+3.7}/_{−3.7} × 10^{6} m^{3} (posterior mean and 95% highest density interval [HDI]) with a peak discharge of 15,600^{+2,000}/_{−1,800} m^{3}⋅s^{−1}. Our estimated GLOF hazard is tied to the rate of historic lake outbursts and the number of present lakes, which both are highest in the Eastern Himalayas. There, the estimated 100-y GLOF discharge (∼14,500 m^{3}⋅s^{−1}) is more than 3 times that of the adjacent Nyainqentanglha Mountains, and at least an order of magnitude higher than in the Hindu Kush, Karakoram, and Western Himalayas. The GLOF hazard may increase in these regions that currently have large glaciers, but few lakes, if future projected ice loss generates more unstable moraine-dammed lakes than we recognize today. Flood peaks from GLOFs mostly attenuate within Himalayan headwaters, but can rival monsoon-fed discharges in major rivers hundreds to thousands of kilometers downstream. Projections of future hazard from meteorological floods need to account for the extreme runoffs during lake outbursts, given the increasing trends in population, infrastructure, and hydropower projects in Himalayan headwaters.

Monsoonal floods are among the most destructive natural hazards in the greater Himalayan region and the adjacent mountain ranges of the Hindu Kush, Karakoram, Nyainqentanglha, and Hengduan Shan (Fig. 1). Regional projections for the lower Indus, Ganges, and Brahmaputra rivers hold that flood frequencies will rise noticeably in the 21st century (1, 2), putting the livelihoods of 220 million people at risk (3). In Himalayan headwaters, such prognoses have disregarded episodic, but potentially destructive, floods from the sudden emptying of moraine-dammed lakes. Such glacier lake outburst floods (GLOFs) occur largely independently of hydrometeorological floods, but can surpass their peak discharges by orders of magnitude in the upper river reaches (4⇓–6). Glacier lakes dammed by abandoned moraines are susceptible to outburst, triggered by ice or debris falls, strong earthquake shaking, internal piping, or overtopping waves that exceed the shear resistance of the dam (7⇓–9). These triggers mostly happen unrecorded in remote terrain, eroding the impounding barriers within minutes to hours, and releasing sediment-laden floods that may travel >100 km downstream (7). With little to no warning, communities and infrastructure downstream have often been caught unprepared, suffering loss of human lives and livestock, and damage to roads, buildings, and hydropower facilities (10⇓–12). An objective and reproducible hazard assessment of such dam-break floods is key to human safety and sustainable development, and is repeatedly emphasized in research and media coverage of atmospheric warming, dwindling glaciers, and growing meltwater lakes (3, 13, 14).

GLOFs have gained growing attention in the Himalayas (15, 16), where these disasters have had the highest death toll worldwide (17). While the distribution and dynamics of moraine-dammed lakes have been mapped extensively in recent years (12, 18⇓⇓–21), objectively appraising the current Himalayan GLOF hazard has remained challenging. The high-alpine conditions limit detailed fieldwork, leading researchers to extract proxies of hazard from increasingly detailed digital topographic data and satellite imagery. These data allow for readily measuring or estimating the geometry of ice and moraine dams, the possibility of avalanches or landslides entering a lake, or the water volumes released by outbursts (9, 18, 20). Ranking these diagnostics in GLOF hazard appraisals has mostly relied on expert judgment (18, 22), because triggers and conditioning factors are largely unknown for most historic GLOFs (23).

Guided by studies on earthquakes, landslides, wildfires, and floods, we use the magnitude and frequency of GLOFs in a given area and period as an objective metric of hazard. GLOFs have occurred at an average rate of ∼1.3 y^{−1} in the Himalayas in the past 3 decades (23), although of differing size, which is conventionally described by the released flood volume *V*_{0} or the associated peak discharge *Q*_{p}. These parameters are difficult to measure during dam failure but can be estimated eventually from the surface areas of glacier lakes. We thus scaled the manually mapped areas of all 5,184 Himalayan moraine-dammed lakes (>0.01 km^{2}) to volumes using a Bayesian robust regression of maximum lake depth versus area from 24 bathymetrically surveyed lakes (*Materials and Methods* and *SI Appendix*, Fig. S2 and Table S1). Our approach assumes that all lakes are equally susceptible to outburst and that any incision depth is possible, which is consistent with data from Himalayan dam breaks (*SI Appendix*, Fig. S3). We used a Bayesian variant of a physical dam-break model (24, 25) to predict peak discharge *Q*_{p} from the product of *V*_{0} and the breach erosion rate *k* (*Materials and Methods* and *SI Appendix*, Figs. S4 and S5). In summary, we obtained 4.6 × 10^{8} and 4.9 × 10^{9} scenarios of *V*_{0} and *Q*_{p}, respectively, for all moraine-dammed lakes in the Himalayas. These numbers ensure that we have sufficiently explored the physically plausible space for fitting an extreme-value distribution to these key flood diagnostics. We further stacked our simulations, accounting for contemporary mean annual GLOF rates between 0.03 and 0.71 y^{−1} in 7 subregions of the Himalayas, and estimated GLOF return periods in terms of the 100-y flood volume *V*_{100} and discharge *Q*_{p100} (*Materials and Methods*) (26, 27).

## Results

The predicted flood volumes and peak discharges span more than 7 orders of magnitude (Fig. 2). Based on a mean posterior rate of 1.26 GLOFs y^{−1} over the past 3 decades (23), we estimate a contemporary *V*_{100} of 33.5^{+3.7}/_{−3.7} × 10^{6} m^{3}⋅s^{−1} and *Q*_{p100} of 15,600^{+2,000}/_{−1,800} m^{3}⋅s^{−1} in the entire study area (Fig. 3). Regionally differing GLOF rates cause variation in *V*_{100} and *Q*_{p100}. For example, GLOF hazard in the Eastern Himalayas is 350% (*V*_{100}) and 338% (*Q*_{p100}) higher than in the next highest region (Fig. 3 and *SI Appendix*, Fig. S6).

Correcting for regionally varying GLOF rates, we can explore how—for a fixed rate of one GLOF y^{−1}—hazard changes as function of the lake size distribution (Fig. 4). We find that GLOF hazard in the Eastern Himalayas, which currently have the largest lakes and have had the highest GLOF activity in the past 3 decades (23), is on average more than 87% (*V*_{100}) and 57% (*Q*_{p100}) higher than the estimate from the entire Himalayas. In contrast, the Western Himalayas had no known outburst from moraine-dammed lakes in the past 30 y, while the rate-adjusted *V*_{100} and *Q*_{p100} there are still about 76% and 58% below that of the overall study area. In summary, these figures robustly support a contemporary hotspot of GLOF hazard in the southern Himalayas, mainly due to the larger abundance of glacier lakes and GLOFs in that part of the mountain belt.

## Discussion

We offer a consistent and reproducible estimate of present GLOF hazard in the Himalayas. Our Bayesian estimates explore the parameter space of plausible flood volumes and associated peak discharges with roughly a million outburst scenarios for any given lake. Our approach expands previous hazard appraisals by explicitly accounting for regionally varying GLOF rates. The estimated GLOF return periods are consistent with the frequency of reported flood volumes and discharges since 1935 (Fig. 3 and *SI Appendix*, Fig. S6 and Table S2). The 3 largest reported flood volumes of 71.6, 19.5, and 17.2 × 10^{6} m^{3} from the lakes Sangwang Tsho (1954), Sabai Tsho (1998), and Lugge Tsho (1994), respectively, had return periods of 237, 56, and 49 y, respectively, according to our predictions. The highest reported peak discharge (15,920 m^{3}⋅s^{−1}) drained from Lake Zhangzangbo in 1981; we estimate this as a 100-y event in the greater Himalayan region. The lakes Sangwang Tsho and Sabai Tsho had similar high peak discharges of ∼10,000 m^{3}⋅s^{−1}, corresponding to a return period of 60 y. Yet reported values of *V*_{0} and *Q*_{p} rely chiefly on point estimates from empirical rating curves, eyewitness accounts, flood marks, and transported clasts in river channels, or measurements several kilometers downstream of the failing dams (28), thus compromising detailed, site-specific validation of our regional predictions. These uncertainties cause some of the scatter in the data that we based our models on (*SI Appendix*, Figs. S4 and S5). However, our approach explicitly propagates these uncertainties and robustly estimates *V*_{0} and *Q*_{p} from the distribution of meltwater lake areas in the entire Himalayan region. The size of glacier lakes will remain the key determinant of GLOF magnitude, regardless of whether we use alternative discharge rating curves (*SI Appendix*, Fig. S7), more elaborate numerical dam breach simulations (29), or any other metric of flood potential (20).

Our study objectively enriches previous qualitative appraisals, which contrasted in the degree and number of “potentially dangerous” Himalayan lakes (22). The tails of the frequency−size distributions (Fig. 2) contain valuable information for practitioners regarding the physically possible margins of *V*_{0} and *Q*_{p} in the Himalayas. We note that half of the total Himalayan lake volume (20.1^{+2.4}/_{−2.1} km^{3}) is in the largest 1% of glacier lakes (*SI Appendix*, Fig. S8). For example, Sangwang Tsho—the largest Himalayan moraine-dammed lake today—accounts for 5% of the total glacier lake volume (∼1 km^{3}) and could release a flood with a return period of 5,000^{+1,400}/_{−1,000} y, should this volume drain entirely (*SI Appendix*, Fig. S9). Such flood magnitudes and associated return periods are consistent with reports on Holocene flood volumes up to 3 orders of magnitude larger in the eastern Himalayas (30).

Himalayan glacier lakes are currently limited in their storage capacity to produce floods of these dimensions, particularly in the Hindu Kush, Karakoram, and Western Himalayas, where lakes are smallest today (Figs. 1 and 2 and *SI Appendix*, Fig. S8). Why glacier lakes occur in clusters and have become largest in the Central and Eastern Himalayas since 1990 (+7.7 to +23% in area) remains debated (19). Glaciers have had negative mass balances in all regions except the Karakoram since the 1970s at least (31), while the highest ice losses in the 21st century were in the Western Himalaya and the Nyainqentanglha (32). The Central and Eastern Himalayas have a high proportion of lakes in direct contact with debris-covered parent glaciers, possibly accelerating glacier mass loss and lake growth (33). The increasing trend of air temperature is consistent over the entire Himalayas, whereas changes in precipitation were mostly insignificant (34), so that the clustering of lakes today may be tied to local, less well-known, controls such as glacier bed topography or valley geometry (35).

Yet the distribution of glacier lakes will likely change in the future, given that a global temperature rise of 1.5 °C could melt half of the Himalayan glacier mass by the end of the 21st century (36). Models of subglacial topography suggest that a complete melting of all glaciers might provide accommodation space for another 16,000 meltwater lakes (>10^{4} m^{2}) with a maximum total volume of 120 km^{3} (37), hence possibly increasing the current lake volume 6-fold. According to these projections, most of these futures lakes could form in the Karakoram and Western Himalayas, where contemporary glacier cover is highest (13). Yet the errors of the inferred bathymetries can involve tens of meters, and many lakes might fill bedrock depressions that are less prone to catastrophic outbursts. Future flood volumes may increase nevertheless, especially in regions with large ice volumes and few lakes, even if only a fraction of the total meltwater will be trapped by moraines. Regional GLOF hazard could thus change commensurately, should annual GLOF frequencies rise with a growing lake abundance in the future (16). Whether lake outbursts will become more frequent with atmospheric warming remains an open question, however. The average rate of GLOFs in the greater Himalayan region has remained unchanged in the past 3 decades (23), showing that rapid growth of glacier lakes alone is an unsuitable predictor for GLOF activity, let alone GLOF hazard (38, 39).

Both the magnitude and frequency of outburst floods are straightforward to change in our model, motivating a more dynamic assessment with regular updates of outburst frequency and lake size distribution. We also encourage future work to learn more about the triggers and drivers of historic GLOFs. With more of this information becoming available, we can improve our Bayesian framework by assigning more informed (and appropriately weighted) prior probabilities of outburst to each glacier lake. In this spirit, we offer a practical tool that is extensible and compatible with flood routing models, design codes, and hazard mitigation. Our appraisal of GLOF hazard complements projections of meteorological flood hazards in a warming climate for the sparsely instrumented Himalayan drainage network. Atmospheric warming is projected to increase mean daily and annual discharges in the Indus, Ganges, and Brahmaputra rivers only by a few percent in this century (40, 41), although the return periods for a given flood stage might drop by 50 to 90% in these rivers in the 21st century. Extreme runoff in headwaters (1), and GLOFs in particular (12), has eluded such projections. We show that flood hazard due to glacier melt requires urgent attention in light of the rapidly growing population, infrastructure, and hydropower projects in the Himalayas (3, 12). Our appraisal adds to recent calls toward identifying regions that are most prone to GLOF impacts, and complements integrative concepts for climate risk adaptation on local, national, or regional levels (42).

## Materials and Methods

### Inventory of Moraine-Dammed Glacier Lakes.

We use *glacier lake* and *meltwater lake* synonymously for water bodies fed by glacier meltwater, snow, and ice as defined in ref. 21. This database is the only consistent, highly accurate glacier lake inventory for our study region, containing a total of 25,614 manually mapped lakes from Landsat images between 2003 and 2007. We extracted all (*n* = 5,184) moraine-dammed glacier lakes >0.01 km^{2} for our 7 subregions, which follow the glacier regions in the Randolph Glacier Inventory. We slightly modified the outlines of the Central and Eastern Himalayas to include glaciers and associated glacier lakes facing toward the Tibetan Plateau.

### Estimating Depths, Volumes, and Flood Volumes of Glacier Lakes.

Global samples from glacier lakes have suggested that glacier depths can vary by one order of magnitude for a given lake area (43). To estimate glacier lake volumes in our study region, we compiled data from bathymetric surveys of 24 Himalayan glacier lakes with reported lake areas a and maximum depths d (*SI Appendix*, Table S1). We used a Bayesian robust linear regression with a normally distributed target variable (lake depth) to account for possible effects of the small sample size and outliers in our data (*SI Appendix*, Fig. S1*A*): *SI Appendix*, Fig. S2).

We assumed that the manually mapped glacier lake areas are circular, so that the radius r of any lake is described as a function of the lake area a:

### Estimating Peak Discharge.

Walder and O’Connor (24) proposed a physically motivated model that predicts peak discharge *Q*_{p} during natural dam failure, based on the breach rate, breach depth, and volume of water released. They compiled data from 63 observed natural dam breaks in diverse settings, and found that the largest values of *Q*_{p} arise from dam breaches that fully develop before the lake level drops substantially. Such large impoundments involve either large lake volumes relative to the breach depth or very rapid dam failure. In contrast, mountain valleys mostly store smaller volumes even behind high dams, so that peak discharge often occurs prior to complete dam incision (24). These 2 end members make up a constant response of dimensionless peak discharge *Q*_{p}* when plotted against the dimensionless product *η* of lake volume and breach rate (24) (*SI Appendix*, Fig. S5). From equations for critical flow, Walder and O’Connor inferred a model with the 2 key components dimensionless peak discharge

and

where *SI Appendix*, Fig. S4). The capped values of *Q*_{p}* express the model’s idea that breach conditions primarily control peak discharge, so that outflow deceleration or downstream ponding limit any further increase in discharge for a given breach geometry.

The empirical data support a piecewise regression model of the form *SI Appendix*, Fig. S5). To better constrain this transition and to explicitly treat all other uncertainties, we learned a Bayesian piecewise robust model from these empirical data (*SI Appendix*, Fig. S1*B*). We assumed a Student’s t distributed noise to make our regression robust against outliers, and specified the likelihood of observing the data

with

where we encoded our prior assumptions about the regression parameters as *SI Appendix*, Fig. S5). We further specified a half-Cauchy prior on the scale parameter *t*-distributed piecewise model at a low value of ν = 10 to impose taller tails on the regression coefficients, hence making it more robust against outliers. We used a numerical sampling scheme implemented in the Stan programming language that relies on MCMC to approximate the marginal posterior distributions of

We integrated out these parameters to simulate draws from the predictive posterior distribution for unobserved inputs, so that we can predict *V*_{0} and 100 physically plausible values of the breach rate k based on a log-normal fit to reported breach rates (24) (*SI Appendix*, Fig. S4). We allowed for twice the range of documented breach rates to account for the effects of failure mechanisms that may have remained unobserved. Multiplying 10^{5} samples of *V*_{0} with 100 samples of k finally resulted in a total of 10^{7} scenarios of

### Estimating GLOF Return Periods.

We used the posterior estimates of *SI Appendix*, Fig. S1*C*). The Generalized Pareto (GP) distribution has been widely used to characterize the distribution of quantities recorded above a high threshold u (also known as the *peak-over-threshold* or *partial duration series* method), and has the probability density

where u is a location parameter, κ ≥ 0 is a scale parameter, and γ is a shape parameter for

We further assumed that GLOFs occur in time according to a Poisson process with rate

where the terms on the right-hand side containing λ denote the time component, and

We obtained the mean annual GLOF rate λ in a given region from remote sensing-based analyses covering the past 3 decades (23, 45). These studies detected 38 GLOFs in 30 y (1988–2017), and demonstrated that the mean annual rate of GLOFs remained unchanged over this period, both in the entire Himalayas and in its 7 subregions. We thus assumed constant posterior annual GLOF rates (brown numbers in Fig. 3 and *SI Appendix*, Fig. S6) in any region, motivating a Poisson-distributed occurrence of GLOFs in time. To estimate GLOF return periods, we created a mixture distribution

We used the mean regional GLOF rate λ (i.e., adjusted by the historic GLOF count from each region in 30 y) to simulate 10,000 y of data from the mixed posterior predictive distribution of *SI Appendix*, Fig. S1*C*). We fitted a GP model to these time series with an MCMC sampler implemented in the package extRemes (46) in the statistical programming language R. We found predictions to be robust for upper quartile thresholds and used the 80th percentiles of the data as thresholds eventually, adjusting regional estimates by the appropriate regional values of λ. We repeated this sampling scheme 1,000 times for all Himalayan regions and report the mean from all average return periods and their 95% HDI. By design of our extreme-value model, changes in annual GLOF rates directly change return periods (*SI Appendix*, Figs. S10 and S11): For example, a doubling of GLOF frequency in a given period shortens the return period of this flood volumes or discharge level by one-half.

### Data and Code Availability.

All input data are available at Zenodo (http://doi.org/10.5281/zenodo.3523213), and model codes are available at GitHub (https://github.com/geveh/GLOFhazard).

## Acknowledgments

This study was funded by Deutsche Forschungsgemeinschaft within the graduate research training group NatRiskChange (Grant GRK 2043/1) at the University of Potsdam (https://www.uni-potsdam.de/natriskchange). We performed our spatial analysis using the System for Automated Geoscientific Analyses (SAGA V. 2.2.0 and 7.3.0) and ArcGIS (V. 10.5), and performed the statistical analyses with R (http://www.r-project.org/) and the probabilistic programming language Stan (https://mc-stan.org/). We thank Kristin Vogel for comments on an early version of this manuscript.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: georg.veh{at}uni-potsdam.de.

Author contributions: G.V. and O.K. designed research; G.V. and O.K. performed research; G.V. and O.K. analyzed data; and G.V., O.K., and A.W. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

Data deposition: All input data are available at Zenodo (http://doi.org/10.5281/zenodo.3523213), and model codes are available at GitHub (https://github.com/geveh/GLOFhazard).

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1914898117/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- R. R.Wijngaard et al

- ↵
- ↵
- P. Wester,
- A. Mishra,
- A. Mukherji,
- A. B. Shrestha

- ↵
- K. L. Cook,
- C. Andermann,
- F. Gimbert,
- B. R. Adhikari,
- N. Hovius

- ↵
- ↵
- R. Osti,
- S. Egashira

- ↵
- S. D. Richardson,
- J. M. Reynolds

- ↵
- ↵
- A. Emmer,
- V. Vilímek

- ↵
- S. K. Allen,
- P. Rastner,
- M. Arora,
- C. Huggel,
- M. Stoffel

- ↵
- T. Watanbe,
- D. Rothacher

- ↵
- ↵
- T. Bolch et al

- ↵
- ↵
- Y. Nie et al

- ↵
- S. Harrison et al

- ↵
- J. L. Carrivick,
- F. S. Tweed

- ↵
- X. Wang et al

- ↵
- Y. Nie

- ↵
- K. Fujita et al

- ↵
- S. B. Maharjan et al

- ↵
- D. R. Rounce,
- D. C. McKinney,
- J. M. Lala,
- A. C. Byers,
- C. S. Watson

- ↵
- G. Veh,
- O. Korup,
- S. von Specht,
- S. Roessner,
- A. Walz

- ↵
- J. S. Walder,
- J. E. O’Connor

- ↵
- D. M. Burr,
- P. A. Carling,
- V. R. Baker

- J. E. O’Connor,
- R. A. Beebee

- ↵
- G. Veh,
- O. Korup,
- A. Walz

- ↵
- G. Veh,
- O. Korup,
- A. Walz

- ↵
- O. Korup,
- F. Tweed

- ↵
- M. J. Westoby et al

- ↵
- D. R. Montgomery et al

- ↵
- T. Bolch et al

- ↵
- F. Brun,
- E. Berthier,
- P. Wagnon,
- A. Kääb,
- D. Treichler

- ↵
- C. Song et al

- ↵
- J. M. Maurer,
- J. M. Schaefer,
- S. Rupper,
- A. Corley

- ↵
- O. King,
- A. Dehecq,
- D. Quincey,
- J. Carrivick

- ↵
- P. D. A. Kraaijenbrink,
- M. F. P. Bierkens,
- A. F. Lutz,
- W. W. Immerzeel

- ↵
- A. Linsbauer et al

- ↵
- ↵
- D. I. Benn et al

- ↵
- S. Nepal,
- A. B. Shrestha

- ↵
- ↵
- S. K. Allen et al

- ↵
- S. J. Cook,
- D. J. Quincey

- ↵
- M. Plummer

- ↵
- G. Veh,
- O. Korup,
- S. Roessner,
- A. Walz

- ↵
- E. Gilleland,
- R. W. Katz

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Environmental Sciences