Earthquake-induced debris flows at Popocatépetl Volcano, Mexico

The 2017 M_w 7.1 Puebla–Morelos intraslab earthquake (depth: 57 km) severely hit Popocatépetl Volcano, located ∼ 70 km north of the epicenter. The seismic shaking triggered shallow landslides on the volcanic edifice, mobilizing slope material saturated by the 3 d antecedent rainfall. We produced a landslide map based on a semi-automatic classification of a 50 cm resolution optical image acquired 2 months after the earthquake. We identified hundreds of soil slips and three large debris flows for a total affected area of 3.8 km². Landslide distribution appears controlled by the joint effect of slope material properties and topographic amplification. In most cases, the sliding surfaces correspond with discontinuities between pumice-fall and massive ash-fall deposits from late Holocene eruptions. The largest landslides occurred on the slopes of aligned ENE–WSW-trending ravines, on opposite sides of the volcano, roughly parallel to the regional maximum horizontal stress and to volcano-tectonic structural features. This suggests transient reactivation of local faults and extensional fractures as one of the mechanisms that weakened the volcanic edifice and promoted the largest slope failures. The material involved in the larger landslides transformed into three large debris flows due to liquefaction. These debris flows mobilized a total volume of about 10⁶ m³ of material also including large wood, were highly viscous, and propagated up to 7.7 km from the initiation areas. We reconstructed this mass wasting cascade by means of field evidence, samples from both landslide scarps and deposits, and analysis of remotely sensed and rainfall data. Although subduction-related earthquakes are known to produce a smaller number of landslides than shallow crustal earthquakes, the processes described here show how an unusual intraslab earthquake can produce an exceptional impact on an active volcano. This scenario, not related to the magmatic activity of the volcano, should be considered in multi-hazard risk assessment at Popocatépetl and other active volcanoes located along volcanic arcs.

1 Introduction

Earthquakes can induce large slope instabilities in tectonically active regions, resulting in a relevant source of hazard and damage. Earthquake magnitude (M) and the resulting intensity of ground vibration control the extent of the area where landslides may occur. One of the first comprehensive historical analyses of earthquake-induced landslides was done by Keefer (1984), who showed that the maximum area likely to be affected by landslides during a seismic event increases with M following a power law scaling relationship. In the following years, a growing number of studies started focusing on the impact of landsliding caused by large-magnitude earthquakes along relatively shallow crustal faults. In particular, it was observed that the fault rupture mechanism strongly influences the distribution of landslides, which are usually more abundant on the hanging wall in the case of reverse or normal faults (Sato et al., 2007) and present a symmetric distribution in the case of strike-slip faulting (Xu and Xu, 2014). In the case of co-seismic landslides related to earthquakes in subduction zones, very few data and inventories are available (LaHusen et al., 2020; Schulz et al., 2012; Serey et al., 2019; Wartman et al., 2013). Best examples are the landslides induced by the 2010 Chile megathrust earthquake (Serey et al., 2019) and by the 2011 Tohoku earthquake in Japan (Wartman et al., 2013). In both cases, thousands of shallow landslides were identified, but the main conclusion of these works is that the number of landslides generated by megathrust earthquakes is smaller than the number of events triggered by shallow crustal earthquakes by at least 1 or 2 orders of magnitude. Today, we know that the spatial distribution of earthquake-induced landslide is also a function of geological parameters (e.g., contrast in rock coherence, permeability), topography (slope and shape), land cover and land use, and ground-motion characteristics such as amplification and shaking frequency (Fan et al., 2019; Von Specht et al., 2019). Concerning the impact of earthquake on sediment-related hazards, a dramatic increase in sediment yield has been documented after large earthquake-induced landslides (Pearce and Watson, 1986; Dadson et al., 2004; Marc et al., 2019). Progressively, source areas on highlands can become quickly stable as fine material is removed and new vegetation grows and stabilizes the slope (Domènech et al., 2019), but debris flows can still occur due to the remobilization of deposited material along the channel (Fan et al., 2021).

On active volcanoes, a large variety of factors can promote slope instability and failure; these are magma intrusions, hydrothermal activity, gravitational spreading of the basements, climate fluctuations, and regional tectonics (Capra et al., 2013; Mcguire, 1996; Norini et al., 2008; Roberti et al., 2017; Roverato et al., 2015). In particular, earthquakes are recognized to play one of the most important roles in the initiation of slope failures on volcanoes (Kameda et al., 2019; Sassa, 2005; Siebert, 2002). Volcanic slopes that are close to a critical state can be particularly susceptible to perturbations produced by regional earthquakes. Volcanic landslides include a wide spectrum of instability phenomena, from small slope failures to large sector collapse evolving into catastrophic debris avalanches. Intermediate processes such as shallow landslides and debris flows are common in the case of an earthquake, but they are relatively poorly documented for past events. Debris flows, often called lahar in volcanic environments, are usually associated with eruptions that induce ice melt or snowmelt or with intense rainfalls occurring during intra-eruptive phases (e.g., Capra et al., 2018; Major et al., 2016; Manville et al., 2009). Few examples of long-runout debris flows triggered by earthquakes have been described on active volcanoes (Schuster et al., 1996; Scott et al., 2001). In Mexico, a M_w 6.5 earthquake that occurred in 1920 induced several landslides in the Pico de Orizaba–Cofre de Perote volcanic chain that transformed into debris flows with catastrophic effects for villages along the Huizilapan ravine (Camacho, 1920; Flores, 1922). More recently, several thousands of shallow landslides were triggered by the Tecomán earthquake of 21 January 2003 (M_w 7.6) in the volcanic highlands north and northwest of Colima City (Keefer et al., 2006).

