# Less extreme and earlier outbursts of ice-dammed lakes since 1900

Georg Veh, Natalie Luetzow, Jenny Tamm, Lisa V. Luna, Romain Hugonnet, Kristin Vogel, Marten Geertsema, John J. Clague, Oliver Korup

<figure><img src="image_1.png" /></figure>

## To cite this version:

Georg Veh, Natalie Luetzow, Jenny Tamm, Lisa V. Luna, Romain Hugonnet, et al. Less extreme and earlier outbursts of ice-dammed lakes since 1900. Nature, 2023, 10.1038/s41586-022-05642-9. hal-04720433

HAL Id: hal-04720433

https://hal.science/hal-04720433v1

Submitted on 4 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Less extreme and earlier outbursts of ice-dammed lakes since 1900

https://doi.org/10.1038/s41586-022-05642-9

Received: 3 March 2022

Accepted: 7 December 2022

Published online: 15 February 2023

Open access

<figure><img src="image_3.png" /></figure>

Check for updates

Georg Veh<fnref n="1" />,²,³, Natalie Lützow<fnref n="1" />, Jenny Tamm<fnref n="1" />, Lisa V. Luna<fnref n="1" />,²,³, Romain Hugonnet⁴,⁵,⁶, Kristin Vogel⁷, Marten Geertsema⁸, John J. Clague⁹ & Oliver Korup<fnref n="1" />,²

Episodic failures of ice-dammed lakes have produced some of the largest floods in history, with disastrous consequences for communities in high mountains<fnref n="1" />⁻⁷. Yet, estimating changes in the activity of ice-dam failures through time remains controversial because of inconsistent regional flood databases. Here, by collating 1,569 ice-dam failures in six major mountain regions, we systematically assess trends in peak discharge, volume, annual timing and source elevation between 1900 and 2021. We show that extreme peak flows and volumes (10 per cent highest) have declined by about an order of magnitude over this period in five of the six regions, whereas median flood discharges have fallen less or have remained unchanged. Ice-dam floods worldwide today originate at higher elevations and happen about six weeks earlier in the year than in 1900. Individual ice-dammed lakes with repeated outbursts show similar negative trends in magnitude and earlier occurrence, although with only moderate correlation to glacier thinning⁸. We anticipate that ice dams will continue to fail in the near future, even as glaciers thin and recede. Yet widespread deglaciation, projected for nearly all regions by the end of the twenty-first century⁹, may bring most outburst activity to a halt.

One of the major hazards caused by ongoing glacier retreat are sudden floods resulting from the failure of unstable glacier dams⁴⁻⁷. These ice dams temporarily and periodically impound rainfall and meltwater that can be released abruptly following dam flotation, subglacial tunnel enlargement or mechanical collapse¹⁰. Such glacier lake outburst floods (GLOFs) initiate hazard cascades with numerous impacts downstream. Historic ice-dammed GLOFs have travelled hundreds of kilometres, entraining sediment along their runout paths; widening river channels by bank erosion; triggering landslides from undercut and destabilized hillslopes; and covering valley floors with metre-thick sheets of sediment, ice and woody debris¹¹⁻¹⁴. Extreme GLOF discharges can raise local river levels by several metres, increasing bed shear stresses and stream power values sufficiently to incise bedrock¹⁵. Floods from the largest lake outbursts can last for several days and destroy local flood control systems²,¹⁶. In many mountain regions, settlements, tourism, forestry and mining advanced into high mountain valleys during the twentieth century¹⁷ and thus closer to the source areas of GLOFs. As a consequence, ice-dam failures have caused many hundreds of fatalities, damaged infrastructure and farmland worth hundreds of millions of dollars, and disrupted transportation and communication routes<fnref n="1" />.

Despite the hazard posed by ice-dammed lakes, temporal trends in the magnitude and timing of their failures remain poorly understood in most mountain regions. We attribute the limited knowledge of contemporary GLOF activity to the largely unsystematic and regionally inconsistent record of GLOFs. How strongly, or whether at all, regional GLOF

activity responds to accelerating glacier melt therefore remains an open question. Melting glacier margins create new space for ice-dammed lakes to form and grow¹³,¹⁴,¹⁸. Thus, in the future, ice dams might store an increasing amount of meltwater, causing larger floods if they fail. However, the few detailed case studies of ice-dammed lakes, particularly in Iceland and northwest North America⁶,¹¹,¹⁴,¹⁹, suggest the opposite. There, individual ice-dammed lakes fill and empty repeatedly, releasing gradually smaller flood volumes over time. The edges of glacier dams might begin to float on the adjacent lake once a critical lake level is reached, suggesting that a gradually thinning glacier dam requires a lower threshold to initiate subglacial drainage⁶,²⁰,²¹. Within this cycle, the timing of GLOFs, that is, the day of the year when a lake begins to drain, is likely to shift, given that thinner dams may require a shorter period of time to refill to their maximum storage capacity¹⁴,²². This cycle of GLOFs from a given lake may end when the ice dam has decayed to the point where it can no longer impound water²⁰,²¹. Yet, new cycles of GLOFs might originate from higher elevations with ongoing exposure of ice-free areas at glacier margins²³.

## A database of historic ice-dam failures

Here we study the core components of the GLOF cycle by examining historical changes in GLOF magnitude, timing and source elevation, and their relationships to glacier thinning. We focus on the six most glaciated mountain regions on Earth, that is, northwest North America,

<footnote n="1">¹Institute of Environmental Science and Geography, University of Potsdam, Potsdam-Golm, Germany. ²Institute of Geosciences, University of Potsdam, Potsdam-Golm, Germany. ³Potsdam Institute for Climate Impact Research, Potsdam-Telegrafenberg, Germany. ⁴LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, Toulouse, France. ⁵Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zürich, Zürich, Switzerland. ⁶Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), Birmensdorf, Switzerland. ⁷Federal Institute for Materials Research and Testing (BAM), Berlin, Germany. ⁸Ministry of Forests, Prince George, British Columbia, Canada. ⁹Department of Earth Sciences, Simon Fraser University, Burnaby, British Columbia, Canada.

