Comparative analysis on seismicity and stress triggering of strong earthquakes sequence in Central Tibet

Central Tibet constitutes part of the central part of the Qinghai–Tibetan Plateau, which is one of the highest seismically active areas in China. This paper discusses the regularity of seismic activity in this area. Based on a stratified viscoelastic earth model, we calculated the Coulomb stress changes imparted from four strong earthquakes (M ≥ 6.3) along the Bengco–southeastern piedmont of the Nyainqentanglha mountain fault zone in this region. The result shows that the study area may be entering a new active period from 2020. There was a trigger between the strong earthquakes (M ≥ 6.3) on the Bengco fault zone. The post-seismic viscous relaxation effect of a strong earthquake had a significant impact on subsequent earthquakes (M ≥ 6.3). In the next 100 years, the Coulomb stress loading is more than 1.0 MPa in the northwest section of the Bengco fault and the central part of Nimu segment of the southeastern piedmont of the Nyainqentanglha mountain fault. Thence, strong earthquakes may occur along these fault segments. The maximum magnitude of the earthquake could be M6.7 in the next 100 years.


Introduction
Central Tibet has strong tectonic activity and frequent strong earthquakes (Bgure 1). The coseismic dislocation of an earthquake will change the stress state of the surrounding faults or adjacent areas, resulting in Coulomb stress changes (Xie et al. 2010). In recent years, the spatiotemporal relationship between Coulomb stress changes caused by earthquakes and subsequent earthquakes has attracted widespread attention from seismologists (Harris 1998;Gomberg and Felzer 2008;Yin et al. 2018;Yang et al. 2019;Shan et al. 2020). Studies have shown that small Coulomb stress changes (over a threshold of 0.01 MPa) may trigger earthquakes (Harris 1998;Stein 1999), leading to a change of future seismic activity in the region.
The crustal deformation of the study area mainly develops a series of near south-north tectonic belts (Jiao et al. 2015), and northwest and northeast tectonic belts. The Yangbajain-Damxung basin is located in the north-central part of the Yadong-Gulu Rift Valley system, and is 120 km long and 10-25 km wide (Qiao et al. 2009).
The southeastern piedmont fault strike of Nyainqentanglha mountain (SPFNM) generally shows a north-east trend, which is mainly distributed in the northwest of the basin and partly in the basin and east of the basin. The boundary fault zone of the basin is very active since the Quaternary, and it is a strong seismic activity zone in the interior of the Qinghai Tibet Plateau (Wu et al. 1992). The epicenter-concentrated zone mainly distributes on the active faults in the NW, NE, and nearly N-S directions. East-west sutures (such as the Yarlung Zangbo suture) have no obvious epicenter-dense zones (Wu et al. 1992).
In this paper, the impact of the four earthquakes on the surrounding area will be analyzed by calculating the Coulomb stress changes after each earthquake. This study will be of great scientiBc significance for understanding seismic activity in the region and estimating future seismic hazards of fault zones.

