The Evolution of Stress and Strain around the Bayan Har Block in the Tibetan Plateau

With the constraint of GPS observation, the tectonic deformation of the Bayan Har block and its periphery faults is investigated based on an elastoplastic plane-stress finite element model. The results show that the elastic model cannot explain the current GPS observation in the Bayan Har block.When East Kunlun fault and Yushu-Xianshuihe fault are under plastic yield state or high strain localization, the calculated velocities fit well with the observation values. It indicates that most of the current shear deformations or strain localizations are absorbed by these two large strike-slip faults. In addition, if the recurrence intervals of large earthquakes are used to limit the relative yield strength of major faults, the order of entering the plastic yield state of the major faults around Bayan Har block is as follows. The first faults to enter the yield state are Yushu-Xianshuihe faults and the middle segment of East Kunlun faults. Then, Margaichaka-RolaKangri faults (Mani segment) and Heishibeihu faults would enter the yield state. The last faults to enter the yield state are the eastern segment of East Kunlun faults and Longmenshan faults, respectively.These results help us to understand the slip properties of faults around the southeastward moving Bayan Har block.


Introduction
The Tibetan Plateau is the main seismic active area in China.From the south to the north, it is composed of the Himalayan block, Lasha block, Qiangtang block, Hoh Xil-Bayan Har block, Qaidam block, and Qilian Shan block.The Bayan Har block is located in the northern Tibetan Plateau [1,2].Its boundary faults include the Kunlun fault in the north, Mani-Yushu-Xianshuihe fault in the south, and Longmen Shan in the east [3,4].A series of strong earthquakes have occurred around the Bayan Har block in the Tibetan Plateau since 1997, including Mani earthquake (Mw7.5 on November, 8, 1997 UTC) in Xinjiang province along the Mani fault [5], Kunlun Shan earthquake (Mw7.8 on November, 14, 2001 UTC) in Qinghai province along the Kunlun Shan fault [6][7][8], Yutian earthquakes (Mw7.1 on March, 20, 2008 UTC and Mw6.9 on February, 12, 2014 UTC) in Xinjiang province along the Altyn Tagh fault [9][10][11], Wenchuan earthquake (Mw7.9 on May, 12, 2008 UTC) in Sichuan province along Longmen Shan fault [12][13][14], Yushu earthquake (Mw6.9 on April, 13, 2010 UTC) in Qinghai province along Yushu-Xianshuihe fault [15,16], and Lushan earthquake (Mw6.6 on April 20, 2013 UTC) [17,18].These earthquakes have composed the main seismic activity zone in the Tibetan Plateau in the past decades years [3].The study of the Coulomb stress among strong earthquakes around the Bayan Har block suggested that it was limited to the fault activity caused by the strong earthquakes.The maximum static Coulomb stress change triggered by the large earthquakes was around 0.01 MPa [19,20].Focal mechanism of these earthquakes implied that the Bayan Har block has been moving southeastward (Figure 1).Yutian earthquake is caused by a normal fault, which implied extension setting [21,22].Wenchuan earthquake indicated the intense activation of Longmen Shan thrust-fold belt, which is a typical compressive tectonic environment [23][24][25].Kunlun Shan earthquake in the north of the block and Yushu earthquake in the south of the block suggested the strike slip of Kunlun Shan fault and Yushu-Xianshuihe fault [26][27][28].Thus, what is the relation between these strong earthquakes and the movement of Bayan Har block?Solving this issue will help us to understand the relation between the block movement and occurrence of strong earthquakes around it.
To describe the large-scale and long-term deformation, the continental lithosphere is regarded as a continuum, obeying a Newtonian or a power law rheology [29][30][31].However, elastic constitutive equations are usually used to study the short-term deformation of the crust or lithosphere [32,33].In plate interiors, deep faults always play an important role during the lithospheric deformation [34,35].For instance, the slip and rupture of faults are elastoplastic deformation process.The recurrence of a large earthquake gets to be more than 1000 years.The recurrence of Wenchuan earthquake is estimated to be 2000-6000 years [14].Thus, elastoplastic rheology with the continuous assumption is adopted in this study, which is effective and convenient to consider the effect of strain localization which is equivalent to the effect of faults [36][37][38][39][40]. Large faults always control the overall movement of the block.Previous studies suggested that large active faults also control the uplift and expansion of the Tibetan Plateau [35,41,42].The slip of faults generally obeys the plastic deformation.However, the strong blocks present rigid or elastic features.Thus, the overall movement of blocks and its peripheral faults should be treated differently in the numerical study.GPS observation is one of the most direct reflection for the short-term tectonic active [33].It can be used as the boundary conditions in our numerical model.
In this study, based on the GPS observation, we constructed a plane-stress numerical model to investigate the evolution of stress and strain around the Bayan Har block.In this process, we try to investigate the relation between the movement of block and occurrence of strong earthquakes.

