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: σ + η 2 k 1 + η 1 k 1 + η 2 k 2 σ ˙ + η 1 η 2 k 1 k 2 σ ¨ = η 2 ε ˙ + η 1 η 2 k 1 ε ¨ (1) 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.

The formula for calculating the Coulomb failure stress change on the fault plane is as follows (King et al., 1994; Harris, 1998): ∆ C F S = ∆ τ + μ ′ ∆ σ n (2) 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).

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): log R L D = 0.58 M w − 2.42   4.8 < M w < 7.6 (3) log R W = 0.41 M w − 1.61   4.8 < M w < 7.6 (4)

The strike–slip earthquake equations are as follows (coefficients are the average values): log R L D = 0.62 M w − 2.57   4.8 < M w < 8 (5) log R W = 0.27 M w − 0.76   4.8 < M w < 7.6 (6) M 0 = μ · R L D · R W · D ¯ (7) M W = 2 3 log ⁡ M 0 − 9.1 (8) 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; 2011; Wang et al., 2011).

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 ( ∆ C F S ) 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 ∆ C F S of the 2008 Wenchuan earthquake on the 2022 Lushan earthquake was 2.36 × 10 4 Pa, and the ∆ C F S at the time of the 2022 Lushan earthquake occurrence was 3.45×10⁴ 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 ∆ C F S (Figure 2I), so the 2013 Lushan earthquake had an unloading effect on the 2022 Lushan earthquake. The coseismic ∆ C F S of the 2013 Lushan on the 2022 Lushan earthquake was −4×10⁴ Pa, and this negative value then decreased to −4.53×10³ Pa. These data are consistent with the coseismic and postseismic stress change evolution of the 2013 Lushan earthquake (Figures 2D–I). The viscoelastic lower crust and upper mantle contributed to the occurrence of the 2022 Lushan earthquake.

The postseismic ∆ τ , ∆ σ n , and ∆ C F S contour maps show that the 2021 Maduo earthquake had loading effects on the 2022 Lushan earthquake that were similar to those of the Wenchuan earthquake, as shown in Figures 2J–L. The coseismic ∆ C F S of the Maduo earthquake on the 2022 Lushan earthquake was 9.26×10¹ Pa, and the postseismic ∆ C F S increased to 1.3×10² Pa at the time of the 2022 Lushan earthquake.

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, ∆ C F S was positive (Figures 2M–O). The overall ∆ C F S 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 ∆ C F S of historical earthquakes on the 2022 Lushan earthquake alone. As shown by the solid blue line (12 km) in (A) and (C), the ∆ C F S remained at a small positive value due to the great distance from the previous historical earthquakes. However, after the 2008 Wenchuan earthquake, the ∆ C F S 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 ∆ C F S 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.

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 cross-border conversion between positive and negative ∆ σ n , and its value was negative (Figure 4E). The hypocenter of MEK1 was in the ∆ C F S 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 ∆ C F S 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 ∆ C F S 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 ∆ C F S caused by the MEK1 and MEK2, indicating that the MEK1 and MEK2 may have jointly triggered the MEK3.

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.

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¹⁶∼10²¹ Pa s. Shao et al. (2011) considered the viscoelastic relaxation effect of the lower crust and upper mantle after the 2008 Wenchuan earthquake. They found that the best long-term viscosity coefficient of the lower crust and upper mantle in the Songpan–Ganzi Terrane was 5×10¹⁷ Pa s. Huang et al. (2014) used GPS and InSAR data to constrain the 1.5-year postseismic deformation of the 2008 Wenchuan earthquake and found that the long-term viscosity coefficient of the upper mantle of the Songpan–Ganzi Terrane was 1×10¹⁸ Pa s. Zhao et al. (2021) determined that the transient viscosity and steady viscosity of the lower crust and upper mantle in Tibet were 5×10¹⁸ Pa s and 4×10¹⁹ Pa s, respectively, through crustal deformation simulation after the 2001 Kekexili earthquake. Wei et al. (2020) estimated the viscosity range beneath the lower crust in the West Qinling–Songpan Terrane and its surrounding areas by a geomorphological analysis method and channel flow model, and the result was 10 18 ∼ 10 20 Pa s. Wang et al. (2021) strictly constrained the long-term viscosity of the lower crust and upper mantle in the Songpan–Ganzi Terrane to (4.3–5.7)×10¹⁸ Pa s and (1–1.6)×10¹⁸ Pa s according to the GPS data after the Wenchuan earthquake. Currently, the viscosity value of the lower crust and upper mantle in the eastern Tibetan Plateau remains controversial, as the viscosity at different time scales changes, but it is concentrated in the range of 10 18 ∼ 10 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 ∆ C F S 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 ∆ C F S . Therefore, the viscosity conditions considered in our study were appropriate.

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 ∆ C F S , but cannot affect the polarity of ∆ C F S .

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 ∆ C F S 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 ∆ C F S changes with the effective friction coefficient. However, the change in the effective friction coefficient in our study does not alter our conclusion.

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 ∆ C F S at the epicenter of the 2022 Lushan earthquake was calculated for depths of 12 km, 16 km, and 20 km. The ∆ C F S 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 ∆ C F S value changed from positive to negative as the depth increased. Combined with Figure 3A, the ∆ C F S 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, ∆ C F S may vary at different depths on the fault plane. Thus, the influence of different hypocenter depths requires further study.

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.

We calculated the Coulomb stress change along the faults around the Songpan–Ganzi Terrane after the Lushan–Maerkang earthquake sequence in 2022. The fault slip data were from previous publications (Deng, 2007; Xu et al., 2016, https://www.activefault-datacenter.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 ∆ C F S 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 ∆ C F S value, similar to the ∆ C F S 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 ∆ C F S changes sign, therefore potential earthquake risk may be high.