Repetitive small seismicity coupled with rainfall can
trigger large slope instabilities on metastable volcanic
edifices

Virginie Durand, Anne Mangeney, Pascal Bernard, Xiaoping Jia, Fabian

Bonilla, Claudio Satriano, Jean-Marie Saurel, El Madani Aissaoui, Aline

Peltier, Valérie Ferrazzini, et al.

► To cite this version:

Virginie Durand, Anne Mangeney, Pascal Bernard, Xiaoping Jia, Fabian Bonilla, et al. Repetitive small seismicity coupled with rainfall can trigger large slope instabilities on metastable volcanic edifices. Communications Earth & Environment, 2023, 4, 10.1038/s43247-023-00996-y. http://dx.doi.org/10.1038/s43247-023-00996-y

HAL Id: hal-04401982

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

Submitted on 18 Jan 2024

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

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

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

# Repetitive small seismicity coupled with rainfall can trigger large slope instabilities on metastable volcanic edifices

Virginie Durand<sup>1,2✉</sup>, Anne Mangeney<sup>2</sup>, Pascal Bernard<sup>2</sup>, Xiaoping Jia<sup>3</sup>, Fabian Bonilla<sup>1,2,4</sup>, Claudio Satriano<sup>1,2</sup>, Jean-Marie Saurel<sup>1,2</sup>, El Madani Aissaoui<sup>2</sup>, Aline Peltier<sup>1,2,5</sup>, Valérie Ferrazzini<sup>1,2,5</sup>, Philippe Kowalski<sup>2,5</sup>, Frédéric Lauret<sup>2,5</sup>, Christophe Brunet<sup>2,5</sup> & Clément Hibert<sup>1,6</sup>

Quantifying the effect of external forcings like seismicity or rain on slope destabilization is a long-standing and challenging issue. To investigate the respective roles of these forcings, we analyze an unprecedented 10-year long catalog of rockfalls occurring in the crater of the Piton de la Fournaise volcano (La Reunion Island), using statistical tools originally developed for earthquakes. Our analysis reveals the predominant effect of low amplitude repetitive seismicity in the triggering of rockfalls located at a few kilometers from the source, due to progressive damaging of the slope. Moreover, we show that the efficiency and time-delay of this dynamic triggering is controlled by the stability state of the slope, i.e. its closeness to the failure, as observed with lab-experiments on metastable granular slopes. Our results show the need to account for long-term swarm-type seismic activity that can affect the stability of geological structures like slopes and faults, but also buildings.

<footnote n="1">¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, Sophia-Antipolis, France. ²Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, Paris F-75005, France. ³Institut Langevin, ESPCI Paris Tech, CNRS UMR, 7587 Paris, France. ⁴Geotechnical Engineering, Environment, Natural Hazards and Earth Sciences Department, Université Gustave Eiffel, Marne-La-Vallée 77447, France. ⁵Observatoire Volcanologique du Piton de la Fournaise, Institut de Physique du Globe de Paris, La Plaine des Cafres, F-97418 La Reunion, France. ⁶Institut Terre et Environnement de Strasbourg, ITE, CNRS UMR 7063, University of Strasbourg/EOST, Strasbourg, France. ✉email: virginie.durand@geoazur.unice.fr</footnote>
Unstable slopes in mountainous, coastal or volcanic area are subject to different forcings, acting at different time-scales, that accumulate and may lead the slopes to fall down and generate gravitational flows such as landslides or rockfalls<sup>1-3</sup> (Fig. 1). Slope instability essentially depends on the balance between external forcings and the resistance force due to cohesion and friction within the material prone to fall. When the rock strength (cohesion and friction) decreases and the system gets close to the conditions for mechanical instability, a state called metastable, any small perturbation can trigger macroscopic failure.

Gravity is the primary permanent and therefore long-term force that may prime a metastable fractured slope for failure by progressively initiating and propagating cracks within the system and thus decreasing its cohesion and internal friction. This driving process (imposed force) can be compared to the loading of faults by relative plate motion (imposed deformation inducing increasing stress) that leads the metastable fault to yield stress and macroslip. Intermediate-term (a few years) to long-term weath-ering processes such as dissolution and alteration notably by rain may also induce the rock walls to progressively lose their stability by reducing rock cohesion or creating unstable reliefs and therefore changing the stress distribution. On the contrary, che-mical effects may also help heal the system and consolidate the material on a few years time-scale as observed for example during post-seismic recovery<sup>4,5</sup>.

At shorter time scales, from seconds to weeks, meteorological (e.g. rainfall, storms) and seismic activity induce changes that may trigger failure in a metastable system (Fig. 1). Rainfall increases local pore pressure for a duration that depends on the drainage capacity of the slope (e.g. porosity, fracture density). Strong winds as well as shaking from seismic events may apply stress fluctuations to the system. Furthermore, shaking may decrease the rock strength by slightly expanding the material (increasing the joint space between grains) and/or developing cracks thus reducing friction and cohesion.

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Sketch summarizing the mechanisms destabilizing the slopes in the case of the Piton de la Fournaise volcano. This sketch highlights how the different mechanisms impact the stability of the slopes. μ is the friction coefficient, p is the pore pressure, c is the cohesion of the material. The dashed line on the volcano profile (right figure) indicates the crater geometry before the 2007 crater floor collapse. Light yellow lines note the schematic overall fracturation of the massif. Thick yellow lines highlight larger fractures. Red lines show the development of fractures due to seismicity.</figcaption></figure>

The resulting triggering of instability, instantaneous or delayed, may depend on the forcings intensity and duration, but also on the stability of the slope, i.e. the distance to failure. A very unstable system having for example experienced a catastrophic event (e.g. caldera collapse) will spontaneously destabilize. In turn, falling down of the unstable material may relax and stabilize the slopes (by purge or annealing) making them less sensitive to triggering processes. On the other hand, in a metastable system, failure may be triggered more or less rapidly by the decrease of the resistance strength due to local changes in the yield strength distribution induced by external forcings. Essentially, triggering effects may locally shorten the time to failure, a process that strongly depends on the slope characteristics.

In the case of the triggering of slope instability by rainfalls, delays of 1 h to 2 days are observed, with rockfall activity potentially lasting several days<sup>2,6-8</sup>. Earthquakes, on the other hand, have been shown to trigger slope instabilities with delays up to several days, with the activity potentially lasting several months, by producing a permanent decrease in the resistance strength or yield stress of faults<sup>9,10</sup>. Finally, rainfall and seismicity can have a joint effect and lead to slope collapses several days after the seismic activity<sup>1,11</sup>. In volcanic settings, the slow deformation of the edifice due to the rising of magma is also expected to play a role by changing the applied stresses as well as decreasing rock cohesion or propagating cracks. The associated time-scale is not clearly defined and may be longer than the time-scale for seismicity as suggested for example at the Piton de la Fournaise volcano<sup>12,13</sup>.

Understanding how the different forcings interplay and impact slopes (as summarized in Fig. 1) is important to better assess landslide and rockfall hazards, in particular in earthquake prone regions. In the study presented here, we will focus on the factors triggering rockfalls at a time-scale of several days to weeks.

No clear triggering threshold is proposed in the literature for the rainfall effect. For earthquakes, it is commonly assumed that $M < 4$ earthquakes are not potential triggers<sup>3,14</sup>. However, some observations suggest that smaller single events combined with rainfall, as well as swarms of smaller earthquakes that cause repetitive shaking, may trigger slope instabilities<sup>1,12,14,15</sup>. The small number of statistically significant studies focusing on the influence of low magnitude seismicity ($M < 4$) on the triggering of slope instabilities prevent us from understanding well the respective role of forcings on slope instability<sup>3,11,15</sup>.

