Causes, consequences and implications of the 2023 landslide-induced Lake Rasac glacial lake outburst flood (GLOF), Cordillera Huayhuash, Peru

Glacierised Peruvian mountain ranges are experiencing accelerated, climate-change-driven glacier ice loss. Peru's second highest mountain range, the Cordillera Huayhuash, has lost about 40 % (∼ 34 km²) of its glacier cover since the 1970s. Newly exposed landscapes are prone to a number of hazard processes including the formation and evolution of glacial lakes, changing stability conditions of mountain slopes, and rapid mass movements. In this study, we integrate the analysis of meteorological data, remotely sensed images and field observations in order to document the most recent (February 2023) large mass-movement-induced glacial lake outburst flood (GLOF) from the moraine-dammed Lake Rasac. The GLOF was triggered by a mass movement from the failure of a part of an arête ridge (i.e. narrow ridge separating parallel glacial valleys) with an estimated volume of 1.1×10⁶ m³–1.5×10⁶ m³. This occurred in a rock zone where climate information – primarily from reanalysis data – indicates cold, deep permafrost and was preceded by several small-magnitude precursory rockfall events. The reduced stability of the frozen rocks in the detachment zone is most likely related to deep warming but not to especially critical conditions of warm permafrost with higher amounts of unfrozen water. Further, we describe the surprisingly short-distance process chain (attenuated by Lake Gochacotan, located 3.5 km downstream from the detachment zone) and analyse the transport of large boulders with the use of hydrodynamic modelling, revealing that flow velocities >5 m s⁻¹ must have been reached in the case of translational motion and >10 m s⁻¹ in the case of rotational motion of the largest transported boulders (diameter >3.5 m). In addition, we analyse climate trends over the past 8 decades and meteorological conditions in the days prior to the GLOF, revealing a statistically significant temperature increase trend over the past 8 decades, as well as an atmospheric temperature rise and positive thermal anomaly in different soil layers before the event. We argue that climate change effects (warming air and permafrost temperatures) served to hasten the failure. This study helps us to understand (i) mechanisms, amplification and attenuation elements in GLOF process chains and (ii) frequency–magnitude relationships of extreme geomorphic processes that undergo alteration due to the rapidly changing high-mountain environments on a regional scale (both large-magnitude rockfalls and GLOFs). This study supports earlier work that indicated an increasing frequency of large mass-movement-induced GLOFs originating from the warming cryosphere in recent decades.

1 Introduction

The world's high-mountain regions are losing glacier ice and permafrost at an accelerated pace in response to climate change (Hugonnet et al., 2021), with important implications for downstream areas (Clague et al., 2012; Knight and Harrison, 2013; Huss et al., 2017; Immerzeel et al., 2020). There are growing concerns about potentially disastrous, far-reaching process chains originating from a warming cryosphere (Haeberli et al., 2017; Ding et al., 2021), as exemplified by recent disasters such as the 2020 Lake Jinwuco outburst (Zheng et al., 2021), the 2021 rock and ice avalanche in Chamoli in India (Shugar et al., 2021), the 2022 Marmolada glacier collapse in Italy (Bondesan and Francese, 2023; Olivieri and Bettanini, 2023), and the 2023 South Lhonak Lake cross-border outburst impacting vast areas along the Teesta River in India and Bangladesh (SANDRP, 2023). Among the hazard processes of greatest concern for mountain communities and foothill areas are glacial lake outburst floods (GLOFs) – sudden releases of water from glacial lakes (Costa and Schuster, 1988; Korup and Tweed, 2007; O'Connor et al., 2013; Emmer et al., 2021; Zhang et al., 2024).

The Peruvian Andes have witnessed >100 GLOFs in the past century, some of which have had catastrophic consequences (Carey, 2005; Emmer et al., 2022b). Strikingly, the two cordilleras located in northern Peru (Cordillera Blanca and Cordillera Huayhuash) exhibit a higher concentration of major GLOFs (Kinzl et al., 1954; Zapata, 2002), including the 1941 Lake Palcacocha disaster that destroyed part of Huaraz (Huggel et al., 2020; Wegner, 2024). The most recent catastrophic Peruvian GLOF was the 2020 GLOF triggered by an ice–rock avalanche at the moraine-dammed Lake Salkantaycocha in Cordillera Vilcabamba in 2020 (Vilca et al., 2021). Emmer et al. (2020) showed that while major GLOFs induced by moraine dam failure in the late 1930s to early 1950s were predominantly triggered by ice avalanches and calving processes, more recent events such as the 2020 Salkantaycocha GLOF were more frequently triggered by landslides and rockfalls. Despite recent advances in compiling GLOF inventories (Emmer et al., 2022a; Veh et al., 2022; Lützow et al., 2023), forecasting and predicting GLOF occurrence in space and time is challenging due to (i) a wide range of possible combinations of lake and dam characteristics favouring GLOF occurrence, triggers and mechanisms (Emmer and Cochachin, 2013); (ii) the often non-periodic, one-off nature of GLOFs (especially GLOFs originating from bedrock- and moraine-dammed lakes; Lützow et al., 2023); and (iii) dynamic conditions, including the evolution of proglacial lakes driven by glacier retreat (Shugar et al., 2020), and changing slope stability conditions driven by glacier recession and permafrost warming (Stoffel and Huggel, 2012; Haeberli et al., 2017). Therefore, the frequency–magnitude relationship and derived recurrence interval approaches traditionally used in flood hazard assessments on a catchment scale cannot be applied to GLOFs; the estimates of GLOF probabilities are instead derived from detailed understanding of past events at a regional scale, with consideration for how future conditions will evolve (Haeberli et al., 2022; Emmer, 2024).

This study aims to examine the preconditions, triggers and impacts of the 2023 Lake Rasac GLOF triggered by large ice–rock failure (Cordillera Huayhuash) and to discuss whether and, if so, to what extent this event can be attributed to climate change. The importance of documenting such events lies in better evidence-driven understanding of (i) mechanisms as well as amplification and attenuation elements in GLOF process chains and (ii) altering frequency–magnitude relationships of extreme processes in rapidly changing high-mountain environments on a regional scale (both high-magnitude rockfalls and GLOFs). Such outcomes can serve the scientific community (in particular geomorphologists, natural hazard scientists and process chain model developers) as well as disaster risk reduction authorities and practitioners.

2 Study area and previous lake and GLOF research

2.1 The Cordillera Huayhuash