In this paper, we investigate the exceptional mass wasting episode triggered by the 19 September 2017 M_w 7.1 Puebla–Morelos intraslab earthquake along the eastern and western sides of Popocatépetl Volcano. The seismic shaking mobilized pre-existing ash and pumice-fall deposits, producing hundreds of co-seismic soil slips. The largest ones had a total volume of about 10⁶ m³ and transformed into debris flows that traveled up to 7.7 km on the western side of the volcanic edifice. This phenomenon, never studied before at Popocatépetl Volcano and probably unique on an active stratovolcano along a continental volcanic arc, has important implications for hazard assessment, as the actual hazard map only includes the impact of lahars related to volcanic activity (Martin Del Pozzo et al., 2017). In the following, we provide a general introduction to the geomorphology of Popocatépetl Volcano and to its recent volcanic activity. Then, we describe the impact of the M_w 7.1 Puebla–Morelos earthquake on the volcano slopes in terms of ground vibrations and landslide activity. Finally, we reconstruct the transformation of the major landslides into long-runout debris flows, and we discuss the hazard implication for an active volcano.

2 Background

2.1 Popocatépetl Volcano

Popocatépetl Volcano (19^∘03^′ N, 98^∘35^′ W; elevation 5450 m a.s.l.) is located in the central sector of the Trans-Mexican Volcanic Belt (TMVB) and represents the active and southernmost stratovolcano belonging to the Sierra Nevada volcanic chain (Pasquaré et al., 1987) (Fig. 1a, b). The Popocatépetl is a composite volcano, and its present shape is the result of eruptive activity that rebuilt the modern cone after the 23.5 ka flank collapse (Siebe et al., 2017). During the Last Glacial Maximum (20–14 ka) glacier activity resulted in extensive moraines and glacial cirques (Vázquez-Selem and Heine, 2011). The lower part of the cone features a gentle slope (10–15^∘) and a dense vegetation cover up to approximately 3800 m a.s.l. (Fig. 2a), where pine trees became scattered and surrounded by dense tropical alpine grasslands (zacatonal alpino, Almeida et al., 1994) that can measure up to 1 m in height. Then, the cone becomes progressively steeper (20–30^∘) and unvegetated up to the summit. In the upper portion of the cone, the slopes are covered by abundant unconsolidated ash named “Los Arenales” from the recent vulcanian eruptions (Fig. 2a).

Figure 1 (a) Plate tectonic settings of Central Mexico (CP – Cocos Plate; NAP – North America Plate; RP – Rivera Plate; PP – Pacific Plate; MAT – Middle American Trench) and location of Popocatépetl Volcano (PV) in the Trans-Mexican Volcanic Belt (TMVB). (b) Details of the area affected by the M_w 7.1 Puebla–Morelos earthquake and location of the seismic stations CU and PPIG and of the rain gauge ALTZ. Popocatépetl Volcano (PV) is the southernmost edifice of the Sierra Nevada volcanic chain along with the Iztaccihuatl (IzV), Telapón (Te), and Tláloc (Tl) volcanoes. The main distribution areas of the pumice-fall deposits from the last two Plinian eruptions of the PV are indicated with white dashed lines (background image Landsat/Copernicus from Google Earth: © Google 2020, © INEGI 2020). (c) Strong motion recorded at station PPIG during the Puebla–Morelos earthquake on 19 September 2017 (channel HLZ refers to the vertical component, channels HLE and HLN to horizontal components).

Figure 2 Optical images of Popocatépetl Volcano acquired before and after the M_w 7.1 Puebla–Morelos earthquake (images © DigitalGlobe) and landslide map. In the pre-event image (a) acquired on 23 March 2017, the 3800 m line and Los Arenales (LA) deposits are indicated in white. In the post-event image (b) acquired on 11 December 2017, the red arrows indicate the debris flows that occurred in the Huitzilac, Hueyatlaco, and Xalipilcayatl ravines. In the landslide map (c), the black polygons corresponding to landslide scars and deposits were extracted from a Pléiades 1A image acquired on 13 November 2017. Main volcano-tectonic lineaments are reported in the map (red lines). The white triangle indicates the PPIG seismic station where a PGA of 106.83 cm s⁻² was measured (Fig. 1c). The black squares indicate the areas zoomed in Fig. 3. Background: map of the PGA distribution (data source: USGS National Earthquake Information Center, Worden et al., 2010) and 12.5 m DEM © JAXA/METI ALOS PALSAR 2008 (coordinate system WGS84 UTM Zone 14Q).

Quaternary volcanic activity of Popocatépetl Volcano has been characterized by catastrophic episodes including sector collapses and Plinian eruptions that emplaced pyroclastic density currents and thick pumice-fall deposits, predominantly toward the east and northeast (Fig. 1b) (Siebe and Macías, 2006). Based on its Holocene eruptive record, Plinian eruptions at Popocatépetl have occurred with a variable recurrence time of about 1–3 kyr (Siebe et al., 1996). Since 1994, the volcano entered a new eruptive phase, which includes dome growths that are subsequently destroyed during strong vulcanian eruptions with columns up to 8 km in height, accompanied with ash fall that has been affecting populations in a radius of approximately 100 km. Eruptive activity played the primary role in accelerating the glacier retreat on the northern slope of the volcano (Julio-Miranda et al., 2008).

