- Email: [email protected]
Migration trend of strong earthquakes in North China from numerical simulations
Journal of Asian Earth Sciences 50 (2012) 116–127
Contents lists available at SciVerse ScienceDirect
Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes
Migration trend of strong earthquakes in North China from numerical simulations Yuanzhong Lu a, Shuxin Yang a, Lianwang Chen a, Jianshe Lei a,⇑, Ping He b a b
Key Laboratory of Crustal Dynamics, Institute of Crustal Dynamics, CEA, Beijing 100085, China Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China
a r t i c l e
i n f o
Article history: Received 12 August 2011 Received in revised form 2 January 2012 Accepted 14 January 2012 Available online 17 February 2012 Keywords: Stiffness reduction Numerical simulation Distant migration Strong earthquake North China
a b s t r a c t Lu et al. (2011) established a finite element model around the entire Chinese continent to simulate the distant migration of strong earthquakes by killing the elements. In this study a new finite element model of North China has been established in the existing Chinese continental model using the GPS observations and the finite element sub-model technique through taking into account geological structure, active faults, surface topography, Moho discontinuity and 3-D velocity structure of the crust and upper mantle in the region. Furthermore, the stress field adjustment was realized by reducing the element stiffness rather than by killing the elements. We calculated the displacement of the Chinese continent model, and took the displacement located at the boundary of North China as the boundary loadings, then computed the initial stress field by considering the gravity effect. The migration of the strong earthquake regions was simulated by performing a number of tests with various kinds of short-term boundary loadings. Locations of strong earthquakes were determined with the highest G (the risk factor of the fault) value and adjustment of the stress field. Our results showed that high stress disturbance regions could be possible regions of future strong earthquakes, and that the initial stress field in North China can be well estimated, if an optimal boundary loading configuration is chosen in the Chinese continental model. These results suggest that the strong earthquake forecast in North China may strongly depend on the stress field and strong earthquake activities in a much wider area, such as the Chinese continent. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Through analyzing earthquake activities and precursory phenomena, including crustal stress and strain, geomagnetical, geoelectrical anomalies, and ground fluids, a few moderate-strong earthquakes were once successfully predicted and natural hazards were mitigated significantly (Ding et al., 2000). However, the failure of the prediction to some destructive earthquakes, such as the 1976 Tangshan (M 7.8) and 2008 Wenchuan (M 8.0), China, earthquakes, renders us to ponder the feasibility of earthquake prediction and related scientific investigations. A common consensus is that due to the complexity and nonlinearity of seismogenic structures, it is extremely hard to obtain useful information on the earthquake pregnancy and occurrence, and nowadays more difficult to predict the exact location, magnitude, and origin time. It is almost impossible to achieve a significant breakthrough in the earthquake prediction if we only depend on precursor observations. One wants to know whether there exists a possibility to predict the trend of strong earthquakes based on the deficient understanding of crustal structure, the property of fault zones, and crustal movements. Wang et al. (1982) once simulated strong earthquakes (M P 7.0) over the past 700 years, and found that most ⇑ Corresponding author. Tel./fax: +86 10 6284 6760. E-mail addresses: [email protected], [email protected] (J. Lei). 1367-9120/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.jseaes.2012.01.006
earthquakes occurred around the region with an increased value of the risk factor G = s=ðrn tan /Þ, where s is the shear stress along the fault plane, rn is the normal stress on the fault plane, and / is the friction angle. Zhang et al. (1990) jointly inverted the 1966 Xingtai (M 7.2), 1969 Bohai (M 7.4), 1975 Haicheng (M 7.3), 1976 Tangshan (M 7.8), China, earthquakes, and found that the occurrence of strong earthquakes could be resulted from the activity and extension of the fault system. Tao et al. (2000) established a viscoelastic finite element model in the Chinese continent, and obtained the distribution of strain and stress field. Chen et al. (2001) investigated a dynamic evolution process of the 1966 Xingtai earthquake, and concluded that the occurrence of the 1966 Xingtai earthquake could play a key role in the stress variation around the 1978 Tangshan source area. Bai et al. (2003) simulated the dynamic process of the Tangshan earthquake, and found that the earthquake significantly affected its adjacent blocks and border faults and led to various sizes of ground movement and surface deformation around it. Based on these studies, Lu (2004) proposed a seismic dynamic prediction model attributable to a remarkable process in the computer technique, geodynamics, deep seismic structure, GPS observations, and continental deep drillings. In the seismic dynamic prediction model, we do not only take into account one or several faults and the interaction between the block and its surrounding fault system, but also dynamic processes of stress concentration and fault slipping in the brittle and
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
brittle-ductile mid-upper crust away from the plate boundaries. Assuming a sudden weakening of a localized weak lower crust embedded in an elastic crust, Kenner and Segall (2000) concluded that weakening in the lower crust can causes stress amplification in the upper crust and result in a sequence of ruptures of the fault zone. Li et al. (2009) investigated different cases of fault weakening after a large fault rupture, and concluded that fault weakening can lead to repeated earthquakes. Lu et al. (2011) made an attempt to establish a finite element model of the entire Chinese continent (Fig. 1a) by killing the elements (no loading capability) after the earthquake, but this method cannot predict where the next earthquake occurs. However, the purpose of the dynamic prediction at the current stage is to forecast the migration trend of an incoming strong earthquake region so as to provide an important scientific evidence for capturing the seismic precursors to mitigate the seismic hazard.
117
Recently, in North China many investigations have being conducted using various approaches, such as seismic tomography (Lei et al., 2008, 2011; Huang et al., 2009), active faults (Deng et al., 2002), and GPS observations (Jiang et al., 2000), and obtained many significant results. These results allow us to model the migration trend of the strong earthquake region in North China by introducing these previous results into our model. To achieve practical and meaningful forecast of strong earthquakes with the purpose of the mitigation of seismic hazard, the region of the predicted earthquakes should not be too large. Therefore, in the present study we took North China as an example to explore the migration trend of follow-up strong earthquakes. Because the loading configuration could not be changed significantly within an active period of about 10 years, the migration of strong earthquakes could be mainly controlled by the source areas without load-bearing capacity and induced stress adjustment. Our simulation steps
Fig. 1. (a) Location of the North China sub-model (open polygon) in the entire Chinese continent finite element model of Lu et al. (2011). Black lines denote major active faults in China, while other lines denote the model parameterization. The digital numbers denote the model boundary. (b) A finite element model of North China. (c) Distribution of the faults used in the finite element model of North China. The numbers denote the faults as shown Table 1.
118
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Table 1 Some parameters of the faults (Deng et al., 2002) used in the North China sub-model. Sequence number
Name
Type Dip-slip Unknown Dip-slip Strike-slip
7 8 9 10 11
Helan Mountain Eastern Foot Fault Luo Mountain Eastern Foot Fault Qinling Mountain Northern Foot Fault Huashan Mountain Piedmont Fault Wentang Fault Zhongtiaoshan Mountain Piedmont Fault Luoyun Mountain Piedmont Fault Huo Mountain Piedmont Fault Shangu Fault Jiaocheng Fault Jizhou Mountain Piedmont Fault
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
Wutai Mountain Western Foot Fault Heng Mountain Northern Foot Fault Kouquan Fault Tangxi Fault Tangdong Fault Handan Fault Yuanshi Fault Zhangjiakou Fault Xinbaoan-Shacheng Fault Shizhuang Fault Sunhe-Nankou Fault Jiyunhe Fault Anyang Fault Heze Fault Cangni Fault Tangshan Fault Cangdong Fault Xiadian Fault Sunzhuangzi-Wulonggou Fault Liaocheng-lankao Fault Cixian-daming Fault Dongming-chengwu Fault Xinhe Fault Huanghe-Lingwu Fault Qinmu Fault Ordos Fault Pikou Fault Jinzhou Fault Gaixian Fault Liaodong Fault
1 2 3 4 5 6
Dip-slip Unkown Dip-slip Strike-slip Dip-slip Strike-slip Dip-slip Dip-slip Unknown Right-lateral strikeslip Dip-slip Dip-slip Strike-slip Unknown Unknown Unknown Unknown Unknown Unknown Unknown Strike-slip Unknown Unknown Unknown Unknown Unknown Strike-slip Dip-slip Unknown Unknown Unknown Unknown Unknown Dip-slip Strike-slip Dip-slip Strike-slip Dip-slip Unknown Unknown Unknown Unknown
are as follows. Firstly, a finite element model of North China (Fig. 1) was established based on the entire Chinese continental model (Fig. 1a, Lu et al., 2011). In the Chinese continental model they did not take into account the undulated Moho discontinuity and faults, while in present study we added both the undulated Moho discontinuity and 41 faults in the North China sub-model to obtain more accurate results. In the sub-model the stiffness diminution by changing the element’s material properties including Young’s modulus and Poisson’s ratio within the faults was designed after the occurrence of strong earthquakes so as to calculate a new stress field. Moreover, the same amount of stiffness is reduced for all the earthquake faults and the stress adjustment is time independent. Then the possible criterion was explored to determine the occurrence of subsequent simulated earthquakes. Finally, the migration trend of strong earthquakes was estimated after a multi-step cycle calculation. 2. Establishment of the finite element model in North China Lu et al. (2011) set up a finite element model in the entire Chinese continent. In order to better investigate the mitigation of strong earthquakes in North China, here we added 41 active faults with an assumed 2 km width in the North China sub-model. The
location and trending of these faults were extracted from Deng et al. (2002) (Table 1 and Fig. 1c). The boundary of the North China model (Fig. 1a) was directly adopted from Han et al. (2003). In the model the 98,433 grid nodes and 86,840 elements were set up. The elements on the faults are dense, where the grid spacing is about 2 km, while those away from the faults are relatively sparse, where the grid spacing is about 60 km. The material on the faults is very weaker than that away from the fault. The model was divided into eight layers. The first layer is from 5 km depth to the ground surface (the topography from www.usgs.gov) (Fig. 2a), while the bottom layer is between 45 km and the Moho discontinuity (Chen et al., 2002) (Fig. 2b). Other layers are evenly distributed between the first and bottom layers. The material is assumed to be viscoelastic at the bottom layer of the mantle, while those at other layers are assumed to be linear elastic. The commercial software ANSYS 12.1 was used to establish the model with SOLID185 hexahedral elements and finish some related calculations. The ANSYS software is suitable for numerical simulations in a variety of physical problems, such as solid mechanics, thermal, electromagnetic and fluid dynamics, but recently it is popularly used in Earth sciences (e.g., Parsons, 2002; Seyferth and Henk, 2002; Zhu et al., 2004; Lu et al., 2011). The SOLID185 element was defined by eight nodes with three freedom degrees at each node. Some medium parameters used in the present study, such as Young’s modulus, Poisson’s ratio and density, were inferred from the tomographic model of Huang et al. (2009). These parameters are not constant but variable in North China. Young’s modulus range from 1.77 1010 to 1.61 1011 Pa, and Poisson’s ratios change from 0.25 to 0.39, while density varies between 2384 and 4100 kg/m3. All the parameters of the material properties are listed in Table 2. The classifications 1–85 of the material properties were assigned in the corresponding elements in the model, while those 86–95 were mainly used to reduce the element stiffness of the elements where the strong earthquake occurred. Note that the classification 107 of the material properties was assigned to all the faults, and that 108 was assigned to the mantle elements. The parameter values at each grid node can be obtained by interpolating those close to this node. Fig. 3 shows the distributions of Young’s modulus in the layers 1, 3, and 8 of the model. In order to examine the reasonability of all the parameters assigned in the model, the annual displacements at all grid nodes on the surface were obtained by applying the least square technique to the GPS observations in the three time intervals of 1999, 2001, and 2004. These annual displacement values on the surface were also assigned to those nodes that are in the same locations but at different depths. Using these constraints of the observed GPS velocity at the boundary of the North China sub-model, the displacement field at all grid nodes of the model was computed. Fig. 4 compares the surface displacements inferred from the GPS observations and calculated by using the finite element method. It is evident that the size and direction of the displacements at most grid nodes from these two ways are very close to each other, except for those in the Ordos interior. Because there is no GPS observations in the Ordos block, they are inferred by interpolating those in the Ordos margins. The mean amplitude difference between observed and calculated GPS values and its rms (the root mean square) are about 2.6 mm and 5.1 mm, respectively, while the mean angle difference and its rms are about 11° and 15°, respectively. All these results suggest that the element geometries and physical parameters assigned in our present finite element model of North China are generally reasonable. 3. The initial stress field of the North China model In order to exactly reveal the activity trend of strong earthquakes in North China, we cannot separate completely it from
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
119
Fig. 2. (a) Illustration of the top layer in the North China model (100 times vertical exaggeration). (b) Illustration of the bottom layer in the North China model (25 times vertical exaggeration).
the Chinese continent. In the 20th century there existed five periods of earthquake activities in the Chinese continent, and the strong earthquakes in North China were only concentrated in the fourth active period of 1966–1976 (Ma et al., 2003), which corresponds to the strong earthquake migrations in the North–South Seismic Zone. Lu et al. (2011) once simulated the strong earthquake occurrence applying the killed element technique in the entire Chinese continent, and suggested that the simulated earthquakes corresponded to each other in North China and the North–South Seismic Zone by simultaneously using the 1st–5th types of the border loadings. In the present study, although a stiffness reduction technique was used, instead of the killed element Table 2 Model settings. Material property
Young’s modulus (Pa)
Poisson’s ratio
Density (kg/m)
Viscosity coefficient (Pa S)
1–85
1.77– 9.01 1010 1 108 5 108 6 108 7 108 8 108 9 108 1 109 2 109 3 109 4 109 1.77 109 1 1011
0.25–0.39
2384– 3229 2300 2300 2300 2300 2300 2300 2300 2300 2300 2300 2700 3200
3 1019
86 87 88 89 90 91 92 93 94 95 107 108
0.35 0.34 0.33 0.32 0.31 0.30 0.29 0.28 0.27 0.26 0.30 0.35
method as shown in previous study of Lu et al. (2011), a similar result has been obtained (Fig. 5a). The earthquakes 1 and 2 in the southern part of the North–South seismic zone were followed by the earthquake 3 in North China. When the earthquakes 10 and 12 occurred in the southern part of the North–South seismic zone, subsequently the earthquake 13 occurred in North China. The earthquake 13 in North China was followed by the earthquakes 14 and 15 in the central part of the North–South seismic zone and the earthquake 17 in the southern part of North China. Fig. 5b and c illustrates the changes of the equivalent stress in the Chinese continent, when the earthquakes 2 and 13 occurred, respectively. It is clearly seen that the incoming distant earthquakes right occurred in the region of the increased equivalent stress. In order to better analyze the strong earthquake migration in North China, we adopted the sub-model technique of the software ANSYS 12.1 to establish a sub-model of North China within the model of the Chinese continent. The displacement at the sub-model boundary is that calculated from the entire Chinese continent model at the same place. Following the Saint–Venant principle, stress and strain only change near the loading location, if the actual loadings are replaced with equivalent loadings, which suggests that stress concentration only occurs where the loading is concentrated. Therefore, we conclude that more accurate results can be obtained in the sub-model when they are far from the stress concentration. In the Chinese continent model, the displacement at each node at the boundary of North China was calculated and then applied to the North China model as the boundary displacement loadings. As an example, Fig. 6a shows the displacement loading at the boundary nodes of the sub-model at 15 km depth. This
120
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 3. Distribution of Young’s modulus (Pa) in the first, third and eighth layers. Red and blue colors denote high and low modulus, respectively. The modulus scale is shown at the bottom.
Fig. 4. Comparison of the surface node displacements between the GPS least square fitting (red color) and finite element simulations (blue color).
displacement ranges from 10 to 24 mm. The numbers in Fig. 6a denote the boundaries of the sub-model. Using such displacement boundary loadings and taking into account the gravity effect, we calculated the initial stress field of the North China region, and simulated the strong earthquake occurrence and the stress adjustment induced. Fig. 6b shows the distribution of equivalent stress at 15–20 km depths. It can be seen that, in addition to the sub-model boundary, there are stress concentration at the end of the faults and a continuous high-stress region in the central part of Shanxi Province.
4. Relationship between simulated earthquakes and boundary loadings The initial stress field was applied to constrain the basic stress environment in North China within a certain period (here means an active period) so as to better investigate the migration of strong earthquakes and main active areas in North China. In order to achieve the goal, there are two issues required to be explored. One is to define the condition of strong earthquakes and how to simulate the disturbance of the displacement and stress and strain
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
121
Fig. 5. (a) Mutual migration of strong earthquakes from North China to the central and south parts of the North–south seismic zone under the 1st–5th Chinese continent boundary loadings, as shown in Table 2 of Lu et al. (2011). The numbers around solid dots denote the sequence codes of simulated earthquakes. (b) The equivalent stress distribution around the Chinese continent induced by the simulated earthquake 2. (c) The equivalent stress distribution around the Chinese continent induced by the simulated earthquake 13. Solid lines denote major active faults. Red and blue colors in (b) and (c) denote high and low stress. The stress scale is shown at the bottom of (b) and (c).
122
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 6. (a) The boundary displacement loadings (unit: m) at 15 km depth of the North China model. The numbers around the boundaries denote the boundary loading configurations as shown in Table 3. (b) Distribution of the simulated equivalent stress (unit: Pa) at 15–20 km depth. Red and blue colors denote high and low stress. The stress color scale is shown at the bottom.
fields, and the occurrence of following-up strong earthquakes after one strong earthquake, while the other is whether there are some significant effects on the distribution of strong earthquakes when the boundary loadings of the North China model are changed in a short time. Some researchers simulated the strong earthquakes with the highest equivalent stress groups and calculated the stress disturbance field by killing these highest equivalent stress groups (Lu et al., 2011), while others used the risk factor (G) of the fault to define the strong earthquake occurrence (Zhang et al., 2005). When G is larger than 1.0, the fault at the boundary of the block slips locally to generate the earthquake, while when G is less than 1.0, the fault does not slip and no earthquake occurs. In the present study, we estimated the location of strong earthquakes using the highest equivalent stress elements, the highest stress intensity elements, and the highest G value elements. Our
results show that the strong earthquake areas inferred from the above three methods are generally consistent with each other, but the high-G value areas are relatively smaller and more close to the location of strong earthquakes. Therefore, in this study we evaluated possible locations of strong earthquakes using the highest G value method. After a simulated earthquake occurred, the Young’s modulus of the corresponding element would lower and Poisson’s ratio would increase. The material property of the element group of the simulated strong earthquake was set to be the classification 86 (Table 2). Except for the decreased Young’s modulus, there were no any other changes in the model. In such a model the highest G value element was recalculated and was set to be a new simulating strong earthquake, where the material property was assigned to be the classification 86, and in the simulated strong earthquake element the material property was changed
123
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127 Table 3 Some representative boundary displacement loadings (Unit: mm). Boundary loading configuration
Loading direction
Boundary 1
1
Eastward Northward Eastward Northward Eastward Northward Eastward Northward Eastward Northward Eastward Northward
0.0025 0.002
2 3 4 5 6
0.0025 0.002 0.0025 0.002 0.0025 0.002 0.0025 0.002
Boundary 2
0.003 0.0025 0.003 0.0025 0.003 0.0025 0.003 0.0025 0.003 0.0025
from the classification 86–87, suggesting that Young’s modulus in the previous element were gradually recovered. Similarly, such simulation continued for 10 times because there were at least 10 strong earthquakes occurred in North China during the active period of 1966–1976 (Ma et al., 2003). The classification of the material property in the first group of the simulated strong earthquakes was changed to 95, and the one in the last group was set to 86. The classifications of the material property in other groups were between 94 and 87, respectively. Because we cannot know the trending and dip of the faults in North China before our simulation, all the faults were assumed to be vertical and a trending angle of NE 40° during the calculation of the G values. These assumptions were based on most strong earthquakes basically dominated by strike-slip and NE trending rupture in North China. All the friction coefficients were set to be 0.4 to better calculate the G value. However, to demonstrate the reasonability of our assumption, the case with the trending angle of NE 15° was also tested, but the results are similar. A variety of short-term displacement boundary loading modes were designed in North China (Table 3 and Fig. 7). The location criterion of simulated strong earthquakes was determined using the highest G value, and the follow-up stress adjustment induced by strong earthquakes was simulated using the element stiffness reduction. Through these techniques we explored the migration of the strong earthquakes in North China. During the simulation, under each loading configuration, after the 10 loop calculations we obtained 10 simulated strong earthquakes and the corresponding stress and strain perturbations. Our results show that, in the same initial stress state, most of the simulated strong earthquakes are concentrated on the seismic zones in North China Plain, and some are around the Ordos block, and a few are around the Tanlu seismic zone. We did not simulate the strong earthquakes in the Bohai Sea due to the scare of the reliable data there, such as active faults and seismic structure. Table 3 lists the six representative loading configurations, and Fig. 7 illustrates the distribution of the simulated strong earthquakes using these loading configurations. It can be seen that all the simulated earthquakes occurred near the fault zones. However, it is worth noting that the boundary loading mode may have significant impacts on the spatial distribution of simulated earthquakes. When there is only a unilateral boundary loading, such as the first and second loadings, most simulated earthquakes occur near the loading boundary (Fig. 7a and b). When the two boundary loadings are simultaneously considered, such as the first and second loadings or the third boundary loading, it is evident that the simulated strong earthquakes move eastward (Fig. 7c). When there are three or more than three boundary loadings, the simulated earthquakes are concentrated on various areas (Fig. 7d–f). For example, the forth loading can cause simulated strong earthquakes to concentrate on the North China Plain seismic zones and the southern part
Boundary 3
Boundary 4
0.002 0.0015 0.002 0.0015 0.001 0.001
Boundary 5
2E4 5E4
Boundary 6
Boundary 7
0.001 0.0003 0.001 0.0005
0.001 0.0005 0.001 0.0005 0.001 0.001
of the Tanlu zone (Fig. 7d). The fifth loading can result in the earthquakes to concentrate on the northern part of the North China Plain (Fig. 7e). The sixth loading can lead the earthquakes to the north and south regions of the North China Plain (Fig. 7f). Therefore, we conclude that if we can determine the boundary loading configuration in North China during a certain period of time, then it is possible to obtain the main active region of strong earthquakes using the present numerical simulation method. To further investigate whether the numerical simulation is useful to obtain some prediction instructions on the location of future strong earthquakes, the classification of the material property of the elements near the faults where the earthquakes (M > 6.0) occurred from 1900 to 1965 was set to be 107 to explore what boundary loadings can make the locations of simulated earthquakes closer to the actual earthquakes from 1966 to 1998. Fig. 8 shows the spatial distribution of the 10 simulated earthquakes under two kinds of different boundary loading configurations. The red circles denote the epicenters of strong earthquakes (M > 6.0) that occurred from 1966 to 1998. As seen from Fig. 8, there are 8–9 out of 10 simulated earthquakes that are close to the actual earthquake locations. Under the first loading configuration, there are 9 out of 10 simulated earthquakes near the hypocenters of the actual earthquakes (within the range of 110 km) (Fig. 8a). For example, the distances from the simulated earthquakes to the 1966 Xingtai, 1969 Bohai, 1976 Helingeer, 1976 Tangshan, 1983 Heze earthquakes are about 120 km, 120 km, 109 km, 62 km, and 150 km, respectively. Under the second loading configuration, there are 7 out of 10 simulated earthquakes close to the locations of the actual earthquakes (Fig. 8b), but they are closer to each other than those (less than 110 km) under the first kind of loading configuration. For example, the distances from the simulated earthquakes to the 1969 Bohai, 1976 Helingeer, 1976 Tangshan, 1983 Heze, and 1998 Hebei-Shangyi earthquakes are about 107 km, 109 km, 53 km, 108 km, and 108 km, respectively. These results indicate that useful instructions could be obtained for the earthquakes larger than 6.0 during 1966–1998 using numerical simulation that considered stiffness reduction near the historic earthquakes and joint boundary loadings from the Indian and Pacific plates over the past decades. However, it is worth noting that in the present study the friction coefficients and trending of the fault planes were set to be 0.4 and NE40°, respectively, which are not completely consistent with those actual earthquakes. This may explain why simulated earthquakes cannot correspond well to the actual earthquakes one by one. In order to explore whether there exist some agreements between simulated and actual earthquakes in time order, the sequence number of our simulated earthquakes in North China is marked in Fig. 7 and 8. It is found that there are some differences between simulated and actual earthquakes in time order, suggesting that stiffness reduction and gradual recovery have an
124
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 7. Spatial distribution of simulated earthquakes inferred from the sixth displacement boundary loadings (arrows) as shown in Table 3. Solid circles with 10 digital numbers denote simulated earthquakes, while open circles denote historic strong earthquakes. Black lines denote major active faults in North China.
accumulating influence on stress field and G values after a simulated earthquake. To further demonstrate if there do exist some inconsistences between simulated and actual earthquakes in time order, one more simulation test has been conducted. In this test, the element of the first simulated earthquake was assigned to locate in an actual earthquake, such as the 1966 Xingtai or 1975 Haicheng earthquake, and the results still show that some inconsistences between them in time order. These results suggest that
it is hard to estimate the time order of strong earthquakes in North China using the present numerical simulation method.
5. Possible areas of future strong earthquakes The location and sequence of simulated earthquakes were mainly controlled by the spatial distribution of stress perturbation
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
125
Fig. 8. Spatial distribution of simulated earthquakes by assuming the 107th material property around the earthquake hypocenters (M > 6.0) (open circles) from 1900 to 1965. The numbers around the arrows (boundary loading) denote the boundary displacements. Other symbols are the same as shown in Fig. 7.
induced by stiffness reduction and the G values inferred from the fault trending and friction coefficients, because the model and boundary conditions remained unchanged during the simulation. In the present study, the trending of the fault and friction coefficients were assumed to be constant, thus the stress field perturbation induced by stiffness reduction may be one of the main factors determining the location and time order of simulated earthquakes. Fig. 9a illustrates the distribution of the equivalent stress in a three-dimensional model, which was caused by the simulated earthquake 8 in the source area of the 1966 Xingtai earthquake with the boundary loading as shown in Fig. 8a. Fig. 9b shows the equivalent stress distribution of all the elements in the fourth layer (at 15–20 km depth) of our model. It is visible that the high stress disturbance area extends northeastward, and at the boundary of the highest G value where the simulated earthquake 9 occurred is 120 km away from the 1969 Bohai earthquake. Although this
stress disturbance is not the highest value, it has amounted to 1.2 MPa. Other simulated earthquakes show a similar feature. Table 4 lists the stress disturbance values around the next simulated earthquake induced by present earthquake. Except for the simulated earthquakes 1 and 3, the stress disturbance values induced by other eight earthquakes are all larger than 0.1 MPa. Note that the stress disturbance value at the location of the simulated earthquake 2 induced by the earthquake 1 was very low, but the earthquake 2 was located at the end of the NE 40° fault and the G value was greater than 1.0. Therefore, it is concluded that the possible region could be predicted by the highest stress disturbance field, although the trending and friction state of the fault zones were not given exactly. In order to investigate the variations of the equivalent stress with depth, we plotted two vertical cross sections 13 and 23 passing through the simulated earthquakes 8 and 9, respectively
126
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 9. (a) Cubed view of the equivalent stress (unit: Pa) induced by the 8th simulated earthquake exaggerated in depth. (b) Map view of the equivalent stress (unit: Pa) at 15– 20 km depth. (c and d) The same as (a) but for two vertical cross sections 13 and 23 as shown in (a), respectively Red and blue colors denote high and low stress. The stress perturbation scale is shown on the left. Open circles denote the simulated earthquakes 8 and 9.
Table 4 Stress perturbation (Unit: MPa). Simulated earthquake (i)
1
2
3
4
5
6
7
8
9
Stress perturbation around the simulated earthquake (i + 1) induced by the earthquake (i)
0.0056
0.59
0.087
0.95
0.17
0.19
2.2
1.2
2.2
(Fig. 9c and d). This stress state is caused by the simulated earthquake 8. It is evident that all the equivalent stress increases from the southwest to the northeast along these two cross sections, suggesting that strong earthquakes would move northeastward. In the cross section 13 passing through the earthquake 8, the
equivalent stress values around the source area become higher from shallow to deep depth, and are relatively smaller above 20 km depth in the southwest and northeast segments, which may explain why many small earthquakes occurred there. In the cross section 23 passing through the simulated earthquake 9, the
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
source areas extend northeastward and are in a high stress state from the shallow to deep depth, which may be a possible stress state after the earthquake 8.
127
instructive comments and suggestions, which improved our manuscript. References
6. Discussion and conclusions In this study a three-dimensional Maxwell viscoelastic model was established for North China using previous geological, geophysical and geodetic observations. In the model, the faults were set to be the weakened material. The initial stress field of North China was calculated using the finite element sub-model technique and considering the gravity and the boundary loadings from the Chinese continent model. Based on this initial model, we considered a variety of displacement boundary loading configurations to obtain the spatial distributions of 10 simulated earthquakes from the spatially distributed stress field and G values. After examining and comparing the displacements from our present initial model and the GPS observations, it is inferred that the initial model was considered to be reasonable. The simulated earthquakes are close to the actual ones and they are about 100 km apart from each other. The number of such earthquakes can amount to over 50%. Therefore, we conclude that the possible region of follow-up strong earthquakes could be predicted by obtaining the location of high stress disturbance values. Our results also show that it is critical to choose the approximate boundary loadings around the Chinese continent so as to obtain more reasonable parameters for North China, suggesting that the strong earthquake prediction in North China could depend in certain extent on the stress field and strong earthquake activities in the Chinese continent. The time sequence of simulated earthquakes can hardly correspond to that of actual ones, maybe due to the unknown trending and friction coefficient of the faults used to calculate the G value. Therefore, we infer that the main active region of strong earthquakes could be predicted using the finite element numerical simulation, but it is still hard to conduct an exact dynamic prediction due to the lack of a better understanding of the medium parameters, including the fault property and the crust and upper mantle structure. Acknowledgments We thank Academician Q. Deng, Profs. Z. Jiang and Z. Huang for providing the active faults, GPS velocity field and seismic velocity model in North China, respectively. This work was partially supported by the Exploration Technology Deep and Experimental Study (SinoProbe-06-03) and Special grants of Basic Scientific and Research (ZDJ2007-1 and ZDJ2009-1) and National Natural Science Foundations (40774044 and 40974021). We thank Prof. Bor-ming Jahn and two anonymous reviewers for providing many
Bai, W., Lin, B., Chen, Z., 2003. Numerical simulations on the block movement and deformation in North China after the 1976 Tangshan earthquake. Science in China (Series D) 33 (Suppl.), 99–107. Chen, L., Lu, Y., Liu, J., et al., 2001. Three-dimensional viscoelastic numerical simulation on dynamic evolution of stress field in North China induced by the 1966 Xingtai earthquake. Acta Seismologica Sinica 23 (5), 480–491. Chen, X., Zhang, J., Tang, R., et al., 2002. The Distribution of Moho Discontinuity in China and Surrounding Areas. Earthquake Press, Beijing, p. 1. Deng, Q., Zhang, P., Ran, Y., et al., 2002. Basic characteristics on the active tectonics in the Chinese continent. Science in China (Series D) 12, 1020–1030. Ding, J., Liu, J., Yu, S., 2000. Exploration and practice of earthquake prediction in China. Earthquake 20 (Suppl.), 13–17. Han, Z., Xu, J., Ran, Y., et al., 2003. Active blocks and strong earthquake activities in North China. Science in China (Series D) 33 (Suppl.), 108–118. Huang, Z., Li, H., Zheng, Y., Peng, Y., 2009. The lithosphere of North China Craton from surface wave tomography. Earth and Planetary Science Letters 288, 164– 173. Jiang, Z., Zhang, X., Chen, B., Xue, F., 2000. Characteristics of recent horizontal movement and strain-stress field in the crust of North China. Chinese Journal of Geophysics 43 (5), 657–665. Kenner, S.J., Segall, P., 2000. A mechanical model for intraplate earthquakes, application to the New Madrid seismic zone. Science 289, 2329–2332. Lei, J., Xie, F., Lan, C., Xing, C., Ma, S., 2008. Seismic images under the Beijing region inferred from P and PmP data. Physics of the Earth and Planetary Interiors 168, 134–146. Lei, J., Zhao, D., Xie, F., Liu, J., 2011. An attempt to detect temporal variations of crustal structure in the source area of the 2006 Wen—an earthquake in North China. Journal of Asian Earth Sciences 40, 958–976. Li, Q., Liu, M., Stein, S., 2009. Spatial–temporal complexity of continental intraplate seismicity: insights from geodynamic modeling and implications for seismic hazard estimation. Bulletin of the Seismological Society of America 99, 52–60. Lu, Y., 2004. Suggestions on promoting earthquake prediction in China. Recent Development in World Seismology 5, 81–84. Lu, Y., Yang, S., Chen, L., Lei, J., 2011. Mechanism of spatial distribution and migration of strong earthquakes in China inferred from numerical simulation. Journal of Asian Earth Sciences 40, 990–1001. Ma, H., Zhang, G., Liu, J., et al., 2003. Correlation between earthquake activity and activity crustal-block in China mainland and its adjacent regions. Earthquake Science Frontiers 10 (Suppl.), 75–80. Parsons, T., 2002. Post-1906 stress recovery of the San Andreas fault system calculated from three-dimensional finite element analysis. Journal of Geophysical Research 107, 1–13. Seyferth, M., Henk, A., 2002. Coupling of PFC2D and ANSYS – concepts to combine the best of two worlds for improved geodynamic models. In: Konietzky, H. (Ed.), Numerical Modeling in Micromechanics Via Particle Methods. Balkema, pp. 283–290. Tao, W., Hong, H., Liu, P., 2000. The numerical simulation of the strong earthquake region around the Chinese continent and surrounding regions. Acta Seismologica Sinica 22 (3), 271–277. Wang, R., Sun, X., Cai, Y., 1982. Numerical simulation of earthquake sequence over the past 700 years. Science in China (Series B) 8, 745–753. Zhang, Z., Deng, Y., Wang, B., 1990. Joint inversion of fracture system models of Ms P 7.0 earthquake sequences and rupture processes in North China. Acta Seismologica Sinica 12 (4), 335–347. Zhang, R., Wei, F., Qiao, C., et al., 2005. Numerical simulations for the 1975 Haicheng and 1999 Xiuyan earthquake processes by using the DDA + FEM method. Acta Seismologic Sinica 27 (3), 163–170. Zhu, Y., Wang, Q., Zhao, X., et al., 2004. Simulation of welding residual stress based on ANSYS. Journal of Wuhan University 26 (2), 69–72.
Contents lists available at SciVerse ScienceDirect
Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes
Migration trend of strong earthquakes in North China from numerical simulations Yuanzhong Lu a, Shuxin Yang a, Lianwang Chen a, Jianshe Lei a,⇑, Ping He b a b
Key Laboratory of Crustal Dynamics, Institute of Crustal Dynamics, CEA, Beijing 100085, China Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China
a r t i c l e
i n f o
Article history: Received 12 August 2011 Received in revised form 2 January 2012 Accepted 14 January 2012 Available online 17 February 2012 Keywords: Stiffness reduction Numerical simulation Distant migration Strong earthquake North China
a b s t r a c t Lu et al. (2011) established a finite element model around the entire Chinese continent to simulate the distant migration of strong earthquakes by killing the elements. In this study a new finite element model of North China has been established in the existing Chinese continental model using the GPS observations and the finite element sub-model technique through taking into account geological structure, active faults, surface topography, Moho discontinuity and 3-D velocity structure of the crust and upper mantle in the region. Furthermore, the stress field adjustment was realized by reducing the element stiffness rather than by killing the elements. We calculated the displacement of the Chinese continent model, and took the displacement located at the boundary of North China as the boundary loadings, then computed the initial stress field by considering the gravity effect. The migration of the strong earthquake regions was simulated by performing a number of tests with various kinds of short-term boundary loadings. Locations of strong earthquakes were determined with the highest G (the risk factor of the fault) value and adjustment of the stress field. Our results showed that high stress disturbance regions could be possible regions of future strong earthquakes, and that the initial stress field in North China can be well estimated, if an optimal boundary loading configuration is chosen in the Chinese continental model. These results suggest that the strong earthquake forecast in North China may strongly depend on the stress field and strong earthquake activities in a much wider area, such as the Chinese continent. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Through analyzing earthquake activities and precursory phenomena, including crustal stress and strain, geomagnetical, geoelectrical anomalies, and ground fluids, a few moderate-strong earthquakes were once successfully predicted and natural hazards were mitigated significantly (Ding et al., 2000). However, the failure of the prediction to some destructive earthquakes, such as the 1976 Tangshan (M 7.8) and 2008 Wenchuan (M 8.0), China, earthquakes, renders us to ponder the feasibility of earthquake prediction and related scientific investigations. A common consensus is that due to the complexity and nonlinearity of seismogenic structures, it is extremely hard to obtain useful information on the earthquake pregnancy and occurrence, and nowadays more difficult to predict the exact location, magnitude, and origin time. It is almost impossible to achieve a significant breakthrough in the earthquake prediction if we only depend on precursor observations. One wants to know whether there exists a possibility to predict the trend of strong earthquakes based on the deficient understanding of crustal structure, the property of fault zones, and crustal movements. Wang et al. (1982) once simulated strong earthquakes (M P 7.0) over the past 700 years, and found that most ⇑ Corresponding author. Tel./fax: +86 10 6284 6760. E-mail addresses: [email protected], [email protected] (J. Lei). 1367-9120/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.jseaes.2012.01.006
earthquakes occurred around the region with an increased value of the risk factor G = s=ðrn tan /Þ, where s is the shear stress along the fault plane, rn is the normal stress on the fault plane, and / is the friction angle. Zhang et al. (1990) jointly inverted the 1966 Xingtai (M 7.2), 1969 Bohai (M 7.4), 1975 Haicheng (M 7.3), 1976 Tangshan (M 7.8), China, earthquakes, and found that the occurrence of strong earthquakes could be resulted from the activity and extension of the fault system. Tao et al. (2000) established a viscoelastic finite element model in the Chinese continent, and obtained the distribution of strain and stress field. Chen et al. (2001) investigated a dynamic evolution process of the 1966 Xingtai earthquake, and concluded that the occurrence of the 1966 Xingtai earthquake could play a key role in the stress variation around the 1978 Tangshan source area. Bai et al. (2003) simulated the dynamic process of the Tangshan earthquake, and found that the earthquake significantly affected its adjacent blocks and border faults and led to various sizes of ground movement and surface deformation around it. Based on these studies, Lu (2004) proposed a seismic dynamic prediction model attributable to a remarkable process in the computer technique, geodynamics, deep seismic structure, GPS observations, and continental deep drillings. In the seismic dynamic prediction model, we do not only take into account one or several faults and the interaction between the block and its surrounding fault system, but also dynamic processes of stress concentration and fault slipping in the brittle and
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
brittle-ductile mid-upper crust away from the plate boundaries. Assuming a sudden weakening of a localized weak lower crust embedded in an elastic crust, Kenner and Segall (2000) concluded that weakening in the lower crust can causes stress amplification in the upper crust and result in a sequence of ruptures of the fault zone. Li et al. (2009) investigated different cases of fault weakening after a large fault rupture, and concluded that fault weakening can lead to repeated earthquakes. Lu et al. (2011) made an attempt to establish a finite element model of the entire Chinese continent (Fig. 1a) by killing the elements (no loading capability) after the earthquake, but this method cannot predict where the next earthquake occurs. However, the purpose of the dynamic prediction at the current stage is to forecast the migration trend of an incoming strong earthquake region so as to provide an important scientific evidence for capturing the seismic precursors to mitigate the seismic hazard.
117
Recently, in North China many investigations have being conducted using various approaches, such as seismic tomography (Lei et al., 2008, 2011; Huang et al., 2009), active faults (Deng et al., 2002), and GPS observations (Jiang et al., 2000), and obtained many significant results. These results allow us to model the migration trend of the strong earthquake region in North China by introducing these previous results into our model. To achieve practical and meaningful forecast of strong earthquakes with the purpose of the mitigation of seismic hazard, the region of the predicted earthquakes should not be too large. Therefore, in the present study we took North China as an example to explore the migration trend of follow-up strong earthquakes. Because the loading configuration could not be changed significantly within an active period of about 10 years, the migration of strong earthquakes could be mainly controlled by the source areas without load-bearing capacity and induced stress adjustment. Our simulation steps
Fig. 1. (a) Location of the North China sub-model (open polygon) in the entire Chinese continent finite element model of Lu et al. (2011). Black lines denote major active faults in China, while other lines denote the model parameterization. The digital numbers denote the model boundary. (b) A finite element model of North China. (c) Distribution of the faults used in the finite element model of North China. The numbers denote the faults as shown Table 1.
118
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Table 1 Some parameters of the faults (Deng et al., 2002) used in the North China sub-model. Sequence number
Name
Type Dip-slip Unknown Dip-slip Strike-slip
7 8 9 10 11
Helan Mountain Eastern Foot Fault Luo Mountain Eastern Foot Fault Qinling Mountain Northern Foot Fault Huashan Mountain Piedmont Fault Wentang Fault Zhongtiaoshan Mountain Piedmont Fault Luoyun Mountain Piedmont Fault Huo Mountain Piedmont Fault Shangu Fault Jiaocheng Fault Jizhou Mountain Piedmont Fault
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
Wutai Mountain Western Foot Fault Heng Mountain Northern Foot Fault Kouquan Fault Tangxi Fault Tangdong Fault Handan Fault Yuanshi Fault Zhangjiakou Fault Xinbaoan-Shacheng Fault Shizhuang Fault Sunhe-Nankou Fault Jiyunhe Fault Anyang Fault Heze Fault Cangni Fault Tangshan Fault Cangdong Fault Xiadian Fault Sunzhuangzi-Wulonggou Fault Liaocheng-lankao Fault Cixian-daming Fault Dongming-chengwu Fault Xinhe Fault Huanghe-Lingwu Fault Qinmu Fault Ordos Fault Pikou Fault Jinzhou Fault Gaixian Fault Liaodong Fault
1 2 3 4 5 6
Dip-slip Unkown Dip-slip Strike-slip Dip-slip Strike-slip Dip-slip Dip-slip Unknown Right-lateral strikeslip Dip-slip Dip-slip Strike-slip Unknown Unknown Unknown Unknown Unknown Unknown Unknown Strike-slip Unknown Unknown Unknown Unknown Unknown Strike-slip Dip-slip Unknown Unknown Unknown Unknown Unknown Dip-slip Strike-slip Dip-slip Strike-slip Dip-slip Unknown Unknown Unknown Unknown
are as follows. Firstly, a finite element model of North China (Fig. 1) was established based on the entire Chinese continental model (Fig. 1a, Lu et al., 2011). In the Chinese continental model they did not take into account the undulated Moho discontinuity and faults, while in present study we added both the undulated Moho discontinuity and 41 faults in the North China sub-model to obtain more accurate results. In the sub-model the stiffness diminution by changing the element’s material properties including Young’s modulus and Poisson’s ratio within the faults was designed after the occurrence of strong earthquakes so as to calculate a new stress field. Moreover, the same amount of stiffness is reduced for all the earthquake faults and the stress adjustment is time independent. Then the possible criterion was explored to determine the occurrence of subsequent simulated earthquakes. Finally, the migration trend of strong earthquakes was estimated after a multi-step cycle calculation. 2. Establishment of the finite element model in North China Lu et al. (2011) set up a finite element model in the entire Chinese continent. In order to better investigate the mitigation of strong earthquakes in North China, here we added 41 active faults with an assumed 2 km width in the North China sub-model. The
location and trending of these faults were extracted from Deng et al. (2002) (Table 1 and Fig. 1c). The boundary of the North China model (Fig. 1a) was directly adopted from Han et al. (2003). In the model the 98,433 grid nodes and 86,840 elements were set up. The elements on the faults are dense, where the grid spacing is about 2 km, while those away from the faults are relatively sparse, where the grid spacing is about 60 km. The material on the faults is very weaker than that away from the fault. The model was divided into eight layers. The first layer is from 5 km depth to the ground surface (the topography from www.usgs.gov) (Fig. 2a), while the bottom layer is between 45 km and the Moho discontinuity (Chen et al., 2002) (Fig. 2b). Other layers are evenly distributed between the first and bottom layers. The material is assumed to be viscoelastic at the bottom layer of the mantle, while those at other layers are assumed to be linear elastic. The commercial software ANSYS 12.1 was used to establish the model with SOLID185 hexahedral elements and finish some related calculations. The ANSYS software is suitable for numerical simulations in a variety of physical problems, such as solid mechanics, thermal, electromagnetic and fluid dynamics, but recently it is popularly used in Earth sciences (e.g., Parsons, 2002; Seyferth and Henk, 2002; Zhu et al., 2004; Lu et al., 2011). The SOLID185 element was defined by eight nodes with three freedom degrees at each node. Some medium parameters used in the present study, such as Young’s modulus, Poisson’s ratio and density, were inferred from the tomographic model of Huang et al. (2009). These parameters are not constant but variable in North China. Young’s modulus range from 1.77 1010 to 1.61 1011 Pa, and Poisson’s ratios change from 0.25 to 0.39, while density varies between 2384 and 4100 kg/m3. All the parameters of the material properties are listed in Table 2. The classifications 1–85 of the material properties were assigned in the corresponding elements in the model, while those 86–95 were mainly used to reduce the element stiffness of the elements where the strong earthquake occurred. Note that the classification 107 of the material properties was assigned to all the faults, and that 108 was assigned to the mantle elements. The parameter values at each grid node can be obtained by interpolating those close to this node. Fig. 3 shows the distributions of Young’s modulus in the layers 1, 3, and 8 of the model. In order to examine the reasonability of all the parameters assigned in the model, the annual displacements at all grid nodes on the surface were obtained by applying the least square technique to the GPS observations in the three time intervals of 1999, 2001, and 2004. These annual displacement values on the surface were also assigned to those nodes that are in the same locations but at different depths. Using these constraints of the observed GPS velocity at the boundary of the North China sub-model, the displacement field at all grid nodes of the model was computed. Fig. 4 compares the surface displacements inferred from the GPS observations and calculated by using the finite element method. It is evident that the size and direction of the displacements at most grid nodes from these two ways are very close to each other, except for those in the Ordos interior. Because there is no GPS observations in the Ordos block, they are inferred by interpolating those in the Ordos margins. The mean amplitude difference between observed and calculated GPS values and its rms (the root mean square) are about 2.6 mm and 5.1 mm, respectively, while the mean angle difference and its rms are about 11° and 15°, respectively. All these results suggest that the element geometries and physical parameters assigned in our present finite element model of North China are generally reasonable. 3. The initial stress field of the North China model In order to exactly reveal the activity trend of strong earthquakes in North China, we cannot separate completely it from
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
119
Fig. 2. (a) Illustration of the top layer in the North China model (100 times vertical exaggeration). (b) Illustration of the bottom layer in the North China model (25 times vertical exaggeration).
the Chinese continent. In the 20th century there existed five periods of earthquake activities in the Chinese continent, and the strong earthquakes in North China were only concentrated in the fourth active period of 1966–1976 (Ma et al., 2003), which corresponds to the strong earthquake migrations in the North–South Seismic Zone. Lu et al. (2011) once simulated the strong earthquake occurrence applying the killed element technique in the entire Chinese continent, and suggested that the simulated earthquakes corresponded to each other in North China and the North–South Seismic Zone by simultaneously using the 1st–5th types of the border loadings. In the present study, although a stiffness reduction technique was used, instead of the killed element Table 2 Model settings. Material property
Young’s modulus (Pa)
Poisson’s ratio
Density (kg/m)
Viscosity coefficient (Pa S)
1–85
1.77– 9.01 1010 1 108 5 108 6 108 7 108 8 108 9 108 1 109 2 109 3 109 4 109 1.77 109 1 1011
0.25–0.39
2384– 3229 2300 2300 2300 2300 2300 2300 2300 2300 2300 2300 2700 3200
3 1019
86 87 88 89 90 91 92 93 94 95 107 108
0.35 0.34 0.33 0.32 0.31 0.30 0.29 0.28 0.27 0.26 0.30 0.35
method as shown in previous study of Lu et al. (2011), a similar result has been obtained (Fig. 5a). The earthquakes 1 and 2 in the southern part of the North–South seismic zone were followed by the earthquake 3 in North China. When the earthquakes 10 and 12 occurred in the southern part of the North–South seismic zone, subsequently the earthquake 13 occurred in North China. The earthquake 13 in North China was followed by the earthquakes 14 and 15 in the central part of the North–South seismic zone and the earthquake 17 in the southern part of North China. Fig. 5b and c illustrates the changes of the equivalent stress in the Chinese continent, when the earthquakes 2 and 13 occurred, respectively. It is clearly seen that the incoming distant earthquakes right occurred in the region of the increased equivalent stress. In order to better analyze the strong earthquake migration in North China, we adopted the sub-model technique of the software ANSYS 12.1 to establish a sub-model of North China within the model of the Chinese continent. The displacement at the sub-model boundary is that calculated from the entire Chinese continent model at the same place. Following the Saint–Venant principle, stress and strain only change near the loading location, if the actual loadings are replaced with equivalent loadings, which suggests that stress concentration only occurs where the loading is concentrated. Therefore, we conclude that more accurate results can be obtained in the sub-model when they are far from the stress concentration. In the Chinese continent model, the displacement at each node at the boundary of North China was calculated and then applied to the North China model as the boundary displacement loadings. As an example, Fig. 6a shows the displacement loading at the boundary nodes of the sub-model at 15 km depth. This
120
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 3. Distribution of Young’s modulus (Pa) in the first, third and eighth layers. Red and blue colors denote high and low modulus, respectively. The modulus scale is shown at the bottom.
Fig. 4. Comparison of the surface node displacements between the GPS least square fitting (red color) and finite element simulations (blue color).
displacement ranges from 10 to 24 mm. The numbers in Fig. 6a denote the boundaries of the sub-model. Using such displacement boundary loadings and taking into account the gravity effect, we calculated the initial stress field of the North China region, and simulated the strong earthquake occurrence and the stress adjustment induced. Fig. 6b shows the distribution of equivalent stress at 15–20 km depths. It can be seen that, in addition to the sub-model boundary, there are stress concentration at the end of the faults and a continuous high-stress region in the central part of Shanxi Province.
4. Relationship between simulated earthquakes and boundary loadings The initial stress field was applied to constrain the basic stress environment in North China within a certain period (here means an active period) so as to better investigate the migration of strong earthquakes and main active areas in North China. In order to achieve the goal, there are two issues required to be explored. One is to define the condition of strong earthquakes and how to simulate the disturbance of the displacement and stress and strain
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
121
Fig. 5. (a) Mutual migration of strong earthquakes from North China to the central and south parts of the North–south seismic zone under the 1st–5th Chinese continent boundary loadings, as shown in Table 2 of Lu et al. (2011). The numbers around solid dots denote the sequence codes of simulated earthquakes. (b) The equivalent stress distribution around the Chinese continent induced by the simulated earthquake 2. (c) The equivalent stress distribution around the Chinese continent induced by the simulated earthquake 13. Solid lines denote major active faults. Red and blue colors in (b) and (c) denote high and low stress. The stress scale is shown at the bottom of (b) and (c).
122
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 6. (a) The boundary displacement loadings (unit: m) at 15 km depth of the North China model. The numbers around the boundaries denote the boundary loading configurations as shown in Table 3. (b) Distribution of the simulated equivalent stress (unit: Pa) at 15–20 km depth. Red and blue colors denote high and low stress. The stress color scale is shown at the bottom.
fields, and the occurrence of following-up strong earthquakes after one strong earthquake, while the other is whether there are some significant effects on the distribution of strong earthquakes when the boundary loadings of the North China model are changed in a short time. Some researchers simulated the strong earthquakes with the highest equivalent stress groups and calculated the stress disturbance field by killing these highest equivalent stress groups (Lu et al., 2011), while others used the risk factor (G) of the fault to define the strong earthquake occurrence (Zhang et al., 2005). When G is larger than 1.0, the fault at the boundary of the block slips locally to generate the earthquake, while when G is less than 1.0, the fault does not slip and no earthquake occurs. In the present study, we estimated the location of strong earthquakes using the highest equivalent stress elements, the highest stress intensity elements, and the highest G value elements. Our
results show that the strong earthquake areas inferred from the above three methods are generally consistent with each other, but the high-G value areas are relatively smaller and more close to the location of strong earthquakes. Therefore, in this study we evaluated possible locations of strong earthquakes using the highest G value method. After a simulated earthquake occurred, the Young’s modulus of the corresponding element would lower and Poisson’s ratio would increase. The material property of the element group of the simulated strong earthquake was set to be the classification 86 (Table 2). Except for the decreased Young’s modulus, there were no any other changes in the model. In such a model the highest G value element was recalculated and was set to be a new simulating strong earthquake, where the material property was assigned to be the classification 86, and in the simulated strong earthquake element the material property was changed
123
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127 Table 3 Some representative boundary displacement loadings (Unit: mm). Boundary loading configuration
Loading direction
Boundary 1
1
Eastward Northward Eastward Northward Eastward Northward Eastward Northward Eastward Northward Eastward Northward
0.0025 0.002
2 3 4 5 6
0.0025 0.002 0.0025 0.002 0.0025 0.002 0.0025 0.002
Boundary 2
0.003 0.0025 0.003 0.0025 0.003 0.0025 0.003 0.0025 0.003 0.0025
from the classification 86–87, suggesting that Young’s modulus in the previous element were gradually recovered. Similarly, such simulation continued for 10 times because there were at least 10 strong earthquakes occurred in North China during the active period of 1966–1976 (Ma et al., 2003). The classification of the material property in the first group of the simulated strong earthquakes was changed to 95, and the one in the last group was set to 86. The classifications of the material property in other groups were between 94 and 87, respectively. Because we cannot know the trending and dip of the faults in North China before our simulation, all the faults were assumed to be vertical and a trending angle of NE 40° during the calculation of the G values. These assumptions were based on most strong earthquakes basically dominated by strike-slip and NE trending rupture in North China. All the friction coefficients were set to be 0.4 to better calculate the G value. However, to demonstrate the reasonability of our assumption, the case with the trending angle of NE 15° was also tested, but the results are similar. A variety of short-term displacement boundary loading modes were designed in North China (Table 3 and Fig. 7). The location criterion of simulated strong earthquakes was determined using the highest G value, and the follow-up stress adjustment induced by strong earthquakes was simulated using the element stiffness reduction. Through these techniques we explored the migration of the strong earthquakes in North China. During the simulation, under each loading configuration, after the 10 loop calculations we obtained 10 simulated strong earthquakes and the corresponding stress and strain perturbations. Our results show that, in the same initial stress state, most of the simulated strong earthquakes are concentrated on the seismic zones in North China Plain, and some are around the Ordos block, and a few are around the Tanlu seismic zone. We did not simulate the strong earthquakes in the Bohai Sea due to the scare of the reliable data there, such as active faults and seismic structure. Table 3 lists the six representative loading configurations, and Fig. 7 illustrates the distribution of the simulated strong earthquakes using these loading configurations. It can be seen that all the simulated earthquakes occurred near the fault zones. However, it is worth noting that the boundary loading mode may have significant impacts on the spatial distribution of simulated earthquakes. When there is only a unilateral boundary loading, such as the first and second loadings, most simulated earthquakes occur near the loading boundary (Fig. 7a and b). When the two boundary loadings are simultaneously considered, such as the first and second loadings or the third boundary loading, it is evident that the simulated strong earthquakes move eastward (Fig. 7c). When there are three or more than three boundary loadings, the simulated earthquakes are concentrated on various areas (Fig. 7d–f). For example, the forth loading can cause simulated strong earthquakes to concentrate on the North China Plain seismic zones and the southern part
Boundary 3
Boundary 4
0.002 0.0015 0.002 0.0015 0.001 0.001
Boundary 5
2E4 5E4
Boundary 6
Boundary 7
0.001 0.0003 0.001 0.0005
0.001 0.0005 0.001 0.0005 0.001 0.001
of the Tanlu zone (Fig. 7d). The fifth loading can result in the earthquakes to concentrate on the northern part of the North China Plain (Fig. 7e). The sixth loading can lead the earthquakes to the north and south regions of the North China Plain (Fig. 7f). Therefore, we conclude that if we can determine the boundary loading configuration in North China during a certain period of time, then it is possible to obtain the main active region of strong earthquakes using the present numerical simulation method. To further investigate whether the numerical simulation is useful to obtain some prediction instructions on the location of future strong earthquakes, the classification of the material property of the elements near the faults where the earthquakes (M > 6.0) occurred from 1900 to 1965 was set to be 107 to explore what boundary loadings can make the locations of simulated earthquakes closer to the actual earthquakes from 1966 to 1998. Fig. 8 shows the spatial distribution of the 10 simulated earthquakes under two kinds of different boundary loading configurations. The red circles denote the epicenters of strong earthquakes (M > 6.0) that occurred from 1966 to 1998. As seen from Fig. 8, there are 8–9 out of 10 simulated earthquakes that are close to the actual earthquake locations. Under the first loading configuration, there are 9 out of 10 simulated earthquakes near the hypocenters of the actual earthquakes (within the range of 110 km) (Fig. 8a). For example, the distances from the simulated earthquakes to the 1966 Xingtai, 1969 Bohai, 1976 Helingeer, 1976 Tangshan, 1983 Heze earthquakes are about 120 km, 120 km, 109 km, 62 km, and 150 km, respectively. Under the second loading configuration, there are 7 out of 10 simulated earthquakes close to the locations of the actual earthquakes (Fig. 8b), but they are closer to each other than those (less than 110 km) under the first kind of loading configuration. For example, the distances from the simulated earthquakes to the 1969 Bohai, 1976 Helingeer, 1976 Tangshan, 1983 Heze, and 1998 Hebei-Shangyi earthquakes are about 107 km, 109 km, 53 km, 108 km, and 108 km, respectively. These results indicate that useful instructions could be obtained for the earthquakes larger than 6.0 during 1966–1998 using numerical simulation that considered stiffness reduction near the historic earthquakes and joint boundary loadings from the Indian and Pacific plates over the past decades. However, it is worth noting that in the present study the friction coefficients and trending of the fault planes were set to be 0.4 and NE40°, respectively, which are not completely consistent with those actual earthquakes. This may explain why simulated earthquakes cannot correspond well to the actual earthquakes one by one. In order to explore whether there exist some agreements between simulated and actual earthquakes in time order, the sequence number of our simulated earthquakes in North China is marked in Fig. 7 and 8. It is found that there are some differences between simulated and actual earthquakes in time order, suggesting that stiffness reduction and gradual recovery have an
124
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 7. Spatial distribution of simulated earthquakes inferred from the sixth displacement boundary loadings (arrows) as shown in Table 3. Solid circles with 10 digital numbers denote simulated earthquakes, while open circles denote historic strong earthquakes. Black lines denote major active faults in North China.
accumulating influence on stress field and G values after a simulated earthquake. To further demonstrate if there do exist some inconsistences between simulated and actual earthquakes in time order, one more simulation test has been conducted. In this test, the element of the first simulated earthquake was assigned to locate in an actual earthquake, such as the 1966 Xingtai or 1975 Haicheng earthquake, and the results still show that some inconsistences between them in time order. These results suggest that
it is hard to estimate the time order of strong earthquakes in North China using the present numerical simulation method.
5. Possible areas of future strong earthquakes The location and sequence of simulated earthquakes were mainly controlled by the spatial distribution of stress perturbation
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
125
Fig. 8. Spatial distribution of simulated earthquakes by assuming the 107th material property around the earthquake hypocenters (M > 6.0) (open circles) from 1900 to 1965. The numbers around the arrows (boundary loading) denote the boundary displacements. Other symbols are the same as shown in Fig. 7.
induced by stiffness reduction and the G values inferred from the fault trending and friction coefficients, because the model and boundary conditions remained unchanged during the simulation. In the present study, the trending of the fault and friction coefficients were assumed to be constant, thus the stress field perturbation induced by stiffness reduction may be one of the main factors determining the location and time order of simulated earthquakes. Fig. 9a illustrates the distribution of the equivalent stress in a three-dimensional model, which was caused by the simulated earthquake 8 in the source area of the 1966 Xingtai earthquake with the boundary loading as shown in Fig. 8a. Fig. 9b shows the equivalent stress distribution of all the elements in the fourth layer (at 15–20 km depth) of our model. It is visible that the high stress disturbance area extends northeastward, and at the boundary of the highest G value where the simulated earthquake 9 occurred is 120 km away from the 1969 Bohai earthquake. Although this
stress disturbance is not the highest value, it has amounted to 1.2 MPa. Other simulated earthquakes show a similar feature. Table 4 lists the stress disturbance values around the next simulated earthquake induced by present earthquake. Except for the simulated earthquakes 1 and 3, the stress disturbance values induced by other eight earthquakes are all larger than 0.1 MPa. Note that the stress disturbance value at the location of the simulated earthquake 2 induced by the earthquake 1 was very low, but the earthquake 2 was located at the end of the NE 40° fault and the G value was greater than 1.0. Therefore, it is concluded that the possible region could be predicted by the highest stress disturbance field, although the trending and friction state of the fault zones were not given exactly. In order to investigate the variations of the equivalent stress with depth, we plotted two vertical cross sections 13 and 23 passing through the simulated earthquakes 8 and 9, respectively
126
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
Fig. 9. (a) Cubed view of the equivalent stress (unit: Pa) induced by the 8th simulated earthquake exaggerated in depth. (b) Map view of the equivalent stress (unit: Pa) at 15– 20 km depth. (c and d) The same as (a) but for two vertical cross sections 13 and 23 as shown in (a), respectively Red and blue colors denote high and low stress. The stress perturbation scale is shown on the left. Open circles denote the simulated earthquakes 8 and 9.
Table 4 Stress perturbation (Unit: MPa). Simulated earthquake (i)
1
2
3
4
5
6
7
8
9
Stress perturbation around the simulated earthquake (i + 1) induced by the earthquake (i)
0.0056
0.59
0.087
0.95
0.17
0.19
2.2
1.2
2.2
(Fig. 9c and d). This stress state is caused by the simulated earthquake 8. It is evident that all the equivalent stress increases from the southwest to the northeast along these two cross sections, suggesting that strong earthquakes would move northeastward. In the cross section 13 passing through the earthquake 8, the
equivalent stress values around the source area become higher from shallow to deep depth, and are relatively smaller above 20 km depth in the southwest and northeast segments, which may explain why many small earthquakes occurred there. In the cross section 23 passing through the simulated earthquake 9, the
Y. Lu et al. / Journal of Asian Earth Sciences 50 (2012) 116–127
source areas extend northeastward and are in a high stress state from the shallow to deep depth, which may be a possible stress state after the earthquake 8.
127
instructive comments and suggestions, which improved our manuscript. References
6. Discussion and conclusions In this study a three-dimensional Maxwell viscoelastic model was established for North China using previous geological, geophysical and geodetic observations. In the model, the faults were set to be the weakened material. The initial stress field of North China was calculated using the finite element sub-model technique and considering the gravity and the boundary loadings from the Chinese continent model. Based on this initial model, we considered a variety of displacement boundary loading configurations to obtain the spatial distributions of 10 simulated earthquakes from the spatially distributed stress field and G values. After examining and comparing the displacements from our present initial model and the GPS observations, it is inferred that the initial model was considered to be reasonable. The simulated earthquakes are close to the actual ones and they are about 100 km apart from each other. The number of such earthquakes can amount to over 50%. Therefore, we conclude that the possible region of follow-up strong earthquakes could be predicted by obtaining the location of high stress disturbance values. Our results also show that it is critical to choose the approximate boundary loadings around the Chinese continent so as to obtain more reasonable parameters for North China, suggesting that the strong earthquake prediction in North China could depend in certain extent on the stress field and strong earthquake activities in the Chinese continent. The time sequence of simulated earthquakes can hardly correspond to that of actual ones, maybe due to the unknown trending and friction coefficient of the faults used to calculate the G value. Therefore, we infer that the main active region of strong earthquakes could be predicted using the finite element numerical simulation, but it is still hard to conduct an exact dynamic prediction due to the lack of a better understanding of the medium parameters, including the fault property and the crust and upper mantle structure. Acknowledgments We thank Academician Q. Deng, Profs. Z. Jiang and Z. Huang for providing the active faults, GPS velocity field and seismic velocity model in North China, respectively. This work was partially supported by the Exploration Technology Deep and Experimental Study (SinoProbe-06-03) and Special grants of Basic Scientific and Research (ZDJ2007-1 and ZDJ2009-1) and National Natural Science Foundations (40774044 and 40974021). We thank Prof. Bor-ming Jahn and two anonymous reviewers for providing many
Bai, W., Lin, B., Chen, Z., 2003. Numerical simulations on the block movement and deformation in North China after the 1976 Tangshan earthquake. Science in China (Series D) 33 (Suppl.), 99–107. Chen, L., Lu, Y., Liu, J., et al., 2001. Three-dimensional viscoelastic numerical simulation on dynamic evolution of stress field in North China induced by the 1966 Xingtai earthquake. Acta Seismologica Sinica 23 (5), 480–491. Chen, X., Zhang, J., Tang, R., et al., 2002. The Distribution of Moho Discontinuity in China and Surrounding Areas. Earthquake Press, Beijing, p. 1. Deng, Q., Zhang, P., Ran, Y., et al., 2002. Basic characteristics on the active tectonics in the Chinese continent. Science in China (Series D) 12, 1020–1030. Ding, J., Liu, J., Yu, S., 2000. Exploration and practice of earthquake prediction in China. Earthquake 20 (Suppl.), 13–17. Han, Z., Xu, J., Ran, Y., et al., 2003. Active blocks and strong earthquake activities in North China. Science in China (Series D) 33 (Suppl.), 108–118. Huang, Z., Li, H., Zheng, Y., Peng, Y., 2009. The lithosphere of North China Craton from surface wave tomography. Earth and Planetary Science Letters 288, 164– 173. Jiang, Z., Zhang, X., Chen, B., Xue, F., 2000. Characteristics of recent horizontal movement and strain-stress field in the crust of North China. Chinese Journal of Geophysics 43 (5), 657–665. Kenner, S.J., Segall, P., 2000. A mechanical model for intraplate earthquakes, application to the New Madrid seismic zone. Science 289, 2329–2332. Lei, J., Xie, F., Lan, C., Xing, C., Ma, S., 2008. Seismic images under the Beijing region inferred from P and PmP data. Physics of the Earth and Planetary Interiors 168, 134–146. Lei, J., Zhao, D., Xie, F., Liu, J., 2011. An attempt to detect temporal variations of crustal structure in the source area of the 2006 Wen—an earthquake in North China. Journal of Asian Earth Sciences 40, 958–976. Li, Q., Liu, M., Stein, S., 2009. Spatial–temporal complexity of continental intraplate seismicity: insights from geodynamic modeling and implications for seismic hazard estimation. Bulletin of the Seismological Society of America 99, 52–60. Lu, Y., 2004. Suggestions on promoting earthquake prediction in China. Recent Development in World Seismology 5, 81–84. Lu, Y., Yang, S., Chen, L., Lei, J., 2011. Mechanism of spatial distribution and migration of strong earthquakes in China inferred from numerical simulation. Journal of Asian Earth Sciences 40, 990–1001. Ma, H., Zhang, G., Liu, J., et al., 2003. Correlation between earthquake activity and activity crustal-block in China mainland and its adjacent regions. Earthquake Science Frontiers 10 (Suppl.), 75–80. Parsons, T., 2002. Post-1906 stress recovery of the San Andreas fault system calculated from three-dimensional finite element analysis. Journal of Geophysical Research 107, 1–13. Seyferth, M., Henk, A., 2002. Coupling of PFC2D and ANSYS – concepts to combine the best of two worlds for improved geodynamic models. In: Konietzky, H. (Ed.), Numerical Modeling in Micromechanics Via Particle Methods. Balkema, pp. 283–290. Tao, W., Hong, H., Liu, P., 2000. The numerical simulation of the strong earthquake region around the Chinese continent and surrounding regions. Acta Seismologica Sinica 22 (3), 271–277. Wang, R., Sun, X., Cai, Y., 1982. Numerical simulation of earthquake sequence over the past 700 years. Science in China (Series B) 8, 745–753. Zhang, Z., Deng, Y., Wang, B., 1990. Joint inversion of fracture system models of Ms P 7.0 earthquake sequences and rupture processes in North China. Acta Seismologica Sinica 12 (4), 335–347. Zhang, R., Wei, F., Qiao, C., et al., 2005. Numerical simulations for the 1975 Haicheng and 1999 Xiuyan earthquake processes by using the DDA + FEM method. Acta Seismologic Sinica 27 (3), 163–170. Zhu, Y., Wang, Q., Zhao, X., et al., 2004. Simulation of welding residual stress based on ANSYS. Journal of Wuhan University 26 (2), 69–72.