The Cordillera Huayhuash (10.0–10.5° S, 76.8–77.0° W; Fig. 1), located 200 km north of Lima, is the second highest mountain range in Peru (Nevado Yerupajá, 6617 m a.s.l.), covering an area of nearly 1200 km² (INAIGEM, 2018). The geological setting of Cordillera Huayhuash is complex and comprises faulted and folded Cretaceous sedimentary carbonates of marine origin, Tertiary silicic volcanic rocks, and intrusive rocks (tonalite, granite, granodiorite, diorite), in places covered by Quaternary deposits of glacial and glaciofluvial origin, alluvium and colluvium. The main fault zone is oriented in the north–south direction and is parallel to the main massif of the Cordillera Huayhuash (IGM, 1975; Hall et al., 2009). Situated in the outer tropics, the climate of the Cordillera Huayhuash is characterised by relatively consistent temperatures year-round, with minimal seasonal variation. Air temperatures are about 1–2 °C warmer during the austral summer compared to the austral winter (Rabatel et al., 2013). The Cordillera Huayhuash exhibits a pronounced east–west precipitation gradient, with considerably wetter conditions on the eastern side, and approximately 75 % of the annual precipitation occurs between November and April (Bookhagen and Strecker, 2008; Hall et al., 2009)

Figure 1 Location of the Cordillera Huayhuash in Peru (a); lakes (Bat'ka et al., 2020), glaciers (RGI Consortium, 2017) and topography (ASF, 2024) of the Cordillera Huayhuash (b) and detail of the upper Jahuacocha and Rasac valleys (c).

The high-elevation areas of the Cordillera Huayhuash (above approximately 4500–5000 m a.s.l.) are glacierised and exhibit several generations of moraines from past glaciations (Hall et al., 2009). The last significant glacier advance occurred during the Little Ice Age, which subsided in the second half of the 19th century. According to INAIGEM (2018), the total glacier extent decreased from 86.89 km² in 1975 to 53.06 km² in 2016. The estimated average snow-line altitude (an approximation of equilibrium altitude) was 5051 m a.s.l. in 1986–2005 (McFadden et al., 2011), and this elevation has most likely increased since then. Repeated Holocene advance–retreat cycles of glaciers preconditioned the formation of different types of glacial lakes. Using Landsat imagery with a resolution of 30 m, Wood et al. (2021) mapped 129 lakes with a total area of 6.92 km², within a 1 km buffer distance from the glaciers, while Bat'ka et al. (2020) used high-resolution (sub-metre) images from the Google Earth collection to map 270 lakes with a total area of 6.32 km² in the whole Cordillera Huayhuash. This implies that accounting for small lakes influences the total number of mapped lakes while large lakes control total lake area. One of these lakes – Rasac ($\mathrm{10}\mathit{°}{\mathrm{15}}^{\prime }{\mathrm{50}}^{\prime \prime }$ S, $\mathrm{76}\mathit{°}{\mathrm{56}}^{\prime }{\mathrm{18}}^{\prime \prime }$ W; 4654 m a.s.l.) – and the GLOF it produced in 2023 are the subject of this study (see Sect. 2.2).

2.2 Lake Rasac

Lake Rasac (Figs. 1 and 2) was a moraine-dammed water body located at 4654 m a.s.l. behind the moraine deposited by Rasac glacier in the north-facing Rasac valley, a left bank tributary of the Jahuacocha/Achin stream (Salazar and Valverde, 2022). Considering the elevation of this moraine and its position in the valley system of moraine assemblages, it is assumed they were last reached during the Little Ice Age advance. While the lake did not have surface outflow (dam freeboard was 10–15 m, depending on lake water level), it was drained by a system of subsurface channels through the moraine dam and the water level (area and volume) fluctuated both inter- and intra-seasonally. The area of the lake varied between 53 000 and 64 000 m² (ANA, 2014; INAIGEM, 2018) and the depth was not known. Using the empirical equation developed by Muñoz et al. (2020), the estimated mean lake depth is 5 m and so the estimated lake volume is approximately 300 000 m³. Another lake – Gochacotan ($\mathrm{10}\mathit{°}{\mathrm{14}}^{\prime }{\mathrm{58}}^{\prime \prime }$ S, $\mathrm{76}\mathit{°}{\mathrm{56}}^{\prime }{\mathrm{25}}^{\prime \prime }$ W; 4285 m a.s.l.) – is located approx. 1.5 km downstream from Lake Rasac (see Fig. 2).

Figure 2 Lake Rasac and its surroundings in July 2019. (a) Upstream view towards the head of the Rasac valley (looking south); (b) downstream view towards both lakes Rasac and Gochacotan (looking north) from near Rasac peak (5129 m a.s.l.); (c) a view towards the Rasac arête ridge (looking east), with the 2023 release zone highlighted. Field images: Adam Emmer.

Bat'ka et al. (2020) conducted a comprehensive assessment of mountain range-wide GLOF susceptibility of 46 lakes larger than 20 000 m². They found that for Lake Rasac, the most probable GLOF scenario involved dam overtopping as a result of a landslide into the lake and ranked this lake as the 16th most susceptible from all assessed lakes in the Cordillera Huayhuash (other moraine-dammed lakes reached higher scores because of their surface outflows; see the assessment methodology of Emmer and Vilímek, 2014).

3 Data and methods

3.1 Remotely sensed data and analysis

We analysed satellite images and derived products (such as digital elevation models, DEM; see Table 1) from various sources for event description and geomorphic interpretation. These data were used to (i) characterise mapped features from high-resolution optical images (e.g. area of the lake, areas of landslide release and GLOF impact zones, location of largest transported boulders), (ii) derive topographical and morphometric characteristics (e.g. a slope of the longitudinal profile along the process chain), and (iii) narrow down the timing of the GLOF from images with fine temporal resolution (Planet and Sentinel images). Further, we cross-checked the interpretation of remotely sensed images against field observations (field visit in July 2019; post-event field images and testimony of an eyewitness – a member of the community settled in the Rio Jahuacocha valley a few kilometres downstream). All the analyses of remotely sensed images were undertaken in the open-source QGIS software (https://qgis.org/, last access: 18 March 2025).

Table 1 Remotely sensed data used in this study.

n/a: not applicable.

3.2 Climate data and analysis

3.2.1 Climate trends

Given that meteorological station records in Peru are very sparse and limited (none of which have consistent records going to 2023), we employed the state-of-the-art ERA5 gridded reanalysis products (∼ 31 km horizontal resolution) with records starting in 1940 through the present, from the European Centre for Medium-Range Weather Forecasts (Hersbach et al., 2020). The list of variables analysed is detailed in Table 2. For the climate analysis, we used the average value over a 1° × 1° box (−77.5°, −76.5°, −10.75°, −9.75°) covering Rasac, representing the meteorological conditions of the lake, to investigate (1) the long-term trend over the past 8 decades 1940–2023 and (2) the anomalous states (compared with the climatological mean averaged over 1981–2010) that occurred in the near-surface atmosphere and in the land surface and subsurface of different soil layers for the time period prior (up to 2 years) to the event. Monthly products were used to analyse the climate trends and daily products were used to analyse the anomalous conditions prior to the GLOF event. The atmospheric circulation variables zonal (U) wind, meridional (V) wind and geopotential height at different vertical levels (850, 500 and 250 hPa, representing the lower, middle and upper troposphere) were also extracted from ERA5 covering the region (−85.0, −65.0, −20.0, and 0.0°) to investigate any possible anomalous circulation patterns associated with the event but did not reveal any notable changes prior to the event and are hence not shown here. It should be noted that the variables from ERA5 are not directly assimilated but are generated by the atmospheric components of the Integrated Forecast System (IFS) modelling framework.

Table 2 Climate variables analysed from ERA5 products in this study.

3.2.2 Climate change attribution

We use the meteorological variables averaged over the 1° × 1° box $(-\mathrm{77.5}\mathit{°},-\mathrm{76.5}\mathit{°},-\mathrm{10.75}\mathit{°},-\mathrm{9.75}\mathit{°})$ covering Rasac, taking the February 1 d maximum as the event definition. As outlined in Philip et al. (2020), we assess how the likelihood of occurrence of the event has changed due to anthropogenic climate change. We use a generalised extreme value (GEV) distribution that scales with global mean surface temperature – GMST (i.e. using GMST as covariate, Philip et al., 2020, and Van Oldenborgh et al., 2021). The resulting probability ratio for the event is calculated for the comparison between the climate of 2023 and the pre-industrial climate – 1.2 °C cooler than the current climate. The probability ratio is calculated as the probability of occurrence under the current climate divided by the probability of occurrence under the pre-industrial climate.

We acknowledge that in typical statistical event attribution studies, climate model outputs are also analysed (often taken from the most up-to-date Coupled Model Intercomparison Project global Earth system models or regional model experiments) through model validation or multi-method attribution (combining observed and model results) to arrive at a synthesised attribution statement. In this study, we did not perform analyses based on climate model outputs for two reasons: (1) the Rasac 2023 GLOF was a localised event with no identifiable impact on humans and the built environment and was not directly triggered by meteorological variables but triggered by mass movements (as detailed below) and (2) the attribution of temperature anomalies is mainly to demonstrate the anthropogenic climate change impact on local warming associated with the GLOF rather than establishing a causal chain directly linking meteorological anomaly to the GLOF occurrence.

3.3 Hydrodynamic analysis

If the slope conditions and the composition of the terrain are known and there is clear evidence of deposited boulders and cobbles, the minimum water velocity that moved them can be estimated. Various empirical relationships exist for calculating the critical flow velocity and shear stress required to initiate particle motion (Van Rijn, 2019). However, these formulae are predominantly derived for fine spherical sediments or stones and under boundary conditions that do not align with the impulse waves generated by landslides. To calculate boulder transport, we employ a method considering the balance of forces and torque acting on boulders modelled as simplified rectangular shapes (Bressan et al., 2018). Since there is no available survey detailing the size and shape of the boulders in the impact area for this analysis, we assume that the boulder's dimensions are characterised by a height half the size of their width and length, which are considered equal. This assumption aims to approximate the realistic shape of boulders where, typically, the boulder height does not match the width and length. Like other sediment transport phenomena, the forces acting on the boulder include the following.

$$\begin{array}{}\text{(1)}& {\displaystyle }\text{Gravity force,}G=gV({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{w}}{V}_{\mathrm{w}}/V).\text{(2)}& {\displaystyle }\text{Drag force,}{F}_{\mathrm{D}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{f}}{v}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{d}}.\text{(3)}& {\displaystyle }\text{Lift force,}{F}_{\mathrm{L}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{b}}{v}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{l}}.\text{(4)}& {\displaystyle }\text{Friction force,}{F}_{\mathrm{S}}=\mathit{\mu }(G-L).\end{array}$$

