# Landslide hazard cascades can trigger earthquakes

Zhen Zhang, Min Liu, Yen Joe Tan, Fabian Walter, Siming He, Malgorzata

Chmiel, Jinrong Su

► To cite this version:

Zhen Zhang, Min Liu, Yen Joe Tan, Fabian Walter, Siming He, et al.. Landslide hazard cascades can trigger earthquakes. Nature Communications, 2024, 15 (1), pp.2878. 10.1038/s41467-024-47130-w. hal-04745947

HAL Id: hal-04745947

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

Submitted on 24 Oct 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.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Landslide hazard cascades can trigger earthquakes

Received: 31 July 2023

Accepted: 21 March 2024

Published online: 08 April 2024

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

Check for updates

Zhen Zhang<sup>1, ∗</sup>, Min Liu<sup>1</sup>, Yen Joe Tan<sup>1, †</sup>, Fabian Walter<sup>2</sup>, Siming He<sup>3</sup>,
Małgorzata Chmiel<sup>2,4</sup> & Jinrong Su<sup>5</sup>

While earthquakes are well-known to trigger surface hazards and initiate hazard cascades, whether surface hazards can instead trigger earthquakes remains underexplored. In 2018, two landslides on the Tibetan plateau created landslide-dammed lakes which subsequently breached and caused catastrophic outburst floods. Here we build an earthquake catalog using machine-learning and cross-correlation-based methods which shows there was a statistically significant increase in earthquake activity (local magnitude ≤ 2.6) as the landslide-dammed lake approached peak water level which returned to the background level after dam breach. We further find that ~90% of the seismicity occurred where Coulomb stress increased due to the combined effect of direct loading and pore pressure diffusion. The close spatial and temporal correlation between the calculated Coulomb stress increase and earthquake activity suggests that the earthquakes were triggered by these landslide hazard cascades. Finally, our Coulomb stress modeling considering the properties of landslide-dammed lakes and reservoir-induced earthquakes globally suggests that earthquake triggering by landslide-dammed lakes and similar structures may be a ubiquitous phenomenon. Therefore, we propose that earthquake-surface hazard interaction can include bidirectional triggering which should be properly accounted for during geological hazard assessment and management in mountainous regions.

Mass movements such as landslides, avalanches, and debris flows are regarded as both major hazards and key forms of geomorphic evolution in mountainous regions<fnref n="1" />⁻³. These surface processes can be triggered by many factors including rainfall⁴,<fnref n="5" />, snow melt⁶,⁷, and human activities⁸,⁹. In seismically active regions, the risk is also heightened since large earthquakes can induce widespread mass wasting with devastating effects¹⁰,¹¹. When these mass movements block rivers and form landslide-dammed lakes (LDLs), subsequent dam breaches may cause catastrophic outburst floods which significantly increase the degree and scope of the disaster¹²,¹³. For example, the 1933 M7.5 Diexi earthquake in the eastern Tibetan plateau directly caused ~7000

fatalities and triggered many large landslides, some of which dammed rivers. These LDLs breached in the following days, generating catastrophic outburst floods which then resulted in >2500 additional fatalities¹⁴. Therefore, the response of surface processes to earthquakes has been studied extensively, in particular the failure mechanisms of earthquake-triggered landslides¹¹,¹², the quantification and prediction of the spatial distribution of mass movements¹³,¹⁵, and the geomorphic evolution after earthquakes in mountainous areas¹⁶,¹⁷, to improve our ability to mitigate the impact of these complex hazard cascades.

To date, studies about the interaction between earthquakes and mass movements have focused mainly on how earthquakes trigger

<footnote n="1">¹Earth and Environmental Sciences Programme, Faculty of Science, The Chinese University of Hong Kong, Hong Kong S.A.R., China. ²Swiss Federal Institute for Forest, Snow and Landscape Research, Zürich, Switzerland. ³State Key Laboratory of Mountain Hazards and Engineering Safety, Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu, China. ⁴Géoazur, OCA, Campus Azur du CNRS, Sophia Antipolis, Nice, France.</footnote>

<footnote n="5">⁵Earthquake Monitoring Center, Sichuan Earthquake Administration, Chengdu, China. ✉ e-mail: zhenzhang@cuhk.edu.hk; yjtan@cuhk.edu.hk</footnote>
mass movements and initiate hazard cascades, while a potential earthquake response to mass movements and their hazard cascades has received less attention. However, various studies have shown that processes which change the stress state of the Earth's crust by surface loading and/or fluid diffusion<sup>18,19</sup>, such as earth tides<sup>20</sup>, water storage behind dams<sup>21</sup>, hydraulic fracturing<sup>22</sup>, waste fluid disposal<sup>23</sup>, and lake filling<sup>24-27</sup> can trigger/induce earthquakes when the stress increase on neighboring active faults exceeds a critical threshold<sup>28,29</sup>. Since hazardous mass movements and their complex hazard cascades are often accompanied by sediment redistributions and changes in water storage and associated surface loads<sup>3,13</sup>, it is conceivable that the resulting stress changes can trigger earthquakes<sup>30-33</sup>. This would have important implications for hazard management in seismically active mountain regions.

Here, we use the 2018 Baige landslide hazard cascades (Fig. 1) on
the Tibetan plateau as a case study to investigate earthquake triggering

by LDLs. We quantify the spatiotemporal evolution of seismicity and
stress changes on a surrounding fault system due to this hazard cas-
cade. We further use Coulomb stress modeling based on a global LDL
database to explore the potential for LDLs to trigger earthquakes.

Results

2018 Baige hazard cascades

At 22:05 China Standard Time (CST) on 10 October 2018, a massive landslide occurred at the Baige village on the Tibetan plateau<sup>34-38</sup> (Fig. 1 and S1) and deposited 2.5 × 10<sup>7</sup> m<sup>3</sup> of sediments into the Jinsha River<sup>35</sup>. This created a landslide dam with ~1500 m length, >450 m width, and ~60 m height<sup>35-37</sup> (Fig. S1). With an upstream inflow of 1680 m<sup>3</sup>/s (ref. 35), the water level of the LDL gradually increased from ~2880 m asl to a peak of ~2932 m asl at 00:45 CST on 13 October (Fig. 2c), corresponding to inundation up to ~45 km upstream and a peak volume of ~2.9 × 10<sup>8</sup> m<sup>3</sup> (refs. 35–37). As the water naturally

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

**Fig. 1** | **Study site and regional seismicity.** a Distribution of seismic stations operating since January 2014 (black triangles) and May 2018 (cyan triangles). Red box marks the specific study region where we developed the earthquake catalog. b Seismicity map from January 2014 to January 2023 (gray filled circles) showing strike-slip focal mechanisms (red/white circles) available from the National Earthquake Data Center<sup>39</sup> plus earthquakes that occurred during the week (10 to 16 November 2018) when the second landslide-dammed lake (LDL) reached peak

