Earthquake statistics changed by typhoon-driven erosion

Tectonics and climate-driven surface processes govern the evolution of Earth’s surface topography. Topographic change in turn influences lithospheric deformation, but the elementary scale at which this feedback can be effective is unclear. Here we show that it operates in a single weather-driven erosion event. In 2009, typhoon Morakot delivered ~ 3 m of precipitation in southern Taiwan, causing exceptional landsliding and erosion. This event was followed by a step increase in the shallow (< 15 km depth) earthquake frequency lasting at least 2.5 years. Also, the scaling of earthquake magnitude and frequency underwent a sudden increase in the area where mass wasting was most intense. These observations suggest that the progressive removal of landslide debris by rivers from southern Taiwan has acted to increase the crustal stress rate to the extent that earthquake activity was demonstrably affected. Our study offers the first evidence of the impact of a single weather-driven erosion event on tectonics.

It is well established that tectonics and climate influence the evolution of the Earth's surface by modulating the rate and pattern of erosion and sedimentation 1,2 . In turn, theoretical predictions and numerical models suggest that changes of surface topography by surface processes can promote tectonic deformation over geological times [3][4][5] (1-10 Myr), enhance fault slip over intermediate time scales [6][7][8][9] (1 kyr-1 Myr), and induce sufficient static stress changes over a seismic cycle (1-1,000 year) to trigger earthquakes 10 . However, the influence of ongoing surface processes on tectonics has not been directly observed. Over seasonal timescales, surface (un)loading induced by rainfall or snow can modulate local to regional seismicity 11,12,13 . Here we ask if a single erosional event can have a discernable effect on seismogenic processes, which dominate deformation of the Earth's upper crust. Both large-magnitude earthquakes and extreme rainfall events can trigger instantaneous and widespread landsliding, driving the export of millions of tons of sediment from mountain areas over periods of months to years 2,14 . At these timescales, geophysical methods allow monitoring of changes in earthquake activity associated with erosional perturbations. Compared to large-magnitude earthquakes, rainfall events do not directly trigger aftershock sequences (even if surface (un)loading can modulate seismicity), which could preclude the observation of a seismicity change induced by erosion. We therefore focus in the following on detecting seismicity changes associated with an erosional perturbation induced by a rainfall event. temporal changes in seismicity after Morakot. To determine how the erosional unloading due to typhoon Morakot may have impacted fault dynamics, we analyze the evolution of shallow (< 15 km) seismicity in Taiwan after it made landfall. Because the expected stress change is small compared to background tectonic loading at seismogenic depth 10 , we focus on detecting changes in the statistics of recorded seismicity, such as the earthquake frequency, seismic moment rate and the b-value of the Gutenberg-Richter earthquake size distribution (see "Methods" section), rather than on individual events. We use the seismicity catalogue of the Central Weather Bureau of Taiwan, which includes > 340,000 earthquakes during the period 1995-2015 over Taiwan island.
First, we assess the time evolution of earthquake statistics using a temporal sliding window of 1,001 earthquakes (see "Methods" section). Results show a stepwise increase of the frequency of shallow (< 15 km) earthquakes in the landsliding zone after Morakot (Fig. 2, Supplementary Fig. S3). Although this frequency increase is orders of magnitude lower than after the Chi-Chi earthquake, it is observed both for earthquakes with magnitude above the completeness magnitude (from ~ 0.8 to ~ 2 earthquakes per day) and for all recorded earthquakes (from ~ 5 to ~ 10 earthquakes per day). The increase of earthquake frequency during the 2.5 year after Morakot has a probability of 1 (see "Methods" section, Supplementary Fig. S4). The probability is still significant, at the 95% confidence interval, when considering a frequency increase by a factor ~ 1.2-1.5. Except for the regional seismicity rate increases due to the Chi-Chi earthquake, this is the only significant increase in earthquake frequency, over a period of 2.5 year, observed over the investigated period (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013).
Moreover, the increase in earthquake frequency after Morakot is associated with an increase in the b-value from 1.18 ± 0.1 to 1.28 ± 0.1 (Fig. 3a). Both increases have a step-like shape, which lasts for at least 2.5 years (Supplementary Figs. S5, S6) and do not follow an Omori-type inverse law, which describes the temporal evolution of aftershock sequences. Because the seismometer network used to detect earthquakes remained similar in the time period January 2007 to December 2011 (see "Methods" section, Supplementary Fig. S7), we restrict our comparison to the time period from 2.5 years before to 2.5 years after Morakot. More instruments were added, mainly in North Taiwan, at the beginning of 2012, which explain why the frequency of earthquakes remains high in the landsliding zone after 2012, despite a decrease of the b-value towards its pre-Morakot value. Using a  Supplementary  Fig. S1). (c) Gutenberg-Richter law fits over the distributions of cumulative earthquake numbers in the landsliding zone as a function of earthquake magnitude during the 2.5 years before (yellow) and after (green) typhoon Morakot (see "Methods" section).  Fig. 2, Supplementary Fig. S3), we observed that earthquake frequency decreases to pre-Morakot values after about 2.5. The increase of the b-value in the landsliding zone after Morakot is found for different fitting methods of the Gutenberg-Richter law and sampling methods associated with the sliding time window (see "Methods" section, Supplementary Figs. S5, S6) and therefore deemed robust.
In addition, considering all earthquakes during the 2.5 years before and after Morakot gives b-value estimates of 1.17 ± 0.03 and 1.31 ± 0.04 (Fig. 2c), respectively, similar to the values obtained by averaging the temporal b-value signal (Fig. 3a). We observe that the increase in earthquake frequency and b-value in the landsliding zone after Morakot coincides with an increase in the number of shallow earthquakes at depths < 15 km (Fig. 3b,c). Moreover, the increase is more pronounced closer to the surface and reaches a factor 2 to 3 in between 0 and 5 km of depth. We also observe a secondary peak at greater depth, 15 to 25 km. However, the rate of seismic moment release remains low after typhoon Morakot, and potentially lower than before ( Supplementary Fig. S3). Crucially, earthquakes outside the landsliding zone do not show a significant temporal increase of both their frequency and b-value after Morakot.
Seismicity changes are not associated to earthquake clustering. These temporal changes in earthquake statistics after Morakot are determined from an undeclustered earthquake catalog. However, using a declustered catalog (see "Methods" section) also leads to similar changes in earthquake frequency and in its associated probability, b-value and depth-distribution after Morakot (Fig. 3c, Supplementary Figs. S9, S10). Moreover, this result is not limited to a single parametrization of the declustering algorithm and appears robust for all the acceptable combinations of declustering parameters. For instance, the increase in the number of earthquakes after Morakot at 15 to 25 km of depths is fully removed when considering any of the resulting declustered catalog, while the shallow peak remains (Fig. 3c). This demonstrates that these seismicity changes at shallow depths are not associated with triggering processes by large mainshocks, or by other earthquakes. In addition, based on comparisons between the observed earthquake catalog and synthetic catalogs that share the same average properties, we demonstrate that the observed changes in earthquake frequency and b-value after Morakot, in the landsliding zone, depart statistically from random temporal changes ( Supplementary Fig. S12). www.nature.com/scientificreports/ Spatial changes in seismicity after Morakot. To assess how the chosen delimited area affects our results, we compute the change in the spatial pattern of the frequency and b-value of shallow (< 15 km) earthquakes from before to after Morakot. For this, we use a sliding window in space with a radius of 30 km, which allows us to detect large-scale features not affected by small sub-samples of events (see "Methods" section). It is applied separately to earthquakes in the 2.5 years before and after Morakot, respectively. Consistent with the temporal evolution of earthquake frequency, results show an increase in the number N of earthquakes after the typhoon over the landsliding zone ( Fig. 4, Supplementary Fig. S8). This increase is observed for earthquakes with magnitudes above the completeness magnitude and also for all recorded earthquakes. It is not limited to the vicinity of the landsliding zone, but it also occurred in northeast Taiwan. However, outside the landsliding zone, increases in the b-value appear not to be associated to increases in the number of earthquakes. For instance, North-East Taiwan displays a significant increase in the number of earthquakes associated to a decrease of the b-value. We note that the spatial correlation between earthquake statistics change and the landsliding zone is less well resolved and less robust (see "Methods" section) than the temporal correlation.

