Stress triggering of the 2022 Lushan–Maerkang earthquake sequence by historical events and its implication for fault stress evolution in eastern Tibet

It remains unknown how stress triggering causes earthquakes in the eastern Tibetan Plateau following the Wenchuan Earthquake in 2008. The MS 6.1 Lushan earthquake on 1 June 2022 in a seismic gap between the 2008 MW 7.9 Wenchuan earthquake and the 2013 MW 6.6 Lushan earthquake provided an opportunity to detect stress evolution and seismic activity on the fault in this study. We calculated the Coulomb stress change of the June 2022 Lushan–Maerkang earthquake sequence using a Burgers viscoelastic model and, herein, discuss how the sequence have been triggered by historical earthquakes since 1900 in the eastern Tibetan Plateau. Our results suggest the following trends: 1) the 1955 M 7.6 Kangding and 2008 MW 7.9 Wenchuan earthquakes contributed most significant loading effects on the 2022 MS 6.1 Lushan earthquake; however, the 2013 MW 6.6 Lushan earthquake had an unloading effect on the 2022 Lushan earthquake. 2) The 2021 MW 7.3 Maduo earthquake contributed a loading effect on the 2022 Lushan earthquake, and the 2022 Lushan earthquake may have triggered the subsequent Maerkang earthquake swarm on 10 June 2022. 3) Viscoelastic relaxation of the lower crust and upper mantle contributed significantly to fault stress level, while a fault in a late-phase earthquake cycle may have ruptured via slight stress perturbation near a cross-border conversion between positive and negative stress from a far-field earthquake. We also provide a seismic potential assessment along the faults in the eastern Tibet. Notably, the MW 6.8 Luding earthquake that ruptured the southern segment of the Xianshuihe Fault on 5 September 2022 supports the conclusions of this study.

provided an epicenter of 30.35°N, 103.08°E, with a moment magnitude of M W 5.8 and a source depth of 21.6 km.
Previous studies have indicated that the 2013 Lushan earthquake was located in an area where the Coulomb stress change increased due to the 2008 Wenchuan earthquake and was triggered by that earthquake (Toda et al., 2008;Lei et al., 2013;Jia K et al., 2014;Liu et al., 2014;Wang et al., 2014;Xie et al., 2014;Zhu and Miao, 2015;Zhu, 2016;Lin et al., 2019). If a triggering correlation exists, it is unclear why the 2022 Lushan earthquake occurred 9 years later. The potential seismogenic trend of faults around the Longmenshan fault system is also in question. The Lushan earthquake and Maerkang earthquake swarm that occurred in June 2022 provided a significant opportunity to study the correlation of earthquake triggering. Therefore, we studied the stress triggering of the 2022 Lushan earthquake and the three subsequent Maerkang earthquakes from the events shown in Table 2. Herein, we discuss the subsequent seismic trend of surrounding faults and the triggering pattern of the Lushan-Maerkang earthquake sequence via historical earthquakes.
2 Data and method

Stratified viscoelastic model
We used a stratified viscoelastic model to simulate the Songpan-Ganzi Terrane, where earthquakes mainly occur along the Longmenshan Fault. The wave velocity and density structure of the crust and upper mantle were taken from Jia S et al. (2014) and Xu et al. (2010). We assumed that the thickness of the elastic layer of the crust was 30 km, and that the viscoelastic material was below that depth. The viscosity coefficient mainly refers to the rheological structural parameters of the lithosphere in eastern Tibet obtained by Wang et al. (2021) based on deformation simulation after the 2008 M W 7.9 Wenchuan earthquake, as shown in Table 3.
The Maxwell body and Kelvin body in rheological structures have defects in fitting short-term and long-term deformation, respectively. Therefore, we chose the Burgers body, which is suitable for simulating transient elastic response, short-term exponential decay response, and long-term linear increase response (Pollitz and Wicks, 2001;Pollitz and Sacks, 2002;Shao et al., 2007). The constitutive relation (Malkin and Isayev, 2022) is as follows: where σ is the stress, ε is the strain, and k 1 and k 2 are the stiffness coefficients of the Kelvin model and Maxwell model, respectively. η 1 and η 2 are the Kelvin model viscosity (short-term viscosity) and Maxwell model viscosity (long-term viscosity), respectively. k 2 , the effective stiffness coefficient, is obtained from the S-wave velocity and density of the medium, and the incomplete relaxation stiffness coefficient is k 1 αk 2 /(1 − α), in which α is set to 0.67, according to Ryder et al. (2011). Unless otherwise specified, the Coulomb stress change mentioned in this paper is the sum of coseismic and postseismic viscoelastic relaxation.

