Numerical simulation of fault deformation and seismic activity in the southern Qinghai–Tibet Plateau

The moderate-strong (M ≥ 5, since 1976; or M ≥ 6.5, overall) earthquakes of the north to south trending (NS-trending) faults and their vicinity in the southern Qinghai–Tibet Plateau exhibit prominent spatial concentration distribution characteristics. However, research on the southern Tibet faults is limited. Hence, this study used a two-dimensional viscoelastic finite element model to simulate crustal movement in southern Tibet based on 1991–2015 GPS velocity data. The current deformation field, tectonic stress field, fault slip, and stress accumulation rates distribution were obtained to analyze the relationship between fault activity and earthquake distribution. The results revealed that the simultaneous effected by of the NE-trending compression and uneven EW-trending tension on the study area. Crustal deformation exhibited simultaneous NS-trending compression and EW-trending stretching. The EW- and NWW-trending faults and the NS-trending faults differed in their mechanical properties and movement modes. The NS faults were primarily subjected to extensional stress with normal motion. The remarkable heterogeneity of the fault segments influenced the distribution of moderate-strong earthquakes and types of earthquake ruptures. The concentrated distribution of moderate-strong earthquakes on the NS-trending normal fault and its vicinity depended largely on the high slip rate, strong tensile stress, and geometric strike of the fault segment. This study aids in understanding the heterogeneity in normal fault activity in southern Tibet and is the basis for seismic hazard assessment.

margin, north to south (NS) normal faults and conjugate strike-slip faults are distributed in the interior, and large strike-slip faults in the near east to west (EW) direction are dominant on the northern margin. This plateau was strongly affected by the collision of the Indo-Eurasian plates and has experienced NS compression, rapid uplift, and EW extensional deformation since the Miocene epoch (Amoji, 1986;Xu et al., 2005).
The southern Qinghai-Tibet Plateau exhibits highfrequency and high-intensity seismic activity with a relatively spatial concentrated distribution (Cui et al. 1990;Zheng, 1995;Qiu et al., 2019). Moderate-strong earthquakes in this area are generally distributed along faults . Specially, the orthogonal fault system in the plateau primarily controls shallow intraplate earthquakes (Armijo et al., 1989;Han, 1987;Wu et al., 1994;Li, 2005). Since 1976, earthquakes with a magnitude (M) ≥ 5.0 have mainly been concentrated on a series of normal fault segments ( Fig. 1), such as the southern segment of the PalongCo fault (PLF), northern segment of the DangreYongcuo fault (DRYCF), and middle of the Nyainqentanggula-Yadong-Gulu fault (NQ-TGF) and Jiagang-Dingjie faults (JGF). Most strong earthquakes (M ≥ 6.5) have occurred in the above-mentioned segments, such as the 1411 Dangxiong M8, 1901 Nimu M6.8, and2008 Dangxiong M6.6 earthquakes, which occurred in and around the middle of the NQ-TGF. The above phenomena reveal that the normal fault system significantly controls moderatestrong seismic activity (Li, 2005). Nonetheless, the kinematics and stress characteristics of normal fault segments with concentrated moderate-strong earthquakes require further investigation. Analyzing normal fault activity and segmentation in southern Tibet is necessary for elucidating seismic activity patterns in the area, conducting seismic hazard assessments, promoting earthquake prediction, and researching seismogenic mechanisms.
Most studies on the normal fault kinematics in southern Tibet are based on geological investigations (Wu et al. 2005;Kali et al., 2010;Song et al. 2013;Yang, 2014;Zhao et al. 2019;Ha, 2019), primarily focusing on slip characteristic analysis of partial NQ-TGF segments. Global positioning system (GPS) data in recent years have contributed greatly to the rapid development of crustal movement research on the southern plateau. Using high-precision GPS monitoring data, Wang et al. (2002) and Zhang et al. (2004) estimated the current extension rate of the plateau, and Qu et al. (2021) revealed the crustal tectonic movement and deformation patterns of the plateau over the past 20 years. Results indicate that the horizontal deformation in southern Tibet is much larger than the vertical deformation, with a maximum of approximately 35 mm/yr. Thus, horizontal deformation may be the main factor controlling the earthquakes. However, owing to the sparsity of GPS sites, comprehensively and quantitatively determining the fault activity using GPS data is difficult. Many scholars have used numerical simulation method based on GPS velocity data to study dynamic problems in southern Tibet (Liu and Yang 2003;Ye and Wang 2004;Wang et al. 2006;Xu and Zhao, 2009;Pang et al. 2018), such as the extensional state uplift and formation mechanisms. These studies, which are largely based on three-dimensional rheological models, emphasize the influence of the rheological layering structure of the lower crust on the surface crustal deformation, but rarely consider the influence of faults. However, crustal deformation is mainly concentrated in the fault zone (Ismail-Zadeh et al., 2007). Faults significantly affect the deformation field distribution of the Qinghai-Tibet Plateau (Zheng et al., 2007), and fault activity determines the seismicity distribution of the primary seismic belts (Song, 2010). Additionally, the strong earthquakes (M ≥ 7) usually occur on active faults of block boundaries with high sliding velocities (Wang, 2005). The two-dimensional finite element numerical simulation method has been used to obtain reliable results regarding the correlation between the deformation movement of fault zones and earthquakes in the Qinghai-Tibet Plateau, Sichuan-Yunnan, and North China regions. However, research on the relationship between complex faults (especially the NS central normal faults) in southern Tibet and moderate-strong earthquakes (M ≥ 5) is limited; therefore, further studies are needed.
This study used long time-scale GPS velocity field data (Zheng et al., 2017) and the finite element numerical simulation method to evaluate the present-day tectonic deformation of the southern Qinghai-Tibet Plateau. We present the simulation results of the crustal stress field distribution, and the movement and stress of each fault in the study area. Finally, the correlation between NS-trending fault activity and the concentration of moderate-strong earthquakes (M ≥ 5) since 1976 is discussed in combination with focal mechanism solutions and geological investigations.