Discussion
Previous intense erosional events associated with the Chi-Chi earthquake (1999), and typhoons Herb (1996) and Toraji (2001) did not induce any detectable change of seismicity. This may be because these events had less erosion, the total volume of ~ 0.45 km 3 of landslides triggered by Chi-Chi earthquake 25 being by far the largest, and because most landslides triggered by these events deposited debris distant from rivers 14 . Besides, the high rate of aftershocks after Chi-Chi prevents detection of a change in earthquake frequency or b-value associated with erosion due to Chi-Chi or Toraji, if there was any. We also note that the changes in earthquake frequency and b-value are not perfectly synchronous with typhoon Morakot and occur about 1-2 months after the typhoon. However, the time resolution of our temporal analysis does not allow us to detect with confidence changes occurring with a period less than about 100 days (see Supplementary Fig. S5b). The duration of the transient increase in earthquake frequency and b-value occurring after typhoon Morakot lasts about 2.5 years.
Non-erosional causes of the observed changes in earthquake statistics in SW Taiwan are possible but appear less likely. Earthquakes result from the stresses induced by tectonics, in general, but also by other earthquakes 26 . The 4th of March 2010, the Mw 6.3 Jiashian earthquake occurred within the landsliding zone, close to its western limit, and was followed by many aftershocks. However, this earthquake and most aftershocks are located at 15-25 km depth ( Supplementary Fig. S11), as shown by seismological records 27 and further confirmed by the declustering process, which mostly removes events below 15 km depth, (Fig. 3, Supplementary Fig. S10). Therefore, the Jiashian earthquake does not affect shallow earthquake statistics and cannot be considered responsible for the increase of the b-value after Morakot (Supplementary Fig. S10).
Hydrological triggering of seismicity after the heavy rainfall during typhoon Morakot, either by surface loading or by pore pressure diffusion 28 , could be an alternative mechanism. Indeed, on the east coast of Taiwan, where landslide erosion was limited, borehole water levels rose by 4 m after typhoon Morakot, and a volumetric contractional strain was observed 29,30 . However, both signals decayed in ~ 6 months and their amplitudes and temporal evolutions do not depart from the mean seasonal trends observed from 2006 to 2011. Moreover, pore pressure diffusion along permeable faults after large rainfall events generally leads to episodic increases of seismicity 31 and not to prolonged changes as observed in our case.
The temporal and spatial collocation of intense landsliding triggered by typhoon Morakot and the observed increase of shallow earthquake frequency and b-value suggest a potential mechanistic link. Direct physical modeling of the impact of erosion during and after Morakot on seismicity is beyond our reach because the location, magnitude and rate of sediment export from the landsliding zone are not reliably constrained. Despite this, simple elastic models show that large erosional events with rapid sediment export can induce static stresses at depth, sufficient to modulate tectonic stresses on the shallower (< 5-10 km) parts of faults 10 . Removing 20 to 50% of the landslide volume in 2.5 years, equivalent to about 2 to 5 cm of average erosion over the landsliding zone, would lead to a Coulomb stress increment on a nearby thrust fault of about 0.1 10 -1 to 0.25 10 -1 bar at 5 km depth 10 , using a thrust dip angle of 30° and an effective friction of 0.6. These stress increments are roughly similar to the ones induced by seasonal hydrologic loading in the Himalaya, ~ 0.2-0.4 × 10 -1 bar, which are suggested to lead to a seasonal modulation of earthquake frequency 32 .
In addition, spring-slider models 33 and 2D elasto-dynamic models of seismogenic faults 34 with rate-and-state friction laws 35,36 show that the rate of seismicity can increase linearly or more than linearly due to a positive, step-like stress perturbation. It is also observed that shallow earthquakes generally have smaller magnitudes and larger b-values than deeper earthquakes 37 possibly because they nucleate at lower differential stresses 38 . Our observations suggest that the intense and prolonged sediment export and surface unloading after typhoon Morakot could have acted as a progressive increase of stresses, with a complex spatio-temporal evolution, on the shallow parts of underlying thrust faults 10 , giving rise to a similar increase of shallow earthquake frequency and b-value in Southwest Taiwan.
The proposed erosional-driven mechanism to explain the changes in seismicity after Morakot offers new perspectives on the links between climate, erosion and tectonics 4,5 at the time scale of elementary processes. While numerous studies have shown that earthquakes and storms can trigger landslides 2,14,20,25 , this is to our knowledge the first potential evidence of the short-term and ongoing influence of erosion on seismicity. Because the mechanical link between erosion and stresses is through crustal elasticity 10 , crustal deformation is sensitive not only to extreme weather-driven erosion but also to the cumulative effects of smaller but numerous erosion events. More frequent extreme rainfall, for instance under a warmer climate 39 , could result in accelerated sediment transport 19 and in turn in changes in the rate of shallow earthquakes, similar to the modulation of seismicity by surface water load variations 13 . Over the duration of the transient increase in seismicity, the shallow but smallmagnitude seismicity induced by erosion is not likely to trigger new landslides, and the seismic moment rate would not be significantly affected. Hence, we do not expect a significant feedback of this additional deformation on erosion. However, the impact of this short-term transient signal on the seismicity and erosion over longer time scales, such as the duration of a seismic cycle, remains to be investigated. Our results therefore call for a new generation of process-based models coupling landscape dynamics 21 and fault dynamics 34 at scales relevant to natural hazards and societal issues. In these models, storms, floods, mass wasting, river sediment transport, redistribution of water, elastic stress transfer, seismicity and seismic wave propagation should all be represented to account for the complexity of the links between climate, erosion and tectonics. To cover most of the islands we mosaicked multiple cloud-free pre-event (01/14, 05/08, 05/09, 05/10, 06/06) and post-event (08/17, 08/19, 08/21, 08/28, 08/30, 09/06) images taken in 2009. For parts of the inventory, especially east of the main divide, landslides were first mapped automatically and then edited manually. For both approaches, the scar, runout and deposit areas are not differentiated. We did not consider debris flow transport areas and excluded gentle slopes (< 20°) from mapping to avoid confusion with human activity. Special attention was given to the separation of individual landslides, which had common transport or deposit areas but independent initiation points 41 . The robust and conservative estimation of landslide surface area is especially important for the estimation of landslide volume 42 . estimation of landslide volume. The landslide volume was estimated based on landslide area, following the method of ref. 16 and 41. Briefly, we assumed constant size ratios between scar and deposit areas of 1.1 and 1.9 for mixed and bedrock landslides, respectively 42 . Then, we converted the scar area into volume using a power law with different prefactor ( α ) and exponent ( γ ) for mixed and bedrock landslides, with α = 0.146 and γ = 1.332 for A scar < 1e5 m 2 and α = 0.234 and γ = 1.41 for A scar > 1e5 m 2 , respectively 41 . Recent studies have proposed a regional scaling relationship for southern Taiwan In 2010, suspended sediment discharge decreases to pre-Morakot values, partly due to the absence of large water discharge events. Timescales for bedload sediment export is probably much larger, between a few years to decades 21 and potentially centuries. Most rivers connected to the landsliding zone experienced rapid and intense sediment aggradation, with 10-100 m of sediment thickness, in response to Morakot typhoon. To our knowledge, no study has been published on the regional evolution of this bedload sediment mass. The estimated total volume 1. www.nature.com/scientificreports/ time evolution of earthquake statistical properties. We use a temporal sliding window to subsample the earthquake catalogue and to assess time variations of earthquake statistical properties. Because the uncertainty of the b-value strongly depends on the number of considered events, and because larger samples give better estimates, we use a sliding window of N = 1,001 events to prevent undue statistical bias. The window is centred on a given earthquake and the corresponding b-value and earthquake frequency are determined for the 500 earthquakes that occurred immediately before and immediately after this event. Earthquake frequency is computed by dividing N by the temporal length of the sliding window. Among the 1,001 earthquakes, only those with magnitude above M c are used to estimate the b-value. We assess the effect of changing the sampling method of the temporal window, by considering an a priori value of M c = 2.25 , and considering a fixed number of earthquakes with a magnitude > M c ranging from N = 101 to N = 1,001 ( Supplementary Figs. S5, S6). We find a consistent increase of b-value in the 2.5 years after Morakot for N ranging from 101 to 501, as all these earthquakes occurred within a time window duration lower than 2.5 years. For N = 751 or N = 1,001 the window duration is equal to or greater than 2.5 years (about 3 years for N = 1,001), which leads to over-smoothing of the signal, and prevents detection of changes occurring on shorter time-scales. It is also notable that changing the fitting method from maximum likelihood to least-square, which is generally considered less reliable, does not significantly change the relative variation of the b-value in time (Supplementary Fig. S5). This includes the changes occurring after Morakot.