²e-mail: georg.veh@uni-potsdam.de</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Map of dated GLOF source locations in the six study regions. a, Global overview. b, Iceland. c, European Alps. d, Scandinavia. e, Northwest North America. f, Andes. g, High Mountain Asia, including legend. We created all maps in a–g using the free software QGIS V3.22.6 (https://www.qgis.org/). Generalized country boundaries and locations of cities are available through the online ArcGIS hub hosted by Esri Data and Maps (https://hub.arcgis.com/).</figcaption></figure>

High Mountain Asia, the Andes, Iceland, Scandinavia and the European Alps (Fig. 1). In reviewing 446 sources of information, we found a total of 1,569 dated outbursts from 186 ice-dammed lakes in the period 1900–2021 (Methods). In 64% of all cases, we had access to primary sources such as original research papers, official documentation from local agencies, reports in newspapers and social media platforms, and written correspondence with local eyewitnesses and experts. Secondary sources (36%) include review articles or data in previously published databases, for which we could not access the primary source directly. For each GLOF, we recorded the name of the source lake and its parent glacier; the year and, if available, the month and day of occurrence; the elevation of the lake surface; and any previously gauged or estimated values of outburst peak discharge $Q_p$ and flood volume $V_0$. Using Bayesian hierarchical quantile regression models, we estimated global and

h, Typical setting of an ice-dammed lake: Summit Lake dammed by Salmon Glacier, British Columbia, Canada (see location in e). Oblique view shows the situation on 11 August 1989, when the lake was approximately 1-km wide. Summit Lake has been draining beneath Salmon Glacier (arrow) since 1961, and only a small lake exists today (Fig. 5c). Photograph by John J. Clague. Mts, Mountains.

regional trends in reported median and extreme discharges, defined as GLOFs exceeding the 50th and 90th percentiles of the conditional distributions of reported $Q_p$ and $V_0$. We also used regression models to assess changes in the regional timing and the median source elevation of ice-dam failures (Methods). We assessed the sensitivity of the posterior trends in GLOF magnitudes learned from data since 1900 against the period since 1990, as regional GLOF reporting improved in the late twentieth century²⁴ (Extended Data Fig. 1).

## Less extreme GLOFs since 1900

Our inventory shows that the median reported discharge and volume released by ice-dammed lakes between 1900 and 2021 were $767_{-756}^{+9,761} \text{ m}^3 \text{ s}^{-1}$ and $51_{-50}^{+1,100} \times 10^6 \text{ m}^3$, respectively (median and 95%
<figure><img src="image_5.png" /><figcaption>**Fig. 2 | Trends of ice-dammed GLOF discharges between 1900 and 2021.**</figcaption></figure>

Posterior trends from Bayesian quantile regression of reported $V_0$ (orange) and $Q_p$ (blue) with time. Dark and light colours are trends for the median and 90th percentile (that is, 10% largest), respectively, for 1900–2021; grey densities are

trends for 1990–2021. Black lines are 95% HDIs and circles are medians of the posterior distributions. Trends refer to $\log_{10}$-transformed values of $V_0$ and $Q_p$, such that the posterior densities show the decadal trend in orders of magnitudes. Numbers ($n$) on the left are sample sizes for 1900–2021.

highest density interval (HDI)). Runoff peaks more than ten times higher than this global median occurred in High Mountain Asia where glacier dams collapsed suddenly after phases of advance³,⁷. Dam failures following fjord closures in the Andes²⁵ and Alaska² produced floods with volumes exceeding $10^9$ m³. However, our quantile regression models suggest that extreme GLOFs have weakened over the past 120 yr (Fig. 2 and Extended Data Figs. 2–4). The 90th percentile decreased by at least an order of magnitude in five regions for $Q_p$, and in four regions for $V_0$. The Andes and northwest North America had a decline of almost two orders of magnitude in the extremes of $Q_p$, whereas the highest $V_0$ dropped even more in the Andes and Iceland. Only $Q_p$ in Scandinavia and $V_0$ in the European Alps had unchanged trends in the extremes. The regional medians remained largely stable since 1900: $Q_p$ decreased credibly only in High Mountain Asia and northwest North America, without clear trends elsewhere. Four regions had a credible decrease in median $V_0$, yet in most cases trends in median $V_0$ were smaller than those of the extremes (Fig. 2 and Extended Data Fig. 3). A trend towards smaller floods is evident also for individual lakes with repeat (at least five) outbursts: hierarchical regression models show credible decreases in median $Q_p$ in 23% (5 of 22), and in median $V_0$ in 31% (8 of 26) of all cases. Only one lake had a credible increase in $Q_p$ or $V_0$ over the GLOF cycle (Extended Data Fig. 5).

All regions had more abundant reported GLOFs since 1990 compared with earlier decades (Extended Data Fig. 1). These changes in data density probably reflect a more consistent documentation of GLOFs²⁴, such that an unknown number of lake outbursts may have gone unnoticed early in the study period. Reducing the observation period to the past three decades (1990–2021), a period with greater research activity, yields broader posterior trend estimates as fewer data points enter the model (Fig. 2). However, in almost all regions, the trends obtained from shorter periods are consistent with the long-term trends, even though based on less than two-thirds of the total data (Extended Data Figs. 2 and 3). Trends of extreme $Q_p$ remain credibly negative in the Andes, High Mountain Asia and northwest North America, whereas median trends are unchanged in five of the six regions. Recent outbursts of

Lac de Faverges (Switzerland) may explain the only increase in the 10% highest records of $Q_p$ in the European Alps since 1990 (Fig. 2), although reported values of $Q_p$ in this region are very low (<100 m³ s⁻¹) by global standards²⁶,²⁷. Regional trends in extreme $V_0$ since 1990 mirror those of extreme $Q_p$, and are negative in five of the six regions. Except for a negative trend in High Mountain Asia, all regional trends in median $V_0$ remain unchanged. By explicitly drawing on stratified subsets of the data, we conclude that regional trends in GLOF size over the recent shorter time interval with probably better documentation hardly deviate from the long-term trends. The recent decrease in extremes is striking because these arguably represent the best documented events with large values of $Q_p$ and $V_0$, making a potential increase in extreme GLOFs with time very unlikely.

## Earlier annual timing of GLOFs

Most GLOFs happen in summer, and winter remains largely free of GLOFs (Extended Data Fig. 6). In the Northern Hemisphere, GLOF activity peaks in July and August, and in the Southern Hemisphere in January when monthly average temperatures are highest, thus supplying ice-dammed lakes with large amounts of meltwater. Yet, the average date of outbursts within a given year has shifted earlier globally by about 40 days compared with 1900. The strongest regional shifts in timing are in High Mountain Asia (77±25 days earlier), the European Alps (70±24 days earlier) and northwest North America (52±25 days earlier) (Fig. 3). In contrast, GLOFs in Iceland occur about six weeks later in the year on average, with the physical causes yet to be determined. At the site scale, we find that 68% of all ($n = 40$) ice-dammed lakes for which there is a documented day of the year (doy) burst out earlier in the year on average; in 30% of these cases, trends are credibly negative (Extended Data Fig. 7). The strongest change occurred at Lake Gornersee (Switzerland), which bursts out almost four months (114±19 days) earlier today than at the beginning of the twentieth century²². Successively smaller and earlier floods could signal a local weakening, if not the end, of the GLOF cycle²¹. Indeed, already half (46) of the 101 glaciers
<figure><img src="image_6.png" /><figcaption>**Fig. 3 |** *Changes in the annual timing of ice-dam failures between 1900 and 2021.* *a, Regional trends in the annual date of GLOFs. Thick black lines and light blue shades show, respectively, the median and 95% HDI of changes in the doy when glacier lakes burst out. Numbers in lower left corners are sample sizes.*</figcaption></figure>

*b, Posterior differences in GLOF timing (doy) between 1900 and 2021. Black bubbles are medians and horizontal lines are 95% posterior HDIs. Probability masses of negative trends are dark.*

with two or more failures stopped releasing GLOFs before 2000. Occasionally, human intervention interrupts the natural GLOF cycle, especially in cases where outbursts have repeatedly damaged infrastructure. In the European Alps and in Scandinavia, engineered structures such as drainage tunnels through bedrock or glacier ice stopped the recurrent outbursts of at least six ice-dammed lakes²⁸,²⁹.

## Increasing elevation of ice-dam failures

Our database further shows that outbursts from ice-dammed lakes originated at progressively higher elevations over the past 12 decades (Fig. 4). The largest elevation shifts are in Scandinavia, Iceland and the Andes, where new sources of GLOFs moved 25–50 m higher each decade on average. Outlet glaciers in these regions retreated rapidly in past decades, eliminating dams at their toes²⁵,³⁰,³¹. GLOFs in the European Alps had no credible change in elevation during our study period, although many large dams formed at lower elevations before 1900, when glaciers descended down to major valleys². These dams disappeared with warming and glacier retreat after the mid-nineteenth century². In contrast, High Mountain Asia has had more than 100 valley-blocking glacier advances to lower elevations in the past three decades³,³², explaining the unchanged elevation trend of ice-dam failures in that region. In most regions, outbursts emerged from lakes that were tens to many hundreds of metres lower in elevation than today (Fig. 4 and Extended Data Fig. 8). In three regions, the number of ice-dammed lakes decreased over time: Scandinavia, the European Alps and Iceland have, in total, only 25 ice-dammed lakes >0.01 km² today, whereas three times as many lakes burst out in the past 120 yr (Fig. 4).

## Impact of glacier decay on GLOF size

How does overall glacier decay help explain changes in ice-dam failures? In theory, the drainage of ice-dammed lakes occurs when the water pressure in the lake is close to or exceeds the pressure of the adjacent ice dam³³⁻³⁵. A GLOF may initiate once a critical ratio between lake depth

and glacier thickness is reached²⁰. As glaciers thin due to atmospheric warming⁸, the necessary water depth to trigger a GLOF should lower, resulting in smaller lakes that produce floods of decreasing magnitudes²⁰. We searched for a possible correlation between glacier thinning and GLOF magnitudes by focusing on 12 glaciers (one glacier impounds two lakes), for which we could obtain the cumulative elevation change of the ice dam, maximum lake areas before outbursts and reported GLOF magnitudes ($V_0$ or $Q_p$) between 2000 and 2019 (Methods). Over this period, we observed a thinning of the glacier dams by $36_{-62}^{+19}$ m (median and 95% HDI; Fig. 5c and Extended Data Fig. 9), and that 12 of 13 lakes shrank, with 8 losing more than 20% of their surface areas (Fig. 5c). Yet, rapid glacier thinning appears to correlate only moderately with reported GLOF magnitudes: average trends between glacier thinning and decreases in flood volumes or peak discharges are positive for only half of the glacier dams (Fig. 5a,b). Credibly positive trends occur only at Tulsequah Lake and Lake No Lake (British Columbia), and at Lago Cachet (Patagonia). Most other lakes have no clear trend. Interactions between ice and water might be an important mechanism for removing parts of glacier dams, causing most ice-dammed lakes to shift laterally towards the glacier dam during the GLOF cycle. For example, Lake No Lake and Summit Lake (northwest North America) eroded several hundreds of metres from their ice dams in the past 20 yr (Fig. 5c).

The unclear relationship between changes in glacier thickness and GLOF magnitudes suggests a nonlinear response, or instead that glacier thinning is a poor proxy for GLOF magnitudes. Variances among the estimated local trends point to local drainage mechanisms. Small variances in $Q_p$ and $V_0$ (ref. ²⁷) may indicate recurrent opening and closing of the same subglacial tunnel system that allows the floods to attain high peak flows despite decreases in glacier thickness²². Floating ice can further exert high pressures on the underlying water in lakes, maintaining high $Q_p$ and $V_0$ independent of glacier thinning that reduces the basin's storage capacity³⁶. In contrast, high variances or abrupt changes in the trends may indicate shifts between drainage regimes, for example due to failure of ice buttresses at the point of inflow into the glacier¹⁹. Ice motion during an outburst can also close subglacial
<figure><img src="image_7.png" /><figcaption>**Fig. 4 |** *Elevation changes of GLOF sources.*** a, Changes in elevation of ice-dam failures over time. Blue bubbles mark years of the first reported outburst and the elevations above sea level (a.s.l.) of ice-dammed lakes. Thick black lines and green shades are the median and 95% HDI of the posterior trends, respectively. Blue numbers in lower left corners are the number of lakes that burst out between 1900 and 2021. Orange ticks and numbers on the right are elevations and the number of present-day ice-dammed lakes >0.01 km² in each region,</figcaption></figure>

respectively⁴⁸,⁵⁵⁻⁵⁸. The grey densities estimate present-day glacier volumes across elevation from ref. ⁴⁷, truncated at 95% of the data. The black areas within the grey densities show the loss of glacier volume over the period 2000–2019 from ref. ⁸, relative to the entire glacier volume. b, Posterior decadal trends in annual elevation change of GLOF sources. Black bubbles are medians and horizontal lines are 95% posterior HDIs. Probability masses of positive trends are dark.

tunnels before the lake has emptied completely³⁴. Some lakes such as Lago Cachet (Fig. 5c) have maintained their surface area but gained depth as sediment was removed from the lake bottom between successive outbursts¹². Outburst floods can also drain other water bodies within and beneath the glacier. The volume from these sources can rival or even exceed those of the subaerial lake, although distinguishing the contributions from the individual sources is subject to high uncertainties³⁷,³⁸. Local weather can also add to the observed variability. Water supply to lakes can suddenly increase following heavy precipitation or during heatwaves when there is high glacier and snow melt, resulting in unusual high peak flows¹⁶,³⁹. Accurately determining such local water fluxes could make GLOFs more predictable²²,³⁷,⁴⁰, although it would require a much denser meteorological network to quantify weather conditions in ice-dammed basins⁴¹.

## Implications for GLOF hazard

In our study regions, >70% of all reported GLOFs since 1900 originated from ice-dammed lakes²⁴. However, their present-day number and sizes are small compared with meltwater lakes dammed by moraines or bedrock. In northwest North America, ice-dammed lakes account only for 6% of the total lake abundance¹⁸, and in High Mountain Asia, as few as ~0.06% of total lake area is impounded by glaciers⁴². However, given their high outburst activity, we infer that ice-dammed lakes are among the most hazardous of all meltwater lakes. Still, catastrophic damage and fatalities remain a rare consequence of ice-dammed failures: only 10% of all floods in our database have documented socio-economic impacts, mostly involving destroyed bridges, pipelines, roads and

campgrounds buried in debris, or loss of farmland and cattle affecting local food supply. Agricultural loss prevails in reports from Scandinavia, Iceland, the Andes and the European Alps. Since 1900, three ice dam failures have claimed human lives, although loss of life was reported more frequently in previous centuries. Well documented cases include the sudden failures of an ice-dammed lake in the Dranse Valley, Swiss Alps, in 1590 and 1818 that killed at least 40 and 140 people, respectively⁴³. Such one-time ice-dam failures and those starting a new, unexpected GLOF cycle with flood defence measures not yet in place have caused much damage¹⁶,⁴⁴. Dykes, floodgates, dams, or the temporary or permanent resettlement of mountain communities have reduced the impact of recurrent ice-dam failures, even for exceptionally large floods¹.

The global reduction in glacier volumes⁸ will probably control ice-dam failure in the future. We find that 90% of all outburst floods in our database originated below the elevation band where most glacier mass is stored today (Fig. 4a). Elevations of historical ice-dam failures are largely within glacier ablation areas, where glaciers melt fastest and temporarily trap excess meltwater at their margins. This zone has been rapidly shifting upwards under ongoing atmospheric warming. In the European Alps and the Himalayas, for example, glaciers have been melting throughout their entire elevation range in the twenty-first century⁴⁵,⁴⁶, and ice-dammed lakes are now forming even in the highest parts of ice-covered catchments (Fig. 4a and Extended Data Fig. 8). Progressive retreat at the glacier terminus may thus push the accommodation space for new lakes towards the natural limit set by the remaining glacier mass close to the drainage divide. According to projections, regions with small amounts and volumes of glaciers, such as the European Alps, Scandinavia and British Columbia, could
a

<figure><img src="image_8.png" /></figure>

b

<figure><img src="image_9.png" /></figure>

**Fig. 5 | Trends in flood volume, peak discharge and lake area with cumulative changes in glacier thickness between 2000 and 2019.** a, Map of 15 ice-dammed lakes that produced more than five GLOFs each in that period. We generated the map in the free software RStudio V2022.07.2 (https://posit.co/products/open-source/rstudio/) using the ggmap package, including free data on administrative boundaries from the naturalearth package. b, Posterior regression slopes of $V_0$ (left) and $Q_p$ (right) versus cumulative change in thickness of the ice dam. Bubbles are medians and horizontal lines are 95% posterior HDIs of the posterior trends. Colours in a and b distinguish the study regions. Missing trends for individual lakes arise from lack of data. c, Local changes in lake area

become largely ice free by 2100⁹. We anticipate that the window for future ice-dam formation and failure might close within a few decades in these mountain ranges. Glacier remnants in steeper and higher terrain may be capable of only damming, if at all⁵, smaller lakes, albeit with higher potential energy. In the course of one to several centuries, ice-dam failures may diminish as a cryospheric hazard, although mass flows from other lakes, possibly tied to permafrost thaw, may take over.

This trajectory of changing GLOF hazard is likely to differ in regions with large glacier volumes. Alaska and Patagonia, for example, host ~70% of the total glacier volume in our study regions, about half of which will persist beyond the end of the twenty-first century⁹,⁴⁷. Extensive glacier streams currently terminate at or near sea level, and hundreds of glacier-dammed lakes have been forming in ice-free tributary valleys in recent decades¹⁸,⁴⁸ (Fig. 4a). Yet, statistical models have failed to identify glaciological or meteorological predictors that could explain

c

<figure><img src="image_10.png" /></figure>

and glacier elevation between 2000 and 2019. We mapped lake areas (polygons) from satellite imagery immediately before the outburst using QGIS V3.22.6. Colours show the year of the GLOFs. Blue bubbles show the average percent change in lake area between the first and last reported GLOF for each lake. Orange bubbles show the mean annual elevation change of the glacier dam. Black horizontal scale bar is 1 km in all panels in c. Data on lake and glacier elevation change are available at ref. ²⁷. All background imagery are from 2019 Planet Labs Inc., obtained under a basic license in the Education and Research programme. We generated all maps in c using QGIS V3.22.6.

the growth pattern of ice-dammed lakes in northwest North America⁴⁹. Trends in the magnitude and timing of GLOFs, let alone any drivers, are largely unknown for most of these lakes²³. Newly exposed proglacial areas could be sources of new cycles of GLOFs, posing hazards for potential development in tourism, mining and hydropower generation in these areas. In this regard, time series of satellite images⁵⁰ and stream gauges¹²,¹⁴,²² are invaluable archives to systematically capture geomorphic changes caused by GLOFs. Integrating newly detected and reported cases into a single database is thus key to quantifying and predicting climate-driven changes in GLOF magnitude, frequency and hazard²⁴.

In summary, our analysis reveals a strong decay in the highest GLOF magnitudes and discharges, whereas GLOFs of median size have changed little in the past 12 decades. Given their steep hydrographs, floods caused by ice-dam failures will remain discharge extremes, even during glacier retreat. How and where the remaining glaciers will form
dams in the future will determine the extent to which GLOFs will feed the sediment cascade¹², disturb ecological and geochemical cycles in mountain rivers⁵¹,⁵², and locally inundate floodplains and terraces⁷. To anticipate these consequences, future work will be required to refine the few available physical models that link glacier retreat and lake formation with the size and timing of outbursts¹⁰,³⁴⁻³⁶. Recently released high-resolution data on elevation, glacier change, ice-flow velocities and climate will be key to understanding better the nexus between meltwater accumulation and hazard. Predicting and monitoring the locations and dynamics of future lakes will become increasingly important, especially in steep mountain regions that are difficult to access, but have growing populations and infrastructure downstream. International collaborations such as the United Nations Development Programme project 'GLOF-I/-II' in Pakistan aim to mitigate GLOF hazards, but more effort will be required to manage the episodic GLOF discharges as the variability in precipitation, temperature and river discharge rises⁵³,⁵⁴. Our dataset forms a solid base to help predict changes in GLOF hazard from the catchment to the regional scale to ensure more sustainable development for one of the most sensitive regions on Earth.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05642-9.

1. Carrivick, J. L. & Tweed, F. S. A global assessment of the societal impacts of glacier outburst floods. *Glob. Planet. Change* **144**, 1–16 (2016).

2. Rabot, C. Glacial reservoirs and their outbursts. *Geogr. J.* **25**, 534–548 (1905).

3. Bazai, N. A. et al. Increasing glacial lake outburst flood hazard in response to surge glaciers in the Karakoram. *Earth Sci. Rev.* **212**, 103432 (2021).

4. Häberli, W. Frequency and characteristics of glacier floods in the Swiss Alps. *Ann. Glaciol.* **4**, 85–90 (1983).

5. Liestel, O. Glacial dammed lakes in Norway. *Nor. Geogr. Tidsskr.* **15**, 122–149 (1956).

6. Thorarinsson, S. Chapter IX. The ice dammed lakes of Iceland with particular reference to their values as indicators of glacier oscillations. *Geogr. Anna.* **21**, 216–242 (1939).

7. Hewitt, K. & Liu, J. Ice-dammed lakes and outburst floods, Karakoram Himalaya: historical perspectives on emerging threats. *Phys. Geogr.* **31**, 528–551 (2010).

8. Hugonnet, R. et al. Accelerated global glacier mass loss in the early twenty-first century. *Nature* **592**, 726–731 (2021).

9. Hock, R. et al. GlacierMIP—a model intercomparison of global-scale glacier mass-balance models and projections. *J. Glaciol.* **65**, 453–467 (2019).

10. Roberts, M. J. Jökulhlaups: a reassessment of floodwater flow through glaciers. *Rev. Geophys.* **43**, 21 (2005).

11. Post, A. & Mayo, L. R. *Glacier Dammed Lakes and Outburst Floods in Alaska*. Hydrologic Atlas 455 (US Geological Survey, 1971); https://doi.org/10.3133/ha455.

12. Jacquet, J. et al. Hydrologic and geomorphic changes resulting from episodic glacial lake outburst floods: Rio Colonia, Patagonia, Chile. *Geophys. Res. Lett.* **44**, 854–864 (2017).

13. Mathews, W. H. & Clague, J. J. The record of jökulhlaups from Summit Lake, northwestern British Columbia. *Can. J. Earth Sci.* **30**, 499–508 (1993).

14. Geertsema, M. & Clague, J. J. Jökulhlaups at Tulsequah Glacier, northwestern British Columbia, Canada. *Holocene* **15**, 310–316 (2005).

15. Larsen, I. J. & Lamb, M. P. Progressive incision of the Channeled Scablands by outburst floods. *Nature* **538**, 229–232 (2016).

16. Tvede, A. M. Floods caused by a glacier-dammed lake at the Folgefonni Ice Cap, Norway. *Ann. Glaciol.* **13**, 262–264 (1989).

17. Wang, Y., Wu, N., Kunze, C., Long, R. & Perlik, M. in *The Hindu Kush Himalaya Assessment: Mountains, Climate Change, Sustainability and People* (eds Wester, P. et al.) 17–56 (Springer, 2019); https://doi.org/10.1007/978-3-319-92288-1-2.

18. Rick, B., McGrath, D., Armstrong, W. & McCoy, S. W. Dam type and lake location characterize ice-marginal lake area change in Alaska and NW Canada between 1984 and 2019. *Cryosphere* **16**, 297–314 (2022).

19. Roberts, M. J., Pálsson, F., Gudmundsson, M. T., Björnsson, H. & Tweed, F. S. Ice-water interactions during floods from Grænآلon*g*lacier-dammed lake, Iceland. *Ann. Glaciol.* **40**, 133–138 (2005).

20. Tweed, F. S. & Russell, A. J. Controls on the formation and sudden drainage of glacier-impounded lakes: implications for Jökulhlaup characteristics. *Prog. Phys. Geogr.* **23**, 79–110 (1999).

21. Evans, S. G. & Clague, J. J. Recent climatic change and catastrophic geomorphic processes in mountain environments. *Geomorphology* **10**, 107–128 (1994).

22. Huss, M., Bauder, A., Werder, M., Funk, M. & Hock, R. Glacier-dammed lake outburst events of Gornersee, Switzerland. *J. Glaciol.* **53**, 189–200 (2007).

23. Capps, D. M. & Clague, J. J. Evolution of glacier-dammed lakes through space and time; Brady Glacier, Alaska, USA. *Geomorphology* **210**, 59–70 (2014).

24. Veh, G. et al. Trends, breaks, and biases in the frequency of reported glacier lake outburst floods. *Earths Future* **10**, e2021EF002426 (2022).

25. Dussaultant, I. et al. Two decades of glacier mass loss along the Andes. *Nat. Geosci.* **12**, 802–808 (2019).

26. Raymond, M., Wegmann, M. & Funk, M. *Inventar Gefährlicher Gletscher in der Schweiz*. Report 182 (Versuchsanstalt für Wasserbau, Hydrologie und Glaciologie der Eidenössischen Technischen Hochschule Zürich, 2003).

27. Veh, G. Supplementary data for 'Less extreme and earlier outbursts of ice-dammed lakes since 1900'. *Zenodo* https://doi.org/10.5281/zenodo.7326571 (2022).

28. Häberli, W., Kääb, A., Mühll, D. V. & Teysseire, P. Prevention of outburst floods from periglacial lakes at Grubengletscher, Valais, Swiss Alps. *J. Glaciol.* **47**, 111–122 (2001).

29. Jackson, M. & Ragulina, G. *Inventory of Glacier-Related Hazardous Events in Norway*. Report 83 (Norwegian Water Resources and Energy Directorate, 2014).

30. Guðmundsson, S. et al. Terminus lakes on the south side of Vatnajökull Ice Cap. *SE Iceland. Jökull* **69**, 1–34 (2020).

31. Andreassen, L. M., Elvehøy, H., Kjellmoen, B. & Belart, J. M. C. Glacier change in Norway since the 1960s—an overview of mass balance, area, length and surface elevation changes. *J. Glaciol.* **66**, 313–328 (2020).

32. Vale, A. B., Arnold, N. S., Rees, W. G. & Lea, J. M. Remote detection of surge-related glacier terminus change across High Mountain Asia. *Remote Sens.* **13**, 1309 (2021).

33. Röthlisberger, H. Water pressure in intra- and subglacial channels. *J. Glaciol.* **11**, 177–203 (1972).

34. Nye, J. F. Water flow in glaciers: Jökulhlaups, tunnels and veins. *J. Glaciol.* **17**, 181–207 (1976).

35. Fowler, A. C. Breaking the seal at Grímsvötn, Iceland. *J. Glaciol.* **45**, 506–516 (1999).

36. Jenson, A., Amundson, J. M., Kingslake, J. & Hood, E. Long-period variability in ice-dammed glacier outburst floods due to evolving catchment geometry. *Cryosphere* **16**, 333–347 (2022).

37. Bigelow, D. G. et al. The role of englacial hydrology in the filling and drainage of an ice-dammed lake, Kaskawulsh Glacier, Yukon, Canada. *J. Geophys. Res. Earth Surf.* **125**, e2019JF005110 (2020).

38. Shangguan, D. et al. Quick release of internal water storage in a glacier leads to underestimation of the hazard potential of glacial lake outburst floods from Lake Merzbacher in central Tian Shan Mountains. *Geophys. Res. Lett.* **44**, 9786–9795 (2017).

39. Ng, F., Liu, S., Mavlyudov, B. & Wang, Y. Climatic control on the peak discharge of glacier outburst floods. *Geophys. Res. Lett.* **34**, L21503 (2007).

40. Kingslake, J. & Ng, F. Quantifying the predictability of the timing of jökulhlaups from Merzbacher Lake, Kyryzstan. *J. Glaciol.* **59**, 805–818 (2013).

41. Kienholz, C. et al. Deglacierization of a marginal basin and implications for outburst floods, Mendenhall Glacier, Alaska. *Front. Earth Sci.* **8**, 137 (2020).

42. Maharjan, S. B. et al. *The Status of Glacial Lakes in the Hindu Kush Himalaya Research Report 2018/1* (ICIMOD, 2018).

43. Ancey, C. et al. Hydraulic reconstruction of the 1818 Giérotral glacial lake outburst flood. *Water Ressour.* **55**, 8840–8863 (2019).

44. Mathews, W. H. Two self-dumping ice-dammed Lakes in British Columbia. *Geogr. Rev.* **55**, 46–52 (1965).

45. Fischer, M., Huss, M. & Hoelzle, M. Surface elevation and mass changes of all Swiss glaciers 1980–2010. *Cryosphere* **9**, 525–540 (2015).

46. Maurer, J. M., Schaefer, J. M., Rupper, S. & Corley, A. Acceleration of ice loss across the Himalayas over the past 40 years. *Sci. Adv.* **5**, eaav7266 (2019).

47. Farinotti, D. et al. A consensus estimate for the ice thickness distribution of all glaciers on Earth. *Nat. Geosci.* **12**, 168–173 (2019).

48. Wilson, R. et al. Glacial lakes of the Central and Patagonian Andes. *Glob. Planet. Change* **162**, 275–291 (2017).

49. Field, H. R., Armstrong, W. H. & Huss, M. Gulf of Alaska ice-marginal lake area change over the Landsat record and potential physical controls. *Cryosphere* **15**, 3255–3278 (2021).

50. Veh, G., Korup, O., Roessner, S. & Walz, A. Detecting Himalayan glacial lake outburst floods from Landsat time series. *Remote Sens. Environ.* **207**, 84–97 (2018).

51. Anderson, S. P. et al. Integrated hydrologic and hydrochemical observations of Hidden Creek Lake jökulhlaups, Kennicott Glacier, Alaska. *J. Geophys. Res.* **108**, 6003 (2003).

52. Milner, A. M. et al. River ecosystem resilience to extreme flood events. *Ecol. Evol.* **8**, 8354–8363 (2018).

53. Huss, M. et al. Toward mountains without permanent snow and ice. *Earths Future* **5**, 418–435 (2017).

54. Khanal, S. et al. Variable 21st century climate change response for rivers in High Mountain Asia at seasonal to decadal time scales. *Water Res.* **57**, e2020WR029266 (2021).

55. Nagy, T. & Andreassen, L. *Glacier Lake Mapping with Sentinel-2 Imagery in Norway* NVE Rapport 40-2019 (NVE, 2019).

56. Buckel, J., Otto, J. C., Prasicek, G. & Keuschnig, M. Glacial lakes in Austria—distribution and formation since the Little Ice Age. *Glob. Planet. Change* **164**, 39–51 (2018).

57. Chen, F. et al. Annual 30 m dataset for glacial lakes in High Mountain Asia from 2008 to 2017. *Earth Syst. Sci. Data* **13**, 741–766 (2021).

58. Mölg, N. et al. Inventory and evolution of glacial lakes since the Little Ice Age: lessons from the case of Switzerland. *Earth Surf. Process. Landf.* **46**, 2551–2564 (2021).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_11.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023
### Developing a database of ice-dam failures

We define an ice-dammed lake as any water body at the margin of and impounded by glacier ice. We exclude all other glacier-fed lakes, such as lakes dammed by moraines or landslides, internal water pockets, supraglacial ponds, lakes below glaciers and lakes formed by geothermal activity. Ice dams often trap runoff from catchments and cause loss and damage when the impounded water is suddenly released and reaches human settlements³⁹. We focus on six major glaciated mountain regions that support 99% of the contemporary ice volume outside polar regions⁴⁷. Our evidence for the 1,569 ice dam failures reported between 1900 and 2021 comes from a total of 446 different sources of information. We distinguish between primary sources, where we have direct observations or access to raw data, and secondary sources that describe, interpret or summarize previous findings on historical GLOFs. Primary sources are available for 64% of all reported cases and include original scientific publications (36%), annals and reports published by local and regional authorities (18%), written communications in emails with eyewitnesses and experts (9%), and newspapers and internet resources such as video footage and posts on social media platforms (altogether 1%). Secondary sources (36% of all cases) include second-hand information derived from reviews and summaries of previous publications (19%) and entries in databases with no direct access to the primary source (17%). Half of all sources report only a single event, whereas a few studies include information on dozens of GLOFs. For example, ref. ⁶⁰ identified more than 100 outbursts from lakes impounded by Kennicott Glacier during the twentieth century by studying stream gauge data and historical aerial photographs. Still, decades may elapse between the date of a GLOF and the reporting, especially where researchers reprocessed and reanalysed historical data (Extended Data Fig. 1). In half of all cases, however, GLOFs were reported within a decade after they occurred. A few studies contributed significantly to improving the regional record of GLOFs in the first decades of the twentieth century, when the overall reporting rate was low³,⁷,¹¹,⁶⁰⁻⁶². In a few cases, reporting may have changed due to resettlement in this period. For example, repeated flooding of an ice-dammed lake below the Folgefonni Ice Cap (Norway) forced farmers to abandon their homes in the mid-1960s¹⁶. It is only since the 1970s that research activity on GLOFs has become more consistent²⁴: 63% GLOFs in our databases were first documented in that period (Extended Data Fig. 1).

In collecting historic GLOFs, we only considered cases with either a documented date (at least the year of occurrence) or a time interval of failure, centroidal coordinates of the lake and at least one reference. For each GLOF, we documented the mountain range, the country from which the flood originated, the name of the parent glacier (local name and the ID in the Randolph Glacier Inventory RGI V6.0³⁵) and the name of the burst glacier lake. We extracted the surface elevation of each lake centroid from the Advanced Land Observing Satellite (ALOS) V3.2 digital elevation model (DEM)⁶⁴. We collected the reported flood volume (n = 510, 33% of all cases), peak discharge (n = 506, 33%) and timing (n = 940, 61%). Finally, we documented any impacts, damage, and losses for events in our database.

### Lake mapping

We could verify the origin, timing and area of 53% of all reported lake outbursts that have occurred since the Landsat 5 mission began to continuously obtain satellite images in the late-1980s⁶³. Later generations of satellites, such as RapidEye and the Planet Cubesat constellation (launched in 2009 and 2015, respectively), have higher spatial resolution and thus capture GLOF features in greater detail. Ice-dammed lakes decrease in size, or even empty entirely, during GLOFs, exposing either partly or fully the lake floor, which is often covered with stranded icebergs (Fig. 5c). Ice, clouds, and shadows compromise the automatic detection of water and lakes⁶⁶. Hence, we manually digitized the extents

of the lakes from satellite images in QGIS V3.16 software. To this end, we scanned the archive of satellite images to retrieve the last available image before and the first available image after the lake drained. Mapping lakes from optical images was limited in Iceland because of frequent cloud cover. For each lake with repeat outbursts, we used a simple linear regression model of the mapped lake area before the outbursts versus time to obtain the change in lake area (Fig. 5). We then estimated the mean percent change in lake area between the first and last reported outburst in the time series. For manually mapped lake areas, uncertainties are generally assumed to be between 4% and 10% of the mapped lake area⁴⁸,⁵⁵ or between 0.5 and 1 pixel of the sensor resolution¹⁸,⁴²,⁵⁷, that is, 3–5 m for Planet and RapidEye, and 30 m for Landsat images.

In Fig. 4, we compare elevations of lakes that have had historic outbursts with elevations of ice-dammed lakes that currently (2018–2020) exist according to regional inventories¹⁸,⁴⁸,⁵⁵⁻⁵⁸. The lakes in these databases were either manually digitized or refined after automatic mapping. We chose a minimum mapping unit of 0.01 km² to ensure comparability between inventories. No contemporary lake inventory is available for Iceland. Thus, we manually mapped ice-dammed Icelandic lakes from high-resolution Google Earth imagery obtained between May and September 2021. We added the elevations obtained from the ALOS V3.2 DEM⁶⁴ to both the previous lake inventories and our manually mapped lakes.

### Bayesian hierarchical modelling

To assess temporal changes in GLOF magnitude, elevation and timing, we fitted Bayesian hierarchical models that distinguish between regional (that is, our six study regions) or local (that is, lakes with repeat outbursts) groups. Bayesian models form a compromise between the likelihood of observing the data under the assumption of specific model parameters and prior knowledge about these parameters. We encode our uncertainty in this prior knowledge before obtaining a posterior distribution for all parameters conditioned on the data. The benefit of hierarchical models is that they account for group-level structure in the data within a single model. We can thus estimate regional or local effects with respect to a population-level model learned from all data regardless of their location. Hierarchical models offer a compromise between separate models that overemphasize regional and local effects and that are independently learned for each group, and pooled models that are learned from all data without distinction⁶⁷. Hierarchical models improve parameter estimates for individual groups, especially if sample sizes are low or vary across the groups. This trait is particularly desirable for our study as the number of reported GLOFs differs widely among both regions and lakes. The model learns the population-level parameters from the data, and these parameters are shared among groups. In this way, the groups inform each other, achieving an effective sample size that greatly exceeds that of single groups. The joint posterior density in a hierarchical model is proportional to the product of the likelihood and the joint prior density $\pi(\theta_1, ..., \theta_j, \tau|D) \propto \pi(D|\theta_1, ..., \theta_j, \tau)\pi(\theta_1, ..., \theta_j, \tau) = \pi(\tau)\prod_{j=1}^{J} \pi(D_j|\theta_j)\pi(\theta_j|\tau)$, where D are the observed data and $D_j$ is the subset of D that contains all data in group j. The group-level parameters $\theta_j = (\theta_{j1}, ..., \theta_{jK_j})$ define the model on the group-level with K individual parameters, which vary for j = 1, ..., J groups. These parameters are drawn from distributions $\pi(\theta_j | \tau)$ and specified by population-level (hyper-) parameters $\tau = (\tau_1, ..., \tau_L)$, which consist of L individual parameters and are also learned from the data. The hierarchical model estimates the posterior distribution for both population-level and group-level parameters conditional on the reported GLOFs.

### Assessing regional and local trends of V₀, Qₚ and Z

Continental and global studies of river floods point to climate-driven changes in both median and extreme discharges in past decades⁶⁸⁻⁷¹. Similarly, repeated inventories reveal recent shifts in the regional size
distribution of ice-dammed lakes, but little about trends in lake outbursts¹⁸,⁴⁸,⁴⁹. By using quantile regression, we examined whether and by how much the quantiles in the conditional distributions of reported peak discharge, Qp, and flood volume, Vo, have changed in our six study regions. Extreme (>90th percentile) flood discharges are important for planning flood protection measures and informing flood risk zoning⁶⁹. Ordinary least squares regression estimates the conditional mean, but quantile regression makes no assumptions about the distribution of the response variable and is robust against outliers. In our model specification, we used an asymmetric Laplace likelihood (AL) with fixed quantiles. The AL has density

$$f_{\text{AL}}(y|\mu, \kappa, p) = \frac{p(1-p)}{\kappa} \exp\left(-\rho_p\left(\frac{y-\mu}{\kappa}\right)\right), \quad (1)$$

where y is the response variable, μ is a location parameter, κ is a positive scale parameter, 0 < p < 1 is a percentile, and ρp(x) = x(p - I(x < 0)) with indicator function I(x). The exponent takes the values -|y - μ|/κ for y ≥ μ and -|y - μ|/κ(1 - p) for y < μ. We modelled the conditional distributions of Qp and Vo with decimal year t for two fixed values of p as

$$y_{ji} \sim \text{AL}(\mu_{ji}, \kappa, p), \text{ for } j = 1, \dots, J \text{ and } i = 1, \dots, n_j \quad (2)$$

$$\mu_{ji} = \alpha_j + \beta_j t_{ji}, \text{ for } j = 1, \dots, J \text{ and } i = 1, \dots, n_j \quad (3)$$

$$\begin{bmatrix} \alpha_j \\ \beta_j \end{bmatrix} \sim \text{MVNormal}\left[\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, S\right], \text{ for } j=1, \dots, J \quad (4)$$

$$S = \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{pmatrix} R \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{pmatrix} \quad (5)$$

$$R = \begin{pmatrix} 1 & \zeta \\ \zeta & 1 \end{pmatrix} \quad (6)$$

where yji are reported values of Qp or Vo with yji referring to the ith observation in group j, nj is the number of observations in group j, and J is the number of groups. The parameters αj and βj are the group-level intercepts and slopes, respectively, and α and β are the corresponding population-level parameters. The covariance matrix S is composed of group-level standard deviations σα and σβ, and R, the correlation matrix with correlation c. In line with decadal analyses of median and extreme flood discharges⁶⁸⁻⁷⁰, we considered fixed values of p = 0.5 (median) and p = 0.9 (90th percentile), and estimated all other parameters from the data.

The Bayesian framework demands prior distributions for each parameter that enters at another model level. We chose the following priors

$$\kappa \sim N(0, 2.5) \quad (7)$$

$$\alpha \sim N(0, 2.5) \quad (8)$$

$$\beta \sim N(0, 2.5) \quad (9)$$

$$\sigma_{\alpha} \sim N(0, 2.5) \quad (10)$$

$$\sigma_{\beta} \sim N(0, 2.5) \quad (11)$$

$$R \sim \text{LKJCholesky}(1). \quad (12)$$

These priors refer to standardized, log₁₀-transformed response variables Qp and Vo, and a standardized predictor (decimal year of

observation) with zero mean and unit standard deviation. Our choice of a zero-mean Gaussian with standard deviation of 2.5 admits both negative and positive trends for β, expressing the contrasting trends reported for the size of ice-dam failures in recent decades¹,⁶,¹⁹,²²,⁷². The Lewandowski-Kurowicka-Joe (LKJ) Cholesky correlation distribution prior for R makes all correlation matrices equally likely.

For a fixed quantile, the joint posterior distribution can be thus written as

$$\pi(\alpha_1, \dots, \alpha_j, \beta_1, \dots, \beta_j, \kappa, \alpha, \beta, \sigma_\alpha, \sigma_\beta, R|y, t, p) \propto \pi(\kappa, \alpha, \beta, \sigma_\alpha, \sigma_\beta, R) \\ \left(\prod_{j=1}^{J} \prod_{i=1}^{n_j} f_{\text{AL}}(y_{ji}|\alpha_j + \beta_j t_{ji}, p, \kappa) \pi(\alpha_j, \beta_j|\alpha, \beta, \sigma_\alpha, \sigma_\beta, R)\right), \quad (13)$$

where π(αj, βj|α, β, σα, σβ, R) is the prior distribution for the group-level parameters as defined in equations (4)–(6) and π(κ, α, β, σα, σβ, R) is the joint prior distribution for the independent population-level parameters as defined in equations (7)–(12). We numerically approximate the posterior distribution using a Hamiltonian sampling algorithm implemented in Stan⁷³ that is called via the software package brms⁷⁴ within the statistical programming language R⁷⁵.

We used the same multi-level structure of the quantile regression model to determine how flood discharges have changed for individual lakes. To this end, we conditioned the model on J = 22 lakes (Qp) and J = 26 lakes (Vo) that produced at least five outbursts each between 1900 and 2021. For these models, we focused only on the trends in median flood discharges, given that lakes with few reported GLOFs have few samples that exceed high thresholds.

Finally, to test whether more recent GLOFs originated from progressively higher elevations, we selected for each lake the year of the first reported GLOF, given that many ice-dammed lakes burst out repeatedly. As for our hierarchical models on GLOF magnitudes, we used quantile regression to determine the trend in the median elevation Z of glacier-dammed lakes against the decimal year of GLOF occurrence t, conditioned on the J = 6 regions.

We ran all regional and local quantile regression models on Qp, Vo and Z with three parallel chains that drew 6,000 samples with 2,000 warm-up runs each. We found no numerical divergences after warm-ups, which were supported by the Gelman-Rubin potential scale reduction factor R̂ = 1.0 in all models, indicating that the Markov chains have converged. We report the posterior distributions of the pooled model parameters in Supplementary Tables 1–3, and the group-level regression slopes in Extended Data Figs. 2, 3 and 5. These regression slopes are conventionally deemed 'credibly' different from zero if a select interval of the posterior probability mass, in our case the 95% HDI, excludes zero. This choice of the HDI is arbitrary, and we also report how much of the probability mass of the posterior distribution is on either side of zero (Extended Data Figs. 2–5).

## Assessing trends in doy

Measuring temporal changes in GLOF timing must consider that the response variable doy can recycle through the observation period, given that calendar day 0 is equal to day 365. The von Mises (VM) distribution is a close approximation to a wrapped normal distribution for circular data and allows estimating the angular response d (doy) from the linear predictor year t (ref. ⁶⁴). The VM distribution has density

$$f_{\text{VM}}(y|\theta, \varphi) = \frac{\exp(\varphi \cos(y - \theta))}{2\pi I_0(\varphi)}, \quad (14)$$

where y are the observed data, θ is a location parameter, φ is a positive precision parameter and I₀ is the modified Bessel function of the first kind of order 0.

Following the parameterization in brms, we rescaled the response d to yᵢ ∈ (−π, π), standardized t to a zero mean and unit standard deviation, and estimated the trend in GLOF timing as
$$y_{ji} \sim \text{VM}(\theta_{ji}, \phi_j), \quad \text{for } j=1, \dots, J \text{ and } i=1, \dots, n_j \quad (15)$$

$$\theta_{ji} = \zeta_j + \eta_j t_{ji}, \quad \text{for } j=1, \dots, J \text{ and } i=1, \dots, n_j \quad (16)$$

$$\zeta_j \sim N(0, 2.5) \quad (17)$$

$$\eta_j \sim N(0, 2.5) \quad (18)$$

$$\phi_j \sim \text{Gamma}(2, 0.01), \quad (19)$$

where $y_{ji}$ are reported values of doy, with $y_{ji}$ referring to the ith observation in group j, $\theta_{ji}$ is the corresponding location parameter, and $\zeta_j$ and $\eta_j$ are the estimated intercepts and slopes, respectively. We used the default Gamma prior in brms to estimate the precision parameter $\phi_j$. This model was fit to each group separately without partial pooling of parameters across groups, because the distinct seasonal differences in GLOF timing in the two hemispheres (Extended Data Fig. 6) caused numerical divergences. Using a VM likelihood function, the joint posterior distribution can be written as

$$\pi(\zeta_j, \eta_j, \phi_j | y_{ji}, t_{ji}) \propto \left( \prod_{i=1}^{n_j} f_{\text{VM}}(y_{ji} | \zeta_j + \eta_j t_{ji}, \phi_j) \right) \pi(\zeta_j, \eta_j, \phi_j), \quad (20)$$
for $j=1, \dots, J$

where $y_{ji}$ and $t_{ji}$ are the observations and corresponding years of group j and $\pi(\zeta_j, \eta_j, \phi_j)$ is the joint prior of the regression parameters. We maintained the set-up of the Hamiltonian sampler in brms, using three parallel chains each with 6,000 samples after 2,000 warm-up runs. We found that $\hat{R} = 1.0$ in all runs, indicating convergent Markov chains. We report the posterior distributions of all model parameters in Supplementary Table 4.

## Estimating local trends between glacier elevation change and $V_0$ and $Q_p$

The length and thickness of the glacier dam are important diagnostics for predicting flood magnitudes and timing⁷⁶. Changes in dam geometry could therefore alter the pattern of lake outbursts, in both frequency and magnitude. In QGIS, we manually trimmed the area of all glaciers in the RGI V6.0 that impound a lake to the area below the dam. We define a dam as the area between the glacier terminus and the elevation contour on the glacier directly upstream of the lake. We then extracted the cumulative elevation change of all glacier dams between 2000 and 2019, following the procedure by ref. ⁸. In summary, this method first generates DEMs from all available stereo-pairs of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) satellite instrument intersecting glaciers during the period of 2000 to 2019. We stacked the ASTER DEMs in temporal order and further added high-resolution data from the ArcticDEM⁷⁷ to the time series. Data from ArcticDEM cover glacier dams in Alaska, Iceland and Scandinavia (Supplementary Table 6). Each pixel in the DEM time series is then interpolated using a Gaussian process regression model to estimate the elevation at any date within this period. We extracted the elevation time series for the areas of the glacier dams, resulting in an area-weighted average of 45 independent DEMs between 2000 and 2019 for each dam. Individual glacier dams are covered by at least 95% of measured pixels, drawing on at least 13 independent DEMs per pixel. We calculated the associated elevation differences between successive time steps, interpolated spatial data gaps hypsometrically⁷⁸ and finally aggregated them to the cumulative mean elevation changes with respect to 1 January 2000. On average, our glacier dams thin at a rate of 3 m yr⁻¹ during the period 2000–2019, with an average uncertainty of ±0.23 m yr⁻¹ (95% confidence level). The elevation change uncertainties

were estimated using a framework that accounts for heteroscedasticity (that is, variability in precision), which more reliably describes errors depending on terrain slope and the quality of stereo-correlation. The framework also accounts for spatial correlation of errors using multiple correlation ranges that better describe ASTER instrument noise⁷⁹. These uncertainties were extensively validated using independent high-precision data from NASA’s Ice, Cloud and land Elevation Satellite (ICESat) and IceBridge missions, as well as LiDAR acquisitions⁸. Supplementary Table 6 reports statistics on spatial coverage, the temporal density of DEM data, thinning rates between 2000 and 2019, and associated uncertainties for each glacier dam.

To assess the impact of glacier thinning on GLOF magnitudes, we selected only $J=13$ glacier lakes that burst out repeatedly ($n > 5$) between 2000 and 2019. For each dated outburst within a GLOF cycle, we obtained the closest temporal estimate of the cumulative elevation change of the glacier dam. We then fitted a hierarchical quantile regression model, that is, the same as in equations (1)–(13), to estimate the median of reported $V_0$ and $Q_p$ at a given glacier dam from the associated cumulative elevation change $h$. As we did when calculating temporal trends of $V_0$, $Q_p$ and $Z$, we chose an AL likelihood to estimate the conditional median of $V_0$ and $Q_p$, with the same (hyper-)prior distributions as in equations (7)–(12). We ran three parallel chains with 6,000 samples and 2,000 warm-up runs and found no numerical divergences. We report the posterior distributions of all model parameters in Supplementary Table 5.

## Data availability

The GLOF inventory V2.0 is freely available at http://glofs.geology.uni-potsdam.de, and archived at https://doi.org/10.5281/zenodo.7326571. The ALOS DEM V3.2 is freely available at https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm. The RGI V6.0 is freely available at https://nsidc.org/data/nsidc-0770. Planet and RapidEye images are available under a basic license from the Education and Research programme from Planet Labs Inc. through the Planet-Explorer (https://www.planet.com/). Landsat images were obtained from the EarthExplorer (https://earthexplorer.usgs.gov/). Data on glacier elevation change for the period 2000–2019 are freely available at https://doi.org/10.6096/13.

## Code availability

All codes to reproduce the statistical analysis are available at https://github.com/geveh/IceDamFailures, and are archived at https://doi.org/10.5281/zenodo.7326865.

59. Iturrizaga, L. in Encyclopedia of Snow, Ice and Glaciers (eds Singh, V. P. et al.) 381–399 (Springer, 2011); https://doi.org/10.1007/978-90-481-2642-2_196.

60. Rickman, R. L. & Rosenkrans, D. S. Hydrologic Conditions and Hazards in the Kennicott River Basin, Wrangell-St. Elias National Park Preserve, Alaska Water-Resources Investigations Report 96-4296 (USGS, 1997); https://doi.org/10.3133/wri964296.

61. Bhambri, R. et al. Ice-dams, outburst floods, and movement heterogeneity of glaciers, Karakoram. *Glob. Planet. Change* **180**, 100–116 (2019).

62. Stone, K. H. Alaskan ice-dammed lakes. *Ann. Assoc. Am. Geogr.* **53**, 332–349 (1963).

63. Pfeffer, W. T. et al. The Randolph Glacier Inventory: a globally complete inventory of glaciers. *J. Glaciol.* **60**, 537–552 (2014).

64. Takaku, J., Tadono, T., Doutsu, M., Ohgushi, F. & Kai, H. Updates of ‘AD3D30’ ALOS global digital surface model with other open access datasets. *Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci.* **XLIII-B4-2020**, 183–189 (2020).

65. Wulder, M. A. et al. The global Landsat archive: status, consolidation, and direction. *Remote Sens. Environ.* **185**, 271–283 (2016).

66. Pekel, J.-F., Cottam, A., Gorelick, N. & Belward, A. S. High-resolution mapping of global surface water and its long-term changes. *Nature* **540**, 418–422 (2016).

67. McElreath, R. Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman and Hall/CRC, 2020).

