Introduction

Earthquakes occurring in the Chinese mainland are frequent, of large intensity, and cover a wide range of areas. Events in western China are more numerous and on average stronger than in eastern China. The statistics of historical earthquakes show that the seismicity frequency in Western China accounts for 90% of the country seismic activity. A huge amount of high-precision and high spatial and temporal resolution data has been provided for studying crustal deformation and strong earthquake preparation processes in the western Chinese Mainland. Such data were based on the CMONOC I since 1999, and on the CMONOC II since 20101,2,3,4. GNSS observations include campaign observations and continuous observations. Campaign observations have a good spatial resolution, whereas continuous observations have a higher temporal resolution5.

Since the CMONOC observations (1999) started in Western China, there were observed two great earthquakes with M ≥ 8 (the West Kunlun Mountain Pass M 8.1 earthquake in 2001, and the Wenchuan M 8.0 earthquake in 2008), seven earthquakes with M 7.0–7.9 ( Yutian M 7.3 earthquake in 2008, Yushu M 7.1 earthquake in 2010, Lushan M 7.0 earthquake in 2013, Yutian M 7.3 earthquake in 2014, Jiuzhaigou M 7.0 earthquake in 2017, Maduo M 7.4 earthquake in 2021 and Wushi M 7.1 earthquake in 2024), more than 60 earthquakes with M 6.0-6.9, and more than 300 earthquakes with M 5.0–5.9 (shown in Fig. 1, where the secondary active blocks in the map are divided from the literature7). To deepen the application of GNSS observation data in earthquake forecasting and improve the capabilities and effectiveness of seismic risk assessment, it is necessary to understand the relationship between GNSS crustal deformation characteristics and strong earthquakes.

Fig. 1
figure 1

Spatial distribution of earthquakes above M ≥ 5 in western Chinese Mainland since 1999, the focal mechanism solutions of No.1 to No.9 correspond to the West Kunlun Mountain Pass M 8.1 earthquake occurred in 2001, the Yutian M 7.3 earthquake occurred in 2008, the Wenchuan M 8.0 earthquake occurred in 2008, the Yushu M 7.1 earthquake occurred in 2011, the Lushan M 7.0 earthquake occurred in 2013, the Yutian M 7.3 earthquake occurred in 2014, the Jiuzhaigou M 7.0 earthquake occurred in 2017, the Maduo M 7.4 earthquake occurred in 2021, the Wushi M 7.1 earthquake occurred in 2024, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image

Over the past few decades, numerous studies have explored the relationship between GNSS observation data and major earthquakes. As a result, a multi-scale method for extracting crustal deformation information from GNSS observations has been developed, focusing on its connection to the earthquake preparation process. The results have been widely applied in long, medium and short-term earthquake forecasting.

Murai et al.8 discovered that the occurrence and ___location of large earthquakes occurring in Japan and other Asian regions can be predicted based on pre-signals-changes in the area measured by triangular networks of GPS receivers. However, more research is needed for precisely forecasting the occurrence time of earthquakes. Wu et al.9 studied the characteristics of crustal deformation in the western region of China. They analyzed GNSS velocity fields, principal strain rates, maximum shear strain rates, east-west strain rates, and rotation rates to investigate the relationship between such characteristics and the preparation of strong earthquakes with a magnitude M ≥ 7.0. These results have good indicative significance for the medium and long-term forecasting of strong earthquake activity sites.

The scientific short-term forecasting of earthquakes by GNSS still poses considerable challenges to researchers. Jiang et al.10 proposed that dynamic information of relative crustal movement and deformation can be extracted using GNSS displacement time series, baseline time series, parameter time series 3-station, or multi-station strain. The GNSS time series curves have a quasilinear trend after deducting non-tectonic influences. They have a certain indicative significance for short-term earthquake forecasting when a relevant deviation from the quasilinear trend in crustal movement at a certain scale (deviation exceeding a certain limit)11 is observed. For the GNSS time series curves, there is a trend change from months to years before large earthquakes, such as the 2008 Wenchuan M 8.0 earthquake12, the 2001 West Kunlun Mountain Pass M 8.1 earthquake10, and the 2011 Japan M 9.0 earthquake13. Thus, there is also a good correspondence between the GNSS time series and the occurrence time of strong earthquakes.

Gitis et al.14 evaluated the efficiency of using the time series of displacements of the earth’s surface according to GPS data for the systematic forecasting of earthquakes. In this paper, geological data, seismic activity records, and GNSS observation data from the western region of the Chinese Mainland are collected. Abnormal data, such as those affected by multipath effects, damaged stations, or large mean square errors, are removed from the dataset. The data includes GNSS continuous station records from CMONOC (1999–2023) and GNSS campaign station data from CMONOC (1999 to 2007). This study examines the common characteristics of crustal deformation movements preceding medium and strong earthquakes with different focal mechanism types in the western Chinese Mainland. By combining seismogenic structure types with GNSS velocity field analysis, strain rate field analysis, and time-series analysis of strain parameters from multiple continuous station combinations, the research aims to identify patterns related to these earthquakes. The possible correlation between the dynamic characteristics of regional crustal deformation before medium and strong earthquakes and the late stage of earthquake preparation are preliminarily analyzed. Finally, the role of GNSS crustal deformation in the forecasting of future medium and strong earthquakes in the western region of the Chinese Mainland is discussed.

Data and methodology

Study area

