In this study, 105 earthquakes were selected from the catalog of the China Earthquake Network Center (CENC), including the 21 May 2021 Yangbi mainshock, with 11 foreshocks after 13 May and 39 aftershocks before 1 June, and the Maduo mainshock, with 53 aftershocks before 1 June. These earthquakes were selected according to magnitudes ≥ 3.0 to improve the quality of Lg recordings, and their locations and origin times were obtained from the catalog. We collected more than 42 thousand vertical-component seismograms recorded by 376 permanent seismic stations to investigate the source properties of the 2021 Yangbi and Maduo earthquake sequences. These stations have good azimuthal coverage and were within epicentral distances between 200–3,000 km to ensure well developed Lg waves (Figure 2).

The seismic Lg-wave, developed from a reverberated S-wave in the continental crust, is one of the prominent seismic phases in regional seismograms and often used in amplitude analysis (Furumura and Kennett, 1997; Fisk, 2006; Pasyanos et al., 2009). To analyze the arrival times and frequency contents of different regional phases, illustrated in Figure 3A is a regional waveform, and in Figure 3B is the corresponding spectrogram. The spectrogram was obtained using the short-time Fourier transform, where the window was determined to be consistent with the regular Lg-extraction length, according to the velocity window with a length of 0.6 km/s. Due to the changing length of the window, we regularized the spectrogram at 58 frequencies within 0.05–10.0 Hz, and then normalized the spectrogram at each frequency. The first arrival in the spectrogram was Pn at ∼2–5 Hz, followed by Pg over ∼0.3–4 Hz. The Sn was very weak compared with other regional phases. Lg could be recognized as a prominent patch with a longer duration and dominating energy within ∼0.2–2.0 Hz. Based on the spectrogram, the velocity window is 3.6–3.0 km/s. We followed the procedure of Zhao et al. (2010) to process the Lg amplitudes. First, instrument responses were removed, and Lg waveforms were cut using the velocity window of 3.6–3.0 km/s (Figure 3B). The window can be slightly adjusted to capture the maximum Lg energy. The waveforms immediately before the first p arrivals were defined as noise (Figure 3B). Then, the spectra of Lg and corresponding noise were calculated by applying a cosine taper and fast Fourier transformation and were sampled at 58 frequency points between 0.05 and 10.0 Hz to obtain the spectral amplitudes A o b s and A n o i s e of Lg and noise, respectively (Figure 3C). Next, the signal-to-noise ratio (SNR) was calculated by A o b s f / A n o i s e f at each frequency point f . We set the SNR threshold of 2.0, in consideration of the trade-off between the quality and amount of available data, to select the data to be used (Figure 3D). The quality control resulted in, for example, 18,039, 18,070 and 10,179 valid Lg spectral measurements at 0.5, 1.0 and 2.0 Hz, respectively. Finally, by assuming that the noise is not correlated with the Lg wave, the true Lg spectral amplitudes were obtained by A o b s 2 f − A n o i s e 2 f .

According to Aki (1967), we modeled the Lg spectral amplitude A f , ∆ as A f , ∆ = S f G ∆ Γ f , ∆ P f R f (1) where Δ is the epicentral distance, S is the source spectrum, G is the geometrical spreading term, Γ is the attenuation term, P is the site response and R is the term containing modeling errors. The geometrical spreading term can be written as Δ 0 ∆ − 1 2 , where Δ 0 is the reference distance to be fixed at 100 km (Street et al., 1975). The attenuation term has the form Γ f , ∆ = exp   − π f v ∫ d s Q x , y , f (2) where the symbol ∫ denotes the integral along the source-station ray path and Q can be obtained from the established high-resolution broadband Lg-wave Q models (Zhao et al., 2013a; Zhao et al., 2013b; He X. et al., 2021). To compute the attenuation term, we mapped the ray path onto the Q model and split it into N segments with a length of d s ; then, Eq. 2 can be discretized and rewritten as summation Γ f , ∆ = exp − π f v ∑ n = 1 N d s Q n   (3) where n denotes the n th segment of the ray, and Q n is the average Q within this segment. It is frequency-dependent and can be obtained from the broadband Lg-wave attenuation model. In previous studies, the relative site response has been determined in conjunction with attenuation (He X. et al., 2021) and is available for 375 stations in this study. For the only station without an available site term, we assigned it with a unit amplification. The Q model and site response dataset used are illustrated in Supplementary Figures S1, S2.