water level (red filled circles). Red box marks the region shown in (d). c Magnitude-frequency distribution of earthquakes from May 2018 to January 2023 within 10 km of the LDLs. Magnitude of completeness (Mc) and a- and b-values were calculated using the ZMAP software<sup>40</sup>. d Distribution of earthquakes from May 2018 to January 2023 (gray filled circles) within 10 km of the LDLs plus earthquakes that occurred during the week (10 to 16 November 2018) when the second LDL reached peak water level (red filled circles).
<figure><img src="image_5.png" /><figcaption>**Fig. 2 | Earthquakes triggered by the 2018 Baige landslide-dammed lake.**</figcaption></figure>

a Spatial distribution of Coulomb stress change (ΔCFS) at depth of 4.5 km corresponding to -0.5 km asl on 13 November 2018 and earthquakes from 10 to 16 November 2018 (green filled circles) within 10 km of the landslide-dammed lakes (LDLs). ΔCFS is from the combined effect of direct gravitational loading and pore pressure diffusion (Fig. 2d). b Cumulative number of earthquakes with local

magnitudes greater than or equal to the magnitude of completeness (Mc = 1.1) and weekly seismicity rate over a 5-year period with timings of the two 2018 Baige landslides (black arrows). Gray bar marks the period shown in (c and d).

c Earthquake magnitudes and measured LDL water level³⁶,³⁷ over a fifty-day period. Gray dashed line marks Mc. d Temporal evolution of ΔCFS at point P (a).

overflowed the west of the dam, the landslide dam breached rapidly. The outburst flood's discharge reached a peak of 10,000 m³/s at 06:00 CST on 13 October while the LDL's water level rapidly decreased. Eight hours later at 14:00 CST, the discharge was reduced to that of the upstream lake inflow. However, due to the residual dam, the water level upstream of the dam (-2893 m asl) remained higher than before the landslide.

About three weeks later at 17:40 CST on 3 November 2018, materials at the trailing edge of the first landslide's scar suddenly failed and the resulting landslide deposited an additional 8.7 × 10⁶ m³ of sediments on the residual dam formed by the first landslide³⁵⁻³⁷. This formed a new >96 m high landslide dam and a new LDL³⁵⁻³⁷. To mitigate the potential impact of another outburst flood, the Chinese government trenched an artificial spillway of 220 m length, 15 m depth, and >3 m width on the west side of the dam³⁶,³⁷. Water entered the spillway at 04:45 CST on 12 November marking the initiation of a second outburst flood and the LDL reached a peak water level of 2957 m asl at 13:40 CST on 13 November (Fig. 2c), corresponding to inundation up to -70 km upstream and a peak volume of ~5.8 × 10⁸ m³. At 18:00 CST on 13 November, the outburst flood's discharge reached a maximum of 31,000 m³/s (refs. 35⁻37). Fourteen hours later at 08:00 CST on 14 November, the discharge decreased to the upstream inflow of ~600 m³/s. Like for the first LDL outburst, the water level upstream of the dam (2906 m asl) remained higher than before the landslide due to the residual dam (Fig. 2c).

## Seismicity rate increase after Baige LDLs' formation

Using continuous seismic data recorded by 13 nearby seismic stations from January 2014 to January 2023, we apply machine-learning and

cross-correlation-based methods to build an earthquake catalog which contains ~3970 earthquakes with documented magnitude within a region of ~1° × 1° (Fig. 1b). 244 of these earthquakes within 10 km of the LDLs occurred after May 2018 when a denser seismic network began operation. It is well-established that earthquakes can trigger surface hazards and initiate hazard cascades⁸⁻⁷,¹². However, the Baige landslides were not triggered by earthquakes³⁴⁻³⁸ as no earthquakes with magnitude $M_L \ge 3.5$ were recorded within 50 km in the 3 preceding years³⁹ (Figs. 1 and 2). Instead, in the week (10 to 16 November) when the second LDL approached its peak water level, 61 earthquakes with local magnitude up to ~2.6 occurred within 10 km of the LDLs.

As the second LDL approached peak water level, the number and magnitude of earthquakes started increasing as the water level rised and then peaked together with the water level peak (Fig. 2 and S2). Subsequently, as the LDL water level decreased following the dam breach, the earthquake activity gradually decreased back to the background rate. 45 of the 61 events in this earthquake sequence occurred before the largest magnitude ($M_L$ 2.6) earthquake, and other $M_L$ > 2.6 earthquakes during our observation period all had fewer than 4 aftershocks in the following week (Fig. S3). Therefore, this earthquake sequence is unlikely to be primarily an aftershock sequence of the largest magnitude event. We further verified that only 5 earthquakes (all with $M_L$ < 2.0) occurred within 60 km of the LDLs besides this earthquake sequence as the second LDL approached peak water level (Fig. 1b), and the largest earthquake that occurred within 60 km of the LDLs in the 6 months before this earthquake sequence has $M_L$ 2.7 (Fig. S4). Therefore, the earthquake sequence within 10 km of the LDLs is unlikely to have been triggered by surrounding large earthquakes.
We estimate the magnitude of completeness (Mc) to be 1.1 (ref. 40; Fig. 1c) for the earthquakes within 10 km of the LDLs. There are 62 earthquakes between May 2018 and January 2023 with local magnitude $M_L \ge 1.1$, 16 of which occurred in the week when the second LDL approached its peak water (Table S1). We then perform declustering using the Reasenberg method⁴¹, leaving us with 54 events (Fig S5), 9 of which occurred in the week the second LDL approached its peak water level (Fig. S5, Table S1). In comparison, there were only 8 earthquakes in the previous 27 weeks (seismicity rate of ~0.3 events/week) and 37 earthquakes in the subsequent ~224 weeks (seismicity rate of ~0.2 events/week) in this region (Table S1). For each of $\ge$ Mc and declustered catalogs, we consistently find that the increase in seismicity rate is statistically significant at a > 99% level based on both the statistic P (ref. 42) and statistic Z (ref. 43) tests (Methods). We further confirm that the increase in seismicity rate is also statistically significant at a >99% level (Table S1) based on both the improved statistic based on Poisson probability⁴⁴ and an empirically derived statistic⁴⁵ (Fig. S6) tests (Methods) which were developed to determine statistically significant changes in earthquake rate under small background rate. Hence, the significant increase in seismicity rate starting when the second LDL approached peak water level (Fig. 2 and S2) suggests that the two are related.

## Stress changes caused by Baige LDLs