68. Gudmundsson, L. et al. Globally observed trends in mean and extreme river flow attributed to climate change. *Science* **371**, 1159–1162 (2021).

69. Blöschl, G. et al. Changing climate both increases and decreases European river floods. *Nature* **573**, 108–111 (2019).

70. Rokaya, P., Budhathoki, S. & Lindenschmidt, K.-E. Trends in the timing and magnitude of ice-jam floods in Canada. *Sci. Rep.* **8**, 5834 (2018).
71. Seneviratne, S. I. et al. in *IPCC Climate Change 2021: The Physical Science Basis* (eds Masson-Delmotte, V. et al.) 1513–1766 (Cambridge Univ. Press, 2021).

72. Ng, F. & Liu, S. Temporal dynamics of a jökuthlaup system. *J. Glaciol.* **55**, 651–665 (2009).

73. Stan Development Team *Stan Modeling Language Users Guide And Reference Manual*, version 2.26.1; https://mc-stan.org (2022).

74. Bürkner, P.-C. brms: an R package for Bayesian multilevel models using Stan. *J. Stat. Softw.* **80**, 1–28 (2017).

75. R Core Team. R: a language and environment for statistical computing. *R Foundation for Statistical Computing* https://www.R-project.org/ (2022).

76. Clarke, G. K. C. Hydraulics of subglacial outburst floods: new insights from the Spring-Hutter formulation. *J. Glaciol.* **49**, 299–313 (2003).