Methods
Probability of earthquake frequency change with time. Following ref. 48 , we compute the probability P that the earthquake frequency increases by a factor greater than r between a period 1 and 2, where Ŵ(n, x) = ∞ 0 e −t t n−1 dt is the incomplete Gamma function, N is the number of earthquake over a certain period t and = N/�t is the earthquake frequency. The subscripts 1 and 2 refer to the time periods 1 and 2. We first apply this approach to determine the probability of earthquake frequency change between the 2.5 year before and after Morakot, in the landsliding zone considering only earthquakes with a depth shallower than 15 km (Supplementary Fig. S4a). We find a probability 1 for an increase of earthquake frequency (i.e. with r > 1 ), when considering all the earthquakes or only earthquakes above the completeness magnitude. The 90% confidence interval, 0.05 < P < 0.95 , of an earthquake frequency change is found for 1.23 < r < 1.40 , when considering only the earthquakes above the completeness magnitude, and 1.45 < r < 1.52 for all magnitudes. This demonstrates that the change of earthquake frequency after Morakot is significant with a ratio of 1.23, at least. We then apply the same analysis to the entire catalog by using a double sliding window of period 2.5 years after and before the center time ( Supplementary Fig. S4b). The largest increase in earthquake frequency, with a factor close to 5 for a probability inside the confidence interval, is Chi-chi earthquake. Except for Chi-Chi, the time period including typhoon Morakot and Jiashian earthquakeis the only one associated with a significant and positive change of earthquake frequency over a period of 2.5 years since 1994. However, as demonstrated in this manuscript, Jiashian earthquake did not lead to a significant increase in the number of earthquakes at shallow depths ( Fig. 3; Supplementary Figs. S9, S10, S11, S13). Spatial variations of earthquake statistical properties. We use a common spatial sliding window to subsample the earthquake catalogue and to assess variations of earthquake statistical properties in space 49 between the 2.5 years before and after typhoon Morakot. We use a radius of 30 km for the sliding window, which enables sampling of a sufficiently large number of earthquakes at a length scale that is smaller than that of the landslide zone of 7,000 km 2 . The sampling window effectively corresponds to a disk shape extending to a depth of 15 km. The b-value is determined for the recorded earthquakes in each disk volume by maximum likelihood estimation. This method has an inherent statistical bias as the number of sampled earthquakes changes significantly depending on the local rate of seismicity. The minimum number of events for the determination of the b-value is set arbitrarily at 50. Our method can give rise to small-scale shapes in the maps of b-value, such as disk and rod shapes, that are not the focus of this study. We use a disk shape kernel with a radius of 15 km, convolved with the initial b-value map to blur the mapped patterns and to isolate features with longer wavelengths. More sophisticated methods exist to compute spatial variations of b-value, including a penalized likelihood-based method 50 and a distance-dependent sampling algorithm 51 . Declustering and seismicity changes after Morakot. Several studies investigating potential earthquake frequency changes use declustered earthquake catalogs. However, declustering is an ill-posed problem that does not have a unique solution 52 and that will lead to method-dependent results. Yet, to test the impact of potential earthquake clustering on our results, we have applied the traditional Reasenberg declustering algorithm 53 , obtained from the ZMAP toolbox 54 , to the CWB earthquake catalog. This deterministic algorithm aims to remove earthquake sequences, defined as chains of connected earthquakes in space and time, leaving only the initial earthquake in a given sequence. We use standard parameters 52 adapted for Taiwan with τ min = 1 day the minimum look-ahead time for for building clusters when the first event is not clustered, τ max = 10 days, the maximum look-ahead time for building clusters, p clust = 0.95 , a confidence probability, x meff = 2.0 , the effective lower cutoff magnitude chosen here to be consistent with the CWB catalog, x k = 0.5 , the increase in lower cutoff magnitude during clusters, and r fact = 10 , the number of crack radii surrounding each earthquake within new events considered to be part of the cluster. We emphasize here that declustering has been performed Scientific RepoRtS | (2020) 10:10899 | https://doi.org/10.1038/s41598-020-67865-y www.nature.com/scientificreports/ over the entire catalog, without any regional selection of the seismicity (i.e. not only in the landsliding zone and not only above 15 km of depth), to prevent potential declustering biases associated to earthquake censoring. Supplementary Figure S9 shows the influence of declustering on earthquake frequency, b-value, depth distribution and the probability of earthquake frequency change after Morakot in the landsliding zone (< 15 km depth). We observe that declustering mainly leads to a decrease of earthquake frequency over the 2.5 years before typhoon Morakot, and to very minor changes after. We find that declustering even slightly enhances the probability of an increase in earthquake frequency after Morakot. In addition, declustering does not significantly changes the time variation in b-value nor the depth-distribution of seismicity. We note that the Reasenberg declustering approach is well-suited to remove earthquake sequences that lead to significant changes (e.g. an aftershock sequence after a large mainshock) and not to identify potential triggered events during quiet periods associated with a constant and low frequency of earthquakes 48 . Under this potential limitation, this analysis demonstrates that earthquake clustering is not the reason for the observed changes after typhoon Morakot.
In addition, we are confident that large and deep mainshocks occurring in the landsliding zone, below 15 km, such as Mw 6.3 Jiashian earthquake in 2010, are not the cause for the observed seismicity changes. Indeed, declustering is efficient to remove most of the aftershocks caused by the Jiashian earthquake, which appear to be concentrated at depth between 15 and 30 km (Supplementary Fig. S10). To confirm this inference, we vary the values of the declustering parameters using a range of acceptable values around the reference values 48,55 : τ min = 0.5−2.5 days, τ max = 3−15 days, p clust = 0.90−0.99 , x k = 0−1 and r fact = 5−20 . We use a random sampling of the parameters space within these ranges, using a Monte-Carlo approach, to define 50 declustering parameter combinations. We obtain similar results, in terms of earthquake depth-distribution ( Supplementary  Fig. S9) and temporal variations of earthquake frequency (Supplementary Fig. S10) as in the reference declustering model. This confirms the robustness of the observed changes in shallow seismicity after Morakot that are not related to earthquake interactions. Seismicity temporal variation: random hypothesis versus significant regional seismicity change versus local earthquake interactions. Having a declustered catalog should theoretically imply that all the earthquakes are seismically independent from each other. In turn, earthquakes should be randomly distributed in time, unless a non-seismic process triggers them. Here we test this hypothesis by comparing the observed changes, in earthquake frequency and b-value after Morakot, obtained from the "true" declustered catalog with those from 200 "synthetic" earthquake catalogs over the period 2006-2015. Each synthetic catalog is generated using the exact same earthquakes as in the true declustered catalog, including their magnitudes, but the time of occurrence of each earthquake is randomly sampled over the period of interest (2006-2015) using a time step of 1 s. It results that each synthetic catalog has the exact same magnitude distribution and average earthquake frequency than the true one, but the temporal distribution of earthquakes is randomly distributed. For each catalog (true or synthetic), the probability of frequency change and the change in b-value after Morakot are computed by comparing the earthquakes occurring in the 2.5 years after and before Morakot (see Supplementary Fig. S12). The probability of a frequency change after Morakot for all the synthetic catalogs drops around a ratio 0.8-1.2, centered around 1, meaning there is no significant frequency change. This clearly departs from the frequency change observed using the true catalog that drops around a ratio of 1.3-1.6. In addition, the change in b-value of the true catalog is significantly greater than the changes of all the synthetic catalogs. Overall, these results mean that both the frequency change and the b-value increase after Morakot are robust features, that depart from random changes. In addition, because we performed these tests using the declustered catalog, which should only include independent earthquakes, this means that these robust and non-random changes should not be associated to the occurrence of large mainshocks. We also note that potential earthquake interactions in a local subset of the landsliding zone 56 cannot explain the increase in earthquake frequency after Morakot ( Supplementary Fig. S13).

Data availability
The datasets used in this study are provided along with the manuscript. The earthquake catalogue is also accessible through the Taiwan Central Weather Bureau (https ://gdms.cwb.gov.tw). The rainfall data is also accessible through Data Bank for Atmospheric Research at the Taiwan Typhoon and Floods Research Institute (https :// www.narla bs.org.tw/en). The landslide dataset is also available by request to O.M., P.M. and N.H.

code availability
The Matlab code used to perform the analysis and the figures of this manuscript, as well as the data needed to process the figures of this paper, are provided here: https ://zenod o.org/recor d/37783 55. The Reasenberg declustering algorithm is accessible through the ZMAP set of tools (https ://www.seism o.ethz.ch/en/resea rch-and-teach ing/produ cts-softw are/softw are/ZMAP/ and https ://githu b.com/Celso Reyes /zmap7 ).