We further investigate whether the LDLs could have triggered the earthquakes by modeling the Coulomb failure stress changes (ΔCFS) on the surrounding fault systems. In our study region, there are three local faults Boluo-Tongmai, Jinshajiang, and Gangtuo-Yidun faults

(from west to east) with similar surface traces of ~N50°W (ref. 46; Fig. 1a) though the exact fault dip and slip directions are unknown. Nevertheless, since January 2014, four M > 3.5 earthquakes occurred at distances between 16 and 35 km from the LDLs and have similar left-lateral strike-slip focal mechanisms with strikes that are also ~N50°W and dips of ~75°NE (ref. 39; Fig. 1a). In addition, we find that the first motions recorded by surrounding seismic station of the two largest earthquakes ($M_L$ 2.6 and 1.9) as the second LDL approached peak water level are generally consistent with the left-lateral strike-slip focal mechanisms of these surrounding large earthquakes (Fig. S7). Therefore, we take these as our receiver fault geometry for the ΔCFS modeling (Methods). Furthermore, the receiver faults are assumed to be planar and located at 4.5 km depth which is deeper than 95% of our observed seismicity (Fig. S8) and corresponds to ~0.5 km asl because ~80% of our study area is at elevation between 4 and 5 km asl³⁵.

We find that as the LDLs' water level increased following the landslides, the ΔCFS in the seismic region generally increased due to both the direct gravitational loading of the LDL and pore pressure diffusion, while ΔCFS decreased when the LDLs' water level decreased following the dam breaches (Fig. 2 and S9). The larger and longer-lasting LDL generated greater ΔCFS. The ΔCFS due to the first, smaller, and shorter-duration LDL peaked at 0.007 MPa which is smaller than the stress increase of >0.01 MPa typically associated with stress triggering of earthquakes (refs. 47,48; Fig. 2d). This could explain why there was no significant increase in seismicity rate associated with the first LDL's formation (Fig. 2c). In comparison, about 3 weeks later on 3 November, the second landslide formed a larger LDL. The ΔCFS due to this second LDL peaked at 0.024 MPa when the LDL's water level

<figure><img src="image_6.png" /><figcaption>**Fig. 3** | General properties of landslide-dammed lakes (LDLs). **a** Frequency distribution and cumulative distribution functions of the height of reservoirs documented to have induced earthquakes⁴⁹ and landslide-dammed lakes³⁰ worldwide. **b** Frequency distribution and cumulative distribution function of the duration of landslide-dammed lakes worldwide³⁰. **c** Temporal evolution of ΔCFS and pore</figcaption></figure>

pressure at point P (d) for a 30-m deep LDL. **d** Spatial distribution of ΔCFS at depth of 4.5 km corresponding to ~0.5 km asl, 11 days after the formation of a 30-m deep LDL. ΔCFS is from the combined effect of direct gravitational loading and pore pressure diffusion.
approached its peak, coinciding with maximum earthquake activity (Fig. 2c). Although the subsequent dam breach lowered the LDL's water level which decreased the ΔCFS due to direct gravitational loading, the pore pressure diffusion effect kept the ΔCFS elevated which might explain why earthquake activity continued over the next few days (Fig. 2c).

While direct loading resulted in a decrease in ΔCFS directly below the LDL, it resulted in a slight increase in ΔCFS in surrounding regions (Fig. S9a). On the other hand, fluid diffusion significantly increased the pore pressure around the LDL, and the ΔCFS in regions closer to the LDL were larger (Fig. S9b). Therefore, due to the combined effect of direct loading and pore pressure diffusion, the ΔCFS around the LDL increased, especially on its east side. We find that ~90% of the earthquakes locate within regions of positive ΔCFS (Fig. 2a). The close spatial and temporal correlation between the calculated positive ΔCFS and the statistically significant increase in earthquake activity suggest that the earthquakes were likely triggered by the LDL.

## Earthquake triggering by LDLs

Considering the fundamental mechanism of earthquake triggering by LDLs and reservoir-induced earthquakes is similar, we further explore the potential for other LDLs to trigger earthquakes. Based on HiQuake, currently the most complete and up-to-date freely available database of ~226 cases of reservoir-induced earthquakes spanning the period 1933-2019 globally⁴⁹, we find that the minimum and median heights of ~190 reservoirs with documented height that induced earthquakes are ~13 and 110 m, respectively (Fig. 3a). In comparison, from a global database of ~410 LDLs with dams >1 million m³ in volume reported worldwide spanning the period 1900-2018 (ref. 50), we find that of the ~300 LDLs with documented dam height, ~73% and ~10% were higher than the minimum and median heights of reservoirs which induced earthquakes, respectively (Fig. 3a). Therefore, some LDLs globally are of similar heights as reservoirs which induced earthquakes.

We further calculate the stress response of surrounding faults in our study area to an LDL with a depth of 30 m which is the median dam height in the global LDL database⁵⁰ (Fig. 3b). In this case, the ΔCFS at a given point (Figs. 2a and 3d) would exceed 0.01 MPa after ~11 days due to direct loading and pore pressure diffusion if the dam does not breach (Fig. 3c). The stress responses to this typical LDL at different

sites can vary (Fig. 3d) and depend on factors such as the fault types⁵¹ (Fig. S10), geometries (Fig. S11), and depths (Fig. S11). Nevertheless, our modeling shows that both the fault types and geometries primarily affect stress distribution rather than amplitude of ΔCFS (Figs. S10 and S11), and stress changes caused by the LDLs are larger at shallower depths (Fig. S11). Considering ~64% of the ~145 LDLs in the global database with documented duration lasted more than 11 days before dam breach (Fig. 3b) with no correlation with dam height (Fig. S12), our modeling suggests that other LDLs globally have the potential to trigger earthquakes if there are critically stressed faults nearby.

## Discussion

While gravitational loading can cause Coulomb stress to decrease in some areas, pore pressure always causes an increase in ΔCFS with its amplitude gradually decreasing away from the LDLs (Fig. S9). For earthquakes that occurred in regions where gravitational loading resulted in a decrease in Coulomb stress, it is clear that pore pressure is the primary mechanism triggering these events, though we neglect the coupling effect between pore pressure diffusion and gravitational loading²⁴⁻²⁷. However, for those earthquakes that occurred in regions where gravitational loading resulted in an increase in Coulomb stresses, both gravitational loading and pore pressure may jointly control the triggered seismicity. After the second LDL formed, while the water level rises and seismicity begins to occur from 10 to 13 November (Fig. S2), the relative contribution to ΔCFS of gravitational loading decreases gradually from 0.45 to 0.41 (Fig. S13). Subsequently, after the dam breached on 13 November which is the day of peak seismicity rate, the water level rapidly decreased and the relative contribution to ΔCFS of gravitational loading decreases sharply from 0.41 to 0.19 so pore pressure becomes the dominant triggering mechanism. However, the relative contribution to ΔCFS of pore pressure depends on the assumed hydraulic diffusivity (Fig. S14). In addition, since the pore pressure diffusion is from a long winding river (LDLs) instead of a point source, the seismicity do not display obvious migration pattern with time (Fig. S15).

