Geomechanical Modeling of Stress and Fracture Distribution during Contractional Fault-Related Folding

Jianwei Feng^1,2*, Kaikai Gu^1,2

¹School of Geosciences, China University of Petroleum, Qingdao, China.
²Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao, China.

DOI: 10.4236/gep.2017.511006                

Abstract

Understanding and predicting the distribution of fractures in the deep tight sandstone reservoir are important for both gas exploration and exploitation activities in Kuqa Depression. We analyzed the characteristics of regional structural evolution and paleotectonic stress setting based on acoustic emission tests and structural feature analysis. Several suites of geomechanical models and experiments were developed to analyze how the geological factors influenced and controlled the development and distribution of fractures during folding. The multilayer model used elasto-plastic finite element method to capture the stress variations and slip along bedding surfaces, and allowed large deformation. The simulated results demonstrate that this novel Quasi-Binary Method coupling composite failure criterion and geomechanical model can effectively quantitatively predict the developed area of fracture parameters in fault-related folds. High-density regions of fractures are mainly located in the fold limbs during initial folding stage, then gradually migrate from forelimb to backlimb, from limbs to hinge, from deep to shallow along with the fold uplift. Among these factors, the fold uplift and slip displacement along fault have the most important influence on distributions of fractures and stress field, meanwhile the lithology and distance to fault have also has certain influences. When the uplift height exceeds approximately 55 percent of the total height of fold the facture density reaches a peak, which conforms to typical top-graben fold type with large amplitude and high-density factures in the top. The overall simulated results match well with core observation and FMI results both in the whole geometry and fracture distribution.

Keywords

Kuqa Depression, Deep Tight Sandstone, Fault-Propagation Fold, Geomechanical Model, Fractures Prediction

1. Introduction

Natural fractures influence the performance of many reservoirs around the world, including carbonate reservoir, deep tight sand reservoir and low-per- meability reservoir in the world [1] [2] [3] [4] . Understanding and interpretation of where and when fractures would develop within a geological structure along with their orientation and intensity can be important to both exploration and production planning activities. The geological factors controlling development of the fractures include proximity to faults, position on folds, differential stress, lithology and their combination and layer thickness [5] - [11] . Many studies have found that fractures often develop around fault zones and anticlinal core, and that fracture spacing is positively correlated with regional differential stress [12] [13] [14] . Some researchers have suggested that lithology and combination of fractures may be more significant than regional stress [15] [16] [17] [18] [19] . Some have even found that the natural fracture networks tend to be heterogeneous and messy near formation interfaces, when the stratum consists of hard and soft rocks, in which fractures are clustered in swarms and are irregularly distributed [7] [20] .

As practiced in fractured reservoir exploration and production, fracture prediction is commonly based on geometric and/or kinematic models, such as analyses of fault-related folds and fold curvature [20] [21] [22] [23] [24] or seismic and logging techniques [25] , to acquire interwell fracture networks and the attributes, but far less on geomechanical modeling method. Many studies have found that fractures often develop around fault zones and anticlinal core, and that fracture spacing is positively correlated with regional stress intensity. Numerical geomechanical modeling such as finite element method (FEM), boundary element method (BEM) and discrete element method (DEM) can provide powerful tools for simulating the spatial and temporal development of geological structures [26] - [34] . In comparison, as an accurate and effective engineering numerical analysis method, the BEM can use simple element to simulate the boundary shape and get high precision. However, the BEM has been used to a limited extent, primarily in specialized applications, such as the heterogeneous formation of basins. Similarly, as a dynamic numerical analysis method, the DEM used to simulate the failure process of rocks, permitting large deformation and obvious slide [35] [36] [37] [38] . Without doubt, rock rheology of the DEM should be further considered, some weak points should be improved and perfected, such as treatment of free movable boundary and coupling problem. In current study, the researchers found the application of finite-element-based geomechanical models to have excellent potential for understanding and interpreting natural fractures in geologic structures. Finite element modeling allows complex geometries (e.g., faults and mechanical stratigraphy) to be combined with realistic material models to produce physically realistic and mechanically rigorous forward models. The basic premise of this approach is that the geomechanical-model-derived permanent strains can be acquired in terms of fracture characteristics (i.e. location, orientation, and intensity). This approach captures the kinematic history and further permits tracking the spatial and temporal evolution of stress and strain in the deformed rock or bedding [4] [24] [33] [39] .

Since the 1960s, there have been many studies on the mechanical methods of structural movement generating fractures, including rock failure criterion and strain energy density. Based on laboratory modeling experiments and tests, [5] [40] proposed that the fracture intensity had fine positive correlation with elastic strain energy in rocks. [41] and [42] successfully used tensional failure criterion and shear failure criterion to calculate rock cracking index in predicting fractures development zone and dominant orientation. [43] put forward a quasi-binary method to quantitatively characterize fracture density using based on strength criterion and energy criterion for which their calculated results were in good agreement with practical measurements. [44] and [45] tried to establish a series of formulas between stress-strain and fracture volume density based on such assumption that paleostress field generates cracking and the current stress field only induces some minor changes in fractures size instead of producing new cracks. Based on interpretation of the geological structure, rock mechanics test, and acoustic emission experiment, [46] simulated the three-dimensional paleotectonic stress field during the Yanshan movement, and used rock failure criterion and comprehensive evaluation coefficient of fractures to determine the quantitative development of fractures.

In this study, we use the finite-element-method (FEM) to better simulate when and where paleotectonic differential stress develop within a fault-propagation fold (D gas field) throughout the entire deformational history based on the analysis of the structural evolutional characteristics, rock mechanical test, and acoustic emission experiments. Then we quantitatively predict the development and distribution zones of tectonic fractures based on composite rock failure criterion and geomechanical model between fracture density and strain energy density. Meanwhile, the predicted results can be verified through the fracture density distributions identified from cores, boreholes and the capacity of gas wells. The ultimate goal is to develop several suites of geomechanical models and experiments for evaluating and analyzing how these important mechanical factors (e.g., stress field, slip displacement of fault, mechanical parameters difference among layers and uplift amplitude of fold) control and influence the development and distribution of these fractures in fault-propagation fold. According to the obtained results, the adopted methodology proved successful in predicting tectonic fractures of tight sandstone reservoirs. Such reality shows that the methods for predicting fractures can provide an important technological approach to exploration and development of gas field in the deep tight sandstone reservoirs.

2. Geologic Setting

The Kuqa Depression is located at the northern margin of the Tarim Basin between the South Tianshan Orogenic Belt and the Northern Tarim Uplift to the south (Figure 1). Thrust faults and related folds widely developed in the Kuqa Depression during the Cenozoic time structures. Laterallythe Kuqa Depression can be divided intothree structural belts and two sags, which are the northern monocline belt, Kelasu structural belt, Baicheng sag, Kuqa depression, and Qi- ulitage structural belt from north to south as depicted in Figure 1 [18] [47] . Since the late Cretaceous, the Kuqa Depression has experienced a complex evolutionary history as a consequence of the northward Indian sub-continent and southward thrusting of the South Tianshan, and is recognized as one of the major Cenozoic depocenters along the margin of the Tarim Basin [48] [49] [50] . Within the Kuqa depression, D gas field is situated in the hanging wall of Qiulitage tectonic belt and north of Baicheng fault with an exploration area of approximately 245 km², displaying a fault-propagation anticline with structure amplitude more than 700 m, the top of which was cross-cut by large-scale Dinan fault and subsidiary Dibei fault striking EW to make it more complicated. The dip of southern limb ranges from 16˚ to 30˚ and the dip of northern limb, ranges from 19˚ to 33˚ (Figure 1).

The Qiulitage tectonic belt and D anticline was mainly characterized by three nearly SN or NNW 350˚ compression [51] [52] [53] since the Himalayan movement. At the end of the Miocene (i.e. Middle Himalayan orogenic period) the Tianshan orogen began to uplift and moved southwards. As a result, the action of compressional stress led to an initial development of Dinan fault from northeast to southwest. Later, along with a strong compression from Tianshan Mountain to the south at the end of Pliocene (i.e. the late Himalaya period), the Qiulitage thrust belt experienced intense deformation process, controlled by an N-dipping imbricated fault in southern limb. During the key Neotectonic movement stage (i.e. Pleistocene) the Dfault-propagation fold rapidly uplifted and the tectonic differential stress gradually reduced from north to south, whose displacement is about 480 m after the restoration (Figure 2). Meanwhile, a South-dipping backthrust fault (i.e. Dibei fault) in the northern limb formed instantaneously to adjust the volume of tectonic deformation. Intervals of fracture in D gas field consist of Paleogene delta front rocks, which are dominated by conglomeratic siltstone, fine sandstone, siltstone and mudstone together with limited sandy conglomerate. The depth of these beds in the study area is about 4600 - 5300 m and the average total thickness can reach up to 450 m. From top to bottom the Paleogene stratum can be further divided into two groups, i.e. the Kumugeliemu Group (EII) and the Suweiyi Group (EI), both containing a set of stable interlayer in the middle. The two sets of interlayer can reach a total thickness

Figure 1. Maps show the elementary structural features of the Kuqa Depression (a) Location of the study area in the Kuqa Depression, I-Northern monocline tectonic zone; II-Kelasu tectonic zone; III-Baicheng sag; IV-Qiulitage tectonic zone; V-Yangxia sag; VI-Northern Tarim Uplift; Pre-T: Silurian to Triassic; J: Jurassic; K:Cretaceous; N₁j: Jidike Formation; N₁k: Kangcun Formation; N₂k-Q₁x: Kuqa Formation to Quaternary (from Zhang et al., 2006); (b) Tectonic cross-section of the location shown in (a); (c) Structural map of the top Palaeogene in the D gas field, with strikes of fractures obtained from orientated cores by means of the geomagnetism method.

(a) (b) (c)

Figure 2. Sketch map showing the restoration of the D fault-propagation fold.