Here, g stands for acceleration due to gravity (9.81 m s⁻²), V is the volume of the boulder (m³), V_w is the submerged volume based on the water depth h (m), ρ_s is boulder density (for granodiorite we assumed 2700 kg m⁻³), ρ_w is water density (1000 m³ s⁻¹), A_f represents the front wetted area of the boulder (m²) based on the water depth h (m), A_b represents the bottom area of the boulder (m²) and v is the depth-averaged approaching speed velocity (m s⁻¹). The drag coefficient C_d is assumed to be 2.05 for a partially submerged boulder and 1.05 for a fully submerged boulder of cubic shape (Bressan et al., 2018). C_l is the lift coefficient with a constant value of 0.178, and μ is the friction coefficient with a constant value of 0.5 (Nott, 2003; Lorang, 2011; Bressan et al., 2018). The current literature does not provide comprehensive guidelines for boulder transport. Friction, drag and lift coefficients, as used in our study, are sourced from the general literature (Çengel and Cimbala, 2013). Accurate determination of these coefficients would require specific experimental or mathematical models, which are beyond the scope this study.

By applying a force balance in the direction of flow and a torque balance around the rear edge of the boulder, the required velocity v (m s⁻¹) to disrupt this balance and initiate the incipient motion of the boulder can be calculated. It is assumed that if the velocity disrupts the force balance, translational motion begins, while if the velocity disrupts the torque balance, rotational motion is initiated (Nandasena and Tanaka, 2013; Bressan et al., 2018). It is evident from the above discussion that certain parameters are contingent on water depth. Due to the absence of data, a simplified relationship between wave depth and averaged velocity must be assumed.

3.4 Assessment of permafrost conditions

The assessment of permafrost conditions is based on the analysis of Google Earth images, (sparse) climatic information, experience of borehole temperatures available in Europe, and scientific literature about the slope stability of ice-clad and perennially frozen mountain peaks. They constitute “best estimates” for a site without direct measurements, dealing with thermal conditions, subsurface and surface ice properties, and the resulting hydro-mechanical aspects in view of climate-related stability changes. Large-scale (global) model results relating to patterns of permafrost occurrence are available from Gruber (2012) based on climate data for the second half of the 20th century (1960–1990) and from Obu (2021) based on climate data (freezing–thawing indices) for the beginning of the 21st century (2000–2016) in combination with empirical thermal offsets mainly related to effects from vegetation and winter snow cover; comparison with measured borehole temperatures indicate an uncertainty range of ±2 °C. Due to the low spatial resolution (about 1 km) of these two important approaches, topo- and microclimatic effects at the local Rasac detachment site on an extremely steep slope must be assessed using additional process-related indications as described further below.

4 Results: description and interpretation of the 2023 event

4.1 Surface and subsurface ice, prior rockfall activity, and the 2023 arête ridge failure

4.1.1 Surface and subsurface ice in relation to slope stability