In recent time, only two large lahar events were observed along the Huiloac Gorge (Hg, Fig. 1b), i.e., in 1997 and 2001, associated with eruptive phases (Capra et al., 2004). At those times, the Ventorillo glacier was still present on the northern face of the volcano. Both lahars propagated to the town of Santiago Xalitzintla (SX, Fig. 1b), located ∼ 15 km E of the volcano summit. The 1997 lahar originated after a prolonged explosive activity with emission of ash, which caused the partial melt of the glacier. The rapid release of water gradually eroded the riverbed and triggered a debris flow. The 2001 lahar originated from the remobilization of a pumice flow deposit emplaced over the Ventorillo glacier on the northern side of the volcano. The event occurred ∼ 5 h after the pyroclastic flow emplacement, and the debris flow was characterized by a stable sediment concentration of 0.75 (Capra et al., 2004). In the distal part, the 1997 lahar transformed into a hyperconcentrated flow, while the 2001 one maintained the characteristics of a debris flow due to its apparent cohesion conferred by a silty-rich matrix inherited from the pumice flow deposit. Apart from the Huiloac Gorge, which was characterized by significant geomorphic transformations due to these latter processes (Tanarro et al., 2010), most of the drainage network of Popocatépetl Volcano has a dense vegetation cover and presents stable, low-energy sediment transport conditions (Castillo et al., 2015). These stable conditions suddenly changed during the M_w 7.1 Puebla–Morelos earthquake.

2.2 The M_w 7.1 intraslab Puebla–Morelos earthquake

On 19 September 2017, central Mexico was hit by a M_w 7.1 intraslab seismic event (depth: 57 km) named the Puebla–Morelos earthquake (Melgar et al., 2018; Singh et al., 2018). The epicenter of the earthquake was located ∼ 70 km south of the summit of Popocatépetl Volcano and ∼ 100 km south of Mexico City (Fig. 1b). The focal mechanism corresponds to a normal fault with a dip angle of 44–47^∘ (Melgar et al., 2018). The 2017 M_w 7.1 Puebla–Morelos earthquake produced the most intense ground shaking ever recorded in Mexico City during a subduction-related earthquake and was the most damaging event for this densely urbanized part of the country since the 1985 M_w 8.1 Michoacán interplate earthquake, which occurred exactly 32 years before (Singh et al., 2018). The damage was surprisingly large in the critical frequency range for Mexico City (0.4–1 Hz), where the earthquake severely damaged hundreds of buildings and killed 369 people (Singh et al., 2018). The 2017 intraslab earthquake occurred closer to Mexico City, at greater depth, and involved a higher stress drop than its interplate counterparts, such as the 1985 Michoacán event. The stress drops of intraslab events have been estimated as being ∼ 4 times greater than that of the interplate earthquakes (García et al., 2005) and the ground acceleration of the intraslab earthquakes are expected to be more enriched at higher frequencies than those of the interplate events (Furumura and Singh, 2002; Singh et al., 2018). During the 2017 M_w 7.1 Puebla–Morelos earthquake, the peak ground acceleration (PGA) recorded at station Ciudad Universitaria (CU) was the highest recorded in the last 54 years of observations (57.1 cm s⁻²) (Singh et al., 2018). Station CU is located on the external boundary of the sedimentary basin responsible for the well-known seismic amplification at Mexico City (Fig. 1b). The strong ground motion recorded at PPIG station (Fig. 1c), located at 3980 m a.s.l. on Popocatépetl Volcano's slopes, featured a much higher value of PGA (106.83 cm s⁻², 0.1 g) than the one observed at station CU.

3 Data and methods

We adopted a combined field- and remote-based approach to retrieve information about the earthquake impact on such a difficult to access environment. Semi-automated satellite-image classification is a rapidly developing tool producing reliable landslide maps (e.g., Fan et al., 2019, and references therein). We used optical satellite data to identify the main areas affected by landslides and to constrain the timing of the landslide occurrences with respect to the earthquake event (Fig. 2a, b). We constructed a preliminary landslide map (Fig. 2c) based on the interpretation of an archive Pléiades 1A image (incidence angle of 14.63^∘, resolution of 0.5 m) acquired 2 months after the earthquake. A normalized difference vegetation index (NDVI) was calculated using band 1 (red) and band 4 (infrared). The resulting raster was classified for excluding vegetation cover, roads, and buildings from the analysis and selecting only landslide scars or depositional areas. The final map (Fig. 2c) was validated and refined based on data that we collected in the field. Most landslides are located on the W side, on the E side, and on the SE side of the volcanic edifice (Fig. 3).

Figure 3 Details of the most affected areas of Popocatépetl Volcano by co-seismic landslides. Black polygons correspond to shallow-landslide scars and deposits, and whitish areas indicate debris flows. Location of panels (a), (b), and (c) is reported in Fig. 2. Background: 12.5 m DEM © JAXA/METI ALOS PALSAR 2008 (coordinate system WGS84 UTM Zone 14Q).

Figure 4 (a) Map of three ravines of Popocatépetl Volcano affected by the debris flows; black dots indicate locations of field surveys (background: 12.5 m DEM © JAXA/METI ALOS PALSAR 2008; coordinate system WGS84 UTM Zone 14Q). (b) View in the downstream direction from the top of the sampling point PO06, located on the scarp of the Huitzilac landslide. (c) Topographic profiles of the larger landslides (1 and 2) that occurred at Huitzilac, which are also reported on the Pléiades 1A image acquired on 13 November 2017; the black arrow indicates the topographic barrier overtopped by the debris flow.