Method
The formula for calculating the Coulomb failure stress change on the fault plane is as follows (King et al., 1994;Harris, 1998):

Frontiers in Earth Science frontiersin.org 04
where Δτ is the shear stress change on the fault plane (taking the fault sliding direction as positive), Δσ n is the normal stress change on the fault plane (the tension is positive), and μ′ is the effective friction coefficient of the fault. The effective friction coefficient is 0.4 (King et al., 1994). We used the PSGRN/PSCMP program for this work (Wang et al., 2006).

Coseismic slip model
The stress field induced by an earthquake can be calculated using fault dislocation theory and the source fault slip model. The Coulomb stress change on the receiving fault plane can be obtained using Formula (1). In our study, we used two methods to obtain the slip model of the source fault: one method used the existing accurate slip distribution model of an earthquake, such as that of the 2021 M W 7.3 Maduo earthquake released by the USGS; in the second method, the empirical formula of the relationship between the fault length and magnitude was obtained using statistics (Wells and Coppersmith, 1994), the scalar seismic distance formula, and the moment magnitude definition formula (Knopoff, 1958;Hanks and Kanamori, 1979): The thrust earthquake equations are as follows (coefficients are the average values): The strike-slip earthquake equations are as follows (coefficients are the average values): where RLD is the length of the fault along the strike direction, RW is the width of the fault along the dip direction, M W is the magnitude of the moment, M 0 is the scalar seismic moment, μ is the shear modulus of the medium, and D is the average dislocation of the fault. According to the existing seismic fault parameters, the length, width, and average dislocation of the fault can be calculated using Formulas 3ormulas -Formulas 8. See Table 2 for the specific seismic fault rupture parameters.
Afterslip on historical earthquake faults can affect some simulation results (Nur and Mavko, 1974). However, we did not consider the afterslip effect in the calculation; compared with the main earthquake, the contribution of afterslip is concentrated in the near field, and the stress change it causes is much smaller than the effect of the coseismic slip model (Shen et al., 2009;Wang et al., 2011).