The region surrounding the 2018 Baige LDLs is seismically not very active which allowed us to identify the sudden, sharp increase in earthquake activity and link it to the LDLs. However, this might explain why the number and magnitude of triggered earthquakes are relatively

<figure><img src="image_7.png" /><figcaption>**Fig. 4 | Schematic representation of earthquake-surface hazards interaction.**
Previous studies⁵,⁶,¹³,¹⁶ had shown that earthquakes can trigger landslides which can block rivers and form landslide-dammed lakes. We show that direct gravitational</figcaption></figure>

loading and pore pressure diffusion from landslide-dammed lakes can in turn increase stresses on surrounding faults and trigger earthquakes. The interaction between these hazards is thus bidirectional.
small since earthquake triggering also depends on other factors such as the availability of critically stressed faults. Since cases of reservoir-induced earthquakes are relatively common and have included M > 7 earthquakes⁴⁹, LDLs in seismically more active regions can potentially trigger more and larger earthquakes, assuming they increase ΔCFS given the fault locations and orientations. For example, the 2008 M7.9 Wenchuan earthquake triggered more than 100,000 landslides and formed a few hundred LDLs⁵²⁻⁵⁴. Some landslides were >30 times larger than our landslides³³ and the largest LDL lasted for -1 month and reached a maximum depth of ~82 m before breaching⁵⁴. Therefore, we suggest that such LDLs have greater potential to trigger earthquakes and the earthquake-surface interaction could even form a feedback loop (Fig. 4). Furthermore, while reservoir planning aims to avoid seismically active regions, ~20% of LDLs are estimated to have formed by earthquake-triggered mass wasting³⁰. Consequently, earthquake triggering by LDLs might be a ubiquitous phenomenon in seismically active regions though at least to some extent, this effect may be masked by earthquake-earthquake triggering in such regions⁵⁵.

While we only demonstrated earthquake triggering by LDLs, other surface hazards may similarly alter the stress state of surrounding faults and trigger earthquakes. For example, glacial lakes whose volumes are increasing due to climate change⁵⁶ can alter the stress state of surrounding faults through both direct gravitational loading and pore pressure diffusion. Therefore, our results suggest that such bidirectional interaction between surface hazards and earthquakes should be properly accounted for in future risk assessment of geological hazards and hazard management in mountainous regions.

## Methods

### Earthquake catalog building

We adopt an AI-based workflow to develop a new earthquake catalog from 2014 to 2023 based on the 13 permanent seismic stations⁵⁷ (Fig. 1). We first use the machine-learning-based detector, PhaseNet⁵⁸, to identify P- and S-wave arrival times from the continuous waveforms. Subsequently, these P- and S-wave arrival times are associated into individual earthquakes using a high-throughput seismic phase associate, PyOcto⁵⁹, and a local 1-D velocity model⁶⁰. A threshold of at least two P picks, one S pick, and a total of five P and S picks is adopted during the PyOcto association, resulting in 6692 earthquakes within a region of 2.5° × 3° (latitude: 30° to 32.5°; longitude: 97° to 100°). We then manually inspect the waveforms of 859 earthquakes around our target area (latitude: 30.7° to 31.7°; longitude: 98.1° to 99.2°) to exclude 112 non-tectonic or unreasonably associated events. We then use the absolute location method Hypoinverse⁶¹ to further refine the hypo-centers of the remaining 747 events and calculate their location uncertainties.

We also estimate the local magnitudes for all retained earthquakes based on their S-wave amplitudes and a recently improved national standard magnitude scale that is specific to the Sichuan region, China⁶²,⁶³. The maximum amplitudes of horizontal component wave-forms are measured after deconvolving the instrument response from the raw waveforms and then convolving the obtained signal with the theoretical Wood-Anderson seismometer response. The measured waveform window starts 0.5 s before the P wave arrival and is twice the predicted S-P travel time in length. Note that we only estimate magnitudes for 669 events with at least two stations having signal-noise ratios of S waves larger than 3. The estimated magnitudes range from M_L 0.1 to 4.2.

To further improve the magnitude of completeness, we use a template matching technique, Graphics Processing Unit-based match and locate technique (GPU-M&L)⁶⁴, to detect and locate more earthquakes based on the 669 retained earthquakes with estimated magnitude. P- and S-waves of these events are simultaneously used in GPU-M&L. The length of template window is 6 s starting 1 s before the P- and S-wave arrivals picked by PhaseNet. We filter the continuous and

template waveforms from 2 to 8 Hz. The 1-D velocity model used in the above absolute earthquake location was adopted in GPU-M&L as well. Note that in GPU-M&L, we relocate the epicenters of detections through a grid search method with searching space 1.2 × 1.2 km in the horizontal direction and searching interval of 0.06 km if template and continuous waveforms were recorded by at least three seismic stations. Otherwise, the locations of detections will be defined as that of the templates. By manually inspecting these detections, the detection threshold is defined as correlation coefficient of 0.3 and 12 times median absolute deviation (MAD) if the common component number between template and continuous waveforms is ≥ 9. Otherwise, the detection threshold is defined as correlation coefficient of 0.5 and 15 times MAD. In total, GPU-M&L resulted in 4294 detections including the 669 template events. The magnitudes of the newly detected earthquakes are estimated based on the local magnitudes of the templates and the amplitude ratio between these detections and their corresponding templates⁶⁵. Note that due to the poor local seismic station coverage with available data before May 2018, 22 earthquakes that occurred from March 2014 to May 2018 in the routine catalog, which utilized data from the entire China, were not detected. Thus, these missing earthquakes are manually added to our earthquake catalog.

### Declustering and seismicity rate change analysis

Using the approach from Reasenberg⁶¹, we identified 5 clusters of earthquakes containing a total of 13 events (5 main earthquakes and 8 foreshocks/aftershocks) out of 62 events which occurred over ~5 years, using the following input parameters⁶¹: look-ahead time from 0.5 to 3 days, confidence probability of 0.9, effective min magnitude cutoff of 1.1, the increase in lower cutoff magnitude during clusters of 0, and the number of crack radii surrounding each earthquake of 5.

We then evaluate the statistical significance of the observed seismicity rate change based on statistic P (ref. 42) and statistic Z (ref. 43). The statistic Z is defined as follows:

