http://orcid.org/0000-0002-8395-2247
Figures
Abstract
The interaction of active faults as a factor affecting the mechanisms of large earthquakes has been observed in many places. Most aftershock and clustering earthquake sequences do not recur on the main seismogenic fault but are controlled by fault interactions with adjacent seismic structures. Four groups of conceptual models were generated in this study to determine how the geometry of the seismogenic faults controls the distributions of stress fields and earthquakes. The influences of the fault length ratio, center distance, overlap ratio, echelon distance and fault opening angle were considered in a 2D viscoelastic model. The results indicate that the interaction in the slipping zone is larger when collinear interacting faults are more closely positioned, with one fault lengthening. For noncollinear faults, the interaction is stronger as the inner tips pass each other, which impedes their growth after some degree of overlap. Additionally, fault interaction at the slipping zone becomes stronger as the opening angle approaches 180°. We further generated a 3D viscoelastic model of fault interactions in Central North China Block and applied the finite element method to analyze the relationship between distributions of earthquakes and fault geometry. The calculated results reveal well-matched higher stress and maximum shear strain concentrations in the southern part of the Fen-wei Graben Zone than in other zones in Central North China Block, which can be explained by the longer faults, shorter center distances, shorter overlap lengths and larger opening angles. The stress distributions and fault interactions should be considered in long-term seismic hazard assessment in these zones.
Citation: Shao B, Hou G (2019) The interactions of fault patterns and stress fields during active faulting in Central North China Block: Insights from numerical simulations. PLoS ONE 14(4): e0215893. https://doi.org/10.1371/journal.pone.0215893
Editor: Shuan-Hong Zhang, Chinese Academy of Geological Sciences, CHINA
Received: January 19, 2019; Accepted: April 10, 2019; Published: April 25, 2019
Copyright: © 2019 Shao, Hou. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: The work and the two authors were supported by funds from the National Natural Science Foundation of China (Grant No. 41530207) and State Key Projects (Grant No. 2016ZX05051004) to GH (http://www.nsfc.gov.cn/english/site_1/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Interactions between active faults have been observed in many places, and stress-change calculations for such interactions can reveal information about the dynamics and evolution of earthquakes, e.g., in California[1–4]. The motivation behind this study is to shed light on the dynamics of seismic initiation and migration by analyzing the interaction between different faults with various geometries and kinematics[4].
Problems involving fault interactions have received much attention through various analytical methods and physical and digital simulations over the past decades. Early studies focused on the stress regime and propagation mechanism of single faults[5–8]. Interactions between double fractures were then investigated in analytical solutions[9]. Only two collinear fractures were considered in these studies[10]. However, analytical solutions are too simple for real geologic conditions and are not suitable for noncollinear fractures.
The calculation of the required quantities in a viscoelastic structure is a sophisticated problem that necessitates advanced numerical methods[11]. Simulation of fault interactions has also been extensively employed in the previous literature. Physical simulation has been used to study the propagation of single-fracture and double-fracture models[12]. Additionally, the interaction between double strike-slip faults has been studied in physical simulation experiments[13]. Limited in terms of materials and processes, physical simulation modeling is often relatively simple but difficult to fit with the observed geology.
The finite element method (FEM) has been exclusively employed in fault interactions[1]. FEM analysis has also been performed under smooth contact conditions in both two and three dimensions. The influences of fault length, center distance, overlap ratio, en echelon distance and other factors have been considered in the research. However, all these solutions are incomplete and define the displacement and stress fields of only limited fault complexes around the cracks. Predecessors mainly considered one or two factors controlling fault interactions instead of coupling all factors and did not consider the interactions between double faults with different orientations.
In this contribution, multiple coupling factors that influence fault interaction are considered with a viscoelastic 2D FEM application. We further analyze the relationship between the distributions of earthquakes and faults in the Central North China Block based on fault interactions.
Conceptual modeling
Construction of the conceptual models
In this paper, we utilize four groups of 2D models with the same set of mechanical properties but different fault combinations to determine the roles of fault interactions using the FEM. The models have been developed within the viscoelastic fracture mechanics framework with 2D assumptions. Each time step is set to 100 years, and the entire modeling process consists of 100 steps. The ANSYS 14.0 (University Version 15.0; www.ansys.com) finite element software package is used to calculate the models.
The viscoelastic rheology is controlled by the Maxwell constitutive equation[14]:
(1)
where
is the strain rate, σ is the differential stress,
is the stress rate, E is Young’s modulus reflecting the elastic component, and ηeff is the effective viscosity. In the above viscoelastic model, Young’s modulus E and Poisson’s ratio ν can be calculated using the following equations, which essentially define Hooke’s law in three dimensions[15]:
(2)
(3)
(4)
(5)
where Vp represents the P-wave velocity, Vs represents the S-wave velocity, ρ represents rock density, μ and λ are Lamé parameters, E is Young’s modulus, and ν is Poisson’s ratio.
Most moderately strong earthquakes in North China occur in the gradient belt between different velocity blocks around the middle crust. The average rock mechanical properties of the 2D models could be regarded as those of the middle crust with a density of 2550 kg/m3, Young’s modulus of 30 GPa, Poisson’s ratio of 0.25 and viscosity of 1×1022 Pa•s[16, 17]. The coefficient of friction of the faults is taken to be 0.4 throughout our calculations, as the selection of the value was proven to have little effect on the spatial pattern of the Coulomb stress change[18].
The specific geometries considered in this paper are shown in Fig 1. Each group of models considers one factor that is related to fault interaction and control variables to keep the results credible. Model Group A includes five models of collinear faults interacting, which differ in the center distance. The five models of Model Group B are noncollinear parallel faults, which differ in the overlap ratio. Model Group C includes five models of interacting collinear faults, which differ in the length of the upper faults. Model Group D investigates the influence of the fracture opening angle using round models to eliminate errors caused by different distances between the faults and the boundaries.
Model A includes models differing in the center distance. Model B includes models differing in the overlap ratio. Model C includes models differing in the en echelon distance. Model D includes models differing in the length of the upper fault. Model E includes models differing in the fracture opening angle.
Constraints on the conceptual models
In finite element methods, the application of boundary conditions remains one of the key problems[19–23]. In the modeling, four groups of two-dimensional viscoelastic models of fault interaction are employed. In each group, five models are employed with different parameters of faults (Table 1). In our first four groups of conceptual numerical models (Models A, B, C and D), boundary conditions are applied as follows: the up and down boundaries of the model are set free, i.e., all stress components are nil. The left boundary is vertically and horizontally fixed to velocity zero. A 50 MPa compressive deviatoric stress is set on the right boundary. In the last group of conceptual numerical models on the influence of fracture opening angle (Model E), a 50 MPa compressive deviatoric stress is set on the circumference to simulate the extensional or compressional environment, in order to reduce the difference in results caused by the difference in load direction. To verify the stress field under different stress conditions, we also set a 50 MPa tensile deviatoric stress on the right boundary in a group of comparative trials.
Modeling results
Center distance.
The results for the differential stress of two equal-length collinear faults with different center distances are shown in Fig 2-Model A. The center distance is defined as the distance between the inner tips of the double faults. Fig 3(a) illustrates the variations in the calculated differential stress at points A, B and C as a function of normalized center distance for five different separations of the two adjacent faults[24]. For normalization, the center distance is divided by the length L of a single fault of the same length. As shown in Fig 3, in general, the interaction is larger as the interacting faults are more closely positioned. The differential stresses at points A and C represent the influence of the double faults on the far ends. The differential stress at point B represents the influence of the interaction on the zone between the double faults. The curve illustrates that the interaction is significant only when the slipping zones are closer than half of the fault length, which shows good agreement with the experimental results of Bobet and Einstein [12].
Black circles denote points A, B and C in each model, where the differential stresses are shown in the curves in Fig 3.
In (a), the variations in the calculated differential stresses at points A, B and C are illustrated as a function of the normalized center distance. In (b), the variations in the calculated differential stresses at points A, B and C are illustrated as a function of the normalized overlap length. In (c), the variations in the calculated differential stresses at points A, B and C are illustrated as a function of normalized en echelon distance. In (d), the variations in the calculated differential stresses at points A, B and C are illustrated as a function of the normalized length of the shorter fault. In (e), the variations in the calculated differential stresses at points A, B and C are illustrated as a function of the opening angles.
Overlap length.
The results for the differential stress of two equal-length collinear faults with different overlap ratios are shown in Fig 2-Model B. Overlap length is defined as the length of the overlapping portion of the double faults. Fig 3(b) illustrates the variations in the calculated differential stresses at points A, B and C as a function of normalized overlap length for five different separations of the two adjacent faults. For normalization, the overlap length is divided by L of a single fault of the same length[24].
As shown in Fig 3, the differential stresses at points A and C represent the influence of the double faults on the far ends, which remains invariant as the overlap ratio changes. The differential stress at point B represents the influence of the interaction on the zone between the double faults. The results from the curve indicate that fault interaction is stronger as the inner tips pass each other and impedes their growth after some degree of overlap. The calculated differential stress increases sharply as the inner tips approach each other and decreases after the inner tips pass each other to form a larger overlap ratio, which shows good agreement with studies on the San Andreas Fault[1].
Echelon distance.
The results for the differential stress of two equal-length collinear faults with different echelon distances are shown in Fig 2-Model C. Echelon distance is defined as the orthogonal distance between two echelon faults with some degrees of overlap. Fig 3(c) illustrates the variations in the calculated differential stresses at points A, B and C as a function of normalized echelon distance for five different separations of the two adjacent faults. For normalization, the echelon distance is divided by L of a single fault of the same length[24].
As shown in Fig 3, the differential stresses at points A and C represent the influence of the double faults on the far ends, which remains invariant as the overlap ratio changes. The differential stress at point B represents the influence of the interaction on the zone between the double faults. The results from the curve indicate that fault interaction is stronger as the echelon distance becomes shorter. Only when the echelon distance is approximately one-fifth of the fault length is the calculated differential stress in the overlap area larger than the stress on the far ends, which shows good agreement with studies on the San Andreas Fault[1].
Fault lengths.
The results for the differential stress of two collinear faults with different fault lengths are shown in Fig 2-Model D. Fig 3(d) illustrates the variations in the calculated differential stresses at points A, B and C as a function of the normalized fault length[24]. For normalization, the upper fault length is divided by L of the invariable lower fault length. As shown in Fig 3, the differential stress at point B represents the influence of the interaction on the zone between the double faults. The results from the curve indicate that fault interaction is stronger as one fault becomes longer. The result emerges that a large preexisting slip zone can reactivate slip and movement along a nearby smaller slip zone. In other words, segmentation of fault slip zones can substantially affect the interpretation of observations, which shows good agreement with experiments[9].
Opening angle.
The results for the differential stress of two faults with different fault opening angles are shown in Fig 2-Model E. The opening angle is defined as the minimum angle between the two faults. Fig 3(e) illustrates the variations in the calculated differential stresses at points A, B and C as a function of the opening angle[24]. The differential stresses at points A and C represent the influence of the double faults on the far ends, which remains invariant as the opening angle changes. The differential stress at point B represents the influence of the interaction on the zone between the double faults. For normalization, Model Group D uses round models to eliminate errors caused by different distances between the faults and the boundaries. The results from the curve indicate that fault interaction at the slipping zone is stronger as the opening angle is closer to 180°. This result implies that collinear faults are more likely to coalesce than faults with an opening angle.
In short, the results from the curves indicate that fault interactions are stronger with longer fault length, shorter overlap length, shorter center distance and larger opening angle. The coupling of the four factors should be taken into account in the stress distribution and long-term seismic hazard assessment of seismic zones and gaps.
Fault interaction modeling in Central North China Block
Seismotectonic setting
The Central North China Block (CNCB) is an old tectonic unit in mainland China, which is also seismically one of the most active intracontinental regions in the world[25–27]. Many previous studies have shown that North China has experienced significant lithospheric thinning since the late Mesozoic[28–31]. The CNCB underwent intensive compression-shear faulting during the earliest Cenozoic[32, 33]. Many NNE- and NE-trending structures passing through the interior of the CNCB developed or were reactivated then[34, 35]. During the early Paleocene–Eocene interval, the NNE- to NE-trending faults were inverted to normal faults, and numerous new normal faults developed. The main seismic belts are oriented in a NNE direction, including the Tangshan–Cixian seismic belt, the Shanxi seismic belt, and the Tan-lu seismic belt. However, the Zhangjiakou-Bohai seismic zone is oriented in a ENE direction.
Setup of fault interaction model in the CNCB
To reveal how fault interactions controlled stress evolution and earthquakes in the CNCB, a 3D finite element model[36, 37] was applied.
The model was developed within a viscoelastic fracture mechanics framework. Each time step was set to 100 years, and the entire modeling process underwent 100 steps. The model covered a region from 32° to 42°N and 106° to 124°E. The depth of the model was approximately 100 km in accordance with the estimated average lithospheric thickness in the CNCB. The lithosphere was modeled using a layered rheology structure, which included five types of rocks (sediment, upper crust, middle crust, lower crust, and mantle; Table 2). During the calculation process, the selection of the parameters was based on the crustal velocity and viscous structure from previous studies[38] (Table 1).
The 3D model includes the main seismogenic faults associated with devastating earthquakes (M≥6.5) in Central North China Block since A.D. 1303 [38] (Fig 4), since the accuracy of the seismic parameters before 1303 is poor.
The black lines depict the distribution of major seismogenic faults in the model (after Shen et al., 2004). Red circles mark the locations of earthquakes (M≥6.5) since 1303; white arrows depict displacement directions of the model boundaries. The surface of the model is set free, the bottom horizontally and vertically fixed to velocity zero.
Earthquake location, magnitude and intensity data are from the Chinese Earthquake Catalog[39, 40]. The length of the seismogenic faults is mainly controlled by the empirical intensity-fault length equation, which is derived from the statistics for historical earthquakes in North China[38] and this equation is:
(6)
For historical earthquakes with magnitude only and no seismic intensity distribution data, the empirical magnitude-aftershock area length equation of North China is used:
(7)
where L is the length of the seismogenic fault, T is the length of the seismic intensity zone VIII and M is the Richter-scale magnitude. The main fault details are recorded in Table 3.
In the 3D model, faults were set as discrete planes of weakness cutting the FE model. These planes were described by so-called contact elements, which were defined at opposite sides of preassigned faults[41, 42]. The faults were set to propagate to the bottom of the upper crust without considering the fault occurrence in the model. The coefficient of friction of the faults was taken to be 0.4 throughout our calculations.
Boundary conditions
Boundary conditions were applied based on GPS data and the inverse method[23, 43], as follows: the surface of the model was free; the bottom of the model was horizontally and vertically fixed to velocity zero; SE-directed displacements were set on the left and right boundaries; and ESE-directed displacements were set on the north and south boundaries. The magnitude of the displacements varied according to the GPS data[44] (Fig 4).
Modeling results
Through the 3D viscoelastic model, investigations were made into the distribution of the tectonic stress field controlled by active fault interactions in the CNCB[23]. The 3D model is to reveal how the laws from the 2d fault interaction models controlled stress evolution and earthquakes in the CNCB. Also, the results from the 3D model can be a verification of the 2d fault interaction models.
Figs 5 and 6 show the maximum and minimum principal strain calculated with the 3D model in the CNCB. On the maximum principal strain figure, an approximately N-S-trending extension is located in the southeast part of the research area, while an approximately N-S-trending shortening is located in the northwest part (Fig 5). On the minimum principal strain figure, an approximately E-W-trending extension is located in the east-central part of the research area, while an approximately N-S-trending shortening is located in the west part (Fig 6). The modeling results are consistent with the phenomenon that the North China Plate escapes to the southeast and with the topography of the eastern basin and western plateau, which is due to the Indian Plate subducting beneath the Eurasian plate and blocking from the Siberian Plate.
Black lines denote the seismogenic faults. An approximately N-S-trending extension is located in the southeast part of the research area, while an approximately N-S-trending shortening is located in the northwest part.
Black lines denote the seismogenic faults. An approximately E-W-trending extension is located in the east-central part of the research area, while an approximately N-S-trending shortening is located in the west part.
The maximum shear strain and differential stress intensity triggered by stick slip and steady-state slip along the seismogenic faults are presented in Figs 7 and 8. The differential stress and maximum shear strain at a depth of 10 km are calculated because most large earthquakes occur at this depth in the CNCB. Higher differential stress and maximum shear strain could result in new active faults and more earthquakes, which indicates the range and tendency of stress changes controlled by fault interactions[17].
Black lines denote the seismogenic faults. The highest strain value is located in the southern part of the Fen-wei Graben Zone. The lowest strain value is located in the northern part of the Fen-wei Graben Zone. Moderate strain values are located in the Bohai-Zhangjiakou Fault Zone and the Tan-lu Fault Zone.
Black lines denote the seismogenic faults. The highest stress value is located in the southern part of the Fen-wei Graben Zone. The lowest stress value is located in the northern part of the Fen-wei Graben Zone. Moderate stress values are located in the Bohai-Zhangjiakou Fault Zone and the Tan-lu Fault Zone.
The laws from the model results of the conceptual model can be verified in the 3D model. For instance, the “fault length” model can be used to explain the higher differential stress and maximum shear strain values in the southern part of the Fen-wei Graben Zone south of Taiyuan. The zone with the highest differential stress and maximum shear strain values in the CNCB can be explained by the longer faults, shorter center distances, shorter overlap lengths and larger opening angles of F1 and F2 (Fig 8). The moderate value between F3 and F4 in the Tan-lu Fault Zone can be explained by the “fault length” model as well. Although F3 is the longest fault in the CNCB and the opening angle between F3 and F4 is nearly 180°, the short length of F4 is the main reason for the moderate stress value in the zone between F3 and F4. The echelon distance model can be used to explain the relatively high stress values in the Tangshan Fault Zone, in which moderate faults with short echelon distances, such as F5 and F6, are distributed. The “center distance” and the “fault length” model can be used to explain the low stress values in the northern part of the Fen-wei Graben Zone north of Taiyuan (Fig 8), in which short faults with relatively long center distances and longer echelon distances, such as F8, F9 and F10, are distributed. The “overlap length” model can be used to explain the relatively moderate stress values in the slipping zones between echelon faults in the Bohai-Zhangjiakou Fault Zone, such as the zone between F11 and F12, which are relatively long faults with a large overlap length and short normalized center distance (Fig 8).
Discussion
Factors affecting the mechanisms of large earthquakes in the CNCB mainly include the coupling of the velocity anomaly body and active faults[45–47], the constraints on the boundaries of different active blocks[48, 49] and the anisotropy of the lithosphere[50]. The current seismic hazard assessment relies on the characteristic earthquake recurrence model, i.e., the detailed locations of the active faults[51, 52]. However, most of the aftershock and clustering earthquake sequences do not recur on the main seismogenic fault but are controlled by the fault interactions of the adjacent seismic structures [53]. Additionally, the mechanism of newly formed active faults is the result of interactions, such as the formation of the San Andreas Fault at the Mendocino triple junction [54].
Problems involving fault interactions have received much attention through various analytical methods and physical and digital simulations over the past decades. However, uncertainty exists in the way that the geometry of the seismogenic faults controls the distribution of stress fields and earthquakes. We generated four groups of conceptual models to reveal the influences of the fault length ratio, center distance, overlap coincidence ratio and opening angle in our study.
In the 3D viscoelastic model of the CNCB, seismogenic faults are used instead of shallow active faults to study fault interactions. The distribution of shallow active faults, which has been used in some models of the CNCB[55, 56], is the shallow manifestation of deep seismogenic structures. Unlike the active faults, the seismogenic faults have magnitudes that are more proportional to the intensity of the earthquake and the degree of influence on the stress field.
The factors of fault interaction should be considered in long-term seismic hazard assessment in these zones. The calculated results reveal a well-matched highest stress concentration between Xi’an and Taiyuan, which can be explained by the longer faults, shorter center distances, shorter echelon distances, shorter overlap lengths and larger opening angles of the faults in the zone. Between Beijing and Tianjin, the shorter fault lengths, longer center distances, longer overlap lengths and smaller opening angles lead to lower stress concentrations. The further work need to do is quantificationally analysis fault interaction, such as works from Liu and Konietzky [57, 58].
Conclusions
The conceptual models considered the influences of the fault length ratio, center distance, overlap coincidence ratio and opening angle, and the following primary conclusions are drawn:
- For two equal-length collinear faults, the interaction is larger as the interacting faults are more closely positioned. Additionally, the modeling makes it clear that the interaction is significant only when the slipping zones are closer than half of the fault length.
- The interaction in the slipping zone increases sharply as the inner tips approach each other and decreases after the inner tips pass each other to form a larger overlap ratio. Additionally, compared with Model Group A, fault interaction is stronger as the inner tips pass each other, which impedes their growth after some degree of overlap. Moreover, fault interaction is stronger as the echelon distance becomes shorter. Only when the echelon distance is approximately one-fifth of the fault length is the calculated differential stress in the overlap area larger than that on the far ends.
- For two unequal-length collinear faults, the interaction is stronger as one fault becomes longer. The result emerges that a large preexisting slip zone can reactivate slip and movement along a nearby smaller slip zone. In other words, segmentation of fault slip zones can substantially affect the interpretation of observations[9].
- The fault interaction at the slipping zone becomes stronger as the opening angle approaches 180°. This result indicates that collinear faults are more likely to coalesce than faults with an opening angle. In other words, when the interval is constant, the risk of a seismic gap is greater when the faults on both sides are collinear. This factor should be taken into account in long-term seismic hazard assessment of a seismic block zone, background gap and preparatory gap.
Based on geological evidence and previous research in North China, the stress field is calculated by the finite element method for a 3D viscoelastic model to reveal the mechanism of seismic fault interactions. The interactions of seismic faults in Central North China Block can be explained by conceptual models. The calculated results reveal well-matched highest stress and maximum shear strain concentrations for Central North China Block located in the southern part of the Fen-wei Graben Zone, which can be explained by the longer faults, shorter center distances, shorter overlap lengths and larger opening angles of the faults in this area. The lowest stress concentration for Central North China Block located in the northern part of the Fen-wei Graben Zone can be explained by the short fault lengths and long center distances. In the Bohai-Zhangjiakou Fault Zone, longer center distances, longer overlap lengths and smaller opening angles lead to moderate stress concentrations. The relatively moderate stress values in the Tan-lu Fault Zone can be explained by interactions between short and long faults.
The factors of fault interactions and calculated stresses under the present-day stress field distribution should be considered in long-term seismic hazard assessment in North China. The southern part of the Fen-wei Graben Zone may have a relatively high seismic risk in the future.
Acknowledgments
We express our sincere gratitude to Dr. Jie Li and Ph.D. student Zhifang Zhou, who provided linguistic assistance. We also thank our editor for the substantial help.
References
- 1. Aydin A, Schultz RA. Effect of mechanical interaction on the development of strike-slip faults with echelon patterns. Journal of Structural Geology. 1990; 12(1): 123–129.
- 2. Jónsson S. Stress interaction between magma accumulation and trapdoor faulting on Sierra Negra volcano, Galápagos. Tectonophysics. 2009; 471(1): 36–44.
- 3. Sorlienet CC, Bennett JT, Cormier MH, Campbell BA, Nicholson C, Bauer RL, et al. Late Miocene–Quaternary fault evolution and interaction in the southern California Inner Continental Borderland. Geosphere. 2015; 11(4): 1111–1132.
- 4. Delorey AA, van der Elst NJ, Johnson PA. Tidal triggering of earthquakes suggests poroelastic behavior on the San Andreas Fault. Earth and Planetary Science Letters. 2017; 460: 164–17.
- 5. Williams ML. The stresses around a fault or crack in dissimilar media. Bulletin of the seismological society of America. 1959; 49(2): 199–204.
- 6. Chinnery MA. The stress changes that accompany strike-slip faulting. Bulletin of the Seismological Society of America. 1963; 53(5): 921–932.
- 7. Tan CL, Fenner RT. Elastic fracture mechanics analysis by the boundary integral equation method. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 1979; 369: 243–260.
- 8. Koller MG, Bonnet M, Madariaga R. Modelling of dynamical crack propagation using time-domain boundary integral equations. Wave Motion. 1992; 16(4): 339–366.
- 9. Rudnicki JW, Kanamori H. Effects of fault interaction on moment, stress drop, and strain energy release. Journal of Geophysical Research: Solid Earth. 1981; 86(B3): 1785–1793.
- 10. Golecki JJ, Ornai D. Extension fracture due to collinear cracks. Engineering fracture mechanics. 1996; 54(1): 25–33.
- 11. Brocks W, Müller W, Olschewski J. Experiences in applying ADINA to the analysis of crack tip fields in elastic-plastic fracture mechanics. Computers Structures. 1985; 21(1–2): 137–158.
- 12. Bobet A, Einstein HH. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. International Journal of Rock Mechanics and Mining Sciences. 1998; 35(7): 863–888.
- 13. Ma J, Liu L, Liu P, Ma S. Thermal Precursory Pattern of Fault Unstable Slip: An Experimental Study of En Echelon Faults. Chinese Journal of Geophysics. 2007; 50(4): 995–1004.
- 14. Christensen RM, Freund LB. Theory of viscoelasticity. Journal of Applied Mechanics. 1971; 38: 720.
- 15. Zhan Y, Hou G, Kusky T, Gregg PM. Stress development in heterogenetic lithosphere: Insights into earthquake processes in the New Madrid Seismic Zone. Tectonophysics. 2016. 671; 56–62.
- 16. De Franco R, Govers R, Wortel R. Numerical comparison of different convergent plate contacts: subduction channel and subduction fault. Geophysical Journal International. 2007; 171(1): 435–450.
- 17. Lu Y, Yang S, Chen L, Lei J, He P. Migration trend of strong earthquakes in North China from numerical simulations. Journal of Asian Earth Sciences. 2012; 50: 116–127.
- 18. Deng J, Sykes LR. Evolution of the stress field in southern california and triggering of moderate-size earthquakes: a 200-year perspective. Journal of Geophysical Research Solid Earth. 1997; 102(B5): 9859–9886.
- 19. Liu J, Shi Y, Zhang G. Seismicity model based on the pseudo-3d finite element method. Chinese Journal of Geophysics. 2001; 44(6): 800–810.
- 20. Zhao Z, Xu J. Basin structure and numerical simulation for the mechanisms of seismic disasters. Acta Geologica Sinica (English edition). 2004; 78(4): 954–959.
- 21. Xu W, Hou Z, Yan Z, Shi D, Meng Y, Zeng P. Numerical Simulation of Fluid Migration during Ore Formation of Carboniferous Exhalation-Sedimentary Massive Sulfide Deposits in the Tongling District, Anhui Province. Acta Geologica Sinica (English edition). 2005; 79(1): 98–105.
- 22. Hou G, Wang C, Li J, Qian X. Late Paleoproterozoic extension and a paleostress field reconstruction of the North China Craton. Tectonophysics. 2006; 422(1): 89–98.
- 23. Dai L, Li Q, Li S, Guo L, Somerville ID, Xu L, et al. Numerical modelling of stress fields and earthquakes jointly controlled by NE-and NW-trending fault zones in the Central North China Block. Journal of Asian Earth Sciences. 2015; 114: 28–40.
- 24. Kotousov A, Berto F, Lazzarin P, Pegorin F. Three dimensional finite element mixed fracture mode under anti-plane loading of a crack. Theoretical and Applied Fracture Mechanics. 2012; 62: 26–33.
- 25. Huang Q, Ma Y, Song X. Seismic Activity and Fuzzy Comprehensive Evaluation of Crustal Stability in the Beijing Region. Acta Geologica Sinica (English Edition). 1994; 7(1): 1–13.
- 26. Wang Y, Cao J, Zhu J. Seismic velocity structure and composition of the continental crust of eastern China. Acta Geologica Sinica (English Edition). 2004; 78(1): 291–297.
- 27.
Huang J, Zhao D. Lithosphere Detachment beneath North China Craton imaged by joint inversion of regional and teleseismic data. In AGU Fall Meeting Abstracts; 2009.
- 28. Xu YG. Thermo-tectonic destruction of the Archaean lithospheric keel beneath the Sino-Korean Craton in China: Evidence, timing and mechanism. Physics and Chemistry of the Earth. Solid Earth and Geodesy. 2001; 26(9): 747–757.
- 29. Ren J, Tamaki K, Li S, Junxia Z. Late Mesozoic and Cenozoic rifting and its dynamic setting in Eastern China and adjacent areas. Tectonophysics. 2002; 344(3): 175–205.
- 30. Hou G, Qian X, Cai D. Cretaceous-Paleogene rifting in the Bohai Sea and Western Shandong: geochemical evidence from volcanic rocks. Chinese Journal of Geology. 2003; 38(1): 13–21.
- 31. Ma L, Jiang S, Albrecht WH, Xu Y. A Possible Mechanism to Thin Lithosphere of the North China Craton: Insights from Cretaceous Mafic Dikes in the Jiaodong Pensinsula. Acta Geologica Sinica (English Edition). 2016; 90(S1):106–108.
- 32. Hou G, Liu Y, Li J. Evidence for∼ 1.8 Ga extension of the Eastern Block of the North China Craton from SHRIMP U–Pb dating of mafic dyke swarms in Shandong Province. Journal of Asian Earth Sciences. 2006; 27(4): 392–401.
- 33. Xu L, Li S, Cao X, Somerville ID, Suo Y, Liu X, et al. Holocene intracontinental deformation of the northern North China Plain: evidence of tectonic ground fissures. Journal of Asian Earth Sciences. 2016; 119: 49–64.
- 34. Yin A, Nie S. An indentation model for the North and South China collision and the development of the Tan-Lu and Honam fault systems, eastern Asia. Tectonics. 1993; 12(4): 801–813.
- 35. Liu M, Yang Y, Shen Z, Wang S, Wang M, Wan Y. Active tectonics and intracontinental earthquakes in China: The kinematics and geodynamics. Geological Society of America Special Papers. 2007; 425: 299–318.
- 36.
Zienkiewicz OC, Taylor RL. The finite element method: solid mechanics (Vol. 2). Butterworth Heinemann; 2000.
- 37. Zhu S, Zhang P. FEM simulation of interseismic and coseismic deformation associated with the 2008 Wenchuan Earthquake. Tectonophysics. 2013; 584: 64–80.
- 38. Shen Z, Wan Y, Gan W, Li T, Zeng Y. Crustal stress evolution of the last 700 years in North China and earthquake occurrence. Earthquake Research of China. 2004; 20: 211–228.
- 39.
State Seismological Bureau. Historical Strong Earth-quake Catalogue of China (23rd Century B.C. to 1911A.D.). Seismological Press; 1995
- 40.
China Seismological Bureau. Strong Earthquake Catalogue of China (23rd Century B.C. to 1999 A.D.). Seismological Press, Beijing; 1999.
- 41. Smart KJ, Ferrill DA, Morris AP. Impact of interlayer slip on fracture prediction from geomechanical models of fault-related folds. Amer Assoc Petroleum Geologist Bulletin. 2009; 93(11): 1447–1458.
- 42. Fischer K, Henk A. 3-D geomechanical modeling of a gas reservoir in the North German Basin: Workflow for model building and calibration. Solid Earth Discussions. 2013; 5: 767–788.
- 43. Wang Q, Zhang PZ, Freymueller JT, Bilham R, Larson KM, Lai XA, et al. Present-day crustal deformation in China constrained by global positioning system measurements. Science. 2001; 294(5542): 574–577. pmid:11641493
- 44. Liu X, Ma J, Fu R, Yang G, Shao Z, Zheng Z. Primary study on the dynamics of the present-day crustal motion in North China Region. Chinese Journal of Geophysics. 2010; 53(6): 1418–1427.
- 45. Li S, Suo Y, Deng L. Development of the Bohai Bay Basin and destruction of the North China Craton. Earth Science Frontiers, 2010; 17(4): 64–89.
- 46. Li S, Zhao G, Dai L, Zhou L, Liu X, Suo Y, et al. Cenozoic faulting of the Bohai Bay Basin and its bearing on the destruction of the eastern North China Craton. Journal of Asian Earth Sciences. 2012; 47(1): 80–93.
- 47. Shi W, Cen M, Chen L, Wang Y, Chen X, Li J, et al. Evolution of the late cenozoic tectonic stress regime in the Shanxi rift, central North China plate inferred from new fault kinematic analysis. Journal of Asian Earth Sciences. 2015; 114: 54–72.
- 48. Chen L, Lu Y, Guo R, Xu G, Zhang J. Evolution of 3d tectonic stress field and fault movement in North China. Acta Seismologica Sinica. 2001; 14(4): 371–383.
- 49. Zhang P, Deng Q, Zhang G, Ma J, Gan W, Min W, et al. Active tectonic blocks and strong earthquakes in the continent of China. Science in China Series D: Earth Sciences. 2003; 46: 13–24.
- 50. Shen J, Shao B, Yu X, Yu Y, Qi G, Deng M, et al. Seismogenic tectonics of the Qian-gorlos earthquake in Jilin Province, China. Earthquake Science. 2016; 29(2): 93–103.
- 51.
Petersen MD, Bryant WA, Cramer CH, Cao T, Reichle MS, Frankel AD, et al. Probabilistic seismic hazard assessment for the state of California No. 96–706. California Dept. of Conservation Division of Mines and Geology; 1996.
- 52. Meletti C, Galadini F, Valensise G, Stucchi M, Basili R, Barba S, et al. A seismic source zone model for the seismic hazard assessment of the Italian territory. Tectonophysics. 2008; 450(1): 85–108.
- 53. Dieterich J. A constitutive law for rate of earthquake production and its application to earthquake clustering. Journal of Geophysical Research: Solid Earth. 1994; 99(B2): 2601–2618.
- 54. Hole JA, Beaudoin BC, Klemperer SL. Vertical extent of the newborn San Andreas fault at the Mendocino triple junction. Geology. 2000; 28(12): 1111–1114.
- 55. Liu C, Shi Y, Zheng L, Zhu B. Relation between earthquake spatial distribution and tectonic stress accumulation in the North China Basin based on 3d visco-elastic modelling. Chinese Journal of Geophysics. 2012; 55(12): 3942–3957.
- 56. Liu L, Hao T, Lü C, You Q, Pan J, Wang F, et al. Crustal structure of Bohai Sea and adjacent area (North China) from two onshore–offshore wide-angle seismic survey lines. Journal of Asian Earth Sciences. 2015; 98: 457–469.
- 57. Liu Y, Konietzky H. Particle-Based Modeling of Pull-Apart Basin Development. Tectonics, 2018; 37: 343–358.
- 58. Liu Y, Konietzky H. Particle-based modeling of transtensional pull-apart basins. Tectonics, 2018; 37:4700–4713.