We conducted four field campaigns from October 2017 to November 2019 to investigate the morphology and stratigraphy of the source area of main landslides, to map and measure faults and fractures caused by the earthquake, and to define the extension, thickness, and textural characteristics of the larger debris flows (Fig. 4). The stratigraphy on the main landslide scars was reconstructed to determine texture and physical properties of the tephra layers involved in the mass wasting process. We selected a soil sample for radiocarbon analysis to identify the age of the stratigraphic sequence and to define its distribution. The ¹⁴C age was obtained through accelerator mass spectrometry dating (BETA Analytic Laboratory) and calibrated with the IntCal20 calibration curve (Reimer et al., 2020). We mapped and sampled main debris-flow deposits, and grain size analyses was performed by dry-sieving for the sand fraction and by means of a laser particle sizer (Analysette 22) for silt and clay fractions.

We analyzed two Sentinel-1 SAR images (synthetic aperture radar, Copernicus program) to define the timing between the earthquake and the observed mass wasting processes (Fig. 9). The analyzed images were acquired before and after the earthquake (17 and 23 September 2017) in 1A level ground-range-detected, ascending-orbit, interferometric wide-sensor mode and dual polarization. A radiometric calibration was applied to extract the most significant amount of backscattering information from the ground linked to the surficial roughness. As a second step, a change detection technique named log ratio was applied to detect pixel values directly related to radar backscattered and correlated to superficial processes; this is an algorithm used to detect changes using a mean ratio operator between two images of the same area but taken at different times (Mondini, 2017; Singh, 1989). Finally, we analyzed rainfall data gathered at the Altzomoni rain gauge station (ALTZ, Fig. 1b), located approximately 10 km north of the volcano summit.

4 Results

4.1 Landslide mapping

The earthquake triggered hundreds of shallow landslides on the volcano slopes covering a total area of 3.8 km² (Fig. 2). Soil slips affected the modern soil and part of the unconsolidated volcaniclastic cover. The five largest slope failures occurred in the basins of Hueyatlaco and Huitzilac on the west side of the volcanic edifice and in the basin of Xalipilcayatl on the east (Fig. 4). The scarps of these landslides were generated at elevations of about 3400–3800 m a.s.l. on the internal faces of ravines or glacial cirques, where slopes are > 20^∘ (Fig. 4c). Sharp rectilinear extensional fractures and small normal faults parallel to the valley slopes were observed in the Hueyatlaco basin after the earthquake (Fig. 5a). These faults and fractures have a maximum length of about 1 km, show displacements of up to 40–50 cm, and are located on the valley flanks (Fig. 5b), suggesting a correlation with local gravitational instability triggered by seismicity. A cluster of smaller shallow landslides is visible on the southwestern side of the volcanic cone (Fig. 2c). These landslides were produced by the collapse of the steep slopes of hummocky hills (Figs. 2, 3b) corresponding to the debris avalanche deposit of the last major flank failure that occurred at 23.5 ka (Espinasa-Perena and Martín-Del Pozzo, 2006; Siebe et al., 2017).

Figure 5 (a) Rectilinear extensional fractures and small normal faults opened parallel to the Hueyatlaco ravine (white arrows); background image: Pléiades 1A image acquired on 13 November 2017. (b) Detail of the normal displacement of about 50 cm (white arrow). Coordinate system WGS84 UTM Zone 14Q.

Figure 6 View of the three main landslides scarps: (a, b) Huitzilac, (c) Hueyatlaco, (d) Xalipilcayatl. In panel (e) the stratigraphic sections of the scarps and the grain size distributions of their strata (data in Appendix A) are reported; upper classic Plinian eruptions (UCPES, pink), lower classic Plinian eruptions (LCPES, orange), and the dated layer are indicated; see text for more details. (f) Geographic location of sampling points; background image: Pléiades 1A image acquired on 13 November 2017.

The scarps or the larger landslides located on the western slope of the volcano show a similar stratigraphy, with the intercalation of pumice- and ash-fall deposits (Fig. 6a–c). Pumice-fall deposits consist of open-framework, clast-supported units composed of gravel–sand-sized fragments of pumice embedded in a matrix of fine material (from silt to clay). Two main layers of pumice-fall deposit were observed at the Hueyatlaco and Huitzilac landslide scars (layers B and D, section PO1906; layers C and E, section PO1927; Fig. 6e). Another pumice-fall deposit crops out at the base of the pyroclastic succession. The fallout deposits are intercalated with massive or stratified ash layers, with variable thicknesses up to 4 m. They mainly consist of sand (71 %–93 %), silt (16 %–1 %), and less than 1 % of clay (see Appendix A). A sample from layer C (section PO06) was dated by using ¹⁴C, giving a calibrated age of 537–643 CE (1500 ± 30 BP conventional radiocarbon age) (Fig. 6e). Based on this age, the two younger pumice-fall deposits are here correlated with the upper and lower classic Plinian eruptions (UCPES and LCPES) of the late Holocene, which had a main dispersal axis towards the east and northeast (Fig. 1b) (Siebe et al., 1996). The thicker deposits of these eruptions crop out on the eastern flank of the volcano, as observed at section PO11, and correspond to the scar of the Xalipilcayatl landslide (Fig. 6d). Here, a main unit of pumice-fall deposit (C in Fig. 6d) features a total thickness of 3.5 m and consists of a massive, clast-supported unit dominated by coarse fragments barren of any silt and clay fractions. This latter unit lies on a 10 cm thick sandy layer (B in Fig. 6d). In all the studied sections, the upper ash unit corresponds to the products accumulated from the frequent vulcanian explosions that characterize the modern eruptive activity of the volcano.