The NW–SE-oriented Rasac arête ridge reaches elevations exceeding 6000 m a.s.l. and is partly glacier-covered. Surface ice on very steep slopes (generally > 45°) occurs in two main forms: ice aprons (ice faces; for additional info, see Ravanel et al., 2023) and hanging glaciers. Both forms, together with ice crests and small summit ice caps, are abundant in the region. Ice aprons or ice faces are thin (a few metres), smooth, near-static, cold ice covers frozen onto steep perennially frozen rock faces. Like perennial ice patches in less inclined topography, such ice covers cool the ground because they have a high albedo and because they buffer against the penetration of above-zero temperatures whilst allowing penetration of cold waves to the ground below. The ice itself may be thousands of years old, and its disappearance due to atmospheric temperature rise induces freeze–thaw cycles within the uncovered frozen rocks, followed by minor-to-major rock destabilisation (Ravanel et al., 2023). Hanging glaciers are much thicker and move from a less inclined accumulation area towards the front, where ablation tends to be by dry calving. Even with strongly negative mean annual air temperatures, they can be polythermal in that meltwater percolating from the surface and refreezing at depth of the firn constitutes a strong source of latent heat, often inducing temperate firn/ice with high water pressure (Margreth et al., 2017). Two important consequences of such conditions are that (i) the stability of such glaciers in steep terrain is mainly related to the cold, near-vertical ice front frozen to bedrock in a narrow zone at the lower margin and (ii) the temperate firn area with percolating meltwater can introduce strong local thermal anomalies (Haeberli et al., 1997) and complex thermo-hydraulic conditions in the underlying permanently frozen rocks (Haeberli et al., 2004). Macroscopic subsurface ice in perennially frozen bedrock usually fills cracks and fissures. As the original formation of such ice involves efficient frost cracking and frost wedging, such ice fillings can be massive, as is sometimes visible in permafrost tunnels (Haeberli et al., 1979; Mamot et al., 2021) or uncovered immediately after detachments of rock avalanches (Gruber and Haeberli, 2007).

One of the areas that developed frequent rockfall activity is located upstream from Lake Rasac at an elevation between 5500 and 5750 m a.s.l. (see Fig. 2c). The detachment zone of the 2023 event was covered by ice aprons rather than hanging glaciers. As a consequence, rather homogenous spatial patterns of subsurface temperatures can be expected. The existence of ice aprons on the WSW-oriented detachment side of the Rasac arête ridge and the lack of any ice cover on its other, ENE-oriented, side is likely to enhance the temperature differences between the two sides and the related horizontal flux of heat across the ridge towards the detachment site.

4.1.2 Prior rockfall activity and the 2023 GLOF-triggering failure

The Rasac arête ridge is formed of tilted layers of sedimentary and metamorphic rocks inclined towards the Rasac valley (i.e. parallel with the slope inclination; see Fig. 3a) that favour disintegration and that precondition mass movements. Available high-resolution satellite images with sub-metric resolution going back to 2010 (Fig. 3b–e) reveal that the release zone gradually migrated to higher elevations, continuously undermining an overhanging block that eventually failed in February 2023. Small rockfalls from this release zone have occurred frequently since 2018, i.e. 5 years before the major failure, as documented by the Sentinel time series 2017–2023 (Fig. 4; see also Fig. 2c) and the analysis of Planet images.

Figure 3 Release zone of the GLOF-triggering landslide. Panel (a) shows a field image of the hanging block as seen from the opposite slope in July 2019 (image: Adam Emmer); panels (b) to (e) show pseudo-3D views of the release zone between 2010 and 2023 in © Google Earth (Maxar images, CNES/Airbus images). Note the very little rockfall activity (change of the release zone) between 2010 and 2014.

Figure 4 Sentinel image time series depicting rockfall activity and development of the release zone between 2017 and 2023. Note gradually developing head scarp and release zone (R) and increasing amount of material deposited on the glacier and slopes above the eastern bank of the lake (D).

According to a member of the community settled in the Jahuacocha valley, the GLOF-triggering failure of an overhanging section of a very steep (>60°) part of the Rasac arête ridge occurred around noon on 12 February 2023. The head scarp was located at an elevation between 5750 and 5800 m a.s.l.; the approximated planar extent of the release zone is 250 m wide and 150 m high on average (Fig. 3a). Field images taken in 2019 reveal that the thickness of the hanging block was approximately 30 to 40 m, resulting in a volume estimation of 1.1×10⁶ m³–1.5×10⁶ m³. This means that the landslide would be considered “large” under the landslide classification scheme of McColl and Cook (2024). The 2023 landslide material was rock mixed with ice. In view of the thin (a few metres), co-detached ice aprons, the relative amount of ice in the 2023 rock–ice avalanche was probably on the order of 10 % to 20 %. The vertical drop between the head scarp and Lake Rasac is 1200 m (Fig. 5).

Figure 5 Topographical profile of the whole process chain (a) and schematic sketch of the 2023 Rasac arête ridge failure (b). Note the x:y ratio of 2:1. DEM: 12.5 m ALOS PALSAR (ASF, 2024).

Large slope failures are often triggered by earthquakes or extreme hydro-meteorological conditions. According to the USGS earthquake catalogue (USGS, 2024), 71 earthquakes with magnitude > 4.5 occurred within the 200 km buffer distance from the landslide site between 2015 and 2024, with the highest magnitude reaching M 5.6 (95.8 km deep near Lima in June 2022). However, if the buffer is decreased to 100 km, the search returns only five earthquakes (highest magnitude M 4.9). No earthquakes occurred within the 200 km buffer distance from the landslide site near the time of its occurrence, and so we exclude the possibility of direct earthquake triggering as well as preconditioning. Instead, we examine climate change preconditions and evaluate the role of subsurface ice and permafrost warming (see Sect. 4.2).

4.2 Climate change and thermal and permafrost conditions

4.2.1 Changing climate in the study area between 1940 and 2023

From a longer-term climate change perspective over several decades, the study region has been warming significantly, receiving less snowfall and experiencing a rise in the 0 °C elevation. Over the studied region, all temperature metrics showed a substantial increase over the 84 years prior to 2023 (Fig. 6 and Figs. S1–S9 in the Supplement). The linear regression of the 2 m air temperature trend revealed an increase of +1.97 °C for February temperatures and +2.06 °C for January temperatures between 1940 and 2023. February 2023 was warmer compared to the historical mean but does not stand out as anomalously warm compared to other individual years (Fig. 6a), and it is also worth noting that in the previous 2 years, the temperature was also slightly warmer than the long-term mean. Although there is no significant trend in total rainfall (Fig. 6d), snowfall in February and January has shown a consistent decrease, with February experiencing a stronger decrease (Fig. S1). As a result, the ratio of snowfall to rainfall has also been decreasing over the study region in January and February. Interestingly, but not surprisingly, the 0 °C elevation has also been rising significantly in January and February, with February experiencing a stronger increase (Fig. 6e and f). The linear regression of the 0 °C elevation trend between 1940 and 2023 revealed an increase of +258 m for February and +248 m for January.

