Present-day crustal deformation based on GPS measurements in southeastern Tibetan Plateau: implications for geodynamics and earthquake hazard

Abstract

GPS-derived strain rates can provide a tight constraint in understanding tectonic evolution of the Tibetan Plateau and are helpful in seismic hazard assessment. For this purpose, we at first verify the method, proposed by Zhu et al. in 2005 and 2006, for the computation of the strain rates, and found that the approach is reasonable and accurate. Then, based on the updated GPS data, we compute strain rates in southeastern Tibetan Plateau (SETP). The computed results showed that the strain rates in the Sichuan Basin and South China Block are very small, and the high values are largely concentrated along the Sagaing Fault and Xianshuihe-Anninghe-Xiaojiang fault system. In addition, the highest principal strain rates are located around the eastern Himalayan syntaxis (EHS), with the compressive orientations perpendicular to the Main Thrust Belt. In particular, the calculated extensive deformation basically matches the normal faulting earthquakes. In general, the characteristic of the strain rate distribution is in good agreement with the tectonic structures from the geological and geophysical investigations. At last, according to the strain rates, seismicity, and tectonic structures, we estimate five most likely zones for future major earthquakes in SETP. These areas include the Zemuhe Fault, northwestern segment of the Red River Fault, and at the intersection between the Ganzi-Yushu Fault and the Xianshuihe Fault, and at the one between the Zhongdian Fault and Jinshajiang Fault. Consequently, this study is significant in deep understanding of the geodynamics in the Tibetan Plateau and earthquake hazard assessment.

Highlights

• The newly-developed method for computing strain rates is validated;

• We calculated strain rates based on updated GPS data by new method;

• Five earthquake-prone areas are estimated in southeastern Tibetan Plateau.

Keywords:

• GPS data
• strain rates
• seismicity
• earthquake-prone area
• southeastern Tibetan plateau

1. Introduction

The continuous Indian-Eurasian collision in the past ∼50 Ma, led to the uplift of the Tibetan Plateau, which manifested itself in broadly distributed active deformation and high seismicity (Figure 1). Overall, the ground surface of the Tibetan Plateau is moving toward northeast and southeast with respect to the stable Eurasia plate (Zhu and Shi 2011). Particularly, southeastern Tibetan Plateau (SETP), located between the backland of the Tibet to the west and the South China Block to the east and covering most of Sichuan and Yunnan provinces in southwestern China (Figure 1), has played a pivotal role in the tectonic evolution of the Tibetan Plateau.

