Modulation of Popocatépetl’s activity by regional and worldwide earthquakes

Volcanoes switching from quiescence to eruption shortly after catastrophic earthquakes have raised interest for volcanic triggering and the influence of earthquakes on volcanic activity. Its influence on already active systems and especially at open-vent volcanoes is more difficult to apprehend. A number of recent observations suggest an influence of tectonic earthquakes on Popocatépetl’s activity, the importance of which remains unknown. To further investigate this, we introduce an index, based on the near-field concept, identifying the earthquakes with the highest potential to promote volcanic activity (hereafter termed “significant earthquakes”). The time series of significant earthquakes is compared with the intensity of the volcanic activity, as characterized by the number and energy of volcano-tectonic earthquakes, the number of dome extrusions, the intensity of thermal and degassing fluxes, and ash production. Three main periods with contrasting activity stand out showing that Popocatépetl presents intense activity when significant tectonic earthquakes are frequent. Enhanced extrusion apparently follows significant earthquakes quickly with pulses of dome extrusion that peak after 1.3 ± 0.3 years. Conversely, extrusive activity vanishes when significant seismicity disappears, as during the period 2003–2011, which coincides with a 12-year-long significant seismicity gap. Hence, we propose that the 1994–2022 open-vent activity at Popocatépetl is in part modulated by the repetitive occurrence of significant earthquakes that periodically promote volcanic activity.


Introduction
It has long been observed that regional or volcano-distal earthquakes can affect volcanic activity (e.g. Darwin 1840;Rockstroh 1903;Yokoyama 1971;Nakamura 1975;Carr 1977;Gudmundsson and Saemundsson 1980;Marzocchi et al. 1993;Linde and Sacks 1998;Hill et al. 2002;Watt et al. 2009;Eggert and Walter 2009;Prejean and Haney 2014). Such earthquakes are now widely acknowledged as potential triggers for nearby and remote volcanic reactivation or fluctuations in eruptive activity at open-vent systems (e.g. Cigolini et al. 2007;Delle Donne et al. 2010;Pritchard et al. 2014;Coppola et al. 2015). Most of these studies use one or several parameters to characterize volcanic activity. This mostly involves the use of various seismic parameters such as the appearance, cessation, or rate of changes in volcanic seismic signals including tremor, [very] long period signals, hybrid signals, and volcano-tectonic earthquakes (Hill et al. 1993;Ukawa et al. 2002;Carniel et al. 2003;Moran et al. 2004;Sanchez 2004;West et al. 2005;Hill and Prejean 2005;Okubo and Wolfe 2008;Pritchard et al. 2014;Mora-Stock et al. 2014;Prejean and Hill 2018;Butcher et al. 2021). It also includes consideration of changes in gas and/or thermal fluxes (Benson and Follet 1986;Harris and Ripepe 2007;Cigolini et al. 2007;Delle Donne et al. 2010;Coppola et al. 2015;Avouris et al. 2017), extrusion rates (Cigolini et al. 2007;Coppola et al. 2015;Carr et al. 2018), deformation (Pritchard et al. 2013(Pritchard et al. , 2014Takada and Fukushima 2014), and/or explosions frequencies (De la Cruz-Reyna et al. 2010;Coppola et al. 2015). Fewer cases describe the decrease or even cessation of volcanic activity after an earthquake, with Kimura (1978), Ebmeier et al. (2016), and Farías and Basualto (2020) being exceptions. A few studies have also reported changes in eruptive style following a regional earthquake (Carniel et al. 2003;Ortiz et al. 2003;Carr et al. 2018).
However, to prove a relation between a trigger event (i.e. a distant earthquake) and an eruptive response (i.e. a change in activity), recourse to a statistically robust dataset for, and index to describe, cause and effect is needed. Many studies have attempted this in a global sense (Carr 1977;Linde and Sacks 1998;Marzocchi 2002;Marzocchi et al. 2004;Manga and Brodsky 2006;Walter and Amelung 2007;Harris and Ripepe 2007;Eggert and Walter 2009;Delle Donne et al. 2010;Avouris et al. 2017;Nishimura 2017;Sawi and Manga 2018). We here attempt this using an extensive and robust data set for a single volcano characterized by open-vent activity, i.e. Popocatépetl, and test an index to describe the earthquake trigger.
The response of such an open-vent persistently degassing system to tectonic earthquakes ranges from small variations of one or several parameters (cf. De la Cruz-Reyna et al. 2010) to significant changes in the level of activity (cf. Delle Donne et al. 2010). The occurrence of trigger events can thus promote the modulation of its level of activity. Here, we focus on significant changes in activity (i.e. switching from passive degassing activity to sustained extrusion of domes following tectonic earthquakes) to evaluate the influence of volcanic triggering by tectonic earthquakes on the ongoing eruption at Popocatépetl.

Earthquake-triggering volcanic activity
Single-case and short-term reactivation (i.e. within a few days of the trigger event) has received most attention in the literature, these interactions being more easily identified in volcanic records. Such immediate responses suggest that reactivated volcanoes had to be poised to erupt in order to respond to the triggering event (Barrientos 1994;Marzocchi 2002;Hill et al. 2002;Manga and Brodsky 2006;Walter and Amelung 2007;Watt et al. 2009;Eggert and Walter 2009;Bebbington and Marzocchi 2011;Nishimura 2017). Using a reliable seismic catalogue spanning 1964 to 2016, Sawi and Manga (2018) recently showed that: (1) ~ 4% of the explosive eruptions (VEI ≥ 2) occurred shortly (within 5 days) after nearby (< 800 km) earthquakes (M w ≥ 6) and were thus potentially triggered.
(2) This proportion is not significantly different from what would happen without the earthquake trigger.
When considering longer-term time windows, Sawi and Manga (2018) found a 5-12% increase in eruption number was apparent within 2 years of the triggering earthquake. This is in agreement with other studies that suggest increased volcanic activity 1-5 years following the trigger (Marzocchi 2002;Walter and Amelung 2007;Watt et al. 2009;Nishimura 2017). Such durations are also consistent with timescales of mechanisms proposed to account for buildup of pressure within the magma chamber prior to eruptions that follow large earthquakes (Manga and Brodsky 2006;Seropian et al. 2021). Although eruptions delayed by a few years (Bautista 1996;Nishimura 2017) and up to 30-35 years (Marzocchi 2002) have been described, delayed responses may not be distinguishable for chance coincidence between the large number of earthquakes and volcanic eruptions that occur every year (Sawi and Manga 2018). Thus, earthquake-triggering of volcanic activity can only be proved through statistically robust analyses of longterm multiparametric time series for both tectonically related seismic events and volcanic activity (cf. Delle Donne et al. 2010).
Open-system volcanoes are usually in a delicate state of balance, which makes them very sensitive to small perturbations (Rose et al. 2013). Although their activity is persistent and relatively constant in the long-term, they frequently undergo short-term variations in activity and/or eruptive style in the short-term (cf. Vergniolle and Métrich 2021;Edmonds et al. 2022). A number of studies have highlighted the sensitivity of open-vent systems to external perturbations in their magmatic system (Patanè et al. 2007;Sottili and Palladino 2012;Carey et al. 2012;Eychenne et al. 2015). This includes the system response to tectonic earthquakes (Linde and Sacks 1998;Ortiz et al. 2003;Harris and Ripepe 2007;Watt et al. 2009;Delle Donne et al. 2010;Avouris et al. 2017;Nishimura 2017;Carr et al. 2018). However, volcanic reactivation due to an earthquake trigger are often delayed by a few days, weeks, months, and even years (Eggert and Walter 2009;Bonali et al. 2013;Nishimura 2017;Sawi and Manga 2018).
Here, we use the open-vent eruption of Popocatépetl volcano (Mexico) to assess earthquake-triggering. Popocatépetl eruption began in December 1994, and the system is located in the actively seismic Central American subduction zone. Thus, there is an extensive and detailed record for tectonic seismicity and volcanic activity available for analysis. Although Popocatépetl does not meet the "sensitive volcano" criteria of Sawi and Manga (2018), i.e. the volcano needs to have erupted at least once within 5 days and 800 km of a M6 or greater earthquake since 1964, several observations suggest that this volcano might actually be sensitive to earthquake-triggering. Two cases have previously been described, including the volcano-tectonic earthquake swarm that followed the 1999 M w 7.0 Tehuacán earthquake within hours (De la Cruz-Reyna et al. 2010) and the early 2005 reactivation following the 2004 M w 9.2 earthquake of Sumatra (Delle Donne et al. 2010). More recently, interest in earthquake-triggering of activity at Popocatépetl was renewed by a series of short-to midterm responses: (1) The coincidence of exhalations with some of the aftershocks of the 20 March 2012 M w 7.4 Ometepec ( Fig. 1) earthquake ; Ferres, personal communication) (2) An ash emission occurring 20s after the arrival of the P wave of the 19 September 2017 M w 7.1 of Puebla-Morelos earthquake, followed 8 days afterwards by a 5-h-long isolated ash-venting episode (where such events had not occurred in several months before the earthquake) (3) The paroxysm of July 2020, 29 days after the 23 June 2020 M w 7.3 earthquake of Oaxaca In this study, we aim to assess the effect of significant earthquakes through the analysis of the 26-year-long records for the tectonic and volcanic activity available for Popocatépetl. In doing so, we set up and test a data processing methodology and index, which provides means of assessing the triggering potential of an earthquake at a recipient volcano. We scaled this index with all global published cases of earthquake-triggered activity and define a "significant earthquake", as an event with an index in the range where > 90% of the published cases occurred. Pulses of sustained dome extrusion shortly followed the earthquakes with the strongest potential to trigger an eruptive response. Instead, no change to the persistent open-vent activity was apparent when no "significant earthquakes" occurred.