77. Porter, C. et al. ArcticDEM. https://doi.org/10.7910/DVN/OHHUKH (Harvard Dataverse, 2018).

78. McNabb, R., Nuth, C., Kääb, A. & Girod, L. Sensitivity of glacier volume change estimation to DEM void interpolation. *Cryosphere* **13**, 895–910 (2019).

79. Hugonnet, R. et al. Uncertainty analysis of digital elevation models by spatial inference from stable terrain. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **15**, 6456–6472 (2022).

80. Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. *Sci. Data* **7**, 109 (2020).

**Acknowledgements** We thank Planet Labs Inc. for a basic license in the Education and Research programme that granted free access to high-resolution Planet and RapidEye

images. We thank many experts for sharing regional data on ice-dammed lakes: B. Einarsson, J. Carrivick and F. Tweed (Iceland); M. Huss and M. Werder (European Alps); M. Jackson (Scandinavia); P. Iribarren Anacona (Andes); R. Johnston (northwest North America); J. Steiner and V. Round (High Mountain Asia); and V. Foord and J. Jacklin for internal reviews.

**Author contributions** G.V. designed and led the study, analysed the data and produced all figures. N.L. and J.T. helped collect the data. R.H. and O.K. contributed to the data analysis. K.V., O.K. and L.V.L. contributed to formulating the Bayesian models. J.J.C., O.K., M.G., R.H. and L.V.L. contributed to the discussion of the results and to the writing of the manuscript.