4.2 Characterization of debris flows and associated deposits

The five largest landslides described in Sect. 4.1 (one each at Hueyatlaco and Xalipilcayatl and three at Huitzilac) mobilized a total volume of about 1.35 × 10⁶ m³ of ash- and pumice-fall deposits (Table 1). Landslide scarps measured 640 m in length and 4 m in depth at Hueyatlaco, 740 m in length and 3 m in depth at Huitzilac, and 400 m in length and 3 m in depth at Xalipilcayatl (Fig. 7). We calculated the volume of the landslides by assuming a constant depth (with an uncertainty of ±0.5 m) over the area of detachment. We measured the depth of the main scars in the field while the area of the main scars was inferred from field surveys and from the inspection of post-event optical images (Fig. 7).

Table 1 Main morphometric data of the landslides that occurred in the headwaters of the Hueyatlaco, Huitzilac, and Xalipilcayatl ravines. The area of the main scars was inferred from the inspection of post-event optical images (see Fig. 8). The depth of the scars was measured in the field. The volume of the landslides was calculated assuming a constant depth (with an uncertainty of ±0.5 m) over the area of detachment.

Figure 7 Comparison of 3D views on the detachment areas in the headwaters of the Hueyatlaco (a, b), Huitzilac (c, d), and Xalipilcayatl (e, f) ravines before and after the larger landslides. Images from Google Earth (© Google 2018, © INEGI 2018, and © DigitalGlobe 2019).

Figure 8 Debris-flow deposits in the upper (a–c), intermediate (d–f), and lower reaches (g–h) of the Huitzilac, Hueyatlaco, and Xalipilcayatl basins: (a) scarp of landslide A-1 at Huitzilac (view from point PO06), (b) main channel of the Huitzilac ravine (PO17), (c) main channel of the Hueyatlaco ravine (PO1701), (d) large wood deposits at Hueyatlaco (PO03), (e) overbank deposits at Hueyatlaco (PO02), (f) mud trace on lateral terraces at Huitzilac (HWM: height of watermark) (PO11), (g) evidence of dewatering at Huitzilac (PO19), (h) detail of the lower deposit at Xalipilcayatl (POE04), and (i) grain size distribution of the samples of the deposits.

The landslides transformed into three long-runout debris flows (Fig. 8). At the Huitzilac ravine, the main landslide body (landslide A-1, Fig. 7c) impacted the opposite side of the valley, partly overtopping it (Figs. 4c and 8a). Two other soil slips (landslide A-2 and A-3, Fig. 7c) contributed to forming the subsequent debris flow, which extended up to 7.7 km from the source before diluting into a streamflow. The total observed thickness of the deposit measures up to 3 m, but mud traces on standing trees and on lateral terraces measure up to 10 m on proximal reaches (PO17, Fig. 8b) and up to 1.5 m in distal reaches with horizontal surfaces at benches (PO11, Fig. 8f). In distal reaches, where the channel was shallow, the flow inundated large plains (PO15 and PO19). The deposit is massive and dark-gray in color and mainly consists of sand (77 %–86 %) with a relevant gravel proportion (15 %) due to pumice fragment enrichment in proximal reaches (Fig. 8i). Clay content is less than 1 %. The lower unit consists of coarse-to-medium ash with evidence of dewatering (Fig. 8g). At Hueyatlaco, the debris-flow runout extended up to 6.4 km (Fig. 4). The deposit appears as a main unit, dark-gray in color, massive, and homogeneous, with a sand fraction consisting of 70 % in proximal reaches (PO01) to 87 % in distal reaches (PO05), with up to 15 % of silt and less than 1 % of clay (Fig. 8i; see also Appendix A). Overbank deposits show sharp edges up to 10 cm thick (PO02, Fig. 8e). The total observed thickness is up to 50 cm (Fig. 8d; erosion was only incipient at the time of the observation) but watermarks of up to 5 m were observed in proximal reaches (PO1701, Fig. 8c). Finally, the deposit in the Xalipilcayatl ravine extended up to 1.5 km (Fig. 7f) and is clearly composed of two main units. The lower unit is massive and dark-gray in color and mostly consists of sand fraction (88 %, POE03-lower; Fig. 8i), up to 1.2 m in thickness, while the upper one is massive and pumice-enriched and represents up to 40 % of the total unit (POE04, Fig. 8h and i).

We estimate a total entrainment of about 205 000 m³ along both hillslopes and a channel network assuming 0.5 m of erosion over the area located downstream from the main scars (Table 2). Large wood (LW) elements entrained by the initial landslides and the subsequent debris flows contributed to the final bulk deposits of about 1.632 × 10⁶ m³. The volume of LW was calculated considering a mean tree height of 25 m (measured in the field, with an uncertainty of ±5 m), a mean trunk diameter of 0.4 m (observed in the field, with an uncertainty of ±0.1 m), and a mean distance of two trees of 10 m (estimated by using the post-event optical images; see Fig. 7). The amount of LW recruited in the Huitzilac basin results was 60 000 m³ (±3000 m³), far more than the sum of wood recruitment estimated for the Hueyatlaco (10 000 ± 500 m³) and Xalipilcayatl (7000 ± 350 m³) basins. The recruited LW stemmed from the combination of hillslope and channel processes originating from the earthquake-induced landslides. In general, these landslides were the dominant recruitment processes in headwaters. In contrast, LW recruitment from lateral bank erosion became significant in the intermediate reaches of the channels. The slope area that collapsed into the Xalipilcayatl basin contained most of the LWs that was later transported by the flow (86 %). In the Huitzilac basin, the LW recruitment mainly occurred on the slopes located right below the collapses (62 %), while in the Hueyatlaco basin most came from the channel banks (75 %). Most of the transported LWs remained trapped by natural obstacles in the main channel (i.e., standing vegetation) and clogged in the flat reaches of the channel (Fig. 8d). In the Xalipilcayatl ravine, most of LW was transported for the whole runout distance into the main landslide deposit.