Regional and tectonic setting
Popocatépetl is located within the Trans-Mexican Volcanic Belt (TMVB) (Fig. 1), which results from the subduction of the oceanic Cocos plate under the continental North America plate (Ferrari et al. 2012). This 1000-km-long volcanic arc comprises thousands of monogenetic edifices, some caldera complexes, stratovolcanoes, and geothermal fields (Gómez-Tuena et al. 2007;Ferrari et al. 2012Ferrari et al. , 2018. Tectonic deformation within the TMVB is mainly extensional (Ego and Ansan 2002;Suárez et al. 2019).
The Cocos plate consists of a flat slab which is up to 250 km long in its central section, terminating with a steep 75° dipping distal segment (Pardo and Suárez 1995;Pérez-Campos et al. 2008;Manea et al. 2017), which allows melt production beneath the TMVB (Manea and Manea 2011). This flat slab segment is laterally connected to a moderate dipping zone associated with normal subduction (Pardo and Suárez 1995), possibly through slab tearing at depth (Dougherty and Clayton 2014). Trenchward migration of volcanism suggests that following the tear, the slab is rolling back (Gómez-Tuena et al. 2003;Ferrari et al. 2012).
In this context, interplate earthquakes are dominated by thrust mechanisms down to a depth of 35 km. At greater depths, the intraslab mechanisms are normal (see Fig. 3b of Sawires et al. 2021). Due to the steep dip of the slab termination, subduction seismicity is absent beneath the TMVB (Pardo and Suárez 1995;Suárez et al. 2019;Sawires et al. 2021). To the west, the rift of the Gulf of California and the East Pacific Rise generates earthquakes with strike-slip mechanisms.

Recent activity of Popocatépetl and interplay between dome extrusion and seismic and degassing activity
The Popocatépetl volcanic complex began its construction < 730 kyr ago (Conte et al. 2004;Sosa-Ceballos et al. 2015;Sunyé-Puchol et al. 2022) and comprises at least three edifices (Espinasa-Pereña and Martín-Del Pozzo 2006). Its history features destructive events including flank collapses and at least nine Plinian eruptions between 23 and 1.2 ky ago (Siebe et al. 1996(Siebe et al. , 2017Arana-Salinas et al. 2010). About 15 eruptions occurred during the last 500 years (De la Cruz-Reyna and Tilling 2008), the penultimate of which ended in 1927. Following 67 years of quiescence, and concluding a progressive increase of hydrothermal activity since late 1991 (Gerlach et al. 1997), a series of phreatic explosions initiated the ongoing eruption on 21 December 1994.
Since March 1996, Popocatépetl's activity has been characterized by cycles of extrusion and destruction of andesitic to dacitic lava domes within the summit crater, with residence time lasting from a few days to 2 years (Gómez-Vazquez et al. 2016). Dome extrusion is usually spasmodic, with individual extrusive episodes lasting several hours to several days and being associated with "trains of exhalations" (cf. Gómez-Vazquez et al. 2016). They consist of ejection of incandescent blocks to distances of more than 500 m and ash-venting episodes (Taquet et al. 2017(Taquet et al. , 2019Campion et al. 2018). Such events are associated with continuous low-frequency and spasmodic high-frequency tremor episodes Arámbula-Mendoza et al. 2016). The domes are usually partially or completely destroyed by subsidence and/or explosive activity that last several days (Gómez-Vazquez et al. 2016;Campion et al. Fig. 1 Location of Popocatépetl volcano (red symbol) and of the earthquakes that occurred during the study period (1992-2020) with a triggering index (d/√S) of < 40. The names, date, and magnitude (M w ) of the main earthquakes described in the text are indicated. The green zone represents the Trans-Mexican Volcanic Belt following Ferrari et al. (2012). Scales are in km 2018). Rare isolated moderate explosions also occur during inter-cycle activity Campion et al. 2018). Gómez-Vazquez et al. (2016) identified changes in the extrusive behaviour of Popocatépetl and defined five regimes of dome emplacement alternating between high and low emplacement rates. The regime with the highest rates, between December 2000 and April 2003, preceded a long period, from April 2003 to August 2011 during which just a few domes were emplaced with low extrusion rates. From September 2003 to December 2004, the level of activity strongly decreased with only 1-3 exhalations per day, so that the authorities reduced the volcanic alert (Espinasa-Pereña 2012). We refer to these 15 months as an "eruptive pause".
How may earthquakes affect volcanic activity at Popocatépetl? Seropian et al. (2021) introduced a classification of earthquake-triggered volcanic events based on three parameters: (1) Open-system degassing (2) Magma viscosity (3) Presence of a hydrothermal system They also proposed a list of possible processes through which the activity could be affected following triggering earthquakes. We apply this threefold classification here: (1) Popocatépetl volcano is considered an open-conduit system because of its high and persistent gas emissions (Goff et al. 1998;Carn et al. 2017;Campion et al. 2018). These result in the excess degassing problem, where the emitted amount of lava is less than the mass of degassed magma (Andres et al. 1991;Christopher et al. 2010;Coppola et al. 2022). Witter et al. (2005) proposed a degassing model for Popocatépetl in which the excess is produced through convection of magma in the conduit (cf. Kazahaya et al. 1994;Shinohara 2008;Moussallam et al. 2015). (3) Prior to the current eruption, an intensification of the hydrothermal activity was observed, with the presence of an acidic crater lake (Armienta et al. 2000), sulphur deposits, and high-temperature fumarole fields (Goff et al. 1998  ) -Dynamic stress changes with, e.g. bubble nucleation, sloshing or resonance (Namiki et al. 2016(Namiki et al. , 2019 The dynamic stress processes would, however, probably affect the deeper portions of the magma system but not in the conduit because of the higher viscosities in the shallowest portions of the plumbing system (Witter et al. 2005).
Method: an index for the triggering potential Following Hill et al. (1993), it is commonly assumed that local stress changes during and/or following earthquakes may alter the volcanic activity (Nostro et al. 1998;Walter 2007;Bonali 2013). Earthquakes generate seismic waves, including transient dynamic seismic waves (P and S far-field body waves, near-field waves, Love and Rayleigh far-field surface waves) and near-field static displacement. The latter can be considered a zero-frequency seismic wave. These waves generate dynamic and static stress fields, which are commonly invoked to explain the triggering of volcanic eruptions (Hill et al. 2002;Manga and Brodsky 2006;Eggert and Walter 2009;Seropian et al. 2021). The dynamic stress field is generated by the passage of dynamic waves that propagate in the near-and far-fields. Instead, the static stress field is generated by permanent dislocation over the fault plane (the static displacement field), which attenuates with distance following 1/d 3 , such that it is only significant in the near-field. The concept of near-and far-field is not an absolute value of distance to the earthquake but instead needs to be described by a relative distance with respect to the length L of the rupture plane generated by the earthquake (Hill et al. 2002). The ratio d/L is thus classically used in seismology to define the near-field, corresponding to distances that scale with the order of magnitude of the rupture length (Aki and Richards 1980;Vidale et al. 1995;Hill et al. 2002), i.e. d < 10 × L.
The potential for an earthquake to provoke a change in volcanic activity thus mainly depends on two parameters: the distance between earthquake epicentre and the volcano and the magnitude of the earthquake. As the length of the rupture is a function of the magnitude, these two parameters are implicitly included in the nondimensional ratio: The smaller the d/L ratio, the stronger the ground acceleration at the point of measurement and thus also the greater is the earthquake potential to trigger a response at the recipient volcano. Whether a response is triggered depends on the internal state of the volcano (Bebbington and Marzocchi 2011;Hamling and Kilgour 2020;Seropian et al. 2021).
Length is not a fully satisfactory parameter because of its dependence on earthquake magnitude and focal mechanism. For example, given two earthquakes with the same magnitude, a strike-slip earthquake will have a longer L than a subduction earthquake (Blaser et al. 2010). The size of the rupture area is also linearly proportional to the moment magnitude (Blaser et al. 2010). To take these effects into account, we thus use √S (where S = L × W , i.e. the surface area of the rupture plane) instead of L and define an index using: To apply this index, the distance between each earthquake and the receiving volcano is calculated in spherical coordinates and converted into km, and S is estimated as follows. First, we used the published values for L and W, or S, when available for large (M w > 5.9) earthquakes (Mendoza and Hartzell 1999;Pollitz et al. 2011;Melgar et al. 2021). When no data was available from the literature, the source parameters were calculated from the scaling relationships of (1) d∕L Blaser et al. (2010). This was applied to earthquakes within the M w focal mechanism sorted ranges defined by Blaser et al. (2010). That is: -4.8-9.5 for a thrust mechanism -5.1-8.4 for a normal case -4.6-8.1 for a strike-slip event These regressions are designed for subduction environments and thus are well suited to our case. For earthquakes with magnitudes less than the lower M w range of the Blaser et al. (2010) scaling relationships, we calculated L using the following relation ): In applying this, we assume a square fault (i.e. W = L , Scholz 1982).
There is no benchmark for our new index, and we propose to calibrate it for the earthquake-triggered volcanic activity using the cases already described in the literature. For this, we calculated the index values, as given above exposed, for all published cases of tectonic earthquakes that apparently triggered an eruptive response at volcanoes globally (Table S1).

Regional and worldwide catalogue of earthquakes with the potential to favour volcanic activity at Popocatépetl
To characterize the triggering potential of regional and/or distal earthquakes at Popocatépetl, we compiled an earthquake catalogue between 1992 and 2020. It is based on the Mexican National Seismological Service (SSN; UNAM, 2021-www. ssn. unam. mx) and the Advanced National Seismic System (ANSS), Comprehensive Catalogue of Earthquake Events (US Geological Survey 2021; https:// earth quake. usgs. gov) for the distal earthquakes (M w ≥ 7 for worldwide earthquakes, M w ≥ 6 for Central American earthquakes). When regional earthquakes were in both catalogues, we gave our preference to the locations determined by the local SSN network because they are better determined than those obtained by a worldwide network (cf. Legrand et al. 2021). The locations of distal earthquakes are those of the ANSS or of local networks. Magnitudes determined by the broader ANSS network are more accurate than those calculated by a local seismic network, so they were always used when available. For specific earthquakes studied in detail in the literature, the published magnitudes and, when available, re-locations were used.
The local volcano-tectonic earthquakes located at less than 35 km from the volcano are not considered as significant earthquakes because they are sourced to the volcano itself (Prejean and Hill 2018). Nevertheless, we consider them as an indicator of volcanic activity. This distance is quite large but is justified by the fact that we should account for potentially large location errors when the earthquakes are located with only a few (and sometimes distant) seismometers. In total, our catalogue contains 95 earthquakes with d/√S values less than 40, spanning 1992 to 2020. A list of the earthquakes with d/√S < 20 is provided in Table S2, along with the details of the parameters used for the calculations.

Characterization of volcanic activity and unrest at Popocatépetl
The influence of a significant earthquake on activity at Popocatépetl can be characterized by a change in intensity of the eruption. A few parameters are available and consistent over the entire study period to evaluate such changes. CENAP-RED (Centro Nacional por la Prevencion de Desastres, institution in charge of monitoring Popocatépetl) changed the definition of exhalations and explosions on various occasions over the years (e.g. , so that these parameters are not totally consistent through time. In the literature, the most commonly described seismically triggered types of unrest at other volcanoes, and especially at openvents, are seismicity and heat fluxes (Carniel et al. 2003;Harris and Ripepe 2007;Delle Donne et al. 2010;Pritchard et al. 2014;Coppola et al. 2015). Cycles of dome construction and destruction are the typical activity at Popocatépetl. The number of new domes is thus one of the best available parameters to characterize its activity and is a good proxy of the intensity, together with the persistence of ash in the atmosphere as given by the number of VAAC (Volcanic Ash Advisory Center) reports.

Popocatépetl volcano-tectonic earthquakes catalogue
The volcano-tectonic earthquakes (VTs) are events corresponding to discrete fractures often related to locally elevated fluid pressure (Chouet and Matoza 2013) and may thus be indicative of: -Transient local perturbations of the hydrothermal system by dynamic stress changes -Increased pressure conditions around magma chambers The former is likely to also affect magmatic systems far from the critical state (Seropian et al. 2021), explaining why this type of event is more frequently observed to trigger near-surface perturbations. We thus used the VT earthquakes as an indicator of volcanic activity as they are associated with fluid flow (cf. Alparone et al. 2010). Their magnitudes can be converted to energy, and cumulative energy can be calculated. Given that the magnitudes of other volcanic seismic signals (e.g. volcanic tremor) are not constrained at Popocatépetl (cf. Arámbula- Mendoza et al. 2016), we prefer not to use this information. The VT catalogue is a compilation of different sources: (1) The CENAPRED daily (available at http:// www. cenap red. unam. mx/ repor tesVo lcanG obMX/) and annual reports ( To account for the evolution of the seismic network through time, we only selected earthquakes with M d ≥ 2.

MIROVA thermal radiation dataset
MIROVA (Middle Infrared Observation of Volcanic Activity) (see Coppola et al. 2016) is a volcano-dedicated hotspot detection system, based on Middle Infrared (MIR) data acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS). These sensors are mounted on NASA's Terra and Aqua satellites (launched in February 2000 and May 2002, respectively). MODIS images (1 km 2 resolution in the IR bands) are automatically processed by a hybrid algorithm that combines spectral and contextual approaches to detect and quantify the radiant heat associated with volcanic activity (volcanic radiative power (VRP), in watts (Coppola et al. 2016)).
From 24 February 2000 to 31 December 2020, MIROVA analysed 13,420 nighttime images, of which 5511 (about 41%) had at least one alerted pixel at Popocatépetl. Since thermal data is often strongly attenuated by clouds and affected by satellite viewing geometry (Francis et al. 1996;Wooster et al. 1998;Coppola et al. 2016), here we used the weekly average of VRP values as a proxy of the heat flux radiated over discrete time windows of 7 days (cf. Harris et al. 1997;Wright and Flynn 2004). The time integration of the heat flow is then used to calculate the cumulative energy radiated by Popocatépetl throughout the 2000-2020 period (VRE in joules) (Coppola et al. 2016).

OMI SO 2 emissions
OMI (ozone monitoring instrument) is an imaging UV-Vis spectrometer operating since 2004 on the EOS-Aura platform in a sun-synchronous orbit. It allows measurement of trace gases produced by natural and anthropic processes such as O 3 , SO 2 , NO 2 , OClO, and BrO, with a spatial resolution of 12 × 25 km at nadir (Levelt et al. 2006). Monthly SO 2 fluxes are calculated from the monthly SO 2 masses retrieved from OMI data, using the methodology described in Campion (2014). First, all OMI images for any given month are stacked over a matrix with a 0.05° grid interval, with a value of zero attributed to its cells for which no data are available, either due to cloud cover or absence of coverage. Stacks are summed to obtain the monthly cumulative SO 2 matrix, which is then divided by a matrix containing the number of data points acquired over the month. A box is defined around the volcano to isolate the monthly averaged SO 2 anomaly corresponding to gas plumes, and the monthly averaged SO 2 mass is calculated as the sum of the SO 2 column density of every grid element multiplied by its area.
To convert the monthly averaged SO 2 mass into a flux, a calibration, which is specific for a given volcano, is necessary (cf. Campion 2014). This is done by performing a linear regression between a series of monthly masses and its corresponding series of monthly averaged fluxes computed with the traverse methods (cf. Campion et al. 2012). The regression plot for Popocatépetl, obtained over 55 months, for the two series is shown in Supplementary Figure S1. The regression equation is and the R 2 is 0.759. For Popocatépetl, a seasonal effect was found in the monthly mass time series, which was corrected by multiplying the time series with a square sinusoid function so as to smooth the 12-month periodical fluctuation. Gómez-Vazquez et al. (2016) published a catalogue of Popocatépetl's domes between 1994 and 2014 based on seismic signals and available overflight photographs. To extend our dataset over the period 1994-2020, we used the CENAPRED's daily and annual reports. According to Gómez-Vazquez et al. (2016), both catalogues are only slightly different between 2006 and 2010. We thus chose the CENAPRED's catalogue as a reference for the timing of the beginning of the dome extrusion and destruction cycles to keep consistency in our dataset.

VAAC reports
We counted the number of Volcanic Ash Advisory Center (VAAC) reports issued by the Washington office of the National Oceanic and Atmospheric Administration (NOAA). This is a proxy of how often and how long-lasting the ashproducing activity is. The maximum number of daily VAAC reports during 1999-2020 was between eight and nine and corresponded to ash-venting episodes or dome destruction phases. We integrated the number of reports over a monthlong period. Sustained dome extrusion and destruction cycles activity generate several tens of VAAC reports per month, while passive degassing phases (cf. Campion et al. 2018) would only generate a few.

Scaling of the d/√S index using published records
To scale our index, we calculated the d/√S values for all published cases of earthquake-triggering globally (Table S1).  (Manga and Brodsky 2006;Walter and Amelung 2007). The d/√S of 36 calculated for Popocatépetl for this event indicates that the earthquake was too distant to have induced significant static stress changes at Popocatépetl, but dynamic stress would have promoted renewal of volcanic activity. This led us to conservatively extend to 40 the range of d/√S values considered as potential triggers for Popocatépetl. This case may be indicative of how "sensitive" Popocatépetl might be to remote earthquake-triggering through dynamic stresses. Figure 3 displays the relationship between the earthquake Popocatépetl distance and the earthquake magnitude, for earthquakes in our catalogue with an index value d/√S < 40. The M w -distance definitions for triggering earthquakes of previous studies (Linde and Sacks 1998;Marzocchi 2002;Delle Donne et al. 2010;Nishimura 2017;Sawi and Manga 2018) are also represented. About half of the index values of these earthquakes have been calculated using the published surfaces of the rupture plane (Table S2); the other half, as explained in the methods section, was calculated using the scaling relationships of Blaser et al. (2010).

Candidate earthquakes for earthquake-triggering of Popocatépetl
Eighty-two percent of these earthquakes arose from the Mexican subduction, involving distances of less than 850 km to the Popocatépetl and were either interslab thrust or intraslab normal earthquakes. Most of the events were moderate in magnitude (M w 5.5-7.5) and located relatively close to Popocatépetl (250-650 km), which complies with the earthquake characteristics of the definition of Sawi and Manga (2018), i.e. an M6+ earthquake at < 800 km. Larger earthquakes (M w ≥ 7.5) with d/√S < 40 were located at distances of ≥ 350 km and predominantly show inverse mechanisms (Fig. 3). The farthest events were the 2004 Sumatra M w 9.2 megathrust earthquake at a distance of 17,083 km and the Chilean events of Maule (M w 8.8, 2010, distance = 6712 km) and Illapel (M w 8.4, 2015, distance = 6314 km) (Fig. 1). The closest (distance < ~ 300 km) earthquakes have only normal focal mechanisms (Figs. 1 and 3). This is due to the flat-slab geometry of the Cocos plate that is subducting under Popocatépetl ( Fig. 1 in Singh et al. 2018). In our catalogue, these earthquakes had magnitudes up to M w 7.1, as reached by the 19 September 2017 Puebla-Morelos earthquake. This event represents the lowest d√S in our catalogue, having a value of 2.9. The majority of the other nearby normal earthquakes have 10 < d√S < 20.
Only two earthquakes with d/√S < 20 occurred at more than 800 km from Popocatépetl. These are the 2009 M w 7.4 Honduras and 2010 M w 8.8 Maule (Chile) earthquakes (Fig. 1). Although various studies concluded that seismically triggered eruptions are not expected to occur at such distances (Fig. 3) (Linde and Sacks 1998;Nishimura 2017), a number of such distant reactivation have been proposed, especially for very large earthquakes. For example Krenitzyn or Sarychev volcanoes reactivated after the M w 9 Kamchatka earthquake in 1952 (Walter and Amelung, 2007) and Barren Island experienced an eruptive phase after the 2004 Sumatra M w 9.2 earthquake (Walter and Amelung 2007) (Table S1). Note that the 2010 M w 8.8 Chilean earthquake triggered non-volcanic (i.e. tectonic) tremor in Guerrero, Mexico (Zigone et al. 2012), but was not followed by an immediate or drastic change in volcanic activity at Popocatepétl.
None of the earthquakes with the lowest d/√S values at Popocatépetl lie within Nishimura's (2017) range of  Table S1 distance-magnitude for a statistically significant increase in the number of eruptions triggered. The triggering of activity at already active volcanoes, with their magmatic systems at or close to a critical state (such as Popocatépetl since late 1994), would however require less impulse (Hill et al. 2002). This is consistent with the average d/√S values 18 for the published cases of triggering at already active systems in Table S1. This compares with the average value of 4 for triggering new eruptions. Figure 4a presents the temporal distribution of earthquakes with d/√S < 40 between 1992 and 2020. The lower this index, the higher the potential of an earthquake to promote an activity change at Popocatépetl. There is an inhomogeneous temporal distribution for earthquakes with d/√S < 20 (Table S2). We distinguish three main periods: 1992-1999 and 2012-2020 during which seismic activity with low < 10 indices was frequent and 2000-2012 when such earthquakes were absent. Earthquakes with indices 10 < d/√S < 20 were slightly more numerous than indices < 10 between 1992 and 2020 (20 versus 11). They also delineate a wide gap, even if shorter (6.3 versus 12.5 years), between the 2003 M w 7.4 Colima earthquake (d/√S = 15) and a small (M w 5.9) local event in 2009 (d/√S = 13) ( Fig. 1) (Table S2). Earthquakes with indices between 20 < d/√S < 40 were more frequent (75) and present a more regular temporal distribution overall. Nonetheless, their number was also smaller between 2002 and 2006.

Temporal distribution of potential trigger earthquakes at Popocatépetl
There are also three short periods (< 20 months) without earthquakes of d/√S < 40 (Fig. 4a)

Temporal evolution of activity at Popocatépetl
Based on the variations in the parameters considered here as a measure of activity (i.e. VTs, thermal emissions, SO 2 emissions, domes extrusion, ash emissions), we split the activity of Popocatépetl into three periods: Each wit h a roughly similar duration of 8.7 years ± 6 months (Table 1; Fig. 4). Figure 4b gives the temporal distribution of the daily number of VTs with M d ≥ 2 and the cumulative energy. The average daily number of VTs with M d ≥ 2 is 0.37, meaning roughly one event every 3 days. The main features of this time series are the increments in activity that often follow tectonic earthquakes in Mexico. These consisted either of short (lasting a few days) microseismic swarms such as those after the 1999 M w 7 Tehuacan earthquake (De la Cruz-Reyna et al. 2010) or longer (week to month long) seismic sequences, such as those after the 2017 M w 8.2 Tehuantepec and M w 7.1 Puebla-Morelos earthquakes (Figs. 1 and 4). VT swarms are relatively frequent at Popocatépetl and have variable cumulative seismic energies. They often follow tectonic earthquakes with the lowest index values. For example, between 1995 and 2020, six such swarms followed all 10 earthquakes with d/√S < 10. They thus represent a classic feature of seismic triggering at volcanic systems (cf. Hill et al. 1993;Ukawa et al. 2002;Moran et al. 2004;Prejean et al. 2004;West et al. 2005;Pritchard et al. 2014;Prejean and Hill 2018;Butcher et al. 2021). There were fewer VTs between 2003 and 2013 (Fig. 4b), which corresponds to the gap in significant earthquakes defined above (Fig. 4a).

Volcano-tectonic earthquakes
The number of VTs is not completely representative of the volcanic activity because a large number of small earthquakes can release less energy than a few large VTs. It is thus important to also consider VTs energy. Cumulative seismic energy rate reveals four main periods (Fig. 4b).
These were: (1) Between 1996 and 2003, the average slope was 750 MJ/ day. This VT activity corresponds to the beginning of the eruption and started about 1 month after the growth of the first dome. It represents 15% of the total energy released during 1994-2020.
(2) Between 2003 and February 2013, the rate of increase (220 MJ/day) was the lowest for the entire period. This period again corresponds to the seismic gap observed in Fig. 4a and confirms that the small number of earthquakes ( Fig. 4b) was also small in magnitude.  (Fig. 4).

Thermal energy
In Fig. 4c, we give the weekly averaged value of VRP and the cumulative thermal energy. The radiated heat flux at such a dome-forming open-vent system can be associated with: The spikes can be related to intense dome extrusion episodes and/or explosions coinciding with the time of the satellite acquisition as are apparent from the CENAPRED reports and previous description by  using GOES data. We emphasize that these spikes are almost absent from the record between 2003 and 2012, coinciding with the period of low extrusion rate described by Gómez-Vazquez et al. (2016).
The cumulative radiated energy (VRE) trend (Fig. 4c) shares some features with the VT energy curve (Fig. 4b) and helps to define the three periods of activity. However, the linear behaviour in Fig. 4c is consistent with overall stability, with perturbations always returning to the general trend (on average ~ 4 × 10 8 MJ/year). This is typical of long-lasting open-vent systems (Wright et al. 2015;Barrière et al. 2017;Laiolo et al. 2018). Figure 4d gives the monthly SO 2 fluxes for Popocatépetl for 2004-2020. SO 2 fluxes are commonly associated to the magma supply rates above the level of exsolution (Shinohara 2008). At Popocatépetl, SO 2 fluxes have been shown to increase during episodes of dome extrusion, reflecting the supply of undegassed magma to sustain the ascent of the magma column (Taquet et al. 2017(Taquet et al. , 2019Campion et al. 2018). A period of sustained high SO 2 emissions starting in September 2011, peaking in May 2012, and decaying through the end of 2014 apparent in Fig. 4d is similar to the trend in the thermal flux (Fig. 4c). This period corresponds to the magma recharge episode proposed by Espinasa-Pereña (2012) and Taquet et al. (2019). As shown by Campion et al. (2018), most of the magma involved in the recharge degassed at depth with only a small amount being extruded. This is similar to the intrusion that fuelled the first phase of the eruption in 1994 (Witter et al. 2005;Roberge et al. 2009). In both cases, less than 2% of the intruded volume reached the surface. SO 2 emissions during PII and most of PIII are stable at values < 2000 tons/day (Fig. 4d). This shows that the magma supply from depth was either stable in the long-term or inaccessible through the SO 2 fluxes, because the injection depth lies below its exsolution level. A few outliers exceed 2000 tons/day and correspond to extrusion of particularly voluminous domes as occurred early 2015 or during February to March 2019. Figure 4e gives the annual number of cycles of dome extrusion and destruction and the number of domes extruded throughout the eruption. The first dome grew in March 1996 (Sheridan et al. 2001;Gómez-Vazquez et al. 2016), about a month after the 25 February 1996 M w 7.1 of Pinotepa Nacional earthquake (d/√S = 10) (Fig. 1). Cycles of dome extrusion and destruction then mainly took place during PI and PIII, as opposed to a period with low extrusion rates, i.e. PII (cf. Gómez-Vazquez et al. 2016), in agreement with the absence of spikes in thermal emissions described above. Seventy-seven domes out of the 85 counted by CENAP-RED (> 90%) were extruded during these two periods. The remaining 10% were extruded during 2003-2011.

Cycles of dome extrusion and destruction
In both PI and PIII, the temporal distribution of dome extrusion and destruction cycles is irregular (Fig. 4e) with regimes of higher rates of dome emplacement recognized by Gómez-Vazquez et al. (2016). We emphasize that such short-term regimes of high dome emplacement rates can be considered as pulses of dome extrusion and will refer to them as such. Based on Fig. 4e  Each of these pulses is ~ 1.5-2.5 years long (cf. Gómez-Vazquez et al. 2016) and has a rate of dome emplacement higher than the average over the period in which it takes place (cf. Table 1). Another less intense pulse is observed in 2017-2019, but the emplacement rate (4.8 domes per year) is lower than that of the PIII. These pulses peak during the following paroxysmal eruptive sequences: Fig. 4 Time series of a earthquakes with d/√S ≤ 40; b daily number (black lines) and cumulative energy (green curve) of volcano-tectonic earthquakes (M d > 2); c weekly average (in blue), 6-month moving average (orange curve) of MIROVA VRP and VRE (cyan curve), the black dashed line is the slope of the whole dataset; d OMI-derived SO 2 flux and cumulative seasonally corrected SO 2 masses (gold curve) emitted since 2004; e annual (grey bars) and cumulative (red line) number of domes extruded; and f monthly number of VAAC reports. Vertical lines: the black dashed line marks eruption onset, the red dotted lines mark the timing of earthquakes with d/√S < 10 in each time series, and the grey bold lines split the time series into the three periods defined here (PI, PII, PIII)

◂
(1) The 30 June 1997 explosion (Sheridan et al. 2001) (2) The December 2000 cycle of dome construction and destruction (Sheridan et al. 2001;Macías and Siebe 2005), which remains the most voluminous of the eruption at the time of writing (May 2022) (3) The May to July 2013 sequence (cf. Campion et al. 2018) that terminated with the first crater collapse episode (

Main periods of volcanic activity at Popocatépetl
Popocatépetl underwent intense volcanic activity during periods I (1994-2003) and III (2011III ( -2020. These two periods were characterized by many dome extrusion and destruction cycles and frequent VTs (Table 1, Fig. 4). The system was less active during period II (2003-2011), which included a 1.5-year eruptive pause.
The temporal distribution of tectonic earthquakes likely to promote volcanic activity at Popocatépetl follows a similar trend. During periods I and III, six and seven tectonic earthquakes with a 10 < d/√S < 20 occurred (Fig. 4a). Instead, during period II, only three of them were recorded. If we consider the < 10 class, then 4-6 events occurred during periods I and III and none during period II (Table 1).
The periods of intense volcanic activity (PI and PIII) are characterized by a high rate of dome extrusion (26 and 51 dome extrusion and destruction cycles, respectively, i.e. 3-5.5 per year). This contrasts with the low rate (1 per year) observed during period II, when only eight domes were reported. Such a drop in activity is consistent with the low number of VAAC reports. There was an average of 32 reports per year during the PII, but four-to-seven Table 1 Volcanic and tectonic activity during the three main periods of volcanic activity between 1994 and 2020. The number of significant regional and worldwide tectonic earthquakes (with d/√S < 20) occurring during these periods is listed. Italic values are given for the period beginning after the date specified in the description. * average daily number of VTs  Fig. 4e). VT events also show a twofold decrease during PII (decreasing from 0.45 to 0.22 per day). In terms of energy, the cumulative energy released by VTs during PII represents just 9% of the total energy liberated during the eruption. Radiated power shows the same variations (Fig. 4c). Whereas PI and PIII exhibit similar annual values of 4.7 and 4.4 × 10 8 MJ/yr, PII has a value of 2.6 × 10 8 MJ/yr. In terms of VRP, the three periods have values of 14 MW, 8 MW, and 14 MW, respectively. The background SO 2 emissions have been very similar over 2004-2020, and the time series is dominated by the strong degassing pattern related to the 2011 recharge and 2012-2013 eruptive crisis.

Triggering potential index
In an oceanic subduction context such as in Mexico, hundreds of tectonic outer-rise, interplate, intraplate, and crustal earthquakes occur every day with magnitudes M w ≥ 3. It is thus of prime importance to identify which of them are likely to trigger any kind of change in volcanic activity. This is the purpose of our new index, d/√S.
In seismology, the near-field corresponds to distances of the order of the rupture length (Aki and Richards 1980;Vidale et al. 1995), i.e. d/√S values < 10. It applies to distances where the near-field waves and the static grounddisplacement are significant with respect to the far-field waves. However, the value of ~ 10 is not well defined and cannot be considered as a unique and absolute value (cf. Cummins 1997). This is because it will depend not only on the relative distance to the earthquake, but also on many other parameters such as the distribution of the slip over the fault plane, the focal mechanism, and the position of the volcano with respect to the fault plane. This is why we here determine the range of index values for published reactivation cases, independently from static or dynamic stress changes.
The distribution of our d/√S values obtained for published cases associated with eruptive activity triggering (Fig. 2) shows that for most of the published reactivation cases, triggering occurs at the lowest values (i.e. > 90% of cases have d/√S indices < 20). Nevertheless, the distribution also shows cases of high index values (> 20) corresponding to reactivation by distant (5000-17,000 km), strong (M w > 8), or regional (< 500 km) low-magnitude (M w < 6.5) earthquakes. These can be explained by dynamic stress triggering Delle Donne et al. 2010). Low index (< 10) values, however, are not sensitive to the type of triggering since: (1) Strong transient dynamic stress changes can accompany the passage of seismic waves in the Near-Field ).
(2) The eruption-triggering mechanism will be different depending on the earthquake scenario and volcano type (cf. Seropian et al. 2021).
Our index cannot inform either on the intensity of the static stress change at the volcano or whether it is increasing or decreasing. Such an assessment requires accurate knowledge of the earthquake and the volcanic system with which it is interacting (Nostro et al. 1998;Walter 2007;Chesley et al. 2012;Bonali et al. 2013Bonali et al. , 2015Bonali 2013). Depending on the state and level of criticality of the magmatic system before a significant earthquake, the amount of seismic stress to be transferred to trigger an eruption can be highly variable (Manga and Brodsky 2006;Bonali et al. 2015;Kriswati et al. 2019;Hamling and Kilgour 2020) so that defining an absolute and globally applicable threshold cannot be achieved.
Our index allows objective discrimination and sorting of large earthquake catalogues according to their potential to promote an eruptive response at a recipient volcano, without using arbitrary magnitude-distance thresholds. In the case of Popocatépetl, we reduced the catalogue from 734 (187 613 with the complete catalogue of small local earthquakes) to 31 earthquakes with a d/√S value < 20. This reflects the benefit of our approach in building a comprehensive and usable database.
An illustration of the importance of this index to the study of system triggering is the case of the 20 December 2001 explosion which was part of the destruction phase of Popocatépetl's dome #17 (Espinasa-Pereña 2012). De la Cruz-Reyna et al. (2010) argued that this event was triggered by a static stress change, generated by a succession of lowmagnitude (between 3.5 and 4.8) regional earthquakes at relatively large distances (between 280 and 575 km) from Popocatépetl. Our calculated d/√S for these earthquakes are greater than 193 and up to 575. These values are too high for the static stress fields of such earthquakes to significantly affect Popocatépetl. Based on our catalogue, such earthquakes occur on a weekly-monthly basis, and, counting only contiguous events on equally short time periods, 10 similar earthquake sequences with an equal or lower average d/√S value (387 for the sequence of De la Cruz-Reyna et al. (2010)) occurred on at least five occasions in the same month (December 2001). From the point of view of our index, this earthquake sequence is therefore very common, and its relationship with this explosion is probably casual.
The main outcome of our new methodology analysis was to thus properly identify those earthquakes with the potential to affect volcanic activity and the need for large catalogues for both seismic and eruptive activity. Figure 4a shows that earthquakes with index values of 20 ≤ d/√S ≤ 40 occurred continuously throughout the 26-year-long dataset, even during the period of low levels of activity between 2004 and 2011. This class of earthquakes is thus less likely to promote volcanic activity at Popocatépetl. Instead, the time series of earthquakes with d/√S of < 20 is irregular and in particular the d/√S < 10 class with a decade long gap between 1999 (Oaxaca) (Fig. 1, Table S2) and 2012 (Ometepec). However, before 1999 and after 2012, they occur rather frequently (Fig. 4a). These two periods coincide with periods when the activity is more intense. Our index thus suggests that these are "significant earthquakes" with the potential to promote activity at Popocatépetl. The d/√S value of < 10 (− 20) value is not a threshold and shall not be used as such. It also has no predictive potential because volcanoes have to be poised to erupt to be triggered by earthquakes (Barrientos 1994;Marzocchi et al. 2002;Manga and Brodsky 2006;Eggert and Walter 2009;Prejean and Haney 2014). Instead, we here stress this value for screening earthquake databases and identifying events with the highest potential to promote volcanic activity at a recipient volcano, as a preliminary step for more detailed studies of triggered volcanic activity.
The gap in significant seismicity between 2000 and 2012 implies that the regional (< 650 km) and high magnitude (M w ≥ 6) events associated with subduction may be missing. Around this gap, such events account for all d/√S < 10 events and ~ 75% of the 10 < d/√S < 20 earthquakes. However, for such magnitudes, seismic catalogues are complete (Sawires et al. 2019;Ammon et al. 2021), so that the lack of large earthquakes is real during the "gap". An analysis of the SSN "large earthquakes catalogue" (with M w > 5.5; http:// www2. ssn. unam. mx: 8080/ sismos-fuert es/) over our study period reveals that between 2004 and 2009, less earthquakes occurred and average magnitudes were lower and that during 2009-2010, large earthquakes (M w > 5.5) were on average farther (> 400 km more) from Popocatépetl than in other years. Explanations of this gap are beyond the scope of this paper but would merit further study. We herein focus on the interaction between earthquakes and activity at Popocatépetl, and the relations revealed through application of our index.

Modulation of Popocatépetl's activity by regional and worldwide earthquakes
Time series of all studied parameters (Fig. 4) delineate three main periods, each defined by characteristic behaviour. The relationship with our triggering potential index shows that Popocatépetl undergoes enhanced activity when significant tectonic earthquakes are frequent, as during 1995-2003 and 2012-2020. Such enhanced activity is characterized by more frequent extrusive and explosive activity, as well as an increase in heat flux (Fig. 4). More intense microseismicity is also observed, the seismic energy release during PI and PIII being a factor of five higher than during PII, likely indicating increased pressure conditions in the magma chamber. Such increased activity may start shortly after a significant earthquake, especially the swarms of VTs (De la Cruz-Reyna et al. 2010), and can sometimes continue for months, as after the September 2017 earthquakes (Fig. 4b).
The four main pulses of dome extrusion and destruction occurred shortly after: 1  Table S2), including the event with the lowest d/√S of the whole study period (Fig. 4a). This apparent contradiction is detailed in the next section. The lag between the earthquakes and the paroxysm of each pulse is on average 1.3 ± 0.3 years. This is consistent with the statistically significant increase in number of eruptions at "sensitive volcanoes" (mainly represented by subduction andesitic stratovolcanoes) within 2 years of a M6+ earthquake observed by Sawi and Manga (2018).
During the period preceding the onset of the Popocatépetl eruption on 21 December 1994, hydrothermal activity progressively increased from late 1991 (Gerlach et al. 1997) and VT seismicity from mid-1993 (De la ). The eruption then initiated just 10 days after the M w 6.4 Zihuatanejo earthquake (d/√S = 20). This earthquake is the last of a series of 12 with d/√S < 40 that took place during the years 1993-1994, including five events with a value < 20 (Table S2). This sequence represents the maximum density of significant regional earthquakes in our whole dataset with the potential to promote volcanic activity at Popocatépetl (Fig. 4a). We infer that such repetitive seismicity might have favoured pressure buildup and initiation of the eruption, and ending ~ 70 years of dormancy.
The corollary to intense volcanic activity accompanying frequent significant seismicity is the low activity of the volcano during periods of seismic calm. Following the second pulse of cycles of dome construction and destruction (2001)(2002)(2003), the activity almost ceased for 1.5 years, where thermal emissions dropped by almost an order of magnitude (Fig. 4c) and the alert level was decreased (Espinasa-Pereña 2012; De la Cruz-Reyna et al. 2017). This corresponded to the beginning of the 6.33-year-long gap without significant earthquakes (d/√S < 20) observed during PII. Dome growth activity resumed in April 2005 (Fig. 4e), about 4 months after the 2004 Sumatra earthquake, but Popocatépetl maintained relatively low levels of activity until 2009, although with a few cycles of dome construction and destruction (five in four years). Following another year (28 April 2008 to 27 April 2009) without earthquakes with index value < 40, a 7-month-long episode with low thermal fluxes was recorded (Fig. 4c). Two very different earthquakes with d/√S < 20 within a week of each other ended this year. These were a regional (d ~ 100 km) low-magnitude (M w 5.6) event and the 2009 M w 7.4 Honduras earthquake (Fig. 1, Table S2). The activity and thermal flux returned to their previous levels by November 2009 (Fig. 4c). The activity at Popocatépetl remained at low levels until late 2011 when the magma recharge episode proposed by Espinasa-Pereña (2012)  Popocatépetl had already passed through a period of low-level volcanic activity during 1998-2000, which followed ~ 2 years without earthquakes with index value < 40 (Fig. 4). In this timespan, dome extrusion activity was reduced, and in 1999, there was no extrusive activity. Following the 1999 (d/√S value of 5.6 and 6.1) earthquakes (Fig. 4, Table S2), activity progressively returned (De la Cruz-Reyna and Tilling 2008) to remain relatively high until 2003 (cf. Gómez-Vazquez et al. 2016) and included some of the most intense extrusive phase, such as the December 2000 dome (Sheridan et al. 2001;Macías and Siebe 2005;Gómez-Vazquez et al. 2016). This shows that following a decline in the recurrence of significant earthquakes, the activity at Popocatépetl declines and can remain low for years, until a new impulse promotes enhanced activity.
Furthermore, during periods of reduced extrusive activity, the system apparently conserves its potential to erupt. Although PII shows lower extrusive activity (with eight domes extruded between 2005 and 2010), the thermal and SO 2 fluxes during most of this period were similar (lower by a factor of two and four for the thermal and SO 2 fluxes, respectively) to the average fluxes during 2012-2016, i.e. after the recharge episode of 2011 ( Fig. 4c and d) (Espinasa-Pereña 2012; Carn et al. 2017;Taquet et al. 2019). The fact that both the thermal and degassing fluxes remained comparable during periods with such contrasted extrusion rates suggests that during PII the system: -Maintained its open-vent behaviour -Had an internal state comparable with that of most of PIII, at least from the thermal and degassing point of view -Could have been fueled with equivalent magma supply rates The lack of cycles of dome extrusion and destruction during PII thus seems to indicate that seismic triggering could be a viable mechanism to modulate extrusion at Popocatépetl. Among such mechanisms, increased magma supply rates and intrusion events usually hold a major role, contributing to increase the pressure in the magma chamber and/or revitalize the convection, allowing maintaining the magmatic system in a critical state.
At . This is around 7 months after both the thermal and SO 2 fluxes started to increase ( Fig. 4c and d), i.e. the end of PII to beginning of PIII. This activity reached its maximum with the May to July 2013 dome extrusion events (Donovan et al. 2018) and ended in July 2013 with a crater collapse episode (Taquet et al. 2019). Following this collapse, the SO 2 emission rate decreased until October 2014 (Fig. 4d) (Taquet et al. 2019) and was only a factor of 1.5 higher than during the period 2005-2010 (Fig. 4d). This suggests that sustained high magma supply rates cannot account for the extrusive activity. Instead, the 18 April 2014 M w 7.3 of Papanoa earthquake (d/√S = 6) occurred about 1 month before the extrusion of domes #50-52 during the last 3 months of this period.
We thus find multiple lines of evidence for the role of significant earthquakes in promoting eruptive activity at an open-vent system, in this case Popocatépetl. Although magma supply from depth maintains the system in a persistently eruptive state, large earthquakes have a role in driving the system towards extrusion and/or facilitating the transfer of magma towards the surface (cf. Bonali 2013; Seropian et al. 2021).

Modification of eruptive style after the 2017-2018 earthquakes
Sustained ash productivity began immediately after the 19 September 2017 M w 7.1 of Puebla-Morelos earthquake, as illustrated by the increased number of emitted VAAC reports following this event (Fig. 4f). An increased number of exhalations was also detected by CENAPRED so that the increase documented by the VAAC cannot be ascribed to increased detection capacities since the end of 2016 and with the launch of GOES16 and GOES17 (Engwell et al. 2021). The number of reports distributed after the earthquake (2.8 per day) was greater than during the 2015-2016 pulse of extrusion (0.8 per day). This change was accompanied by a decrease in the number of dome extrusions between 2018 and 2021, there having been 11 domes extruded (3.3 per year), as opposed to 17 during 2015-2016 period (8.5 per year). This suggests that the 2017-2018 earthquakes (d/√S values < 8.4) ( Table S2) may have modified the eruptive style at Popocatépetl, pushing it towards a more explosive regime. This transition involved the beginning of low intensity but long-duration ash emissions, as well as medium to strong explosions unrelated to the construction of domes (Caballero and Valderrama 2021). Immediately following the Puebla-Morelos earthquake, several short episodes (sometimes < 1 h) of ash-venting or "Strombolian" activity occurred. Such events had previously been associated with episodes of dome extrusion but after September 2017 occurred without extrusion .
This modification of eruptive style could have been driven by several different processes. For example, an increased volatile content in the ascending magma column could result in higher levels of explosivity (Parfitt and Wilson 1995;Cashman 2004). However, such an increase is not supported by the SO 2 emissions (Fig. 4d) that have instead decreased since 2014. Alternately, increased ash-producing activity could be due to increased magma viscosity, which can result either from its dehydration or a drop in its temperature (Cashman 2004;Dingwell 2006). We thus suggest that propagation of seismic waves or the stress field change following the 2017-2018 earthquakes favoured the exsolution of volatiles (cf. Brodsky et al. 1998;Manga and Brodsky 2006;Davis et al. 2007). H 2 O in particular would have been the major dissolved volatile component, whose exsolution would also have increased viscosity (cf. Giordano and Dingwell 2003). This would have favoured crossing of a viscosity threshold and the subsequent switch from extrusive to explosive regime (cf. Slezin 2003).

Conclusion
The influence of earthquakes on volcanic activity has long been debated and linked with changes in gas flux, seismicity, heat flux, and/or eruption rate. However, such cases are usually based on a single-case, a single-event, or over a limited time scale, thus not being statistically robust. In this study, we report how tectonic earthquakes may have contributed to both the onset and changes in activity at an open-vent system using a database spanning three decades, i.e. at Popocatépetl for the period 1994-2021.
We introduce a new index, which we scale with published cases of triggering. Our index allows objective discrimination of earthquakes with the potential to trigger an eruptive response at a recipient volcano, without recourse to an arbitrary distance-magnitude threshold. As we demonstrate, such a simple index is useful in volcanic hazard assessment, allowing potential changes in activity to be forecast.
At Popocatépetl, our index allowed a catalogue of subduction earthquakes comprising 187 613 events to be reduced to just 31 significant events, i.e. events that have the potential to trigger a change in activity and with index (d/√S) values < 20. We show that earthquakes with the lowest (< 10) index values are absent during the period with the lowest extrusion rates. Instead, the four pulses of extrusive activity all reached their peak 1.3 ± 0.3 years after an earthquake with an index of < 10. At open-system volcanoes, times for delayed responses on the order of months to years can raise doubts on the causal link between the earthquake trigger and the change in volcanic activity. Our study, based on a 26-year-long time series, not only provides a statistically robust base for our observations, but more importantly captures an episode of seismic and volcanic quiescence. This 12-year-long gap in "significant seismicity" is indeed coincident with almost 9 years of low extrusive and volcanic activity during which the magmatic system though maintained its open-vent character (with comparable thermal and gas fluxes). Instead, and although magma supply rates were comparable, extrusive activity increased during periods of frequent "significant tectonic earthquakes". This highlights the influence of "significant earthquakes" at Popocatépetl in modulating its eruptive activity. We thus propose that the recurrence of significant earthquakes over years partly controls and contributes sustaining the extrusion at the volcano during this long eruption.
In short, we find that although the magma supply from depth maintains the system in a persistently eruptive state, changes (quasi-static or dynamic) of the volcano stress field imposed by "significant earthquakes" (low d/√S index values) modulate its output rate.
Funding DL is funded by the PAPIIT project number # IN106019 and SV by the PAPIIT project number #IA102221 and RC by the CONA-CYT-SEP project number A1-S-30127.

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