China has the largest number of strong earthquakes in the world. 33% of the world’s continental strong earthquakes have occurred in China, whereas its land area covers only 7% of the planet15. The western part of the Chinese Mainland contains the two active blocks of Qinghai-Xizang Plateau and the Western Regions7. Large-scale crustal deformations are generated due to the compression, uplift, strike-slip, and thrust generated by the convergence of the Indian and Eurasian plates16,17. The northward subduction and compression of the Indian Plate is the source of block folding movement and deformation in the western Chinese Mainland. The Qinghai-Xizang Plateau is compressed, uplifted, and thickened in the near NS direction and extended in the EW direction under the action of the Indian plate, thus causing the Qinghai-Xizang Plateau and its marginal areas to become high-incidence areas of large M ≥ 7 earthquakes in the Chinese Mainland. The Secondary active blocks in the study area mainly including the areas of Lhasa, Qiangtang, Bayan Har, Qaidam, Qilian, Sichuan-Yunnan, Western Yunnan, Southern Yunnan, Tarim, Tianshan, Zhungeer, Alxa, Ordos and South China (see Fig. 1). Among them, Lhasa, Qiangtang, Bayan Har, Qaidam, Qilian, Sichuan-Yunnan, Western-Yunnan, and Southern-Yunnan jointly constitute the Qinghai-Xizang Plateau. The large-scale fault zones mainly include the Himalayan main thrust belt, the Karakorum-Jiali fault, the Mani-Yushu-Xianshuihe fault, the Kunlun-Maqin fault, the Altyn-Haiyuan fault, and the Jinshajiang-Red River fault7.

Strong earthquakes

