Lithospheric Equilibrium and Anisotropy around the 2021 Yangbi Ms 6.4 Earthquake in Yunnan, China

0. INTRODUCTION

On May 21, 2021, a strong Ms 6.4 earthquake occurred in Yangbi County, Dali Prefecture, Yunnan Province, China. The epicenter was 25.67°N, 99.87°E and the focal depth was 8 km (Fig. 1), resulting in three deaths and serious damage to some houses. The earthquake was a shallow strike-slip event rupturing on a concealed fault (Yang et al., 2021). Fourty-two earthquakes with Ms ≥ 3.0 occurred in a week before and after the main earthquake, of which Ms 5.6 and Ms 4.5 foreshocks occurred continuously within half an hour before the main earthquake, and Ms 5.0 and Ms 5.3 strong aftershocks occurred on the day after the earthquake. The earthquake sequence showed typical foreshock-mainshock-aftershock characteristics (Zhang et al., 2021). The foreshocks that began on May 18, 2021 migrated and expanded step by step from the middle of the aftershock area to the northwest, showing a long and narrow strip structure. The foreshock sequences were dense and triggered each other, which was in line with the cascade mode of foreshock occurrence (Long et al., 2021; Su et al., 2021). The seismic sequence is 3–10 km away from the Weixi-Qiaohou fault (F17) on the northeast side. The long axis of the aftershock area is distributed in the NW-SE direction, about 22 km long, and generally presents the distribution characteristics of narrow in the northwest and wide in the southeast. The dominant focal depth is 4–8 km. The earthquakes in the northwest section are shallower and more intensive, and the events in the southeast section are deeper and more scattered, with southward deflection (Fig. 1c). The main shock is located at the northwest end of the aftershock intensive area, showing obvious dextral strike-slip and tension characteristics. The earthquake series has the characteristics of one-sided rupture in southeast direction. The seismogenic structure can be preliminarily determined as a parallel associated hidden fault of Weixi-Qiaohou fault (Long et al., 2021; Zhang et al., 2021).