**Funding** Open access funding provided by Universität Potsdam.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-022-05642-9.

**Correspondence and requests for materials** should be addressed to Georg Veh.

**Peer review information** *Nature* thanks Guoxiong Zheng and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
a

<figure><img src="image_12.png" /></figure>

b

<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 1** | GLOF reporting between 1900 and 2021. a, Histograms of all reported GLOFs (top), reported GLOFs with documented peak discharge $Q_p$ (middle), and reported GLOFs with documented flood volume $V_0$. Numbers are the total number of GLOFs every three decades. b, Plot of year of GLOF

c

<figure><img src="image_14.png" /></figure>

occurrence versus year of first available document. White circles are individual GLOFs; colours are aggregated cases in 5-year intervals. c, Histogram showing the year of first available documentation for all cases in our database.
a

<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 2** | Regional trends of GLOF peak discharge $Q_p$ with time obtained by hierarchical quantile regression. a, Temporal trends of $Q_p$ for the periods 1900–2021 and 1990–2021. Colours distinguish between the median (blue) and 90th (orange) percentile of reported $Q_p$, including the posterior 95% HDI (shades). Numbers (n) and percentages in lower left of panels are, respectively, the sample sizes and shares of GLOFs used in the quantile

b

<figure><img src="image_16.png" /></figure>