Triggering effect on Coulomb stress change in the 2022 Lushan earthquake
We calculated the coseismic and postseismic viscoelastic Coulomb stress change effects of historical earthquakes on the 2022 Lushan earthquake using the seismic rupture parameters, as shown in Table 2. We show the Coulomb stress change at a depth of 12 km, which is also the result of waveform fitting, and from the USGS. The focal mechanism parameter of the receiving fault is the result from the USGS in Table 1.
The postseismic shear stress change (Δτ) at the hypocenter of the 2022 Lushan earthquake (30.395°N, 102.958°E) was positive (consistent with the direction of fault slip), which contributed to fault rupture (Figure 2A). The normal stress change (Δσ n ) at the hypocenter was in the negative shadow zone ( Figure 2B), indicating that the 2008 Wenchuan earthquake caused a significant extrusion change in the fault normal direction, which was not prone to fault instability. The Coulomb stress change (ΔCFS) at the hypocenter was in the positive shadow zone ( Figure 2C), so the 2008 Wenchuan earthquake had a very obvious loading effect on the 2022 Lushan earthquake. According to Table 2, the coseismic ΔCFS of the 2008 Wenchuan earthquake on the 2022 Lushan earthquake was 2.36 × 10 4 Pa, and the ΔCFS at the time of the 2022 Lushan earthquake occurrence was 3.45×10 4 Pa, indicating that the loading effect of the 2008 Wenchuan earthquake on the 2022 Lushan earthquake increased with ongoing loading of the viscoelastic lower crust and upper mantle.
The hypocenter of the 2022 Lushan earthquake was near the cross-border conversion between positive and negative Δτ, but had positive values ( Figure 2G). The hypocenter of the 2022 Lushan earthquake was near the Δσ n cross-border conversion between positive and negative, but had negative values ( Figure 2H). The hypocenter of the 2022 Lushan earthquake was in the negative shadow zone of ΔCFS ( Figure 2I   Frontiers in Earth Science frontiersin.org Under the combined effects of historical earthquakes, the Δτ at the hypocenter of the 2022 Lushan earthquake was positive, and the Δσ n was negative. Due to their combined effects, all of which were near the cross-border conversion between positive and negative stress change of the three, ΔCFS was positive ( Figures 2M-O). The overall ΔCFS value was high, and the accumulated seismic moment was at a high level. Therefore, under the influence of the 2021 Maduo earthquake, the 2022 Lushan earthquake occurred along the southwest section of the Longmenshan Fault. Figure 3C shows the evolution of ΔCFS of historical earthquakes on the 2022 Lushan earthquake alone. As shown by the solid blue  Frontiers in Earth Science frontiersin.org 08 line (12 km) in (A) and (C), the ΔCFS remained at a small positive value due to the great distance from the previous historical earthquakes. However, after the 2008 Wenchuan earthquake, the ΔCFS rapidly increased to exceed the threshold value of 0.1 bar (King et al., 1994;Harris, 1998); the coseismic effect of the 2013 Lushan earthquake then caused it to decrease to a negative value. Under the continuous viscoelastic relaxation of the lower crust and upper mantle, the ΔCFS went from negative to positive, resulting in the M S 6.1 earthquake that occurred on 1 June 2022. This pattern may also explain why the M S 6.1 earthquake gap did not rupture when the 2013 Lushan earthquake occurred. In summary, Table 2 and Figure 3 show that the loading effect of the 2008 Wenchuan earthquake on the 2022 Lushan earthquake was dominant in these historical earthquakes, while the unloading effect was dominant in the 2013 Lushan earthquake.

Triggering effect on Coulomb stress change in the Maerkang earthquakes
After the Lushan M S 6.1 earthquake on 1 June 2022, the M W 5.6 (MEK1), M W 5.9 (MEK2), and M W 4.9 (MEK3) earthquakes occurred within a few hours on 10 June 2022 in Maerkang County. The three successive earthquakes were very close in time and space. We studied the triggering effect on Coulomb stress change using the earthquake rupture parameters in Table 2.
We considered a stress change at a depth of 10 km, consistent with the hypocenter depth given by the USGS. The focal mechanism solutions of the three earthquakes in Maerkang were equivalent and similar, and these solutions can be regarded as a single event (MEK1) under the action of historical earthquakes.
The Δτ and ΔCFS at the MEK1 hypocenter was negative, which was not conducive to fault slip ( Figures 4A, C). The Δσ n was positive, indicating that the normal stress change of the Maerkang earthquakes caused by the Wenchuan earthquake changed significantly ( Figure 4B), which contributed to the instability of the fault.
The hypocenter of MEK1 was near the cross-border conversion between positive and negative Δτ, and its value was positive ( Figure 4D). The hypocenter of MEK1 was also near the crossborder conversion between positive and negative Δσ n , and its value was negative ( Figure 4E). The hypocenter of MEK1 was in the ΔCFS positive shadow area ( Figure 4F), so the 2013 Lushan earthquake may have triggered the 2022 Maerkang earthquake. Remarkably, the focal mechanism solutions of the 2021 Maduo earthquake and the 2017 Jiuzhaigou earthquake, similar to that of the Maerkang earthquake, were also located near the positive and negative cross-border conversion of shear stress change and normal stress change caused by the 2013 Lushan earthquake. This finding indicates that the 2013 Lushan earthquake might also have triggered the 2017 Jiuzhaigou earthquake and the 2021 Maduo earthquake.
The triggering effect of the 2022 Lushan earthquake was similar as shown in Figures 4J-L, the Δτ and Δσ n at the hypocenter of the 2022 Maerkang earthquake included one positive and one negative value, but the ΔCFS was entirely in the positive shadow area. Therefore, the 2021 Maduo earthquake also had a certain loading effect on the 2022 Maerkang earthquakes to that of the 2013 Lushan earthquake Figures 4G-I.
According to Figure 5, the MEK2 was located in the ΔCFS positive shadow area caused by the MEK1, which triggered the MEK2, and the MEK3 was located near the cross-border conversion between positive and negative ΔCFS caused by the MEK1 and The hypocenter of the Lushan-Maerkang earthquake sequence seems to have been near the cross-border conversion between positive and negative stress change. It appears that the stress change at the hypocenter experienced a transform from negative to positive when the hypocenter finally became a starting point of rupture and the earthquake occurred. This phenomenon has been confirmed by Freed (2005), whose work showed that the hypocenter of the 1999 Hector Mine earthquake was very close to the border between positive and negative Coulomb stress change caused by the 1992 Landers earthquake. In addition, Xie et al. (2022) found that the hypocenter of the 2011 M S 9.0 Tokyo earthquake was located near the border between positive and negative normal stress change. We assumed that the positive and negative changes in stress were more likely to cause the rock at the hypocenter to plastically fail and thereby cause an earthquake.

Effect of model viscosity parameters on the results
The viscosity of the lower crust and upper mantle in the Tibetan Plateau and its surrounding areas has been studied over the past two decades. Clark and Royden (2000) used the topographic gradient method to estimate the viscosity beneath the lower crust around the Tibetan Plateau and obtained approximate results of 10 16~1 0 21 Pa s. Shao et al. (2011) considered the viscoelastic relaxation effect of the Plateau remains controversial, as the viscosity at different time scales changes, but it is concentrated in the range of 10 18~1 0 20 Pa s. Therefore, testing of different viscosity coefficients is necessary. We used the latest results from Wang et al. (2021) and tested the influence of two groups of viscosities (one order of magnitude lower and one order of magnitude higher) on the Coulomb stress change calculation results, as shown in Figures 6A-C. The ΔCFS evolution at the hypocenter of the 2022 Lushan earthquake using the three groups of viscosities is shown in Figure 3A. Variation in viscosity can induce different Coulomb stress changes, but it has little influence on the distribution range and evolution trend of ΔCFS. Therefore, the viscosity conditions considered in our study were appropriate.

Effect of the effective friction coefficient
The effective friction coefficient in Formula (2), involving the pore fluid and medium of the fault plane, is an uncertain parameter in Coulomb stress change calculations (King et al., 1994;Harris, 1998;Tang et al., 2023), and it is important to discuss the influence of its sensitivity on the results. King et al. (1994) proposed that the effective friction coefficients are generally 0. 2-0.8. Shen (2003) showed that the effective friction coefficient of most faults can change only the relative magnitude of ΔCFS, but cannot affect the polarity of ΔCFS.
To ensure the reliability of our results, we calculated the influence of the four groups of effective friction coefficients, i.e., 0.2, 0.4, 0.6, and 0.8, on the Coulomb stress change. As the effective friction coefficient increased from 0.2 to 0.8, the distribution range of Coulomb stress change in the 2022 Lushan earthquake area showed little variation (Figure 7). The ΔCFS at the hypocenter of the 2022 Lushan earthquake was still near the border between positive and negative stress change, but the value was positive. This result shows that the ΔCFS changes with the effective friction coefficient. However, the change in the effective friction coefficient in our study does not alter our conclusion.

Effect of depth on the results
Previous studies focused on the triggering effect of Coulomb stress change between earthquakes with stress fields at the depth of the hypocenter (Freed and Lin, 2001;Cheng, 2018;Li et al., 2022;Stein, 1999). However, hypocenter depth is difficult to estimate. The USGS suggested two hypocenter depths of 12 km and 19.5 km for the 2022 M S 6.1 Lushan earthquake, and the GFZ and the CENC suggested hypocenter depths of 10 km and 17 km. respectively. For strike-slip faults, the change in the Coulomb stress field with depth is small (Wan et al., 2007), but the Longmenshan Fault is under oblique-strike motion. Thus, the contribution of the Coulomb stress change on fault planes with different depths is discussed in further detail here. As shown in Figure 8, the ΔCFS at the epicenter of the 2022 Lushan earthquake was calculated for depths of 12 km, 16 km, and 20 km. The ΔCFS distribution range of the Longmenshan area varies significantly with depth. The 2022 Lushan earthquake occurred near the cross-border conversion between positive and negative stress change, but the ΔCFS value changed from positive to negative as the depth increased. Combined with Figure 3A, the ΔCFS evolution curves of the 2022 Lushan earthquake epicenter at different depths were consistent and coincided before 2013. However, the 2013 Lushan earthquake caused great differences at different depths. On the fault of the 2022 Lushan earthquake, the historical earthquakes had a strong unloading effect at a depth of 20 km, but a slight unloading effect at depths of 16 km and 12 km. Therefore, ΔCFS may vary at different depths on the fault plane. Thus, the influence of different hypocenter depths requires further study.

Effect of the focal mechanism solution of the receiving fault on the results
We evaluated focal mechanism solutions from the USGS. As shown in Figure 3B, the fluctuations in the curves of different focal mechanism solutions were similar. However, the effects of the 2013 Lushan earthquake on the 2022 Lushan earthquake using the GFZ solution were quite different from those of the other four groups. The large uncertainties in the strike and dip angle in the GFZ inversion results may be the source of the discrepancy. Therefore, the collection of focal mechanism solutions considered in our study greatly influenced the calculated results. We suggest that, when a similar calculation of fault stress change is carried out to study the interaction mechanism between strong earthquakes, focal mechanism solutions as similar as possible to actual earthquakes should be considered.

Influence of the Lushan-Maerkang earthquake sequence and historical earthquakes on the faults around the Songpan-Ganzi Terrane
We calculated the Coulomb stress change along the faults around the Songpan-Ganzi Terrane after the Lushan-Maerkang

Frontiers in Earth Science
frontiersin.org earthquake sequence in 2022. The fault slip data were from previous publications (Deng, 2007;Xu et al., 2016, https://www.activefaultdatacenter.cn/map). Figure 9 shows the Coulomb stress change induced by historical earthquakes on the faults in the Songpan-Ganzi Terrane at a depth of 10 km. The eastern section of the East Kunlun fault zone and the central section of the Longriba fault zone are in a continuous zone of high Coulomb stress change, which is consistent with the results of Cheng (1983) and Wang and Xu (2017). The western section of the Songgang Fault lies along an area in which ΔCFS changes sign, which corresponds to the occurrence of the Maerkang earthquake swarm in 2022. Additionally, the southern section of the Xianshuihe Fault has a higher ΔCFS value, similar to the ΔCFS distribution along the Xianshuihe Fault in Liu e al. (2014). The M S 6.8 Luding earthquake that occurred on 5 September 2022 also confirmed our prediction.
Furthermore, the middle section of the Yuke Fault, the southern section of the Xianshuihe Fault, the eastern section of the Dari Fault, the Bayankala Fault, the Wudaoliang-Changshagongma Fault, and the Chuan-Zang railway are near areas in which the ΔCFS changes sign, therefore potential earthquake risk may be high.

Conclusion
We constructed a stratified viscoelastic model to simulate postseismic stress on the 2022 Lushan-Maerkang earthquake sequence from historical strong earthquake data in eastern Tibet. We further assessed the seismic potential of faults in eastern Tibet. Our conclusions are as follows.
(1) The 1955 M 7.6 Kangding and 2008 M W 7.9 Wenchuan earthquakes are the most significant loading events on the 2022 M S 6.1 Lushan earthquake. But, the 2013 M W 6.6 Lushan earthquake is the most significant unloading event. This is a potential explanation for the occurrence of the 2022 Lushan earthquake 9 years later. Effects of viscoelastic relaxation of the lower crust and upper mantle play a significant role in fault activity. (2) The M W 5.6, M W 5.9, and M W 4.9 Maerkang earthquakes, which occurred after the 2022 Lushan earthquake, could have been due to both the 2013 and 2022 Lushan earthquakes. The 2021 Maduo earthquake had a loading effect on the 2022 Lushan earthquake, and the 2022 Lushan earthquake might have triggered the Maerkang earthquakes. Furthermore, the M W 5.6 Maerkang earthquake may have triggered the M W 5.9 Maerkang earthquake, and the M W 5.6 and M W 5.9 Maerkang earthquakes may have jointly triggered the M W 4.9 Maerkang earthquake. (3) The tectonic stress of the Songpan-Ganzi Terrane has accumulated to a high level. The hypocenters of the 2022 Lushan-Maerkang earthquake sequence all occurred near the cross-border conversion between positive and negative stress change. Rupture of a fault in a late phase of the seismic cycle near the cross-border conversion between positive and negative stress change may have been induced by a slight stress disturbance from a far-field earthquake. Therefore, study of the the middle section of the Yuke Fault, the southern section of the Xianshuihe Fault (which has been confirmed by the 5 September 2022, M W 6.8 Luding earthquake), the eastern section of the Dari Fault, the Bayankala Fault, and the Wudaoliang-Changshagongma Fault should be emphasized in future research, as they all correspond to areas of cross-border conversion.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, and further inquiries can be directed to the corresponding author.

Author contributions
DT and WG conceived the study and wrote the manuscript. XC helped archive fault trace data. All authors discussed, commented on, and edited the manuscript.