Figure 6 Selected climate variables, averaged over the 1° × 1° box $(-\mathrm{77.5}\mathit{°},-\mathrm{76.5}\mathit{°},-\mathrm{10.75}\mathit{°},-\mathrm{9.75}\mathit{°})$ covering Rasac, showing their seasonal cycles and historical trends. Seasonal cycle of 2 m air temperature (°C) for 1940–2023 shown in panel (a), and the historical trend for February and January shown in panel (b). Seasonal cycle of total precipitation (rainfall plus snowfall, in mm d⁻¹) for 1940–2023 shown in panel (c), and the historical trend for February and January shown in panel (d). Seasonal cycle of 0 °C level (m) for 1950–2023 shown in panel (e), and the historical trend for February and January shown in panel (f). The 2 prior years – 2021 and 2022 – are also highlighted in panels (a), (c) and (e).

4.2.2 Thermal conditions, the role of permafrost and hydro-mechanical aspects

Ice-filled clefts and negative temperatures have strong hydro-mechanical effects on the stability of steep icy peaks. The ice-filling of cracks and fissures reduces the hydraulic permeability of the frozen rock mass to essentially zero. New cracks without ice may, however, form through slow rock deformation. Where water enters these open cracks, efficient advective heat transport, local rock wedging and high water pressures can result (see Hasler et al., 2012; Walter et al., 2019). The strength of negative-temperature rocks with ice-filled clefts is high but decreases with increasing temperature (Fig. 7). Critical conditions develop especially in “warm” permafrost with increasing amounts of unfrozen water at temperatures between about −1.5 and 0 °C. This is due to the fact that not only do the rock and ice components lose strength, but so do the rock–rock and rock–ice contacts, where temperatures and unfrozen water contents increase (Krautblatter et al., 2013). Complete permafrost thawing with elimination of ice then leads to a slight recovery of stability but still at a much lower level than with cold permafrost (Fig. 7).

Figure 7 Stability of frozen rocks with ice-filled clefts as a function of rock temperature from centrifuge experiments for two sets of slope angle α and inclination of discontinuity β (redrawn from Davies et al., 2001). Red arrows point to best guesses of permafrost temperatures at the Rasac detachment site.

The detachment site is in a very steep, WSW-oriented slope at an altitude of 5600 m a.s.l. reaching up to the ridge top at 5780 m a.s.l. (Fig. 1c). As such steep faces cannot accumulate a thick, thermally insulating snow cover, their mean surface temperature is close to the mean annual air temperature. Schauwecker et al. (2017) estimate the altitude of zero mean annual air temperature in the region to be about 5000 m a.s.l. Applying an environmental lapse rate of 6 ± 2 °C km⁻¹ (Stewart-Jones and Gruber, 2023) provides a rough estimate of the mean annual air and surface temperature at the altitude of the detachment site of around −4 ± 2 °C. This indicates cold and deep permafrost.

The N–S-oriented Rasac arête ridge induces pronounced thermal contrasts between its E- and W-facing slopes. Due to atmospheric humidity on average increasing during the day as a consequence of convection and cloud formation, E-oriented slopes, which receive morning sun, tend to be considerably warmer than W-exposed slopes receiving evening irradiation reduced by increased atmospheric humidity and clouds. This effect is independent of latitude. In the European Alps, the difference between permafrost occurrence on E- and W-oriented slopes is around 500 m as noted by Haeberli (1975) and implemented by Keller (1992); see later model approaches by Boeckli et al. (2012), or Kenner et al. (2019). This translates into a temperature difference between the two orientations of about 3 °C. Correspondingly adjusting the estimated mean surface temperature for W exposure by −1.5 °C leads to a best estimate of about −5.5 ± 2 °C at the W-oriented detachment site (Fig. 7). Applying a corresponding adjustment of +1.5 °C to the E-oriented slope of the ridge results in a mean surface temperature at the same altitude of −2.5 ± 2 °C. It follows that the temperature pattern is strongly asymmetric. The isotherms between the two sides with a temperature difference of about 3 °C must be near-vertical, causing heat flow inside the ridge to have a correspondingly strong horizontal component across the ridge from the warmer E- to the colder W-oriented face (see Noetzli et al., 2007). Such thermal asymmetries are not uncommon in detachment zones with permafrost of rock avalanches: Cathala et al. (2024) provide a detailed analysis of slightly warmer but comparable NW–SE asymmetric thermal conditions at the detachment site of the Étache rock avalanche in the French Alps, using in situ measurements and numerical model calculations. In the uppermost part of the Rasac detachment, where the ridge width is around 100–200 m, temperature gradients (around 1–3 °C per 100 m) indicate a horizontal heat flow component, which is comparable to the geothermal heat flow in flat topography under conditions of thermal equilibrium. Such horizontal asymmetries and related relatively weak heat fluxes nevertheless have the potential to become critical for slope stability in connection with warming-induced water infiltration from the warm side of ridges.

Due to atmospheric warming since the Little Ice Age, subsurface temperature gradients at depth, as documented by thermal monitoring of deep mid- and high-latitude boreholes, are generally reduced roughly within the uppermost 100 m below the surface (Etzelmüller et al., 2020). The observation by Thompson et al. (1995), that the coldest temperature of −5.2 °C in their borehole on the cold firn–ice saddle of Huascarán (6048 m a.s.l.) was measured at the ice–rock interface at 82.5 m depth, is in general agreement with the above-estimated surface temperatures and the warming-induced reduction of subsurface temperature gradients and related (vertical) heat flow as documented in mountain borehole temperatures at higher latitudes (Etzelmüller et al., 2020). The uppermost part of the Rasac arête ridge, including the detachment site, is most probably frozen throughout. Subsurface heating due to ongoing atmospheric temperature rise comes from both sides of the ridge. Only small remaining parts unaffected by colder historical/Holocene times may still exist deep inside the mountain (see transient 3D model calculations by Noetzli and Gruber, 2009).

4.2.3 Meteorological conditions before the GLOF