regression model for each region. b, Posterior trends in $Q_p$ for the 50th and 90th percentiles. Posterior 95% HDI (lines) and median (dot) regression trend of $Q_p$ are grouped by the periods 1900–2021 (dark green) and 1990–2021 (light green). Distributions overlapping with zero (thick black line) indicate the lack of a credible trend in $Q_p$. Percentages show posterior mass smaller (left) and larger (right) than zero.
a

<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 3** | Regional trends of GLOF volume $V_0$ with time $t$ obtained by hierarchical quantile regression. **a**, Temporal trends of $V_0$ for the periods 1900–2021 and 1990–2021. Colours distinguish between the median (blue) and 90th (orange) percentile of reported $V_0$, including the posterior 95% HDI (shades). Numbers ($n$) and percentages in lower left of panels are, respectively, the sample sizes and shares of GLOFs used in the quantile

b

<figure><img src="image_18.png" /></figure>

regression model for each region. **b**, Posterior trends in $V_0$ for the 50th and 90th percentiles. Posterior 95% HDI (lines) and median (dot) regression trend of $V_0$, grouped by the period 1900–2021 (dark green) and 1990–2021 (light green). Distributions overlapping with zero (thick black line) indicate the lack of a credible trend in $V_0$. Percentages show posterior mass smaller (left) and larger (right) than zero.
<figure><img src="image_19.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_24.png" /></figure>