With known geometrical spreading, attenuation, and site response, we linearized Eq. 1 to relate source perturbations with the residuals between the observed and predicted Lg spectra (Zhao et al. (2010). The least squares QR (LSQR) factorization algorithm (Paige and Saunders, 1982) was used to obtain the source perturbations by minimizing the L2 norm of the residual. The source terms, starting with the unit source function, were iteratively updated until convergence. Based on the attenuation, site response and restored source terms, the synthetic Lg spectral amplitudes fit well to the observations. For example, the mean absolute error between synthetic and observed Lg amplitudes at 1 Hz decreased from 1.78 to 0.81 and the deviation from 1.27 to 0.84 (red and blue histograms in Figure 4). We also tested a constant Q to restore the source terms (yellow and green histograms in Figure 4), and found that Lg spectra can also be fitted to a certain degree, but not as good as that based on both laterally-varying Q model and site response.

After retrieving the source spectra S f for each earthquake at 58 individual frequencies, we estimate the scalar seismic moment M 0 and the corner frequency f c based on the source model S f = M 0 4 π ρ ν 3 1 + f f c γ n 1 γ (4) where n is the high-frequency falloff rate, γ controls the sharpness of the spectral corner, and ρ and ν are the density and shear wave velocity in the source region, respectively. Typical values of ρ = 2.7   k g / m 3 and ν = 3.5   k m / s were used. Generally, γ can be set to 1, although a larger value of γ = 2 was also used for some studies (Boatwright, 1980). Letting γ = 1 and n = 2 yields the commonly used Brune (1970) model. However, in this study, the high-frequency spectra often fall faster than 2. Therefore, we treated n as a variable to be estimated as well. By fitting the theoretical source model to the retrieved source spectrum using the annealing simulation algorithm (Kirkpatrick et al., 1983), M 0 , f c and n can be determined for each earthquake. The above results were also compared to those using other source models, including the Brune model ( γ = 1 ; n = 2 ), the ω − 3 model ( γ = 1 ; n = 3 ) and the Boatwright model ( γ = 2 ; n as variable) (refer to Supplementary Figure S3)

The ideal scenario is that the source spectrum has a flat low-frequency asymptote and a high-frequency descending portion, permitting to accurately estimate M 0 and f c . In practice, our analysis showed that for small earthquakes, the Lg source spectra at frequencies lower than ∼0.2 Hz were sometimes raised, possibly by ambient noise or surface waves generated by preceding earthquakes. Although not apparent for larger earthquakes because their source spectra were significantly higher than the noise level, the contaminated step-like spectrum would possibly result in an underestimated corner frequency and overestimated seismic moment for small earthquakes, introducing large errors into the following analysis. To suppress these effects and improve the reliability of parameter estimation, we set the fitting band truncations for small earthquakes and then screened problematic estimations based on comparisons of magnitude measurements. First, for M ≤ 3.5 earthquakes, we fitted the retrieved source spectra within only 0.3–10.0 Hz because the observed noise level is high at approximately 0.2 Hz and the corner frequencies for earthquakes with magnitudes between 3.0 and 3.5 are generally higher than 0.3 Hz (e.g., Madariaga, 1976). The fitting band of 0.3–10.0 Hz could generally exclude the noise spectrum and retain the high-frequency portion of source spectra as much as possible to estimate the corner frequencies. Such a magnitude-dependent fitting band strategy was also used in previous source analysis studies (Abercrombie, 1995; Allmann and Shearer, 2007; Allmann and Shearer, 2009). Although effective in finding M 0 and f c for small earthquakes, such rough truncation would sometimes omit all of the flat portion, leading to a problematic estimation. Therefore, next, we screened these problematic results based on the obtained seismic moment M 0 . The moment magnitude M W can be calculated from M 0 using the China Earthquake Administration scaling relationship M W = 2 3 log ⁡ M 0 − 9.1 ( M 0 is in the unit of N m). Magnitude scaling analysis based on earthquakes in Switzerland demonstrated that the difference between M W and the local magnitude M L , M W − M L generally ranged from −1.0 to 0.5 (Goertz-Allmann et al., 2011). Considering the geological settings and algorithm used to calculate M 0 , we set a broader threshold and discard earthquakes with M L − M W > 1.0 (Figure 5). In total, 17 earthquakes were discarded due to abnormally high M 0 values. The remaining earthquakes yielded a linear relationship M L = 1.31 ± 0.02 M W − 1.64 ± 0.12 . The slope was in agreement with previous studies (e.g., Bethmann et al., 2011), whereas slight differences could result from geological settings, algorithms to calculate the M W or the uncertainty in the catalog magnitude.

Based on a circular fault model (Eshelby, 1957), the stress drop can be calculated using Δ σ = 7 M 0 16 f c k v 3 (5) where the constant k is 0.21 by assuming a rupture velocity of 0.9 shear wave velocity (Madariaga, 1976). The constant k was also suggested to be 0.29 and 0.37, depending on specific source model used (Brune, 1970; Sato and Hirasawa, 1973). We also tested the effects of using different k values, and found that the absolute stress drop values may vary with k, but their relative relationships are kept unchanged among events (Supplementary Figure S4). For simplicity, stress drops for all events were calculated based on same assumptions, including using the circular fault model, although it may be improper for large earthquakes with excessive long ruptures. Therefore, the mainshock stress drops will not be compared with aftershock stress drops.

To refine the spatial distribution of the stress drop, we relocated the 2021 Yangbi and Maduo earthquake sequences for all earthquakes with local magnitudes greater than 3.0 using the HypoDD program, which allows for simultaneous accurate relocation of clusters of earthquakes based on double-difference travel times (Waldhauser and Ellsworth, 2000). Earthquake relocation and stress drop calculation based on consistent catalog avoid possible errors of assigning the calculated stress drops to incorrect events. We first visually inspected first P-wave arrivals automatically picked by the deep-learning-based PickNet program (Left panels in Figure 6) (Wang et al., 2019). This program can offer better picks at regional distances and was previously trained by the travel-time dataset from earthquakes in the Myanmar region (Wang et al., 2019). Then, we further processed the picked arrival times to calculate the double-difference first-P travel times of earthquake pairs at their common stations and only those pairs observed by more than eight stations were used. For the Yangbi earthquake sequence, a total of 1,493 first arrival P-wave travel times with epicentral distances less than 5° generated 20,124 P-phase pairs. The earthquake pairs had a separation of less than 10 km. For the Maduo earthquake sequence, a total of 2,594 first arrival P-wave travel times at 214 stations were used to generate 3,872 P-phase pairs. When processing the Maduo earthquake sequence, the performance of PickNet significantly reduced the amount of available data, and therefore, we allowed larger observation distances of 10°. Due to its ∼150-km-long aftershock zone, we also allowed a larger spatial separation of 20 km within each earthquake pair so to include as many earthquakes as possible into the relocation clusters. The station distributions have relatively good azimuthal coverage for both earthquake sequences (Supplementary Figure S5). Proper parameters were chosen through several trials to control the initial conditions, clustering of earthquakes and data weighting in the HypoDD program. During relocation, local 1D crustal structures for the source regions of Yangbi and Maduo sequences were extracted from the USTClitho1.0 model (Right panels in Figure 6) (Xin et al., 2018). The relocation results were obtained using the LSQR factorization algorithm (LSQR, Paige and Saunders, 1982) with eight iterations. For convergence purpose, in each iteration, the data outside of 3 times the standard deviation of the RMS residual were treated as outliers and dropped. The iterations updated epicenter locations. Meanwhile, the rms travel time residuals were reduced from 0.49 s to 0.28 s and from 0.67 s to 0.25 s for Yangbi and Maduo sequences, respectively (Supplementary Figure S6). To assess relocation errors, we relocated the events for 100 times, in each time, with only 80% randomly-selected double-different travel times being used. For the Yangbi sequence, the mean errors were 245, 200 and 558 m in the east-west, north-south and vertical directions, respectively (Supplementary Figure S7). For the Maduo sequence, the mean errors were 812, 667 and 1,053 m in the east-west, north-south and vertical directions, respectively (Supplementary Figure S8). The errors were comparable to those of previous studies (Long et al., 2021; Xu et al., 2021; Zhang J. Y. et al., 2022; Zhu et al., 2022). The accuracy and reliability of earthquake relocation could be further improved by incorporating the cross-correlation double-different travel times, which is, however, out of the scope of this study.

The stress drop of the Yangbi mainshock is 22.59 ± 4.85 MPa. For the foreshocks, the stress drops range from 0.44 to 11.56 MPa, with a median value of 0.59 MPa, and for aftershocks ranged from 0.08 to 5.55 MPa, with a median value of 0.75 MPa (Figures 8A,B). The foreshocks emerged on 18 May, with their stress drops exhibiting a general decrease from the initial value, until foreshocks became active again on 19 May. After that, a gradual increase in stress drop was visible, followed by two earthquakes, including the largest Ms 5.6 one, characterized by a significantly higher stress drops of ∼10 MPa occurring ∼25–27 min before the mainshock. In the aftershock sequence, the stress drops were relatively fluctuated and decreased gradually along with the decrease of seismicity. The minimum stress drop appeared ∼26 min after the mainshock, indicating of a great amount of stress was released by the mainshock. The stress drop was again increased along with several strong aftershocks, including the largest Ms 5.2 aftershock occurring ∼43 min after the mainshock, with the highest stress drop among all aftershocks. After strong fluctuation lasting for approximately 1 day, the stress drop values remained around the median value, which could indicate the stress adjustment following the mainshock.

For the Maduo earthquake sequence, the stress drop of the mainshock is 28.31 ± 3.92 MPa, and the stress drops of aftershocks range from 0.05 to 3.21 MPa, with a median value of 0.51 MPa (Figure 8C). The values of the stress drop for the Maduo aftershock sequence are generally lower than those for the Yangbi sequence in terms of their maximum, minimum and median values, suggesting a possibly lower stress level on the Maduo seismogenic fault after the mainshock. The Maduo aftershock sequence is characterized by more scattered stress drop values undulated for ∼9 days after the mainshock. For example, several large aftershocks occurred ∼8 h and ∼8 days after the mainshock, with their stress drops greater than 1.5 MPa. Compared with the Yangbi sequence, the prolonged undulation in the stress drop and late occurrence of the minimum value nearly 3 days after the mainshock probably indicate a relatively long stress repartition process, which could be related to the longer seismogenic fault involved in the Maduo sequence.

Based on the first P arrivals picked by the PickNet program, we used the double-difference earthquake relocation program HypoDD to refine the locations of the calculated stress drop. The results were generally consistent to previous studies (e.g., Long et al., 2021; Xu et al., 2021; Zhang et al., 2021; Zhang J. Y. et al., 2022). For the May 2021 Yangbi earthquake sequence, we relocated 51 epicenters (Figure 9; Supplementary Table S1). The Yangbi earthquake sequence was generally distributed along a northwest‒southeast trend and extended for approximately 21 km. Although roughly parallel to regional tectonics, e.g., the northwest‒southeast-trending Red River fault and Jinshajiang fault, the distribution of the Yangbi earthquake sequence appeared not directly linked to any known faults, including the Weishan–Qiaohou fault ∼10 km to the northeast. The foreshock activities migrated back and forth within the epicentral area, including the largest M5.6 foreshock near the center and the M4.5 foreshock near the southeastern end occurring ∼2 min later. The Yangbi mainshock was located at 99.8635 ° E and 25.7039 ° N , with a depth of 9.090 km at the northwestern end of the epicentral area, suggesting a unilateral rupture toward the southeast. For the aftershocks, their epicenters extended further southeast and gradually scattered to a wider area, possibly indicating secondary structures with variable orientations. An M4.1 event was located notably off the primary aftershock trend to the northeast, which may indicate a cluster of earthquakes on an activated secondary fault (Zhang et al., 2021). The sources of the Yangbi earthquake sequence were mostly located in the upper crust at approximately 10 km depths.

By relocation, we obtained the fine distribution of the stress drop. The highest stress drop in the foreshock sequence was observed for a M5.6 event located near the center of the foreshock zone, where foreshock activities had been initiated since 18 May. Another significantly high stress drop was linked to the M4.5 foreshock, occurred at the southeast end of the foreshock area. The mainshock and the largest aftershock were located at the northwestern and southeastern ends of the epicentral zone, respectively, both with notably high stress drops.

We relocated 54 earthquakes in the Maduo sequence (Figure 10; Supplementary Table S1). The results demonstrate that the mainshock was located at 98.3641 ° E and 34.6656 ° N , with a depth of 21.909 km. The epicenter of the mainshock was ∼5 km northeast to the mapped trace of the Kunlunshankou–Jiangcuo fault (KJF). The aftershocks were nearly linearly distributed on both sides of the mainshock, extending approximately 170 km along a northwest–west trends and generally consisting with the KJF. Gaps with absent or sparse aftershocks are visible, for example, at sections ∼40–60 km (gap 1), ∼15–20 km northwest (gap 2), and ∼30–50 km southeast (gap 3) of the mainshock (Figure 10B). The aftershock zone intersected with the Maduo–Gande fault (MGF) near its southeast end, with the strike turning from northwest–west to west–east. The relocated focal depths in the Maduo sequence mainly ranged between ∼5–15 km, which were generally deeper than those in the Yangbi sequence and appeared more scattered. Except for the mainshock, which has the highest stress drop, other events with significantly high stress drops (e.g., >10^3.5 kPa) in the Maduo sequence can be generally associated with the aftershock gaps, especially the cluster ∼30 km northwest of the mainshock (gap 1). The cluster of high stress drops ∼60 km southeast of the mainshock may also be related to the intersection with the MGF.

The earthquake stress drop can be related to the localized stress variability on a seismogenic fault. Larger stress drop variations could result from fault plane heterogeneities (Kirkpatrick et al., 2020), e.g., a rougher fault plane may cause higher stress drops due to its additional resistance to slip (Fang and Dunham, 2013). Allmann and Shearer (2009) observed higher overall stress drops for strike-slip earthquakes at a global scale, which is consistent with previous findings of increased apparent stress for strike-slip events. However, at local scales, the dependence on the focal mechanism is less significant, and stress drops appear more dominated by the local tectonics (Shearer et al., 2006; Goebel et al., 2015; Boyd et al., 2017). Therefore, investigating stress drops in the Yangbi and Maduo earthquake sequences may reveal stress variations in the related seismogenic zones.

In the Yangbi sequence, the stress drops of foreshocks may suggest the accumulation of stress and hint the stress variation that led to the nucleation of the mainshock. The pre-slip mechanism is one of the popular explanations for how an earthquake initiates, in which concentrated and repeated seismicity is expected (Zhu et al., 2022). However, the relocated foreshock sequence revealed a back-and-forth migration pattern. In addition, high-stress-drop foreshocks, including the first and the largest ones, have their stress drops significantly higher than the median value of the entire sequence, implying the potential to trigger earthquake. By modeling the cumulative static stress change from large foreshocks, Zhu et al. (2022) found that the local shear stress at the Yangbi mainshock epicenter was increased by major foreshocks. Therefore, our results of foreshock relocation and stress drops were consistent with the cascade triggering process during the nucleation of the Yangbi mainshock (Long et al., 2021; Zhang et al., 2021; Zhu et al., 2022). For the aftershock sequence, fluctuations of stress drops including many significantly high and low stress drop events may be related to the stress adjustment and continually release of the local stress field disturbed by the mainshock. Similar stress development pattern were also suggested by the observed b-values, which demonstrated long-term low values before 20 May 2021, showed abrupt reductions before the Ms 5.6 foreshock and the Ms 6.4 mainshock, then demonstrated large fluctuations around their mean value after the mainshock (Duan et al., 2021; Jiang et al., 2021). The b-value measures the distribution of earthquakes with different sizes in the recurrence relation of Gutenberg and Richter (1941) and is found to be reversely related to the differential stress within the crust (Schorlemmer et al., 2005). A low b-value and high stress were also observed preceding the M6.4 and M7.1 Ridgecrest, California earthquakes in 2019 (Nanjo, 2020).

In the Maduo sequence, the aftershocks were distributed on both the northwest and southeast sides of the mainshock along the Kunlunshankou–Jiangcuo fault. For aftershocks to the east and west of the mainshock, our results showed their median stress drops were 0.49 MPa and 0.62 MPa, respectively. The difference probably suggested differentiated rupture processes between the east- and west-sections. Previous studies have suggested that the eastward rupture of the Maduo mainshock propagated at a speed exceeding 3.5 km/s S-wave velocity, i.e., supershear rupture (Zhang X. et al., 2022; Li et al., 2022; Lyu et al., 2022; Yue et al., 2022). The back-projection imaging technique is commonly used to obtain large earthquake rupture kinematics. With this method, the Maduo mainshock was found ruptured toward the east and west, with velocities of 3–5 km/s and 2–3 km/s, respectively (Zhang X. et al., 2022; Li et al., 2022; Yue et al., 2022). Zhang X. et al. (2022) further confirmed the eastward supershear rupture by recognizing far-field Love Mach waves. Joint analysis using Global Positioning System, Interferometric Synthetic Aperture Radar, and teleseismic data determined the westward and eastward rupture velocities to be 2.0 and 4.6 km/s, respectively (Yue et al., 2022). Near-field Global Navigation Satellite Systems recordings demonstrated significant sensitivity to determine rupture velocities and confirmed sub- and supershear ruptures on the west and east sections of the mainshock fracture, respectively (Lyu et al., 2022). We suggest that the lower aftershock stress drops on the east section can be related to the larger stress release due to the supershear rupture. Yue et al. (2022) pointed out that the coseismic slip east of the mainshock was concentrated in seismogenic layers, in which rate-weakening friction law might support instability failure to facilitate the supershear rupture. Since a large amount of energy was consumed to feed the supershear rupture, the remaining stress was low, resulting in low stress drops for the aftershocks to the east of the mainshock. The relatively high aftershock stress drops observed on the west section of the mainshock fracture suggest there were more strain stored within this part of the crust. This energy can be further released by later aftershocks or interseismic creeping.

In this study, stress drops from the relocated earthquake sequences demonstrated that high stress drops were often related to fault junctions or earthquake gaps, suggesting stress concentrations at these locations.

The relocated sources of the Yangbi sequence suggested a complicated fault system with a northwest‒southeast-trending major fault, which may split into multiple branches in the southeast section. This is consistent with previous earthquake relocations; however, a detailed fault structure can be suggested by relocation for a larger number of small earthquakes. Long et al. (2021) relocated more than 1,000 Yangbi events with local magnitudes greater than 1.0 based on the seismic phase bulletin, suggesting a northwest‒southeast-trending primary fault, with northeast‒southwest-trending conjugate faults in its northwestern segment, and multiple secondary faults along the primary fault and gradual turning toward the south in the southeast. The multiscale en echelon structures were further revealed by high-resolution catalogs built using a deep-learning phase picker (Su et al., 2021; Liu et al., 2022). Combining the fault structures with our relocation results, high-stress-drop events in the Yangbi sequence were generally related to fault junctions and tips. For example, the M4.5 and M5.6 foreshocks with significantly high stress drops were located around the middle of the epicentral zone, where the south-trending en echelon faults intersect with the primary fault; several aftershocks with high stress drops could be found near fault junctions and tips, including the M5.2 and M4.1 aftershocks. Simulations suggested significantly increased stress localized at the junctions of branching faults (e.g., Ando et al., 2018). Therefore, the high stress drops in the Yangbi sequence may result from the stress concentrations near the fault junctions and tips and may in turn indicate the precise locations of sophisticated fault branches.

In the Maduo sequence, the stress drops demonstrated a more fluctuating pattern than in the Yangbi sequence. The prolonged stress drop variations could be related to the fault properties. The cross-section view reveals that the clusters of earthquakes with significantly high-stress-drop values were generally associated with the earthquake gaps along the fault. These gaps should be robust features despite a magnitude threshold 3.0 was used for earthquake relocation in this study. In previous studies that relocated more than 1,300 earthquakes with magnitudes greater than 0, reduced seismicity was still pronounced at these locations, and in particular, few aftershocks occurred ∼50 km northwest and ∼30 km southeast of the mainshock (Wang W. L. et al., 2021; Xu et al., 2021) (corresponding to the gaps 1 and 3 in Figure 10B). The correspondence between high stress drops and aftershock gaps probably indicated redistributed stress following the mainshock. Geodetic observations revealed coseismic slip on the seismogenic fault following the mainshock. They varied significantly along the fault strike and concentrated at several localized segments with maximum slips of ∼3–4 m (He K. et al., 2021). These high-slip patches may represent the asperities where stresses were released by significant rupture during the mainshock, resulting in localized reduced or lack of aftershock activities. Since the asperities are generally related to variations in strength or friction properties due to the changes of fault geometry, and result in stress concentration (Mildon et al., 2019), we propose that the stress remaining on the fault was redistributed and concentrated near the edges of fault asperities, thus leading to the occurrence of larger aftershocks with high stress drops. A similar complementary correspondence between the stress drop and mainshock slip was also suggested for the 1992 Landers earthquake in southern California (Shearer et al., 2006).

The stress drops of the Yangbi, Yunnan and Maduo, Qinghai earthquake sequences in May 2021 demonstrated different temporal and spatial patterns, which reflected varying localized crustal stress variability and deformations across the eastern Tibetan Plateau. A deformation field recently constructed from GPS data revealed predominantly continuous strains inside the eastern Tibetan Plateau, with significant shear associated with large strike-slip faults to accommodate the eastward extrusion of geological blocks (Wang and Shen, 2020). On the southern flank of the extrusion, the Chuandian block represents the majority of the escaping crustal material from the interior plateau. Wang and Shen (2020) found that the translation of the Chuandian block was accompanied by notable clockwise rotation at a rate of ∼2.9°/Myr around the Eastern Himalayan syntax and widespread dilatation at 10–20 nanostrain/yr. The dilatation should be responsible for the normal component obtained in the focal mechanism solutions of the Yangbi earthquake sequence. The complicated combination of southeastward translation, clockwise rotation and interior deformations requires a large number of active faults to accommodate, possibly resulting in the multiscale en echelon faults and conjugate faults that inherited motions of the associated large strike-slip faults. Therefore, these small faults could also generate large earthquakes, including the 2021 Yangbi Ms 6.4 earthquake, with the repartitioned stress concentration subject to the geometric complexity of the fault systems, for example, the fault junctions and tips, giving rise to several high stress drops in the Yangbi earthquake sequence. The local stress field inverted from focal mechanisms suggested a change in the orientation of the horizontal maximum principal stress from near north–south to nearly east‒west in the southernmost epicentral zone, coincident with the regional stress and possibly reflecting the tectonic rotation of the Chuandian block (Wang Y. et al., 2021).

In contrast, the focal mechanism of the Maduo earthquake sequence demonstrated a consistent northwest‒southeast orientation of maximum principal stress, corresponding to the movement of the Bayan Har block (Xu et al., 2021). The Bayan Har block is bounded in the north by the Kunlun fault, with significant left-lateral strike-slip, to accommodate the northeastward expansion of the eastern Tibetan Plateau. The slip rate on the Kunlun fault shows a west–east decrease from ∼12.5 mm/yr to less than 3 mm/yr, which was suggested to be decomposed to slips on several secondary faults subparallel to and south of the Kunlun fault (Kirby et al., 2007), including the Kunlunshankou–Jiangcuo fault, which was the seismogenic fault of the 2021 Maduo earthquake sequence. The slip around the middle segment of the Kunlun fault, which is faster than that near the northern segment of the Red River fault, is more likely to result in a concentration of shear strain (Zheng et al., 2017; Wang and Shen, 2020). Localized variations in crustal strength could occur due to shear heating and increased crustal temperature (Leloup et al., 1999; Bao et al., 2012). As a consequence, large earthquakes would rupture the high-strength asperities, and the perturbed stresses were likely to be repartitioned around the asperities and lead to larger variability in the stress drops of aftershocks. Therefore, the stress drops of the Yangbi and Maduo earthquake sequences were controlled by the eastward expansion of the Tibetan Plateau, and the differences in their patterns showed varying impacts on the seismogenic environment.