To investigate whether anomalous meteorological conditions contributed to triggering the event on 12 February, we investigated the evolution of different temperature and precipitation metrics over the study region, focusing on both the short term (a few days up to 2 months prior to the GLOF) and long term (2 years prior to 2023). In Fig. 8, the anomaly patterns compared to climatology (30-year average over 1981–2010) are shown because in extreme events it is important how anomalous the conditions are. The anomaly is calculated with respect to each month, i.e. the January 2023 anomaly is calculated by subtracting the January climatology from the January 2023 absolute values, to take away the influence of seasonality. Figure 8a presents temperature anomalies from January to February for variables including 2 m air temperature, skin temperature and soil temperatures at different depths. Before 12 February, temperature anomalies for 2 m air temperature, skin temperature and soil layers exhibit notable fluctuations. In early January, all variables are below the climatological average, with 2 m air temperature and skin temperature showing negative anomalies of approximately −1.5 °C, while soil temperature anomalies at all depths range between −1.0 and −0.5 °C. From late January to 10 February, all variables increase, with 2 m air temperature and skin temperature peaking at anomalies above +1.5 °C around 10 February. Soil temperature positive anomalies are more dampened, with surface soil (layer 1) reaching a maximum anomaly of about +0.8 °C, while deeper soil layers (3 and 4) exhibit a smoother, lagged and subtle increase up to ∼ +0.4 °C on 12 February. By 15 February, 2 m temperature and skin temperature anomalies drop to approximately −1.5 °C, with soil temperature anomalies also declining but more gradually. Surface soil (layer 1) anomalies decrease to −1.0 °C, while deeper layers (3 and 4) show smaller and lagged declines. This highlights the thermal buffering effect of deeper soil layers, which respond more slowly to atmospheric changes. Rainfall anomalies remain consistently negative, fluctuating around −7 mm, with a drop to −8 mm at the beginning of February. Snowfall anomalies, in contrast, remain near 0 mm throughout the period, suggesting limited deviation from the snowfall climatology. In the previous 2 years, the anomaly patterns align with broader seasonal cycles, with 2 m temperature and skin temperature peaking in summer and reaching lows in winter. Soil temperature anomalies are more dampened, with deeper soil layers (3 and 4) exhibiting even more subdued variations, rarely exceeding ±0.5 °C, reflecting the thermal inertia of deeper soil. In 2021 and 2022, precipitation anomalies are dominated by persistent negative anomalies. Periodic reductions in rainfall deficits are observed, also reflecting seasonal increases of precipitation in this region. Snowfall anomalies are negligible across the entire period, fluctuating around 0 mm, indicating minimal deviation from climatological average snowfall conditions. In summary, the most salient feature is that almost all temperature metrics investigated here (apart from deep soil layer temperature; see the Supplement) started to become anomalously warmer compared to the long-term mean (Fig. 8).

Figure 8 Anomalies in temperature and precipitation compared to the 1981–2010 mean climatology. Panel (a) shows the daily temperature anomalies for January and February of 2023, showing variables across 2 m air temperature, surface skin temperature and soil temperatures at different soil column depths. Panel (b) shows the corresponding precipitation anomalies during the same period. Panel (c) expands the temperature anomalies over 2 years prior to 2023, and panel (d) shows the corresponding precipitation anomalies over the same 2-year period. Dashed vertical line indicates the timing of the Rasac landslide and GLOF event on 12 February 2023.

4.3 GLOF process chain, hydrodynamics and transport of large boulders

4.3.1 GLOF process chain

The 2023 rockfall into Lake Rasac created a displacement wave that overtopped the crest of the lake's moraine dam at two locations and impacted the section of valley between Lake Rasac and Lake Gochacotan (Fig. 9c). The impact area between the crest of the dam of the GLOF-producing Lake Rasac and the downstream Lake Gochacotan covers 142 200 m². The reach of geomorphologically effective GLOF impacts is approximately 1.4 km with a vertical drop of 220 m (average slope of the trajectory 10°). The impact area is characterised by the presence of several large boulders with diameters > 3 m (Fig. 9d). These boulders were transported by the GLOF from the dam for a distance of > 600 m. Further geomorphic impacts include the erosion of existing stream networks as well as the creation of new branches of the stream network. The rockfall also hit a small bedrock-dammed lake located east of the moraine dam of Lake Rasac. High-resolution Planet Labs images reveal minor geomorphic impacts on the east-facing side of the ridge (part of the released material was deposited there, i.e. in the direction towards Yerupajá glacier). Interestingly, and fortunately, only very limited impacts are observed further downstream from Lake Gochacotan, and the GLOF did not reach the community settled near Lake Jahuacocha (3 km downstream from Gochacotan; see also Sect. 5.1).

Figure 9 Impact area of the 2023 process chain. Panels (a) and (b) show the location before (August 2019) and after (2023); (c) shows the disappeared lake and entire GLOF impact area and (d) shows detail of the impact area, where both erosional and depositional features associated with the 2023 GLOF are visible. B1 to B5 refer to large boulders analysed in Sect. 4.3.2. Maxar and CNES/Airbus satellite images available from © Google Earth collection.

4.3.2 Transport of large boulders

Large boulders (diameter > 3 m) were transported within the GLOF impact area, indicating that a transient, high-energy flow mobilised them. The boulders probably originated from the moraine dam, and their transport was facilitated by the steepness of the first section of the downstream face of the moraine (up to 30° in places) after wave impact. There are no apparent indications within the impact zone suggesting the occurrence of hyper-concentrated flows or debris flows, characterised by densities exceeding 1800 kg m⁻³ (Costa, 1988). According to Google Earth Pro imagery (July 2023) and in situ photographs of the post-GLOF impact area, no substantial morphological changes can be associated with any such high-density flows. Such flows typically induce substantial morphological alterations, resulting in the stripping of the entire surface layer of the landscape (Jakob and Hungr, 2005). Nevertheless, clear evidence of erosion and deposition processes in regions 2 and 3 (Fig. 9c) indicates a possibility of stratified flows. Given the above, it is reasonable to classify this event as a regular flood with a combination of water and sediment, displaying Newtonian fluid characteristics.

The minimum velocities required for translational and rotational motion and corresponding wave height are illustrated in Fig. 10. The results encompass a range of boulder dimensions (corresponding to width and length respectively) from 0.5 to 6 m, with observed boulder dimensions estimated from aerial imagery. Each curve represents a distinct region within our study site. It is assumed that the boulders originate from the dam crest with zero slope. For the largest observed boulder with a dimension of 5.6 m, the minimum velocity required for incipient translational motion is shown to be 8.8 m s⁻¹, with a corresponding wave height of approximately 7.9 m (see Fig. 10). Additionally, the minimum velocity for incipient rotational motion is determined to be 14.6 m s⁻¹, resulting from an approximate wave height of 21.6 m.

Figure 10 Flow velocities needed to transport the mapped large boulders.

5 Discussion

5.1 The limited GLOF impacts

The 2023 Rasac GLOF is characterised by limited reach and impacts (see Figs. 9 and 11). This observation is in contrast with the 2020 Salkantaycocha GLOF that was induced by an ice–rock avalanche of similar volume and characteristics (estimated volume 1.0×10⁶ m³–2.0×10⁶ m³, vertical drop 1000 m), but the reach of impacts exceeded tens of kilometres (Vilca et al., 2021). This could be explained by (i) dissipation of triggering rockfall energy on the bedrock slope before reaching Rasac Lake; (ii) possible gradual (not one-off) failure of the overhanging block (this may explain the two dam overtopping locations visible in Fig. 9c); (iii) the impact of a rockfall at the frontal part of the lake, diverting the primary direction of displacement wave towards the rear part of the lake; (iv) relatively small lake volume (0.3×10⁶ m³ in Rasac vs. 0.8×10⁶ m³ in Salkantaycocha); and (v) the attenuation role of Gochacotan Lake.

