According to Chinese Seismic Network, the large Nepal earthquake occurred on 25 April, 2015, the magnitude is measured as MS8.1, the epicenter is located nearby Pokhara of Nepal (84.7°N, 28.2°E), the distance is about 80km from the epicenter to Nepal's capital Kathmandu, the hypocenter depth is about 15km, the earthquake occurred in the boundary zone of India Plate and Eurasian Plate. The Maximum intensity is Ⅹ degrees in Pokhara, heavy intensity area extended eastward. The earthquake caused more than 7, 000 deaths, more than 14, 000 injuries and huge property losses. The earthquake fault belongs to thrust type, the fault strike is about 290 degrees and the dip angle of fault is small. InSAR and GPS data show that the main co-seismic deformation distribution range is 150km×100km (Shan Xinjian et al., 2015). The northern deformation zone is wider than the southern on both sides of the fault, and its maximum slip distance is more than 5m. The earthquake caused southward movement of the Qinghai-Tibetan block, We may analyze and understand the effects of the earthquake on the seismic activity in the Qinghai-Tibetan Pleatau through finite element numerical simulations.
Due to the development of finite element method to the field of geology, some valuable research results have been obtained in the Qinghai-Tibetan area (Tapponnier P. et al., 1976; Fu Rongshan et al., 2000; Lu Shikuo et al., 2004; Deng Zongce et al., 1990; Dai Liming et al., 2011). The finite element method has become one of the important methods to study and explore the dynamic process of this special area. In the Qinghai-Tibetan Plateau the tectonic stress field simulation results showed that Indian Plate in a certain speed is continuing to squeeze into the north, and this is the main driving force of the Qinghai-Tibetan Plateau's uplift, and further more the uplift process is an inhomogeneous evolution process (Zheng Hongwei et al., 2006). Zhang Dongning et al. (1995) applied the viscoelastic power-law creep constitutive relation to explore the possible mechanism of stress distribution in the tectonic plateau, and their result show that the upper crust of the southern Qinghai-Tibetan Plateau can still be in tensional stress state, which may mainly exist in the upper crust of the southern Qinghai-Tibetan Plateau. Cao Xuefeng et al. (2003) have shown the correlation between seismic activities on the North-South Seismic Belt, which may be the result of regional stress field equilibrium, and the overall performance of the function of the block's boundary force (external load). After a strong earthquake occurs, the seismogenic fault is in the stage of strong stress unloading, while the coulomb rupture stress loading effect is dominant in other potential seismogenic faults. The interaction between the two strong earthquakes is mainly due to the loading of the stress. The occurrence of a strong earthquake may accelerate the occurrence of the next strong earthquake, causing a strong earthquake or multiple earthquakes until the regional strain energy reaches to a low level state, then the region will enter a new low seismicity period.1 THE VISCOELASTIC MODEL ESTABLISHMENT OF THE FINITE ELEMENT METHOD 1.1 Entity Model and Meshing
The model area included 2 plates (India Plate and Eurasian Plate), 12 sub-blocks (in Eurasian Plate) and 15 boundary zones of blocks and plates (Zhang Peizhen et al., 2003; Zhang Guomin et al., 2005) (Fig. 1). According to Zhang's results, a first-order boundary is a boundary between active tectonic-regions whick include several active tectonic-blocks and a second boundary is a boundary between two active tectonic-blocks. Among these boundaries, 8 were first-order boundaries, and 7 were secondary boundaries. From the CRUST 1.0 model(Artemieva I. M. et al., 2001) we obtained the P- and S-wave velocities, medium density, and geometric parameters of the upper, middle, and lower crust in a 1° × 1° grid and used these parameters to calculate the Young's modulus and Poisson's ratio (Yao Xiuyun et al., 1989). For the upper crust, a fully elastic model was used, while a viscoelastic model was used for the middle and lower crust.The coefficients for the middle and lower crust were 1019 and 1020Pa ·s, respectively, and the viscosity of each layer was uniform. Based on the establishment of entity model, we use Prony series form for describing the shear relaxation kernel function G(t) and the volume relaxation kernel function K(t). We selected solid 187 elements (3-D 10-node tetrahedral solid structure element) in the meshing, the average length of element grids is 20km. The model grid is divided into 488, 707 elements with 740841 nodes.
The description of the fault zone took two forms: (1) the number of curves and lines (i.e., a fracture is connected with multiple lines), and (2) different width zones were composed of a series of quadrilateral areas(i.e., a description of the actual width of the fracture zone). The width of our block boundary is determined by the results of Zhang Guomin et al. (2005) study on the width of the boundary bandwidth. To weaken the material parameters of all elements located in the boundary zones of these fault zones or blocks (Chen Lianwang et al., 2011), fault parameters were set to: (1) Young's modulus=original Young's modulus/3, and a Poisson's ratio=original Poisson's ratio +0.02 for the upper crust elastic layer; (2) Young's modulus=original Young's modulus/3, Poisson's ratio=the original Poisson ratio +0.02, and viscosity coefficient=the original viscosity coefficient/2 for the viscoelastic middle and lower crustal layers (the original value is calculated from the CRUST1.0 model).1.3 Boundary Conditions and Seismic Fault Dislocation Loading of the Model
GPS data were selected for the model area and these data were interpolated in a 2-D plane, and the annual displacement of data points was obtained. The displacement of the side boundary was restrained, and the vertical displacement condition was equivalent to the side boundary. We assume that the surface of the model was free, and the bottom of the crust was not displaced perpendicularly to the Earth's surface. A horizontal displacement constraint was applied to the model structure and the horizontal displacement component was zero. According to field investigation and coseismic data of the Nepal MS8.1 earthquake, the total length of the rupture is about 120km, and the maximum surface displacement of USGS fault slip model is 3.5m (Zhang Bei et al., 2015). Since the earthquake occurred within the boundary zone of the block, the width of the boundary zone was not less than 60km. The loading width of seismic dislocation fault is 60km (30km on both sides of the fault) and the depth of the fault was 25km (assuming the focal depth is 15km). For the thrust of positive segment dislocation, only the hanging wall was loaded. The seismic source depth was uniformly loaded to the top of the upper crust, and other directions were linearly attenuated. The minimum loading dislocation is 0.1m.2 CALCULATION AND RESULTS
The deformation process of viscoelastic material after loading is a time-varying process, and the recovery process after unloading is a delayed process. Therefore, the stress in viscoelastic material is not only related to the strain at that time, but also to the whole history of strain change.Therefore, a reference time needs to be selected for the calculation of the viscoelastic model. The calculation of each earthquake is divided into two steps. The first step is to calculate the displacement and strain before the earthquake occurred. The second step is to calculate the corresponding results after the earthquake is loaded. The displacement field and the stress-strain field of the earthquake can be obtained by calculating the difference between the results of seismic dislocation loading and the results before the earthquake. The following is an analysis of the Nepal MS8.1 earthquake.2.1 The Distribution Characteristics of Displacement and Strain Induced by the Earthquake
Coseismic displacement: The Nepal MS8.1 earthquake, which occurred on April 25, 2015, located in the middle of the Himalayan zone. The seismogenic fault is thrust fault. The maximum dislocation value exceeded 3.5m, the actual calculation is close to 3.6m. The distribution area with displacement larger than 0.3m is mainly located in the central area of the Lhasa block in the north of the seismic rupture zone. When it is far away from the epicenter, there is an obvious attenuation of the displacement as long as it through a boundary zone. In the vicinity of the southern boundary of the meizoseismal area, the displacement direction is NNE, which is opposite to the north. The direction of the north of the meizoseismal area is southward or points to the epicenter (Fig. 1). The average displacement distribution induced by this earthquake is as follows: The Lhasa block is 0.5m, the Qiangtang block is 15cm, the Bayan Har block is 4cm, the south of the Tarim block and the western part of the Qaidam block are about 1cm, the displacement of the northeastern, eastern and southeastern edges of Qinghai-Tibetan Plateau is about 1mm. This earthquake mainly caused the southward movement of the Qinghai-Tibetan Plateau. The western part of the Bayan Har block and the west of the Qaidam block form a region with a small amount of coseismic displacement.
Characteristics of coseismic strain: There are two main positive and negative strain areas caused by the earthquake. Himalayan zone and some region of the southern part of the Lhasa block are the increase region of strain. The equivalent strain near the epicenter increased by 0.06×10-5-0.09×10-5 (Fig. 2(a)). Maximum shear strain increased by 0.01×10-5-0.065×10-5 (Fig. 2(b)). There is an attenuating discontinuous area in the east of the epicenter. There are three strong aftershocks with magnitude exceeding 7 in this area. And the maximum aftershock is MS7.5, which occurred on May 12, 2015. The unloading region of strain are mainly distributed in the north of the epicenter with a width of about 10 degrees along the direction of NNE. And it is located in the midwest of the Qinghai-Tibetan Plateau. The reduction range of equivalent strain in this region is from 0.03×10-5-0.12×10-5 (Fig. 2(a)). The maximum shear strain decreases in the range of 0.01×10-5-0.07×10-5 (Fig. 2(b)). Similar to the distribution characteristics of the displacement field, the midwest of Bayan Har block and the western region of the Qaidam block constituted an area of positive coseismic strain. And this plays a loading role. The negative strain region is the unloading area of strain energy. Except for the fracture zone, the positive strain area is the loading area of strain energy.
In order to study the influence of the Himalayan earthquake on the Qinghai-Tibetan Plateau, we choose the relatively complete period of the earthquake (MS≥5.0) catalogue (Chinese Seismic Network has been improved since 1970). So far, there are 6 earthquakes with MS≥7.0 in Himalayan zone (Table 1). Three of them are strong aftershocks of the Nepal MS8.1 earthquake on April 25, 2015. And aftershocks mainly occurred in the east of the main shock. The epicenter of the MS7.5 aftershock is about 200km from the epicenter of the main shock, which occurred in the stress and strain loading area caused by the main shock. The seismogenic faults of these earthquakes are thrust faults or mainly composed by thrust faults. A MS7.1 earthquake occurred on August 21, 1988 is locate at about 140km southeast of the MS7.5 aftershocks. The earthquake is close to the southern boundary of the Himalayan zone, and most of the seismic faults are thrust faults. During the period of about 3 years after the earthquake (August 21, 1988 to July 2, 1991), the western and southern Tibetan block are very quiet and basically no earthquakes of MS≥5.0 occurred (Fig. 4). Which may have resulted in the release of partial stress in the northern part of the MS7.1 earthquake. Its seismic quiet areas is located in the stress-unloading area. On October 2, 1991, the MS7.0 earthquake occurred in the northwest of the MS8.1 earthquake in 2015. The distance between the two earthquakes is about 640km. Both belong to the Himalayan zone. The seismogenic fault is a thrust fault. Within 3 years after the earthquake (actually 7 years), the earthquake activity is calm in the southwestern part of the Qinghai-Tibetan block. There is no MS≥5.0 earthquakes(Fig. 5). The seismic quiet area is located in the stress unloading area.
The displacement field and strain field generated by the 2015 Nepal MS8.1 earthquake were used to analyze the changes of seismic activity around the Qinghai-Tibetan Plateau after the earthquake. The earthquake occurred in the suture zone of the Qinghai-Tibetan block and the Indian plate. The strike of the seismogenic fault is NWW and the fault is dominated by thrust and has a right-handed component. The stress and strain generated by the earthquake are weakly four-quadrant distribution and the stress and strain at the two ends of the fault are increasing, which acts as a stress loading. The other two areas are the upper and lower plates of the fault and the stress and strain are decreasing, which acts as a stress unloading. In the southern part of the seismic fault, the displacement field direction is NNE, the northern part of the fault is basically pointing to the epicenter, the displacement direction of the west end of the fault is WS, and the displacement direction of the east end is ES (Fig. 1). The northern part of the earthquake fault is a stress-strain unloading area, mainly involving the western part of the Qinghai-Tibetan Plateau. Due to the unloading effect of the 2015 Nepal MS8.1 earthquake on the stress and strain in the western part of the Qinghai-Tibetan Plateau, since then the seismic activity is very quiet in the western part of the Qinghai-Tibetan block (Fig. 3 and Fig. 6), and this state may continue for much longer time.
The three earthquakes are all located in the Himalayan zone, border of the India plate and the Eurasian plate, and the seismogenic faults are all belong to the thrust type or mainly composed by thrust faults. Each of them can produce a stress relief process for the Qinghai-Tibetan block, especially to the western regions after they occurred. And it had unloaded some part of stresses. Then the activity of strong earthquakes can remain stable for a period of time. Although the magnitudes of the three earthquakes used are difference, but its influence on Qinghai-Tibetan block's seismic activity are basically the same, the large earthquakes in the area where had obviously unloaded will have a quiet period in a time scale of years.3 CONCLUSIONS AND DISCUSSION
The following conclusions can be obtained by the finite element simulating calculation and seismic activity analysis:
(1) After each thrust type earthquake in the Himalayan zone, there is an impact on a certain range of seismic activity in the north of the epicenter of Qinghai-Tibetan block, and it mainly plays the role of unloading stress. After that, seismic activity has a relatively stable stage.
(2) After the occurrence of the large earthquake, its surrounding stress and strain field will be changed, and the obvious loading and unloading areas will be formed at the same time. The unloading area of the stress-strain is the quiet area of the earthquake activity and there are differences in the calm time scale.
(3) Compared to the average three-dimensional model used in the finite element simulation with the predecessors (generally the parcel is treated as a uniform entity, that is, the same parcel is described by a set of mechanical parameters), but actuallythe non-uniformity between the blocks is obvious. In this paper, the global crustal CRUST 1.0 model is used to divide the crust into upper, middle and lower layers. The model is more refined and the simulation results are closer to the actual situation.
(4) This paper only calculates and analyzes the relationship between deformation field characteristics and seismic activity in a time phase after an earthquake. However, the change of the deformation field is a slow process of time accumulation. The analysis result of a single object has certain limitations. In the future research, it is expected to conduct a comprehensive analysis of the deformation field caused by multiple earthquakes in a certain period of time.
After any major earthquake, the stress in the source region is greatly released. If the stress released completely, a major earthquake occurs again in the source area requires a long time's stress accumulation. There is still some change in the stress field within a certain range beyond the earthquake rupture zone. The most important feature is the loading and uploading of stress. The unloading area is a region where the earthquake is calm, while the opposite loading area is a region where the earthquake is more active, and its active degree is closely related to the tectonics. In a word, it is necessary to know the distribution of theloading and unloading areas generated by the large earthquake. Thus we can make a judgement about the future trends of seismic activity and the activity areas.ACKNOWLEDGEMENTS
The authors thank GAN Weijun (Institute of Geology, China Earthquake Administration), who provided valuable GPS data. The authors also thank anonymous reviewers for providing feedback that has greatly enhanced this manuscript.
Artemieva I.M., Mooney W.D.Artemieva I.M., Mooney W.D. Thermal thickness and evolution of Precambrian lithosphere: a global study[J]. Journal of Geophysical Research, 2001, 106(B8): 16387-16414. DOI:10.1029/2000JB900439
Cao Xuefeng, Yin Jingyuan, Yang LimingCao Xuefeng, Yin Jingyuan, Yang Liming. Study of simulation for the evolution of stress field on north-south seismic belt in China[J]. Northwestern Seismological Journal, 2003, 25(3): 221-225, 245 (in Chinese with English abstract).
Chen Lianwang, Zhang Peizhen, Lu Yuanzhong, Chen Huaran, Ma Hongsheng, Li Li, Li HongChen Lianwang, Zhang Peizhen, Lu Yuanzhong, Chen Huaran, Ma Hongsheng, Li Li, Li Hong. Numerical simulation of loading/unloading effect on Coulomb failure stress among strong earthquakes in Sichuan-Yunnan area[J]. Chinese Journal of Geophysics, 2008, 51(5): 1411-1421 (in Chinese with English abstract).
Chen Lianwang, Zhan Zimin, Ye Jiyang, Li YanChen Lianwang, Zhan Zimin, Ye Jiyang, Li Yan. Numerically modeling the influence of rheological properties on tectonic deformation of Tibet plateau[J]. Journal of Geodesy and Geodynamics, 2011, 31(3): 8-14 (in Chinese with English abstract).
Dai Liming, Li Sanzhong, Tao Chunhui, Li Xishuang, Liu Xin, Suo Yanhui, Zhou Shuhui, Zhou Yonggong, Gao WubinDai Liming, Li Sanzhong, Tao Chunhui, Li Xishuang, Liu Xin, Suo Yanhui, Zhou Shuhui, Zhou Yonggong, Gao Wubin. Three-dimensional numberical modeling of activity of the Longmenshan fault zone driven by the India Plate[J]. Progress in Geophysics, 2011, 26(1): 41-51 (in Chinese with English abstract).
Deng Zongce, Cui Junwen, Feng Xiaofeng, Wu Changde, Zhu HongDeng Zongce, Cui Junwen, Feng Xiaofeng, Wu Changde, Zhu Hong. The analysis of tectonic dynamics and deformation mechanics of Qinghai-Xizang plateau by elastic finite element method[J]. Bulletin of the Chinese Academy of Geological Sciences, 1990, 11(21): 65-75 (in Chinese with English abstract).
Fu Rongshan, Xu Yaomin, Huang Jianhua, Li LigangFu Rongshan, Xu Yaomin, Huang Jianhua, Li Ligang. Numerical simulation of the compression uplift of the Qinghai-Xizang plateau[J]. Chinese Journal of Geophysics, 2000, 43(3): 346-355 (in Chinese with English abstract).
Lu Shikuo, Cai YongenLu Shikuo, Cai Yongen. Review of numerical simulation on the dynamics of Qinghai-Xizang plateau[J]. Acta Seismologica Sinica, 2004, 26(5): 547-559 (in Chinese with English abstract).
Shan Xinjian, Zhang Guohong, Wang Chisheng, Li Yanchuan, Qu Chunyan, Song Xiaogang, Yu Lu, Liu YunhuaShan Xinjian, Zhang Guohong, Wang Chisheng, Li Yanchuan, Qu Chunyan, Song Xiaogang, Yu Lu, Liu Yunhua. Joint inversion for the spatial fault slip distribution of the 2015 Nepal MW7.9 earthquake based on InSAR and GPS observations[J]. Chinese Journal of Geophysics, 2015, 58(11): 4266-4276 (in Chinese with English abstract).
Tapponnier P., Molnar P.Tapponnier P., Molnar P. Slip-line field theory and large-scale continental tectonics[J]. Nature, 1976, 264(5584): 319-324. DOI:10.1038/264319a0
Yao Xiuyun, Zhao Debin, Zhao HongruYao Xiuyun, Zhao Debin, Zhao Hongru. Relations between longitudinal and shear wave velocities and other physical parameters[J]. Progress in Geophysics, 1989, 4(2): 7-20 (in Chinese).
Zhang Bei, Cheng Huihong, Shi YaolinZhang Bei, Cheng Huihong, Shi Yaolin. Calculation of the co-seismic effect of MS8.1 earthquake, April 25, Nepal[J]. Chinese Journal of Geophysics, 2015, 58(5): 1794-1803.
Zhang Dongning, Xu ZhonghuaiZhang Dongning, Xu Zhonghuai. An alternative interpret ation of earthquake activity caused by normal fault of upper crust in southern Qinghai-Xizang (Tibetan) plateau[J]. Acta Seismologica Sinica, 1995, 17(2): 188-195 (in Chinese).
Zhang Guomin, Ma Hongsheng, Wang Hui, Wang XinlingZhang Guomin, Ma Hongsheng, Wang Hui, Wang Xinling. Boundaries between active-tectonic blocks and strong earthquakes in the China mainland[J]. Chinese Journal of Geophysics, 2005, 48(3): 602-610 (in Chinese with English abstract). DOI:10.1002/cjg2.693
Zhang Peizhen, Deng Qidong, Zhang Guomin, Ma Jin, Gan Weijun, Min Wei, Mao Fengying, Wang QiZhang Peizhen, Deng Qidong, Zhang Guomin, Ma Jin, Gan Weijun, Min Wei, Mao Fengying, Wang Qi. Active tectonic blocks and strong earthquakes in the continent of China[J]. Science in China (Series D: Earth Sciences), 2003, 46(S2): 13-24.
Zhang Yuansheng, Zheng Xiaojing, Wang LanminZhang Yuansheng, Zheng Xiaojing, Wang Lanmin. The distribution characteristics of deformation field caused by three great earthquakes in the Qinghai-Tibetan Plateau and its vicinity since 2001[J]. Chinese Journal of Geophysics, 2016, 59(10): 3637-3645 (in Chinese with English abstract).
Zheng Hongwei, Li Tingdong, Gao Rui, He RizhengZheng Hongwei, Li Tingdong, Gao Rui, He Rizheng. The advance of numerical simulation in geodynamics[J]. Progress in Geophysics, 2006, 21(2): 360-369 (in Chinese with English abstract).