Source Model for the 2022 Qinghai Menyuan Ms 6.9 Earthquake Based on D-InSAR

Introduction

At 01:45 on 8 January 2022, a magnitude 6.9 earthquake occurred in Menyuan County, Haibei Prefecture, Qinghai Province (37.77 N, 101.26 E) with a focal depth of 10 km. The epicenter is 54 km away from Menyuan County, with an average altitude of about 3,674 m within 5 km (http://www.cenc.ac.cn/). The area is sparsely populated. Nine people were injured in the earthquake. Many buildings, including the Lanxin Railway, were severely damaged with a maximum intensity of IX (China Earthquake Administration; Lu et al., 2022). As of April 1, two aftershocks of magnitude 5 or above were recorded, which occurred 26 min and 4 days after the earthquake, respectively, and the Ms 5.2 earthquake was the largest aftershock on January 12.

This earthquake occurred on the Qilian-Haiyuan fault zone on the northeastern margin of the Qinghai-Tibet Plateau (Figure 1). Since the Cenozoic, the northeastern margin of the Qinghai-Tibet Plateau has expanded into the mainland under the double effects of the continuous northeast compression of the Indian plate and the resistance of the adjacent strong Alashan block (Zhang et al., 2003; Shi et al., 2018). The GPS velocity field shows that the crustal displacement of the northeastern margin of the Qinghai-Tibet Plateau is NE-EW, and the continuous northward expansion is the major tectonic cause for the frequent occurrence of strong earthquakes in this area (Wang and Shen, 2020). Qilian-Haiyuan fault is a large active fault in the study area, which is dominated by strike-slip. It starts from the Muli-Jiangcang fault in the west, and stretches eastward as the Tuolaishan fault, the Lenglongling fault, the Jinan fault, the Qianghe fault, the Maomaoshan fault, the Laohushan fault, the Haiyuan fault, the Liupanshan fault. Since 1900, two M ≥ 8 earthquakes have occurred on this fault—Haiyuan M 8^1/2 earthquake in 1920 and the Gulang M 8 earthquake in 1927 (Zhu et al., 2022). The Menyuan Ms 6.9 earthquake occurred in a region with a risk code of A3 in the middle part of Qilian Mountains (Xu et al., 2017). It is located at the intersection of the Lenglongling fault and the Tolaishan fault (Fan et al., 2022; Xu et al., 2022; Yang et al., 2022). The Lenglongling fault is connected to the Jinqianghe fault in the east and the Tolaishan fault in the west. It is a Holocene sinistral strike-slip and thrust fault with a strike of NWW and a total length of about 150 km (Xu et al., 2016). On 26 August1986 and 21 January2016, two historical earthquakes of M 6.4 occurred on a blind thrust fault developed on the northeast of Lenglongling fault, at 28 and 33 km from this epicenter, respectively (Liu et al., 2019; He et al., 2020).

FIGURE 1. Tectonic backgrounds of the Menyuan Ms 6.9 earthquake. (A) The green and blue dashed rectangles represent the data coverage of ascending and descending orbits of Sentinel-1 satellites, respectively; the black rectangle represents the scope of the study area. The black and red solid lines constitute together the Qilian-Haiyuan Fault, and the latter is seen as the Tianzhu Gap (Gaudemer et al., 1995). The yellow circles represent historical earthquakes with surface wave magnitude (Ms) larger than 7 since 1900 (data from Sichuan Earthquake Administration). Black beach balls represent focal mechanisms of M ≥ 5.0 events (Global CMT, since 1976). Brown arrows represent the GPS velocity field with respect to the Eurasian frame (Wang and Shen, 2020). The red star represents the Ms 6.9 earthquake. The background image is DEM from SRTM (Farr and Kobrick, 2000). (B) Topography of the study area [black box in Figure (A)] and aftershock distribution. The mainshock is marked by the red star. Pink circles represent two Ms ≥ 5.0 aftershocks (Yang et al., 2022); black beach balls represent the focal mechanisms of two Ms 6.4 Menyuan earthquakes in 1986 and 2016. The solid black line represents the surface fault trace (Xu et al., 2016). (C) Geographical location of the study area. The red box marks the extent of figure (A). MJF: Muli-Jiangcang Fault, TLSF: Tuolaishan Fault, LLLF: Lenglongling Fault, JQHF: Jinqianghe Fault, MMSF: Maomaoshan Fault, LHSF: Laohushan Fault, HYF: Haiyuan Fault, LPSF: Liupanshan Fault, SNQLF: Sunan-Qilian Fault, MYDMYF: Minyue-Damaying Fault, HCSTF: Huangcheng-Shuangta Fault, MYF: Menyuan Fault, DBSF: Dabanshan Fault, RYSF: Riyueshan Fault.

After the Menyuan Ms 6.9 earthquake, some scientific research teams responded quickly, and made significant research achievements in the field geological surveys, precise location of aftershocks, focal mechanism solutions, fault slip distribution, regional Coulomb stress changes, regional deep structures and so on (Feng et al., 2022; Liang et al., 2022; Pan et al., 2022; Peng et al., 2022; Wang et al., 2022; Yang et al., 2022; Zhu et al., 2022). Field investigation reveals that the surface rupture mainly includes linear shear cracks, oblique tension cracks, tension shear cracks, extrusion bulge and other types. It is the rupture of the left-order part of the left-lateral strike-slip fault under the tensile stress regime, and it is considered that the fault is ruptured bilaterally initiating from the compressional bend (Li et al., 2022a; Pan et al., 2022). The aftershock location results show a distinct spatiotemporal migration exist from the west to the east (Fan et al., 2022; Xu et al., 2022; Yang et al., 2022). Focal mechanism results show that the earthquake is a strike-slip type. Further, the source characteristics of aftershocks vary from the mainshock zone to the east side. Focal mechanisms of aftershocks near the mainshock are similar to that of the mainshock, and the reverse component increases to the east side obviously (Liang et al., 2022). The finite fault inversion results show that the fault is dominated by left-lateral strike-slip with a maximum slip of 3.5 m. The Coulomb stress change results show that the earthquake is located at the positive stress zone of the 2016 Menyuan earthquake (Li et al., 2022b; Peng et al., 2022). The historical strong earthquakes on the northeastern margin of the Qinghai-Tibet Plateau have a stress loading effect on the Muli-Jiangcang and Tuolaishan faults in the western segment of the Qilian-Haiyuan fault zone, as well as the Jinqianghe-Laohushan fault in the middle segment (Zhu et al., 2022). In addition, the thickness of crust changes dramatically in the eastern part of the Qilian block, and this earthquake is located at the transition zone where physical properties of the crustal medium rapidly change (Wang et al., 2022).

Results above show that the fault structure is complex. Some problems still remained unclear. For example, the influence of the earthquake on the adjacent faults is not resolved. Previous studies mainly focused on quick disaster emergency response after the earthquake. The ascending orbit images used in the InSAR-based deformation observation only covered part of the rupture area, so the interferogram was incomplete. Would the results be improved by using ascending orbit data that completely covered the rupture zone?

In this paper, the co-seismic displacement field is obtained by D-InSAR technology, using the data of synthetic aperture radar (SAR) in ascending and descending orbits covering the seismic rupture area. At the same time, based on Pixel Tracking technology, the range and azimuth offset maps of the ascending and descending orbits are obtained, and the surface fault traces are also extracted. Then, the linear inversion of the distributed slip model is carried out. Based on the finite fault model, the Coulomb stress change caused by the earthquake are calculated. Finally, and the stress condition of the adjacent faults after the earthquake is analyzed.

InSAR Data and Processing

According to the principle of minimum time interval of SAR data coverage, the Level 1 Single-Look Complex data (SLC) in Sentinel-1 Interferometric Wide (IW) mode is used. The ascending orbit (T128A) and the descending orbit (T33D) each contain two images to form an interferometry pair. The detailed parameters of the image pairs are shown in Table 1. The time baselines of the satellite images of the ascending and descending orbits are both 12 days. The ascending orbit image pair were collected 3 days before and 9 days after the earthquake. The descending orbit image pair were collected 10 days before and 2 days after the earthquake. So, more afterslip is included in the ascending orbit pair. In addition, the corresponding spatial perpendicular baselines of the two pairs of images are both less than 60 m, indicating that the topographic phase has little influence on the interferometry phase, increasing the sensitivity of the deformation phase.

TABLE 1. Basic parameters of the SAR data.

In this paper, the open-source software ISCE (InSAR Scientific Computing Environment, Rosen et al., 2012) is used for two-pass differential interference processing of the SAR data in ascending and descending orbits. Firstly, the software achieves accurate image registration through enhanced spectral diversity (ESD) algorithm (Fattahi et al., 2017), which makes the azimuth matching accuracy higher than 0.001 pixels, avoiding the phase jump caused by adjacent image bursts, and greatly improving the reliability of interference results. The precise orbit (POEORB) data of European Space Agency (ESA) is used to eliminate orbit errors (https://scihub.copernicus.eu/). The 30-m resolution SRTM elevation data product published by NASA was used to simulate the terrain phase, and the interferogram is obtained by the second difference. In order to improve the signal-to-noise ratio of the interferograms, the multi-look ratio of range direction and azimuth direction is set as 5:1. A weighted power spectrum adaptive algorithm for phase filtering (Goldstein and Werner, 1998). Then, the phase unwrapping is carried out based on the minimum cost flow algorithm, to obtain the geocoded Line Of Sight (LOS) displacement. The atmospheric errors in InSAR measurement mainly come from the influence of the ionosphere and troposphere, and the ionosphere has less influence on the C band due to the relatively short wavelength (Gray et al., 2000). Therefore, the tropospheric atmospheric delay model provided by the General Atmospheric Correction Online Service (GACOS) of Newcastle University was used for atmospheric correction (Yu et al., 2018). The residual error of orbit will produce long wavelength noise, so a linear trend is removed from the line of sight (LOS) displacement result. Finally, the co-seismic LOS deformation field after atmospheric correction and detrending is obtained. (Supplementary Figure S1).

Coseismic Deformation

In this paper, D-InSAR technology was used to obtain the interferograms of the ascending and descending orbit SAR data (Figure 2). Except for a small incoherent area in the deformation center caused by deformation jump on the fault trace, strong coherence is preserved among the most area, showing continuous interference fringes. Both ascending and descending interferograms show the typical “butterfly shape” of strike slip rupture (Figures 2A,C). The deformation zone is about 50 × 40 km², mainly along the Lenglongling fault and the eastern section of Tuolaishan fault.

FIGURE 2. Interferogram and co-seismic deformation field (LOS direction) of Menyuan Ms 6.9 earthquake. (A,C) are the interferograms of T128 ascending orbit and T33 descending orbit, respectively. Focal mechanism solutions from GCMT, GFZ, and USGS are displayed with the red, blue, and green focal spheres. The red star marks the epicenter of the mainshock. (B,D) are the LOS co-seismic deformation fields of T128 ascending orbit and T33 descending orbit, respectively. The dashed lines labeled by (A–D) are four cross-fault profiles shown in Figure 3, and the red solid lines represent the seismogenic faults. Black arrows indicate satellite flight direction (AZI) and right-looking observation direction (LOS).

For an earthquake, the co-seismic deformation fields of ascending and descending orbits are different, which is mainly caused by the difference of incidence angle and flight azimuth (Hu et al., 2010). The relation between three component of deformation and the LOS displacement is:

Where $u_{los}$, $u_e$, $u_n$, $u_z$ are the LOS, east, north and vertical displacement respectively; $\theta$ are the incidence angle of radar pulse; $\alpha$ represents the azimuth angle of satellite flight. According to the average incidence and azimuth angle of ascending and descending tracks (Table 1), the contribution values of the east-west, north-south and vertical displacements to the LOS displacement of ascending and descending orbits are (−0.5770, −0.1065, 0.8097) and (0.6072, −0.1095, 0.7869), respectively. It can be seen that the Sentinel-1 satellite is most sensitive to vertical deformation, but insensitive to north-south deformation. And the opposite characteristics and great differences between the ascending and descending interferograms are mainly due to east-west displacement. The results of T128 ascending orbit show that the displacement at the northern side of seismogenic fault are mainly positive with a maximum value of 42.72 cm (LOS direction), and are mainly negative at the southern side with a minimum value of -59.65 cm (LOS direction). On the contrary, the displacement results of T33 descending orbit show that negative values are mainly observed at the northern side (−71.52 cm) and positive values is mainly observed at the southern side (94.65 cm) of the fault. Considering the geometric characteristics of satellite flight and right-sided look and LOS changes, the earthquake is associated with sinistral strike-slip faulting. The absolute LOS displacement on the south side of the fault is greater than that on the north side, which may result from slip on dipping fault plane as suggested by focal mechanism solutions. Meanwhile, the focal mechanism results show that the earthquake contains a small amount of thrust component (Table 2). Thus, it could be inferred that the south side of the fault is the hanging wall, and the north side is the footwall, slightly inclined to SW.

TABLE 2. Focal mechanism solutions of the Menyuan Ms 6.9 earthquake.

Figure 3 shows four displacement profiles in the ascending and descending tracks. The displacement gradually decreases away from the fault. Compared with the AA’ and DD’ sections far from the epicenter, BB “and CC” sections near the epicenter have obviously larger relative displacements, indicating that the ground deformation caused by the earthquake decreases from the center to both sides. The surface displacement caused by the Menyuan Ms 6.9 earthquake gradually decreased from the center of the rupture zone to both ends, and the displacement on both sides of the cross-fault was opposite in direction but similar in size, indicating that the seismogenic fault has a high dip Angle.

FIGURE 3. Co-seismic displacement profiles of the Menyuan Ms 6.9 earthquake (LOS direction). (A) Cross-fault displacement profiles of the T128 ascending orbit, with profile location shown in Figure 2; (B) Cross-fault displacement profiles of the T33 descending orbit. Gray circles are the data within 1 km along the profile. The vertical gray dotted line indicates the location of the faults, of which the southern branch fault (fault2) is located about 1 km south from the main fault.

Fault Geometry Constraints and Slip Distribution Inversion

After the Menyuan Ms 6.9 earthquake, different institutions and researchers quickly conducted focal mechanism inversion, and the results were somewhat different (Table 2). Xu et al. (2022) and Liang et al. (2022) have used more waveform data to obtain a more reliable solution of focal mechanism, which shows that the focal depth was shallow, both at 4 km. Compared with the teleseismic records, InSAR observations can provide more near-field constraints for the inversion of fault geometry and slip distribution. In order to determine the fine motion characteristics of seismogenic faults, the down-sampling process is carried out on the basis of co-seismic surface displacement, and then the slip distribution of the fault is inverted based on the classical elastic dislocation model.

InSAR co-seismic deformation field usually contains millions of points, and adjacent observation points have strong spatial continuity (Hanssen, 2001; Ji et al., 2017). In order to improve the calculation efficiency and reduce the influence of unreliable values in inversion, filter and downsampling preprocess are applied to the LOS displacement data. Firstly, the interferograms of ascending and descending orbit are masked according to the coherence. Sampling points with coherence greater than 0.5 are extracted to ensure the high quality co-seismic deformation data (Yu et al., 2020). Then, an adaptive quadtree method based on gradient variation is applied to downsample the co-seismic deformation data (Decriem et al., 2010; Bagnardi and Hooper, 2018). For each point, independent incidence and azimuth angles are used. This method retains the main characteristics of the original deformation field completely, and reduces the number of sampling points in the far field. Finally, 872 LOS displacement samples of the ascending orbit and 908 LOS displacement samples of the descending orbit are obtained (Supplementary Figure S2).

Spatial Characteristics of the Fault

Fault traces on the surface are extracted from the range and azimuth offset map. Amplitude pixel tracking method is used to obtain the range and azimuth offset of the Menyuan Ms 6.9 earthquake (Figure 4). The basic principle of this method is to calculate the pixel offsets in the range direction and azimuth direction according to the cross-correlation of amplitudes on the basis of sub-pixel registration of two SAR images (Hu et al., 2010). This method greatly compensates the insensitivity of InSAR to the north-south ground displacement, and helps to explain the incoherent phenomenon caused by the local large rupture near the fault. The azimuth offset results of T128 ascending and T33 descending orbits are opposite on both sides of the fault, showing a sharp color discontinuity (Figures 4B,D). Utilizing the offset maps with the boundary between positive and negative offset values, the surface traces of seismogenic fault can be extracted by outlining the discontinuity (Zhang et al., 2011; Hua et al., 2021). The finally obtained surface trace is composed of two faults (Figures 2B,D), of which the average strike of the small south branch is 77°, and the main north fault is 108°, with a total length of 49.5 km.

FIGURE 4. The results of the Pixel-Tracking method for the ascending and descending orbit amplitudes. (A,C) are the range offset of the T128 ascending orbit and the T33 descending orbit, respectively, and the positive and negative values indicate the decreasing and increasing of the horizontal surface displacement along the distance direction from the satellite and to the target point (the range direction); (B,D) are the azimuth offset of the T128 ascending orbit and the T33 descending orbit, respectively, and the positive and negative values represent the forward and backward displacement of the horizontal displacement along the flight direction of the satellite, respectively. The red star indicates the epicenter. Black arrows indicate satellite flight direction (AZI) and right-looking observation direction (LOS).

In the addition, the four profiles of the co-seismic deformation fields of the ascending and descending orbits (Figure 3) show that the displacements on both sides of the cross fault are approximately complementary. The existing focal mechanism solutions show that the dip angle of the seismogenic fault is close to vertical (Table 2), so the range of dip angle is set to 70°–90°. In order to select the optimal dip angle, a grid search with step of 1° was conducted in the inversion (Li et al., 2018), and the optimal dip angle was obtained as 79° (Figure 5A).

FIGURE 5. Determination of the optimal fault dip angle and the smoothing factor. (A) The misfit-dip curve. The light blue star indicates the optimal dip angle of 79°. (B) The model roughness and the fitting residual curve. The star indicates the best smoothing factor of 0.05.

Distributed Slip Model

According to the uniform isotropic elastic half-space dislocation model, the slip on the seismogenic fault plane can be linearly projected to the earth’s surface by Green’s function (Okada, 1985). The relationship between the slip distribution of the fault and the co-seismic displacement field is as follows:

where $u$ represents the co-seismic displacement measured by InSAR; G(m) represents Green’s function related to fault parameters m; $\mathit{s}$ represents the slip vector on the fault plane; $\epsilon$ represents the data observation error.

Generally speaking, the inversion of fault parameters includes the location, length, width, depth, strike and dip angle of the fault, and the slip vector on the fault, which is a highly nonlinear problem. Usually “two-step” inversion is adopted. Firstly, geometric parameters of seismogenic faults are resolved using a nonlinear inversion method (Bagnardi and Hooper, 2018). Then, according to the linear relationship between $u$ and $\mathit{s}$, the least square solution is obtained (Feng and Li, 2010). Especially, the Menyuan Ms 6.9 earthquake produced obvious surface ruptures. So the geometric shape along strike of the fault can be determined according to the surface rupture (Fukahata and Wright, 2008). Further, the dip angle of faults can be estimated according to the aftershock distribution, focal mechanisms and tectonic settings. Therefore, the slip distribution of faults can be directly solved by a linear inversion method.

In this paper, the Steepest Descent Method (SDM) inversion method (Wang et al., 2013) was used to solve the rake and slip size of each discrete sub-fault in the distribution slip model. This method has been widely used in GPS and InSAR co-seismic displacement inversion (Chen et al., 2018; Shi et al., 2019; Yu et al., 2020; He et al., 2022). SDM is an iterative optimization algorithm with better convergence and less dependence on the initial value. In order to minimize the objective function in each calculation, a step size proportional to the negative gradient is used in each iteration.

Meanwhile, in order to ensure smoothness and continuity between each sub-faults, a smoothing factor is added in the inversion process, and the objective function is as follows:

where ${\alpha }^2$ denotes the smoothing factor used to control the roughness of the slip distribution; $H$ denotes the finite-difference approximation of the Laplace operator; $\tau$ denotes the shear stress drop.

In the process of inversion, Poisson’s ratio is set to 0.25, the top depth of the fault is 0 km, and the dip angle is 79° (Figure 5A). Referring to the focal mechanisms of the mainshock (USGS, GCMT, and GFZ), the slip angle is empirically constrained to -15°–25° to accelerate the convergence. According to the relocation results of aftershocks, the aftershock concentrates at 5–10 km depth, and a few of aftershocks occurred between 11 and 20 km depth (Xu et al., 2022; Yang et al., 2022). Thus, the width of the fault was set to 20 km corresponding to the depth of 19.6 km. Then, the whole fault plane is divided into a 2 × 2 km rectangle, with 250 sub-faults in total. The maximum slip is set at 10 m. At the same time, the weight of ascending and descending orbit deformation was set to 0.8:1, taking into account the impact of afterslip in the ascending displacement. Independent incidence and azimuth angles are used for each sampling point (Zhang et al., 2011).

In order to further ensure the stability and reduce the uncertainty of the inversion, a grid search method was adopted to obtain the L-curve between model roughness and fitting residual for a series of smoothing factor (Figure 5B). Consequently, the optimal smoothing factor 0.05 was selected.

Figure 6 shows the fault slip distribution of the Menyuan Ms 6.9 earthquake. On the whole, the joint inversion results of ascending and descending orbits are similar to those of Yang et al. (2022), but different from those of Li et al. (2022a). This may be related to the difference in the geometry of faults. It can be seen from the slip distribution results that the slip is mainly concentrated in the strike range of 12–30 km, and 0–10 km along the dip. The maximum slip is 3.12 m at a depth of about 5 km, and the average slip of the main fault is 0.48 m (Supplementary Figures S3–S5).

FIGURE 6. Fault slip distribution of the Menyuan Ms 6.9 earthquake. (A,B) are the slip distributions of the main fault and the minor branch fault along the fault strike, respectively; black arrows represent the slip vector, and its direction represents the movement direction of the hanging wall relative to the foot wall; the blue star indicates the location of the mainshock determined by Yang et al. (2022).

Slip directions of fault patches show that the rupture is mainly sinistral strike-slip. In the slip center on the main fault, a small amount of thrust slip component exists.The inversion results show that the seismic moment is 1.28 × 10¹⁹ N·m, and the equivalent moment magnitude is 6.67, consistent with the GCMT result.

In the joint inversion results of ascending and descending orbits, the correlation coefficient between the observed and the simulated deformation is 98.87%. The root mean square error (RMSE) ascending descending orbits are 5.5 and 5.6 cm, respectively.

As can be seen from the spatial distribution of residuals (Figures 7C,F), the fitting residuals are small in the whole. The largest residuals appear in the area near the surface rupture zone. This may be caused by displacement errors from strong incoherence in the zone and the complexity of the very shallow rupture.

FIGURE 7. Fitting results in distributed slip model inversion. (A–C) represent the observed deformation, simulated deformation and the fitting residual of the T128 ascending orbit, respectively; (D–F) are same for the T33 descending orbit. The background image is DEM.

Figure 8 shows the histogram of fitting residuals. The distribution of ascending fitting residual is more concentrated (Figure 8A). Except for a few large residuals in the descending orbit, most of the residuals are close to zero, showing approximate Gaussian distribution. The ascending average fitting residual error is 1.52 cm with the variance of 2.8 mm, and the descending average fitting residual error is −0.93 cm with the variance of 3.1 mm.

FIGURE 8. Histograms of the fitting residuals. (A,B) represent the fitting residual histograms of the inversion results of the T128 ascending orbit and the T33 descending orbit. The red solid line indicates the Gaussian distribution curve fitted by the distribution.

Discussion

Regional Seismogenic Environment and Source Characteristics

From the perspective of historical earthquakes, there were two Ms 6.4 earthquakes in 1986 and 2016 around the Lenglongling fault (Figure 1). Coulomb stress analysis results show that the 2016 Menyuan Ms 6.9 earthquake has a loading effect on the 2022 Menyuan Ms 6.9 earthquake (Li et al., 2022a; Peng et al., 2022). Judging from the trace of shallow fault, Menyuan Ms 6.9 earthquake occurred at the intersection of Tolaishan fault and Lenglongling fault (Figure 1B), which belongs to the compression-bending zone on the Qilian-Haiyuan fault zone, and is prone to stress concentration (Pan et al., 2022). Geological survey shows that although the Tuolaishan fault and Lenglongling fault are mainly left-lateral strike-slip faults, their slip rates are different: the eastern segment of the former is about 3–4 mm/a, while the latter is about 6.4 ± 0.7 mm/a (Guo et al., 2017). P-wave velocity structure shows that there is an obvious low-velocity zone below the earthquake source (Wang et al., 2022; Xu et al., 2022), and the eastern part of the Qilian Mountains is located at the region where crustal thickness changes dramatically. This makes the source region favorable to stress concentration (Wang et al., 2022).

Based on the co-seismic deformation field of the ascending and descending orbits obtained by InSAR, the fault slip distribution is inverted. Its focal mechanism shows that the Ms 6.9 earthquake was a sinistral strike-slip event with a small amount of thrust component. The average strike of the seismogenic faults was 108°, and the dip angle was 79°. The average rake is 4° and there is only one rupture center (Figure 6). Combining the spatial position of the maximum slip (Figure 6A) and the deformation coverage shown in the co-seismic deformation field (Figure 2D) and azimuth offset (Figure 4D) of T33 descending orbit, we can conclude that the seismogenic fault of this earthquake is the west part of the Lenglongling fault. It is belonged to a left-stepped echelon fault system, forming a tensional zone with the adjacent Tolaishan fault. Stress transfer between the faults is likely to occur, which may be the cause of surface rupture at the eastern end of Tolaishan fault (Pan et al., 2022).

To sum up, under the continuous northward expansion of the northeastern margin of the Tibetan Plateau, the shallow faults in the area around Menyuan are characterized by compressing and bending with inconsistent slip rates. In the deep crust, great differences of physical properties exist. In addition, the historical earthquake had a stress loading effect on the rupture. These conditions may jointly promote the strain accumulation in the western segment of Lenglongling fault and ultimately lead to the occurrence of this earthquake.

The Relationship Among Aftershock, Surface Rupture Traces and Co-seismic Slip Distribution

Figure 9 is the surface rupture trace delineated according to the obvious amplitude azimuth offset and the boundary between positive and negative displacement in the co-seismic deformation field (Figures 2B,D). Its overall shape is similar to the surface rupture distribution map shown by Peng et al. (2022). Coupled with the relocation of aftershocks, it is revealed that this earthquake produced obvious surface rupture. According to the aftershock distribution in Figures 9A,B the aftershocks in 2–9 days after the earthquake continued to expand in the SE direction along the surface rupture trace and are concentrated in the eastern section of the fault. Aftershocks are evenly distributed on both sides of the fault in the eastern segment. However, aftershocks are mainly distributed on the south side of the in the middle and western segment, indicating a possible deflection of the fault.

FIGURE 9. Surface rupture traces and aftershock distribution. (A) The background image is the azimuth offset of the ascending orbit, and the small white circles represent the aftershocks in 9 days after the mainshock (Yang et al., 2022); (B) The background image is the azimuthal offset of the descending orbit, the white small circles represent the aftershocks of 2 days after the earthquake (Yang et al., 2022). Green lines denote the seismogenic faults; the yellow star represents the epicenter of the mainshock; pink circles represent two aftershocks with magnitude≥ 5.0, which occurred on January 8 and 12 from left to right, respectively.

At the same time, there is an obvious varying dipping angles from east to west in the aftershock profiles (Figure 10). The fault plane is determined by surface rupture and distribution of aftershocks. We tested a new fault geometry that each section has a different dip angle (Supplementary Figure S6). However, the results are similar to those obtained with uniform 79° dip (Supplementary Figures S7, S8) and the correlation coefficient between the observed and the simulated deformation is also over 98%. So 79° may be the average of these dips.

FIGURE 10. Dip angles determined from aftershock profiles of different fault sections. (A) Green lines, the yellow star and pink circles denote the seismogenic faults, the epicenter of the mainshock and two aftershocks with magnitude ≥ 5.0, respectively. (B–F) represent the five aftershocks profiles (aa’, bb’, cc’, dd’, ee’). Fault positions are marked by the green rectangles.

The aftershock (purple solid circles) is projected onto the seismogenic fault surface (Figure 11). The results show that the aftershock extended along the east and west sides from the source location within the depth range of 2–12 km. Meanwhile, the main shock depth (7.05 km) was close to the depth (5 km) where the maximum slip of the fault was 3.12 m. In addition, most of aftershocks occurred in the position where the co-seismic slip was small.

FIGURE 11. Fault slip model of the main fault and aftershocks. Red lines represent the contour of the slip; the blue star represents the mainshock location determined by Yang et al. (2022); the cyan circles denote two aftershocks with magnitude ≥ 5.0.

Therefore, the stress release associated with the aftershocks was either triggered by the mainshock or it was a compensation to the stress release associated with the mainshock.

Regional Seismic Hazard Assessment

In general, after an earthquake of moderate magnitude or above, the co-seismic dislocation of the fault will cause static co-seismic permanent deformation, which will change the Coulomb stress in the near and far field. the Coulomb failure stress change $\left(\text{Δ}CFS\right)$ can be used to analyze the triggering of aftershock, and the interaction between earthquakes and faults (King et al., 1994; Xu et al., 2018). According to Coulomb instability criterion, the Coulomb stress change can be expressed as:

where $\text{Δ}{\tau }_s$ represents the change in shear stress; ${\mu }^{\prime }$ represents the effective coefficient of friction on the fault; $\text{Δ}{\sigma }_n$ represents the change in normal stress (dilation for positive); $\text{Δ}CFS$ represents the coulomb stress change on the receiver faults caused by the slip of the source fault, with a positive value indicating stress loading and a negative value indicating inhibition of rupture of the receiver fault. In this paper, Coulomb 3.3 software (Toda et al., 2011) was used to calculate the co-seismic Coulomb stress changes caused by Menyuan Ms 6.9 earthquake (Figure 12).

FIGURE 12. Regional co-seismic static Coulomb stress changes at different depths. (A–D) represent the changes of static Coulomb stress at depths of 0, 5, 15, and 10 km, respectively. White circles represent aftershocks in 10 days after the earthquake, including two M ≥ 5 aftershocks (pink circles) (Yang et al., 2022), which are located at the specified depth within 2 km up and down; the gray solid lines indicate the fault traces in this area.

The receiver fault was set as seismogenic fault itself (strike = 108°, dip angle = 79°, dip angle =4°) with a friction coefficient of 0.4 (King et al., 1994). The co-seismic Coulomb stress changes at the depths of 0, 5, 10 and 15 km are calculated respectively (Figure 12). The results show that the coulomb stress loading region is mainly distributed in the eastern segment of Tuolaishan fault and the western segment of Lenglongling fault, and a few parts of the Minyue-Damaying fault, Minyue-Yongchang fault and Daban Mountain fault in the NS direction. The static Coulomb stress change at 10 km shows that the largest aftershock (Ms 5.2) occurred in the positive area of the Lenglongling fault 4 days after the mainshock. The aftershock mainly migrated towards to SE direction (Fan et al., 2022; Yang et al., 2022), which may be caused by stress loading yielded by the mainshock.

There are many faults in this area, and the stress disturbance at the depth of 10 km on six faults is also calculated (Figure 13). The results show that some parts of the Lenglongling fault, the Tolaishan fault, the Sunan-Qilian fault, and the Minyue-Damaying fault are loaded, indicating high seismic risk in the future. At the same time, the whole Qilian-Haiyuan fault (Figure 1A), the Jinqianghe-Maomaoshan-Laohushan fault is located in the “Tianzhu Gap” between the present earthquake and the Haiyuan earthquake in 1920 (Gaudemer et al., 1995; Li et al., 2017), so the risk of these faults is also worth noting.

FIGURE 13. Static Coulomb stress changes on the surrounding faults at depth of 10 km. The colored lines represent the changes of the static Coulomb stress on the surrounding faults caused by the dislocation of the seismogenic fault, and the positive represent the stress loading. Fault parameters are from Xu et al. (2016). The calculation parameters are: Menyuan Fault (reverse fault, dip 45°, rake 90°); Lenglongling Fault (left-lateral, dip 90°, rake 0°); Tolaishan Fault (left-lateral, dip 90°, rake 0°); Sunan-Qilian Fault (reverse fault, dip 45°, rake 90°); Minyue-Damaying fault (reverse fault, dip 45°, rake 90°); Huangcheng-Shuangta Fault (reverse fault, dip 45°, rake 90°). Grey lines indicate the fault traces in this area; black line indicates the “Tianzhu Gap”, including the LLLF and the JQHF, MMSF, LHSF to its east. The background image is DEM.

Conclusion

In this paper, the co-seismic deformation field (LOS direction) of Menyuan Ms 6.9 earthquake in 2022 is obtained by D-InSAR technology. The co-seismic slip distribution of the fault is inverted by the Steepest Descent Method, and the disturbance of Coulomb stress changes on the surrounding faults is calculated and analyzed. This provides a reference for further research on the seismogenic property and seismic risk of the Qilian-Haiyuan fault. The main conclusions are as follows:

1) The seismogenic fault of the Menyuan Ms 6.9 earthquake is the western segment of Lenglongling fault, and surface deformation filed caused by the co-earthquake is 50 km by 40 km approximately. According to the distribution of ground deformation in the LOS direction and azimuth offsets of the ascending and descending orbits. The Menyuan Ms 6.9 earthquake can be identified as a sinistral strike-slip event.

2) The co-seismic slip distribution shows that the main seismogenic fault strike 108°, dip 79° and rake 4°. The slip is concentrated at a depth of 0–10 km, with an average slip of 0.48 m and the maximum slip of 3.12 m, which is located at a depth of about 5 km. The seismic moment is 1.28 × 10¹⁹ N·m, and the corresponding moment magnitude is Mw 6.67.

3) The results of Static Coulomb stress changes show that some parts of Lenglongling fault, Tolaishan fault, Sunan-Qilian fault, and Minyue-Damaying fault are loaded by the Menyuan Ms 6.9 earthquake. Risk of strong earthquake on these faults in the future deserves attention.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://search.asf.alaska.edu/#/.

Author Contributions

JL processed the SAR data and wrote the manuscript; CL contributed to the conception of the research and result analysis; CW contributed to invert the slip distribution and calculate the Coulomb stress changes; FC contributed the related code and beautified the figures; CY and YY helped the result analysis with constructive discussions.

Funding

This research was funded by the National Natural Science Foundation of China (Grants 4217040570, 41674059).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Sentinel-1A SAR data are provided by the European Space Agency. The historical earthquake catalogue is provided by Zhang Zhiwei, Sichuan Earthquake Administration. SDM inversion program is provided by professor Wang Rongjiang. Aftershocks data is provided by Yang Hongfeng, The Chinese University of Hong Kong. GMT software is used to create most figures in the paper.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.948661/full#supplementary-material

References