Numerical Equation.
Faults activities present stick-slip characteristic [44].In this study, we give high yield strength for the blocks and low yield strength for the faults.Thus, the faults would get to the plastic state firstly under the same increasing stress.High strain rate would be localized along the faults, which causes large slip and stress drop.Corresponding to the seismic activity, large slip or decrease of stress implies the occurrence of earthquake.In this study, a plane-stress model is adopted to study the deformation of block and its surrounding faults.For the plane-stress, where V is Poisson's ratio,   is the normal stress with  direction, and   and   are two components of shear stresses with  direction.  is the normal strain with  direction, and   and   are two components of shear strains with  direction.
The constitutive equations in this study are the incremental relation between strain  and stress  [45].In the Cartesian coordinate system, the constitutive equations can be expressed as where The terms in (3) are expressed as [ e ] and [ p ] are the elastic and plastic stiff matrixes, respectively.[ ep ] is the elastoplastic stiff matrix. is the yield function. is the parameter indicating the loading history, which can be taken as the plastic work  p or the accumulated plastic strain  p .{ e } and { p } are the elastic and plastic strain increments, respectively.{} is the total strain increment.The same definition is used for the stress and the yield function.E and V are Young's modulus and Poisson's ratio, respectively.Superscript (e) indicates the elastic strain increment and subscript (p) indicates the plastic strain increment.
The yield function, or can be called loading function, satisfies Drucker-Prager yield criterion is used in this calculation: 1 is the first invariant of stress tensor and  2 is the second invariant of deviatoric stress tensor.
For the plane-stress state, where   ,   , and   present the three independent components of stress in the Cartesian coordinate system. and  denote the material parameters, which can be expressed as where  and  are cohesion and angle of internal friction, respectively.

Finite Element Mesh and Boundary Conditions.
The finite element model covers the whole Bayan Har block and its peripheral blocks and faults.The boundaries of blocks are regarded as faults, which have low yield strength.To facilitate the establishment of the numerical model, the faults around the Bayan Har block are approximately divided into 6 sections (Figure 2).The total number of quadrilateral elements and nodes are 32397 and 32959, respectively.The boundary conditions are imposed with the recent GPS observation [46][47][48].We assume that the GPS velocity is not changed within the calculation time.According to the incremental method for elastoplasticity, an annual displacement from the GPS observation is regarded as an incremental time step.
In the given calculation time, the final strain or stress is the accumulation of incremental strain or stress.In the elastic state, the strain and stress are having a linear relationship.However, if the deformation gets to the plastic state, the accumulation of strain and stress presents strong nonlinear relationship. of GPS, the calculation result reflects the deformation of blocks and faults during this time.Active faults present plastic deformation when they slip.But, during the interval of slip, the locking faults present elastic state.Thus, the faults along which large earthquakes occurred can be considered as being under the plastic state.Firstly, we try to get the present strain and stress for the Bayan Har block and its peripheral blocks and faults under the constraint of GPS observation.Four cases are designed to investigate the present strain and stress for the Bayan Har block and its peripheral blocks and faults (Table 1).In case 1, both the blocks and faults are under the elastic state, and none of them gets to the plastic yield.In case 2, F5 is given low yield strength, which would get to plastic state firstly.In case 3, F2 and F5 are given low yield strength; other block and faults are given high yield strength.In case 4, F2, F4, and F5 are given low yield strength.We will compare the calculated results with GPS observation in the interior of the model for each case.Root mean square error (RMSE) would be calculated and the minimum error model would be selected as the optimization model to calculate the present stress and strain for Bayan Har block.

