High-rate GNSS-based rapid determination of coseismic deformation and source characteristics for the 2023 M6.2 Jishishan earthquake

1.   Introduction

According to the China earthquake Network Center, an M6.2 earthquake struck Jishishan Bonan, Dongxiang, and Salar Autonomous County (Jishishan County), Linxia Hui Autonomous Prefecture, Gansu (35.70°N, 102.79°E) at 23:59:30 (UTC+8) on December 18, 2023, at a focal depth of 10 km. The epicenter was only 8 km from the urban area of Jishishan County. The earthquake generated violent ground motions within a 100 km radius of the epicenter, with a maximum intensity of VIII. A significant number of houses were damaged, some even collapsed, and overall infrastructure sustained substantial damage. Additionally, the earthquake triggered secondary geological disasters such as landslides, slope failures, soil liquefaction, and several localized mudflows (Huang GW et al., 2023; Chen B et al., 2024; Xu Q et al., 2024). As of 08:00 on December 22, 2023, the earthquake had resulted in 117 deaths and nearly a thousand injuries (https://news.cctv.com/2023/12/22/ARTI5QGeiSwnxCdb1UGmzabn231222.shtml). It was the deadliest earthquake in China in 2023.

The epicenter was located in the eastern part of the Qaidam-Qilian Block, near the Lajishan-Jishishan Fault. The Lajishan Fault Zone comprises two northeastward protruding arc-shaped tectonic belts: the northern and southern margin faults. This fault zone serves as a transition zone, controlling both the right-lateral strike-slip movement of the Riyueshan Fault in the northwest and the left-lateral strike-slip movement of the northern margin fault of the West Qinling in the southeast (Yuan DY et al., 2013). The Lajishan Fault is characterized by compressive tectonics, with local left-lateral strike-slip movement. Historically, there have been more than 20 moderate earthquakes of magnitude ~5 (Yuan DY et al., 2005). Field investigations and trench excavations have revealed at least two moderate-to-strong or severe paleo-earthquakes since 3700 BP, events that may be closely associated with the formation of the Lajia Site (Li ZM et al., 2009). The northwest-striking Jishishan Fault, located south of the Lajishan Fault, primarily exhibits thrusting (Lease et al., 2011). Current geodetic research suggests that the Jishishan Fault demonstrates significant thrusting as well as right-lateral slip movement, indicating the possibility of moderate-to-strong earthquakes (Zhuang WQ et al., 2023). Moment tensor solutions, based on the inversion of broadband seismic waveform data from China and the Global Seismographic Network (centroid moment tensor (CMT) products for large earthquakes from the China earthquake Networks Center, reports from the Institute of Geophysics, China earthquake Administration, and W-phase MT as well as global (G)CMT from the United States Geological Survey (USGS)), all indicate significant thrusting. Specifically, the strike/dip/rake of nodal planes 1 and 2 are 156°–170°/28°–57°/93°–123° and 300°–333°/45°–62°/50°–88°, respectively. The macroscopic tectonic characteristics of the region and the motion characteristics of the Lajishan-Jishishan Fault suggest a west-dipping thrust seismogenic mechanism. However, the accurately plotted aftershock distribution illustrates that the fault is east-dipping (Figure 1) (Wang SG et al., 2024; Zuo KZ and Zhao CP, 2024). This indicates complex source characteristics for the Jishishan earthquake.

Figure  1.  Tectonic background and aftershock distribution. The red star denotes the epicenter. (a) Location of the study area. (b) Regional tectonic background. The beach ball illustrates the focal mechanism (adopted from the reports released by the Institute of Geophysics, China earthquake Administration). RYSF: Riyueshan Fault, DTH-LXF: Daotanghe-Linxia Fault, LJSSF: Lajishan Southern Margin Fault, LJSNF: Lajishan Northern Margin Fault, WQLNF: West Qinling Northern Margin Fault. (c) Distribution of aftershocks and high-rate GNSS stations (aftershock data adopted from Zhang L et al., 2024; Wang AJ and Gao Y, 2024). The arrows represent the surface displacements caused by the earthquake. JSSWF: Jishishan Western Margin Fault, JSSEF: Jishishan Eastern Margin Fault. AB denotes the cross-section of aftershocks (d), and the dashed line is the fitted fault plane based on the spatial distribution of aftershocks.

DownLoad: Full-Size Img PowerPoint

The high-rate global navigational satellite system (GNSS) directly records displacements based on the absolute position reference. It is advantageous because there is no range saturation and no accumulation of measurement error over time. Thus, it has been widely applied in static dislocation extraction, tsunami early warning systems, real-time moment magnitude determination, near-field hypocenter information acquisition, and quasi-real-time inversion of fault rupture processes (Bock et al., 2000; Crowell et al., 2009; Colosimo et al., 2011; Fang RX et al., 2013; Li XX et al., 2013, 2019; Li XX, 2017; Geng T et al., 2016; Shan XJ et al., 2019). Li ZC et al. (2024) utilized high-rate GNSS to perform post-event data resolution for the Jishishan earthquake and obtained distinct coseismic deformation waveforms. In this study, high-rate data from near-field continuous GNSS stations were utilized to simulate data processing under real-time conditions. This approach enabled the rapid determination of coseismic static and dynamic deformation caused by the event and the estimation of the empirical magnitude. High-rate GNSS data and far-field body waves were combined to constrain the source rupture process and analyze potential seismogenic fault structures. The ultimate goal of the study was to provide references for further identification of seismogenic faults and their tectonic characteristics.

2.   Data and methods

2.1   GNSS and teleseismic data

High-rate (1 Hz) GNSS data were collected from four continuously operating reference stations (CORS) in Gansu and Qinghai within 100 km of the epicenter and the Delingha Station of the China Continental Tectonic Environment Monitoring Network (Figure 1). Three approaches, namely dynamic differential positioning, precise point positioning with ambiguity resolution (PPP-AR), and single-point velocity determination, were applied to the high-rate GNSS data. (1) Quasi-real-time dynamic deformation (RTDD) was determined through differential positioning. The Delingha Station, which is relatively far from the epicenter, was selected as the reference station to determine the double-difference models with the four stations (LXJS, GUTI, XUNH, and DUOW) (Figure 1b and 1c) and eliminate receiver and satellite clock errors. With the help of the IGS Ultra-rapid products (IGU) available in real-time, single-epoch differential dynamic resolution was performed using GAMIT/TRACK to obtain the relative displacements of the stations concerning the reference station (Fan SJ et al., 2013). (2) Precise point positioning (PPP) was adopted to accurately estimate the dynamic displacements after the event. Products released by Wuhan University, including rapid satellite ephemerides at 5-min intervals and clock error products at 30-second intervals (with a delay of approximately 15–20 h) were adopted to obtain float solutions by precisely modeling all errors of single-station GNSS data. Next, code/phase bias products were utilized to remove the ambiguity and acquire the single-epoch absolute displacements during the earthquake (with the help of PRIDE PPP-AR) (Geng JH et al., 2019). (3) The single-point velocity was determined rapidly based on the difference between epochs (RTV). The GNSS broadcast ephemeris available in real-time was used in combination with the variometric approach (Colosimo et al., 2011; Grapenthin et al., 2018; Fang RX et al., 2021; Zang JF et al., 2022), which is based on differences between GNSS epochs, to rapidly estimate the velocity series during the earthquake. In this study, the displacements obtained via PPP-AR based on high-precision products and the velocity results calculated using the displacement differences were regarded as the true values. They were compared to the displacements obtained through dynamic differential positioning in real-time and the single-point velocity determination results to verify the feasibility of rapid or real-time resolution GNSS data after an earthquake occurrence.

At epicentral distances ranging from 30° to 90°, far-field broadband seismic waves propagate in the relatively homogeneous mantle, with distinct and clean seismic phases. They are insignificantly affected by the lateral heterogeneity of the medium. Direct P-wave data with high signal-to-noise ratios were adopted from 41 stations of the Global Seismographic Network. More specifically, 30-second data were extracted starting 10 before the direct P-wave arrival. The extracted data were then band-pass filtered at 0.01–0.5 Hz and used to constrain the source model after resampling.

2.2   Moment magnitude estimation based on GNSS dynamic waveforms

Typical ground motion indices in the time domain include peak ground acceleration (PGA), effective peak ground acceleration (EPGA), peak ground velocity (PGV), and peak ground displacement (PGD). Here, the PGD and PGV recorded by the high-rate GNSS (1 Hz) during the strong earthquake were adopted as the key parameters for magnitude estimation. The formulas for extracting PGD from three-component displacement waveforms recorded by high-rate GNSS during the earthquake are as follows:

where $d_{\mathrm{N}}(t)$/$v_{\mathrm{N}}(t)$, $d_{\mathrm{E}}(t)$/$v_{\mathrm{E}}(t)$, and $d_{\mathrm{U}}(t)$/$v_{\mathrm{U}}(t)$ are the GNSS displacement/velocity records during the earthquake in the north-south, east-west, and vertical directions, respectively. The noise in vertical displacements recorded by GNSS is three to five orders of magnitude higher than that of the horizontal ones (Melgar et al., 2015). Furthermore, no noticeable seismic signals were observed for the Jishishan earthquake in the vertical displacement records (Figure 2). A remarkable vertical velocity increase was only noted for the GUTI station (Figure 3). To avoid errors introduced by the vertical components, the horizontal peak displacement and velocity were used for magnitude estimation. The peak ground displacement scaling law proposed by Crowell et al. (2013) for high-rate GNSS (1 Hz) data records was applied. The scaling law considers ground motion attenuation associated with the magnitude to obtain the relative strength of the near-, intermediate-, and far-field ground motion radiation, as shown below:

Figure  2.  High-rate GNSS displacement waveforms. PPP (AR) results in red; RTDD (Track) results in black; 0 denotes the origin time

DownLoad: Full-Size Img PowerPoint

Figure  3.  High-rate GNSS velocity waveforms. Velocities based on displacement differences (Diff vel) in red; SPV (Vel) results in black; 0 denotes the origin time

DownLoad: Full-Size Img PowerPoint

where A, B, and C are regression coefficients; M_W is the moment magnitude; and R denotes the focal distances. The epicentral location and depth were adopted from the values measured by the China earthquake Network Center. R is the focal distance instead of the epicentral distance. This is particularly important to magnitude estimation for stations with epicentral distances similar to or smaller than the focal depth. In reported studies, multiple attempts have been made to determine and update regression coefficients A, B, and C based on high-rate GNSS (1 Hz) displacement and velocity records of actual earthquakes (Crowell et al., 2013, 2016; Melgar et al., 2015; Ruhl et al., 2019; Fang RX et al., 2021). The specific values of the coefficients are listed in Table 1. In this study, the regression coefficients reported by Crowell et al. (2016) were adopted for PGD: A = −6.687, B = 1.047, and C = −0.138. For PGV, the values proposed by Fang RX et al. (2021) based on 22 moderate-to-strong earthquakes were used: A = −5.025, B = 0.741, and C = −0.111.

Table  1.  Regression coefficients and standard deviations of empirical moment magnitude equations

ID	Parameter	A	B	C	Standard deviation	Unit	Source
1	PGD	−5.013±0.211	1.219±0.046	−0.178±0.01	0.224	cm	Crowell et al., 2013
2	PGD	−4.434±0.141	1.047±0.022	−0.138±0.003	0.270	cm	Melgar et al., 2015
3	PGD	−6.687	1.500	−0.214	0.170	cm	Crowell et al., 2016
4	PGD	−5.919	1.009	−0.145	0.210	m	Ruhl et al., 2019
5	PGV	−5.025±0.084	0.741±0.017	−0.111±0.003	0.389	m/s	Fang RX et al., 2021

 | Show Table

DownLoad: CSV

2.3   Slip distribution inversion

Multiple-window linear inversion was used to determine the source rupture process. The fault plane was divided into several sub-faults along its strike and dip directions. Each sub-fault was approximated as a point source, and its source time function was approximated by the linear superposition of several isosceles triangular functions with fixed durations. In this way, the linear relation between the displacement of the sub-fault and the observed data could be established (Liu G et al., 2017; Li Q et al., 2022). To resolve the displacement of the sub-fault, both the fitting degree of the theoretical calculation results with the observations and the roughness of the slip distribution should be considered:

where M denotes the total number of sub-faults; K is the number of time windows for each sub-fault; x is the position vector of station n; G presents the Green’s function matrix; and g_nj is the displacement response at station n to a unit slip (1 m) of the jth sub-fault. For far-field body waves, dynamic Green’s functions were calculated based on spherical layered media, which could be completed with the help of QSSP (Wang RJ et al., 2017). For near-field GNSS, the corresponding static Green’s functions were calculated based on planar layer media using PSGRN/PSCMP (Wang RJ et al., 2006). O denotes the vector of the observed data; u_n(t) is the observed value at station n; m represents the slip vector to be resolved; m_ij is the slip of the jth sub-fault in ith time window; G·m is the simulated value; s(t) denotes the source time function; T₀ is the origin time; D_j denotes the rupture propagation distance, i.e. the distance from the initial rupture point to the sub-fault; v is the mean rupture velocity; d represents the time width of the isosceles triangle; T_ij is the rupture time of the jth sub-fault in ith time window; 1 and 2 indicate two slip directions at angles typically equal to rake±45° to realize non-uniform rake angles of sub-faults; L is the Laplace second-order differential operator; β denotes the smoothing factor, the optimum of which can be obtained from the difference curve of the model roughness || L·m||² with the fitting residuals of observed data ||G·m−O||²; and ω is the weight of the observations. The same weight, 1, has been assigned to GNSS and far-field seismic waves.

3.   Results

3.1   High-rate GNSS coseismic displacements and moment magnitudes

On the day of the earthquake, the displacement and velocity time series at GNSS stations from 15:00–17:00 (UPS time) were acquired. After the original time was aligned with UTC (+18 s), the coseismic displacement and velocity waveforms 10 s before and 120 s after the origin were plotted, as illustrated in Figures 2 and 3. The Jishishan (LXJS) and Guanting (GUTI) stations showed significant coseismic responses. At LXJS, the horizontal PGD reached approximately 3.7 cm 8 s after the origin, whereas the horizontal PGV was 4.9 cm/s 9 s after the origin. At GUTI, a horizontal PGD of 6.3 cm and a horizontal PGV of 6.1 cm/s were observed at 12 and 13 s after the earthquake, respectively. Compared to LXJS, which is situated south of the epicenter and has a shorter epicentral distance, GUTI, located 23 km north of the epicenter, has captured more significant transient deformation amplitudes. This implies predominantly northward rupture propagation. Strong ground motions propagated to the Xunhua (XUNH) station, located approximately 30 km from the epicenter, with a horizontal PGD of 2.3 cm and horizontal PGV of 1.6 cm/s. Subsequently, ground motions weakened and became undetectable in the noise at the relatively distinct Duowa (DUOW) station, with an epicentral distance of 74 km, as the displacement amplitude decreased to below 1 cm. No significant coseismic signals were detected in the vertical displacements resolved by the quasi-RTDD and PPP approaches at all stations. Furthermore, the vertical velocity at GUTI obtained by the SPV method increased noticeably 10 s after the origin, and the peak reached 2.63 cm/s.

The mean epoch coordinates under the pre- and post-earthquake steady states in the displacement time series were obtained. Next, based on the coordinate differences, the permanent coseismic surface displacement caused by this earthquake was calculated. It was approximately 1.9 cm at LXJS (−1.3 and 1.4 cm for the north-south and east-west directions, respectively) and approximately 0.5 cm at GUTI (−0.4 and 0.1 cm for the north-south and east-west directions, respectively). However, no permanent deformation was noted at XUNH and DUOW. These findings agree with the principle that permanent displacement rapidly attenuates with distance. The differences between displacement waveforms acquired using the quasi-RTDD and post-earthquake PPP resolution approaches were calculated. At GUTI, the root mean square (RMS) values of the east-west and north-south differences were 0.44 and 0.75 cm, respectively (Figure 4a). The differential velocities estimated based on the SPV and PPP methods were also compared. The RMS values of the differences in the east-west and north-south velocities were 0.22 and 0.31 cm/s, respectively (Figure 4b). These RMS values are lower than the error level of the single-epoch resolution strategy (Li XX et al., 2019).

Figure  4.  Differences between estimated displacements and velocities during the earthquake. RMS in east-west, north-south, and vertical directions in red, blue, and black, respectively

DownLoad: Full-Size Img PowerPoint

With the help of the high-rate GNSS magnitude scaling equation, the magnitude was estimated in real-time using the horizontal PGD and PGV. At LXJS, the PGD-based magnitude (M_PGD) peaked and reached 5.71 at 8 s after the origin time, whereas the maximum PGV-calculated magnitude (M_PGV) of 5.96 was attained at 9 s. Similarly, at GUTI, the maximum M_PGD and M_PGV of 6.24 and 6.51 were observed 12 and 13 s after the earthquake, respectively. At XUNH, the maximum M_PGV and M_PGD of 5.69 and 6.03 were observed only 14 and 19 s after the origin time, respectively. The mean moment magnitude determined at the three stations (M_PGD 5.99 and M_PGV 6.05) is consistent with the value (M_W5.9) released by the USGS (Table 2). The calculated M_PGD at LXJS, the closest station to the epicenter, is significantly underestimated.

Table  2.  PGD, PGV, and estimated magnitudes of the stations

ID	Station	Focal distance (km)	PGD (cm)	PGV (cm/s)	MPGD3	MPGD3 (Detrend)	MPGV	MPGD1	MPGD2	MPGD4
	Standard deviation	-	-	-	±0.17	±0.17	±0.389	±0.224	±0.27	±0.21
1	LXJS	11.73	3.68	4.85	5.71	5.7	5.96	5.42	5.56	5.25
2	GUTI	25.34	6.3	6.08	6.24	6.03	6.51	6.0	6.13	5.86
3	XUNH	36.06	2.26	1.62	6.03	5.99	5.69	5.7	5.75	5.53
4	MEAN	-	-	-	5.99	5.91	6.05	5.71	5.81	5.55
Note: The number after MPGD corresponds to the ID in Table 1, and MPGD3 is the result based on the scaling law adopted in this study.

 | Show Table

DownLoad: CSV

3.2   Source rupture process of the Jishishan earthquake

Based on the hypocenter location (102.8320°E, 35.7423°N, 15.1 km) and focal mechanism released by the Institute of Geophysics, China earthquake Administration, the Jishishan earthquake likely occurred on the southwest-dipping nodal plane with a strike/dip of 156°/44° or the northeast-dipping nodal plane with a strike/dip of 307°/50°. The rakes of the two nodal planes are 112° and 70°, respectively. Owing to multiple factors including the timing of data acquisition, it is often difficult to immediately determine which nodal plane corresponds to the actual seismogenic fault, particularly for blind earthquakes that produce no surface rupture. In this study, two scenarios were analyzed by considering the two nodal planes as the seismogenic fault separately. The source rupture process was constrained jointly by high-rate GNSS-based static and dynamic deformation results and far-field body waves. Based on the back-projected rupture dimension and duration, the rupture velocity was set to 1.6 km/s (Tang XW et al., 2024).

Based on its parameters, the southwest-dipping nodal plane was established with a length and depth of 52.5 km × 35.0 km. Along the strike and dip directions, the rectangular plane was divided into 21 × 14 = 294 sub-faults, each of which was 2.5 km × 2.5 km. The source-time function of each sub-fault was fitted using the linear superposition of three isosceles triangle functions, each with a duration of 2 s. The optimal model shows that the rupture propagates from the hypocenter to shallow regions along the dip direction. The maximum slip of approximately 0.44 m was located northwest of the epicenter, and the slip occurs in the shallow regions. The total seismic moment release was approximately 2.0 × 10¹⁸ N·m (M_W 6.1), and the energy was primarily released in the first 5 s. The rupture was predominated by thrusting, with a left-lateral strike-slip component. Aftershocks were mainly located at greater depths where slip is observed.

The fitting of optimal simulation results illustrates the significant correlation between the observed and simulated far-field P-wave measurements. Similarly, the observed and simulated GNSS-based static deformation values are consistent in terms of magnitude and direction, whereas the simulated values are slightly larger. Furthermore, the observed and simulated GNSS dynamic waveforms are poorly coordinated. Poor-fitting is noted for waveform details, whereas baseband signals (permanent deformation) are well simulated. This is because the geometric parameters of the seismogenic fault defined based on the focal mechanism deviate from the actual location, strike, and dip angle of the fault. Near-field GNSS deformation waves are more sensitive to these parameters than far-field seismic waves.

Figure  5.  Source rupture of the southwest-dipping nodal plane. (a) Source slip distribution and measured and simulated GNSS values. (b) Source time function. (c) Locations of the far-field seismic stations. (d) Observed and simulated far-field body waves. (e) Observed and simulated high-rate GNSS values.

DownLoad: Full-Size Img PowerPoint

Similarly, based on the parameters of the northeast-dipping nodal plane, a finite fault was constructed. The fault plane division and source time function parameterization procedures are identical to those described for the southwest-dipping plane model. The results suggest that the rupture propagates from the hypocenter to shallow regions, whereas aftershocks are primarily distributed in shallow parts of the rupture area. The maximum slip was approximately 0.57 m, and a small amount of slip was observed in the shallow regions. The total seismic moment release was approximately 1.66 × 10¹⁸ N·m (M_W6.1), and the energy release primarily occurred in the first 5 s. The rupture was dominated by thrust, with a right-lateral strike-slip component. The observed and simulated far-field P-wave measurements demonstrate a significant correlation. The observed and simulated GNSS-based static deformation values are less consistent in terms of magnitude and direction than those for the southwest-dipping nodal plane. The fitting of the GNSS dynamic waveforms is similar to that of the southwest-dipping plane. Based on the fitting results, it is impossible to identify the seismogenic fault from the two nodal planes using high-GNSS. This is related to the sparse distribution of GNSS stations and the complexity of hypocenter geometric parameters.

Figure  6.  Source rupture of the northeast-dipping nodal plane. (a) Source slip distribution and measured and simulated GNSS values. (b) Source time function. (c) Observed and simulated far-field body waves. (d) Observed and simulated high-rate GNSS values.

DownLoad: Full-Size Img PowerPoint

4.   Discussion and conclusions

As the capability to process high-rate GNSS data in real-time has been enhanced, many studies have demonstrated through simulations that this technology can measure strong ground motions in real-time with centimeter-level accuracy. Thus, it is possible to determine first-order characteristics of the source rupture rapidly or in real-time (Fang RX et al., 2013; Geng T et al., 2016; Shan XJ et al., 2021; Gao ZY et al., 2021, 2022; Zang JF et al., 2022, 2024). In this study, high-rate GNSS data from the Jishishan earthquake were processed using three strategies: RTDD, SPV, and PPP. The consistency of the results proves the reliability and accuracy of the GNSS-based real-time determination of dynamic coseismic displacement (Li ZC et al., 2024). Nevertheless, whether this technique can be further utilized to reliably identify source characteristics is highly dependent on the spatial distribution of continuous GNSS networks.

In this study, PGD and PGV results from three stations and an empirical equation for calculating the moment magnitude were adopted to estimate the mean moment magnitude: M_PGD 5.99 and M_PGV 6.05. The values are reasonably consistent with the published moment magnitude (Table 2). The regression coefficient combinations of the empirical PGD equations in Table 1 were compared. The combinations proposed by Crowell et al. (2016) and Melgar et al. (2015) give results closest to the actual moment magnitude of the Jishishan earthquake, with differences of ±0.09. Conversely, the coefficient combination proposed by Ruhl et al. (2019) leads to the largest deviation (−0.35) from the actual value (Table 2). These results suggest that magnitude estimation results based on Equation (3) are highly dependent on the regression coefficients. According to Fang RX et al. (2021), the mean absolute deviations of M_PGD and M_PGV from M_W are 0.25 and 0.26, respectively. The main cause for the significant underestimation using the coefficients proposed by Ruhl et al. (2019) may be the sparse and uneven distribution of GNSS stations in the region and the insufficient sample space for magnitude simulation. Therefore, given uneven GNSS network distribution, deviations are observed when applying empirical equations for magnitude estimation. In addition, although the empirical scaling law (Equation (3)) is derived from global observations involving various tectonic regimes, it is a one-order approximation only (Melgar et al., 2015). As more earthquake events will be recorded by regional GNSS networks in the future, empirical equations and regression coefficients specific to particular regions or seismogenic mechanisms can potentially be formulated.

Regardless of the coefficient combination used, the M_PGD values estimated at LXJS, 11 km from the hypocenter, are smaller than the moment magnitude (Table 2). This has been reported in literature. Yao WM et al. (2019) utilized a surface-wave magnitude formula to estimate the magnitude of the 2016 M_W7.8 Kaikoura earthquake, New Zealand, using data from 41 high-rate GNSS stations in the area. The earthquake involved northwestward rupture, and the magnitude is estimated to be only M6.65 at a station perpendicular to the rupture direction located approximately 20 km from the epicenter. Zheng JW et al. (2023) used another surface-wave magnitude formula to determine the magnitude of the 2022 M_W6.6 Luding earthquake using GNSS data from seven stations. At a station approximately 48 km from the hypocenter, the calculated magnitude equals M6.49. Gao ZY et al. (2021) and Wang YB et al. (2024) utilized the coefficients proposed by Ruhl et al. (2019) and Melgar et al. (2015), respectively, to estimate the magnitude of the 2021 M_W7.4 Maduo earthquake. The corresponding results were 7.8 and 7.13 approximately 30 km from the hypocenter (Table 3). The empirical magnitudes calculated based on near-field waveforms are systematically underestimated because the geometric dimension of the hypocenter or the high heterogeneity of the rupture invalidates the point source assumption. The near-field waveforms and peak values observed at stations are controlled by local asperity ruptures, leading to magnitude underestimation. To verify this hypothesis, assuming the fault is east-dipping and using InSAR data as constraints, the source-slip distribution was obtained via inverse modeling and the three-component coseismic displacement waveforms via forward modeling at 40 stations uniformly distributed in azimuth 5, 10, 20, 30, and 40 km from the epicenter. The magnitude determination methods presented in this study were applied for analysis. Along the rupture direction, the estimated M_PGD values at the stations increase first but then decrease as the epicentral distance increases (Figure 7). The acquired mean magnitude is M_PGD5.53. The mean magnitude obtained via simulation is underestimated because only the low-frequency components of the simulated waveforms were used, instead of the complete frequency band. This is consistent with the observation that the estimated magnitudes are smaller than the actual values after low-pass filtering (considering the site amplification effect of the Yellow River alluvial plain where the study area is located) (Table 2: M_PGD3(Detrend)).

Table  3.  Magnitude estimation results based on near-field stations for different earthquakes

ID	earthquake	Date	Mechanism	MW	Equation	Estimated magnitude	Focal distance (km)	Source
1	Kaikoura earthquake, New Zealand	2016-11-13	strike-slip	7.8	$M={\mathrm{log}}A+1.66{\mathrm{log}}D+2.0$	6.65	20*	Yao WM et al., 2019
2	Luding earthquake	2022-09-05	strike-slip	6.6	${M}_{{\mathrm{S}}}={\mathrm{log}}\left(A/T\right)+1.66{\mathrm{log}}\left(\Delta \right)+3.5$	6.49	48	Zheng JW et al., 2023
3	Maduo earthquake	2021-05-21	strike-slip	7.4	Ruhl et al., 2019	6.8	30	Gao ZY et al., 2021
4	Maduo earthquake	2021-05-21	strike-slip	7.4	Melgar et al., 2015	7.13	35	Wang et al.,2024
Note: *The epicentral distance is adopted here.

 | Show Table

DownLoad: CSV

Figure  7.  Forward modeled magnitude results. The red star denotes the epicenter. The concentric circles, radiating outward from the epicenter, represent station locations at distances of 5, 10, 20, 30, and 40 km.

DownLoad: Full-Size Img PowerPoint

Given both the inaccurate empirical magnitude estimation and invalid point source assumption, finite fault inversion constrained by near-field GNSS can provide useful references for seismic characteristics. Near-field GNSS-based static and dynamic deformation observations and far-field P-wave measurements were combined to perform inversion modeling of the slip distributions on the two nodal planes. The common characteristics are asperity ruptures concentrated at the hypocenter. The rupture dimension is approximately 15 km × 15 km, which is consistent with the results constrained by InSAR and regional seismic records (Huang X et al., 2024; Liu ZJ et al., 2024; Yang JY et al., 2024). The maximum slip is larger than the inversion result solely based on seismic waves (Wang AJ and Gao Y, 2024); therefore, the corresponding moment magnitude is slightly larger. The rupture process demonstrates that upward and northward rupture propagation from the hypocenter is limited and lasts for approximately 10 s. Based on the data fitting results, although the near-field GNSS-based permanent deformation has a high degree of fitting for the southwest-dipping plane scenario, the fitting improvement is insignificant compared to that in the northeast-dipping plane scenario. Thus, it is impossible to distinguish the actual seismogenic faults from the two nodal planes. Nevertheless, the results still provide insights into the interrelationship between the rupture propagation process and aftershock distribution: Aftershocks propagate from the deep regions in the east to the shallow parts in the northwest. This agrees with the rupture propagation direction of the northeast-dipping plane from the hypocenter to the shallow regions and the northwest. However, the rupture direction of the southwest-dipping plane is the opposite. Given the horizontal accuracy of the precise location of aftershocks (Wang SG et al., 2024; Zuo KZ and Zhao CP, 2024), if the seismogenic fault is southwest-dipping (Tang XW et al., 2024), a large number of aftershocks should be concentrated at depths below the focal depth. Conversely, a northeast dip does not align with the distribution direction of the Lajishan and Jishishan Faults (Yuan DY et al., 2005, 2013). After considering the regional tectonic background, the seismogenic fault may be a local blind fault (Gao Y et al., 2024; Wang et al., 2024; Yang JY et al., 2024). Although the interrelationship between the rupture direction and aftershock distribution support a northeast-dipping fault, the northward distribution shift of aftershocks is poorly explained by the relatively uniformly distributed asperities concentrated at the hypocenter, indicating the complexity of the source rupture.

For both southwest- and northeast-dipping planes, the focal mechanism-based rupture model results show limited data fitting with the GNSS dynamic waveforms. Compared to far-field seismic waves, near-field GNSS-based dynamic deformation is more easily affected by the geometric parameters and crustal medium structure of the seismogenic fault. The poor fitting may imply that a single-fault structure cannot describe the complex geometric parameters of the hypocenter. Relatively speaking, the proposed GNSS-based static deformation and existing InSAR source models illustrate that both the southwest- and northeast-dipping planes alone can explain the permanent deformation caused by the Jishishan earthquake (Huang X et al., 2024; Liu ZJ et al., 2024; Yang JY et al., 2024). Tang JY et al. (2023) and Fang N et al. (2024) both adopted InSAR data as constraints to examine the geometric parameters of the seismogenic fault. The former believed that the fault is southwest-dipping, whereas the latter argued that the fault is northeast-dipping, with a strike/dip of 325.2°/32.2°. Significant differences exist in studies partly because InSAR data have low signal-to-noise ratios and cannot provide sufficiently strong constraints on earthquakes due to blind thrust faults. Furthermore, these remarkable differences suggest the complexity of the ruptures involved in the Jishishan earthquake. For example, ruptures may occur on multiple faults with different geometric parameters. Based on the source rupture characteristics, the slip distribution constrained by InSAR data is relatively concentrated (Huang X et al., 2024; Liu ZJ et al., 2024; Tang XW et al., 2024), with insignificant rupture details. Fang N et al. (2024) integrated far-field P-waves and InSAR data to conduct inversion modeling. The obtained source rupture slip distribution is also relatively concentrated. No significant northeastward rupture propagation characteristics have been noted. Furthermore, the spatial distribution of the rupture differs considerably from that of aftershocks.

Based on the precise location of the aftershocks, the Jishishan earthquake may be caused by conjugate faults with opposite dip directions (Wang SG et al., 2024). Aftershocks can be divided into two categories based on their spatial distribution characteristics to obtain two cross-sections: in the south and the north. The following was obtained through fitting: The strike/dip of the first (northeast-dipping) fault is 327°/ 60°, and those of the second (west-dipping) fault are 170°/70°. The inversion results based on these conjugate faults (Figure 8) reveal a significant correlation between observed and simulated far-field P-waves and GNSS-based static and dynamic deformation values. In particular, the fitting of GNSS dynamic waveforms has been considerably improved. The source model demonstrates a significant slip on the west-dipping fault in the north, forming a bilateral rupture parallel to the northeast-dipping fault in the south. This shows high spatial consistency with the northward shift in the aftershock distribution. The maximum slip equals 0.67 m, which is similar to the result by Liu ZJ et al. (2024). However, the moment magnitude obtained via inverse modeling (6.3) is slightly overestimated, likely because of the inaccurate crustal medium structures of the source area and GNSS stations. Considering both regional tectonics and geological survey results, the seismogenic fault is determined to be a local northeast-dipping blind thrust fault, and northward rupture propagation may have caused the slip of conjugate faults. Given the limited number of aftershocks on the northern side, the rupture of the conjugate faults requires further study.

Figure  8.  Source rupture of the conjugate fault model. (a) Source slip distribution and measured and simulated GNSS values. (b) Source time function. (c) Observed and simulated far-field body waves. (d) Observed and simulated high-rate GNSS values.

DownLoad: Full-Size Img PowerPoint