Finite element model construction and parameter selection
Given that fault kinematics rather than dynamics were mainly studied, and factors such as gravity and mantle drag were not considered, we developed a two-dimensional viscoelastic finite element model, using Ansys finite element analysis software. This model was constructed based on the active blocks, major active faults, and moderate-strong earthquakes in the region shown in Fig. 1. Most differences in the physical properties of the Lhasa, Qiangtang, and Himalayan blocks in the model were observed for the Young's modulus parameter. The model covered 12 faults distributed on the edge and inside of the Lhasa block, the shapes of which were described in detail according to the active tectonic map of China (Deng, 2007). The southern boundary faults of the Lhasa block were the Karakorum (KKF), Himalayan (HMLYF), and Mangion (MGF) faults; the Himalayan fault zone included the main mountain front rupture and the North Himalayan fault (NHMLYF). The northern boundary faults were the BangongCo (BGCF), GerzedongCo (GZ-DCF), GyaringCo (GRCF), BengCo (BCF), and Jiali (JLF) faults. Inside the model, the PLF, DRYCF, JGF, and NQ-TGF in the NS-trend were arranged in parallel. These faults were regarded as weak zones, with widths of 6 km. Table 1 listed block and fault physical parameters referring to related studies (Liu and Yang, 2003;Wang, 2005;Xu and Zhao, 2009;Jiang et al., 2009), wherein Young's modulus of the fault was 1/10 th of the blocks (Zhu et al., 2016). The model was meshed using triangular surface element of plane182, with a total of 6619 elements, as shown in Fig. 2. The size of the element was controlled at a maximum of 10 km, and the element of the fault is about 6 km.