of up to 40 m in the study area and act as the role of the slip layer. The porosity of Paleogene reservoir ranges from 5% - 15% through core tests and permeability lies in the range of (0.05 - 1.0) × 10⁻³ μm². However, the fracture permeability can reach (1.00 - 10.00) × 10⁻³ μm². On the whole, the physical property of upper EII reservoir is slightly better than that of lower EI reservoir, Finally, all the above evidences prove that the D gas field belongs to tight sandstone reservoir with low porosity and low permeability.

The most frequently encountered fractures in the reservoir of the study area are planar discontinuities which are sub-perpendicular to the bedding. The fractures can be further divided into three basic types, namely shear fracture type, tenso-shear fracture type, and tension fracture type (Figure 3(a)), respectively accounting for 54.2%, 27.7%, 18.1% of the total fracture volume. Each type is also characterized by a distinct fracture shape. For instance, the shear fracture always has more straight fracture plane and longer extended distance than the

Figure 3. Types of tectonic fractures (a) and statistics of fracture linear density (b) observed in the drilling cores and interpreted by imaging logging (FMI). The fracture development geological model (c) is established on the basis of fault interpretation, stratigraphic division, experimental measurement and fracture observation. The paleostress results were measured by acoustic emission experiments [18] [79] .

other two types, and, in most cases, it can cut through rock grains. The extension fracture, however, exposes dendritic structure and frequently bypasses rock grains with a relatively shorter distance [54] , wherever these types are observable by drill cores and microscope. Based on the orientation, the fractures present in the Kumugeliemu Group and the Suweiyi Group can be subdivided into four distinct, mutual abutting fracture sets oriented NW-SE (set I), nearly SN (Set II), NE-SW (set III) and nearly EW (set IV), among which the former two sets are mainly distributed in the fold limbs, while the later two sets are mainly located in the top of fold. Set I striking 310˚ and Set II striking 10˚ and present in the hinge and limbs, are interpreted as a regional conjugate fracture system that was interpreted as related to the NNW-oriented Himalaya compression just prior to and during initial anticline growth. Fractures striking 45˚ (set IV) and striking 95˚ (set III)approximately parallel to the fold axial trend, are found mainly within the top area and are interpreted to have formed in response to local flexural stresses and tensile stresses during folding. Additionally, some scattered tenso-shear fractures striking NE, found throughout the fold and perpendicular to bedding, may have initiated as shear fractures or deformation bands, or as joints that subsequently were sheared. In the top of fold, these fractures were reactivated as late-folding normal faults. Observations from cores show that dip of fractures mainly ranges from 75˚ - 90˚ (i.e. vertical fractures), followed by 45˚ - 75˚ (i.e. high-angle fractures) and 15˚ - 45˚ (i.e. low-angle fractures). The fracture density generally contains linear density (1/m⁻¹) and volume density (m²/m³), and the former is an important parameter to illuminate the development and distribution characteristics of fractures [55] . Observations made in the cores of the tight sandstone in the limbs and top of the anticline show obvious relationship between fractures and structural positions. Vertical fractures and echelon fractures mainly distribute in the fold hinge area, with linear density ranging from 1.2/m⁻¹ to 1.5 m⁻¹ and length and vertical extent often greater than 10 m (Figure 3(b)). Structures in two limbs are mainly composed of high-angle conjugated fractures with linear density ranging from 0.65 m⁻¹ to 1.1 m^{−1, that is characteristic of systematic joints as defined in} [56] . On the whole, the fracture density achieves the highest values in hinge area of D anticline, then has the higher values near the thrust fault zone, especially around the boundary faults. In the backlimb, the fracture density widely is higher than that in the forelimb probably because of strong back thrust faulting in the key Neotectonic movement stage. All the evidences indicate that the fracture development model in D anticline conforms with a conceptual fold-evolution mode with graben faults and tensile fractures on the hinge area [57] (Figure 3(c)).

3. Fractures and Strain-Energy Density

In order to obtain fracture density and occurrence from geomechanical modeling approach, and further analyze how the mechanical factors control the development and distribution of fractures in fault-propagation fold, a connection between the fracture density and geomechanical modeling results must be established. Accordingly, we used improved “Quasi-Binary Method”, i.e., composite failure criterion inducing strain-energy release as an indicator for fracture development, and establish a relationship between fracture density and strain energy density. And the details of this method could be found in research results of [58] . Ultimately, we numerically calculated paleostress distributions and used them to infer the distribution of fractures in the geomechanical models for comparison with core data.

When a rock mass is subjected to spatial principal stresses σ₁, σ₂ and σ₃, the strain-energy density at a given point can be expressed as [59] :

$\varpi =\frac{1}{2}\left({\sigma }_1{\epsilon }_1+{\sigma }_2{\epsilon }_2+{\sigma }_3{\epsilon }_3\right)$ (1)

where $\varpi$ is strain energy density (J/m³); and ${\epsilon }_1$ , ${\epsilon }_2$ , ${\epsilon }_3$ respectively are the strains corresponding to three principal stresses.

Recalling the Hooke’s Law relations and substituting for the strain components into [60] , we have the most general state of stress.

$\varpi =\frac{1}{2E}\left[{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2-2\mu \left({\sigma }_1{\sigma }_2+{\sigma }_2{\sigma }_3+{\sigma }_1{\sigma }_3\right)\right]$ (2)

where $\mu$ is the Poisson’s ratio.

Deep tight sandstone is generally characterized by strong brittleness, high elasticity modulus and low Poisson’s ratio. According to brittle fracture mechanics theory and maximum tensile stress theory, when elastic strain energy release rate accumulating in brittle rocks equals to the energy per unit volume for generating fractures, rocks will break down. When the surrounding three-dimensional stress state reaches the rock strength, macro brittle fractures will occur with strain energy releasing, part of which will offset the surface energy of new fractures, and the rest of which will offset in form of elastic waves [61] . However, for fractures or fractures, the elastic wave energy is so weak that it can be completely ignored. We assumed that all of fractures in rocks are caused by the releasing energy. If ${\varpi }_f$ is regarded as the residual strain energy density by current strain energy density per unit volume subtracting elastic strain energy density for new fractures( ${\varpi }_e$ ), a formula keep to the principle of energy conservation will be given by:

$D_{vf}=\frac{S_f}{V}=\frac{{\varpi }_f}{J}=\frac{\varpi -{\varpi }_e}{J}=\frac{1}{J}\varpi -\frac{{\varpi }_e}{J}=a\varpi +b$ (3)

where $D_{vf}$ is fracture volume density per unit volume (m²/m³); ${\varpi }_f$ is strain energy density for new fractures (J/m³); V is characterization of unit cell volume (m³); S_f is surface area of new fractures (m²); $J$ is the required energy per unit area for fractures (J/m²), i.e., fracture surface energy (here the energy is different from and far less than theoretical value of molecular dissociation); ${\varpi }_e$ is the necessary elastic strain energy to be overcome for new fractures (J/m³); and a and b are the relative coefficients.

Furtherly, relationships between fracture volume density and strain energy density under uniaxial stress state, triaxial stress state will be established after formula transformation, the final fracture linear density is derived as following:

1) Commonly, tectonic fractures are divided into tensile and shear fractures based on the forming mechanism, and they can be discriminated with the Griffith Criterion and Coulomb-Navier Criterion (also called the Coulomb-Mohr Criterion), respectively [11] [27] [42] [44] [62] [63] [64] [65] [66] . If there exist only compressive stresses ( ${\sigma }_3\ge 0$ ), the Mohr-Coulomb criterion will be selected, that is

$\frac{{\sigma }_1-{\sigma }_3}{2}\ge C_0\cos \phi +\frac{{\sigma }_1+{\sigma }_3}{2}\sin \phi$ (4)

The Coulomb-Mohr Criterion suggests that a shear fracture only forms if the internal strength or cohesion of the rock (C₀) is exceeded, and depends on the magnitude of normal stress along a fracture plane. The relationships between fracture density, aperture and strain energy density, stress-strain are written as follows