To quantify the potential of external forcings to trigger mass destabilization, it is necessary to consider and analyze the activity of the numerous existing small events. The difficulty is that slope instabilities are mainly identified using field observations and aerial or satellite images<sup>8,16</sup>. This leads to high uncertainties in the estimation of their occurrence time and incomplete detection of small and moderate events. However, advanced monitoring techniques have been developed using continuous imaging instruments<sup>2,16</sup> or the seismic waves generated by the slope instabilities<sup>17-19</sup> to detect and locate them and estimate their volumes.

Here, we analyze the temporal evolution of rockfalls occurring in the crater of the Piton de la Fournaise volcano (La Reunion Island) using 10 years (2007–2017) of continuous seismic data. The seismic network<sup>20</sup> used to monitor the volcanic seismicity (Fig. 2) allows us to precisely detect, locate, and estimate the volumes of the rockfalls from the waveform and energy of the generated seismic signal, leading to a 10-year high-precision catalog of rockfalls<sup>21</sup>. To get insight into the response of the slopes of the crater to external forcings, i.e. rain and seismic shaking as well as deformation due to the rising of magma, we analyze this catalog using statistical tools and methods originally developed for earthquake analysis.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Eruptions and monitoring stations of the Piton de la Fournaise Volcano used in this study. Digital elevation model of the Piton de la Fournaise volcano and location of the seismic, Global Navigation Satellite System (GNSS), and rain gauge stations used in this study, along with locations of the eruptions that occurred from 2008 to 2017: (1) 21 September to 2 October 2008, (2) 27 and 28 November 2008, (3) 14 December 2008 to 4 February 2009, (4) 5 and 6 November 2009, (5) 14 and 15 December 2009, (6) 2-12 January 2010, (7) 14-31 October 2010, (8) 9-10 December 2010, (9) 20 and 21 June 2014, (10) 4-15 February 2015, (11) 17-30 May 2015, (12) 31 July to 2 August 2015, (13) 24 August to 31 October 2015, (14) 26 and 27 May 2016, (15) 11-18 September 2016, (16) 31 January to 27 February 2017, and (17) 13 July to 27 August 2017. The lines next to the eruption numbers show the locations of the eruptive fissures. Insert: Piton de la Fournaise on La Reunion Island.</figcaption></figure>

The Piton de la Fournaise volcano is a unique site to study the dynamics of rockfall activity, including its response to external factors. Indeed, this volcano is located in a tropical region and is consequently exposed to heavy rains (possibly more than 800 mm in a day). It is also one of the most active basaltic volcanoes in the world, with an eruption every 8 months on the average²². All the eruptions are accompanied by numerous small seismic events (volcano-tectonic type), with an average magnitude of 0.6 and a maximum magnitude of 3.8 for an event that occurred in 2007. During a major eruption, the Dolomieu crater floor collapsed in April 2007²³, exposing steep rock walls, 320 meters high, prone to a continuous rockfall activity involving a broad range of volumes from m³-scale individual rocks to 100,000 m³ granular flows. Furthermore, the summit of the volcano is densely instrumented²⁰,²⁴, with seismic and GNSS stations (Global Navigation Satellite System), rain gauges (Fig. 2), cameras, and regular photogrammetric campaigns¹². The 2007-2017 time-span includes periods with different forcings (Fig. 3). First, directly after the crater floor collapse, the seismic activity decreases in the first months (Period 1, 2007/05/02 - 2009/08/31) while exhibiting long-duration high activity episodes before subsequent eruptions (Period 2, 2009/09/01 - 2011/12/31). This is followed by a 2.5-year period with no eruptive activity (Period 3, 2012/01/01 - 2014/05/31) and consequently low seismic activity, during which the main forcing is rain, making the analysis of its impact easier. Finally, 2014 corresponds to the renewal of eruptive activity, with short and burst-like seismic activity preceding the eruptions (Period 4, 2014/06/01 - 2017/12/31).

We will quantify here the spatio-temporal change of rockfall triggering in response to these short term seismic and rainfall forcings (time-scale of days) in a context of long-term evolution (over years) of this volcano edifice, and show that the slow deformation of the Dolomieu crater induced by the rising of the magma has a much smaller impact at this short time scale.

## Results

The temporal evolution of the rockfalls shows the global relaxation of the crater. The temporal evolution of the cumulative number and the cumulative volume of rockfalls, that begins after the crater floor collapse, is shown in Fig. 3a (red and orange curves, respectively). As already highlighted by Hibert et al. 2017¹³, we clearly see intense rockfall activity during the first 2 months after the 2007 crater floor collapse, followed by a decrease in the activity towards a more stable state during the following years. We can nicely fit an Omori-type law to the first 2.5 years of rockfall activity (red dashed line in Fig. 3a):

$$n(t) = \frac{K}{(t + c)^p} \quad (1)$$

with $n(t)$ the rate of events, $t$ the time after the main event (here the collapse), and constants $K$, $p$, and $c$, being the productivity of events, the power law exponent, and the time delay before the onset of the power-law event decay rate, respectively. This law is classically used to describe the evolution of the aftershock activity following the occurrence of an earthquake, characterizing the relaxation of the system. The fit to the rockfall data is performed by varying the parameters $K$, $c$ and $p$. The $c$ and $K$ values used here are equal to 1.9 days and 197.4, respectively. The $p$ value we use to fit the rockfall data is equal to 0.8, which is in the range of the values observed for earthquake aftershocks, generally close to ¹²⁵. It describes the rapidity of the aftershock rate decay. This similarity suggests that rockfalls happening during this period are mainly the “aftershocks” of the collapse, generated by the relaxation of the slopes. It also confirms the observations of a relaxation phase lasting several years after the collapse (comparable to the one observed after a large earthquake²⁶), making it possible to quantify the duration of this phase at 2.5 years, a period during which the relaxation law fits the observations.

At the end of 2009, the rockfall activity deviated from the Omori-law (Fig. 3a, black arrow), suggesting it is no longer controlled by the crater floor collapse. Comparing the temporal evolution of the rockfalls to the local seismicity (Fig. 3a, b), we note an increase in the rockfall activity correlated to an increase of the seismic activity before and during three eruptions (indicated by the black arrows in Fig. 3a, b).

**Spatio-temporal distribution of rockfalls on unstable to more stable slopes.** Right after the crater floor collapse, the whole crater is unstable, and large and moderate rockfalls ($V > 5000 m^3$) are homogeneously distributed over all the slopes exposed by the collapse (Fig. 4). Even though the spatial distribution of small rockfalls is homogeneous over the whole time-span investigated here (Supplementary Fig. S1), spatial concentrations of larger rockfalls is observed after Period 1 when the crater walls become more stable. In particular, in 2011, the largest rockfalls were clustered on the steepest (around 70 degrees) north-west slopes (Fig. 4b). According to Derrien et al. 2019²⁷, this area exhibited a drastic fracture opening between 2007 and 2009. After several months, during which the unstable slopes were purged by the occurrence of many rockfalls (~500 events with $V > 100 m^3$ on the NW slopes) and the seismic activity weakened, the number of large rockfalls decreased, and we observed a gap of activity in this area during the following years (Fig. 4c). Finally, after a few years
<figure><img src="image_4.png" /></figure>

**Fig. 3** Temporal evolution of rockfall (RF) activity compared to seismicity and rainfall. The vertical gray lines and rectangles show the start times and durations of the eruptions. **a** Cumulative number (red curve) and volumes (orange curve) of rockfalls (V > 100 m³). The thick dashed red line is the Omori fit to the whole catalog over the first 2.5 years. The black arrow highlights the end of the relaxation period. **b** Temporal evolution of shallow volcanotectonic (VT) events plotted as cumulative number (black) and rainfall recorded at the summit (light blue). The gray arrow highlights an important increase in seismicity that is not followed by an increase in rockfall activity. On the contrary, the black arrow shows another increase in seismicity clearly related to an increase in rockfall activity. **c** Deformation of the summit of the terminal cone is shown as distance changes between two pairs of GNSS stations. The locations of GNSS stations are indicated in Fig. 2.