$$Z = \frac{N\Delta t_M - M\Delta t_N}{\sqrt{N\Delta t_M^2 + M\Delta t_N^2}} \quad (1)$$

Where N and M represent the number of earthquakes after and before an event in time windows Δt_N and Δt_M, respectively. When N and M are sufficiently large and in the null hypothesis of no seismicity rate change, Z follows a Gaussian distribution with zero mean and unit standard deviation. Statistic Z was first proposed by Habermann⁴³ to evaluate the statistical significance of seismicity rate changes. After declustering, there were 9 earthquakes in the week (10 to 16 November) when the second LDL approached its peak water level. There were only 8 earthquakes in the previous 6 months (~27 weeks) which translates to a seismicity rate of 0.3 events/week. Based on Eq. 1, we obtain a Z value of 2.90 taking N, M, Δt_N, and Δt_M as 9, 8, 1, and 27, respectively. A change is considered statistically significant if |Z| > 2 (ref. 43), and our observed seismicity rate increase after declustering⁴¹ is statistically significant at a 99.62% level based on statistic Z.

In the null hypothesis that the seismicity rate in the week after (Poisson process with mean rate λ_N) and the week before (Poisson process with mean rate λ_M) an event is the same, the statistic P reads:

$$P(\lambda_N > 1) = 1 - \frac{1}{N!M!} \int_0^\infty e^{-x} x^M \Gamma(N+1, x) dx \quad (2)$$

Where Γ(n,x) = ∫_0^x e^x t^n dt is the incomplete Gamma function and follows a uniform distribution between 0 and 1 (ref. 42). 1-P gives us the probability that an observed value could be obtained by chance if the null hypothesis of no rate change is true⁴². While there were no earthquakes in the -1 month before the second LDL, we still take the
seismicity rate for the week before the second LDL approached its peak water level as 0.3 events/week based on the average seismicity rate in the previous 6 months to be conservative. Based on Eq. 2, we obtain a *P* value of 0.9982 taking *N* and *M* as 9 and 0.3, respectively. This implies that the probability that our observed seismicity rate increase occurred by random chance is 0.18%. Therefore, our observed seismicity rate increase after declustering⁴¹ is statistically significant at a >99.82% level based on statistic *P*.

We further evaluate the statistical significance of the observed seismicity rate change based on both an improved statistic based on Poisson probability⁴⁴ and an empirically derived statistic⁴⁵ which were developed to determine statistically significant changes in earthquake rate under small background rate. We assumed that earthquakes occur independently at a constant rate, following a Poissonian distribution⁴⁴. To determine the statistical significance of seismicity rate increase, we compare the number of expected earthquakes in a 1-week window to the number of earthquakes in the 1-week window after the trigger and calculate the Poisson probability of obtaining the number of earthquakes (*ν*) in the 1-week window after the trigger given the expected number of events (*μ*) in the 1-week window.

$$P_{\mu}(\nu) = e^{-\mu} \frac{\mu^{\nu}}{\nu!} \qquad (3)$$

For this case of low background seismicity rate, similar to the idea of ref. 44, we conservatively regard the largest number of earthquakes (3) in a 1-week window over our observation period, except for this earthquake sequence, as the expected number of events (*μ*) in a 1-week window. Based on Eq. 3, we obtain a *P* value of 0.0027 taking *μ* and *ν* as 3 and 9, respectively. This implies that our observed seismicity rate increase after declustering⁴¹ is statistically significant at a >99.73% level.

We also use the empirical statistical method of ref. 45 to determine the significance in the seismicity rate change. First, we count the number of earthquakes in each 1-week time window incorporating a 1-day sliding window. There are >1700 windows from May 2018 to January 2023. We then assign the number of earthquakes *N*<sub>count</sub> within each window to a timestamp at the start of the window. We next build a histogram of *N*<sub>count</sub> values (Fig. S6). Finally, using percentile measurements, the >99% statistically significant level is 3 earthquakes in a 1-week window. We observed 9 earthquakes (after declustering) as the second LDL approached the peak water level, hence this seismicity rate increase is statistically significant at a >99% level.

We also tested the impact of different input parameters on the declustering based on the approach from Reasenberg⁴¹. We identified 4 clusters of earthquakes containing a total of 16 events (4 main earthquakes and 12 foreshocks/aftershocks) out of 62 events which occurred over ~5 years, using other suggested input parameters⁴¹: look-ahead time from 1 to 10 days, confidence probability of 0.95, effective min magnitude cutoff of 1.5, the increase in lower cutoff magnitude during clusters of 0.5, and the number of crack radii surrounding each earthquake of 10. While the declustered earthquake catalogs with different input parameters are different, we further confirm that this seismicity rate increase is statistically significant at a > -90% level based on the four above statistical tests (Table S1).

## Coulomb stress change calculation

Coulomb stress change (ΔCFS) on the surrounding fault systems due to the LDL is a result of the combined effect of direct gravitational loading and pore pressure diffusion. We model ΔCFS due to both direct loading and pore pressure diffusion using the GeoTaos software⁶⁶. Since the direct gravitational loading effect can be regarded as the result of point forces acting vertically on the surface of a homogenous elastic half-space, the loading from the LDL can be calculated by the convolution of the Green's function and the distributed

surface forces. The distributed surface forces can be estimated based on the spatiotemporal evolution of the LDL's water level. To quantify the spatiotemporal evolution of the LDL's water level, the LDL is mapped as a series of square cells with a size of 200 by 200 m. The elevation of each square cell is estimated from the GDEM V3 30 m. Based on the elevation of each cell, the water depth of each cell can then be estimated. Therefore, the distributed surface forces and thus the direct gravitational loading effect of each cell can be calculated using the water depth of each cell. Note that we ignore the ΔCFS resulting from gravitational loading of landslide sediments, since these landslide sediments are not a new source of mass but only the very local redistributions (~1–2 km) of sediments along the landslide sliding path.

## Uncertainty analyses

Before more seismic stations became operational in May 2018, no earthquakes were detected within 10 km of the LDLs due to limited station coverage and/or low seismicity rate. Nevertheless, in the ~250 weeks of observation period since May 2018 for which we have an improved earthquake catalog, 61 earthquakes occurred within 10 km of the LDLs during the week when the second LDL's water level almost peaked while there were only 183 earthquakes over the ~247 weeks outside the time period of the two LDLs i.e., a background seismicity rate of only ~0.7/week over ~5 years (Table S1). Therefore, our catalog shows that the seismicity rate during the LDL is significantly higher than the background seismicity rate of the region. We chose to focus on earthquakes that occurred within 10 km of LDLs since the region with stress increase of >0.01 MPa is always located within 10 km of the LDLs (Fig. 2a). Nevertheless, the earthquake activity within larger and smaller regions shows a similar accelerated trend during the landslide hazard cascades (Fig. S16).