Model
Secondly, we try to get the evolution of strain and stress around Bayan Har block in the future.It is difficult to estimate the actual yield strength.However, under specific boundary conditions and the geometry of fault, the yield strength of faults is corresponding to the earthquake recurrence intervals.Therefore, to investigate the evolution of strain and stress around the Bayan Har block, we give the yield strength for different faults according to the recurrence intervals of earthquakes (case 5) [14,[50][51][52].All the blocks are under the elastic state in the calculation.
The elastic modulus and Poisson's ratio are obtained by the seismic velocity and crustal density [53][54][55].Plastic parameters are more complicated.The internal friction angle  and cohesion  have close relation with lithology, confining

Present Deformation around Bayan Har Block.
In case 1, we calculate the velocity in one incremental displacement because velocity is the same in every incremental step under elastic state.For other cases, incremental displacement is loaded on the boundary of the model until the faults which have low yield strength get to plastic state.Root mean square errors (RMSE) [58,59] have been calculated for each case (Table 2).By comparing the calculated velocity and GPS observation, we select the optimization case (Figure 3).For case 1, the slip among the blocks is relatively small because there are not plastic faults among them.Actually, large earthquakes occurred along the large strike-slip faults, such as Kunlun Shan fault (F2) and Yushu-Xianshuihe fault (F5), which caused the movement of the blocks.Thus, the fitting is not very well along the large strike-slip faults (Figure 3(a)).
Case 3 has the minimum fitting error.Thus, we select it as the optimization model to calculate the present stress and strain rate for Bayan Har block (Figure 4).In this case, the Kunlun Shan fault and Yushu-Xianshuihe fault have high accumulation strain (Figure 4(b)).Large displacement along these strike-slip faults adjusts the substance movement from the strong Indo-Asian collision [58].High strain rate along the faults allows the crust to move southeastward in the eastern Tibetan Plateau [69].In the block interior, the strain rate is relatively uniform and small.Corresponding to the high strain rate, the accumulation stress along the strike-slip faults is relatively low because the plastic yield causes stress drop.The stress directions show that the extension is dominant in the northwest of Bayan Har block (Figure 4(c)), for instance, nearly the westeast extension in the location of 2008 Yutian earthquake.However, compression is dominant along the Longmen Shan fault.This results are consistent with the focal mechanism (Figure 4(d)).

The Possible Evolution of Stress and Strain around
Bayan Har Block.Although there are many factors that influence the deformation results, the evolution of stress and strain for the mainly faults under the present velocity constraint can be investigated with specific assumption.In our model, elastic state implies that the fault is locked and plastic state indicates that the fault is rupture or slip fault.Thus, with the incremental displacement constraint, the faults with different yield strength would enter plastic state gradually.This process would help us learn which fault is easy for the stress accumulation and which is not.Thus, in case 5, we imposed the yield strength for different faults according to the recurrence intervals of the large earthquakes.The values given in the case may not be equal to the real values.But the order by which the faults get to the plastic yield helps us to understand the deformation process.Based on the surface rupture of earthquakes and strain rate, the recurrence intervals for the  ≥ 6.7 earthquakes along Xianshuihe fault are estimated to be 100-150 years [52].The recurrence intervals of large earthquakes along the East Kunlun Shan fault, such as 2001 Kusai Lake Ms8.1 earthquake, are 250-350 years [27,50].For the Altyn Tagh fault, the recurrence intervals of earthquakes, such as 2008 Kangxiwa Mw7.4 earthquake, are 370-500 years [50], while, for the Longmen Shan fault, which is characterized slow strain accumulation, such as 2008 Wenchuan earthquake, the recurrence intervals are 3000-6000 years [14,51].With these studies, we give a relative values of the yield strength for different faults (case 5 in Table 1).With the constraint of boundary velocity, we will investigate the evolution of stress and strain around the Bayan Har block.But the initial stress is 0 for all the blocks and faults.The order of fault rupture or slip only has relation with the distribution of faults, boundary conditions, and plastic yield strength.
We calculate the effective strain (Figure 5) and effective stress (Figure 6) for case 5 as follows: When the fault enters the plastic state, its effective strain is very large, indicating large slip along it (Figure 5).However, the effective stress is decreased, indicating that the faults release the accumulative stress (Figure 6).Stress is accumulated at the end the rupture fault, which caused the strain and stress increase along the surrounding faults.The process of accumulation and transfer of stress lasts until all the faults enter the plastic state.Because of the large yield strength, the blocks remain in the elastic state and have little deformation.The evolution of stress and strain shows that the first fault to get to the plastic state is the Kunlun Shan fault (F2) and Yushu-Xianshuihe fault (F5) because they have low yield strength.Then, the Heishibeihu fault (F1) and Margaichaka-RolaKangri faults (F6) get to plastic state.The last yield faults are the East Kunlun Shan fault (F4) and Longmen Shan fault (F3).This order has relation with not only the plastic strength, but also the distribution of faults and boundary conditions.The movement of block depends on the state of surrounding faults.It indicates that the overall movement of Bayan Har block is shown sectionally and periodically.

Discussion and Conclusions
In this study, the values of yield strength imposed on the blocks and faults are not equivalent to the realistic values because we cannot obtain the actual value of the faults at the present time.Therefore, if we change the values of them, the last time of entering the plastic state may change accordingly.However, the order of entering the plastic yield state for different major faults around Bayan Har block is almost the same.Actually, slip and locked state are alternatives to the activation of faults.We have not considered the stick-slip of the faults.More complex factors should be considered in the 3D model, such as heterogeneity, geometric distribution of fault.But the evolution of stress and strain in this 2D model help us to understand the process of stress and strain accumulation for different faults.
In this study, a 2D plane-stress numerical model was used to investigate the present deformation of Bayan Har block and its possible evolution of stress and strain.Primary conclusions include the following.figures are generated by the Generic Mapping Tools (GMT) package [70].

Figure 1 :
Figure 1: Tectonic setting and distribution of recent large earthquakes around Bayan Har block.Mechanisms of earthquakes come from Harvard CMT catalogue (http://www.globalcmt.org/CMTsearch.html).The blue arrows show the movement of block in the Tibetan Plateau according to GPS observation.Red solid lines indicate the major active faults around Bayan Har block [3, 12, 43].F1: Heishibeihu fault; F2: the middle segment of East Kunlun fault; F3: eastern segment of East Kunlun fault (Maqu-Heye segment); F4: Longmen Shan fault; F5: Yushu-Xianshuihe-Anninghe fault; F6: Margaichaka-RolaKangri faults.The black dashed lines show the boundary of different block in the Tibetan Plateau.The red dashed lines show the area of the finite element model (Figure 2).

Figure 2 :
Figure 2: Finite element model used in this study.Different colors in the model show the blocks and faults around Bayan Har block.B1: Qaidam-Qilian Shan block; B2: Bayan Har block; B3: Qiangtang block; B4: Sichuan basin.F1-F6 indicate the major faults (Figure1).The enlarge region shows the grid mesh of finite element model.

Figure 3 :
Figure 3: Velocity comparison between calculation and GPS observation for cases 1-4.The calculative velocity in each figure shows the results after 1020 year (1020 time steps) in the elastoplastic state.

Figure 4 :
Figure 4: The F2 and F5 faults are under the yield state, the distribution of velocity (a), strain rate (b), and cumulative stress (c) for case 3. Mechanism of earthquakes (d) comes from Harvard CMT catalogue (Mw > 4.0, 1976-2011) (http://www.globalcmt.org/CMTsearch.html).Red and blue lines indicate the tension and pressure axes from the mechanism of earthquake, respectively.

Figure 5 :
Figure 5: The effect of periphery faults on the accumulative strain of Bayan Har block.The left lower diagram in each figure indicates the state of periphery faults.The red color marked on the fault shows that it is under the yield state. (a)-(d) show the accumulative strain distribution after loading with GPS observation for 264 years, 274 years, 280 years, and 542 years, respectively.The parameters used in this calculation are shown in Table 2 (case 5).

Figure 6 :
Figure 6: The effect of periphery faults on the accumulative stress of Bayan Har block.The left lower diagram in each figure indicates the state of periphery faults.The red color marked on the fault shows that it is under the yield state. (a)-(d) show the accumulative stress distribution after loading with GPS observation for 264 years, 274 years, 280 years, and 542 years, respectively.The parameters used in this calculation are shown in Table 2 (case 5).

Table 1 :
Cases and material parameters.
Design and Parameters.The GPS data used in this study is from 1991 to 2005 [46, 48, 49].Based on the constraint

Table 2 :
[56,57]r each model.In this study, we select them based on the low strain rate test results for marble and sandstone[56,57].