of low seismic activity, we observe the re-occurrence of large rockfalls, that may be linked to the renewal of seismic activity (Fig. 4d). From 2014, the largest rockfalls cluster again on the west and north-west steepest slopes as already observed using a structure-from-motion method²⁷. This is the place where the N-S fractures were mainly opening E-W and SE-NW²⁷, and so were strained in the long-term by the dominant E-W deformation of the edifice (DERG-BORG baseline in Fig. 3c), leading to a faster opening during this period.

This qualitative spatio-temporal analysis of the rockfalls suggests the existence of three different behaviors of the slopes that are linked to their stability state: (1) after the crater floor collapse in 2007 (Period 1) the slopes are very unstable and quasi-spontaneous large rockfalls occur; (2) during the following period (Period 2), the slopes are more stable and need a triggering mechanism like seismicity or rainfall to collapse; (3) in the absence of seismic activity (Period 3), the slopes continue to evolve under the long-term competitive effects of gravity, erosion, and chemical healing, generating sparse large spontaneous or rain-triggered rockfalls; (4) finally, when the seismic activity increases (Period 4), it contributes to the collapse of

numerous rockfalls, by moving forward the time of the slope ruptures.

**Relative impacts of seismicity and rain.** To better identify and quantify the link between the 2007-2017 rockfall activity and seismicity and rainfall forcing, we use statistical methods that allow us to check the significance of our observations (see Methods for details). The combination of various methods, as much as possible, frees the results from the pitfalls of each individual method. By calculating the change of seismicity and rainfall rates preceding and following each rockfall (Fig. 5a, b), we investigate whether the forcing rates are higher before the occurrence of the rockfalls, and if this observation is statistically significant.

In the tests shown in Fig. 5, we consider all rockfall events with V ≥ 100 m³, independently of their size. To confirm the role of seismicity and rainfall in slope destabilization, we performed these tests considering all detected rockfalls (see Supplementary Fig. S2 and S3), and considering only the background rockfall activity, therefore discarding the correlated interactions between
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Spatio-temporal evolution of the largest rockfalls (V > 5000 m³). The four panels show the spatio-temporal distribution of the largest rockfalls during (a) Period 1, (b) Period 2, (c) Period 3 and (d) Period 4. The size of the circles scales with the volumes of the rockfalls and the colors represent the time from the beginning of each period.</figcaption></figure>

rockfalls (Fig. 5). For this, we applied an earthquake declustering method²⁸ to our rockfall catalog (Supplementary Fig. S4). However, it should be noted that not all the clustered rockfalls might be removed, due to potential errors on rockfall locations. In the tests, for each time period (Periods 1–4 indicated in Fig. 5a), we sum up, for the whole rockfall set, the cumulative number of volcano-tectonic events (VT) and rainfall height before and after each event (black curves in Fig. 5c–j). An asymmetry with a higher slope before the rockfall occurrence (left side of the vertical red line) than after means that, on average, the rockfalls are preceded by higher and thus probably triggering seismic or rainfall rates.

In order to see the degree of significance of the observed asymmetry, we compare these observations to simple statistical simulations performed by generating random rockfall catalogs based on the statistical characteristics of the observed catalog (gray dashed lines in Fig. 5c–j and see Supplementary Fig. S5 for the parameters used to generate the random catalogs). The gray dashed lines before and after the rockfalls represent simulated rates estimated without consideration of when rockfalls occur, thus representing rates when no forcing occurs. As expected, these lines are almost symmetric for all periods. A higher rate (i.e. steeper slope) of averaged VT or rainfall is considered significant if the observed black curve is above the interval between the two simulated gray lines. For this assumption to be valid, black curves after rockfalls should be roughly within the simulated dashed lines, as observed approximately in all cases.

We show here the results for the declustered catalog because most of the potential correlated interactions between rockfalls are

removed (Fig. 5). However, the general conclusions are the same for both catalogs (see Supplementary Fig. S2 and S3 for the results with the whole catalog). The main difference is observed during the first months after the collapse where the background rockfall activity does not exhibit any dependency on external forcings (Fig. 5). This suggests that the dependency found for the whole catalog during these first months (Supplementary Fig. S2) is artificial and simply results from interactions between rockfalls, as aftershocks of the crater collapse. Another way of interpreting the difference between the background and the whole catalogs is to consider that during Period 1, rockfalls and volcano-tectonic events are both related to the edifice damage following the crater collapse. As to whether a rockfall may trigger other rockfalls in the surrounding area, also submitted to seismicity forcing, is beyond the scope of this paper.

In addition to this analysis, we also compare the forcings during the preceding days in the cases of small and large rockfalls (Supplementary Fig. S7 and S8). We reach the same conclusions for all the rockfall subsets. However, we note a higher impact of rain on the destabilization of larger volumes (greater than 1000 m³) during Period 2 (Supplementary Fig. S7). This higher sensitivity during the period following the relaxation of the crater might be due to the fact that the slopes are still unstable, following the crater collapse. Let us discuss below the rate asymmetry observed in Figure 5c–j in the four different periods.

The apparent correlation between volcano-tectonic and rockfall activities observed during the early stage of Period 1 (until the beginning of 2008, Fig. 3a) can be explained either by the fact that the slopes are very unstable and therefore produce a high number
<figure><img src="image_6.png" /></figure>

**Fig. 5** Comparison of background rockfall activity to volcano seismicity and rainfall at the summit. **a** Number of volcano-tectonic events (VT, black) vs. the number of rockfalls (RF, red), both summed over 15 days. The rockfall catalog used is the declustered catalog. The vertical gray lines and rectangles show the start times and durations of eruptions. The horizontal black lines delimited by dots mark the different time spans used in (c) to (j). **b** Rainfall amount (light blue) vs. the number of rockfalls (red), both summed over 15 days. The light blue rectangle indicates the period over which no rain data is available. **c–f** Statistical significance of the effect of seismicity on rockfall triggering. The plots are centered around the time of the rockfall occurrence (red vertical line). The dashed gray lines show the 95% confidence intervals for the simulations (based on a randomly generated catalog), i.e. the zones where the link between seismicity and rockfall occurrences is not significant. The black lines represent the seismicity data stacked for all rockfalls of the given period. The periods used for each plot are indicated in (a). For more details on how these plots are made, see Materials and Methods. **g–j** Same as (c–f), for rainfall.

of rockfalls in response to the strong seismic excitation, or by the fact that rockfalls and seismic events at this time both respond to a common forcing, i.e. the crater floor collapse. On the contrary, from early 2008 to the end of 2009, comparing the temporal evolution of the rockfalls (red curves in Figs. 3a, 5a, b) with that of the seismicity (black curves in Figs. 3b and 5a) and of the rain

(light blue curves in Figs. 3b and 5b), we note that neither the short term variation in rain rate nor that of seismicity produce any visible change in the rockfall rate. Figure 5c and g does not highlight any clear correlation with seismicity or rainfall. The rock mass is not particularly responsive to these stimulations: at the end of 2008, the highest seismic activity over the 10 years
(indicated by the gray arrow in Fig. 3) had no strong impact on the rockfall activity. As discussed in results (section The temporal evolution of the rockfalls shows the global relaxation of the crater), the rockfall activity during this period is dominated by “aftershocks” of the crater floor collapse.