In western China, there have been two great earthquakes with M ≥ 8 and seven earthquakes with M 7.0–7.9 since 2001. Table 1 provides the numbers and information of these earthquakes. The West Kunlun Mountain Pass M 8.1 earthquake occurred in the Hoh Xil region at the border of Qinghai and Xinjiang. The macroscopic epicenter was located in the East Kunlun sinistral strike-slip active tectonic zone at the boundary between the Bayan Har block and the Qaidam block18. According to the CENC, the epicenter of the earthquake was 36.2° N, 90.9° E. The strike and dip values of the obtained two nodal planes were 94°/61° and 190°/ 80°, with a focal depth of 15 km (https://www.globalcmt.org/CMTsearch.html). The type of focal mechanism solution for this earthquake was a sinistral strike-slip movement with thrust19,20.

Table 1 Statistical information of strong earthquakes in the study area.
Full size table

The Wenchuan M 8.0 earthquake occurred in Wenchuan County, Sichuan Province, with a macroscopic epicenter located on the NE-trending Longmenshan fault Zone at the junction of the Sichuan Basin and the Bayan Har Block21,22,23. According to the CENC, the epicenter of the earthquake was 30.95° N, 103.40° E. The strike and dip values of the obtained two nodal planes were 231°/35° and 357°/68°, with a focal depth of 12.8 km (https://www.globalcmt.org/CMTsearch.html). The Wenchuan earthquake was a thrust movement with a dextral strike-slip24.

The Yutian M 7.3 earthquake occurred in 2008 had a macroscopic epicenter located at the intersection of the Altyn Tagh Fault Zone and the West Kunlun Fault Zone25. According to the CENC, the epicenter of the earthquake was 35.60° N and 81.60° E. The strike and dip values of the obtained two nodal planes were 358°/41° and 203°/52°, with a focal depth of 12 km (https://www.globalcmt.org/CMTsearch.html). The mechanism of this earthquake was a normal movement with a slight strike-slip.

The Yushu M 7.1 earthquake had a macroscopic epicenter located on the Ganzi-Yushu fault at the tectonic junction of the Qiangtang block and the Bayan Har block26,27,28,29. According to the CENC, the epicenter of the earthquake was 33.22° N, 96.59° E. The strike and dip values of the obtained two nodal planes were 296°/65° and 205°/89°, with a focal depth of 16 km (https://www.globalcmt.org/CMTsearch.html). In this case, the focal mechanism solution type was a sinistral strike-slip.

The Lushan M 7.0 earthquake was located in the middle and southern sections of the Longmenshan fault zone at the junction of the Bayan Har block and the South China block30,31,32. According to the CENC, the epicenter of the earthquake was 30.30° N, 103.00° E. The strike and dip values of the obtained two nodal planes were 212°/42° and 19°/49°, with a focal depth of 22 km (https://www.globalcmt.org/CMTsearch.html). The focal mechanisms of the Lushan earthquake result from a pure thrust and are consistent with the characteristics of the Longmenshan fault structure.

The Yutian M 7.3 earthquake occurred in 2014 was located in the intersection area of the West Kunlun fault zone and the Altyn Tagh fault zone33,34,35. According to the CENC, the epicenter of the earthquake was 36.10° N, 82.50° E. The strike and dip values of the two nodal planes were 332°/85° and 242°/86°, with a focal depth of 18 km (https://www.globalcmt.org/CMTsearch.html). Its focal mechanism was a sinistral strike-slip.

The Jiuzhaigou M 7.0 earthquake occurred on the eastern boundary of the Bayan Har block36. According to the CENC, the epicenter of the earthquake was 33.20° N, 103.82° E. The strike and dip values of the two nodal planes were 151°/79° and 243°/82°, with a focal depth of 16 km (https://www.globalcmt.org/CMTsearch.html). Its focal mechanism was a sinistral strike-slip.

The Maduo M 7.4 earthquake, with a macroscopic epicenter located within the Bayan Har block, occurred near the northeast boundary of the East Kunlun Fault Zone37,38. According to the CENC, the epicenter of the earthquake was 34.59° N and 98.34° E. The strike and dip values of the two nodal planes were 283°/87° and 13°/84°, with a focal depth of 12 km (https://www.globalcmt.org/CMTsearch.html). Its rupture mode showed a sinistral strike and tension.

The Wushi M 7.0 earthquake had a macroscopic epicenter located in the Tianshan seismic zone (https://www.sohu.com/a/753697575_162758). According to the CENC, the epicenter of the earthquake was 41.19° N, 78.57° E. The strike and dip values of the two nodal planes were 236°/48° and 110°/57°, with a focal depth of 14 km (https://www.globalcmt.org/CMTsearch.html). Its rupture mode showed pure thrust.

GNSS data

GNSS measurements can be divided into continuous observations and campaign observations. The continuous observation mode refers to the installation of GNSS equipment on ground permanent stations while the data is continuously stored in the receiver storage device or transmitted back to the data center through a network. With this setup, the main output will be GNSS time series results. Differently, campaign observations adopt the mode of multi-stage retesting at certain time intervals (generally longer than a year) on ground permanent stations, where each observation phase is approximately 96 h. With this setup, the main output will be GNSS velocity field results. Campaign observations provide high spatial resolution for studying crustal deformation. While continuous observations have lower spatial resolution, they compensate for the limited temporal resolution of campaign observations by offering a more accurate reflection of the dynamic temporal processes of crustal deformation12,39,40,41. To better analyze the dynamic change of crustal deformation before strong earthquakes in the western Chinese Mainland since GNSS observations, the results of CMONOC I and CMONOC II observations are selected for velocity field analysis, strain rate field analysis, and strain time series analysis.

Methodology

The GNSS velocity field can directly reflect the crustal relative movement and deformation42,43. When the crustal movement and deformation are analyzed in a region, the rigid motion without deformation will generally be removed in the study area. The South China block is the largest relatively stable block on the Chinese Mainland. Therefore, it is considered the best choice to take it as reference datum when the movement and deformation of the western Chinese Mainland are used by GNSS data44,45,46.

The GNSS strain rate field is directly calculated using the partial derivative relationship between strain and displacement. Melchier47 and Savage et al.48 mentioned the differential expression of spherical displacement and strain, as shown in formula (1). In the equations, φ is latitude, λ is longitude, and R is the average radius of the earth.

$$\left\{ \begin{gathered} {\varepsilon _\varphi }=\frac{1}{R}\frac{{\partial {u_\varphi }}}{{\partial \varphi }} \hfill \\ {\varepsilon _\lambda }=\frac{1}{{R\cos \varphi }}\frac{{\partial {u_\lambda }}}{{\partial \lambda }} - \frac{{{u_\varphi }}}{R}\tan \varphi \hfill \\ {\varepsilon _{\lambda \varphi }}=\frac{1}{2}\left[ {\frac{1}{{R\cos \varphi }}\frac{{\partial {u_\varphi }}}{{\partial \lambda }}+\frac{{{u_\lambda }}}{R}\tan \varphi +\frac{1}{R}\frac{{\partial {u_\lambda }}}{{\partial \varphi }}} \right] \hfill \\ \end{gathered} \right.$$
(1)

For the GNSS strain rate field, the single-period results can be analyzed from the characteristics of strain rate magnitude, strain rate direction, shear/tension, and compression properties. The results are used to assess the current deformation intensity and its relationship to the tectonic deformation background, along with its connection to tectonic fractures. Multi-phase strain rate field analysis allows for a deeper investigation into the strengthening or weakening of regional tectonic deformation. However, the calculation of strain rates directly from station velocities is influenced by the uneven distribution of GNSS sites, with site density playing a key role49. The least squares collocation method, in contrast, offers significant advantages in terms of robustness, mitigating edge effects, error distribution, and stability50,51. Because the relationship between GNSS strain rates is established using the least-squares collocation method, the strain rates εφ, ελ, and ελφ can be calculated directly from formula (1), and their errors can also be estimated. Additional strain rate components were calculated using formula (2)52, in which ε1 is the principal tensile strain rate, ε2 is the principal compressive strain rate, γmax is the maximum shear strain rate, γ1 is the first shear strain rate, εareal is the plane strain rate, and τ2inv is the second invariant of strain rate.

$$\left\{ \begin{gathered} {\varepsilon _1}=\frac{{{\varepsilon _\lambda }+{\varepsilon _\varphi }}}{2} - \frac{1}{2}\sqrt {4\varepsilon _{{\lambda \varphi }}^{2}+{{\left( {{\varepsilon _\lambda } - {\varepsilon _\varphi }} \right)}^2}} \hfill \\ {\varepsilon _2}=\frac{{{\varepsilon _\lambda }+{\varepsilon _\varphi }}}{2}+\frac{1}{2}\sqrt {4\varepsilon _{{\lambda \varphi }}^{2}+{{\left( {{\varepsilon _\lambda } - {\varepsilon _\varphi }} \right)}^2}} \hfill \\ {\gamma _{\hbox{max} }}=\frac{{{\varepsilon _1} - {\varepsilon _2}}}{2} \hfill \\ {\gamma _1}={\varepsilon _\lambda }-{\varepsilon _\varphi } \hfill \\ {\varepsilon _{areal}}={\varepsilon _1}+{\varepsilon _2} \hfill \\ {\tau _{2inv}}=\sqrt {\varepsilon _{\lambda }^{2}+\varepsilon _{\varphi }^{2}+2\varepsilon _{{\lambda \varphi }}^{2}} \hfill \\ \end{gathered} \right.$$
(2)

The magnitude of the principal strain rate (the principal tensile strain rate ε1, the principal compressive strain rate ε2) reflects the corresponding change in principal stress. The value of the maximum shear strain rate γmax reflects the change strength of the maximum shear stress. The first shear strain γ1 reflects the shear deformation of fault in the northeast and northwest directions. The EW strain rate ελ represents the longitudinal strain rate, and the NS strain rate εφ represents the latitudinal strain rate. Negative values of the plane strain rate εareal represent the contraction rate of a region, while positive values represent the expansion rate of a region. The second invariant of strain rate mainly characterizes the total deformation of the study area.

The GNSS strain time series refers to the time series of continuous deformation parameters of small regional deformation units composed of multiple (3 or more) GNSS continuous stations in the seismogenic area calculated using the observation data of GNSS reference stations and the strain model53,54. The strain time series calculation method based on multi-station combination can effectively study the allocation characteristics of the local strain39,54,55,56.

Results

The velocity field

The intensity (magnitude) of earthquakes is significantly related to the scale of fracture development. Due to the large number of earthquakes of M 5.0–5.9 and M 6.0–6.9, and the fact that the spatial distribution is relatively scattered, the movement and variation characteristics of the velocity field before these earthquakes are not obvious. Earthquakes with M ≥ 7, except the Wushi M 7.1 earthquake, occurred around the Bayan Har block on the Qinghai-Xizang Plateau in the western Chinese Mainland57. The seismogenic structure was generally dominated by typical deep and large faults and had a relatively obvious relationship with the characteristics of velocity field movement. Therefore, the deformation characteristics of the velocity field before earthquakes with M ≥ 7 in the western region of the Chinese Mainland since 2001 are analyzed. It should be noted that the West Kunlun Mountain Pass M 8.1 earthquake occurred after the 1999–2001 velocity field observation, which can be used to study the crustal movement before the earthquake. Differently, the other earthquakes with M ≥ 7 occurred after 2007. The 1999–2007 velocity field is a relatively stable background field of crustal movement in the Chinese Mainland. Since it was not only eliminated the effect of the 2001 West Kunlun Mountain Pass M 8.1 earthquake, but also not affected by the 2008 M 8.0 Wenchuan earthquake, it can be used to study the crustal movement before the aforementioned earthquakes with M ≥ 7. Therefore, the 1999–2001 and 1999–2007 velocity fields are presented in this study based on the South China block45,46,58 (see Fig. 2).

Fig. 2
figure 2

GNSS velocity field in western Chinese mainland (South China datum). The black arrow in the figure represents the velocity field from 1999 to 2001; The red arrow represents the velocity field from 1999 to 2007. The focal mechanism solutions of No.1 to No.9 correspond to the West Kunlun Mountain Pass M 8.1 earthquake, the Yutian M 7.3 earthquake in 2008, the Wenchuan M 8.0 earthquake, the Yushu M 7.1 earthquake, the Lushan M 7.0 earthquake, the Yutian M 7.3 earthquake in 2014, the Jiuzhaigou M 7.0 earthquake, the Maduo M 7.4 earthquake, the Wushi M 7.1 earthquake, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image

Few stations are present in the central and western parts of the Bayan Har and Qaidam blocks for the velocity fields of 1999–2001 and 1999–2007. The number of velocity field stations in the Lhasa and Qiangtang blocks in the period 1999–2001 is significantly less than that in the period 1999–2007. The characteristics of the velocity fields in the two periods are gradually decreasing from southwest to north, northeast, and southeast. Considering data at the same latitude, the velocity values in the west are significantly greater than those in the east. Compared with the Lhasa and Qiangtang blocks of the Qinghai-Xizang Plateau, the velocity values in other regions are significantly smaller. This becomes particularly evident in the northwest, northeast, east, and southeast boundaries of the Qinghai-Xizang Plateau. These areas with large far-field velocity and small near-field velocity can easily accumulate strain energy, thus making them relatively prone areas to experience earthquakes.

The comparison of velocity fields at the same ___location in the two periods shows that the value of the velocity field in the 1999–2007 period is significantly larger than that in the 1999–2001 period moving from the southern part to Lhasa, Qiangtang, Bayan Har, and the Qaidam blocks in the northeastern of the Qinghai-Xizang Plateau, while its direction change is also not obvious. However, it is not apparent in the relatively dense areas of the northwest, northeast, and southwest edges of the Qinghai-Xizang Plateau. In the northwestern margin of the Qinghai-Xizang Plateau (the southern part of the Tarim Basin), the rate value of the period 1999–2007 is significantly larger, and the direction is deflected to the NW direction when compared with the previous period of NWW. In the northeastern margin of the Qinghai-Xizang Plateau (the Qilian block, the suture zone between the Qilian block and the Alxa block, and the eastern part of the Qaidam block), the changes in the value of the velocity field between the two periods are not obvious, although they significantly differ in direction. During the 1999–2001 period, the main movements were in the NNW and near EW directions. However, from 1999 to 2007, there was a clockwise rotation of about 45° west of 97° N and a clockwise deflection of approximately 90° east of 97° N, based on the earlier period’s data. In the eastern margin of the Qinghai-Xizang Plateau (the eastern region of the Bayan Har block), the rate values slightly increased in the 1999–2007 period compared to the 1999–2001 period, while the rate values near the eastern boundary of the Bayan Har block near the Wenchuan M 8.0 earthquake decreased. The SE and SSE movements have a counterclockwise deflection in the direction of motion. To the southwestern margin of the Qinghai-Xizang Plateau, the velocity field of the Sichuan-Yunnan, western Yunnan, and southern Yunnan blocks generally has a clockwise rotation of SE to SW from north to south.

There is no GNSS station around the epicenter of the West Kunlun Mountain Pass M 8.1 earthquake in the 1999–2001 velocity field. Based on the 1999–2007 velocity field containing the earthquake and the post-earthquake movement, it was found that the earthquake may have caused the unlocking of the boundary between the Bayan Har block and the Qaidam block9, resulting in increased values of the velocity field in a wide range, especially in Lhasa, Qiangtang, Bayan Har and Tarim blocks. However, the velocity values of the eastern boundary of the Bayan Har block, where the epicenters of the Wenchuan M 8.0 earthquake and the Lushan M 7.0 earthquake are located, become smaller. This indicates that the region may have a high degree of pre-earthquake locking. The occurrence of the Wenchuan M 8.0 earthquake might have caused the unlocking of the eastern boundary of the Bayan Har block59. In addition, the occurrence of the two large earthquakes with magnitude 8 has accelerated the rate of strain accumulation around the Bayan Har block, resulting in multiple earthquakes M ≥ 7 occurring at the boundary of the block.

The strain rate field

The least-squares collocation method is adopted to calculate the principal strain rate, the maximum shear strain rate, the first shear strain rate, the EW linear strain rate9, the NS linear strain rate, the second invariant of strain rate16,60 and the plane strain rate44,61,62 from GNSS velocities of 1999–2001 and 1999–2007. The results of the maximum shear strain rate field, the first shear strain rate, the EW linear strain rate field, the NS linear strain rate field, the plane strain rate field, and the second invariant of strain rate field are shown from Figs. 3, 4, 5, 6, 7 and 8. The results of the 1999–2001 period (Figure (a)) for each strain parameter simultaneously provide the corresponding principal strain rate field and the focal mechanism solution of all earthquakes M ≥ 5 after the observation period from January 1, 2002 to January 31, 2024. The 1999–2007 period results (Figure (b)) also give the corresponding principal strain rate field and the focal mechanism solutions of all earthquakes M ≥ 5 after the observation period from January 1, 2008 to January 31, 2024. Because the West Kunlun Mountain Pass M 8.1 earthquake occurred after the 1999–2001 observation, and it was mainly a sinistral strike-slip earthquake, the strain rate field characteristics before the earthquake can be analyzed by the maximum shear strain rate field and the first shear strain rate field of the 1999–2001 period. Table 2 shows the number of M 5-M 8 earthquakes in the 1999–2001 period (Figure (a)) and 1999–2007 period (Figure (b)) of each parameter strain rate field.

Fig. 3
figure 3

The maximum shear strain rate (colored contours) and the principal strain rate (black arrows) field of the Western Chinese mainland, (a) the maximum shear strain rate and the principal strain rate field of the 1999–2001 period; (b) the maximum shear strain rate invariant and the principal strain rate field of the 1999–2007 period, the focal mechanism solutions of No.1, No.4, No.6, No.7, and No.8 correspond to the West Kunlun Mountain Pass M 8.1 earthquake, the Yushu M 7.1 earthquake, the Yutian M 7.3 earthquake in 2014, the Jiuzhaigou M 7.0 earthquake, and the Maduo M 7.4 earthquake, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image
Fig. 4
figure 4

The first shear strain rate (colored contours) and the principal strain rate (black arrows) field of the Western Chinese mainland, (a) the first shear strain rate and the principal strain rate field of the 1999–2001 period; (b) the first shear strain rate and the principal strain rate field of the 1999–2007 period, the focal mechanism solutions of No.1, No.4, No.6, No.7, and No.8 correspond to the West Kunlun Mountain Pass M 8.1 earthquake, the Yushu M 7.1 earthquake, the Yutian M 7.3 earthquake in 2014, the Jiuzhaigou M 7.0 earthquake, and the Maduo M 7.4 earthquake, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image
Fig. 5
figure 5

The EW orientation strain rate (colored contours) and the principal strain rate (black arrows) field of the Western Chinese mainland, (a) the EW orientation strain rate and the principal strain rate field of the 1999–2001 period; (b) the EW orientation strain rate and the principal strain rate field of the 1999–2007 period, the focal mechanism solutions of No.2, No.3, No.5, and No.9 correspond to the Yutian M 7.3 earthquake in 2008, the Wenchuan M 8.0 earthquake, the Lushan M 7.0 earthquake, and the Wushi M 7.1 earthquake, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image
Fig. 6
figure 6

The NS orientation strain rate (colored contours) and the principal strain rate (black arrows) field of the Western Chinese mainland, (a) the NS orientation strain rate and the principal strain rate field of the 1999–2001 period; (b) the NS orientation strain rate and the principal strain rate field of the 1999–2007 period, the focal mechanism solutions of No.2, No.3, No.5, and No.9 correspond to the Yutian M 7.3 earthquake in 2008, the Wenchuan M 8.0 earthquake, the Lushan M 7.0 earthquake, and the Wushi M 7.1 earthquake, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image
Fig. 7
figure 7

The plane strain rate (colored contours) and the principal strain rate (black arrows) field of the West Chinese mainland, (a) the plane strain rate and the principal strain rate field of the 1999–2001 period; (b) the plane strain rate and the principal strain rate field of the 1999–2007 period, the focal mechanism solutions of No.2, No.3, No.5, and No.9 correspond to the Yutian M 7.3 earthquake in 2008, the Wenchuan M 8.0 earthquake, the Lushan M 7.0 earthquake, and the Wushi M 7.1 earthquake, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image
Fig. 8
figure 8

The second invariant strain rate (colored contours) and the principal strain rate (black arrows) field of the Western Chinese mainland, (a) the second invariant strain rate and the principal strain rate field of the 1999–2001 period; (b) the second invariant strain rate and the principal strain rate field of the 1999–2007 period, the focal mechanism solutions of No.1 to No.9 correspond to the West Kunlun Mountain Pass M 8.1 earthquake, the Yutian M 7.3 earthquake in 2008, the Wenchuan M 8.0 earthquake, the Yushu M 7.1 earthquake, the Lushan M 7.0 earthquake, the Yutian M 7.3 earthquake in 2014, the Jiuzhaigou M 7.0 earthquake, the Maduo M 7.4 earthquake, the Wushi M 7.1 earthquake, respectively. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image
Table 2 Statistics number of earthquakes in different periods in the study area.
Full size table

The maximum shear strain rate (Fig. 3) and the first shear strain rate (Fig. 4) reflect the strike-slip state of the faults. Therefore, they can be used to determine the relationship between the strike-slip-dominated earthquakes with these two parameters. The EW strain rate (Fig. 5) and the NS strain rate (Fig. 6) reflect the compression and tension states of the crust in different directions. The plane strain rate (Fig. 7) reflects the state of regional contraction or expansion, while all the three parameters together can reflect the dip-slip motion of faults and can be used to describe the relationship between the dip-slip earthquakes. The second invariant of strain rate (Fig. 8) mainly characterizes the total deformation of the study area, so the relationship between all types of earthquakes with this parameter can be found.

Characteristics of the corresponding relationship between the GNSS strain rate field and earthquakes with M ≥ 5 in the western Chinese Mainland are obtained according to the results of each GNSS strain rate field, as shown in Figs. 3, 4, 5, 6, 7 and 8.

9 earthquakes with M ≥ 7 have occurred in the study area since 2000, 5 of which were dominated by strike-slip ruptures, while the other 4 were dominated by dip-slip ruptures. The 5 strike-slip earthquakes are the West Kunlun Mountain Pass M 8.1 earthquake, the Yushu M 7.1 earthquake, the Yutian M 7.3 earthquake in 2014, the Jiuzhaigou M 7.0 earthquake, and the Maduo M 7.4 earthquake. The 4 dip-slip earthquakes are the Yutian M 7.3 earthquake in 2008, the Wenchuan M 8.0 earthquake, the Lushan M 7.0 earthquake, and the Wushi M 7.1 earthquake.

The relation between M 7–8 strike-slip earthquakes and the GNSS strain rate field is obtained based on the analysis of Figs. 3 and 4. The maximum shear strain rate results (Fig. 3) show that 5 M 7–8 strike-slip earthquakes occurred at the intersection of the edge areas of the high-value of the maximum shear strain rate and the large strike-slip fault zone, with a range of 0.6 × 10− 8/year ≤ γmax ≤ 2.0 × 10− 8/year. The results of the first shear strain rate (Fig. 4) show that the sinistral strike-slip earthquakes No.4, No.7, and No.8 occurred at the edge areas of the high negative value of the first shear strain rate, with a range of − 1.2 × 10− 8/year ≤ γ1 ≤ − 1.6 × 10− 8/year. This indicates that the high-value or edge areas of the shear strain rate consistent with the tectonic activity background of the main fault is prone to large earthquakes.

The relationship between M 5–6 strike-slip earthquakes and GNSS strain rate field is also investigated. The maximum shear strain rate result (Fig. 3a,b) shows that the high value and its marginal areas were prone to strike-slip earthquakes, with a range of 0.8 × 10− 8/year ≤ γmax ≤ 3.0 × 10− 8/year. The relationship between the first shear strain rate and earthquakes below M 7 is not significant. Therefore, the first shear strain rate parameter has weak earthquake forecasting significance for M<7.

The relationship between M 7–8 dip-slip earthquakes and the GNSS strain rate field is drawn by analyzing Figs. 5 and 6. The results of the EW strain rate (Fig. 5) show that the Yutian M 7.3 (No.2) earthquake occurred with a tension-type on the western edge areas of the high positive value of EW strain in the central and western Qinghai-Xizang Plateau, that is, the western edge of areas with high-value tension, ελ = 0.1 × 10− 8/year (Fig. 5b). The two thrust earthquakes of the Wenchuan M 8.0 (No.3) and the Lushan M 7.0 (No.5) occurred in the east edge areas of the high negative value of the EW strain rate, that is, the eastern edge of the high-value areas of east-west compression, with a range of − 2.2 × 10− 8/year ≤ ελ ≤ − 2.4 × 10− 8/year (Fig. 5b). The results of the NS strain rate (Fig. 6) show that the Wushi M 7.1 thrust earthquake (No.9) occurred at the edge area of the high negative value of the NS strain rate, that is, the eastern edge of the NS direction compression high-value area, εφ = − 2.8 × 10− 8/year (Fig. 6b). This is consistent with the fact that the Wenchuan M 8.0 earthquake (No.3) is located in the area of the eastward extrusion movement of the crust, and the Wushi M 7.1 (No.9) earthquake is located in the area dominated by the near north-south extrusion movement, as shown by the red arrow of 1999–2007 velocity field in Fig. 2. The results of the plane strain rate (Fig. 7) showed that the three thrust earthquakes, including the Wenchuan M 8.0 (No.3), the Lushan M 7.0 (No.5), and the Wushi M 7.1 (No.9), all occurred at the eastern edge areas of the negative high-value of plane strain, with a range of − 1.4 × 10− 8/year ≤ εareal ≤ − 3.4 × 10− 8/year (Fig. 7b). Differently, the Yutian M 7.3 tension earthquake (No.2) occurred at the edge of the positive high-value area of plane strain rate, εareal = − 0.1 × 10− 8/year (Fig. 7b).

The relationship between dip-slip earthquakes of M 5–6 and GNSS strain rate field was also investigated. The crustal movement in the Qinghai-Xizang Plateau region is mainly east-west (Fig. 2). The results of EW strain rate (Fig. 5) shows that the M 5–6 tension-dominated earthquakes were prone to occur in the positive high-value and its edge areas of EW strain rate of the Qinghai-Xizang Plateau region, with a range of 0 × 10− 8/year ≤ ελ ≤ 2.6 × 10− 8/year. The thrust earthquakes were prone to occur at the edge areas of negative high-value, with a range of 0 × 10− 8/year ≤ ελ≤-2.1 × 10− 8/year (Fig. 5a,b). The results of the NS strain rate (Fig. 6) show that the M 5–6 dip-slip earthquakes occurred in the negative high low-value areas in the Tarim, Tianshan, and Zhungeer regions, which are mainly characterized by north-south movement, with a range of 0 × 10− 8/year ≤ εΦ≤-3.8 × 10− 8/year (Fig. 6a,b). Figure 7 shows that the M 5–6 thrust earthquakes are located in the high negative value and its edge areas of the plane strain rate, with a range of 0 × 10− 8/year ≤ εareal≤− 3.4 × 10− 8/year. The M 5–6 tensile earthquakes are located in the high positive value and its edge areas of the plane strain rate, with a range of 0 × 10− 8/year ≤ εareal≤− 3.4 × 10− 8/year (Fig. 7a,b).

The second invariant of strain rate, τ2inv, mainly characterizes the total deformation of the study area, which is different from the maximum shear strain rate and the plane strain rate. The key reason is that the principal tensile strain rate and principal compressive strain rate will not be offset by differences in symbols, and always reflect the total deformation strength of the study area. The results of the τ2inv (Fig. 8) show that almost all the M ≥ 6 earthquakes occurred in the high-value area and the edge of the high-value area of the second invariant of strain rate in the study period of the western Chinese Mainland, except for the Pishan M 6.5 earthquake on July 3, 2015, the Hutubi M 6.2 earthquake on December 8, 2016 and the Shaya M 6.1 magnitude earthquake on January 30, 2023. Ranges moved from 0.8 × 10− 8/year ≤ τ2inv ≤ 3.6 × 10− 8/year at the period of 1999–2001(Fig. 8a), and from 1.0 × 10− 8/year ≤ τ2inv ≤ 4.4 × 10− 8/year at the period of 1999–2007 (Fig. 8b), with a proportion above 93%. The M 5.0–5.9 earthquakes in the period above, except for several earthquakes in the Tarim Basin and the middle and eastern part of the Tianshan Mountains and one earthquake in the intersection area of the Ordos and Alxa blocks, occurred in the high-value area and the edge of the high-value area of the second invariant of strain rate ranges. The latter occurred from 0.8 × 10− 8/year ≤ τ2inv ≤ 3.8 × 10− 8/year at the period of 1999–2001(Fig. 8a), and from 1.0 × 10− 8/year ≤ τ2inv ≤ 4.4 × 10− 8/year at the period of 1999–2007 (Fig. 8b), with a proportion above 90%.

At the same time, we found that the color contour’s high value, its transition areas of the second invariant of strain rate (Fig. 8), and the maximum shear strain rate (Fig. 3) in the 1999–2001 and 1999–2007 periods are relatively consistent with the high value and its transition areas of the principal strain in each period. Besides, the maximum shear strain better matched with the principal strain. The M ≥ 5 earthquakes are located in the high-value area and the edge of the high-value area of these three strain parameters. Such areas are also the ones structurally relatively active and with relatively large deformation.

The strain time series

The 25 continuous stations established in the CMONOC I have been observed since 1998, while the 260 continuous stations established in the CMONOC II have been observed since August 2010. Two great earthquakes with M ≥ 8 and seven earthquakes with M 7.0–7.9 occurred in western China since 2001. However, it is difficult to extract anomalies before the earthquakes with M < 7 due to the sparsity of continuous stations in the CMONOC I of the network. Thus, 25 earthquakes with M ≥ 6 (as of September 5, 2022, including the M 6.8 Luding earthquake, and excluding aftershocks of earthquakes with M ≥ 7) and 108 earthquakes with M ≥ 5 (as of June 1, 2022, excluding aftershocks earthquakes with M ≥ 6 other earthquakes in the swarm) after 2013 were counted. Such approach covers two years of continuous observation in CMONOC II. In order to reduce the interference of periodic signals on the identification of time series anomaly information, the effects of seasonal terms, seismic events, instrument replacement, and steps caused by unknown reasons are removed from the single station time series before calculating the strain time series.

Take the Milin M 6.9 earthquake on November 18, 2017 as example. Figure 9 shows a small area enclosed by three stations including LHAZ, XZCD, and YNZD, covering the epicenter surrounded. Figure 10 shows the time series curves of the EW strain rate and the principal compressive strain rate in the area before Milin M 6.9 earthquake. The slopes of the two curves are negative, indicating that the enclosed area is in a squeezed state and has a significant slowdown in the downward trend one more years before the earthquake. The strain rate values before the turning are − 0.51 × 10− 8/a and − 0.62 × 10− 8/a, respectively, whereas the values after the turning are 0.02 × 10− 8/a and 0.04 × 10− 8/a, respectively. This abnormal phenomenon deviating from the linear trend reflects that the strain accumulation rate near the epicenter decreased significantly one more years before the earthquake. Furthermore, this reflects that there is a sign of deformation deficit in the seismogenic area before the earthquake. Consequently, a high degree of strain energy may have been accumulated before the earthquake. Therefore, there are background anomalies in the epicentral region before the earthquake.

Fig. 9
figure 9

Spatial distribution of GNSS strain time series in small areas near the epicenter of the Milin M 6.9 earthquake in 2017. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image
Fig. 10
figure 10

Time series of the strain parameters in the epicenter of the 2017 Milin M 6.9 earthquake in the region enclosed by stations LHAZ-XZCD-YNZD, (a) The EW strain rate time series, (b) The principal compressive strain rate time series. The figure was created using GMT 6 (Generic Mapping Tools) software6.

Full size image

Table 3 shows the characteristics of the pre-earthquake strain time series of an earthquake in each of the four magnitude levels of 8, 7, 6, and 5 earthquake38,39,63. Through the analysis of the characteristics of GNSS strain time series before earthquakes with M ≥ 5 in the western Chinese Mainland since 2013, it is found that abnormal changes were deviating from the linear trend before two earthquakes with M ≥ 8 and seven earthquakes with M ≥ 7 since 2001. 6 earthquakes of the 25 earthquakes with M ≥ 6 that occurred since 2013 cannot be used for calculating the strain time series due to missing numbers or lack of stations. The other 19 earthquakes had the anomalies of deviating from the linear trend in strain time series before the earthquake. 20 earthquakes among the 108 earthquakes with M ≥ 5 had the anomalies of deviating from the linear trend in strain time series before the earthquake. The abnormal form deviating from the linear trend is mainly characterized by a slowing trend. The relative variation of the strain time series deviating from the linear trend is greater than 18%, as shown from the statistical results of Table 3. Therefore, it can be used as a threshold for judging the pre-earthquake anomalies. The forecasting indexes of the strain parameters time series of earthquakes above M 5 in the study area as follows.

Table 3 Statistical table of abnormal characteristics of GNSS strain time series before earthquakes.
Full size table
  1. (1)

    The strain rate parameters of earthquakes with M ≥ 8 deviate from the linear rate with a trend slowing down or accelerating. The Advantage duration of abnormal was 6–36 months, where the advantage corresponding parameters are the principal tensile strain rate, the maximum shear strain rate, the first shear strain rate, and the plane strain rate.

  2. (2)

    The strain rate parameters of earthquakes with M 7.0–7.9 deviate from the linear rate with a decreasing trend. The Advantage duration of abnormal was 20–24 months. The advantage corresponding parameters are the first shear strain rate and the plane strain rate.

  3. (3)

    The strain rate parameters of earthquakes with M 6.0–6.9 deviate from the linear rate with a decreasing trend. The Advantage duration of abnormal was 10–24 months. The advantage corresponding parameters are the first shear strain rate, the plane strain rate, the principal tensile strain rate, and the EW strain rate.

  4. (4)

    The strain rate parameters of earthquakes with M 5.0–5.9 deviate from the linear rate with a decreasing trend. The Advantage duration of abnormal was 6–24 months. The advantage corresponding parameters are the principal tensile strain rate, the first shear strain rate, the EW strain rate, the maximum shear strain rate, the NS strain rate, and the plane strain rate.

Discussion

In this paper, several parameters and a relatively low magnitude were applied to study the crustal deformation characteristics in the western Chinese Mainland and its relationship with the preparation of M ≥ 5 strong earthquakes. The study has been performed through the analysis of GNSS velocity field, principal strain rate, maximum shear strain rate, first shear strain rate, EW strain rate, NS strain rate, second strain invariant rate and plane strain rate. We found that the motion state of the seismogenic faults of earthquakes with M ≥ 7 is consistent with the geological tectonic background, according to the results of the velocity fields 1999–2001 and 1999–2007 with South China datum (Fig. 2).

It can be seen from Figs. 3 and 4 that the relationship of the M 7–8 strike-slip earthquakes with the maximum shear strain rate and the first shear strain rate is relatively clear, both of which are located at the edge areas of the high-value of the two parameters. Moreover, most strike-slip earthquakes with M 5–6 are located in the high-value and edge areas of the maximum shear strain rates.

Figures 5 and 6, and 7 show that the relationship between the M 7–8 dip-slip earthquakes with the EW strain rate, NS strain rate, and plane strain rate is more obvious. The thrust earthquakes are prone to occur at the edge zones of the EW compression high value. This is consistent with the direction of compression motion reflected by the velocity field (such as the No.3 and No.5 earthquakes in Fig. 5), as well as the edge areas of the NS extrusion high value (such as the No.9 earthquake in Fig. 6), and the edge areas of the high positive value of plane strain rates (such as earthquakes No.3, No.5 and No.9 in Fig. 7). The tensioning earthquakes are prone to occur at the edges of areas of high positive values of EW strain rates or plane strain rates (such as the No. 2 earthquake in Fig. 5 and the No. 2 earthquake in Fig. 7). The thrust earthquakes with M 5–6 are mainly located in the edge areas of the high negative value of EW strain rates and the high negative value and its edge areas of plane strain rates or NS strain rates (west of the Tarim block). Differently, the tension earthquakes are located in the high positive value areas and their edges of EW strain rates and plane strain rates.

Figure 8 shows that 93% of M 6.0–6.9 earthquakes and 90% of M 5.0–5.9 earthquakes occurred at the high value and edge regions of the second invariant of strain rate.

The junction between the high-value edge area of the maximum shear strain rate and the large strike-slip fault, the intersection of the high-value edge area of the tensile strain rate and the large normal fault, and the intersection of the high-value edge area of the compressive strain rate and the large thrust fault, are the places that should be carefully monitored for M 7–8 strong earthquakes in the future. This conclusion is consistent with the viewpoint given by Wu9 with the results of the maximum shear strain rate and the EW strain rate. However, this paper does not provide the relationship between the characteristics of each strain rate parameter and earthquakes below M 7.0.

Regarding the relationship between the M 7–8 earthquakes and the strain rate field, the mechanism could be attributed to the strong locking effect of the causative fault during the late seismogenic period. This results in a wide distribution of deformation on both sides of the fault, with more gradual deformation occurring near the fault itself. This phenomenon has also been observed by GNSS in the Italian region64.

Differently from the M ≥ 7 earthquakes, the seismogenic structures of M 5.0–5.9 and M 6.0–6.9 earthquakes are not all dominated by typical deep faults, and their dynamic causes exhibit diversity. All earthquakes with M ≥ 5 had the following characteristics: (1) the high-value areas and their transition zones of the maximum shear strain rate were prone to strike-slip earthquakes. (2) There were relatively more extensional earthquakes that occurred in the high positive value areas and their edge areas of EW strain rate, whereas more thrust earthquakes occurred in the high negative value areas and their edge areas. (3) There were relatively more extensional earthquakes that occurred in the high positive value areas and their edge areas of plane strain rate, while more thrust earthquakes occurred in the high negative value areas and the edge areas of plane strain rate.

In general, strain variables are determined by the relative motion between observation points and do not depend on the reference datum. The magnitude and tension or compression, sinistral strike-slip, or dextral strike-slip shear of the strain are not affected by the datum of the displacement or velocity fields62. This paper systematically extracts the abnormal characteristics of GNSS strain time series before earthquakes with M ≥ 5 in the study area through small-scale area strain time series analysis. The overall performance is that the GNSS strain time series curve deviates from the long-term background trend in the 6 months to 5 years before the earthquake, where the relative variation of the strain time series deviating from the linear trend is greater than 18% (Table 3). Despite the promising results, the accuracy of this threshold will need to be verified by more earthquakes in the future. Due to the influence of short-period non-structural noise in GNSS observation time series data65, and since the baseline time series reflects the shortening or elongation of the crust along the line connecting two stations, there are limitations in characterizing the regional deformation properties and intensity. Therefore, this article did not analyze the abnormal characteristics of GNSS single station and baseline time series before strong earthquakes.

Conclusion

In this paper, the strain accumulation state in the epicenter area before strong earthquakes was obtained through the analysis of the three methods of velocity field, strain field, and strain parameter time series of earthquakes with M ≥ 5 since 2001 in the western Chinese Mainland10. The preliminary forecasting indexes and time-___location-magnitude comprehensive indicators of the three methods above were obtained through the further refinement of the characteristics of the velocity field, strain field, and strain time series parameters. The conclusions of this study are as follows.

  1. (1)

    The M 7–8 strike-slip earthquakes are prone to occur at the intersection zone between the edge of the positive high-value area of the maximum shear strain rate and the large-scale strike-slip fault zone, as well as at high shear strain rate areas or at the edge areas consistent with the active background of the main fault tectonic activity. The M 7–8 dip-slip earthquakes tend to occur at the edge of the high negative or positive value areas or the intersection areas of EW strain and NS strain, that is, at the intersection zone between the edge of the high-value areas or the intersection areas of the strain accumulation and the large-scale strike-slip fault zone, where the far-field velocity of the fault is larger and the near-field velocity is smaller.

  2. (2)

    The M 5–6 earthquakes are prone to occur in the high-value areas or at the edge of the second strain invariant rate, maximum shear strain rate, and principal strain rate. The high-value or edge regions of the maximum shear strain rate are prone to strike-slip earthquakes. The positive high-value and edge areas of the EW strain rate and the plane strain rate are prone to tensile earthquakes, while the edge regions of the negative high-value of the EW strain rate and the negative high-value and edge areas of the plane strain rate are prone to thrust earthquakes.

  3. (3)

    The time series deviation from the long-term background trend (0.5–5 years) based on continuous data can be used as a medium-term forecasting indicator.

  4. (4)

    Earthquakes with M ≥ 7 are prone to occur in the intersection of large strike-slip and dip-slip fault zones and high-strain accumulation areas.