Constitutive relationship for the viscoelastic model
To consider the crustal brittleness and rheological characteristics in the study area, the deformation of the lithosphere in geological time can be regarded as a Maxwell viscoelastic model composed of a linear elastic part and an ideal viscoelastic part in series. The elastic and viscous parts obey Hooke's law and Newton's flow law, respectively. With stress σ, the constitutive relationship is expressed as: The stress rate ( ̇) and strain rate ( ̇ ) are approximated by Eq. (2). The elastic coefficient is Young's Δt modulus E; the stress of the viscous part has a linear relationship with the strain rate, and the proportionality coefficient is the viscosity coefficient η. At time t, the strain increment includes the elastic d e and viscous {d v } strain increments in the form of: In the formulas: {} represents the component of the strain tensor; e and v represent the elastic and viscous strains, respectively; σ t is the stress at time t; dσ is the stress increment; [ E e ] is the elastic material matrix related to the elastic modulus; and η v is the viscous material matrix related to the viscous coefficient.
The Maxwell viscoelastic model has the advantage of being able to handle the elasticity-stickiness transition more naturally. It behaves as an elastic body on a short time scale, and when subjected to a long-term external stress, the potential energy stored in the elastic body will gradually disappear in the viscous body, and then stress relaxation occurs. Viscoelastic relaxation time is generally defined as the ratio of viscous hysteresis coefficient and Young's modulus. When the rheological time is much larger than the relaxation time, the stress is relaxed and stabilized (Cao et al., 2009). The viscoelastic relaxation time of the model in this study is approximately 3,500 years, while the rheological properties of the brittle crust can take over 10,000 years to manifest; thus, the simulation calculation must be loaded for an adequate time.

Model loading and calculation
The model adopted boundary loading. According to the 1991-2015 GPS velocity data of mainland China (Zheng et al., 2017), the loading displacement of boundary nodes was obtained by performing doublespline interpolation on the velocity of GPS stations around the model boundaries (Fig. 2). In this process, the data of stations near the fault or inconsistent with the surrounding movement were excluded. In cases where there were multiple GPS stations around the boundary nodes, the average value was obtained before interpolation. The load step time was set to 1.5 Ma, and hydrostatic analysis was conducted to simulate crustal motion and stress field evolution under the driving environment provided by GPS observations. When loading was approximately 0.8 Ma, the displacement and stress fields tended to be stable. After evaluating the reasonableness of the calculation results, the simulation results with a loading time of approximately 0.9 Ma were selected for further analysis.

Reasonable analysis of model results
The results of simulation and observation at GPS stations within the model showed the same movement trends and generally small deviations (Fig. 3a).
The residual values at the GPS stations, representing the difference between the observed value and the modeled value of the GPS station, were small and disordered (Fig. 3b), indicating that the residual distributions did not contain comprehensive motion information. A total of 440 east and north velocity component residuals were counted, in which the residuals < 1, 2, and 2-3 mm/yr corresponded to 52.5, 82.7, and 89.5% of the total, respectively (Fig. 3c). According to the statistics of the 220 directional angle errors between the simulation and actual measurements (Fig. 3d), errors of < 5, 10, and 15° corresponded to 80.1, 82, and 98.4%. The above results show that the results calculated by the 2D finite element model are reasonable.
The simulation results of the fault properties and slip rate values were consistent with the geological results (Table 2), indicating that the simulation results generally reflect the kinematics of regions and faults, meaning the finite element model is suitable. Notably, the simulation rate of some faults was lower than the geological result, which may be because the model was a continuous-medium model and the faults were simplified to easily deformable strips, which inevitably weakened the role of faults in the tectonic evolution of the Tibetan Plateau (Song, 2010). Moreover, the width of some faults in the model was narrower than the actual faults, and the absorbed deformation was relatively small. In addition, the geological data was only local results of the fault section, did not represent the movement rate of the entire fault.
3.2 Regional deformation field, stress field, and fault influence The velocity field and its speed contour provided by the simulation (Fig. 3a) showed a fast regional velocity field of 13-36 mm/yr with a highly uneven spatial distribution. The velocity gradually decreased from south to north and diverged from the center to the east and west. Specifically, the velocity field of the Himalayan tectonic belt was approximately 36 mm/yr, and the northern border of the research area exhibited a motion of 13-25 mm/yr, which was shortened by 11-23 mm/yr in the NS direction. Bounded by 90° E (near the NQ-TGF), the eastern region velocity was 36 mm/yr in the middle and 18.6 mm/yr in the eastern limit, with a maximum decrease close to 17 mm/ yr. The direction rotated clockwise from NE to near EW, with the western region gradually decreasing to 36 mm/yr in the middle and 15.1 mm/yr in the northwest, showing a maximum decrease of 21 mm/yr and shifting the direction from NE to near NS. These results reflect the NS-trend of compressional deformation and the EW-trend of extensional deformation in the research area. The EW-trending extension was particularly strong in the area east of the NQ-TGF, and the plateau material was extruded rapidly eastward, consistent with the numerical simulation results of Wang et al. (2006). The calculated stress field was obtained and presented inhomogeneity in the following aspects.   Figure 4 showed the direction of principal compressive stress changed from NS in the west to NE in the east. The magnitude of the principal extensional stress in the area was generally lower than the principal compressive stress, indicating that the compressional deformation was predominant and accompanied by tensile deformation. The ratio of principal extensional stress/principal compressive stress varied from region to region. The ratio was about 1/3 in the south and about 1/1 in the north. These results were consistent with the principal strain field results obtained using GPS Wu et al., 2020;Qu et al., 2021). The effect of the faults is evident in the velocity contour distribution (Fig. 3a). The velocity contours were bent on both sides of the EW-and NW-trending faults, indicating that the NEE movement was slower on the north side of the fault than on the south side. On both sides of the internal NS-trending normal faults, particularly the DRYCF and NQ-TGF, the velocity contour curve showed that movement was slower on the east side than on the west side. It is observable that the fault can absorb deformation and impede northward movement and east-west extensional deformation in the study area, changing the movement pattern. Additionally, the stress was significantly higher on the south side of the fault at the southern boundary of the Lhasa block than on the north side (Fig. 4), indicating that the fault hinders stress transfer and promotes local high-stress accumulation.