During the following years, when the slopes were no longer under the influence of the collapse (Period 2), we observed a very strong relationship between the increases in seismic and rockfall activities (Fig. 5a). Until 2011, all the increases in seismicity rate are concomitant with an increase of the number of rockfalls (Fig. 5a). This is also visible in the statistical tests, where the curves of the seismicity observations are clearly asymmetric and the pre-rockfall slope is higher than those of the simulations (Fig. 5d), indicating a clear effect of seismicity on rockfall occurrences. The second period presenting high seismic activity, beginning in mid-2014 (Period 4), is characterized by impulsive seismic episodes, with numerous events (up to 1200 per day) over a very short time period (a few days, generally less than a week) that also trigger rockfalls (Fig. 5a, f). However, we note a change in the response of the rockfalls to the seismic excitation compared to the 2009–2011 period: the increases in rockfall activity are not concomitant with the increases of seismicity, but rather occur after the end of the seismic crisis, with a delay of 5 to 10 days after the seismicity peak (Period 4 in Fig. 5a). This is also shown in Fig. 5f by the change in the slope between 5 and 10 days before the rockfalls (black arrow). We performed the same analysis using the seismic moment of the volcano-tectonic events, instead of their number, and we observed a similar pattern, i.e. a more intense seismic activity, with a higher cumulated moment preceding the rockfalls during Periods 2 and 4 (Supplementary Figs. S9C and S9E). Interestingly, Figure Supplementary S9E does not show any delay in the response of rockfalls to seismicity contrary to Fig. 5f. Overall, the cumulative seismic moment underlines the role of large magnitude volcano-tectonic events (Supplementary Fig. S9E) contrary to the cumulative number of events, that is dominated by small events. Thus these two figures suggest that relatively large VT may trigger rockfalls while smaller ones will have to accumulate and will induce a delayed triggering of rockfalls. During Period 4, we also observe that some seismicity increases are not followed by activation of the rockfall activity (Fig. 5a). Some of these lack of triggering might be a detection artifact: during active periods in terms of eruptions, the seismic waveforms are very noisy, with rockfall and volcano-tectonic event signals as well as tremor signals generated by the magma intrusions. The rockfall signals being less impulsive than the seismic events, they are more likely to be drawn into the noise generated by the tremors, and therefore to be undetected.

During periods of high seismic activity, rainfall does not appear to have a clear impact on rockfall triggering (Fig. 5g, h, and j), suggesting that the main triggering factor during these periods is seismicity. During the 4 years of non-eruptive activity, from 2011 to 2014 (Period 3), seismic activity was low, with more than 10 times fewer events than during the other periods (1.2 seismic events per day on the average, compared to 15.7 per day for the other periods) and lower magnitudes (maximum magnitude below 2). This is when the impact of rain on the triggering of rockfalls is more visible (Fig. 5b), as we note in particular a large increase of rockfall activity in 2013 (black arrow in Fig. 5b) that correlates with a rainfall episode. The higher rainfall activity preceding the rockfalls is slightly visible in Fig. 5i, showing an asymmetry with a higher rain rate before the rockfall time than after. The effect of rainfall is more visible when using the whole rockfall catalog. If the sensitivity of background rockfalls to rainfall is considered to be significant (Fig. 5i), Supplementary Figure S2i may either reflect that rainfall triggers clusters of rockfalls or that a large rockfall may trigger other rockfalls around

in areas also submitted to rainfall-induced increase of pore fluid pressure.

The slopes of the crater are also exposed to the global deformation of the summit measured by GNSS stations (Fig. 3c). We checked the impact of this deformation on the slope destabilization using statistical tools (see Supplementary Figs. S10 and S11 and Supplementary Methods for details), and showed that the changes in the inflation rates are not a direct triggering factor for rockfalls. However, they contribute, with gravity and erosion, to the long term destabilization of the slopes. Since in this study, we focus on the short term triggering of rockfalls, in the following discussion we will consider only the impact of seismicity and rainfall.

**Statistical significance of the results.** To further confirm the significance of our results, we first summarize the results on the seismicity and rain rates for the 15 days long pre- and post-rockfall periods, as presented in Fig. 5, by averaging the values of these rates (i.e., the slopes in Fig. 5) over sliding periods of 4 months. The results are presented in Fig. 6a and b, which show a non-symmetric distribution, with more volcano-tectonic activity before rockfalls than after them, and, similarly but less clearly, with higher rainfall levels before rockfalls than after them.

For quantifying the degree of significance of these asymmetries, we consider another statistics: the difference between the rates during the 15 days preceding and the 15 days following each rockfall, normalized by the square root of the number of events (volcano-tectonic events, or rainfall data points) in the one month window centered on the rockfall time. This normalization allows the mixing of periods of very different activity levels. It provides a bell-shaped, Gaussian-type distribution considering all the rockfalls, as presented in Fig. 6c. For the seismicity rates, more positive differences are counted than negative ones, i.e., the pre-rockfall rate of seismicity is larger than the post-rockfall rate, suggesting some triggering capability of rockfalls by the volcano-tectonic events. This is also the case for the rainfall rate (larger before the rockfalls), but it is less prominent (Fig. 6d). In the Methods section, we show that the larger number of positive differences in Fig. 6c, i.e., of larger seismicity rates before rockfalls than after them, is statistically significant: seismicity does trigger rockfalls. The same test for rainfall does not allow us to reject the hypothesis of no correlation between rainfall and subsequent rockfalls, but it may not be appropriate, owing to the noisier structure of the corresponding distribution. Finally, one finds an excess of 4.8% of rockfalls having a larger rate of volcano-tectonic events in the two weeks before them than in the two weeks after. In short, about 5% of the rockfalls are triggered by seismic events.

The stars in Fig. 7a, shown as significantly abnormal in the distributions of Supplementary Figure S6, indicate the time periods when the impacts of seismicity (black) and rainfall (blue) on rockfall production are the strongest. This supports and confirms the previously presented observations. Thus, the seismicity appears to have no impact on rockfall triggering during the relaxation phase (up to the second half of 2009). On another hand, rain seems to trigger some rockfalls during the most active period after the collapse, but has no impact in the following years until 2011. When the slopes are no longer under the influence of the collapse, but the edifice is seismically active (Period 2), Fig. 7 confirms the significance of the impact of the volcano-tectonic events on rockfall production. Finally, as stated before, the effect of rainfall is clearly visible only when the seismic activity is weak (Period 3). Our observations and the statistical tests suggest an influence of rainfall on rockfall production, but this influence is not systematic and is independent of the amount of rain. For example, in May 2013, there was a major increase in
<figure><img src="image_7.png" /></figure>

**Fig. 6** Significance of the observations of triggering by seismicity and rainfall. **a** Slope of the stacked seismicity before (R<sub>b</sub>) vs. after (R<sub>a</sub>) the background rockfalls. The black line is where the two slopes are equal (i.e. same cumulative number of volcano-tectonic events before and after the rockfalls). The slopes are computed on stacked data in 4-month long moving windows, each shifted by 15 days from the previous one (i.e. one point per 4-month period). The stars show when the slope before the rockfalls is significantly larger than the slope after (determined from Supplementary Fig. S6). The dots and stars are color-coded to show time. **b** Same as (a) for rainfall. See Materials and Methods for more details on the computation of the slopes. **c** Distribution of the normalized difference between the number of volcano-tectonic events preceding and following the rockfalls. The vertical blue and red arrows indicate the mean values of the symmetric distribution (i.e. if no influence of the seismicity on the rockfall occurrences) and the observed distribution, respectively and the horizontal blue and red arrows indicate the +/−σ range (see the Methods section for the values). **d** Same as (c) for the rain. The rockfall catalog used for this figure is the declustered catalog.

rockfall activity not directly explainable by higher rainfall activity (thickest black arrow in Fig. 5b), in the sense that the same amount of rain in 2012 was not followed by a significant increase of the rockfall activity (Fig. 5b). This is one of the highest rockfall activities during the 10 years (Figs. 3 and 5), involving the whole crater (Supplementary Fig. S12). This episode is not preceded by any noteworthy seismic activity (it occurs during the seismic quiescence period) and the rainfall occurring just before is not the highest of the study period: as stated earlier, six months before, an even higher rain episode occurred without triggering rockfall activity. A potential explanation for this destabilization in 2013 might be the combination of rainfall with the occurrence of strong winds generated by two tropical cyclones at the beginning of 2013. These cyclones exhibited winds of 180 km/h and 133 km/h at the volcano (January and February 2013, https://reunion-

extreme.re/cyclones-a-la-reunion.html, Météo France). A recent study<sup>29</sup> suggests that strong winds might play a role in the triggering of rockfalls, particularly in a seismically active region where seismicity can generate dynamic cracks in the slopes, thus weakening them.