As the second LDL approached peak water level, ~50% and ~95% of the hypocenters in the earthquake sequence were at depths of <1.2 and <4.5 km, respectively (Fig. S8). Similarly, for the earthquake catalog before template matching, the median depth of these triggered earthquakes is ~1.2 km (Fig. S8). However, the mean horizontal and vertical location uncertainties of our cataloged earthquakes are ~4 and ~8 km, respectively. To further constrain the likely depths of these events, we look at the larger earthquakes which are more accurately located since they are recorded by more stations. We find that the median depth of all $M_L > 1.1$ earthquakes during the second LDL is ~1.5 km with average depth uncertainty of ~6 km. Our modeling shows that ΔCFS down to a depth of 7 km can reach levels that trigger earthquakes (Fig. S11d), and the ΔCFS is larger at shallower depths (Fig. S11). Furthermore, ~60% of reservoirs with documented earthquake depths have induced earthquakes up to ~8 km depth (Fig. S17). Therefore, we conclude that the earthquake sequence can be triggered by stress changes due to the LDL.

The stress change on surrounding faults due to the LDLs partly depends on the fault geometry. In contrast to the well-documented fault strike, the dip of the fault is less well constrained. We assume a fault dip of 75°NE based on the focal mechanisms of four $M > 3.5$ earthquakes nearby and the consistency of the first motions of the two largest earthquakes (ML 2.6 and 1.9) as the second LDL approached peak water level with these focal mechanisms (Fig. S7), though we cannot confirm that the focal mechanisms of the smaller earthquakes necessarily share similar geometries. Nevertheless, our modeling shows that the fault dip mainly affects the spatial distribution and amplitude of ΔCFS but not its sign, and the seismicity generally locates within regions of positive ΔCFS for a range of different assumed fault dips (Fig. S11e–h). Hence, the dip uncertainty has minimal effect on our conclusion that landslide hazard cascades can trigger earthquakes.

The stress change on surrounding faults caused by LDL is also related to the fault friction coefficient and the fluid diffusion
coefficient. We assumed a friction coefficient of 0.5 for our ΔCFS modeling. However, the chosen friction coefficient has minimal effect on both the spatial distribution and amplitude of ΔCFS (Fig. S18). In addition, due to differences in rock properties and depth-dependent confining pressure, crustal hydraulic diffusivity has been found to range from 0.01 to 5 m²/s (refs. 67,68). Both the geometry and closure of faults and fractures also cause significant spatial variation in hydraulic diffusivity⁶⁹,⁷⁰. These factors make it challenging to obtain a precisely constrained hydraulic diffusivity for the calculation of ΔCFS. Nevertheless, higher fluid diffusivity causes the pore water pressure to increase faster, which increases the amplitude of ΔCFS over a wider area. Our chosen value of 0.3 m²/s is on the lower end of observed values and is a typical value observed in fluid-driven seismic swarms in different regions⁷¹⁻⁷³. While hydraulic diffusivity strongly influences the time at which ΔCFS exceeds a certain threshold, we further confirm that ΔCFS exceeds 0.01 MPa for hydraulic diffusivities ranging from 0.02 to 2 m²/s. On the day of peak seismicity rate, the relative contribution to ΔCFS of pore pressure increases from -0.1 to 0.9 as the assumed hydraulic diffusivity increases from 0.02 to 2 m²/s (Fig. S14). We conclude that our main conclusion that LDLs can trigger earthquakes is robust to changes in the discussed parameter choices.

## Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

The reservoir data which induced earthquakes are available at the Human-Induced Earthquake Database (HiQuake) http://inducedearthquakes.org/. The database for landslide-dammed lakes is available from https://doi.org/10.1016/j.earscirev.2020.103116. The earthquake catalog produced in this study is available at https://doi.org/10.5281/zenodo.10807902. The GDEM data are provided by ASTER https://asterweb.jpl.nasa.gov.

## Code availability

The software PhaseNet⁵⁸ and PyOcto⁵⁹ used in this study are publicly available at https://doi.org/10.1093/gji/ggy423 and https://doi.org/10.26443/seismica.v3i1.1130, respectively. The software Hypoinverse⁶¹ and GPU-M&L⁶⁴ are publicly available at https://www.usgs.gov/software/hypoinverse-earthquake-location and https://doi.org/10.1785/0220190241, respectively. The software Zmap⁴⁰ and GeoTaos⁶⁶ are publicly available at http://www.seismo.ethz.ch/en/research-and-teaching/products-software/software/ZMAP and http://bemlar.ism.ac.jp/lxl/Taos/Download.htm, respectively.

## References

1. Allen, S. K., Rastner, P., Arora, M., Huggel, C. & Stoffel, M. Lake outburst and debris flow disaster at Kedarnath, June 2013: Hydro-meteorological triggering and topographic predisposition. Landslides 13, 1479-1491 (2016).

2. Hovius, N., Stark, C. P., Tutton, M. A. & Abbott, L. D. Landslide-driven drainage network evolution in a pre-steady-state mountain belt: finisterre mountains, Papua New Guinea. Geology 26, 1071-1074 (1998).

3. Korup, O. et al. Giant landslides, topography, and erosion. Earth Planet. Sci. Lett. 261, 578-589 (2007).

4. Segoni, S., Piciullo, L. & Gariano, S. L. A review of the recent literature on rainfall thresholds for landslide occurrence. Landslides 15, 1483-1501 (2018).

5. Keefer, D. K. et al. Real-time landslide warning during heavy rainfall. Science 238, 921-925 (1987).

6. Moreiras, S., Lisboa, M. S. & Mastrantonio, L. The role of snow melting upon landslides in the central Argentinean Andes. Earth Surf. Process. Landf. 37, 1106-1119 (2012).

7. Kawagoe, S., Kazama, S. & Sarukkalige, P. R. Assessment of snowmelt triggered landslide hazard and risk in Japan. Cold Reg. Sci. Technol. 58, 120-129 (2009).

8. Bozzano, F., Cipriani, I., Mazzanti, P. & Prestininzi, A. Displacement patterns of a landslide affected by human activities: insights from ground-based InSAR monitoring. Nat. hazards 59, 1377-1396 (2011).

9. Alexander, D. On the causes of landslides: human activities, perception, and natural processes. Environ. Geol. Water Sci. 20, 165-179 (1992).

10. Keefer, D. K. Landslides caused by earthquakes. Geol. Soc. Am. Bull. 95, 406-421 (1984).

11. Meunier, P., Hovius, N. & Haines, J. A. Topographic site effects and the location of earthquake induced landslides. Earth Planet. Sci. Lett. 275, 221-232 (2008).