Table 2 Main morphometric data of the debris flows that were observed in the Hueyatlaco, Huitzilac, and Xalipilcayatl basins. The entrained volume was calculated assuming 0.5 m of erosion over the area located downstream from the main scars where the vegetation was destroyed. The volume of large wood (LW) recruitment was calculated considering a mean tree height of 25 m (with an uncertainty of ±5 m), a mean trunk diameter of 0.4 m (with an uncertainty of ±0.1 m), and a mean distance of two trees of 10 m based on field observations and inspection of post-event optical images.

4.3 Timing of the events

Results of Sentinel-1 SAR image processing clearly indicate that both landslides and debris flows occurred between 17 and 23 September 2017. A binary image was produced where pixels values are linked to spatial change that occurred in this time span (Fig. 9a). Their distribution corresponds with the deposits of the larger debris flows that occurred in the Huitzilac and Hueyatlaco basins, as is easily observable in a later optical Sentinel-2 image (Copernicus program) acquired on 18 October 2017 (Fig. 9b).

Figure 9 (a) RGB (R, post-earthquake; G, pre-earthquake; B, ratio between post- and pre-event) representation of the 17 and 23 September Sentinel-1 (© Copernicus data) change in amplitude analysis. (b) RGB composition of post-event Sentinel-2 image (© Copernicus data). Coordinate system WGS84 UTM Zone 14Q.

Figure 10 Rainfall measurements at rain gauge ALTZ from 1 August to 4 October 2017. A total accumulated rainfall of 200 mm was recorded during the 30 d preceding the earthquake, 19.7 mm of which was on 17 September 2017 (red bar). On 4 October 2017, the population of San Juan Tehuixtitlán noticed the passage of sediment-laden flow in the Hueyatlaco ravine.

A total of 200 mm of accumulated rainfall were recorded during the 30 d preceding the earthquake, with the accumulation of 19.7 mm 2 d before the earthquake (Fig. 10). Thus, we expect that the slope material was wet at the time of the earthquake. Based on the remote-sensing analysis and considering that between 19 and 23 September only a few millimeters of rainfall accumulated (Fig. 10), it is thus highly probable that both slope failures and debris-flow emplacement were co-seismic. Witnesses from the town of Atlautla, which is located at the outlet of the Huitzilac ravine (Fig. 1b), also confirmed this information. During the following weeks, rainfall remobilized fine material from the landslide deposits reaching the town of San Juan Tehuixtitlán (Fig. 4a). On 4 October 2017, the population of San Juan Tehuixtitlán noticed the transformation of the shallow water flow of the Hueyatlaco ravine into a hyperconcentrated flow. It was the first time that this local community located on the western volcano slope observed such a phenomenon. Rainfall measurement at Altzomoni rain gauge station (ALTZ, Fig. 1b) shows an accumulation of 35.7 mm of rainfall over 12 h beginning at 10:00 UTC on 4 October, with a peak between 20:00 and 21:00 UTC (Fig. 10). The rainfall event of 4 October only remobilized fine material from the landslide deposits reaching the town of San Juan Tehuixtitlán; the debris flows along the Huitzilac and Xalipilcayatl were never reported since they never extended out to any populated area in 2017. During the 2018 and 2019 rainy seasons, the fine sediment remobilized from the debris-flow deposit in the Huitzilac ravine reached the road connecting San Juan Tehuixtitlán to Atlautla (Fig. 1b).

5 Discussion

5.1 Predisposing factors to slope instabilities

The Popocatépetl area is tectonically characterized by a Quaternary, roughly NE–SW- or ENE–WSW-trending maximum horizontal stress regime, responsible for arc-parallel E–W-striking transtensive faults and NE–SW or ENE–WSW arc-oblique normal faults (Arámbula-Mendoza et al., 2010; García-Palomo et al., 2018; Norini et al., 2006, 2019). This stress regime generated ENE–WSW extensional fracturing and faulting of the volcanic edifice (Fig. 11), controlling the orientation and propagation by magmatic overpressure of dikes within the volcanic cone and recent eruptive fissures on its flanks (Arámbula-Mendoza et al., 2010; De Cserna et al., 1988).

Figure 11 Simplified tectonic setting of the Popocatépetl area and location of the main earthquake-induced debris flows that occurred on 19 September 2017. Background image: © DigitalGlobe/ESRI 2019; coordinate system WGS84 UTM Zone 14Q.

The size of the slope failures triggered by the 2017 M_w 7.1 Puebla–Morelos earthquake is highly variable although (i) the epicenter of the earthquake is far from the volcano, with seismic shaking expected to be of similar intensity all over the symmetric volcanic cone, and (ii) soil and recent pyroclastic cover is quite homogeneous on the edifice flanks. Small shallow landslides occurred all over the volcano flanks, while the few larger landslides described in our work are limited to the eastern and western sides of the volcanic cone (Fig. 2). Thus, seismic shaking originating from the earthquake triggered large (volume > 10⁵ m³) landslides only in specific sectors of the volcano flanks.