**Lab-scale experiments on the acoustic triggering of granular flows.** To provide insight into the physical processes behind our field-scale observations we show some similarities with lab-scale experiments on cohesionless (dry) granular layers where flows are triggered by small-amplitude ultrasound (continuous wave) (Fig. 8a). In these experiments, the applied ultrasound vibrations mimic the small-magnitude seismic forcing and the triggered granular flows may be compared to triggered rockfalls. In the case
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Temporal evolution of the triggering effect of seismicity and rainfall on background rockfalls. **a** Temporal evolution of the difference between the slopes before and the slopes after the background rockfalls for seismicity (black circles) and rainfall (light blue diamonds). The slopes are computed on stacked data in 4-month long moving windows, each shifted by 15 days from the previous one. As in Fig. 6, the stars denote the times when the slope before the rockfalls is significantly larger than the slope after, i.e. when rockfalls are likely to have been triggered by seismic activity (black) or rainfall (blue). **b** Number of volcano-tectonic events (VT, black) and rainfall amount (light blue) vs. the number of rockfalls (RF, red) summed over the same 4-month long moving windows as in panel **a**. This plot shows the number of seismic events, the amount of rainfall, and the number of rockfalls (using the declustered catalog) in the time-windows used to compute the slope differences.</figcaption></figure>

of dry granular packings, these vibrations have been shown to reduce the maximum angle of slope stability $\theta_m$ (avalanche angle) due to the decrease of the static coefficient of interparticle friction $\mu_p^{30}$. Such an ultrasound-induced triggering is also observed in granular packings saturated by water, confirming the principal effect of ultrasound on the acoustic lubrication of grain contacts (preliminary results not yet published). When the slope is close to the avalanche angle $\theta_m$, triggered grain motions quickly reach fast inertial flows (Fig. 8a), whose dynamics are mostly independent of vibration amplitude. In contrast, when the slope is below the angle of repose, 2°–3° smaller than $\theta_m^{31}$, the ultrasonic vibration triggers slow creep-like motions of grains where the flow stops upon the removal of ultrasound. Our observation of inertial granular flow triggered by ultrasound is consistent with the observed high sensitivity of rockfall activity to seismicity both in Period 2 (2009–2011) and Period 4 (2014–2017) for stable slopes. The underlying physics of frictional granular flows triggered by small-amplitude acoustic vibrations can be explained by the competition between rejuvenation (damage due to slipping) and healing (growth or strengthening) of the contacts³². As healing dominates for stable slopes another rejuvenation mechanism such as acoustic fluidization of the medium via the lubrication of the contacts³³ is necessary to drive the frictional system to failure. In the case of a low-amplitude acoustic input (continuous) or burst-like excitation (low energy), delayed triggering of a large avalanche is expected, which may explain partly those observed in Period 4 as discussed below.

For real cohesive granular materials such as the volcanic slopes, the failure may not occur spontaneously by vibration-induced material weakening or slipping but rather through progressive damage nucleation (or fatigue) to failure when the damaged zone is large enough or/and the medium is sufficiently weakened to

rupture (with a longer delay) as shown in Fig. 8c³⁴,³⁵. The time scale of the failure mechanisms driving the rockfall activity has thus to be compared to the duration of the seismic forcing. In the case of the very unstable slopes after the crater collapse, the growth of the fractures leading to the rockfalls is expected to be much faster than later, when the slopes were more stable (Supplementary Fig. S13). In other words, the duration of damage propagation, from initiation to the final rupture, is expected to be longer for more stable slopes, leading to a longer effective time window during which the stress perturbations generated by the external forcings such as seismicity could interact with the damage growth and advance the final rupture time. This hypothesis leads to a higher probability that, on stable slopes, low level perturbations from seismicity and rain may find many mature rock masses with an initial strength close to the threshold, thus triggering their advanced failure, and statistically out-numbering the slowly self-rising ruptures. On the contrary, for unstable slopes, the window to failure is too narrow to be effectively influenced (or triggered) by additional, sparse stresses from forcing processes (see below).

## Discussion

Our results provide insight into the response of a metastable volcanic edifice to external forcings. A main finding is that even if volcano-tectonic seismicity and rainfall both trigger rockfalls, seismicity is the main triggering factor in this volcanic environment. Rainfall seems to have a limited impact, detectable only when the seismic activity is negligible, and it becomes the principal triggering mechanism. This could be partly due to the highly fractured slopes of the Piton de la Fournaise volcano that may easily drain water, leading to a fast relaxation of the associated
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Schematic illustration of granular flows (rockfalls-like) triggered by small-amplitude ultrasonic (seismic) vibration. **a** Granular layers (rockfalls) on an inclined plane at angle $\theta$ below the maximum angle of stability $\theta_m$ ($= \tan^{-1}\mu_s$)(see text). The inertial flow triggered by ultrasonic vibration is mostly uniform and continuous (shown by arrows). **b** Typical plot of shear stress $\tau$ vs strain $\gamma$ in a sheared granular gouge (inset) with $\tau_{s,d} \sim \mu_s$, $d\sigma_n (\mu_{s,d} \sim \tan\theta_{m,r})$ being the static and dynamic coefficients of friction. Small-amplitude vibrations can progressively weaken the shear modulus (rigidity) from $G_0$ to $G_f$ and lower the shear strength (yield stress) from $\tau_s$ to $\tau_s^*$, triggering finally advanced failure for an imposed shear (slope) $\tau < \tau_s$. **c** Schematic illustration of abrupt drops (and then slower decreases) in the shear modulus and strength of granular gouges or rocks due to strong pulsed vibrations (ultrasonic or seismic), followed by thermally-activated healing of frictional contacts or fractured bonds. The Slow Dynamic (SD) recovery may have a function in delayed dynamic triggering where repeated or burst-like ultrasound (seismicity) (i-iv) lead progressively the systems to rupture (failure).</figcaption></figure>

increase in pore fluid pressure. The global deformation of the volcano also might play a role in the long term destabilization, making the slopes more prone to rupture, but will not act as the final trigger at the time scale of one or two weeks. At an even longer time scale (of several months), healing related to chemical processes also affects the slope by increasing its stability. For example, in Period 3 with almost no seismicity, one could expect that healing may be more efficient, decreasing porosity and draining and thus increasing the triggering effect of rain. These combined effects may explain why we cannot identify a clear triggering threshold for rainfall or seismicity and why some relatively small seismicity episodes are followed by an increase in rockfall activity.

Another key point is that even low magnitude seismicity (average magnitude of 0.6, maximum magnitude 3.8) can destabilize the slopes. This observation is all the more important as only a few studies have looked at the influence of low magnitude seismicity ($M<4$) on the triggering of slope instabilities. Furthermore, the results of these previous studies are often not statistically significant given the small amount of data, and so do not make it possible to establish a clear link between low magnitude seismicity and the triggering of rockslides³,¹¹,¹⁵. Our exceptionally long period of observations (ten years of continuous data) has made it possible for us to infer that in the case of low magnitude seismicity, one earthquake may not be enough to destabilize the slopes, even if it is located very close (less than 2 km below), while the accumulation of vibrations produced by repeated seismic activity will trigger instabilities. Since it is the shaking generated by the seismic events that trigger the unstabilities, the distance-

magnitude relation plays an important role in the triggering potential of the seismicity. However, in the particular case of the Piton de la Fournaise, the distance between volcano-tectonic events and the destabilized slopes is roughly constant, since all the seismic events occur approximately at the same depth. Consequently, the variation of the shaking amplitude is mainly controlled by the magnitude of the volcano-tectonic events. Our analysis of the impact of the seismic moment does not highlight a clear magnitude threshold for the triggering of the rockfalls. Nevertheless, the fact that the test on the cumulative seismic moment during Period 4 (Supplementary Fig. S9E) does not show any delay in the triggering of rockfalls suggests that relatively large VT events may trigger rockfalls, while smaller ones will need to accumulate.