**Extended Data Fig. 4** | Global trends of GLOF volume $V_0$ and peak discharge $Q_p$ with time obtained by hierarchical quantile regression. a, Pooled global trends in median and extreme $V_0$ for the period 1900–2021. b, Pooled global trends in median and extreme $V_0$ for the period 1990–2021. c, Posterior pooled trends in $V_0$ for the 50th and 90th percentiles for both periods. d, Pooled global trends in median and extreme $Q_p$ for the period 1900–2021. e, Pooled global trends in median and extreme $Q_p$ for the period 1990–2021. f, Posterior pooled trends in $Q_p$ for the 50th and 90th percentiles for both periods. In a, b, c, and

e, colours distinguish between the median (blue) and 90th (orange) percentile of reported $V_0$ and $Q_p$, including the posterior 95% HDI (shades). In c and f, lines are the posterior 95% HDI and dots are the median of the global mean regression trends, grouped by the period 1900–2021 (dark green) and 1990–2021 (light green). Distributions overlapping with zero (dashed black line) indicate the lack of a credible trend in $V_0$ and $Q_p$. In c and f, $\pi_{post} < 0$ denotes a posterior mass smaller than zero.
a

<figure><img src="image_25.png" /></figure>

b

<figure><img src="image_26.png" /></figure>