Movement and stress characteristics of the NWand EW-trending fault zones on the southern and northern boundaries of the Lhasa block
Overall, the simulation results of the tensile-compressive and strike-slip rates, which were the slip rate components perpendicular to and along the fault strike, and corresponding normal and shear stress of faults (Fig. 5) indicated variability at multiple faults along the southern boundary. The KKF was subjected to remarkable right-lateral shear and compression, dominated by right-lateral strike-slip with thrust movement. The strike-slip rate was 5.5-9.0 mm/ yr, and the compressive rate was 4.3 mm/yr. As the western border of the Qinghai-Tibet Plateau, the KKF played an essential role in the eastward plateau movement. Furthermore, the middle segment of the HMLYF absorbed a large amount of convergent deformation, and the compression rate gradually decreased from 10.1 mm/yr in the middle to 2 mm/yr in the west and east ends. The strike-slip rate decreased from 9.9 mm/yr in the middle to 5 and 3 mm/yr in the west and east ends, respectively. The most significant compression rates were observed in the central and western sections. Contrastingly, the largest strike-slip rate was in the eastern segment.
These findings indicate that the northern Himalayan orogenic belt is now dominated by extensional structures, consistent with previous geological findings (Zhang and Ding, 2003). Regarding stress, the middle and eastern segments were subjected to relatively strong compression and strong right-lateral shear, respectively, indicating deformation coexistence with NS-trending shortening and EW-trending tension in the Himalayan orogenic belt. The tension activity was strongest in the eastern segment of the fault. The focal mechanism solution (Xu et al., 2005) showed that the central segment of the Himalayan fault is dominated by thrust-type earthquakes, whereas the eastern segment is dominated by strike-slip earthquakes. Nonetheless, relevant geological research on the Magnion fault zone was relatively scarce. The simulated strike-slip and tensile rates were 0.7-5.0 mm/yr and 5.0-7.7 mm/yr, respectively, indicating considerable left-lateral strike-slip and tensional normal fault characteristics, which is consistent with the numerical simulation obtained by Bao (2012) and others. The fault activity of the northern boundary of the Lhasa block is closely related to its strike. For instance, the NWW-and EW-trending faults showed right-lateral strike-slip and compression movements, whereas the NW-trending faults showed right-lateral strike-slip and tension movements. Meanwhile, the faults exhibited significant differences in segmentation. For example, in the BGCF and GZ-DCF, an interphase distribution dominated by compression and strike-slip was Fig. 3 GPS velocity field, simulation results and residual distribution. a Horizontal velocity contours and vectors of simulation and GPS velocity field, b residual distribution, c statistical distribution of residuals, d statistical distribution of direction angle ◂ observed along the fault, and the fault segment with significant compression had a weak shear. These results reveal that the northern fault zone of the Lhasa block is an important adjustment zone for NS-trending compression and EW-trending extension. However, the slip rate values differed significantly in different fault sections. For instance, simulated strike-slip rates of the northwest, middle, and southeast segments of the JLF were 3.0-6.0, 0.5-1.5, and 3.0-4.0 mm/yr, respectively; the strike-slip rate decreased, and then gradually increased from west to east along the fault. Previous studies have also demonstrated these movement features (Ren et al., 2000;Shen et al., 2001;. The slip rate of the JLF was mostly lower than the 10-20 mm/yr suggested by the "extrusion-escape model" (Armijo et al., 1989), and higher than the 2-3 mm/yr suggested by the "horizontal shortening thickening model" (England and McKenzie, 2010). The present tectonic model of the Qinghai-Tibet Plateau is complicated and may be between the "extrusion-escape" and "horizontal shortening thickening" models; this is consistent with the views of Li et al. (2021).  The NW-trending GRCF and BCF are large-scale right-lateral strike-slip faults and strong seismic activity belts in southern Qinghai-Tibet (Wu et al., 1994;Yang et al., 2011). These two faults are located at the intersection of the NS rifts and the EW fault, and the pull-apart basin is the main tectonic landform in this region (Wu et al., 1994;Ning et al., 2006). The simulation revealed that right-lateral strike-slip dominated the GRCF and BCF, wherein their eastern segments exhibited a 6-7-mm/yr and 3-6-mm/yr strike-slip, respectively. Similar results were obtained in the geological survey performed by Yang et al. (2011).
Finally, the dextral and normal motion characteristics of the fault were closely related to its background stress field and strike characteristics. Compared with other near-EW-trending faults on the northern boundary of the Lhasa block, NW-trending faults were more affected by NS-trending compression and NStrending rifts.

Movement and stress characteristics of the NS-trending normal faults and their relationship with earthquakes
The NS-trending fault zones in the central part of the Lhasa block were mainly normal faults, with an extension rate of 1.1-7.2 mm/yr (Fig. 5). The total extension rate of the four normal faults was 11-18 mm/yr, indicating their strong ability of absorbing EW-trending deformation. The principal stress of NS-trending normal faults (Fig. 6) showed that principal tensile stress along most segments was significantly higher than the principal compressive stress. The direction of principal tensile stress was perpendicular to the fault orientation, and the magnitude along the eastern faults was the most prominent. The modeled results of motion and stress of faults were generally consistent with the focal mechanism solutions of strong earthquakes in this region, which are primarily normal fault types. The simulation results (Fig. 5) showed that the movement and stress distribution of the NQ-TGF zone were related to fault strike. Regarding the fault properties and activity intensity, the northern segment of NQ-TGF was NE-trending, with tensional and left-lateral strike-slip rates of 2.2-5.7 mm/ yr and 2.9-5.4 mm/yr, respectively. Moreover, the middle segment was near-NS-trending, with tensional and right-lateral slip rates of 5.5-9.1 mm/yr and 1.1-4.1 mm/yr, respectively; the southern segment was NNE-trending, with tensional and rightlateral slip rates of 2.1-4.5 mm/yr and 4.9-7.0 mm/ yr, respectively. Overall, the strike-slip rate first decreased and then increased from north to south, while the extension rate first increased and then decreased, consistent with the geological results (Armijo et al., 1986) and the GPS observation . Furthermore, the NE-and NStrending sections of the northern segment exhibited tensile stress accumulation rates of 200-400 Pa/ yr and 600-900 Pa/yr, respectively. The middle and southern segments had rates of 800-1600 Pa/yr and  -500 Pa/yr, respectively. That indicate that tensile stress accumulates faster in the NS-trending sections of the middle and northern segments of the fault than in the southern sections. Notably, the direction of principal stress changed with the fault strike. The maximal principal stress was the principal tensile stress in the middle and southern segments, with the near-EW direction perpendicular to the fault strike, while it was the principal compressive stress in the NE-trending section of the north segment, with a direction angle of 30-60° relative to the fault, which was conducive to promote left-lateral slip and weaktensile motion. Overall, in the middle NS-trending sections of the NQ-TGF, the rates of extensional movement and tensile stress accumulation were relatively large, the principal tensile stress was dominant, and the tensile force and fault activity intensity were highly remarkable. Similarly, the DRYCF also showed characteristics of fault movement related to the strike. The tension rates in the NNW-, NS-, and NEE-trending fault segments were 1.1-2.3, 3.5-4.9, and 2.6-3.8 mm/yr and the tensile stress accumulation rates were 200-630, 800-1200, and 400-800 Pa/yr, respectively. Obviously, among the segments of DYCF, the extensional movement and tensile stress were the most significant in the NS-trending northern segment. Compared with other normal faults, the JGF as a whole had a relatively large tensile rate of 3.9-6.1 mm/yr and had a rapid accumulation of tensile stress, with a cumulative rate of 680-1300 Pa/yr. According to historical earthquake distribution, moderate-strong earthquakes were evenly distributed across and around the fault. In addition, the overall tensile rate and tensile stress of the PLF were not high, but the segmental characteristics were significant. Compared with the northern segment, the central and southern segments of PLF were more active, with left-lateral strike-slip and tensile rates of 3.0-4.2 mm/yr and 3.3-5.0 mm/yr, respectively.
Overall, the fault movement and stress simulation showed a noticeable difference in the mechanical properties and movement patterns between faults and fault segments. Some fault segments, such as the central-southern segment of PLF (NNE-trending), NS-and NNE-trending fault segments of DRYCF, and the NS-trending middle segment of the NQ-TGF and JGF faults, often had high tension movement and tensile stress accumulation rates. The direction of the principal tensile stress of above-mentioned segments was near-EW-trending and perpendicular to the near-NS-trending fault strike, which was conducive to the occurrence of normal fault-type earthquakes. Therefore, moderate-strong earthquakes (M ≥ 5, since 1976; or M ≥ 6.5, overall) were concentrated on above-mentioned fault segments. In addition, those NS-trending segments aggregated by earthquakes mostly intersect with near-NNW-trending faults, the special geometric trend distribution may also be beneficial to stress accumulation and earthquake breeding. The above results indicate that the fault segments with high tensile rates and near-NS direction partially controlled moderate-strong earthquakes.

Conclusion
This study used a two-dimensional viscoelastic finite element model to simulate the deformation and stress fields of the main faults in the southern Qinghai-Tibet Plateau. Combining geological data, focal mechanism solutions, and other data, the correlation between the intensity of fault activity and the concentrated distribution of earthquakes was analyzed.
The results of the velocity and stress fields in the southern Tibetan region indicated that the research area is simultaneously subjected to the combined effects of NE-trending compression and inhomogeneous EW-trending stretching, producing NS-trending compressional deformation and EW-trending extensional deformation. Furthermore, the tensile stress magnitude gradually increased from south to north.
Faults have changed the pattern of crustal deformation and movement in the interior of the Qinghai-Tibet Plateau, as they can absorb deformation, regulating the deformation and stress inhomogeneity of the region.
The simulated movement and stress accumulation rates of fault revealed evident activity variability and heterogeneous of the fault segments, which affected and even controlled the regional crustal deformation and rupture types of moderate-strong earthquakes. The concentrated distribution of moderate-strong earthquakes (M ≥ 5, since 1976; or M ≥ 6.5, overall) on the central normal fault largely depended on the slip rate and geometric trend of the fault. The moderate-strong earthquakes were generally concentrated in fault segments with high slip rates and structurally complex parts, such as the intersection of multiple faults in different directions.