Interestingly, the triggering potential of seismicity and rain is shown to depend on the stability state of the slopes, but not in the most intuitive manner. During the relaxation stage following the crater floor collapse, when the crater walls were very unstable, the rockfalls were found not to be particularly responsive to seismicity and rainfall. Following the stabilization of the crater, the rockfalls become more sensitive to seismicity due to the longer duration of damage propagation up to failure on more stable slopes, providing a longer period of potential sensitivity to seismicity (Supplementary Fig. S13). Subsequently, during the second seismically active period (Period 4 from mid-2014), not all the seismicity increases were followed by intensification of rockfall activity. This can be explained by the fact that the very short seismic crises (burst-like), despite the high number of events, were insufficient in energy and less efficient in triggering
rockfalls on these stable slopes than the longer seismic crisis occurring before 2011. A complementary explanation could be that the slopes were more stable and purged, and so less prone to generate numerous rockfalls. Comparing the two seismically active periods 2 and 4, the response delay of the slopes (i.e. destabilization) increased by several days as the crater stabilized. This observation of delayed triggering can be compared with observations reported for earthquakes in tectonic areas$^{4,36}$.

The complex behavior of the lab- and field-scale measurements described in this paper, also observed for earthquakes, has been shown to be well described by the rate and state variable friction (RSF) laws$^{37}$. In the framework of such laws, we infer that the repetitive vibrations generated by the volcano-tectonic earthquakes weaken the shear strength $\tau_s$ (which corresponds to $\tan^{-1}\theta_m$ in the lab-experiments of granular flows described above) and/or the cohesion of fractured rock masses by cumulative breaks or losses of contact (Fig. 8c). The age of contact (healing) or the history of cohesion on the fractured rock masses play the role of the state variable$^{4,32,37}$ and determine the delay of dynamic triggering up to failure, for given initial conditions$^{38-40}$. This delay may thus represent the time needed for the fractures to be developed to their full length across the system. Consequently, on more stable slopes (Period 4, 2014–2017) with smaller energy of seismicity (burst-like), this growth takes more time than on less stable ones (Period 2, 2009–2011). Furthermore, the potentially higher initial friction and strengthening when the slopes were more stable, i.e. from 2014 to 2017, may also contribute to explaining the longer time delays to rupture observed during this period. Therefore, we propose a scenario for the global evolution of the system stability schematically illustrated in Fig. 8c where the abrupt decrease of the shear strength and elastic modulus (rigidity) is caused by breakage of contacts due to the vibrations, followed by logarithmic-like, thermally-activated healing of damaged frictional bonds$^{34,35}$ or chemical activation effects$^{4}$, also called Slow Dynamic recovery (SD). Such a fracture process is somehow akin to the fatigue process of a material under cyclic loading up to failure$^{4}$. More specifically, these seismicity-induced rigidity reduction and micro-damage recovery have been monitored by elastic wave velocity measurements in lab experiments$^{33,35}$ and field observations by ambient noise detections$^{5,41}$. Nevertheless, the driving processes are complex and will need to be addressed in further studies, including the long-term impact of mechanisms such as static deformation and interactions with fluids on the destabilization potential of the slopes.

The observation of the destabilization of the slopes by low but repetitive burst-like seismicity has implications for risk management linked to slope instabilities in inhabited areas, suggesting that the existence of seismicity swarms might have to be taken into account in mountainous and especially volcanic regions at distances of a few km from the low magnitude seismic sources. It also has a strong implication for monitoring and planning of human activities generating seismicity swarms such as repetitive explosions in mining and underground exploitation for extraction or storage of CO2, nuclear, or natural resources.

We want to emphasize that the study we report here is focused on observations, and prepares the ground for more systematic studies of slope instabilities by showing the relevance of applying statistical tools originally developed for earthquake catalog analysis, thus suggesting similarities in the failure processes. Our observations are also of interest to the seismology and civil engineering communities which are concerned by the action of cumulative seismic activity on the destabilization of faults or on the dynamic response of buildings, in particular during periods of aftershocks.

## Materials and methods

**Data.** The Piton de la Fournaise volcano is densely instrumented with, among other instruments, 39 seismic and 24 Global Navigation Satellite System (GNSS) stations. For the 2007–2017 study period, we consider only the ten seismic stations and the four GNSS stations indicated in Fig. 2 (white and black dots, respectively). We use the seismic signals generated by rockfalls recorded at the seismic stations to locate the rockfalls. To compute their volumes (for details on the method see Hibert et al. 2011$^{18}$), we use only the stations closest to the crater (BON, BOR, and SNE). We discard DSO station because of its low saturation level. BOR station has been recorded continuously since 2007. BON has been active only since 2011 and before 2011, SNE station was a short-period seismometer. Consequently, the volumes from 2007 to 2011 are the volumes computed at BOR station, whereas the volumes from 2011 to 2017 are the average volumes at the three stations.

GNSS stations are used to quantify the deformation of the summit. GNSS data were post-processed with the GAMIT/GLOBK software package$^{42}$. GAMIT uses (i) the precise ephemerides of the International GNSS Service (IGS); (ii) a stable support network of 20 IGS stations not located on La Reunion Island but scattered elsewhere in the Indian Ocean; (iii) a tested parameterization of the troposphere; and (iv) models of ocean loading and Earth and lunar tides. Plate motion used to correct data is deduced from the REUN IGS station located 15 km to the west of the summit and assumed not to be affected by any volcanic deformation. The final daily accuracy of the measurements is of the order of millimeters.

Daily rain data are obtained by averaging the measurements at three rain gauges located close to the summit of the cone (gray dots in Fig. 2). Furthermore, three cameras continuously record the activity in the crater. They are distributed all around it.

**Statistical analysis.** To highlight potential changes in volcanic tectonic activity, rain, and deformation, we stack the corresponding time-series preceding and following each observed rockfall ($V \ge 100 \text{ m}^3$) as follows:

$$ \sum_{i=1}^{N} f_{\text{obs}}(t_i - \Delta t : t_i + \Delta t) \qquad (2) $$

with $f_{\text{obs}}$ the symmetrical time series defined by

$$ f_{\text{obs}} = \sum_{j=1}^{k} g_{\text{obs}}(t_j), \text{ for } t > t_i, \text{ and} $$

$$ f_{\text{obs}} = \sum_{j=1}^{-k} g_{\text{obs}}(t_j), \text{ for } t < t_i, \text{ with } k = 1, \dots, m, \text{ and } t_m = \Delta t $$

with $g_{\text{obs}}$ being the volcano-tectonic, rain, or deformation time series, $t_i$ the time of the $i^{th}$ rockfall, and $\Delta t$ the considered time period before and after the rockfall occurrences (15 days), and N the total number of rockfalls.

For Fig. 5, we consider the time-series (seismicity and rainfall) preceding and following the rockfalls contained in the periods indicated in Fig. 5a.

To statistically test the significance of our observations, we compare them to the results obtained for randomly generated rockfall catalogs. The inter-event times of rockfalls follow a gamma distribution (see Supplementary Fig. S5):

$$ p(\Delta t) = \Delta t^{(\gamma-1)} \frac{e^{-\Delta t/\beta}}{\beta^{\gamma}\Gamma(\gamma)} \qquad (3) $$

where $\Delta t$ is the inter-event time and $\Gamma(x)$ is the gamma function. $\beta$ and $\gamma$ are computed using the mean $\Delta t$ and the variance $\sigma_{\Delta t}^2$ of
the inter-event times: $\beta = \sigma_{\Delta t}^2/\bar{\Delta}t/\beta$ and $\gamma = \bar{\Delta}t/\beta$ (from Aki and Richards 1980$^{43}$). The values for $\beta$ and $\gamma$ are indicated for each period in Supplementary Fig. S5.