12. Wasowski, J., Keefer, D. K. & Lee, C. T. Toward the next generation of research on earthquake-induced landslides: current issues and future challenges. Eng. Geol. 122, 1-8 (2011).

13. Kargel, J. S. et al. Geomorphic and geologic controls of geohazards induced by Nepal's 2015 Gorkha earthquake. Science 351, aac8353 (2016).

14. Ren, J. et al. Surface rupture of the 1933 M 7.5 Diexi earthquake in eastern Tibet: implications for seismogenic tectonics. Geophys. J. Int. 212, 1627-1644 (2018).

15. Fan, X. et al. Earthquake-induced chains of geologic hazards: patterns, mechanisms, and impacts. Rev. Geophys. 57, 421-503 (2019).

16. Marc, O., Hovius, N., Meunier, P., Uchida, T. & Hayashi, S. I. Transient changes of landslide rates after earthquakes. Geology 43, 883-886 (2015).

17. Korup, O., Densmore, A. L. & Schlunegger, F. The role of landslides in mountain range evolution. Geomorphology 120, 77-90 (2010).

18. Amos, C. B. et al. Uplift and seismicity driven by groundwater depletion in central California. Nature 509, 483-486 (2014).

19. Borsa, A. A., Agnew, D. C. & Cayan, D. R. Ongoing drought-induced uplift in the western United States. Science 345, 1587-1590 (2014).

20. Cochran, E. S., Vidale, J. E. & Tanaka, S. Earth tides can trigger shallow thrust fault earthquakes. Science 306, 1164-1166 (2004).

21. Bell, M. L. & Nur, A. Strength changes due to reservoir-induced pore pressure and stresses and application to Lake Oroville. J. Geophys. Res. Solid Earth 83, 4469-4483 (1978).

22. Schultz, R., Atkinson, G., Eaton, D. W., Gu, Y. J. & Kao, H. Hydraulic fracturing volume is associated with induced earthquake productivity in the Duvernay play. Science 359, 304-308 (2018).

23. Healy, J. H., Rubey, W. W., Griggs, D. T. & Raleigh, C. B. The denver earthquakes. Science 161, 1301-1310 (1968).

24. Gough, D. I. & Gough, W. I. Stress and deflection in the lithosphere near Lake Kariba—I. Geophys. J. Int. 21, 65-78 (1970).

25. Kebeasy, R. M. & Gharib, A. A. Active fault and water loading are important factors in triggering earthquake activity around Aswan Lake. J. Geodyn. 14, 73-85 (1991).

26. Simpson, D. W., Stachnik, J. C. & Negmatoullaev, S. K. Rate of change in lake level and its impact on reservoir triggered seismicity. Bull. Seismol. Soc. Am. 108, 2943-2954 (2018).

27. Hill, R. G., Weingarten, M., Rockwell, T. K. & Fialko, Y. Major southern San Andreas earthquakes modulated by lake-filling events. Nature 618, 761-766 (2023).

28. Gomberg, J. & Davis, S. Stress/strain changes and triggered seismicity at The Geysers. Calif. J. Geophys. Res. Solid Earth 101, 733-749 (1996).

29. Keranen, K. M., Savage, H. M., Abers, G. A. & Cochran, E. S. Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 5.7 earthquake sequence. Geology 41, 699-702 (2013).

30. Calais, E., Freed, A. M., Van Arsdale, R. & Stein, S. Triggering of New Madrid seismicity by late-Pleistocene erosion. Nature 466, 608-611 (2010).
31. Steer, P., Simoes, M., Cattin, R. & Shyu, J. B. H. Erosion influences the seismicity of active thrust faults. *Nat. Commun.* **5**, 5564 (2014).

32. Steer, P. et al. Earthquake statistics changed by typhoon-driven erosion. *Sci. Rep.* **10**, 10899 (2020).

33. Jeandet Ribes, L., Cubas, N., Bhat, H. S. & Steer, P. The impact of large erosional events and transient normal stress changes on the seismicity of faults. *Geophys. Res. Lett.* **47**, e2020GL087631 (2020).

34. Zhang, Z. et al. Source characteristics and dynamics of the October 2018 Baige landslide revealed by broadband seismograms. *Landslides* **16**, 777–785 (2019).

35. Zhang, L., Xiao, T., He, J. & Chen, C. Erosion-based analysis of breaching of Baige landslide dams on the Jinsha River, China, in 2018. *Landslides* **16**, 1965–1979 (2019).

36. Zhong, Q., Chen, S., Wang, L. & Shan, Y. Back analysis of breaching process of Baige landslide dam. *Landslides* **17**, 1681–1692 (2020).

37. Gao, Y., Zhao, S., Deng, J., Yu, Z. & Rahman, M. Flood assessment and early warning of the reoccurrence of river blockage at the Baige landslide. *J. Geogr. Sci.* **31**, 1694–1712 (2021).

38. Chen, F. et al. Kinematic process and mechanism of the two slope failures at Baige Village in the upper reaches of the Jinsha River, China. *Bull. Eng. Geol. Environ.* **80**, 3475–3493 (2021).

39. Guo, X., Jiang, C., Han, L., Yin, H. & Zhao Z. Focal mechanism data set in Chinese mainland and its adjacent area (2009–2021). *Natl. Earthquake Data Center*. https://data.earthquake.cn. https://doi.org/10.12080/nedc.11 ds.2022.0004 (2022).

40. Wiemer, S. A software package to analyze seismicity: ZMAP. *Seismol. Res. Lett.* **72**, 373–382 (2001).

41. Reasenberg, P. Second-order moment of central California seismicity, 1969-82. *J. Geophys. Res.* **90**, 5479–5495 (1985).

42. Marsan, D. & Nalbant, S. S. Methods for measuring seismicity rate changes: a review and a study of how the Mw 7.3 landers earthquake affected the aftershock sequence of the M w 6.1 Joshua Tree earthquake. *Pure Appl. Geophys.* **162**, 1151–1185 (2005).

43. Habermann, R. Precursory seismicity patterns: stalking the mature seismic gap. Earthquake prediction. *Int. Rev.* **4**, 29–42 (1981).

44. Alfaro-Diaz, R., Velasco, A. A., Pankow, K. L. & Kilb, D. Optimally oriented remote triggering in the Coso geothermal region. *J. Geophys. Res. Solid Earth* **125**, e2019JB019131 (2020).

45. Pankow, K. L. & Kilb, D. Going beyond rate changes as the sole indicator for dynamic triggering of earthquakes. *Sci. Rep.* **10**, 4120 (2020).

46. Wu, F., Jiang, L. & Zhang, G. Discussion on quaternary activity characteristics of northern section of Jinshajiang fault zone along Sichuan-Tibet Railway. *High. Speed Railw. Technol.* (Chin.) **10**, 23–28 (2019).