The location of the larger slope failures defines a sharp ENE–WSW unstable sector, crossing the volcano summit and parallel to many deep rectilinear valleys carved in the volcanic cone (Fig. 11). In this ENE–WSW elongated sector of the volcano, some faults and extensional fractures were generated by the 2017 earthquake in the same basins where the larger landslides occurred (Fig. 5). This configuration suggests strongly localized site effects and/or a structural control on the location of the slope instability. Indeed, the unstable sector is roughly parallel to the ENE–WSW maximum horizontal stress, where local volcano-tectonic structural features are recognized on the volcano (Arámbula-Mendoza et al., 2010; De Cserna et al., 1988). The remobilization of larger quantities of material in this sector with respect to other areas of the volcano flanks may be correlated to the presence of ENE–WSW-striking faults and fractures that progressively weakened the volcanic edifice. Some of these volcano-tectonic structures may also have undergone transient reactivation by seismic shaking, increasing local slope deformation by the opening of fractures that promoted the largest slope failures triggered by the earthquake.

5.2 Initiation of co-seismic landslides

Slopes collapse when the shear stress across a potential failure plane exceeds the substrate strength. Earthquakes reduce the slope stability and can cause landslides through the perturbation of the normal and shear stresses in the slope. In the case of soft, saturated soils, the coalescence of cracks during earthquakes may results in liquefaction due to the increase in substrate permeability. At Popocatépetl Volcano, a combination of these two mechanisms produced the soil slips observed in the headwaters of the Hueyatlaco, Huitzilac, and Xalipilcayatl basins. Shapiro et al. (2000) already noticed that a large earthquake occurring in the vicinity of the volcano may result in flank instability because of the seismic waves traversing the poorly consolidated material composing the volcanic edifice. The ground motion during the 2017 earthquake was anomalously large in the frequency range 0.4–1 Hz, as intraslab earthquakes involve a higher stress drop than their interplate counterparts (Singh et al., 2018). Consequently, the ground motion is relatively enriched at high frequencies as compared with that during interplate earthquakes, which is dominated by lower-frequency waves (f < 0:5 Hz), and this effect can contribute to explaining the high value of PGA measured on the volcano slope.

Unexpectedly large peak accelerations have been recorded along crests of mountain ridges during several earthquakes (Davis and West, 1973; Meunier et al., 2008). Topographic amplification of ground vibrations is primarily due to the reflection or diffraction of seismic waves, which are progressively focused upwards (Bouchon et al., 1996; Davis and West, 1973). The constructive interference of reflections and the associated diffractions of seismic waves increase towards the ridge crest also due to local geologic factors, giving rise to enhanced ground accelerations on topographic highs (Del Gaudio and Wasowski, 2007; Meunier et al., 2008; Von Specht et al., 2019). Geli et al. (1988) show that topographic complexity (presence of neighboring ridges) may be responsible for large crest or base amplifications resulting in complex amplification–deamplification patterns and significant differential motions along the slopes. The amplification at the crest of a mountain can be as large as or larger than the amplification normally caused by the presence of near-surface unconsolidated layers (Davis and West, 1973). It is well-known that shallower earthquakes may cause large landslides (e.g., Marc et al., 2019), but the Puebla–Morelos earthquake was moderately deep (i.e., 57 km). The PGA produced by the 2017 earthquake at station PPIG (106.83 cm s⁻²) was about 2 times higher than the PGA observed at CU (57.1 cm s⁻²). Indeed, the distance epicenter–PPIG (68 km) is about half the distance epicenter–CU (111 km), and this partially explains the difference in PGA observed at the two stations. However, during the earthquake the headwaters of the Hueyatlaco, Huitzilac, and Xalipilcayatl ravines could have experienced even higher values of PGA due to the effect of the topographic amplification of seismic waves. The PGA map produced by the USGS Seismic Hazard Program shows values between 0.28 g in the southern sector of the cone and up to 0.18 g closer to the vent (Fig. 2c). The spatial interpolation of PGA clearly shows the interaction between the energy distribution and the topography, which played an important role in the location of landslides. The cluster of smaller landslides located on the southwestern side of the volcanic cone, closer to the epicenter, is likely due to the combination of large ground motion and high slopes that consist of debris avalanche hummocks (Fig. 3b).

The complex topography of Popocatépetl Volcano, characterized by neighboring ridges and valleys, probably produced local amplification values, which makes it difficult to explain why larger soil slips did not occur in other similar locations in terms of elevation, slope, and stratigraphy. However, the deposits located along the ENE–WSW unstable sector of the volcano (see Sect. 5.1), at an elevation ranging from 3400 to 3800 m and characterized by a slope > 20^∘, appear the most likely to suffer collapse in the case of an earthquake. This sector of Popocatépetl Volcano consists of a mantle of loose volcaniclastic material with the intercalation of silty–sandy ash layers and gravel–sand pumice-fall deposits (up to 5 m thick; see Fig. 6), covered by a modern soil with thick alpine grassland. At higher altitude, the steeper slopes are unvegetated and consist of unconsolidated pyroclastic granular material where superficial granular flows can be easily observed. The largest landslides occurred at the boundary of the vegetation line, where pine trees become scattered but grassland is still abundant (Fig. 7a, c). The intercalation of layers with different grain size and the soil coverage probably promotes water accumulation. Indeed, one mechanism that can possibly explain the collapse of this material is liquefaction through the disruption of internal, suspended aquifers. A similar observation was recently made at Nevado del Huila Volcano, Colombia, during 2007 when lahars originating from large fractures formed across the summit area of the volcano as a consequence of a strong hydromagmatic explosion that drained small, perched aquifers (Johnson et al., 2018). On the unvegetated portion of the cone, mass remobilization processes such as raveling and superficial granular flows likely occurred, but without leaving any scarp, because of the lack of a compacted soil.