We fit the inter-event time distributions of rockfalls using a different gamma distribution for each period (see Supplementary Fig. S5). The power ($\gamma - 1$) provides the constant slope in the log-log plot of Supplementary Fig. S5. $\gamma$ getting closer to 1 indicates getting closer to a Poisson distribution, characteristic of the absence of clustering (i.e. independent events). Supplementary Fig. S5 thus suggests that the rockfalls are more clustered in Period 1, than in Periods 2 and 4, while rockfalls in Period 3 (with almost no seismicity) are more independent from each other. Then we use the $\beta$ and $\gamma$ parameters determined from the observed catalog to generate random catalogs following the same distribution for each period. We generate 500 random catalogs for each period. For each random catalog, we stack the external factor time-series (volcano-tectonic events, rain, and deformation) preceding and following each randomly generated rockfall. For clarity, we represent in Fig. 5 the 95% confidence interval for the 500 simulations. If rockfall occurrences are independent of the seismic activity, rainfall, and/or deformation, the observed stacked temporal evolution plots should be symmetric (same activity before and after the rockfall occurrence) and fall inside the 95% confidence interval of the simulations. An asymmetry with a higher rate before the rockfall occurrence (left side) means that most of the rockfalls are preceded by higher seismic activity, rainfall, and/or deformation. However, these simulations present some limits related to the gamma distribution used to generate random rockfall catalogs (Supplementary Fig. S5). This distribution does not take into account the possible non-stationarity of the gamma distribution, i.e. the variation within a given period of the degree of clustering. Hence the synthetic catalogs are not an exact replica of the observed catalog.

To avoid problems related to the possible over-simplification of the gamma law statistics of the synthetic test above, we directly analyze the distribution of the observed rates, as plotted in Fig. 6. For this, we consider the time-series preceding and following the rockfalls contained in 4-month long moving windows, each shifted by 15 days from the previous one. These 4-month long windows contain enough rockfalls to allow a statistical representation. In Fig. 6, we plot the rates before ($R_b$) as a function of the rates after ($R_a$) for each 4-month period: each point corresponds to a pair ($R_a, R_b$) for a given 4-month period. Under the null hypothesis of independence of rockfalls with respect to seismic shaking and rain, a symmetrical distribution of the points is expected with respect to the diagonal, with a rate before ($R_b$) nearly equal to the rate after($R_a$). For a given forcing, i.e. rain or seismicity, with N pairs ($R_a, R_b$) randomly distributed on each side of the diagonal $R_b = R_a$, we have $N_b$ (respectively $N_a$) points with the difference $R_b - R_a > 2500$ (respectively $< 2500$) and $N = N_a + N_b$. We choose 2500 as a threshold based on the distributions of Supplementary Figure S6 (equivalent to one standard deviation). The pairs with the difference $R_b - R_a > 2500$ potentially highlight periods with a significant impact of the given forcing on the rockfall activity. To check for this significance, we test the probability $P_r(N_b)$ to have $N_b$ pairs with $R_b - R_a > 2500$ under the null hypothesis that $R_b - R_a > 2500$ and $R_a - R_b > 2500$ are equally likely to occur. For this, we perform a binomial test. If $P_r(N_b) < 5\%$, the null hypothesis can be significantly rejected, meaning that volcano-tectonic activity or rainfall has a significant impact on rockfall triggering. The results of the test are described in Supplementary Figure S3.

Let us evaluate the significance of the asymmetry observed in Fig. 6c for the seismicity rates. Taking $N_{b,i}$ (respectively $N_{a,i}$) as the number of volcano-tectonic events in the 15 days before

(respectively after) the rockfall $RF_i$, we consider the normalized statistics of the difference: $D = (N_{b,i} - N_{a,i})/\sqrt{(N_i)}$, with $N_i = N_{a,i} + N_{b,i}$. If there is no correlation between the rockfall occurrences and the seismicity, the sampled distribution of $D$ should be close to symmetric, with a variance of the order of unity, and a mean $D^*$ close to zero. To the first order, under the hypothesis of no correlation, the bell-shaped distribution of $D$ follows a Gaussian-type distribution. Thus, defining $S$ as the standard deviation of $D$, and $S^*$ as the standard deviation of $D^*$, one has: $S^* = S/\sqrt{(N^*)}$, where $N^*$ is the number of independent rockfalls. The comparison of $D^*$ and $S^*$ is our statistical test. The last step is to evaluate $N^*$, which differs from the total number of rockfalls $N$ considered in the distribution. Indeed, the rockfalls are not independent of each other, as clearly evidenced by the clustering deduced from the gamma law of their distribution (See Supplementary Fig. S5). The quantification of the latter shows that $x\%$ of the rockfalls are independent, so that one can approximate $N^*$ by $x/100*N$.

For the seismicity, one has: $N = 4272$, $x = 47$, $N^* = 2007$, $S = 12$, hence $S^* = 0.27$. The mean of $D$, $D^* = 1.19$, is four times larger than its expected standard deviation, leading us to confidently reject the null hypothesis of a symmetrical distribution, i.e., the hypothesis of no correlations between rockfall occurrences and seismicity. The sign of the asymmetry corresponds to a clear - although modest - triggering of rockfalls by volcano-tectonic events.

A similar analysis has been performed for the rainfall, giving $D^* = 0.12$ and $S^* = 0.63$. In that case, the mean is within one standard deviation, so this test does not allow to reject the null hypothesis that the rainfall rate has an effect on the rockfall rate. This test is however less robust, as the distribution of $D$ appears much more irregular than it is for the seismicity, so the Gaussian approximation may not be appropriate.

## Data availability