Figure 1. Geological setting and distribution of earthquakes in southeastern Tibetan Plateau. Major historical earthquakes with M ≥ 6.0 occurred from Jan.1, 1900 to Mar. 30, 2023. (Earthquake data are downloaded from website: https://data.earthquake.cn/index.html).

In SETP, numerous large-scale strike-slip faults such as the Xianshuihe-Aninghe-Zemuhe-Xiaojiang Fault, the Red River Fault and the Sagaing Fault are well developed, and catastrophic earthquakes occur frequently. In the recent 120 years, more than thirty M ≥ 7.0 events have occurred in this region, shown in Figure 1. Moreover, extensive GPS observations have been carried out in SETP over the last three decades (e.g. Shen et al. 2005; Allmendinger et al. 2007; Gan et al. 2007; Zheng et al. 2017; Wang and Shen 2020; Wang et al. 2021). The GPS vectors show southward motion of the crust from the Tibet, with extraordinary clockwise rotation around the eastern Himalayan syntaxis (EHS), as shown in Figure 2.

Figure 2. Distribution of GPS velocities in southeastern Tibetan Plateau with respect to the stable Eurasian plate from 1991 to 2016 (data from Wang and Shen (2020) and Wang et al. (2021)). The distribution of large earthquakes and fault trace is the same as in Figure 1.

Over the past four decades a number of competing models have been developed to interpret the tectonic evolution and uplift of the Tibetan Plateau (e.g. Molnar and Tapponnier 1975; Tapponnier and Molnar 1977; England and McKenzie 1982; Tapponnier et al. 1982; Houseman and England 1986, 1993, 1996; Peltzer and Tapponnier 1988; Molnar et al. 1993; England and Molnar 1997, 2005; Royden et al. 1997, 2008; Clark and Royden 2000; Holt et al. 2000; Flesch et al. 2001; Replumaz and Tapponnier 2003; Zheng et al. 2017; Bischoff and Flesch 2018; Panda et al. 2018; Wang and Shen 2020; Wang et al. 2021; Li et al. Mir et al. 2023). It seems, however, that all models have their limitations and shortcomings. Differences between them mainly focus on the following issues: (1) Whether the crustal deformation behaves block-like or broadly distributed. (2) The north-south shortening of the Tibetan Plateau is mainly accommodated by crustal thickening or by eastward extrusion along numerous large-scale strike-slip faults. (3) Lower crustal channel flow is present in the whole Tibetan Plateau or on some small regions (4) Whether gravitational spreading due to excess gravitational potential energy of the high Tibet can explain the present-day crustal deformation or not. But at any rate, no single model could account for all observed phenomena including fault evolution, crustal deformation, topographic relief and seismicity.

Surface strain and seismicity in SETP can provide important constraints for the examination of different models for the tectonic evolution of the Tibetan Plateau. Unfortunately, different strain fields were obtained by different authors although almost the same GPS observational data were used. For example, by distance weighted approach, Wang and Shen (2020) calculated strain rates in the Chinese continent based on GPS measurements with 2,403 site velocities spanning 25 years. Zheng et al. (2017) also computed strain rates in and around the Chinese continent with 2,576 GPS velocities measured during 1991–2015. If we compare the results of strain rates between them, the differences are evident. The highest maximum compressive principal strain rate from Wang and Shen (2020) is located near the eastern Himalaya syntaxis, but it is ambiguous where the highest maximum compressive principal strain rate is in the result of Zheng et al. (2017). In addition, the maximum shear strain rate is found along the Sagaing Fault in the result of Zheng et al. (2017), but it is at some other places from Wang and Shen (2020）. For detail, refer to supplementary material Figure S1. Another example is the strain rates calculated in the Tibetan Plateau by Ge et al. (2015), in which strain rates varies with the computing parameter Wt (sum of the reweighting coefficients). In Figure S2(a), we could observe that principal strain rates change with Wt, and the results from Ge et al. (2015) are also much different from Zheng et al. (2017) and from Wang and Shen (2020).

The reason why different strain rates are calculated by different authors with almost the same GPS data is mainly because GPS sites are not evenly distributed in space (Zhu and Shi 2011; Zhu 2022). Fortunately, Zhu et al. (2005, 2006) developed the method to calculate strain rate by GPS measurement data, and presented the spatial distribution of strain rates in the Chinese Mainland (Zhu and Shi 2011) and in Ordos Block (Zhu et al., 2022; Zhu, 2022). The modeled results confirmed that the method is robust and reasonable in computation, and the computed strain rates coincide with the corresponded geological structures and geophysical data such as earthquake focal mechanisms and borehole surveying. Moreover, the strain rate results are helpful for assessment of future seismic hazards (Zhu 2022; Zhu et al. 2022).

In order to give tighter constraints for the models to explain the tectonic evolution and uplift of the Tibetan Plateau, and try to estimate the future earthquake hazard in the SETP, in the study we first validate the method for calculation of strain rates, then compute the crustal strain rate field by this new method with up-to-date GPS measurement data. We finally elaborate on geodynamic implications and estimate the future seismic hazard in southeastern Tibetan Plateau.

2. Geodynamic background

Southeastern Tibetan Plateau lies in the special region of the Tibetan Plateau, bordering the South China Block to its east, the Indian-Burma plate to its west, and overlapping with the northern part of the Indochina block, shown in Figure 1. SETP is composed of five major terranes or blocks separated by a complex fault system (Zhang et al. 2003).

The northern part of SETP is composed of the Bayan Har terrane and the eastern edge of the Qiangtang terrane. These two terranes in this region are known as the passageway for the crustal material to exit southeastward from the high plateau.

Further south is the Burma-Yunnan terrane, located to the south of the eastern Himalayan syntaxis and to the north of the Indo-China terrane. The Burma mountain subduction zone forms its western boundary. The oblique subduction results in the right-lateral movement along the Sagaing Fault.

The central region of SETP is the Sichuan-Yunnan block, the northern and eastern parts are bounded by the Xianshuihe-Zemuhe-Xiaojiang fault system and the southwestern part by the Red River Fault. The eastern margin of SETP borders the South China Block including the Sichuan Basin, which is tectonically stable and low seismicity with high seismic velocity in the lithosphere (Jun et al. 2003), having a major impact on the crustal strain in SETP.

From another point of view, southeastern Tibetan Plateau is separated by a series of active faults, among which the Xianshuihe- Anninghe-Xiaojiang fault system is the most active. This fault system, form north to south, is about 1200 km long and a few hundred meters wide. The NW-SE trending Xianshuihe Fault and the nearly NS-trending Anninghe-Zemuhe-Xiaojiang Faults are the principal faults in SETP. These large-scale left-lateral strike-slip faults promote the eastward and southeastward motion of the Tibetan crust (Wang et al. 1998). Furthermore, they are tectonically active and in high seismicity, in which at least 14 M > 7.0 earthquakes were recorded historically since 814, including the 1973 M7.5 Luhuo earthquake and the 1970 M7.3 Tonghai earthquake (Shen et al. 2005), as shown in Figure 1.

Being situated to the south of the Xianshuihe-Xiaojiang fault system, the Red River Fault is more than 1000 km in length from southeastern Tibet to the Hanoi basin, separating the South China block from the Indo-China block. More than 700 km left-lateral offset were accumulated on the Red River Fault during the Oligo-Miocene period (Leloup et al. 1995), exhibiting important evidence for the hypothesis of tectonic extrusion. NE-SW trending secondary faults are located within the Burma-Yunnan terrane, and along the faults the present-day clockwise crustal rotation occurs (Shi et al. 2018b), as displayed by the GPS velocities (shown in Figure 2). Immediately south of the Red River Fault lies the Dien Bien Phu Fault, stretching more than 150 km from Yunnan, China through northwest Vietnam into Laos and maybe linking to the Nan suture in Thailand (Lai et al. 2012).

On the southwestern side of SETP lies the N–S-striking Sagaing Fault, which is a large-scale right-lateral strike-slip continental fault, extending more than 1200 km and linking to the Andaman spreading center at its southern end (Bertrand and Rangin, 2003). Between the Red River and Sagaing faults lies the Shan Plateau, which comprises two groups of active strike-slip faults (Shi et al. 2018). The dominant group consists of the NE-trending, arcuate left lateral faults, such as the Dien Bien Phu Fault, that straddle across areas of China, Vietnam, Myanmar, Thailand and Laos. The subordinate group of NW- to nearly N-trending contains a series of right lateral strike-slip faults.

East of the Xianshuihe Fault, the Longmen Shan thrust fault system is located on the eastern margin of the Tibetan Plateau, trending N-E for over 500 km, where the 2008 Ms8.0 Wenchuan earthquake and 2013 Ms7.0 Lushan event occurred (Zhu and Zhang 2010, 2013; Zhu, 2016).

3. GPS observational data and procedure for calculation of strain rates

3.1. GPS data in southeastern Tibetan Plateau

For the purpose of monitoring crustal deformation and try to capture seismic deformation precursors used to make earthquake forecasts, China has deployed a vast number of GPS measurement sites since the 1980s. In particular, the Crustal Movement Observation Network of China Project (CMONOC) has been carried out in the recent three decades, yielding huge amounts of GPS observational data. Therefore, the GPS data used in the paper were mainly from CMONOC. For the sake of unity, we did not process the GPS original observational signals, and we only employ the previous published GPS velocity results in the study. The GPS vectors with respect to Eurasian refence frame, shown in Figure 2, are the updated data combined the data from both Wang and Shen (2020) and Wang et al. (2021).

It should be noted that the co- and post-seismic deformations on the ground surface induced by major earthquakes, such as the 2004 Sumatra–Andaman earthquake (Mw = 9.3), the 2005 Nias–Simeulue earthquake (Mw = 8.7), the 2011 Tohoku-Oki earthquake (Mw = 9.0), Japan, and the 2008 Wenchuan event (Ms = 8.0), had been eliminated from the GPS measurements (Wang and Shen 2020). Therefore, GPS vectors represent the surface motions of the secular inter-seismic period in southeastern Tibetan Plateau.

It can be seen that the GPS sites densely distributed in some regions and sparsely scattered in other regions such as within the boundaries of Vietnam and Burma, shown in Figure 2. In general, the distribution of the GPS vectors agrees well with neo-tectonic structures. In the northern part of SETP, the GPS vectors direct eastwards, corresponding to eastward extrusion because of continuous India-Eurasia collision. Further south of the northern part, the GPS vectors turn southeastwards on the most part of SETP, possibly due to the resistance of the strong South China Block. Most notably, GPS vectors show remarkable clockwise rotation around the Eastern Himalayan Syntaxis (EHS). Overall, the magnitude of GPS velocity is larger in the west of SETP than that in the east.

3.2. Method and its verification

3.2.1. Method for calculation of strain rates base on GPS data

In order to better interpret geodynamics based on GPS measurements, quite a few authors calculated strain rate field. However, different workers presented different strain rate results even though the same GPS data were utilized. In addition to what has been mentioned above, here are a few more examples. Copley (2008) presented GPS-derived strain rates in southeastern Tibetan Plateau with the method of cubic spline interpolation.   (Jin et al. 2019)derived strain rate fields in SETP by means of the spherical wavelet-based multiscale approach, in which different-scale spherical wavelets are used for different densities of GPS stations based on weighted and damped least squares method. Zhang et al. (2019) obtained strain rates in SETP by means of interpolation method with a “spline in tension” approach. Similarly, Pan and Shen (2017) gained markedly different principal strain rates based on the same GPS data in SETP utilizing (a) the Grid-Distance weighted routine approach, and (b) the Grid-Nearest Neighbor routine technique, respectively. Using a modified least-squares algorithm and spatial weighting functions, Zhao et al. (2022) gained the different strain rates with almost the same GPS data in SETP. Unfortunately, however, these strain rate results are much different from each other. We have no idea what the correct result is. For more details, refer to Figure S2.

In conclusion, different strain rates were obtained in SETP by different authors with nearly the same GPS observational data although it is a deterministic forward computation. Most important of all, we have no idea what the correct answer is at all. We cannot determine whether the computed strain rates are right or not by field observations or physical experiments, or by some others. Therefore, we do not know which result is the best, and which model is reasonable, giving rise to different strain rate results published by different authors. For this purpose, Zhu et al. (2005, 2006) put forward the approach to calculate strain rate field based on GPS measurement data, and presented the spatial distributions of strain rates in the Chinese Mainland and in Ordos Block (Zhu et al. 2006; Zhu and Shi 2011; Zhu 2022), in which strain rates are clearly agreement with tectonic structures and geophysical observations. The main characteristic of the method is the use of Kriging technique (Deutsch and Journel 1997), which is the well-known and the most widely used in interpolation, especially it is suitable for unevenly distributed data points in space. Accordingly, Kriging interpolation is ad hoc lent itself to the processing of spatially unevenly distributed GPS data. In the approach, we at first use the Kriging to interpolate the spatially irregularly scattered GPS vectors to evenly grid sites, and then compute the strain rate in each grid cell. For the GPS vector data shown in Figure 2, some parameters in Kriging interpolation are obtained: Nugget is 0.8 mm², sill 4.6 mm², and range 415.3 km in WE trending velocity, but nugget is 0.4 mm², sill 12.1 mm², and range 706.8 km in NS trending velocity.

In order to demonstrate the accuracy and efficiency of this method, we will offer some comparative results with correct strain rates by numerical simulations.

3.2.2. Validation of the method

We can by no means directly confirm the strain rates computed from GPS measurement although we can obtain the value of strain rate through borehole surveying. Hence, we do not have any knowledge of which approach is reasonable and which result is the best, giving rise to a variety of different results of strain rates by different authors in SETP. In the following, we will corroborate the validity of the calculation method developed by Zhu et al. (2005, 2006).

On the whole, we at first construct an arbitrary function, from which surface velocities (CSV), where the GPS observational sites are located in SETP, can be calculated mathematically, and the corresponding actual strain rate tensors (ASR) are derived by mathematical operations. On the other hand, we can calculate the strain rates (CSR) by means of the method of Zhu et al. (2005, 2006) with the same data of CSV produced above. By means of comparing CSR with ASR through the same velocity data (CSV), we can tell how well the method is.

Here, surface vectors in θ- and φ-direction (east- and north-direction in spherical coordinate system) are given in accordance with the following functions. (1) $$\{\begin{array}{c}{v}_{\theta }=f_u(\theta ,\phi )\\ v_{\phi }=f_v(\theta ,\phi )\end{array}$$(1)

Thus, the components of strain rates on the spherical surface can be expressed as, (2) $$\{\begin{array}{l}{\dot{\epsilon }}_{\theta }=\frac{1}{R}\cdot \frac{\partial v_{\theta }}{\partial \theta }\\ {\dot{\epsilon }}_{\phi }=\frac{1}{R\sin \hspace{0.17em}\theta }\cdot \frac{\partial v_{\phi }}{\partial \phi }+\frac{v_{\theta }}{R}\mathit{\text{ctgθ}}\\ {\dot{\epsilon }}_{\theta \phi }=\frac{1}{2}\left[\frac{1}{R}\right(\frac{\partial v_{\phi }}{\partial \theta }-v_{\phi }\mathit{\text{ctgθ}})+\frac{1}{R\sin \hspace{0.17em}\theta }\cdot \frac{\partial v_{\theta }}{\partial \phi }]\end{array}$$(2) where R stands for the Earth radius.

Then, principal strain rates could be obtained according to equation group (2). (3) $${\dot{\epsilon }}_{1,2}=\frac{{\dot{\epsilon }}_{\theta }+{\dot{\epsilon }}_{\phi }}{2}\pm \sqrt{\frac{({\dot{\epsilon }}_{\theta }-{\dot{\epsilon }}_{\phi }{)}^2}{4}+{{\dot{\epsilon }}_{\theta \phi }}^2}$$(3) where ${\dot{\epsilon }}_{1,2}$ represent maximum and minimum principal strain rates, respectively.

In the same way, the maximum shear strain rate and the second invariant of strain rate (SR) can be written as (4) and (5). (4) $$\dot{\chi }=({\dot{\epsilon }}_1-{\dot{\epsilon }}_2)/2$$(4) (5) $$SR=\sqrt{{{\dot{\epsilon }}_{\theta }}^2+{{\dot{\epsilon }}_{\phi }}^2+2{{\epsilon }^2}_{\theta \phi }}$$(5)

In an actual simulation, equation system (1) can be taken many expressions. In general, equation set (1) is chosen as a spatial periodic function. For example, in the paper, we specify the functions as (6) $$\{\begin{array}{c}{v}_{\theta }=\mathrm{a}1\cdot \sin \hspace{0.17em}(a2\cdot (\theta +a3))+a4\cdot \cos \hspace{0.17em}(a5\cdot (\phi +a6\left)\right)\\ v_{\phi }=b1\cdot \sin \hspace{0.17em}(b2\cdot (\theta +b3))+b4\cdot \cos \hspace{0.17em}(b5\cdot (\phi +b6\left)\right)\end{array}$$(6) where v_θ and v_φ are surface velocities in north- and east-directions, respectively, with the physical unit of mm/yr. θ represents colatitude and φ denotes longitude.

According to Eq. (6)(6) $$\{\begin{array}{c}{v}_{\theta }=\mathrm{a}1\cdot \sin \hspace{0.17em}(a2\cdot (\theta +a3))+a4\cdot \cos \hspace{0.17em}(a5\cdot (\phi +a6\left)\right)\\ v_{\phi }=b1\cdot \sin \hspace{0.17em}(b2\cdot (\theta +b3))+b4\cdot \cos \hspace{0.17em}(b5\cdot (\phi +b6\left)\right)\end{array}$$(6) , we could obtain velocity vectors at the points where the GPS stations are situated in SETP as shown in Figure 2. Based on these same velocity data, we can obtain two strain rate results of ASR and CSR, respectively. To get a better sense of the performance of the model, we set up two sets of parameters to calculate and do some comparisons. In case 1, the parameters in Eq. (6)(6) $$\{\begin{array}{c}{v}_{\theta }=\mathrm{a}1\cdot \sin \hspace{0.17em}(a2\cdot (\theta +a3))+a4\cdot \cos \hspace{0.17em}(a5\cdot (\phi +a6\left)\right)\\ v_{\phi }=b1\cdot \sin \hspace{0.17em}(b2\cdot (\theta +b3))+b4\cdot \cos \hspace{0.17em}(b5\cdot (\phi +b6\left)\right)\end{array}$$(6) is given as following, in which a1=-9.0; a2 = 11.0; a3 = 2.0; a4=-10.0; a5 = 20.0; a6 = 2.0;b1 = 12.0; b2 = 15.0; b3 = 4.0; b4 = 8.0; b5 = 9.0; b6 = 26.0. And in Case 2, the parameters are a1 = 11.0; a2 = 56.0; a3 = 20.0; a4 = 16.0; a5 = 33.0; a6 = 22.0; b1=-12.0; b2 = 35.0; b3 = 40.0; b4=-28.0; b5 = 19.0; b6 = 7.0. The comparison of strain rates between ASR and CSR in Case 1 is shown in Figure 3. By carefully comparing Figure 3(a) with (b), we found that CSR is almost the same as ASR both in orientations and in magnitudes of principal strain rates, demonstrating the method of Zhu et al. (2005, 2006) may be reasonable. Figure 4 displays another comparison of strain rates in Case 2. From Figure 4, we confirm that the computed strain rates by the method of Zhu et al. (2005, 2006) match the corresponded correct strain rates well. In the study, we also performed a lot of comparisons like Case 1 and Case 2, and found the comparison results are satisfactory in any case.

Figure 3. Comparing the actual strain rates (ASR) with the computed counterparts (CSR) based on the same surface velocity data in Case 1. (a) Distribution map of principal strain rates calculated mathematically based on the Eqs. (1)–(6); (b) computed principal strain rates with the method developed by Zhu et al. (2005, 2006). It is clear that CSR in (b) is almost the same as ASR in (a).

Figure 4. Comparing the actual strain rates (ASR) with the computed counterparts (CSR) based on the same surface velocity data in Case 2. (a) Distribution map of principal strain rates calculated mathematically based on the Eqs. (1)–(6); (b) computed principal strain rates with the method developed by Zhu et al. It is evident that CSR in (b) is almost the same as ASR in (a).

In sum, from numerical verifications, we assume that the method of computation strain rates based on GPS data by Zhu et al. (2005, 2006) is reasonable and accurate.

4. Strain rate results

Figure 5 shows the distribution of principal strain rates in southeastern Tibetan Plateau. Obviously, it highlights that the principal strain rates are not homogeneously distributed in space. The strain rates are small in the Sichuan Basin and in the South China Block, whereas they are large in the west of the South China Block. The highest values are located around the EHS, where materials rotate clockwise. Also, we can observe that the orientations of the maximum compressive principal strain rates are perpendicular to the Main Boundary Thrust zone between India and Tibet in EHS, as shown in Figure 5. The next largest magnitude of the principal strain rates is distributed along the Sagaing Fault with the value as high as ∼8.0 × 10⁻⁸/yr, and the angle between orientation of the principal strain rate and the fault strike is nearly 45°, suggesting the dominant strike-slip sense of fault motion. In addition, Xianshuihe-Xiaojiang fault system has a characteristic of large strain rates, where major earthquakes occur frequently (shown in Figure 1). Meanwhile, we can observe that the directions of maximum and minimum principal strain rates along the Xianshuihe-Zemuhe-Xiaojiang fault are exactly opposite those of principal strain rates along the Sagaing Fault, corresponding the fact that one fault is a sense of left-lateral strike-slip, and the other is right-lateral. Additional interesting feature shown in Figure 5 is that one tensional region, located around northwestern end of the Xianshuihe fault at 99–102°E 29–32°N, behaves N-S extension, and the other do ∼ E–W extension at 99–102°E 26–28°N.

Figure 5. Spatial distribution of principal strain rates in southeastern Tibetan Plateau with grid size of 1.0°× 1.0° (the arrow outward stands for tensile, and the inward represents compressive).

The spatial distribution map of maximum shear strain rates is displayed in Figure 6. In the figure, we can clearly see that the high magnitudes of the maximum shear strain rates are located along the Sagaing Fault and the Xianshuihe-Xiaojiang fault system, with the values of 5–6 × 10⁻⁸/yr, and the lower ones are found in the inner part of the Sichuan Basin and in the South China Block and Indo-China Block, which is in line with tectonic structures. Dilatational strain rate (DST) reflects surface deformation, being in tension regime when it is positive, and in compressive state when it is negative. In Figure 7 is illustrated contour map of dilatational strain rates in SETP. The figure presents a complicated style of DST distribution. DST is negative around the EHS due to the collision between India and Eurasia. It is also negative along the Longmenshan thrust fault due to the resistance of eastward motion material by the Sichuan Basin. However, DST appears positive around the Jiangshajiang Fault and the Lijiang-Xiaojinhe Fault in the interior of the Sichuan-Yunnan Block.

Figure 6. Spatial contour map of maximum shear strain rates in southeastern Tibetan Plateau.

Figure 7. Spatial contour map of dilatational strain rates in southeastern Tibetan Plateau.

Figure 8 depicts the spatial contour map of the second invariant of strain rate (SR), which represents the total magnitude of the strain rate. In general, the distribution pattern of SR is roughly the same as that of maximum shear strain rate, although there are some differences between them. The large value of the SR is mainly scattered around the EHS and along the Sagaing Fault and the Xianshuihe-Xiangjiang fault system, possibly implying that there exist fault slips along these faults in the period of inter-seismic phase.

Figure 8. Spatial contour map of second invariant of strain rates (SR) in southeastern Tibetan Plateau.

5. Discussions

5.1. Comparison of the strain rates with previous results

Previous studies (e.g. Shen et al. 2005; Gan et al. 2007; Copley 2008; Pan and Shen 2017; Zheng et al. 2017; Zhang et al. 2019; Jin et al. 2019; Wang and Shen 2020; Wang et al. 2021; Zhao et al. 2022) made use of different approaches to study strain rates in southeastern Tibetan Plateau, which offered helpful clues to the geodynamics in the study area. However, their strain rate results in SETP show huge varieties from each other, and deviated from the basic feature of actual deformation. In the following, I will present some of them.

For example, Pan and Shen (2017) presented the GPS derived principal strain rates by applying (a) the Grid-Distance weighted technique, and (b) the Grid-Nearest Neighbor interpolation method. It is clear the obviously different results with different interpolation methods, shown in Figure S3(a). Particularly, there are evidently high values of principal strain rates scattered in the interior part of the rigid and stable Sichuan Basin where the surface deformation should be small (Zhu and Zhang 2013), but the small values of principal strain rates appear in and around the EHS, which did not agree with the basic tectonic features (Holt et al. 1991; Gupta et al. 2015).

Zhang et al. (2019) attained the principal strain rate tensors derived from GPS horizontal velocities by “spline in tension” technique. We can observe from Figure S3(b) that tensile-dominated strain rates appeared along the Xuanshuihe-Xiaojiang-Zemuhe strike-slip fault system, and the high tensile principal strain rates also occurred around the EHS where the Indo-Eurasia collision is located, which are not in accordance with the basic geological structures (e.g. Holt et al. 1991; Gupta et al. 2015).

Zhao et al. (2022) attained horizontal strain rates from updated GPS measurements using a modified least-squares algorithm and spatial weighting functions. The large values of strain rates are distributed along the Xianshuihe-Xiaojiang fault system, and the very small deformation is observed in and around the EHS (shown in Figure S3(b)), which is also not in agreement with the basic common sense on the deformation in the EHS (Gupta et al. 2015). Similarly, Copley (2008) presented small strain rates around the EHS based on GPS observations (shown in Figure S3(d)), which is evidently violated the basic common sense of deformation on the EHS (e.g. Holt et al. 1991; Gupta et al. 2015).

In a word, the above different results for the strain rates from different authors did not represent the basic deformation pattern of the SETP. The reason is extremely complicated, but it is chiefly due to the different methods they used and the unevenly distribution of GPS receiver sites in space. Additionally, I found that the approaches they utilized were basically according to least square principles, possibly producing large errors in the processing of the unevenly spatially scattered GPS data. In contrast, the method developed by Zhu et al. (2005, 2006) applies the Kriging interpolation, which was specially developed to process unevenly spatially distributed measurements data (Deutsch and Journel 1997).

Furthermore, there are many interpolation methods, e.g. linear, cube, spline, or natural neighbor interpolation. How do these other methods compare with the Kriging technique in interpolating the irregularly distributed GPS vectors? For this purpose, we interpolate GPS vectors in the SETP with linear, cubic, spline, natural neighbor interpolation, and the kriging method with different theoretical models for the semivariogram. In order to know which approach is the best in interpolation of GPS velocities in the SETP, we perform cross-validation for each one, and the square root of the variance between the observed GPS vectors and the predicted ones. Here, the model with the lowest variance is regarded as the best. The modeled results showed that Kriging interpolation is optimal.

In order to evaluate calculated results quantitatively, uncertainties of the strain rates are particularly computed. Figure S4 shows the spatial contour map of uncertainties for the second invariant of strain rate. It is evident that the uncertainties of SR in most of the regions are less than 15 × 10⁻⁹/yr which is much less than the absolute value of SR, as shown in Figure 8.

In conclusion, the principal strain rates computed in the paper are accurate and in good agreement with neo-tectonic structure and geological and geophysical surveys.

5.2. Implications for geodynamics and seismic hazard

In general, mild deformation is located in the Sichuan Basin and the South China Block, corresponding to the strong Sichuan Basin and the rigid South China Block, and large strain rates are scattered mainly along the Xianshuihe-Xiaojiang fault system and Sagaing Fault, matching weak fault zones reconciling most of the deformation. In addition, the largest compressive principal strain rates are distributed on the Eastern Himalayan Syntaxis, with the directions perpendicular to the main boundary thrust belt, as a result of Indo-Eurasia collision. Most of all, a striking feature of the deformation results is the appearance of two extensional regions in the central part of the Sichuan-Yunnan Block, represented by Zone A and B, respectively, in Figure 9. We can observe that the dominant tensile direction is different between A and B, with N-S direction in A and nearly E-W extension in B. Basically, this style of deformation is in agreement with earthquake focal mechanisms shown in Figure 9, in which the normal faulting earthquakes with east–west fault planes in region A, and roughly N–S nodal planes in region B. This peculiar characteristic could be interpreted as gravitational driving forces (Copley 2008). In Zone A, gravitational driving forces are arisen from surface slope with gradient from north to south. But, on the western boundary of the Zone B the surface slopes direct to southwest, and on the eastern part the slope is to the southeast. Therefore, the variation in downslope direction results in extension in the orientation nearly parallel to the relief contours (Copley 2008). Thus, roughly E-W extension deformation and normal faulting earthquakes with north–south nodal planes are distributed on Zone B with the center at ∼27°N 100°E, as shown in Figure 9.

Figure 9. Map of principal strain rates and the spatial distribution of earthquakes with normal sense (Ms ≥ 6.0) from Jan.1, 1976 to Mar. 30, 2023 in southeastern Tibetan Plateau. It is seen that the earthquakes with normal sense are mainly located in Zone a and Zone B. And extensive orientations both in Zone a and Zone B match the fault nodal planes of earthquakes with normal sense. Earthquake data are downloaded from websites (https://www.globalcmt.org/CMTsearch.html).

In general, seismicity is associated with surface deformation to some degree. But there is no one-to-one correspondence between them. High deformation, in other words, does not necessarily correspond to a major earthquake, and minor surface deformation is not necessarily followed by small ones. For example, the 2008 Ms8.0 Wenchuan earthquake occurred on the Longmenshan thrust fault where the deformation is small, as shown in Figures 5–8. Meanwhile, large deformation zone, because of quick fault slip in inter-seismic period, does not accumulate enough elastic strain energy for major earthquakes. This is because small frictional coefficient is on the fault surface when fault surface is slipping (Tang and Zhu 2020; Cui and Zhu 2022; Zhu and Cui 2022). It is the region where faults are locked that huge amounts of energy can be accumulated to produce major earthquakes. Thus, it is evident that the relationship between surface deformation and earthquakes is much complex (Stevens and Avouac 2021).

Based on our previous work (Zhu 2022), we assume that the most likely regions for future major earthquakes are the sites where lower SRs are distributed with fault surface locked in high friction, and adjacent to active faults where SRs are large in which the fault is controlled by velocity strengthening friction law (Cui and Zhu 2022; Zhu and Cui 2022) with low frictional coefficient. Moreover, seismic gaps are the zones associated with high seismic potential in the future (e.g. Mogi 1979; Kagan and Jackson 1991; Rong et al. 2003; Rajendran and Rajendran 2005; Kumar et al. 2013; Srivastava et al. 2015). For this purpose, Figure 10 displays the spatial distribution map of historical strong earthquakes (M ≥ 6.0) from Jan. 1, 1900 to Mar. 30, 2023 and recent small events (M ≥ 0.0) from Jan. 1, 2010 to Mar. 30, 2023. It should be pointed out that a large number of small events outside China’s territory maybe missed in the figure. Furthermore, we can notice from the figure that few earthquakes occurred in the Sichuan Basin and South China Block, as well as in some regions in and around EHS. Most of earthquakes occurred along the Xianshuihe-Xiaojiang fault system and in Sichuan-Yunnan Block. Therefore, based on surface deformation and seismicity in historical and recent years, five most likely earthquake-prone zones were selected in southeastern Tibetan Plateau in Figure 11, marked by red or white dotted circles.

Figure 10. Spatial distribution map of historical strong earthquakes (M ≥ 6.0) from Jan. 1, 1900 to Mar. 30, 2023 and recent events (M ≥ 0.0) from Jan. 1, 2010 to Mar. 30, 2023. Earthquake data are downloaded from websites (https://data.earthquake.cn/index.html).

Figure 11. Spatial distribution of SR and the most likely regions for future major earthquakes marked by A1, A2, A3, A4 and B1 with dotted circles.

It can be seen from Figure 11 that Zone A1 is located on the northwestern segment of the Red River Fault, where SR is small and seismicity is low (shown in Figure 10). Meanwhile, SR is large and seismicity is high in the vicinity of Zone A1. In addition, there has been no major earthquakes with M ≥ 7.0 in Zone A1 for over 100 years. Hence, Zone A1 is the candidate region for future major earthquakes.

Zemuhe Fault lies on the middle part of the Xianshuihe-Xiaojiang fault system, where SR is relatively small. Although the seismicity is high, but no major earthquakes with M ≥ 7.0 have occurred in this area. As a result, this region is assumed to be one of the most likely earthquake-prone area, called Zone A2 as shown in Figure 11.

Zone A3 is located at the intersection of multiple active faults such as Zhongdian Fault, Jinshajiang Fault and Jiali Fault. In particular, SR is small and seismicity is low with a seismic gap in it, shown in Figure 10. Accordingly, Zone A3 is a likely region for future major earthquakes.

Zone A4 is between the active Ganzi-Yushu Fault and Xianshuihe Fault where major earthquakes occur frequently. At the same time, Zone A4 is low in seismicity with small value of SR. Therefore, Zone A4 is regarded as a high risk zone for future major earthquakes.

In addition, Zone B1 is located along the Jiali Fault near the EHS, in which SR is small and seismicity is low, as shown in Figure 10. Most likely, there will be a major earthquake occurred in Zone B1 in the future.

6. Conclusion

In this study, I discussed the complexity of calculations of strain rates based on GPS observations and presented the strain rates in the SETP. Primary conclusions are as follows.

1. By means of constructing a mathematical function, we create spatial velocity data with the same locations as GPS receiver sites in southeastern Tibetan Plateau. From these artificial data, we validated the approach for the computing strain rates proposed by Zhu et al. (2005, 2006), and found that the method is accurate and reasonable in computation of strain rates.

2. The modeled results show that the strain rates calculated based on GPS data in the STEP are consistent with the geological and geophysical investigations. Particularly, the calculated extensive deformation regions are in good agreement with those of the normal faulting earthquakes.

3. We identified the most likely zones for future major earthquakes based on earthquake seismicity and SR. These zones include the Zemuhe Fault, northwestern segment of the Red River Fault, and at the intersection between the Ganzi-Yushu Fault and the Xianshuihe Fault, and at the region between the Zhongdian Fault and Jinshajiang Fault. Additionally, the region on the Jiali Fault near the EHS is a high risk area for future major earthquakes.

Acknowledgments

I appreciate the critical and constructive comments by four reviewers and the Editor (Prof. Ramesh Singh) that helped me to improve the manuscript.

Disclosure statement

The author declares no conflicts of interest relevant to this study.

Data availability

All GPS data used in the paper are from Wang and Shen (2020) and Wang et al. (2021); and all earthquake data are downloaded from websites (https://www.globalcmt.org/CMTsearch.html; https://data.earthquake.cn/index.html). No other new data were used in the current study, and all the model output data are available on request from authors.

Additional information

Funding

This work was jointly supported by the research grant from National Institute of Natural Hazards, MEMC (ZDJ2020-15), National Natural Science Foundation of China (41874060, 41574041), and by National key research and development program of China (2017YFC1500104).