Figure 11 The GLOF impacts and Lake Rasac shortly after the GLOF. Panel (a) shows standpoints of the photographer and captured scenes; panel (b) shows a downstream view of the GLOF impact zone towards Lake Gochacotan, with large boulders analysed in Sect. 4.3.2 visible; panel (c) shows an upstream view of the GLOF impact zone towards the moraine dam of Lake Rasac; panel (d) shows the formed lake basin infilled with deposited material with still-visible water puddles remaining. Satellite image in panel (a) is the 2023 Maxar image available from © Google Earth collection. Field images are reproduced with kind permission of the author, a resident of the local Jahuacocha community, Javier Ramirez.

The Sentinel and Landsat time series reveal that Lake Rasac persisted in a limited areal extent shortly after the GLOF. The 19 March 2023 Sentinel scene (i.e. 4 weeks after the GLOF) shows a clear water body in the frontal part of the former lake basin, while no lake is identifiable on the next cloud-free Sentinel image taken on 23 March 2023. Such post-GLOF lake emptying is noticeably similar to what was observed in the case of the 2020 Salkantaycocha GLOF in Cordillera Vilcabamba (Vilca et al., 2021). The reasons for this unusual post-GLOF complete lake emptying may include (i) enlargement of underground outflow channels caused by increased hydrostatic pressure immediately after the GLOF, (ii) opening of new underground channels associated with erosion of distal face of moraine dam and (iii) decreased water inflow into the lake resulting from landslide deposition in its headwater area. Considering the deposition of rockfall material in the frontal part of the lake, option (iii) is rather unlikely.

5.2 The role of climate change and meteorological conditions in Rasac slope failure

Both the long-term change and the anomalous conditions prior to the event, in particular the deeper soil layer temperature change, suggest the role of permafrost warming in this event (see Sect. 4.2.3). However, the cold temperatures of the permafrost at the detachment site of the 2023 Rasac event mean that permafrost thawing was not likely an important factor. The detachment must have taken place within the deeply frozen rocks. A general stability reduction due to deep permafrost warming, however, must be assumed to have played a partial role in the sense that the geological conditions must already have been in a critical state. Furthermore, penetration of some water from the warmer, E-oriented slope towards the detachment cannot be excluded. Numerous small-magnitude precursory events (see Sect. 4.1 and 4.2) were also initiated from cold and deep permafrost but eliminated the cooling ice aprons, which must have caused additional local warming and frost wedging. These precursory events must have led to the final destabilisation of the main 2023 event despite the fact that this uppermost part of the ridge face remained covered by thin but thermally protective ice aprons.

5.3 Challenges in attributing the 2023 event to climate change

There are several challenges in attributing a single GLOF event such as that of 2023 to climate change. First, typically the initial step in an event attribution study is to define the event in question. The central aim is to link the event definition as closely as possible to the impacts of the event, and quite often the event definition is focused on impacts on humans and/or the built environment. But in this case, the actual impact is mostly felt downstream and with no identifiable “traditionally defined” impacts. Second, event attribution studies are mostly focused on key contributing meteorological variables (e.g. temperature for heatwaves; heavy precipitation for flooding; heat, drought and winds for fires). However, as demonstrated earlier, this sudden 2023 event was triggered by a mass movement, preceded by a few small rockfalls, instead of meteorological triggers. Although local temperature had been anomalously warm prior to the GLOF event, it was not the main trigger. Third, isolating the exact role of anthropogenic climate change for the 2023 event is challenging due to the various timescales of concern here. The anthropogenic influence on climate is undeniable when looking at decadal and centennial timescales, and glacier retreat, lake formation and permafrost degradation are all clear consequences of a warming planet. These long-term changes increase the overall risk of GLOFs. However, for a single and sudden event like Rasac 2023 – with complex cumulative developments and effects – a complex system reaction influenced by ongoing anthropogenic climate change, it is challenging to conduct a meteorological variable-based attribution. Rasac 2023 is a case of climate change interacting with various natural processes (rockfalls, mass movements etc.) creating the conditions for a GLOF.

A detailed and systematic attribution of glacial lake expansion and GLOFs to climate change has not yet been undertaken, and there are a few roadblocks to such efforts which we further discuss here. First and foremost, there is a need to identify the causal physical process pathway linking a GLOF event to weather and climate variables (concurrent with the GLOF and/or antecedent conditions) that trigger the GLOF or make the GLOF more likely. In the case of Rasac 2023, this could be by either exploring the immediate weather conditions to establish whether they serve as triggers or investigating whether the antecedent meteorological conditions contributed to the mass movement processes. The lack of clear understanding and targeted modelling efforts to unpack the role played by weather and/or climate in the physical process chain leading to GLOFs is a major challenge. Currently, there is a great need to fill this gap to better understand glacial lake evolution processes and to develop physical-process-based modelling of such processes. Second, in “traditional” event attribution studies, the intensity, frequency of occurrence, spatial extent and/or persistence (duration) of an event are investigated and compared across the factual scenario (i.e. the current state) and a counterfactual scenario (i.e. removing the influences of climate change) to arrive at an attribution statement of whether and to what extent climate change had made the event more, less or equally likely to occur (and/or intense if event intensity is the focus). For Earth surface processes such as landslides and GLOFs, defining the counterfactual scenario involves defining how the regional land surface and glacial lake evolution processes would have unfolded without the effect of climate change. Currently, there are a lack of efforts modelling and defining what these potential counterfactual scenarios could look like. In this regard, a promising way forward to attribute GLOFs to climate change is to use integrated process chain climate–glacier–land surface process modelling to address if and how much glacial lake expansion and GLOFs are attributable to past total global warming and to an anthropogenic component. Nonetheless, while a single GLOF event might be hard to attribute to climate change, the overall trend of increasing GLOF risk due to climate change is clear, with ongoing climate change leading to glaciers vanishing, lake formation and expansion, deep permafrost warming and degradation, and the resulting reduction in slope stability through glacial de-buttressing and reduced strength of frozen rocks (Emmer, 2024; Zhang et al., 2024).