5.3 Transformation into long-runout debris flows and implications for hazard assessment

Once generated, the earthquake-induced soil slips transformed into debris flows. The two major debris flows that occurred in the Hueyatlaco and Huitzilac basins covered a runout distance of 6.4 and 7.7 km, respectively. In Fig. 12, we show the conceptual model of this transformation at Popocatépetl Volcano: the propagation of an earthquake-induced crack in the saturated slope (1) produces a shallow landslide composed of a mix of ash and pumice. The collapsed material disaggregates and impacts on the opposite side of the valley (2), and rapidly the landslide evolves into debris flows, due to the high water content of the collapsed unconsolidated material (3). The subsequent debris flow is highly viscous due to the high sand and silt fraction of the mixture (Fig. 8i) and contains abundant LW entrained along the channel network, especially along the Huitzilac and Xalipilcayatl channels, which had entire mature trees incorporated, thus leaving abundant log-strewn debris. Even if no direct observations are available to assess whether the collapsed slopes were partially or completely saturated, it is clear that debris flows contained a large amount of water as observed from dewatering features of the deposits and high-water marks along the channels (Fig. 8g). Beginning on 21 August 2017, 138 mm of rainfall accumulated continuously for 2 weeks, with 19.7 mm just 2 d before the earthquake (Fig. 10). This large amount of rainfall was then stored in the open-framework pumice-fall deposits intercalated by meter-thick sandy layers and in the root fabric of the trees in the dense forest cover.

Figure 12 Conceptual model of transformation of earthquake-induced soil slips into debris flows at Popocatépetl Volcano: (1) the earthquake produces the collapse of the saturated slope composed of a mix of ash and pumice; (2) the landslide impacts on the opposite side of the valley entraining a large amount of large wood (LW) and (3) evolves into a debris flow due to the high water content of the material. The simplified stratigraphy in (1) reflects the one observed at the scarp of the Huitzilac landslide (see Fig. 6b).

Volcanoes store or drain water in and through aquifers that can grow and empty as impermeable barriers develop or as they are breached by deformation, respectively (Delcamp et al., 2016). Even if not completely saturated, ground vibrations induced positive pore pressure and triggered liquefaction and slope failure (Kameda et al., 2019; Wang et al., 2019). It is important to note that on 10 August a rainfall of 35 mm, similar to the 4 October event that triggered the sediment-laden flow observed at San Juan Tehuixtitlán town, did not induce any channel response, indicating the stability of the slopes of this sector of the volcano prior to the earthquake. In fact, except for the 2010 lahar that occurred in the Nexplayantla ravine after 100 mm of accumulated rainfall (Zaragoza-Campillo et al., 2020), lahars are related to major eruptions, which explains why the hazard map of Popocatépetl Volcano includes only rainfall-triggered lahars during or after eruptions (Martin Del Pozzo et al., 2017). Detailed field investigations of the role of aquifers on volcanic landslides are very scarce to date (Delcamp et al., 2016). Knowledge of the distribution of perched aquifers and the water content of volcanic deposits can provide precious insights into a complex mass wasting chain like the one experienced at Popocatépetl Volcano in 2017.

Finally, it is worth mentioning that during our last field campaign in November 2019 we observed that the source areas of the larger debris flow that occurred at the Huitzilac ravine were becoming stable as a result of the combined effect of the removal of fine material and of the growth of new vegetation. By contrast, the large amount of material deposited in the channel remains available for remobilization for many years, resulting in a remarkable increase in sediment yield as observed in other locations (e.g., Fan et al., 2021). Indeed, during the 2018 and 2019 rainy seasons, the fine sediment remobilized from the debris-flow deposits in the Huitzilac ravine reached the road connecting San Juan Tehuixtitlán to Atlautla.

6 Conclusions

The catastrophic event of 19 September 2017 at Popocatépetl Volcano is an exemplary case of interrelated multiple hazards in volcanic environments: earthquakes, landslides, and sediment-laden flows. During the M_w 7.1 Puebla–Morelos intraslab earthquake, hundreds of shallow landslides were triggered on the volcano flanks. The combination of strong ground motion due to local amplification with the presence of water-saturated, tephra-rich superficial deposits resulted in large slope failures and subsequent liquefaction of the collapsed material. A total volume of about 10⁶ m³ of volcaniclastic deposits transformed into two large debris flows on the western slope of the volcano and one on its eastern side. While the source areas rapidly stabilized in the months and years following, the fine material deposited in the channels remains exposed to possible remobilization for many years. These observations imply the need to revise the hazard assessment for Popocatépetl Volcano, where multi-hazard risk scenarios should be taken into account, as well as for other volcanic settings. The mass wasting cascade described here may occur in other locations, especially in continental volcanic arcs and mountain chains located in tectonically active regions.

Appendix A

Table A1 Grain size distributions of samples collected in the landslide scarps and deposits; cutoff particle sizes: gravel 64–2 mm, sand 2 mm–64 µm, silt 64–2 µm, clay < 2 µm. Refer to Fig. 4 for sample locations.