$\{\begin{array}{l}{w}_f=w-w_e=\frac{1}{2E}[{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2-2\mu \left({\sigma }_1+{\sigma }_2+{\sigma }_3\right)\\ -{{0.85}^2{\sigma }_p^2+2\mu \left({\sigma }_2+{\sigma }_3\right)0.85\sigma ]}_p\\ {\sigma }_p=\frac{2C_0\cos \phi +\left(1+\sin \phi \right){\sigma }_3}{1-\sin \phi }\\ E=E_0{\sigma }_p\\ D_{vf}=\frac{w_f}{J}\\ J=J_0+\Delta J=J_0+{\sigma }_{3}b\\ D_{lf}=\frac{2D_{vf}{L}_1{L}_3\sin \theta \cos \theta -L_1\sin \theta -L_3\cos \theta }{L_1^2{\sin }^2\theta +L_3^2{\cos }^2\theta }\\ D_{lf}b=\left|{\epsilon }_3\right|-\left|{\epsilon }_0\right|\\ {\epsilon }_0=\frac{1}{E}\left({\sigma }_3-\mu \left(0.85{\sigma }_p+{\sigma }_2\right)\right)\end{array}$ (5)

where $\phi$ is the internal friction angle (˚); $\theta$ is the angle between normal to the newly formed fracture plane and the maximum principal stress (˚); ${\sigma }_1$ is the peak stress (MPa); ${\sigma }_2$ is the intermediate principal stress (MPa); ${\sigma }_3$ is the middle principal stress, (MPa); ${\sigma }_p$ is the rupture stress under action of ${\sigma }_3$ , different from the maximum principal stress; ${\sigma }_t$ is the tensile strength, (MPa); $J_0$ is fracture surface energy with no confining pressure or under unaxial compressive stress (J/m²), $\Delta J$ is the additional surface energy caused by confining pressure ${\sigma }_3$ (J/m²); E is the elasticity modulus with no confining pressure (GPa); b is the fracture aperture (i.e. paleo-aperture, here for reference only), m; $D_{lf}$ is the fracture linear density (1/m); ${\epsilon }_3$ is the tensile strain under current state of stress, dimensionless parameter; ${\epsilon }_0$ is the maximum tensile strain, dimensionless parameter, corresponding to tensile strain when crack beginning to form; $E_0$ is the proportionality coefficient related to lithology; and L₁, L₂, L₃ are side length of the selecting representing element volume (REV) (refer to Figure 4).

2) When there exists tensile stress, for brittle tight sandstone material the Griffith Criterion is used, that is

When ${\sigma }_3\le 0$ and $\left({\sigma }_1+3{\sigma }_3\right)\ge 0$ , the applied failure criterion is given as

${\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\ge 24{\sigma }_T\left({\sigma }_1+{\sigma }_2+{\sigma }_3\right)$ (6)

$\cos 2\theta =\frac{{\sigma }_1-{\sigma }_3}{2\left({\sigma }_1+{\sigma }_3\right)}$

Figure 4. In-situ stress coordinate system representing fracture distribution and differential stress based on REV. (a) An Representing Element Volume (REV) is selected to establish the relationship between stress and fracture parameters under complex stress condition, and some hypotheses are made as follows: it is small enough to be easily cut through, the scattered microcracks inner element can be negligible and the element is supposed as a parallelepiped on with sides L₁, L₂, L₃, thus the σ₁ corresponding to side L₁, the σ₂ corresponding to side L₂, the σ₃ corresponding to the side L₃; (b) Transection in REV perpendicular to σ₂, namely (σ₁ - σ₃) plane, here, θ is the acute angle between fracture surface and maximum principal stress σ₁.

When $\left({\sigma }_1+3{\sigma }_3\right)<0$ , the failure criterion is simplified to ${\sigma }_3=-{\sigma }_T$ , $\theta =0$ .

If the failure criterion is reached, the relationships between fracture density, aperture and strain energy density, stress-strain are expressed as

$\{\begin{array}{l}{w}_f=w-w_e=\frac{1}{2}\left({\sigma }_1{\epsilon }_1+{\sigma }_2{\epsilon }_2+{\sigma }_3{\epsilon }_3\right)-\frac{1}{2E}{\sigma }_t^2\\ D_{vf}=\frac{w_f}{J}\\ J=J_0+\Delta J=J_0+{\sigma }_{3}b\\ D_{lf}=\frac{2D_{vf}{L}_1{L}_3\sin \theta \cos \theta -L_1\sin \theta -L_3\cos \theta }{L_1^2{\sin }^2\theta +L_3^2{\cos }^2\theta }\text{or}{D}_{lf}=D_{vf}\\ D_{lf}b=\left|{\epsilon }_3\right|-\left|{\epsilon }_0\right|\\ \left|{\epsilon }_0\right|=\frac{{\sigma }_t}{E}\end{array}$ (7)

The above fracture mechanics models are all under local coordinates. However subsequent simulation of stress field and calculation of fracture parameters need to be transferred into global coordinates. In three-dimensional space, conformation of fracture strike and dip depend on reasonable projection method, e.g. X-axis agrees with east-west direction, Y-axis agrees with north-south direction, and Z-axis agrees with vertical direction. If direction cosine under global coordinate of normal direction vector of fracture plane is determined as $\bar{n}=\left\{\begin{array}{ccc}l& m& n\end{array}\right\}$ , and $\bar{n}$ is projected to x-o-z plane with angle between projection line and negative direction of z-axis, through the angle conversion ${\alpha }_Z=\arctan \left(-l/n\right)$ , the angle of strike $\alpha$ will be calculated.

If $0\le {\alpha }_Z<{90}^{\circ }$ ，then $\alpha ={90}^{\circ }-{\alpha }_Z$ (8)

If ${90}^{\circ }<{\alpha }_Z<0$ , then $\alpha =\left(-{90}^{\circ }-{\alpha }_z\right)+{360}^{\circ }$ (9)

From the geological point of view, fracture dip ${\alpha }_{dip}$ between x-o-z plane and fracture plane also can be defined as angle between the $lx+my+nz=0$ plane and $y=0$ plane, that is ${\alpha }_{dip}$ ( $0^{\circ }\le {\alpha }_{dip}\le {90}^{\circ }$ ), which is expressed as

$\cos {\alpha }_{dip}=\frac{\left|l\cdot 0+m\cdot 1+n\cdot 0\right|}{\sqrt{l^2+m^2+n^2}\sqrt{0^2+1^2+0^2}}=\frac{\left|m\right|}{\sqrt{l^2+m^2+n^2}}$ , $0^{\circ }\le {\alpha }_{dip}\le {90}^{\circ }$ (10)

4. Geomechanical Simulation of Tectonic Stress Field

The simulation approach used in this study utilizes the Finite Element (FE) technology, and geomechanical models will be run in order to simulate the distribution of tectonic stress for further predicting fracture development and distribution. This powerful tool allows for robust simulations of complex structures with non-linear material behavior or large deformation based on frictional contact mechanics [24] [67] . The basic concept behind FEM is that the geological bodies are discretized into finite continuous elements connected by nodes. Each element is allocated with appropriate geomechanical parameters determined for real rocks. The continuous field function for the region is first transformed into node function that incorporates the basic variables of stress, strain and displacement resulting from applied external forces [68] [69] . All these elements are combined to obtain tectonic stress field over the entire geological bodies [69] [70] . The general-purpose finite element code Ansys (version 15.0) was selected for this study because it is well suited to analyzing geomechanical problems over awide range of scales in one, two and three dimensions [4] [32] [71] . Ansys can accurately handle large strains and rotations as well as complex contact interfaces with frictional behavior where significant sliding can occur. It also has efficient algorithms for solving highly nonlinear problems that result from both geometric complexity and material behavior. Finally, Ansys has a large library of built-in constitutive relationships that are appropriate for simulating the behavior of rock, ranging from simple elastic material models to advanced elastic-plastic and visco-elastic material models [4] .

4.1. Geometry and Boundary Conditions

The purpose of this step is to set up several suites of geomechanical models and experiments on the base of fold evolution and analyze the development and distribution of tectonic fractures in fault-propagation fold. The initial 2D geometric model (Figure 5(a)) was constructed based on the restored geologic section constrained by field measurements during Pliocene Kuqa Period (Figure 2) and incorporated the generalized mechanical stratigraphy of the Kumugeliemu Group (EII) and the Suweiyi Group (EI). The three sand stone members within the FE model are all set to an average thickness of 90 m, they are the first Sandstone Member of EI, the second Sandstone Member of EI, the third Sandstone

Figure 5. Initial structural geometry, boundary conditions and meshing grid of the finite element model. Red frame represents fault, green lines represents frictional sliding interfaces. Segment 1 is the Upper Sandstone Member of EI, Segment 2 is the Middle Mudstone Member of EI, Segment 3 is the Lower Sandstone Member of EI, Segment 4 is the Upper Mudstone Member of EII, Segment 5 is the Lower Sandstone Member of EII.

Member of EII from top to bottom, respectively. The average thickness of the two mudstone interlayers are all set to 70 m, they are the first Mudstone Member of EI and the first Mudstone Member of EII from top to bottom, respectively. The target zone associated with the upper cover layers (i.e. Q. N. Formations) and basement (i.e. K. J. Formations)are considered as an isolated bodies or boundary conditions for simulation. The conformation of boundary conditions is important to demarcate the study area from the rest parts of Kuqa Depression, and eliminate the influence of boundary effect on the simulation results. Considering the Middle Himalaya Period after Miocene Jidike Formation sedimentary and the Late Himalaya Period after Pliocene Kangcun Formation sedimentary as the crucial time of fracture formation in the Paleogene reservoir, the cover layers are set to approximately 1200 m thick and 3400 m thick, respectively. In order to guarantee the accuracy of the simulation, the width of the 2D geometric model is set to 50m according to analogous theory of test model [72] [73] . Given the complexity of this stratigraphy, and the current limitation of FE models, we chose to model a six-layer composite of initially constant thickness that represents presumed layered ductile-brittle sequences. This idealized configuration does not explicitly represent any specific formation or part thereof, but we suggest it does provide an analogous mechanical stratigraphy, while honoring overall structural geometry of the natural example [20] .

And we mainly analyzed the contemporary faulting in this fault-controlled anticline. The southern boundary fault (i.e. Dinan fault) in the geological model began to develop mostly at the end of Miocene, and strongly reactivated at the end of Pliocene associated with stresses redistribution. In the Neotectonics Period, the same time the back thrust fault (i.e. Dibei fault) in the back limb rapidly formed as an accommodation structure along with rapid uplift of the D anticline and tending towards stability. In contrast, most fractures in D anticline were mainly controlled by the boundary Dinan fault in the forelimb, while were partly reactivated or influenced by the late-folding thrust fault (i.e. Dibeibackthrust fault) [53] [74] . According to the geological mechanics theory, these approximate parallel fault sets striking NEE70˚ indicate a NNW330˚ compression perpendicular to the northern boundary frame [10] [48] . In the model, the fault cannot propagate but is allowed to displace under tectonic stress. Bedding-pa- rallel frictional sliding interfaces based on a classic Coulomb friction model were included between each layer so that configurations with and without interlayer slip could be analyzed. For the case where slip is allowed, a friction coefficient of 0.1 was assigned to the fault interfaces and a friction coefficient of 0.25 was assigned to the interlayers. For the case where slip is prevented, the interfaces are tied so that no sliding can occur (Figure 5).

Generally, the vertical stress can be calculated from the bulk density of rocks based on Equation (11), and the horizontal stress component at the end of Miocene caused by the bulk density was less than 3 MPa in the Suweiyi Group and Kumugeliemu Group calculated by Equation (12) with an average depth of 1200 m. In the same way, the horizontal stress component at the ends of Pliocene and Pleistocene caused by the bulk density was less than 5 MPa.

${\sigma }_z=\rho hg$ (11)

${\sigma }_h={\sigma }_z{\left(\frac{\mu }{1-\mu }\right)}^{\frac{1}{n}}$ (12)

where ${\sigma }_z$ is the vertical stress, ${\sigma }_h$ is the horizontal stress caused by the bulk density, μ is the Poisson’ ratio, ρ is the density, g is the gravitational acceleration, and n is a constant related to the nonlinear compression and is 0.67 here [66] [75] [76] . Based on regional analysis and acoustic emission of rocks, the study area has experienced three important thrusting movements during the Himalayan Periodand Neotectonics Period, whose differential stress values are 52 MPa, 104 MPa and 72 MPa, respectively [77] . The average differential stress acting on the D anticline is calculated as about 84 MPa. The maximum/mini- mum principal stress ratio ( ${\sigma }_1/{\sigma }_3$ ) of about 2.1 for shallow crust (<4000 m; Gao, 2011) enables the principal stresses to be calculated in D gasfield. Hence, in the present study, a maximum principal stress of (180.4 + ${\sigma }_h$ ) MPa is applied in the NNW direction (i.e. x-direction)and a minimum principal stress of (66.4 + ${\sigma }_h$ ) MPa is applied in the NEE-SWW direction during modeling (Figure 5). The entire model is subjected to gravity loading in the vertical direction, which could be automatically applied in the Ansys software. Besides, some appropriate displacement constrains are applied to the geological model to prevent it from rotation and undergoing rigid displacement, and to consequently facilitate simulation. The top portion of the geometry model is set as a free surface, whereas its bottom is fixed in vertical and southern edge as fixed in the horizontal direction. Additionally, is the horizontal the maximum principal stress of thrusting is not applied on the basement. It is provided that compressive stresses were positive and tensile stresses are negative in this study.

Note that our intention is to precisely simulate the overall deformation history of the D anticline under horizontal compressive stresses and capture real-time development and distribution of strain energy density and fractures. Based on the sustainability of tectonic movement on the earth, each thrusting movement can be subdivided into two critical stages, so that a total of 6 deformation stages numbered in geological models are conducted through finite element analyses. From beginning to the end, a horizontal stress of 50 MPa, 60 MPa, 80 MPa, 100 MPa, 90 MPa, and 70 MPa from northern Tianshan Mountain is applied to the geological model.

4.2. Material Properties

Material properties are assigned to the elements representing the various lithologies. The finite element models can describe elastic and plastic rock deformation to adequately simulate the folding behavior. The exposed strata and corresponding lithology can be obtained from the stratigraphic data of the Kuqa Depression. To begin the examination of macroscopic effects, the geological model of this study was divided into the following major types of unit: the sandstones, mudstone interlayers, basement, and fault. The mechanical parameters of the EI group and EII group are assigned with the different values because of lithologic differences respectively, which had certain effect on the simulation results. These parameters were determined by the testing of rock mechanics (Table 1). In the present work, a total of 34 samples from the drillcores of EI group and EII group were collected for rock mechanics and other analysis. The various rock densities were determined by density analysis, the Young’s moduli and Poisson’s ratios were acquired from uniaxial compressive tests, the tensile strength values were obtained through splitting tests, and the internal friction angle and cohesion were from triaxial compressive tests. The confining pressure that were employed in the triaxial rock mechanics testing were 0, 10, 20, 30, 40, 50, and 60 MPa, and the role of ground fluid was not considered. The average rock mechanics parameters in the study area are listed in Table 1.

Table 1. Rock mechanical parameters in the EI and EII of D gas field.

The definition of mechanics parameters in fault zones is extremely significant to the modeling results, however exact mechanical parameters for fault zones are unavailable. Based on previous studies [69] [78] [79] , in this FE model, faults were defined as weakness zones with Young’s moduli about 60% of the corresponding sedimentary layers. Commonly, Poisson’s ratios in fault zones were larger than those of the corresponding sedimentary rock stratum, and their differences were typically between 0.02 and 0.10. The parameters of internal friction angle, cohesion and tensile strength for fault zones within the FE model were defined as 60% of the corresponding mechanical parameters of rocks (Table 1).

4.3. Mathematical Model

The Ansys software was used to construct a geological model that was meshed based on area distribution, rock structure, and lithology distribution of the target stratum. The geological model was meshed in the form of hexahedron elements with 8 nodes by an artificial mesh-generation method. In general, the fault zones and targeted layer are divided by finer elements, whereas other layers were meshed by coarse elements. The geological units of the whole model are meshed by setting the element length from small to large. Each layer was meshed and then connected with one another in space for a 2D mesh model. The model contained a total of 12,576 elements and 2796 nodes (Figure 5). The fault zone and target stratum are meshed with a fine grid, and the rest of the area was meshed with a relatively coarse grid. The mesh refinement steps of the model can be described as follows: Firstly, the model is meshed with a coarse grid and hexahedron elements, and then the meshing result could be obtained. Secondly, the accuracy of the solution is computed, and whether the grid can be accepted was determined. Finally, steps one and two are repeated until the accuracy of solution satisfied the inversion standard.

4.4. Calculation of Solution

After establishing several suites of mathematical models, we used the Ansys software to simulate the paleotectonic stresses fields during the Himalayan Period and Neotectonics Period in the D area, Kuqa Depression. The simulation results included the distribution of the maximum principal stress, minimum principal stress intermediate stress and differential stress, which provided stress parameters for fracture prediction (Figure 6). In the Suweiyi Group and Kumugeliemu Group of D gas field, the σ₁, indicative of compression, range from 140 MPa to220 MPa. The σ₂ (i.e. vertical stress) indicative of compression, are generally between 56 MPa and 123 MPa, and the σ₃, indicative of both compression and tension, ranges from −30 MPa to 15 MPa (Figure 6(b)). The differential stress by σ₃ subtracting σ₁, indicative of possible damage zone, mainly ranges from 0 to 85 MPa (Figure 6(a)). The pre-existing fault tip show a higher differential stress, because with the regional contraction and the increasing displacement

Figure 6. Stress intensity (a) and the minimum principal stress (b) in progressive deformation for fault-related folding model with various loading. The interlayer slip and fault plane slip are permitted to facilitate the continuous uplift of fold. The structural stability and the stress concentration with discontinuous contours indicating probable damage zone during the deformation process is shown after (a). The positive of minimum principal stress in (b) reflect locations where extension is occurring. Compressive stresses are considered as positive and tensile stresses are negative.

along the fault plane and uplift of fold amplitude, the stress and strain gradually focus in the thrust fault. Along with the shorting and folding of D anticline, the differential stress gradually increases and the positive values of minimum principal stress migrates from the back limb to the forelimb, indicating a successive transition from regional stress field to local tensile stress field in the core of fold. For comparison, considering the mechanical differences between sandstone member and mudstone interlayers, a series of inconsistent deformation phenomena occur near the lithologic interfaces, accompanied with obvious bed-parallel slip, wide development of tensile stress in sandstones and concentration of differential stress. The distributions of three principal stresses are similar to each other, all of them are mainly fold-controlled and secondly fault-controlled. In addition, both the highest values of tensile stress and differential stress are all located in the hanging wall of fault and top areas of D gas field (Figure 6).

5. Prediction of Tectonic Fractures

The simulation of the 2D tectonic stress field may be used to calculate the tensile stress and compressive stress of rocks during tectonic movement. The rock rupture criteria in Equation (4) and Equation (6) may then be employed to determine whether ruptures will occur in rocks. $D_{vf}$ , $D_{lf}$ and ${\alpha }_{dip}$ are used to express the development degree of the tensional and shear fractures comprehensively. According to Equation (5) and Equation (7), the fracture parameters related to stress, strain and strain energy density are calculated or directly extracted from the above numerical simulation of paleotectonic stress field. Firstly, based on abundant mechanical experiments [58] , it is found that when axial stress reached or exceeded 0.85σ_c (σ_c is compressive strength of rock), the number of fractures explode instantaneously and the number of large-scale fractures increase faster than that of small-scale fractures, which can be considered as damage threshold indicating a key stage of abundant microcracks coalescence and macrocrack initiation. At this key stage, the corresponding strain density energy is close to ${\varpi }_e$ in theory, hence, we assume that ${\varpi }_e$ is strain energy density at $\sigma =0.85{\sigma }_c$ . Recalling the Hooke’s Law the ${\varpi }_e$ and $\varpi$ can be obtained [60] , then $w_f$ is calculated. According to [58] , for low-permeability sandstone the corresponding average $J_0$ (fracture surface energy under unaxial stress state) is 1087.35 J/m², and the average $E_0$ is 112.6. Secondly, putting $J_0$ and $E_0$ into the Equation (5) and combining with formulas of $D_{vf}$ , $D_{lf}$ , b and $J$ , the fracture volume density $D_{vf}$ , the fracture linear density and fracture aperture b under triaxial state of stress can be calculated, in that order. Thirdly, considering appearance of tensile stress, the fracture parameters $D_{vf}$ , $D_{lf}$ and ${\alpha }_{dip}$ are obtained using the same method, which is built into the Ansys platform to calculated spatial distribution of fracture parameters (Figure 7 and Figure 8).

Figure 7. The distribution sections of fracture linear density during key evolutionary stages. At every step, different external stress coming from the Southern Mountain was imposed on Paleogene strata, a series of continuous elastico plastic deformations along fault plane were simulated and the anticline raise continuously along with the constant migration and changing of fracture zones.

6. Results and Discussions

To gain insight into the stress field and damage zone during fault-related fold- ing, we conducted six simulations for different combinations of mechanical stratigraphy (thickness and mechanical properties of the rock layers), continuous contraction of the fold, and strength of bedding surfaces. In sum, the final geomechanical finite element modeling result during the six step matches well with cores observation and imaging logging interpretation results both in the whole geometry and facture distribution (Figure 3), which implies that the Strain Energy Density and Composite Failure Criterion used here are good measures of fracture development and distribution. Therefore, we can further study how the different mechanical factors control and influence the development and evolution of fractures in the fault-propagation fold on basis of the geomechanical finite element models, and probably predict the spatial distribution of fracture parameters (i.e., fracture density, aperture, length, dip and strike).

6.1. Effect of Stresses on Fracture Development and Distribution

Commonly, the development of tectonic fractures is generally controlled by

Figure 8. The distribution sections of fracture dip during key evolutionary stages.

tectonic stress field [80] [81] [82] . Due to remote effects of the Cenozoic India-Eurasia collision and subsequent continuous compression, the D gas field experienced a NNW-SSE compressional regime during the Himalayan period in the Kuqa Depression. During the pre-folding stages, such as Step 1 and Step 2 in Figure 6, the compressive stress and shear stress were predominated in the whole area, a large number of shear fractures probably generated if the internal strength of the rock was exceeded and the Coulomb-Mohr Criterion was met [65] . At the same time, according to the theory of Cosseratt Medium Constitutive Model for laminated rocks [83] [84] , the uncoordinated deformation and parallel-slip tendency will occur along the lithologic interfaces and spread laterally in both directions. If viscoelastic plasticity of mudstone interlayers is stronger, flexural slip interacted between layers is more important for complex deformation. As a surface’s mechanics response, tensile stresses will be derived and concentrated around the interfaces of laminated rocks, which subsequently yielded overall increasing of differential stress. Once the concentration at a point near the peripheral edge of sandstones reaches the cohesive strength, the material will begin rupturing [64] . Therefore, based on the above analysis, the development of shear and tensile fractures can be determined by increasing of differential stress and appearance of tensile stresses. During the rapidly folding stages, such as Step 3, Step 4, Step 5 and Step 6 in Figure 6, the high-value areas of differential stress and tensile stresses gradually migrated from the limbs to the hinge and from the bottom to the top of anticline, leading the generation of large numbers of fractures (Figure 7 and Figure 8).

For all simulations in the geological model, the horizontal pressure of each key construction period is set to 50 MPa, 60 MPa, 80 MPa, 100 MPa, 90 MPa, and 70 MPa in turn. For further study of continuous influences on development of fractures in fault-propagation fold, the horizontal pressure of 30 MPa, 110 MPa, 120 MPa and 130 MPa is respectively added to the geological model as a supplement. The stresses intensity is computed at the center of elements, which is identified in Figure 6 with various color zones. A logarithmic relationship is established between fracture linear density and corresponding differential stress (Figure 9), and the formula is

$y=0.7426\ln x-1.1603\left(R=0.9572\right)$

where y is the simulated fracture density in model (m⁻¹), x is the differential stress (i.e., σ₁ − σ₃), and R is the correlation coefficient.

As is shown in Figure 10(a), it is generally believed that higher differential stress is able to promote the development of tectonic fractures by causing stress concentration in formations and fault tips, or contribute to the creation of associated and derived fractures. However, along with the increasing differential stress the fracture density firstly increases and then begins tending towards stability at an approximate threshold value of 110 MPa, which indicates that the number of fractures in tight sandstones cannot grow without limit as the pressure increases constantly.

Besides, the simulated developed areas of fracture density coincide with the areas with high tensile stress, including the hinge of anticline and the northerns and stones near the lithologic interfaces (Figure 6 and Figure 7), which indicates that the appearance of local tensional stresses promotes the generation of extensional fractures superposing on early fracture networks.

Figure 9. Relationship between the stresses and corresponding fracture development intensity. (a) Fracture linear density vs differential stress; (b) Fracture linear density vs tensile stress.

Figure 10. Relationship between the slip displacement, uplift height and corresponding fracture development intensity, dip. (a) Fracture linear density vs uplift height; here, the data points are extracted from six simulations and when the point location is selected the vertical height between base level of the sandstone formation and the sample point is calculated immediately; the threshold value is approximately confirmed by a transition of fracture density from high values to low values; (b) Fracture linear density vs distance from fault surface; here, the perpendicular distance from sample point in model to the fault surface represents the x-axis; (c) Fracture linear density in hinge position vs slip displacement along fault surface; the filled circles indicate average fracture density of top EI in different folding stages, and the hollow circles indicate average fracture density of bottom EII in different folding stages; (d) Fracture dip variation in hinge position vs uplift height; the filled circles indicate fracture dip of EII in different folding stages, and the hollow circles indicate fracture density of EI in different folding stages.

6.2. Effect of Slip Displacement or Uplift Amplitude on Fracture Development and Distribution

Of interest here is the relation between different slip displacements and uplift amplitude with fault-parallel slip. Generally, vertical amplitudes that develop during folding are affected by the interaction between slip displacements along fault plane and differential stress with different mechanical properties [85] [86] [87] . Given the aforementioned mechanical constraints, the basement remains stable but the overlying formations on upper wall will slide along the fault interface and yield overall horizontal contractions, 40 m, 170 m, 320 m, 650 m, 840 m and 1020 m, respectively. According to the Fault-related Fold Theory [88] [89] including the fault terminal displacement conversion and the front-ent propagation of the upper wall, the uplift amplitudes (i.e., 20 m, 70 m, 140 m, 230 m, 310 m and 430 m respectively) along the hanging wall increase gradually, and the crestal roundness decreases with persistent slip displacements. For comparison, because the slip displacement in the Step 4 is the largest (i.e. 330 m) and the corresponding uplift height is relatively larger (i.e. 90 m), the folding strength plays a more prominent role in the fracture generation. The results agree with previous numerical [90] and experimental results [91] .

The influence of slip displacement and fold amplitude on the development and distribution of fractures in fault-propagation fold is observed and discussed (Figure 10). A logarithmic relationship is established between fracture linear density and corresponding uplift amplitude due to successive stages of a simulation (Figure 10(a)), and the formula is

$y=0.839\ln x-2.6501\left(R=0.8823\right)$

where y is the simulated fracture density in model (m⁻¹), x is the uplift height (m), which is equivalent to curvature in folding process, and R is the correlation coefficient.

As is shown in Figure 10(a), along with increasing of uplift height from the base level of sandstone formation, the fracture density rapidly increases until the peak point (approximately 2.37 m⁻¹) in initial segment of the curve then slowly increases until the end. The corresponding height to this peak point of fracture density is about 240 m, which is regarded as the most critical stage of fracture development. However, as in Figure 10(b), the phenomena of several data points deviating from main trend, probably implies that the fracture development and distribution are results of combination effects by faulting and folding.

At the same time, a power relationship is barely established between the fracture linear density and corresponding distance from fault surface (Figure 10(b)), and the formula is

$y=9.3044x^{-0.482}\left(R=0.7749\right)$

where y is the predicted fracture linear density in model (m⁻¹), x is the distance from fault (m), and R is the correlation coefficient.

As illustrated in Figure 10(b), along with increasing distance from fault zone the fracture density rapidly reduces until reaching a stable value (about 0.3 m⁻¹), and in this process some abnormal points seriously deviate from the curve track indicating a minor quantitative effect of fault activity on development of fracture density.

For comparison, another integrated logarithmic relationship is confirmed between the fracture linear density and corresponding slip displacement as follows (Figure 10(c)).

$y=0.7309\ln x-2.6306\left(R=0.9481\right)$

where y is the predicted fracture linear density in model (m⁻¹), x is the slip displacement along fault during folding (m), and R is the correlation coefficient.

As is shown in Figure 10(c), along with successive increase of slip displacement along fault surface, the average fracture density steadily in EI and EII formations increases until reaching a peak point (about 2.49 m⁻¹). Considering the difference of buried depth, within the first three stages the high-value areas of fracture density are located in upper EI sandstones, however, in the late three stages the high-value areas are transferred to the lower EII sandstones, which is similar to the results of hybrid cellular automata numerical technique [24] [87] .

For the model results described above, the spatial variation of fracture dip during folding is a result of stress state transformation under sustaining loadings. As with an elastic-plastic model, the Mohr-Coulomb and Griffith criterions give the different rupture angles in extension and compression. Therefore, an experimental linear relationship is established between the fracture dip and corresponding uplift height of folding as follows (Figure 10(d)).

$y=0.1684x+11.33\left(R=0.503\right)$

where y is the predicted fracture dip in model (˚), x is the uplift height during folding (m), and R is the correlation coefficient.

It should be emphasized that this obvious linear relationship is only occurred in upper EI formation (Figure 10(d)), whereas in the bottom EI formation there shows a parabola relationship, which indicates the dipping region resulting from interactions among multi-factors (e.g., lithology, depth, stress, fault and fold).

As mentioned above, the fractures observed in cores and identified in FMI can be divided into four sets, and the fracture Set I and Set II mainly presenting in hinge and limbs, are interpreted as a low-angle conjugate fracture system at initial folding stage. Along with increase of depth, the compressional stress state gradually converts to typical strike-slip stress environment, which usually is marked by a series of strike-slip fractures. This result agrees with the simulated distribution characteristics of fracture density and dipping region as shown in first three stages (Figure 7 and Figure 8). At the completion of fold growth, predicted fracture dip on the fold top is about twice that predicted for the fold core. If fold continuously rise along with increasing slip displacement, extensional stress is significantly increasing. Further, the density of high-angle tensile fractures (i.e., set IV and set III) approximately parallel to the fold axial trend begins to get larger than that of shear fractures and eventually will be far more than the later.

6.3. Effect of Lithology on Fracture Development and Distribution

In general, the lithology influences the development of tectonic fractures in reservoirs basically [16] [82] , which can be reflected by the differences in the mechanical parameters. An increase in the proportion of brittle minerals will decrease the tensile strength of rocks and facilitate the generation of tectonic fractures under the actions of external forces [30] . The lithology and corresponding mechanical parameters varied in the Paleogene formations, leading to the difference in the development and distribution of tectonic fractures (Figure 11). As is shown in Figure 7, initially, when the sandstones respond elastically, high-value of facture density are mainly located in Sandstone Formation of EII due to original higher elastic modulus and/or lower Poisson’s ratio in other formations. As the fold grows in association with greater shortening (e.g., slip displacement is 650 m) the total thickness of stratum at fold hinge zone increases, the confining pressure magnifies rupture limit and promotes elasticoplastic

Figure 11. Lithologic effect on development of tectonic fractures. Fractures containing macrocracks and microcracks are observed in drilling cores and thin sections; the y-axis is defined as fracture linear density; microcracks in glutenite and conglomeratic siltstone (i.e., inner-grain fractures and grain boundary fractures) account for more than 38% and 40% of total tectonic fractures, respectively.

deformation [92] , as a result the fracture density in Sandstone Formation of EII decreases rapidly. In contrast, high-value areas, where tensile stress and differential stress rapidly concentrate, gradually transform from the lower sandstones to Upper Sandstone Formation of EI due to later effect of differential stress.

Statistical analysis from cores and thin sections, shows that lithology still has effect on fracture development and distribution to a certain extent in anticline. As is shown in Figure 11, along with the average grain size of rocks becoming coarser and the sorting getting better, the fracture density increases overall although most fractures in glutenite and conglomeratic siltstone are identified as inner-grain fractures and grain boundary fractures. For example, during initial folding stage, density of netted fractures in Sandstone Formation of EII is obviously higher than that of EI, which is primarily attributed to higher content of coarse-grained components and better sorting.

7. Conclusions

Our non-linear 2D FE models of a fault-propagation fold (the D gas field) developed in a mechanically stratified sequence of the Kuqa Depression are analyzed to understand the development and distribution characteristics of fractures during different folding stages. As a result, the overall geomechanical FE modeling results are in agreement with field observation results both in the geometry and in the fracture development and distribution. Clearly, the Strain Energy Density and Fracture Density used here as indicators are not appropriate for all types of fractures, however, if combined with fracture criterions, it might be useful indicators in the present study and the reliability of this premise is acceptable. Certainly, aiming at more complicated structures, the improved and elaborate anisotropic geomechanical models including various lithologies and relatively true structural feature should be constructed for more accurate analysis based on more advanced software platform and theory.

Several important geological factors, such as the differential stress, tensile stress, distance from fault, slip displacement, uplift height and lithology are discussed and analyzed quantitatively based on the simulated results and measured results. A logarithmic relationship between fracture linear density and stress density, fracture density and uplift height, fracture density and slip displacement was established respectively. However, unlike the formers, relationship between fracture linear density and tensile stress values accords with good linear trend, which implies that the derived tension stress field during folding has the most influence on fracture development. Moreover, a negative power relationship between fracture density and distance to major fault surfaces is roughly established through scattered data points, which indicates during the Early Initial Compressional Stage and Mid-term Strong Thrusting and Uplift Stage, most fractures probably are co-controlled by folding faulting activities. Additionally, the potential evolutionary relationship between fracture dip and uplift height during folding is studied. On the whole, along with the increasing of fold amplitude the facture dip becomes steeper, but some fractures with shallower dip exist in the fold core and near the mudstone interlayers, oblique to the fold trend due to widely bed-parallel slip. Lithology as an innate factor influencing the mechanical properties of rock and possibility of fracturing, has non-linear relationship with fracture linear density. In terms of microcosmic, along with the grain size of minerals becoming coarser and the sorting getting better, the fracture density increases although most fractures in glutenite and conglomeratic siltstone belong to inner-grain fractures and grain boundary fractures. Therefore, these geomechanical factors altogether control and influence the development and distribution of tectonic fractures during evolutionary stages of the fault-propagation fold, and the relative height of fold uplift has the largest slope value among all the factors, which implies that the fold uplift in the most important factor for fractures. The second geomechanical factor influencing the development and distribution of fractures in high thrust fault-propagation fold is slip displacement along fault, and the third is lithology, the last is distance to fault plane. Here it must be emphasized that crosscutting relationships lie among these geomechanical factors, such as the derived stress factors are the outcome of fold contraction and uplift.

Finally, this evolutionary model of the structures present to us is a distinct and intuitive top-graben fold with high amplitude and high-density factures in the top, has experienced strong folding and uplift process. From early to late, in profiles both high-density and high-angle areas of fractures migrated from south to north, from deep to shallow and from the limb to the top, and fractures begin to generate in the sandstones, then gradually extend into the mudstones rather than unboundedly clustering in sandstones. Moreover, when the uplift height of fold along with corresponding slip displacement exceeded approximately 55 percent of the total height the early fractures on top parallel to the axis of fold were reactivated until developed into normal faults, which undoubtedly reduced the later density of tensile fractures at present.

Acknowledgements

This research was financially supported by the National Oil and Gas Major Project (2016ZX05047-003,2016ZX05014002-006), the National Natural Science Foundation of China (41572124) and the Fundamental Research Funds for the Central Universities (17CX05010). The authors wish to thank master student Naser lovely for the meticulous review of the manuscript. We also thank Prof. Jusheng Dai, Prof. Shaochun Yang for their insightful reviews and comments.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Nelson, R.A. (2001) Geologic Analysis of Naturally Fractured Reservoirs. 2nd Edition, Gulf Professional Publishing, Boston, 332.
[2]	Florez-Nino, J.M., Aydin, A., Mavko, G., Antonellini, M. and Ayaviri, A. (2005) Fault and Fracture Systems in a Fold and Thrust Belt: An Example from Bolivia. AAPG Bulletin, 89, 471-493. https://doi.org/10.1306/11120404032
[3]	Lonergan, L., Jolly, R.J.H., Rawnsley, K. and Sanderson, D.J. (2007) Fractured Reservoirs: Geological Society (London) Special Publication, 270, 285.
[4]	Smart, K.J., Ferrill, D.A., Morris, A.P. and McGinnis, R.N. (2012) Geomechanical Modeling of Stress and Strain Evolution during Contractional Fault-Related Folding. Tectonophysics, 576-577, 171-196. https://doi.org/10.1016/j.tecto.2012.05.024
[5]	Price, N.J. (1959) Mechanics of Jointing in Rocks. Geological Magazine, 96, 149-167.
[6]	Stearns, D.W. (1968) Certain Aspects of Fracture in Naturally Deformed Rocks. In: Riecker, R.E., Ed., NSF Advanced Science Seminar in Rock Mechanics, Air Force Cambridge Research Laboratories Special Report, Bedford, MA, 97-118.
[7]	Narr, W., Schechter, D.W. and Thompson, L.B. (2006) Naturally Fractured Reservoir Characterization. Society of Petroleum Engineers, Richardson, TX, 112.
[8]	Ackermann, R.V., Schlische, R.W. and Withjack, M.O. (2001) The Geometric and Statistical Evolution of Normal Fault System: An Experimental Study of the Effects of Mechanical Layer Thickness on Scaling Law. Journal of Structural Geology, 23, 1803-1819. https://doi.org/10.1016/S0191-8141(01)00028-1
[9]	Lin, X.B., Chen, H.L., Cheng, X.G., Shen, Z.Y., Yang, S.F. and Xiao, A.C. (2010) Conceptual Models for Fracturing in Fault Related Folds. Journal of Mining Science and Technology, 20, 103-108.
[10]	Zeng, L.B. and Li, Y.G. (2010) Tectonic Fractures in Tight Gas Sandstones of the Upper Triassic Xujiahe Formation in the Western Sichuan Basin, China. Acta Geologica Sinica (English Edition), 84, 1229-1238. https://doi.org/10.1111/j.1755-6724.2010.00293.x
[11]	Ju, W., Hou, G.T. and Hari, K.R. (2013) Mechanics of Mafic Dyke Swarms in the Deccan Large Igneous Province: Palaeostress Field Modeling. Journal of Geodynamics, 66, 79-91. https://doi.org/10.1016/j.jog.2013.02.002
[12]	Ringenbach, J.C., Pinet, N., Stephan, J.F. and Delteil, J. (1993) Structural Variety and Tectonic Evolution of Strike-Slip Basins Related to the Philippine Fault System, Northern Luzon, Philippines. Tectonics, 12, 187-203. https://doi.org/10.1029/92TC01968
[13]	Little, T.A. and Roberts, A.P. (1997) Distribution and Mechanism of Neogene to Present-Day Vertical Axis Rotations, Pacific Australian Plate Boundary Zone, South Island, New Zealand. Journal of Geophysical Research-Solid Earth, 102, 20447-20468.
[14]	Alik, I. Z, Mueller, B. and Schubert, G. (2005) Three-Dimensional Numerical Modeling of Contemporary Mantle Flow and Tectonic Stress Beneath the Earthquake-Prone Southeastern Carpathians Based on Integrated Analysis of Seismic, Heat Flow, and Gravity Data. Physics of the Earth and Planetary Interiors, 149, 81-98. https://doi.org/10.1016/j.pepi.2004.08.012
[15]	Corbett, K., Friedman, M. and Spang, J. (1987) Fracture Development and Mechanical Stratigraphy of Austin Chalk, Texas. AAPG Bulletin, 71, 17-28.
[16]	Hanks, C.L., Lorenz, J., Teufel, L. and Krunhardt, A.P. (1997) Lithologic and Structural Controls on Natural Fracture Distribution and Behaviour within the Lisburne Group, Northeastern Brooks Range and North Slope Subsurface, Alaska. AAPG Bulletin, 81, 1700-1720.
[17]	Wennberg, O.P., Svana, T., Azizzadeh, M., Aqrawi, A.M.M., Brockbank, P., Lyslo, K.B. and Ogilvie, S. (2006) Fracture Intensity vs. Mechanical Stratigraphy in Platform Top Carbonates: The Aquitanian of the Asmari Formation, Khaviz Anticline, Zagros, SW Iran. Petroleum Geoscience, 12, 235-245.
[18]	Zeng, L.B., Wang, H.J., Gong, L. and Liu, B.M. (2010) Impacts of the Tectonic Stress Field on Natural Gas Migration and Accumulation: A Case Study of the Kuqa Depression in the Tarim Basin, China. Marine and Petroleum Geology, 27, 1616-1627. https://doi.org/10.1016/j.marpetgeo.2010.04.010
[19]	Falaknaz, N., Aubertin, M. and Li, L. (2015) Evaluation of the Stress State in Two Adjacent Backfilled Stopes within an Elasto-Plastic Rock Mass. Geotechnical and Geological Engineering, 10, 1-10.
[20]	Pablo, F.S., Pollard, D.D., Allwardt, P.F. and Borja, R.I. (2008) Mechanical Models of Fracture Reactivation and Slip on Bedding Surfaces during Folding of the Asymmetric Anticline at Sheep Mountain, Wyoming. Journal of Structural Geology, 30, 1177-1191. https://doi.org/10.1016/j.jsg.2008.06.002
[21]	Lisle, R.S. (1994) Detection of Zones of Abnormal Strains in Structures Using Gaussian Curvature Analysis. AAPG Bulletin, 78, 1811-1819.
[22]	Hennings, P.H., Olson, J.E. and Thompson, L.B. (2000) Combining Outcrop Data and Three-Dimensional Structural Models to Characterize Fractured Reservoirs. AAPG Bulletin, 84, 830-849.
[23]	Bergbauer, S. (2007) Testing the Predictive Capability of Curvature Analyses. In: Jolley, S.J., Barr, D., Walsh, J.J. and Knipe, R.J., Eds., Structurally Complex Reservoirs, Geological Society (London) Special Publication 292, 185-202.
[24]	Ju, W., Hou, G.T. and Zhang B. (2014) Insights into the Damage Zones in Fault-Bend Folds from Geomechanical Models and Field Data. Tectonophysics, 610, 182-194. https://doi.org/10.1016/j.tecto.2013.11.022
[25]	Masaferro, J.L., Bulnes, M., Poblet, J. and Casson, N. (2003) Kinematic Evolution and Fracture Prediction of the Valle Morado Structure Inferred from 3D Seismic Data, Salta Province, Northwest Argentina. AAPG Bulletin, 87, 1083-1104. https://doi.org/10.1306/02070301102
[26]	Erickson, S.G. and Jamison, W.R. (1995) Viscous-Plastic Finite-Element Models of Fault-Bend Folds. Journal of Structural Geology, 17, 561-573. https://doi.org/10.1016/0191-8141(94)00082-B
[27]	Gudmundsson, A., Simmenes, T.H., Larsen, B. and Philipp, S.L. (2010) Effects of Internal Structure and Local Stress on Fracture Propagation, Deflection, and Arrest in Fault Zones. Journal of Structural Geology, 32, 1643-1655. https://doi.org/10.1016/j.jsg.2009.08.013
[28]	He, J.K., Lu, S.J. and Wang, W.M. (2013) Three-Dimensional Mechanical Modeling of the GPS Velocity Field around the Northeastern Tibetan Plateau and Surrounding Regions. Tectonophysics, 584, 257-266. https://doi.org/10.1016/j.tecto.2012.03.025
[29]	Hou, G.T., Wang, C.C., Li, J.H. and Qian, X.L. (2006) Late Paleoproterozoic Extension and a Paleostress Field Reconstruction of the North China Craton. Tectonophysics, 422, 89-98. https://doi.org/10.1016/j.tecto.2006.05.008
[30]	Ju, W. and Sun, W.F. (2016) Tectonic Fractures in the Lower Cretaceous Xiagou Formation of Qingxi Oilfield, Jiuxi Basin, NW China. Part Two: Numerical Simulation of Tectonic Stress Field and Prediction of Tectonic Fractures. Journal of Petroleum Science and Engineering, 146, 626-636. https://doi.org/10.1016/j.petrol.2016.05.002
[31]	Ostadhassan, M., Zamiran, S., Jabbari, H., Osouli, H., Bubach, B. and Oster, B. (2015) Stability Analysis of Multilateral High Density Pad Wells in the Three Forks Formation. SPE Western Regional Meeting, Society of Petroleum Engineers, Garden Grove, CA.
[32]	Smart, K.J., Wyrick, D.Y. and Ferrill, D.A. (2011) Discrete Element Modeling of Martian Pit Formation in Response to Extensional Fracturing and Dilational Normal Faulting. Journal of Geophysical Research: Planets, 116, 383-392. https://doi.org/10.1029/2010JE003742
[33]	Yin, A. (1994) Mechanics of Monoclinal Systems in the Colorado Plateau during the Laramide Orogeny. Journal of Geophysical Research, 99, 22043-22058. https://doi.org/10.1029/94JB01408
[34]	Rezaei, M., Hossaini, M.F. and Majdi, A. (2015) Determination of Longwall Mining-Induced Stress Using the Strain Energy Method. Rock Mechanics and Rock Engineering, 48, 2421-2433. https://doi.org/10.1007/s00603-014-0704-8
[35]	Fakhimi, A. (2004) Application of Slightly Overlapped Circular Particles Assembly in Numerical Simulation of Rocks with High Friction Angle. Engineering Geology, 74, 129-138. https://doi.org/10.1016/j.enggeo.2004.03.006
[36]	Fakhimi, A. (2009) A Hybrid Discrete-Finite Element Model for Numerical Simulation of Geomaterials. Computers and Geotechnics, 36, 386-395. https://doi.org/10.1016/j.compgeo.2008.05.004
[37]	Potyondy, D.O. (2007) Simulating Stress Corrosion with a Bonded-Particle Model for Rock. International Journal of Rock Mechanics and Mining Sciences, 44, 677-691. https://doi.org/10.1016/j.ijrmms.2006.10.002
[38]	Tarokh, A. and Fakhimi, A. (2014) Discrete Element Simulation of the Effect of Particle Size on the Size of Fracture Process Zone in Quasi-Brittle Materials. Computers and Geotechnics, 62, 51-60. https://doi.org/10.1016/j.compgeo.2014.07.002
[39]	Igor, A.B. (2015) Regularization of Non-Convex Strain Energy Function for Non-Monotonic Stress-Strain Relation in the Hencky Elastic-Plastic Model. Acta Mechanica, 226, 2681-2691. https://doi.org/10.1007/s00707-015-1349-8
[40]	Price, N.J. (1966) Fault and Joint Development in Brittle and Semi-Brittle Rock. Pergamon Press, Oxford, 133-164.
[41]	Wen, S.P. and Li, D.T. (1996) Simulation Technology on Structural Joints of Reservoirs. Journal of China University of Petroleum (Edition of Natural Science), 20, 17-24.
[42]	Song, H.Z. (1999) An Attempt of Quantitative Prediction of Natural Crack on Brittle Rock Reservoir. International Journal of Geomechanics, 5, 76-84. (In Chinese with English abstract)
[43]	Tan, C.X. and Wang, L.J. (1999) An Approach to the Application of 3-D Tectonic Stress Field Numerical Simulation in Structural Fissure Analysis of the Oil-Gas-Bearing Basin. Acta Geoscientia Sinica, 20, 392-394.
[44]	Dai, J.S., Wang, B.F. and Ma, Z.R. (2007) Research on Cracking Principles of Brittle Low-Permeability Sands. Xinjiang Petroleum Geology, 28, 393-395. (In Chinese with English abstract).
[45]	Dai, J.S., Feng, J.W., Li, M., et al. (2011) Discussion on the Extension Law of Structural Fracture in Sand-Mud Interbed Formation. Earth Science Frontiers, 18, 277-283. https://doi.org/10.1016/j.jsg.2016.03.007
[46]	Guo, P., Yao, L.H. and Ren, D.S. (2016) Simulation of Three-Dimensional Tectonic Stress Fields and Quantitative Prediction of Tectonic Fracture within the Damintun Depression, Liaohe Basin, Northeast China. Journal of Structural Geology, 86, 211-223.
[47]	Zhao, W.T., Hou, G.T. and Sun, X.W. (2013) Influence of Layer Thickness and Lithology on the Fracture Growth of Clastic Rock in East Kuqa. Geotectonica et Metallogenia, 37, 603-610.
[48]	Zhang, Z.P., Wang, Q.C. and Wang, Y. (2006) Brittle Structure Sequence in the Kuqa Depression and Its Implications to the Tectonic Paleostress. Earth Science-Journal of China University of Geosciences, 31, 310-316.
[49]	Jia, C.Z. (2004) Characteristics of the Mesozoic and Cenozoic Structures and Petroleum Occurrence in the Tarim Basin. Petroleum Industry Press, Beijing, 229.
[50]	Zheng, C.F., Hou G.T., Lu Y., et al. (2016) An Analysis of Cenozoic Tectonic Stress Fields in the Kuqa Depression. Geological Bulletin of China, 35, 130-138.
[51]	Zeng, L.B. and Wang, G.W. (2005) Distribution of Earth Stress in Kuqa Thrust Belt, Tarim Basin. Petroleum Exploration and Development, 32, 59-60. (In Chinese)
[52]	Zhang, R.H., Zhang, H.L. and Shou, J.F. (2008) Geological Analysis on Reservoir Mechanism of the Lower Cretaceous Bashijiqike Formation in Dabei Area of the Kuqa Depression. Chinese Journal of Geology, 43, 507-517. (In Chinese with English abstract)
[53]	Li, F., Jiang, Z.X. and Li, Z. (2015) Genetic Type and Hydrocarbon Accumulation Mechanism of Dina 2 Gas Reservoir in Kuqa Depression. Journal of Central South University (Science and Technology), 46, 1346-1352.
[54]	Dai, J.S. (2006) Structural Geology (with Tectonics). Petroleum Industry Press, Beijing, 1, 110-112.
[55]	Van Golf-Racht, T.D. (1982) Fundamentals of Jointed Reservoir Engineering. Elsevier Scientific, New York.
[56]	Hegelson, D.E. and Aydin, A. (1991) Characteristics of Joints Propagation across Layers Interfaces in Sedimentary Rocks. Journal of Structural Geology, 13, 897-911. https://doi.org/10.1016/0191-8141(91)90085-W
[57]	Bellahsen, N., Fiore, P. and Pollard, D.D. (2006) The Role of Fractures in the Structural Interpretation of Sheep Mountain Anticline, Wyoming. Journal of Structural Geology, 28, 850-867. https://doi.org/10.1016/j.jsg.2006.01.013
[58]	Feng, J.W., Dai, J.S. and Ma, Z.R. (2011) The Theoretical Model between Fracture Parameters and Stress Field of Low-Permeability Sandstones. Acta Petrolei Sinica, 32, 664-671. (In Chinese)
[59]	Beer, F.P., Johnston, E.R., Dewolf, J.T. and Mazurek, D.F. (2012) Mechanics of Materials. 6th Edition, McGraw-Hill Companies, New York, 758.
[60]	Hoek, E. (1990) Estimating Mohre Coulomb Friction and Cohesion Values from the Hoeke Brown Failure Criterion. International Journal of Rock Mechanics and Mining Sciences, 27, 227-229. https://doi.org/10.1016/0148-9062(90)94333-O
[61]	Barton, N. and Choubey, V. (1977) The Shear Strength of Rock and Rock Joints in Theory and Practice. Rock Mechanics, 10, 1-54. https://doi.org/10.1007/BF01261801
[62]	Jaeger, J.C. (1960) Shear Failure of Transversely Isotropic Rock. Geology Magazine, 97, 65-72. https://doi.org/10.1017/S0016756800061100
[63]	Jaeger, J.C., Cook, N.G.W. and Zimmerman, R.W. (2007) Fundamentals of Rock Mechanics.4th Edition, Blackwell Publishing, London.
[64]	Griffith, A.A. (1921) The Phenomena of Rupture and Flow in Solids. Philosophical Transactions of the Royal Society London, 221, 163-197. https://doi.org/10.1098/rsta.1921.0006
[65]	Fossen, H. (2010) Structural Geology. Cambridge University Press, Cambridge, 463. https://doi.org/10.1017/CBO9780511777806
[66]	Ju, W., Sun, W.F. and Hou, G.T. (2015) Insights into the Tectonic Fractures in the Yanchang Formation Interbedded Sandstone-Mudstone of the Ordos Basin Based on Core Data and Geomechanical Models. Acta Geologica Sinica (English Edition), 89, 1801-1812.
[67]	Fischer, K. and Henk, A. (2013) A Workflow for Building and Calibrating 3-D Geomechanical Models—A Case Study for a Gas Reservoir in the North German Basin. Solid Earth, 4, 347-355.
[68]	Ding, W.L., Fan, T.L., Yu, B.S., Huang, X.B. and Liu, C. (2012) Ordovician Carbonate Reservoir Fracture Characteristics and Fracture Distribution Forecasting in the Tazhong Area of Tarim Basin, Northwest China. Journal of Petroleum Science and Engineering, 86-87, 62-70. https://doi.org/10.1016/j.petrol.2012.03.006
[69]	Jiu, K., Ding, W.L., Huang, W.H., You, S.G., Zhang, Y.Q. and Zeng, W.T. (2013) Simulation of Paleotectonic Stress Field within Paleogene Shale Reservoirs and Prediction of Favorable Zones for Fracture Development within the Zhanhua Depression, Bohai Bay Basin, East China. Journal of Petroleum Science and Engineering, 110, 119-131. https://doi.org/10.1016/j.petrol.2013.09.002
[70]	Zhou, C.M., Zhang, Z.J., Wang, S.W. and Li, X.F. (2009) Exploration of Numerical Simulation on Paleo-Tectonic Stress Field. Procedia Earth and Planetary Science, 1, 875-881. https://doi.org/10.1016/j.proeps.2009.09.136
[71]	Joshi, S., Prashant, A., Deb, A. and Jain, S.K. (2011) Analysis of Buried Pipelines Subjected to Reverse Fault Motion. Soil Dynamics and Earthquake Engineering, 31, 930-940. https://doi.org/10.1016/j.soildyn.2011.02.003
[72]	Fischer, M.P. and Woodward, N.B. (1992) The Geometric Evolution of Foreland Thrust Systems. In: Mc Clay, K.R., Ed., Thrust Tectonics, Chapman and Hall, London, 181-189. https://doi.org/10.1007/978-94-011-3066-0_16
[73]	Fisher, D.M. and Anastasio, D.J. (1994) Kinematic Analysis of a Large-Scale Leading Edge Fold, Lost River Range, Idaho. Journal of Structural Geology, 16, 337-354. https://doi.org/10.1016/0191-8141(94)90039-6
[74]	Zhu, G.Y., Yang, H.J. and Zhang, B. (2012) The Geological Feature and Origin of Dina 2 Large Gas Field in Kuqa Depression, Tarim Basin. Acta Petrologica Sinica, 28, 2479-2492.
[75]	Teeuw, D. (1971) Prediction of Formation Compaction from Laboratory Compressibility Data (SPH Paper 2973). Society of Petroleum Engineers Journal, 11, 263-271. https://doi.org/10.2118/2973-PA
[76]	Wang, P. (2001) Tectonic Mechanic Principles of Oil-bearing Basins. Petroleum Industry Press, Beijing, 156. (In Chinese)
[77]	Liu, H.T. and Zeng, L.B. (2004) Expression of the Himalayan Orogeny in the Kuqa Depression, Tarim Basin—Evidence from the Rock Acoustic Emission Experiment. Geological Bulletin of China, 23, 676-679. (In Chinese)
[78]	Liu, C., Huang, X.B., Fan, T.L., Wang, Z.X. and Zeng, Q.B. (2008) The Simulation of Present Tectonic Stress Field and the Prediction of Tectonic Fractures of Ordovician in Tazhong Area, Tarim Basin. Xinjiang Petroleum Geology, 29, 475-477. (In Chinese with English abstract)
[79]	Zeng, L.B., Qi, J.F. and Wang, Y.X. (2007) Origin Type of tectonic Fractures and Geological Conditions in Low-permeability Reservoirs. Acta Pctrolci Sinica, 28, 52-56.
[80]	Mckinnon, S.D. and Barra, I.G. (1998) Fracture Initiation, Growth and Effect on Stress Field: A Numerical Investigation. Journal of Structural Geology, 20, 1663-1672. https://doi.org/10.1016/S0191-8141(98)00080-7
[81]	Tuckwell, G.W., Lonergan, L. and Jolly, R.J.H. (2003) The Control of Stress History and Flaw Distribution on the Evolution of Polygonal Fracture Networks. Journal of Structural Geology, 25, 1241-1250. https://doi.org/10.1016/S0191-8141(02)00165-7
[82]	Ding, W.L., Li, C., Li, C.Y., Xu, C.C., Jiu, K., Zeng, W.T. and Wu, L.M. (2012) Fracture Development in Shale and Its Relationship to Gas Accumulation. Geoscience Frontiers, 3, 97-105. https://doi.org/10.1016/j.gsf.2011.10.001
[83]	Forest, S., Barbe, F. and Cailletaud, G. (2000) Cosserat Modelling of Size Effects in the Mechanical Behaviour of Polycrystals and Multi-Phase Materials. International Journal of Solids and Structures, 37, 7105-7126. https://doi.org/10.1016/S0020-7683(99)00330-3
[84]	Forest, S., Pradel, F. and Sab, K. (2001) Asymptotic Analysis of Heterogeneous Cosserat Media. International Journal of Solids and Structures, 38, 4585-4608. https://doi.org/10.1016/S0020-7683(00)00295-X
[85]	Suppe, J. and Medwedeff, D.A. (1990) Geometry and Kinematics of Fault Propagation Folding. Eclogae Geologicae Helvetiae, 83, 409-454.
[86]	Shaw, J.H., Bilotte, F. and Brennan, P.A. (1999) Patterns of Imbricate Thrusting. GSA Bulletin, 111, 1140-1154. https://doi.org/10.1130/0016-7606(1999)111<1140:POIT>2.3.CO;2
[87]	Salvini, F., Storti, F. and Mc Clay, K. (2001) Self-Determining Numerical Modeling of Compressional Fault-Bend Folding. The Geological Society of America, 29, 839-842. https://doi.org/10.1130/0091-7613(2001)029<0839:SDNMOC>2.0.CO;2
[88]	Suppe, J. (1983) Geometry and Kinematics of Fault-Bend Folding. American Journal of Science, 283, 684-721. https://doi.org/10.2475/ajs.283.7.684
[89]	Suppe, J. and Chang, Y.L. (1983) Kink Method Applied to Structural Interpretation of Seismic Sections, Western Taiwan. Petroleum Geology of Taiwan, 19, 29-49.
[90]	Erickson, S.G., Strayer, L.M. and Suppe, J. (2001) Initiation and Reactivation of Faults during Movement over Thrust-Fault Ramp: Numerical Mechanical Models. Journal of Structural Geology, 23, 11-23. https://doi.org/10.1016/S0191-8141(00)00074-2
[91]	Kuenen, P.H. and de Sitter, L.U. (1938) Experimental Investigation into the Mechanisms of Folding. Leidshe Geologische Mededeelingen, 10, 217-240.
[92]	Fossen, H. (2010) Structural Geology. Cambridge University Press, Cambridge.