With that being said, it is still valuable to start exploring some initial steps to perform attribution analysis of GLOFs. Therefore, in this study we still perform a simple observational-based analysis, attributing the anomalous temperature leading up to the event (which is the most salient feature through the climate section analysis), averaged over the 1° × 1° box $(-\mathrm{77.5}\mathit{°},-\mathrm{76.5}\mathit{°},-\mathrm{10.75}\mathit{°},-\mathrm{9.75}\mathit{°})$ covering Rasac. Since the anomalous temperature was seen in a few different variables (Fig. 8a) and to make sure the results are robust, the attribution analysis was undertaken for four variables: February maximum daily mean 2 m temperature anomaly from ERA5-Land, February maximum daily mean land surface skin temperature anomaly from ERA5-Land, February maximum daily mean top layer soil temperature anomaly and February 1 d maximum daily maximum temperature anomaly. The detailed methods for observational-based attribution are used according to the World Weather Attribution Protocol, described in Philip et al. (2020), with supporting details found in van Oldenborgh et al. (2021) and Ciavarella et al. (2021). We also briefly outline the methods in Sect. 3.2.2 for the readers' benefit. The results shown in Fig. 12 suggest that the likelihood of events, in terms of February temperature extremes, at least as extreme as what happened in February 2023, is increased due to global warming, with a return period of 1.2, 1.3, 1.3 and 1.3 years (for the four variables respectively) under 2023 climate conditions, a return period of ∼ 11, ∼ 20, ∼ 22 and ∼ 9 years under −1.2 °C climate conditions, and a probability ratio of ∼ 9, ∼ 15, ∼ 17 and ∼ 7. This means that global warming has increased the likelihood of a February temperature anomaly at least as severe as that over Rasac in 2023 February by almost 10–20 times, highlighting the impacts of climate change on local warming.

Figure 12 The GEV fit to the data at two levels of the covariate global mean sea surface temperature: in 2023 (red line) and in a 1.2 °C cooler climate (pre-industrial, blue line). The purple line shows the magnitude of the February 2023 event. The four panels show results for four different event definitions: (a) event defined based on February maximum daily mean 2 m temperature anomaly from ERA5-Land, (b) event defined based on February maximum daily mean land surface skin temperature anomaly from ERA5-Land, (c) event defined based on February maximum daily mean top layer soil temperature anomaly and (d) event defined based on February 1 d maximum daily maximum temperature anomaly from ERA5 (as ERA5-Land does not provide this variable). Data over 1950–2023 are used, and all anomalies are with respect to 1980–2010 climatology period.

5.4 Regional implications of the 2023 Rasac GLOF

Regional-scale insights about possibly changing occurrence patterns of GLOFs in space and time can be derived from systematic analysis of events with specified triggering events (type and magnitude) in a homogenous geographical region. In this study, we consider a “large” (volume >10⁶ m³, using the universal landslide size classification scheme of McColl and Cook, 2024) ice and rock mass failure, triggered in the glacierised cordilleras of the Peruvian Andes. Building on the regional GLOF inventory of Emmer et al. (2022b), we found four examples of such GLOFs that are documented from the Peruvian Andes in 21st century (see Table 3). Apart from the 2023 Rasac GLOF reported in this study, these are the 2002 Lake Safuna Alta rockslide-induced GLOF in Cordillera Blanca (Hubbard et al., 2005), the 2020 Lake Salkantaycocha ice–rock-avalanche-induced GLOF in Cordillera Vilcabamba (Vilca et al., 2021) and the 2022 Lake Upiscocha rockslide-induced GLOF in Cordillera Vilcanota (Vilca et al., 2022). Further, several GLOFs were triggered with similar processes but smaller initial volumes (Emmer et al., 2022b), for example, the GLOF triggered by an ice–rock avalanche at an unnamed lake located beneath the Panta peak in the Cordillera Vilcabamba in rainy season 2012/13, with an estimated landslide volume of >10⁵ m³, and the 2012 Lake No. 513 ice–rock-avalanche-induced GLOF in Cordillera Blanca, with an estimated landslide volume 2×10⁵ m³–5×10⁵ m³ of ice and rock (Carey et al., 2012; Vilímek et al., 2015; Emmer et al., 2022b), although some sources indicate an initial volume exceeding 10⁶ m³ (Valderrama and Vilca, 2010). The timing of these events implies that while there was on average one large GLOF triggered by a ice–rock mass failure per decade in the Peruvian Andes in 2000s, this figure increased 3-fold in the 2020s. Importantly, the rather rare nature of these events underscores the importance of careful reporting and analysis that is necessary to disentangle occurrence patterns and the role of climate change and associated processes (see Sect. 5.2).

Table 3 Large GLOFs triggered by ice and mass failure in the Peruvian Andes in the 21st century (from Emmer et al., 2022b).

We performed the same set of climate analysis for two GLOF events that occurred in the 2020s (2022 Upiscocha and 2020 Salkantaycocha), focusing on the longer-term trend over the last 8 decades and shorter-term anomalous conditions prior to the events, to determine whether there are similar meteorological patterns (results shown in Figs. S10–S36) and whether there are clear signals of weather and/or climate factors serving as triggers. Similar to Rasac, these two lakes have also seen sustained significant warming over the past 8 decades, with the 0 °C isotherm rising significantly over time, exposing more glacier areas to above-freezing conditions. Although the short-term conditions preceding the GLOFs at Lake Upiscocha and Lake Salkantaycocha were also characterised by elevated temperatures and precipitation deficits (similar to Rasac), the short-term variability, intensity and timing suggest distinct meteorological processes for each GLOF. Unlike 2023 Rasac, which had seen a clear peak in temperature anomalies right before the GLOF, 2022 Upiscocha and 2020 Salkantay did not experience such clear peaks right before the GLOF. In the case of 2022 Upiscocha, the positive temperature anomalies peaked by the end of July then dropped continuously to even negative anomalies right before 9 August. This makes it even harder to derive key contributing meteorological metrics to perform attribution analysis in the “traditional way”, further highlighting the need for future research to integrate climate–glacier–land surface process modelling to perform GLOF attribution.

6 Conclusions

The GLOF process chain reported in this study is the most recent large GLOF process chain triggered by mass movement (initial volume >10⁶ m³) in the Peruvian Andes. The triggering mass movement detached from cold, deep-reaching permafrost, and the failure of the main event must have occurred within the frozen rocks. These frozen rocks most likely underwent deep warming and a related stability reduction but not to especially critical conditions of warm permafrost with increased amounts of unfrozen water, which is corroborated by the observed period of 6 years of precursory small-magnitude rockfall activity. Our field observations emphasise the mobilisation of very large boulders (five boulders with diameter >3 m) and surprisingly limited geomorphic impacts downstream at Lake Gochacotan. While hydrodynamic analysis reveals flow velocities exceeding 5 m s⁻¹ for translational motion and 10 m s⁻¹ for rotational motion, understanding the attenuation role of floodplain topography and Lake Gochacotan requires detailed modelling. The 2023 Rasac GLOF in the context of similar events in the Peruvian Andes suggests that the frequency of such events has been increasing since the beginning of the 21st century. We assume that this increase may primarily be attributed to climate-change-induced ice effects, i.e. glacial de-buttressing and deep warming with stability reduction of perennially frozen rock slopes. These two dominant ice-related effects on changing preconditioning must be discussed in view of the challenges tied with attributing short-term triggering factors for individual events. In conclusion, we recommend that future research should especially address (i) quantitative analysis of stability conditions of lake-facing slopes with probable permafrost conditions, (ii) conceptual and modelling-based analysis of amplification and attenuation elements in GLOF process chains, and (iii) frequency–magnitude analysis of GLOF-triggering events in a changing climate.