47. Reasenberg, P. A. & Simpson, R. W. Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake. *Science* **255**, 1687–1690 (1992).

48. Harris, R. A. Introduction to special section: Stress triggers, stress shadows, and implications for seismic hazard. *J. Geophys. Res. Solid Earth* **103**, 24347–24358 (1998).

49. Wilson, M., Foulger, G., Gluyas, J., Davies, R. & Julian, B. HiQuake: The human-induced earthquake database. *Seismol. Res. Lett.* **88**, 1560–1565 (2017).

50. Fan, X. et al. The formation and impact of landslide dams–state of the art. *Earth Sci. Rev.* **203**, 103116 (2020).

51. King, G. C., Stein, R. S. & Lin, J. Static stress changes and the triggering of earthquakes. *Bull. Seismol. Soc. Am.* **84**, 935–953 (1994).

52. Xu, C., Xu, X., Yao, X. & Dai, F. Three (nearly) complete inventories of landslides triggered by the May 12, 2008 Wenchuan Mw 7.9 earthquake of China and their spatial distribution statistical analysis. *Landslides* **11**, 441–461 (2014).

53. Huang, R. & Li, W. Analysis of the geo-hazards triggered by the 12 May 2008 Wenchuan Earthquake, China. *Bull. Eng. Geol. Environ.* **68**, 363–371 (2009).

54. Cui, P., Zhu, Y., Han, Y., Chen, X. & Zhuang, J. The 12 May Wenchuan earthquake-induced landslide lakes: distribution and preliminary risk evaluation. *Landslides* **6**, 209–223 (2009).

55. Hill, D. P., Pollitz, F. & Newhall, C. Earthquake-volcano interactions. *Phys. Today* **55**, 41–47 (2002).

56. Shugar, D. H. et al. Rapid worldwide growth of glacial lakes since 1990. *Nat. Clim. Change* **10**, 939–945 (2020).

57. Liu, M. et al. Complexity of initiation and evolution of the 2013 Yunlong earthquake swarm. *Earth Planet. Sci. Lett.* **612**, 118168 (2023).

58. Zhu, W. & Beroza, G. C. PhaseNet: a deep-neural-network-based seismic arrival-time picking method. *Geophys. J. Int.* **216**, 261–273 (2019).

59. Münchmeyer, J. PyOcto: A high-throughput seismic phase asso- ciator. *Seismica 3*. https://doi.org/10.26443/seismica.v3i1. 1130 (2024).

60. Fang, L. et al. Aftershock observation and analysis of the 2013 M$_{S}$ 7.0 Lushan earthquake. *Seismol. Res. Lett.* **86**, 1135–1142 (2015).

61. Klein, F. W. User's guide to HYPOINVERSE-2000, a Fortran program to solve for earthquake locations and magnitudes (No. 2002-171). *US Geological Survey* (2002).

62. Richter, C. F. An instrumental earthquake magnitude scale. *Bull. Seismol. Soc. Am.* **25**, 1–32 (1935).

63. Tang, L. et al. Comparison of the new national standard magnitude scales with the traditional magnitude scales. *Acta Seismol. Sin.* (Chin.) **40**, 121–131 (2018).

64. Liu, M., Li, H., Zhang, M. & Wang, T. Graphics processing unit-based match and locate (GPU-M&L): An improved match and locate method and its application. *Seismol. Res. Lett.* **91**, 1019–1029 (2020).

65. Shelly, D. R., Ellsworth, W. L. & Hill, D. P. Fluid-faulting evolution in high definition: connecting fault structure and frequency-magni- tude variations during the 2014 Long Valley Caldera, California, earthquake swarm. *J. Geophys. Res. Solid Earth* **121**, 1776–1795 (2016).

66. Lei, X. Possible roles of the Zipingpu Reservoir in triggering the 2008 Wenchuan earthquake. *J. Asian Earth Sci.* **40**, 844–854 (2011).

67. Husen, S., Bachmann, C. & Giardini, D. Locally triggered seismicity in the central Swiss Alps following the large rainfall event of August 2005. *Geophys. J. Int.* **171**, 1126–1134 (2007).

68. Montgomery-Brown, E., Shelly, D. R. & Hsieh, P. A. Snowmelt-triggered earthquake swarms at the margin of Long Valley Caldera, California. *Geophys. Res. Lett.* **46**, 3698–3705 (2019).

69. Kuempel, H. Poroelasticity: parameters reviewed. *Geophys. J. Int.* **105**, 783–799 (1991).

70. Scholz, C. H. *The Mechanics of Earthquakes and Faulting* (Cam- bridge University Press, 2002).

71. Kitagawa, Y., Fujimori, K. & Koizumi, K. Temporal change in per- meability of the Nojima fault zone by repeated water injection experiments. *Tectonophysics* **443**, 183–192 (2007).

72. Parotidis, M., Shapiro, S. A. & Rothert, E. Evidence for triggering of the Vogtland swarms 2000 by pore pressure diffusion. *J. Geophys. Res.* **110**, B05SS10 (2005).

73. Hainzl, S. & Ogata, Y. Detecting fluid signals in seismicity data through statistical earthquake modeling. *J. Geophys. Res.* **110**, B05SS07 (2005).

Acknowledgements

We thank Jinping Zi, Peifeng Wang, and Yiyuan Zhong for valuable dis- cussions. This work was supported by the Croucher Tak Wah Mak Innovation Award (Y.J.T.), CUHK Direct Grant for Research (Grant 4053594; Y.J.T.), CUHK Research Fellowship Scheme (Grant 4200720; Y.J.T.), CUHK Postdoctoral Fellowship Scheme (Grant 3135043; Y.J.T.), the National Key Research and Development program of China (Grant No.2022YFF0800604; S.H.), and the Major Program of the National Natural Science Foundation of China (Grant No.42090051; S.H.).
**Author contributions**

Conceptualization: Z.Z. and Y.J.T.; Investigation: Z.Z., M.L., Y.J.T., F.W., S.H., M.C., and J.S.; Methodology: Z.Z., M.L., and Y.J.T.; Project administration: Y.J.T. and S.H.; Supervision: Y.J.T. and S.H.; Writing-original draft: Z.Z. and Y.J.T.; Writing-review and editing: Z.Z., M.L., Y.J.T., F.W., S.H., M.C., and J.S.

**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/s41467-024-47130-w.

**Correspondence** and requests for materials should be addressed to Zhen Zhang or Yen Joe Tan.

**Peer review information** *Nature Communications* thanks Ya-Ju Hsu, Philippe Steer and the other, anonymous, reviewer for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions 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.

**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) 2024