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

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.


Introduction
Natural fractures influence the performance of many reservoirs around the world, including carbonate reservoir, deep tight sand reservoir and low-permeability 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 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 Journal of Geoscience and Environment Protection reservoirs.

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 Qiulitage 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 2 , 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).   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 −3 μm 2 .However, the fracture permeability can reach (1.00 -10.00) × 10 −3 μm 2 .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 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 Journal of Geoscience and Environment Protection 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 −1 ) and volume density (m 2 /m 3 ), 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 −1 to 1.5 m −1 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 −1 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)).

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 σ 1 , σ 2 and σ 3 , the Journal of Geoscience and Environment Protection strain-energy density at a given point can be expressed as [59]: ( ) where ϖ is strain energy density (J/m 3 ); and 1 ε , 2 ε , 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.
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 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( e ϖ ), a formula keep to the principle of energy conservation will be given by: where vf D is fracture volume density per unit volume (m 2 /m 3 ); f ϖ is strain energy density for new fractures (J/m 3 ); V is characterization of unit cell volume (m 3 ); S f is surface area of new fractures (m 2 ); J is the required energy per unit area for fractures (J/m 2 ), i.e., fracture surface energy (here the energy is different from and far less than theoretical value of molecular dissociation); e ϖ is the necessary elastic strain energy to be overcome for new fractures (J/m 3 ); 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: σ ≥ ), the Mohr-Coulomb criterion will be se- lected, that is The Coulomb-Mohr Criterion suggests that a shear fracture only forms if the internal strength or cohesion of the rock (C 0 ) 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 ( ) where ϕ is the internal friction angle (˚); θ is the angle between normal to the newly formed fracture plane and the maximum principal stress (˚); 1 σ is the peak stress (MPa); 2 σ is the intermediate principal stress (MPa); 3 σ is the middle principal stress, (MPa); p σ is the rupture stress under action of 3 σ , different from the maximum principal stress; t σ is the tensile strength, (MPa); 0 J is fracture surface energy with no confining pressure or under unaxial com- pressive stress (J/m 2 ), J ∆ is the additional surface energy caused by confining pressure 3 σ (J/m 2 ); E is the elasticity modulus with no confining pressure (GPa); b is the fracture aperture (i.e.paleo-aperture, here for reference only), m; lf D is the fracture linear density (1/m); 3 ε is the tensile strain under current state of stress, dimensionless parameter; 0 ε is the maximum tensile strain, di- mensionless parameter, corresponding to tensile strain when crack beginning to form; 0 E is the proportionality coefficient related to lithology; and L 1 , L 2 , L 3 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 3 0 σ ≤ and ( ) ( ) 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 1 , L 2 , L 3 , thus the σ 1 corresponding to side L 1 , the σ 2 corresponding to side L 2 , the σ 3 corresponding to the side L 3 ; (b) Transection in REV perpendicular to σ 2 , namely (σ 1 -σ 3 ) plane, here, θ is the acute angle between fracture surface and maximum principal stress σ 1 .
When ( ) If the failure criterion is reached, the relationships between fracture density, aperture and strain energy density, stress-strain are expressed as ( ) 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 , and n is projected to x-o-z plane with angle between projec- tion line and negative direction of z-axis, through the angle conversion ( ) , the angle of strike α will be calculated.
If 0 90 From the geological point of view, fracture dip dip α between x-o-z plane and fracture plane also can be defined as angle between the 0 lx my nz ), which is expressed as

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].

Geometry and Boundary Conditions
The purpose of this step is to set up several suites of geomechanical models and 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.
where z σ is the vertical stress, 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/minimum principal stress ratio ( 1 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 + h σ ) MPa is applied in the NNW direction (i.e.x-direction)and a minimum principal stress of (66.4 + 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 Journal of Geoscience and Environment Protection 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.

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 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).

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.

Calculation of Solution
After establishing several suites of mathematical models, we used the Ansys software to simulate the paleotectonic stresses fields during the Himalayan Pe-

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.vf D , lf D and dip α are used to ex- press 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 e ϖ in theory, hence, we assume that e ϖ is strain energy density at 0.85 c σ σ = .Recalling the Hooke's Law the e ϖ and ϖ can be ob- tained [60], then f w is calculated.According to [58], for low-permeability sandstone the corresponding average 0 J (fracture surface energy under unaxi- al stress state) is 1087.35J/m 2 , and the average 0 E is 112.6.Secondly, putting 0 J and 0 E into the Equation ( 5) and combining with formulas of vf D , lf D , b and J , the fracture volume density vf D , 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 vf D , lf D and dip α are obtained using the same method, which is built into the An- sys platform to calculated spatial distribution of fracture parameters (Figure 7 and Figure 8).

Results and Discussions
To gain insight into the stress field and damage zone during fault-related folding, 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).

Effect of Stresses on Fracture Development and Distribution
Commonly, the development of tectonic fractures is generally controlled by  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 ( ) where y is the simulated fracture density in model (m −1 ), x is the differential stress (i.e., σ 1 − σ 3 ), 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.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 ( ) where y is the simulated fracture density in model (m −1 ), 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 −1 ) 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 where y is the predicted fracture linear density in model (m −1 ), 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 −1 ), 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)).
( ) where y is the predicted fracture linear density in model (m −1 ), 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 −1 ).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)).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.

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.

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 Journal of Geoscience and Environment Protection 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.

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 .
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 1 j: Jidike Formation; N 1 k: Kangcun Formation; N 2 k-Q 1 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.

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

Figure 3 .
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].

Figure 5 .
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.
riod 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 (Figure6).In the Suweiyi Group and Kumugeliemu Group of D gas field, the σ 1 , indicative of compression, range from 140 MPa to220 MPa.The σ 2 (i.e.vertical stress) indicative of compression, are generally between 56 MPa and 123 MPa, and the σ 3 , indicative of both compression and tension, ranges from −30 MPa to 15 MPa (Figure6(b)).The differential stress by σ 3 subtracting σ 1 , indicative of possible damage zone, mainly ranges from 0 to 85 MPa (Figure6(a)).The pre-existing fault tip show a higher differential stress, because with the regional contraction and the increasing displacement

Figure 6 .
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.

Figure 7 .
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.

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

Figure 9 .
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 . 6 . 2 .
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.
(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 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 (Figure10(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).

Table 1 .
Rock mechanical parameters in the EI and EII of D gas field.
).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 Table1.12 Journal of Geoscience and Environment Protection