The maps presented in this publication have been generated using the GMT tool (https://wwwwatson mapping tools org/). The rockfall catalog used in this manuscript is available through IPGP Research Collection (research-collection.ipgg.fr) at the following link https://doi.org/10.18715/IPGP.2023.ljofdfqx. The volcano-tectonic and rain data are available at the following: https://doi.org/10.18715/REUNION.OVPF, and on request on the volobsis.ipgg.fr portal.

Received: 13 May 2022; Accepted: 5 September 2023;

Published online: 18 October 2023

## References

1. Bontemps, N., Lacroix, P., Larose, E., Jara, J. & Taine, E. Rain and small earthquakes maintain a slow-moving landslide in a persistent critical state. *Nat. Commun.* **11**, 780 (2020).

2. D'Amato, J. et al. Influence of meteorological factors on rockfall occurrence in a middle mountain limestone cliff. *Nat. Hazards Earth Sys. Sci.* **16**, 719–735 (2016).

3. Tatard, L., Grasso, J.-R., Helmstetter, A. & Garambois, S. Characterization and comparison of landslide triggering in different tectonic and climatic settings. *J. Geophys. Res.115* (F04040) https://doi.org/10.1029/2009JF001624 (2010).

4. Scholz, C. *The Mechanics of Earthquakes and Faulting* (3rd Edition, Cambridge, 2019).

5. Brenguier, F. et al. Towards forecasting volcanic eruptions using seismic noise. *Nat. Geosci.* **1**, 126–130 (2008).

6. Dietze, M., Turowski, J. M., Cook, K. L. & Hovius, N. Spatiotemporal patterns, triggers and anatomies of seismically detected rockfalls. *Earth Surface Dynamics* **5**, 757–779 (2017).

7. Helmstetter, A. & Garambois, S. Seismic monitoring of sèchilienne rockslide (french alps): analysis of seismic signals and their correlation with rainfalls. *J. Geophys. Res.115* (F03016) https://doi.org/10.1029/2009JF001532 (2010).
8. Krautblatter, M., Moser, M., Schrott, L., Wolf, J. & Morche, D. Significance of rockfall magnitude and carbonate dissolution for rock slope erosion and geomorphic work on alpine limestone cliffs (reintal, german alps).
*Geomorphology* (2012).

9. Lin, G. W. et al. Effects of earthquake and cyclone sequencing on landsliding and fluvial sediment transfer in a mountain catchment. *Earth Surface Process. Landforms* **33**, 1354–1373 (2008).

10. Marc, O., Hovius, N., Neunier, P., Uchida, T. & Hayashi, S. Transient changes of landslide rates after earthquakes. *Geology* **43**, 883–886 (2015).

11. Keefer, D. K. Investigating landslides caused by earthquakes-a historical review. *Surveys Geophys.* **23**, 473–510 (2002).

12. Durand, V. et al. On the link between external forcings and slope instabilities in the piton de la fournaise summit crater, reunion island. *J. Geophys. Res.: Earth Surface* **123** https://doi.org/10.1029/2017JF004507 (2018).

13. Hibert, C. et al. Spatio-temporal evolution of rockfall activity from 2007 to 2011 at the piton de la fournaise volcano inferred from seismic data. *J. Volcanol. Geotherm. Res.* https://doi.org/10.1016/j.jvolgeores.2017.01.007 (2017).

14. Keefer, D. K. Landslides caused by earthquakes. *Geol. Soc. Am. Bull.* **95**, 406–421 (1984).

15. Del Gaudio, V., Trizzino, R., Calcagnile, G., Calvaruso, A. & Pierri, P. Landsliding in seismic areas: the case of the acquara-vadoncello landslide (southern italy). *Bulletin of Engineering Geology and the Environment* **59**, 23–37 (2000).

16. Lacroix, P. Landslides triggered by the gorkha earthquake in the langtang valley, volumes and initiation processes. *Earth, Planets Space* **68** (46) https://doi.org/10.1186/s40623-016-0423-3 (2016).

17. Dammeier, F., Moore, J. R., Haslinger, F. & Loew, S. Characterization of alpine rockslides using statistical analysis of seismic signals. *J. Geophys. Res.* **116** (F04024) https://doi.org/10.1029/2011JF002037 (2011).

18. Hibert, C., Mangeney, A., Grandjean, G. & Shapiro, N. M. Slope instabilities in dolomieu crater, réunion island: From seismic signals to rockfall characteristics. *J. Geophys. Res.* **116** (F04032) https://doi.org/10.1029/2011JF002038 (2011).

19. Hibert, C. et al. Automated identification, location, and volume estimation of rockfalls at piton de la fournaise volcano. *J. Geophys. Res.: Earth Surface* **119**, 1082–1105 (2014).

20. Observatoire Volcanologique Du Piton De La Fournaise (OVPF) and Institut De Physique Du Globe De Paris (IPGP). Seismic, tiltmeter, extensometer, magnetic and weather permanent networks on piton de la fournaise volcano and la réunion (2008).

21. Durand, V. et al. 2007-2017 Rockfall Catalog of Piton de la Fournaise https://doi.org/10.18715/IPGP.2023.ljofdfqx (2023).

22. Roul, G. et al. A new comprehensive classification of the piton de la fournaise activity spanning the 1985-2010 period. search and analysis of short-term precursors from a broad-band seismological station. *J. Volcanol. Geother. Res.* **241**, 78–104 (2012).

23. Staudacher, T. et al. The april 2007 eruption and the dolomieu crater collapse, two major events at piton de la fournaise. *J. Volcano. Geotherm. Res.* **184**, 126–137 (2009).

24. Institut De Physique Du Globe De Paris (IPGP). Data collection of the volcanological observatory of piton de la fournaise (2021).

25. Utsu, T., Ogata, Y. & Matsuura, R. The centenary of the omori formula for a decay law of aftershock activity. *J. Phys. Earth* **43**, 1–33 (1995).

26. Reasenberg, P. A. & Jones, L. M. Earthquake hazard after a mainshock in california. *Science* **243**, 1173–1176 (1989).

27. Derrien, A., Villeneuve, N., Peltier, A. & Michon, L. Multi-temporal airborne structure-from-motion on caldera rim: Hazard, visitor exposure and origins of instabilities at piton de la fournaise. *Prog. Phys. Geogr.* **43**, 193–214 (2019).

28. Zaliapin, I., Gabrielov, A., Keilis-Borok, V. & Wong, H. Clustering analysis of seismicity and aftershock identification. *PRL* **101**, 018501 (2008).

29. Pei, X., Luo, J. & Huang, R. Failure mechanisms of wind-induced post-seismic rockfall hazard. *Bull. Eng. Geol. Environ.* **78**, 5707–5725 (2019).

30. Léopoldès, J., Jia, X., Tourin, A. & Mangeney, A. Triggering granular avalanches with ultrasound. *Phys. Rev. E102* (042901) (2021).

31. Nagel, S. R. Instabilities in a sandpile. *Rev. Mod. Phy.* **64**, 321–325 (1992).

32. Baumberger, T. & Caroli, C. Solid fraction from stick-slip down to pinning and aging. *Adv. Phys.* **55**, 279–348 (2006).

33. Jia, X., Brunet, T. & Laurent, J. Elastic weakening of a dense granular pack by acoustic fluidization: slipping, compaction, and aging. *Phys. Rev. E* **84**, 020301 (2011).

34. Teencate, J., Smith, E. & Guyer, R. Universal slow dynamics in granular solids. *Phys. Rev. Lett.* **85**, 1020–1023 (2000).

35. Johnson, P. A. & Jia, X. Nonlinear dynamics, granular media and dynamic earthquake triggering. *Nature* **437**, 871–874 (2005).

36. Gomberg, J., Reasenberg, P. A., Bodin, P. & Harris, R. A. Earthquake triggering by seismic waves following the landers and hector mine earthquakes. *Nature* **411**, 462–466 (2001).

37. Marone, C. Laboratory-derived friction laws and their application to seismic faulting. *Annu. Rev. Earth Planet. Sci.* **26**, 643–696 (1998).

38. Belardinelli, M., Cocco, M., Coutant, O. & Cotton, F. Redistribution of dynamic stress during coseismic ruptures: evidence for fault interaction and earthquake triggering. *J. Geophys. Res.* **104**, 14925–14945 (1999).

39. Parsons, T. A hypothesis for delayed dynamic earthquake triggering. *Geophys. Res. Lett.* **32** (L04302) https://doi.org/10.1029/2004GL021811 (2005).

40. Freed, A. Earthquake triggering by static, dynamic and postseismic stress transfer. *Annu. Rev. Earth Planet. Sci.* **33**, 335–367 (2005).

41. Obermann, A., Planès, T., Larose, E. & Campillo, M. Imaging preeruptive and coeruptive structural and mechanical changes of a volcano with ambient seismic noise. *J. Geophys. Res.: Solid Earth* **118**, 6285–6294 (2013).

42. Herring, T. A., King, R. W. & McClusky, S. C. Gamit: Gps analysis at mit, 10.4". Massachusetts Institute of Technology, Cambridge (2010).

43. Aki, K. & Richards, P. G. *Quantitative seismology: Theory and methods* (CA: University Science Books, 2nd ed., 645 pp., Sausalito, 1980).

## Acknowledgements

This work was funded by the ERC project ERC-CG-2013-PE10-617472 Slidequakes. We thank the three anonymous reviewers and the editors for their constructive comments, that largely contributed to improve this paper.

## Author contributions

V.D. and A.M contributed to the Conceptualization, Formal Analysis, Investigation, Methodology, Writing of original draft. P.B and X.J contributed to the Methodology, Writing of original draft. F.B, C.S, J-M.S, A.EM, A.P, V.F, P.K, F.L, C.B and C.H contributed to the Discussion of the results and the Review of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-00996-y.

**Correspondence** and requests for materials should be addressed to Virginie Durand.

**Peer review information** *Communications Earth & Environment* thanks Joan Gomberg, Shengwen Qi and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Teng Wang and Joe Aslin. A peer review file is available.

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

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

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

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

© The Author(s) 2023