Seismic data and integrity analysis
The seismic data are mainly acquired through citing and checking the Catalogue of Chinese Historical Strong Earthquakes (23rd century BC-1911 AD) (China Earthquake Administration 1995), Catalogue of Chinese Modern Earthquakes  AD) (China Earthquake Administration 1999), Catalogue of Chinese Earthquakes (1831 BC-1969 AD) (Gu 1983), the catalogue of earthquakes of China Earthquake Network Center (CENC 2020), and the catalogue of recent and historic world earthquakes provided by the USGS (USGS 2020).
This study area has very few historical records of earthquakes, among which the early seismic data are missing, that of the earthquakes lower than M5 are missing, and the record of the destructive earthquakes originating since 1411 AD. With the gradual establishment of the seismic network, since 1950, there have been instrument monitoring data for strong earthquakes in the study area. Based on this estimate, the missed records for earthquakes above M5 are very low. Judging from the seismic records (China Earthquake Administration 1995, 1999, the earthquake records of M C 5 of this area are basically complete since 1950, and the earthquake records of M C 2 of this area are basically complete since 1970. The seismic data could be used as the reference basis for the study of the seismic activity and the estimation of their parameters in this area.

Characteristics of earthquake temporal distribution
According to available historical seismic data (China Earthquake Administration 1995, 1999, the earliest earthquake to occur in the area was the Duilongdeqing earthquake with a magnitude of M6 in 1264. Figure 2(a) shows the magnitude-time (M-T) diagram of M C 4.7 earthquakes in Central Tibet since 1200. Between 1200 and 1900, only nine earthquakes were recorded, indicating that the seismic data were lacking. However, the Damxung M8 earthquake record in 1411 was more detailed (Wu et al. 1992). It is clearer from Bgure 2(b) that the area has experienced two seismically active periods since 1900, with the Brst active period from 1920 to 1960 and the second active period from 1970 to the present. The study area is at the end of the second active period. From a statistically conservative point of view, at present, the trend of seismic activity in the study area in the next hundred years should be estimated based on the average seismic activity level of the study area.

Statistics of a-value and b-value of seismicity parameters
The Gutenberg-Richter magnitude-frequency relationship is equation (1) (Gutenberg and Richter 1942): where N is the number of cumulative earthquakes with magnitude greater than M, and a and b are statistical constants, where a is the intercept in the magnitude-frequency relationship, which represents the total number of earthquakes in the area or the severity of the frequency of seismic activity. The b value is the slope in the magnitude-frequency relationship, which represents the ratio of different magnitude earthquake frequencies in an area. It is obtained from the actual data statistics. According to the seismic data integrity analysis and time distribution characteristics of the study area, the cumulative annual average incidence (V) of earthquakes in different periods of the study area was calculated (V = cumulative number of earthquakes/statistical time period in years). The calculation results are shown in table 1.
The least squares method is used to Bt the data in each period to form a variety of statistical schemes. After comparative analysis, the best Btting scheme is Bnally selected, then a and b values of the study area are statistically obtained. a = 4.240, and its standard error is 0.038. b = 0.530, its standard error is 0.007, R is the correlation coefBcient, SD is the sum of the squared residuals, and the Bt is shown in Bgure 3.

Research method
According to the stress triggering theory of an earthquake, the stress accumulated on the fault plane is released when the earthquake occurs. Part of the stress will be transferred to the area around the seismic rupture surface, which will change the state of the stress Beld and further aAect the seismicity in the future (Stein 2003). Using elastic dislocation theory, we can obtain the Coulomb stress Beld changes caused by earthquakes in the study area.
The expression of the static Coulomb stress change is equation (2): where DCFS is the Coulomb stress changes on the fault plane, Ds s is the shear stress variation on the fault plane, l 0 is the eAective friction coefBcient of the fault, Dr n is the normal stress variation on the fault plane, and the tensile is positive. In this study, l 0 = 0.4 (Stein et al. 1992;King et al. 1994;Verdecchia and Carena 2015). This calculation uses a PSGRN/PSCMP calculation program (Wang et al. 2006) that is based on a layered gravity viscoelastic Earth model, which can be used to study co-seismic and post-seismic Coulomb stress Beld changes caused by strong earthquakes.

Lithosphere layered structure parameters
Based on the research results of the P-wave velocity and crustal density structure model below the Lhasa block (Christensen and Mooney 1995; Table 1 Note. N is the number of times (frequency). V is the cumulative annual average incidence of earthquakes in different periods. Zhao et al. 2001;Haines et al. 2003), the average crustal velocity ratio V P = 1.73V S model , the crustal velocity, and density models of the study area were determined. In this paper, the Maxwell volume rheological model (Wang et al. 2006) is used to simulate the viscoelastic eAects of the lower crust and mantle over decades to hundreds of years. Due to the high crust thickness and high temperature on the Tibetan Plateau, the lower crust viscosity coefBcient is low. The crustal viscosity coefBcient is 1.0 9 10 20 PaÁs. The mantle viscosity coefBcient is one order of magnitude higher than that of the lower crust, namely 1.0 9 10 21 PaÁs. The middle and upper crust is a completely elastic medium, in which the viscosity coefBcient is inBnite (Shi and Cao 2008;Yin et al. 2018). Table 2 shows the parameters of the multilayer crustal model.

Seismic dislocation model
The dislocation model of the Damxung earthquake M8.0 earthquake in 1411, the Bengco M7.7 earthquake in 1951, and the Gulu earthquake M7.4 in 1952 refer to the fault fracture model proposed by Yang et al. (2019). The strike, dip, rake, and maximum slip of the seismogenic fault of the Damxung earthquake, M6.3, in 2008 are 183.3°, 49.5°, À115°, and 1.30 m, respectively (Gao et al. 2010), and the length is 12 km and the width is 11 km (Qiao et al. 2009). Table 3 shows the data of the four strong earthquake dislocation models.

Stress evolution and earthquake triggering
According to the parameters listed in table 3, we calculated the Coulomb stress changes on the fault plane of the four earthquakes caused by the postseismic dislocation of each earthquake in chronological order (Bgure 4). As illustrated in Bgure 4, the red and white beach ball represents the focal mechanism of the receiver fault in each subgraph. The gray and white beach ball represents the focal mechanism of the seismogenic fault.  post-seismic viscoelastic relaxation of the 1411 M8.0 earthquake and 1951 M7.7 earthquake, the computed depth is 25.0 km, and it is located in the region of increased Coulomb stress. (c) Coulomb stress changes on the 2008 M6.3 earthquake's fault plane associated with the post-seismic viscoelastic relaxation of prior earthquakes before the 2008 M6.3 earthquake, the computed depth is 9.5 km, and it is located in the region of increased Coulomb stress. When calculating Coulomb stress, it is necessary to deBne seismogenic fault and receiver fault. The seismogenic fault is deBned as 'source fault' and the surrounding given fault to be studied is 'receiver fault' (King et al. 1994). The inCuence time of Coulomb stress changes after the earthquake that can last for more than 700 years. Shi et al. (2020) used the Burgers rheological model, which is more in line with the long-term deformation and short-term post-earthquake deformation, to simulate the co-seismic and post-seismic Coulomb stress evolution caused by earthquakes above M6.5 in North China since 1303. In the eastern Qinghai-Tibet Plateau, Xu et al.   Plateau. Therefore, the continuous inCuence of Coulomb stress changes may aAect the time until now. Figure 4(a) shows that the Bengco earthquake, M7.7, in 1951 was located in the area where the post-seismic Coulomb failure stress was significantly increased due to the Damxung earthquake, M8.0 in 1411. The Coulomb stress at the center of the rupture increased 0.143 MPa. Figure 4(b) shows that the Gulu earthquake, M7.4 in 1952, was also located in the area where the post-seismic Coulomb failure stress was significantly increased due to the Damxung earthquake, M8.0, in 1411 and the Bengco earthquake, M7.7, in 1951. The Coulomb stress at the center of the rupture increased by 0.136 MPa. Figure 4(c) shows that the Damxung earthquake, M6.3, in 2008 was also located in the area where the post-seismic Coulomb failure stress was slightly increased due to prior earthquakes before the 2008 earthquake. The Coulomb stress at the center of the rupture increased by 0.279 MPa.
In order to calculate the eAect of four strong earthquakes on a fault zone, the parameters of the receiver fault must be known. The main fault sliding properties and parameters in the study area are shown in table 4. Figure 5 shows the estimated, calculated distribution of Coulomb stress changes on the three fault zones after the Damxung earthquake, M6.3, in 2008 due to the four strong earthquakes. Based on the combined eAects of four earthquakes, this paper analyzes estimated future Coulomb stress changes (from 2033 to 2108 AD) on each fault. The east and west ends of the Yarlung Zangbo fault (F1) will gradually change from stress loading (0-0.01 MPa, 2033-2083 AD) to stress unloading (À0.01 to 0 MPa, 2108 AD). The Coulomb stress in the south of F2-1 will continuously load, which will lead to an increase as yet more than 1.0 MPa near the fault of Damxung earthquake, M6.3, in 2008, but the northern end of F2-1 is under Coulomb stress unloading, and the DCFS can reach À1.0 MPa. F2-2 will be under stress unloading as a whole (the DCFS can reachÀ1.0 MPa), and only the northeast end of F2-2 will be under stress loading (the DCFS can reach 0.1 MPa). The north and south of F2-3 will be under stress unloading as a whole (the DCFS can reach À0.1 * À1.0 MPa), and only the middle of F2-3 will be under stress loading (the DCFS can reach 0.1 MPa). The northwest of F3 will be under stress loading as a whole (the DCFS can reach 1.0 MPa), and the southeast of F3 will be under stress loading (the DCFS can reach À1.0 MPa). Chen et al. (2001) believed that there had been Bve seismically active periods and four seismically quiet periods in China since 1900. The distribution of seismicity from 1900 through the present is shown in table 5. The b-values in the seismically active periods were obviously lower than the bvalues in the quiet periods between earthquakes. The Bfth round of seismic activity in China from 1988 may not have ended as yet.

Seismic activity characteristics and trends
According to the seismicity in the study area, the author preliminarily identiBed the seismic active and quiet periods (interval period). The result was similar to that of Chen et al. (2001) with shortterm quiet periods between the active periods. The division of seismically active and quiet periods in this paper may be able to reCect the characteristics of seismicity in this area.
From the perspective of seismicity, it is mainly divided into active periods and quiet periods. Based on the statistical analysis of the duration, frequency, and energy of the active periods, the development trend of the seismically active period can be determined by analyzing the b value of active periods and stationary periods (Chen et al. 2001;Liu et al. 2015).
The a-value describes a measure of the overall seismicity of the earthquake, which represents the total number of earthquakes or the severity of the frequency of seismic activity in the area. The b value reCects the regional crustal rupture strength and stress state. A low b value can be used as an indicator of high eAective shear stress to determine strong earthquake danger zones. The value of a in the study area is high (a = 4.240), which indicates that the total number of earthquakes in the area is high or the frequency of seismic activity is high (Xie et al. 2019). The b value is relatively low (b = 0.530), meaning that the study area is currently at a relatively high level of stress accumulation. It is more likely that moderate and strong earthquakes will occur in the future (Xie et al. 2019). It can be seen from Bgure 2 that the Brst active period was 1920-1960 AD, and the b value is 0.518. The second active period was 1970-2010 AD, and the b value is 0.544. Each active period is around 40 yrs, and the interval period is at about 10 yrs. Starting from 2022, the study area may enter a new active period.

Analysis of seismogenic fault
Since 1411     Year 1900-19121912-19201920-19381938-19461946-19571957-19661966-19761976-19871988-(M7.4), in 1952, and Damxung earthquake (M6.3), in 2008 in the fault zone occurred under Coulomb stress loading due to a series of prior strong earthquakes. There was a trigger between the strong earthquakes on the fault zone and the increase of Coulomb stress generated after large earthquakes was conducive to the occurrence of subsequent earthquakes. According to the PSGRN/PSCMP Coulomb stress calculation program (Wang et al. 2006), the variation of Coulomb stress on each receiving fault at different times after the four strong earthquakes can be calculated. In this study, two points on the F2-1 fault are used as examples. Table 5 shows the estimated values of the Coulomb stress changes at two points on the F2-1 fault plane in the next 100 yrs. One hundred years after the Damxung earthquake (M6.3), in 2008, the Coulomb stress of the Yarlung Zangbo fault will unload gradually. The overall Coulomb stress of F2-1 is in the loading stage, and the central Coulomb stress loading is more than 1.0 MPa at Point 2. Therefore, earthquakes will occur in this segment. F2-2 is in the stage of Coulomb stress unloading, and the possibility of a future earthquake is little. Segment F2-3 is in the stage of Coulomb stress unloading, and the possibility of large earthquakes is also relatively little. The Coulomb stress in the northwest section of F3 is in the loading stage, and the loading amount of Coulomb stress is more than 1.0 MPa. So, it is estimated that strong earthquakes will occur in the future along this segment.
From table 6, as also observed by King et al. (1994), l 0 controls mainly the magnitude of DCFS, The characteristic two-point Coulomb change value is shown as an example. Figure 6. The map of different probabilities of exceedance-magnitude (the purple square spot is the likely magnitude 10% at 100 years).
rather than the overall pattern of stress-loading lobes and stress shadows. Based on the considerations above, we carried out all our calculations with a value of eAective friction of 0.4.

Seismic analysis on different time scales
The earthquake catalog in this study area was incomplete, which may have a certain impact on the analysis of seismic activity characteristics. In view of the incompleteness of the earthquake catalog, the evaluation of its earthquake disaster parameters can be carried out by the Kijko method (Kijko and Sellevoll 1992). This article introduces the calculation method of Kijko and Sellevoll (1989), Kijko and Sellevoll (1992), Kijko et al. (2016) to analyze the seismic activity of the study area in order to obtain a judgment of the possibility of future earthquakes in the region. Figure 6 shows the possibility of earthquakes of different magnitudes under different probabilities of exceedance. Figure 7 shows the return period [years] of  earthquakes of different magnitudes. Figure 8 shows the magnitude of earthquakes under different exceedance probabilities. It can be seen from Bgures 6-8, when the probability of exceedance in 100 yrs is 10%, an earthquake with a magnitude of *7.8 may occur (the purple square spot in Bgure 6). The return year of an M8 earthquake is at least 700 years. The maximum magnitude of earthquake in this area shall not exceed M9. The maximum estimated magnitude of the earthquake will be *M6.7 in the next 100 yrs.

Conclusion
Through the research and analysis of this paper, the following conclusions and understandings are drawn from the present study: • There was a trigger between strong earthquakes on the Bengco-southeastern piedmont of Nyainqentanglha mountain (B-SPNM) fault zone, and the increase of Coulomb stress generated after a large earthquake was conducive to the occurrence of subsequent earthquakes. • The Coulomb stress in the northwest section of F3 and the middle segment of F2-1 is in the loading stage, and the central Coulomb stress loading is more than 1.0 MPa, so strong earthquakes will occur in the future along these two segments. • The seismicity level of the study area is high, and it is currently at a relatively high level of stress accumulation. It is more likely that moderate and strong earthquakes will occur in the future. The interval period is approximately 10 yrs. Starting from 2022, the study area may enter a new active period. The maximum estimated magnitude of the earthquake will be M6.7 in the next 100 yrs.