Figure  1.  Gravity network in western Yunnan and location of the 2021 Yangbi Ms 6.4 Earthquake. (a) Distribution of gravity stations. Blue dots indicate the observation points. Yellow lines denote the profiles across the focal area. Red pentagram gives the position of the 2021 Yangbi Ms 6.4 Earthquake. Red frame shows the focal area (subfigure c), and the white triangles indicate the towns. Red lines present the faults. F1. Nujiang fault; F2. Nujiang west branch fault; F3. Tengchong Quaternary volcanic fault zone; F4. Sudian fault; F5. Dayingjiang fault; F6. Longchuan fault; F7. Longling-Ruili fault; F8. Wanding-Anding fault; F9. Caojian fault; F10. Lancang River fault; F11. Lanping-Yunlong fault; F12. Yongping fault; F13. Shidian fault; F14. Kejie fault; F15. Yongde fault; F16. Nanting River fault; F17. Weixi-Qiaohou fault; F18. Zhongdian-Longpan-Qiaohou fault; F19. Ailaoshan fault; F20. Honghe fault; F21. Jinggu fault; F22. Wuliangshan fault; F23. Chenghai fault; F24. Doupo-Pupeng fault; F25. Nanhua-Chuxiong fault; F26. Babianjiang fault. (b) Location of the study area. Solid curves are the primary and secondary plot boundaries. Blue frame shows the study area. (c) Earthquakes before and within one month after the 2021 Yangbi Ms 6.4 Earthquake. EQ1–EQ4 are four earthquakes with magnitude > 5 (Duan et al., 2021). Elevation data are from digital elevation model SRTM15+ (Tozer et al., 2019). Fault and active block data are from China Geospatial Dataset (https://github.com/gmt-china/china-geospatial-data/releases). Sta.. Station.

DownLoad: Full-Size Img PowerPoint

Weixi-Qiaohou fault (F18) is one of the branch faults of Honghe fault zone (F20) extending northward, with a length of about 280 km. It is connected with Jinshajiang fault in the north and Honghe fault in the narrow sense in the south. It is an important part of the southwest boundary of Sichuan-Yunnan Rhombic Block (Chang et al., 2016). Weixi-Qiaohou fault shows dextral strike-slip movement as a whole, with obvious activity characteristics in Late Quaternary, which is similar to Honghe fault in kinematics. It is considered that the fault belongs to the northern extension of the Honghe fault, is the main carrier of block movement and deformation, and bears and absorbs the movement energy transmitted by the block in the north of Sichuan and Yunnan (Chang et al., 2016). Weixi-Qiaohou fault and Zhongdian-Longpan-Qiaohou (F18) and other small and medium-sized faults intersected near the epicenter of the 2021 Yangbi Ms 6.4 Earthquake, making the fault structure in the focal area show dispersion distribution characteristics (Fig. 1). Historically, there have been many destructive earthquakes near the source area of the 2021 Yangbi Ms 6.4 Earthquake, including the Midu Earthquake with Ms 7.0 in 1652, Dali Earthquake with Ms 7.0 in 1925 and Maden Earthquake with Ms 6.0 in 1948. Over the past decade, the small and medium-sized seismicity in this region has been continuously enhanced, with Eryuan Earthquake swarm in 2013 and Yangbi Earthquake swarm in 2017 (Long et al., 2021; Pan et al., 2019). In a word, the source area of the 2021 Yangbi Ms 6.4 Earthquake belongs to a zone with strong seismicity.

Deep structure of the focal area can provide important information about the deep seismogenic environment and seismogenic mechanism. Ye et al. (2021) obtained the three-dimensional resistivity structures of the source area of the 2021 Yangbi Ms 6.4 Earthquake through the inversion of magnetotelluric data. Their results show that there is significant electrical lateral heterogeneity in the source area of the great event, and the main earthquake and its sequence occur on the high resistance side near the transition zone of high and low electrical resistance. Li et al. (2021) found significant lateral differences in the crustal medium in the focal area. Thus they deemed that the inhomogeneous distribution of crustal structure is the deep structural factor controlling the distribution of the 2021 Yangbi Ms 6.4 Earthquake. Based on the joint observations of in-situ gravity/GNSS (Global Navigation Satellite System) around the 2021 Yangbi Ms 6.4 Earthquake, this paper carries out the research on lithospheric density structure, lithospheric equilibrium, as well as shear wave splitting, and intends to explore the occurrence mechanism of the great event from the perspective of lithospheric equilibrium and seismic anisotropy.

1. IN-SITU MEASUREMENTS AND GRAVITY ANOMALIES AROUND THE 2021 YANGBI Ms 6.4 EARTHQUAKE

The gravity and elevation data used in this paper are the data of 318 points observed in western Yunnan, China from August to October 2020. The whole observation network covers 97.5°E–101.5°E and 23°N–27°N, reaching Jianchuan in the north, Jinggu in the south, Zhenyuan in the east, and Tengchong and Ruili in the west (Fig. 1). All observation points are selected on relatively flat and hard open land, and the distance between stations is less than 10 km. The 2021 Yangbi Ms 6.4 Earthquake occurred just within the survey area, which provides a data basis for studying the lithospheric equilibrium state in the focal area.

The gravity survey is completed by simultaneous observation of two Burris gravimeters. The resolution of the gravimeter is 0.01 mGal (10^-5 ms^-2) and the observation accuracy is 0.02–0.03 mGal. In the process of measurement, the round-trip series observation in the form of 1-2-3...3-2-1 is adopted to reduce the transmission error of gravity observation. All gravity observations are closed on the same day to reduce the observation error caused by the nonlinear drift of the instrument. The observation results are adjusted by using the gravity adjustment software LGADJ (Liu et al., 1991) recommended by the China Earthquake Administration. Taking the absolute gravity values at Tengchong and Dali as constraints, the adjustment is calculated by the classical adjustment method. The overall accuracy of the adjustment results is about 0.02 mGal. The instruments used for elevation and coordinate observation are Leica GX1230 dual frequency GNSS receiver and LEIAX 1202 dual frequency antenna. The observation time of each station shall not be less than 40 min, and the sampling interval shall be 30 seconds to ensure the elevation accuracy of centimeter level. Finally, GAMIT\GLOBK (Herring et al., 2010) is used for baseline calculation, and the mean square error is about 5 cm. In the process of adjustment, the continuous GNSS data belonging to the nationwide network Crustal Movement Observation Network of China (Liang et al., 2013) in our survey network is used for position control.

Generally speaking, the accuracy of in-situ gravity and elevation data is high, but the number of observation stations is limited and the spatial resolution is low. Our gravity network consists of 318 stations in western Yunnan (Fig. 1). From the point of view of in-situ observation, it has been very dense, but the stations are mainly distributed along roads, so it is still necessary to further improve the resolution. In order to obtain high spatial resolution gravity data of the whole study area, it is necessary to fuse in-situ survey data and model data (Pavlis et al., 2012; Hwang and Parsons, 1995). In this study, the EIGEN-6C4 global gravity model (Förste et al., 2014) and SRTM15PLUS elevation model (https://www2.jpl.nasa.gov/srtm/, accessed on 5 May 2021) are used. EIGEN-6C4 model is a global gravity model released by German Research Center for Geoscience (GFZ) in 2014. It is expanded to 2 190 orders and is the latest model of EIGEN-6C series. The data sources of the model are LAGEOS, GRACE, GOCE and DTU data in-situ(Andersen, 2010). In the land area where the wavelength exceeds order 235, the model is basically consistent with EGM2008 (Pavlis et al., 2012).

Firstly, the normal field correction and elevation correction are carried out for the observed gravity value to obtain the free-air gravity anomaly. Then, the measured and modeled free-air gravity anomalies are fused by using the least square collocation method based on a remove-and-restore prediction algorithm (Hwang and Parsons, 1995). Specifically, the in-situ measurement of gravity is subtracted from the model gravity value at the corresponding position to obtain the systematic deviation between the measured value and the model, and then the systematic deviation of all positions in the study area is estimated by the least square collocation method, and then the systematic deviation is superimposed on the model value to obtain the fused data. The calculation formula (Hwang and Parsons, 1995) is

herein, $ \Delta {g}_{p} $ is the gravity anomaly of the point to be calculated, $ {C}_{pi} $ is the covariance of the system deviation between the point to be solved and the known point, $ {C}_{ij} $ is the covariance between known points, $ \Delta {g}_{oi}–\mathrm{ }\Delta {g}_{mi} $ is the systematic deviation of each observation point, $ \Delta {g}_{mp} $ is the model value of the point to be solved. The covariance between the point to be calculated and the known point can be estimated by using the relationship between the distance between the known points and the covariance. In this paper, it is assumed that the relationship between distance $ d $ and covariance $ c $ conforms to the linear function $ c\left(d\right)=B\times d+c\left(0\right). $ The most appropriate coefficient is selected by the least square fitting method, and the average residual of the fitting is less than 0.01 mGal.

Figure 2a shows the free-air gravity anomalies based on the fused data. It can be seen that the regional free-air gravity anomalies generally present the characteristics of strip distribution, which is significantly correlated with the distribution of main fault zones. From Figs. 2a and 2b we can see that the fusion value and model value of the free-air anomaly field are close in the overall distribution. Along the Nujiang fault (F1), Sudian fault (F4), Dayingjiang fault (F5), east section of Wanding-Anding fault (F8), Yongde fault (F15), and Nanting River fault (F16), they are the low values of free-air gravity anomaly. On the contrary, the values are high in Nujiang West Branch fault (F2), Tengchong Quaternary volcanic fault zone (F3), Longling-Ruili fault (F7), west section of Wanding-Anding fault (F8), Caojian fault (F9), south section of Lancang River fault (F10) and Lanping-Yunlong fault (F11). In general, there is a significant deviation between the observed and modeled free-air gravity anomalies (Fig. 2c), and the maximum difference is more than 100 mGal. Such a difference is due to the limited accuracy of the EIGEN-6C4 model in high degrees. The data of the first 235 degrees of the EIGEN-6C4 model are mainly from the observation results of GRACE and GOCE satellites, and the higher degree results are deduced from the residual terrestrial terrain gravity model (Fu et al., 2021). The high accuracy in-situ gravity data in China is not easy to available for the makers of the EIGEN-6C4 model, therefore accuracy of the EIGEN-6C4 is not high in part areas of our country. This difference not only reflects the limitation of the accuracy of the model gravity anomaly field, but also proves the importance of gravity observation in-situ. The maximum deviation is located between Lushui and Dali, close to the source area of the 2021 Yangbi Ms 6.4 Earthquake (Fig. 2c), which is basically consistent with the profile P1 in this study.

Figure  2.  Gravity anomalies around the 2021 Yangbi Ms 6.4 Earthquake. (a) Free-Air gravity anomalies based on fusion of modeled and observed gravity data; (b) modeled free-air gravity anomalies; (c) difference between the fused and modeled free-air gravity anomalies; (d) bouguer gravity anomalies. Black arrows show the magnitude and directions of the gradient of Bouguer gravity anomaly. Legend and fault are consistent with Fig. 1a.

DownLoad: Full-Size Img PowerPoint

Then, the high precision and high resolution topographic model SRTM15-PLUS (Tozer et al., 2019) is used to obtain the complete Bouguer gravity anomaly in the study area (as shown in Fig. 2c). For terrain correction, the grid terrain data with resolution 5″×5″, 25″×25″ and 50″×50″ are used respectively in the near field (0–2 km), middle field (2–20 km) and far field (20–167 km) to reduce the calculation time and ensure the calculation accuracy at the same time (Fu et al., 2014). In order to eliminate the boundary effect, the grid data used for terrain correction extend 0.5°, 1° and 1.5° outside the study area respectively. In areas with gentle terrain, the terrain correction is less than 5 mGal. However, in Hengduan Mountain area with large topographic relief, the correction is up to 20 mGal.

The elevation observation error of 1 cm leads to an uncertainty of 0.003 mGal in Bouguer gravity anomaly. Our elevation accuracy measured by GNSS at most stations is less than 3 cm, which brings an error within 0.01 mGal. The gravity observation error is about 0.02 mGal. The error of Bouguer correction is within 0.1 mGal. The error of data fusion is related to the terrain and the distance from the observation point (Wang et al., 2007). Speak specifically, in the area with flat terrain, the error of less than 60 km from the observation point is less than 5 mGal. In areas with drastic topographic changes, the error of less than 2 km from the observation point is less than 5 mGal, and the error of 2–5 km is less than 10 mGal. To sum up, the accuracy of complete Bouguer gravity anomaly at each observation point obtained in this study is less than 0.2 mGal, that in flat areas outside the observation point is about 5 mGal, and that in areas with severe topographic changes is about 10 mGal (Table 1).

Table  1.  Errors of Bouguer gravity anomalies in the study area

Distance from observation stations (km)	GNSS observation error (mGal)	Gravity observation error (mGal)	Bouguer correction error (mGal)	Data fusion error (mGal)	Total error (mGal)
0	0.01	0.02	0.1	0	< 0.2
0–2	0.01	0.02	0.1	< 5	5
2–5	0.01	0.02	0.1	< 10	10

 | Show Table

DownLoad: CSV

Figure 2d shows the Bouguer gravity anomaly field in the study area. In general, the Bouguer gravity anomaly varies from -120 to -360 mGal. Generally speaking, the anomaly value in the north is large but that in the south is small. The overall anomaly characteristics should be related to the fluctuation of the Moho. The Bouguer gravity anomaly around Ruili in the southwest of the study area is the smallest, indicating that the crustal thickness is the shallowest, and the material from the Tibetan Plateau mainly migrates to the southeast of the study area. The Bouguer gravity anomaly around the epicenter of the 2021 Yangbi Ms 6.4 Earthquake is about -260 mGal, and there are significant spatial changes in both NS and EW directions (Fig. 2d), that is, there are significant lateral changes in the density structure in the source area of the 2021 Yangbi Ms 6.4 Earthquake.

2. DENSITY STRUCTURE AROUND THE 2021 YANGBI Ms 6.4 EARTHQUAKE

The density structure of lithosphere can be inversed by using Bouguer gravity anomaly data, but gravity inversion has multiple solutions, and the reliability of its inversion results strongly depends on the reliability of the initial model. In this paper, the initial density model is determined with reference to the seismic wave velocity structure. For the area around the 2021 Yangbi Ms 6.4 Earthquake, based on the S-wave velocity structure obtained by Zheng et al. (2016) using the receiver function, we determine the initial density model by applying the empirical relationship between crustal density and S-wave velocity given by the following formula (Laske et al., 2013).

herein, $ \rho $ is the density of crust, $ {v}_{s} $ is the structure of S-wave velocity.

According to the gravity anomaly and initial crustal density structure model, we use the software named GMSYS that developed by NGA to inverse the two-dimensional crustal density structure in the study area. The inversion process is as follows: (1) select a profile to extract the complete Bouguer gravity anomaly; (2) determine the initial model of crustal density structure of the profile according to the previous tomography results and the velocity-density relationship (Eq. 2); (3) calculate the theoretical gravity anomaly of the initial model based on the two-dimensional polygon prism forward modeling method (Talwani et al., 1959); (4) compare the difference between the calculated and observed values. If it is less than a certain threshold, the final result will be output. Otherwise the nonlinear least squares method (Marquardt, 1963) will be used for inversion until the difference between the observed and simulated values is less than the threshold.

Based on the previous researches (Zhang et al., 2017; Deng et al., 2014), we set the threshold to 5 mGal, and then iterate repeatedly according to the initial density model given above to obtain the crustal density structure in which the difference between the observed and modeled gravity anomalies is less than the threshold. Figure 3 is the final inversion result of the crustal density structure of the two profiles passing through the source area of the 2021 Yangbi Ms 6.4 Earthquake. The specific locations of the two profiles are shown in the yellow line of Fig. 1. Figure 3 clearly shows the density spatial distribution of the Crust and upper Mantle of each profile. It can be seen that the crustal density structure in the study area has obvious vertical stratification characteristics. According to the inversion of crustal density structure, the matching degree between the gravity anomaly calculated by forward modeling and the observed value is an important standard for the reliability of the result. The difference between the observed and fitted values on the two profiles does not exceed 4 mGal (P1: 3.57 mGal; P2: 3.82 mGal), which proves the reliability of our inversion results.

Figure  3.  Crustal density structure profiles passing through the focal area of the 2021 Yangbi Ms 6.4 Earthquake. (a) Observed and inverted Bouguer gravity anomalies in near EW direction (Profile P1); (b) terrain of Profile P1; (c) density structure of Profile P1; (d) observed and inverted Bouguer gravity anomalies in near SN direction (Profile P2); (e) terrain of Profile P2; (f) density structure of Profile P2. Red Pentagram denotes the 2021 Yangbi Ms 6.4 Earthquake. Black line is the Moho, which is obtained using surface wave dispersion and receiver function (Zheng et al., 2018).

DownLoad: Full-Size Img PowerPoint

According to Fig. 3, the 2021 Yangbi Ms 6.4 Earthquake is located in an area with significant spatial variation of Bouguer gravity anomaly. The density structure of the source area shows significant lateral differences in both shallow and deep crust, especially in the EW direction (Fig. 3b). Such a result indicates that the crust around the great earthquake does not belong to a solid and stable tectonic unit but a relatively weak tectonic boundary, so the ability to bear crustal stress is relatively poor. This structural feature should be the deep structural reason for the frequent occurrence of moderate and strong earthquakes in this area. The lateral changes of the above density structures in the source area of the great earthquake are also shown in similar documents (Li et al., 2021; Zhao et al., 2021).

3. LITHOSPHERIC EQUILIBRIUM AROUND THE 2021 YANGBI Ms 6.4 EARTHQUAKE

According to the buoyancy principle, if the lithosphere is in complete equilibrium, that is, when the lithosphere is in hydrostatic equilibrium, the gravity of the lithosphere itself is equal to the buoyancy generated by its displaced mantle. The isostatic surface reflects the theoretical position of the Moho in the complete equilibrium state. That is, when the lithosphere is completely balanced, the isotatic surface coincides with the Moho (Fu et al., 2014). However, when there is a vertical external force acting on the lithosphere, such as dynamic effect or material loading, in order to maintain the equilibrium state, the isostatic surface and the Moho will offset, and the buoyancy generated by the residual materials (difference between the densities of crust and mantle) between the offsets is basically the same as the external force, but with opposite direction. In fact the above buoyancy is basically carried by the deflection of the lithosphere. The greater the difference between the isostatic surface and the Moho, the more unbalanced the region is, and the greater the isostatic additional force borne by the lithosphere (Fu et al., 2021; Wang et al., 2018; Gao et al., 2016).

The force borne by the lithosphere before and after lithospheric equilibrium can be expressed as

Herein, $ {\rho }_{\mathrm{M}} $ is the density of Mantle, taken as 3.30 g/cm³; g is taken as 9.8 m/s². $ {h}_{\mathrm{M}\mathrm{o}\mathrm{h}\mathrm{o}} $ and $ {h}_{\mathrm{i}\mathrm{s}\mathrm{o}} $ are respectively the depths of Moho and isostatic surface. $ {\rho }_{i} $ and $ {h}_{i} $ denote the density and thickness of each stratum. $ {P}_{vs} $ is the isostatic additional force borne by the lithosphere. After further derivation we get

We use the method proposed by She et al. (2017) to calculate the depth of the isotatic surface $ {h}_{\mathrm{i}\mathrm{s}\mathrm{o}} $. This method assumes that the earth􀆳s crust is incompressible, completely elastic and conforms to the principle of buoyancy. When calculating the position of each interface of multilayer density structure in equilibrium state in frequency domain, the disturbance of different interfaces will cause the deflection of other interfaces to compensate. $ {h}_{\mathrm{M}\mathrm{o}\mathrm{h}\mathrm{o}} $ takes the depth of the Moho in the inverted crustal density structure, as shown in Fig. 3. Then we can obtain the spatial distribution of the isostatic additional force borne by the lithosphere in western Yunnan through Eq. (4).

The isostatic additional force borne by the lithosphere in western Yunnan shows a gradual overall change trend in the NS direction (Fig. 4). Among them, the isostatic additional force in Jianchuan, Tongdian and other areas in the north is positive, and the maximum value is about 10–15 MPa. The direction of the isostatic additional force in this area is upward (away from the Geo-center). On the contrary, the isostatic additional force in the south is negative, and the maximum negative value is about -40 MPa, which is distributed in Lincang, Jinggu and other places. The isostatic additional force in the source area of the 2021 Yangbi Ms 6.4 Earthquake is small and generally changes near the zero isoline. The research of She et al. (2017) shows that earthquakes frequently occur near the zero isoline of isostatic additional force, which is confirmed again by the 2021 Yangbi Ms 6.4 Earthquake. Figure 5 shows the two profiles of the isostatic additional force across the source area of the 2021 Yangbi Ms 6.4 Earthquake. It can be seen that the lithosphere in the source area of the great event is basically in equilibrium, and the spatial change of the isostatic additional force is very gentle. The results show that the occurrence of the 2021 Yangbi Ms 6.4 Earthquake has little to do with the lithospheric equilibrium in the surrounding areas, which should be mainly caused by the differential movement in the horizontal direction. The above conclusions are consistent with the almost pure strike-slip focal mechanism of the 2021 Yangbi Ms 6.4 Earthquake (Zhang et al., 2021).

Figure  4.  Isotatic additional force borne by the lithosphere around the 2021 Yangbi Ms 6.4 Earthquake. Legend and fault are consistent with Fig. 1a.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Two profiles of the isotatic additional force borne by the lithosphere across the source area of the 2021 Yangbi Ms 6.4 Earthquake. (a) Profile in near EW direction (P1); (b) profile in near NS direction (P2). Red Pentagram denotes the location of the great earthquake.

DownLoad: Full-Size Img PowerPoint

4. ANISOTROPY AROUND THE 2021 YANGBI Ms 6.4 EARTHQUAKE

From 2002 to 2004, we deployed 63 broadband seismic observation stations in western Yunnan and carried out mobile seismic observation. Among the seismic records with epicentral distance of 85°–110° and magnitude greater than 5.0, we checked 9 720 three-component seismic records, and selected out more than 1 800 records with clear SKS seismic phases, involving 103 earthquakes. We use multichannel analysis (Chevrot, 2000) to obtain shear wave splitting parameters. Multichannel analysis is a simple and stable shear wave splitting analysis method. The shear wave splitting parameters are estimated by analyzing the relationship between the radial and tangential relative amplitudes of waveforms and the incident azimuth. Firstly, the splitting intensity related to the incident polarization direction is obtained by using the derivatives of the tangential and radial components of the seismic waveform, and then a sinusoidal curve is fitted according to the splitting intensity in each direction. The phase and amplitude of the sinusoidal curve are the fast wave polarization direction and delay time respectively. In this paper, the splitting intensity of each waveform record is calculated respectively, and then a weighted average value is obtained in the 10° azimuth window as the splitting intensity of the azimuth window, and the weighting factor is the reciprocal of each splitting intensity error.

After the above processing, we obtained the shear wave splitting parameters on 51 stations (red lines in Fig. 6). Detailed numerical results are shown in Table 2. The results show that the polarization direction of fast waves in the lithosphere in southwestern Yunnan is generally 100°–125° north by east, and the delay time of fast and slow waves is generally between 0.3 and 1.0 s. The fast wave polarization direction is generally consistent with the surface movement direction relative to the stable Eurasian Plate observed by GNSS in this area (Wang et al., 2001).

Figure  6.  XKS (SKS, SKKS and PKS seismic phases) shear wave splitting results in the area around the 2021 Yangbi Ms 6.4 Earthquake. Red lines are the results presented by this study. Blue lines are the results of Chang et al. (2015).

DownLoad: Full-Size Img PowerPoint

Table  2.  The SKS split parameters of 51 stations in the study area

Station	Lat. (°E)	Long. (°N)	Φ (°)	δt (s)	Number of event
007YY	23.390 0	101.860 0	80.8 ± 5.48	0.6 ± 0.55	7
009PD	24.060 0	101.950 0	97.5 ± 3.78	0.6 ± 0.36	5
010JX	23.686 1	101.637 8	79.4 ± 4.83	0.4 ± 0.26	5
011ES	24.180 0	102.390 0	166.0 ± 6.45	1.2 ± 1.02	6
013ZX	23.242 4	101.744 2	106.6 ± 5.26	0.4 ± 0.28	4
014TG	23.310 0	101.390 0	116.4 ± 8.05	0.6 ± 0.60	6
015DA	23.320 0	101.150 0	115.5 ± 9.24	0.7 ± 0.63	5
016MX	22.995 2	101.225 6	80.1 ± 4.63	0.7 ± 0.42	6
018GC	23.711 3	101.146 0	109.0 ± 6.48	0.7 ± 0.41	6
019AB	23.857 2	100.894 8	106.7 ± 7.15	0.7 ± 0.42	7
020FS	23.688 1	100.815 4	96.5 ± 3.48	0.6 ± 0.36	6
021JG	23.737 9	100.621 1	85.0 ± 2.21	0.9 ± 0.38	7
022JT	23.501 9	100.735 1	80.6 ± 7.24	0.7 ± 0.56	6
023YP	23.422 4	100.406 4	105.0 ± 4.46	0.7 ± 0.28	7
024MN	23.384 4	100.944 9	113.7 ± 7.77	0.7 ± 0.61	6
026LC	23.876 5	100.071 9	103.6 ± 4.05	0.7 ± 0.22	6
027MT	23.798 0	100.262 6	102.7 ± 4.00	0.6 ± 0.30	7
028QN	23.594 2	100.065 6	74.9 ± 7.06	0.7 ± 0.55	6
029SJ	23.500 7	99.821 8	101.9 ± 7.56	0.7 ± 0.43	7
030MY	24.007 4	99.823 6	103.6 ± 5.12	0.8 ± 0.29	5
032GM	23.546 0	99.392 5	104.3 ± 10.00	0.7 ± 0.48	7
033MJ	23.707 5	99.345 5	104.3 ± 12.22	1.3 ± 1.14	5
034NC	23.543 5	99.068 5	108.0 ± 11.54	0.4 ± 0.27	7
036YD	24.035 7	99.245 4	107.5 ± 6.23	0.7 ± 0.34	7
037BK	24.184 0	99.476 6	41.6 ± 18.34	0.7 ± 0.56	7
038YL	24.232 7	99.571 2	111.1 ± 6.40	0.8 ± 0.49	7
039SC	24.930 0	99.070 0	127.2 ± 13.71	0.6 ± 0.45	7
040SD	24.720 0	99.200 0	116.4 ± 11.16	0.6 ± 0.29	7
041YG	24.577 4	99.314 3	119.1 ± 9.50	0.4 ± 0.24	6
042QZ	24.503 8	98.809 9	127.9 ± 8.51	0.8 ± 0.47	8
043HN	24.353 0	98.913 0	124.2 ± 11.49	0.6 ± 0.34	9
044LJ	24.753 2	98.687 8	123.0 ± 11.38	0.4 ± 0.36	7
045MZ	25.210 0	98.470 0	101.2 ± 4.82	0.4 ± 0.27	7
046QS	25.230 0	98.590 0	116.4 ± 5.90	0.7 ± 0.33	7
047HQ	25.340 0	98.270 0	143.1 ± 7.18	0.5 ± 0.20	6
048SY	25.020 0	98.640 0	133.8 ± 6.40	1.1 ± 0.67	6
049TC	25.030 0	98.520 0	117.0 ± 7.85	0.6 ± 0.29	9
050HX	24.820 0	98.280 0	170.8 ± 6.11	0.5 ± 0.39	8
051TB	24.620 0	97.660 0	179.4 ± 1.01	0.5 ± 0.38	4
052YJ	24.690 0	97.950 0	130.3 ± 7.67	0.5 ± 0.17	8
053KC	24.990 0	97.800 0	142.0 ± 9.64	0.3 ± 0.19	8
055QJ	25.506 2	98.053 8	122.9 ± 5.88	0.5 ± 0.24	6
056JD	24.520 0	98.390 0	100.3 ± 5.73	0.3 ± 0.24	6
057LX	24.420 0	98.590 0	119.1 ± 24.15	0.3 ± 0.30	9
059DX	25.609 3	100.254 0	125.5 ± 18.11	1.0 ± 0.75	6
061DX	25.054 0	100.522 5	114.3 ± 7.92	0.7 ± 0.34	6
062DX	25.148 2	102.747 6	118.8 ± 11.14	0.3 ± 0.28	7
063DX	25.215 0	102.491 5	118.0 ± 9.55	0.4 ± 0.28	7
064DX	25.104 3	101.918 0	122.5 ± 15.02	0.5 ± 0.32	6
065DX	25.029 3	101.535 6	115.4 ± 12.06	0.7 ± 0.46	5
066DX	25.238 9	101.144 7	108.1 ± 6.81	1.0 ± 0.46	6
9999S	24.420 0	98.590 0	117.3 ± 7.16	0.6 ± 0.25	12

 | Show Table

DownLoad: CSV

Predecessors have done a lot of similar work in the study area and adjacent areas. Ruan and Wang (2002) studied the shear wave splitting phenomenon at 23 fixed stations in Yunnan. Chang et al. (2006) studied the shear wave splitting phenomenon at fixed stations and some mobile stations in Yunnan. The polarization direction of the fast wave is mostly north by east. The fast wave direction is basically consistent with the results of this paper, and the delay time is generally greater than that of Ruan et al. (2002). Flesch et al. (2005) studied the shear wave splitting of seven mobile stations and two fixed stations in the region, and its fast wave direction is between 87° and 105°, with an average of 95.7°. Chang et al. (2015) did polarization analysis based on the tele-seismic XKS (SKS, SKKS and PKS seismic phases) waveform data recorded by "ChinArray Exploration-Southern Segment of South-North Seismic Belt" and the broadband fixed stations of the China Seismological Network, and obtained the upper mantle anisotropy image in the southern section of the north-south seismic belt of China. In order to better reveal the anisotropic characteristics around the source area of the 2021 Yangbi Ms 6.4 Earthquake, the XKS shear wave splitting results given by Chang et al. (2015) are also added in Fig. 6 and expressed by blue lines. In general, the results of Chang et al. (2015) are in good coordination with the SKS shear wave splitting results given in this paper, which proves the correctness and reliability of our results to some extent. Our results are a good supplement to the results of Chang et al. (2015) and others, and improve the spatial resolution of shear wave splitting images in southwestern Yunnan to a certain extent.

According to Fig. 6, the lithospheric fast wave direction near the source of the 2021 Yangbi Ms 6.4 Earthquake is generally 97° north by east, and the fast and slow wave delay time is about 1.0 s. The shear wave splitting in the south of the earthquake has a strong regularity, showing obvious EW distribution characteristics, and has a gradual change in space. On the contrary, the distribution of shear wave splitting results is complex in the north of the great event. Both the dispersion of fast wave direction and fast and slow wave delay time are high. Based on the results of near-field small earthquakes, far-field earthquakes and background noise data, Gao et al. (2020) comprehensively gave the regional spatial distribution and vertical hierarchical distribution of seismic anisotropy in the southeastern margin of Tibetan Plateau. It was found that the azimuthal anisotropy presents the characteristics of south and north zoning, and the north-south boundary is near latitude 26°. The fast wave direction is approximately NS direction in the north, however, in the south it is approximately EW direction. The 2021 Yangbi Ms 6.4 Earthquake is located near such a boundary.

Through comparison, it is found that the fast wave direction in western Yunnan (Fig. 6) and the maximum change direction of Bouguer gravity anomaly gradient (Fig. 2d) are approximately vertical in general. The regional Bouguer gravity anomaly shown in Fig. 2d mainly reflects the information of the Moho. The results show that the crustal thickness in western Yunnan is generally deep in the north and shallow in the south. This situation is also shown on the Moho profile in nearly NS direction (Fig. 3f). Figure 4 shows that the lithospheric equilibrium in southwestern Yunnan is generally negative anomaly, indicating that the crust in this area tends to sink under the action of self-gravitation, or the deep mantle material has an obvious supporting effect on the crust, so that the crust is lifted to a certain extent. Therefore, we believe that in the process of escaping to the southeast, the material of the Tibetan Plateau is blocked by the shallow Moho (or high-density mantle material) in the front in southwestern Yunnan, so it disperses and escapes to the near EW direction, which may be the internal cause of the fast wave direction in the near EW direction in southwestern Yunnan, China (Fig. 6).

Both deep and shallow earthquakes should be the result of the accumulation of stress and strain in the whole crust. Revealing the deep structure and dynamic environment around the source is helpful to understand the preparation and occurrence process of the earthquake. That is why we studied the lithospheric equilibrium and seismic anisotropy around the 2021 Yangbi Ms 6.4 Earthquake based our in-situ observations.

5. CONCLUSION

From August to October 2020, we conducted in-situ gravity/GNSS joint observation at 318 stations in western Yunnan and obtained high-precision gravity data. The 2021 Yangbi Ms 6.4 Earthquake happened in our survey area, which provides a good basis for studying the equilibrium state of the focal area. We carried out a series of data processing, optimized the Bouguer and free-air gravity anomalies around the great earthquake, inversed the density structure of the lithosphere, and obtained the lithospheric equilibrium distribution around the focal area. The results show that the free-air gravity anomaly in western Yunnan presents a north-south strip distribution, which has a good correlation with the terrain. The Bouguer gravity anomaly varies from -120 mGal in the south to -360 mGal in the north, and about -260 mGal in the source area of the 2021 Yangbi Ms 6.4 Earthquake. Our inversion results show that the density structure in the source area of the great earthquake presents significant lateral differences, especially in the EW direction, indicating that the crust around the great earthquake does not belong to a solid and stable tectonic unit but a relatively weak tectonic boundary, and it is difficult to accumulate large stress. This feature should be one of the main reasons for the frequent occurrence of moderate and strong earthquakes in this area. The calculated isostatic additional force show that the lithosphere in the source area is basically in equilibrium, indicating that the occurrence of the great earthquake has little to do with the regional lithosphere equilibrium, but mainly related to the differential crustal movement in the horizontal direction that results from the southeast extrusion of Tibetan Plateau.

At the same time, according to the broadband mobile seismic observation carried out in southwestern Yunnan, we analyzed the teleseismic SKS phases recorded at 63 stations, obtained the shear wave splitting parameters at 51 stations, and improved the spatial resolution in the study area. The results show that the polarization direction of fast waves at most stations in southwestern Yunnan is between 100° and 125° north by East, and the delay time of fast and slow waves is between 0.3 and 1.0 s. There is no obvious change in the polarization direction of fast waves at stations on both sides of the Honghe fault. The overall fast wave polarization direction in the focal area of the 2021 Yangbi Ms 6.4 Earthquake is 97°, and the fast and slow wave delay time is about 1.0 s.

In general, the polarization direction of fast wave in western Yunnan is approximately perpendicular to the maximum gradient direction of regional Bouguer gravity anomaly that reflects the change of the Moho. Based on the comprehensive analysis of Bouguer gravity anomaly, lithospheric equilibrium and anisotropy, we believe that in the process of escaping to the southeast, the Tibetan Plateau material is blocked by the shallow Moho (or high-density mantle material) in southwestern Yunnan, and the crust is lifted to a certain extent, thus forming the characteristics of dispersion and escape in the near EW direction.