**Extended Data Fig. 5** | Posterior trends in local median (50th percentile) peak discharge $Q_p$ and flood volume $V_0$ obtained by hierarchical quantile regression. Posterior 95% HDI (lines) and median (dot) regression trend of a, $V_0$ and b, $Q_p$ of ice-dammed lakes. Glacier names are those of the parent glacier (RGI ID). Lakes labelled 'NA' have no name. Distributions overlapping with zero

(dashed black line) indicate lack of credible trend in $V_0$ and $Q_p$. Percentages are posterior masses smaller (left) and larger (right) than zero. Trends refer to log₁₀-transformed values of $V_0$ and $Q_p$. The posterior densities show decadal trends in orders of magnitudes.
<figure><img src="image_27.png" /></figure>

**Extended Data Fig. 6** | Regional seasonality of the monthly number of GLOFs and mean monthly air temperatures. Mean monthly air temperatures at the locations of all lakes that caused a GLOF during the period 1901–2021

(data from CRU TS V4.05, ref. ⁸⁰) were used to calculate the regional monthly average during that period. Months are sorted from 1 (January) to 12 (December).
<figure><img src="image_28.png" /></figure>

**Extended Data Fig. 7 | Posterior differences in GLOF timing (doy) for 40 ice-dammed lakes obtained by Bayesian regression.** Lines and dots show, respectively, the posterior 95% HDI and median of the change in GLOF timing,

expressed in days, between the first and the last reported GLOF from a given lake. Distributions overlapping with zero (thick vertical line) suggest no credible change in the day of year (doy). Lakes labelled 'NA' have no name.
<figure><img src="image_29.png" /><figcaption>**Extended Data Fig. 8** | Regional elevation distribution of ice-dammed lakes. Histograms in top row show all lakes (>0.01 km²) that currently exist in the six regions. Those in the bottom row are all lakes that burst out between</figcaption></figure>

<caption>1900 and 2021. Vertical orange lines are the median elevation of ice-dammed lakes in each region.</caption>
<figure><img src="image_30.png" /></figure>

**Extended Data Fig. 9** | Mean elevation change of glacier dams with reported GLOF magnitudes between 2000 and 2019. Black line and blue shading are, respectively, the mean and standard error of glacier elevation change obtained

by Gaussian Process regression⁸. Black ticks are dated GLOFs. Note that Tulsequah Glacier (RGI60-01.01521) dams two lakes. The reversed trend for Kyagar Glacier (RGI60-14.07098) is due to a glacier surge.