An improved geomechanical model for the prediction of fracture generation and distribution in brittle reservoirs

It is generally difficult to predict fractures of low-permeability reservoirs under high confining pressures by data statistical method and simplified strain energy density method. In order to establish a series of geomechanical models for the prediction of multi-scale fractures in brittle tight sandstones, firstly, through a series of rock mechanics experiments and CT scanning, we determined 0.85 σc as the key thresholds for mass release of elastic strain energy and bursting of micro-fractures. A correlation between fracture volume density and strain energy density under uniaxial stress state was developed based on the Theory of Geomechanics. Then using the combined Mohr-Coulomb criterion and Griffith’s criterion and considering the effect of filling degree in fractures, we continued to modify and deduce the mechanical models of fracture parameters under complex stress states. Finally, all the geomechanical equations were loaded into the finite element (FE) platform to quantitatively simulate the present-day 3-D distributions of fracture density, aperture, porosity, permeability and occurrence based on paleostructure restoration of the Keshen anticline. Its predictions agreed well with in-situ core observations and formation micro-imaging (FMI) interpretations. The prediction results of permeability were basically consistent with the unobstructed flow distributions before and after the reservoir reformation.


Introduction
As an important unconventional resource and clean-burning fuel, tight sandstone gas has been widely distributed in many countries of the world, such as the U.S., Canada, Australia, Mexico, Russia, and China [1], [2], [3]. In general, tight sandstone reservoirs distinctively differ from conventional reservoirs because of their great burial depth, strong diagenesis, strong heterogeneity, abnormal overpressure, low porosity and low permeability, and developed fractures [4], [5], [6], [7]. In these tight low-permeability sandstones, the majority of reservoir spaces and seepage channels for hydrocarbon are primarily provided by widely distributed natural fractures, especially tectonic fractures, which can significantly and effectively improve the permeability of a reservoir and enhance hydrocarbon delivery to wellbores [8], [9], [10], [11], [12], [13]. Therefore, understanding and interpreting where and when tectonic fractures develop within a geological structure, along with their orientation, intensity and porosity, are important in both exploration and production of tight sandstone reservoirs. However, how to effectively predict fractures is still a worldwide challenge. Three-dimensional (3D) characterization and modeling of subsurface tectonic fractures in deep tight gas reservoirs with complex tectonism and diagenesis presents more difficulty [14], [15], [16], [17]. For exploration and production of fractured reservoirs, fracture prediction or modeling is commonly based on geometric and/or kinematic models, such as analyses of fault-related folds and fold curvature [18], [19], [20], [21], seismic techniques or logging methods [22]. These approaches, aiming to acquire the attributes of inter-well fracture networks, are useful in that they stick to the actual measured data and relate deformation to corresponding structural position [23]. Nevertheless, they are specifically tied to ideal geometric models or simplified assumptions that may not completely reflect the multi-phase deformation behaviors and mechanical properties of underground rocks [8], [14], [23].
Generally, tectonic stress field is the most important factor in controlling development and distribution of tectonic fractures in reservoirs [24], [25], [26], [21], [27]. Therefore, one efficient geomechanical modeling strategy used in recent exploration and development of brittle reservoirs is studying concentrations and changes in paleotectonic stress field so as to determine critical process involved in fracture development and combine various rupture criterions to predict favorable zones of fractures [20], [28], [29], [23], [30], [31].
Since the 1960's, many studies have been published on the mechanisms of fracture-generating structural movement, including rock failure criterion, indicator of comprehensive rupture rate, and strain energy density. Price [32] proposed that fracture intensity was positively correlated with the elastic strain energy in rocks based on laboratory experiments. Song [33] successfully applied tensional failure criterion and shear failure criterion to calculate a rock fracturing index that predicted fracture development zones and dominant orientation. Tan and Wang [34] as well as Zhou et al. [35] put forth a similar quasi-binary method that combined failure degree value and energy degree value to quantitatively characterize fracture density. The method was based on Maximum Strain Energy Density Theory that rocks with highstrain energy possessed more fractures than rocks with low strain energy. According to the theory of strain energy density, Dai et al [36] attempted to establish a series of formulas relating stress-strain and fracture volume density and linear density under paleostress field based on comprehensive rock mechanical experiments. Olson and Laubach [14] presented a mode of natural fracture analysis that incorporates fracture mechanics and diagenetic processes to predict fracture network geometry and fracture aperture distribution and preservation in the absence of subsurface rock samples. However, almost all previous geomechanical modeling approaches for fracture prediction have been based on distributions of stress-strain, strain energy density, or rupture rate, but far less directly with further quantitative mechanical models between fracture parameters and strain energy density during different critical tectonic deformation stages. Additionally, there has been little research on the characterization of multi-scale fractures. Based on technologies of FIB-SEM tomography and X-ray μ-CT, Huang et al. [37] and Li et al. [38] conducted multi-scale quantitative characterization of 3-D porefracture networks in bituminous and anthracite rocks and meso-scale fracture modelling. From the perspective of evolution, Ghamgosar and Erarslan [39] used CT scan technique to investigate and evaluate the micro-fracturing samples in fracture process zones under cyclic and static loadings. The following conclusions are drawn from this study: intergranular cracks are formed due to particle breakage under cyclic loading compared with smooth and bright cracks along cleavage planes under static loading, and the macro-scale main crack causing failure is seen in cement without any dust or debris material under monotonic loading.
In this study, we first conducted rock mechanical tests and advanced X-ray CT scanning to accurately determine the key threshold values for multi-scale fracture generation and development in brittle tight sandstones. A quantitative relationship between strain energy density and fracture volume density was established. Applying the classic geological mechanics theories, i.e., the Maximum Tension Theory and Maximum Strain Energy Density Theory, we continued to deduce the associated mechanical models of fracture parameters (fracture linear density, aperture, porosity and permeability) under the complex stress states. Then, on the basis of tectonic evolutionary restoration of Keshen anticline in the Kuqa Depression, we used finite element method (FEM) to better simulate 3-D paleostress field and current stress field during the entire deformation process of faults and folds. Combining with composite fracture criterion and considering the influences of the current stress and filling degree on the aperture of a fracture, we incorporated these mechanical models into a numerical simulator ANSYS to predict the distributions of tectonic fractures in Keshen tight gas field, which were effectively verified by the measured fracture density from drilled cores and imaging logs.

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 (Fig 1). Structures that develop widely in the Kuqa Depression during the Cenozoic Period are dominated by thrust faults and related folds. Laterally, the Kuqa Depression can be divided into three structural belts and two sags, which are, from north to south, the northern monocline belt, the Kelasu structural belt, the Baicheng sag, the Kuqa depression, and the Qiulitage structural belt (Fig 1a) [40], [29]. The Kuqa Depression is recognized as one of the major Cenozoic depocenters along the margin of the Tarim Basin because it has experienced a complex evolution since the Late Cretaceous with the northward movement of the Indian sub-continent and the southward thrusting of the South Tianshan [41], [42], [43]. The Keshen gas field is situated in the footwall of Kelasu tectonic belt, south of Kela fault, and displays a narrow anticline with a structure amplitude of less than 500 m. The top of the Keshen gas field is cross-cut by several faults striking E-W to become complicated, acquiring a dip in the southern wing ranging from 16˚to 20˚and in the northern wing ranging from 19˚to 23˚ (Fig 1b).
Fracture intervals in the Keshen gas field consist of Cretaceous delta front rocks (the Bashijiqike Formation, K 1 bs), which are dominated by fine sandstone, siltstone, and mudstone with limited sandy conglomerate. These beds reach a total thickness of up to 290m within the study area and are divided into three members according to lithological cycles and interbeds (Fig 1c), [29], [7]. The porosity of the Bashijiqike reservoir, as determined by core tests, ranges from 2-7%. Permeability of the reservoir lies in the range of (0.05-0.50) ×10 −3 μm 2 , while fracture permeability can reach (1.00-10.00) ×10 −3 μm 2 . In summary, the above evidence indicates that the Keshen gas field belongs to a typical ultra-deep low-porosity and low-permeability tight sandstone reservoir.
Based on quantitative observation and statistical description in drilling cores and formation micro-imaging (FMI), simultaneously using field outcrops for reference, the fractures most frequently encountered in the reservoir of the area are planar discontinuities that are sub-perpendicular to the bedding and can be divided into three basic types, namely tension fracture, tensoshear fracture, and shear fracture, respectively accounting for 79.5%, 18.7%, 1.8% of the total fracture volume (Fig 2a). The shear fracture has a straight rupture plane, a longer extended distance than the other two types, and in most cases can cut through rock grains. The tensional fracture has a dendritic structure, a relatively shorter distance, and frequently bypasses rock grains, wherever they are observable by drill core and FMI. Observations from cores and analysis from FMI 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) (Fig 2b). Due to the effect of stress unloading after exposure to ground, the physical parameters of fracture network will change slightly relative to underground conditions, such as fracture aperture and porosity. With fine integration of core observation and well-logging data, the detailed parameters of tectonic fractures can be obtained. It is shown that most fracture in K 1 bs are unfilled and slight-filled, which accounts for more than half of the total fracture number (approximately 72%) (Fig 2c). In contrast, only 17% of total fractures are most-filled and full-filled, with calcite and dolomite as the main fillings, followed by small amounts of mud and carbon. The statistical result shows that the fracture aperture in this study area has a bimodal distribution with the main peak range of 0-0.2 mm and 0.4-0.6 mm, and the former dominated in addition (Fig 2d).
The fracture density generally contains linear density (1/m) and volume density (m 2 /m 3 ), and the former is an important parameter to illuminate the development and distribution characteristics of fractures [21]. Observations made in the drilling cores of the tight sandstones 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 (Fig 2e). Structures in two limbs are mainly composed of high-angle conjugated fractures with linear density reaching up to 2.33/m, that is characteristic of systematic joints as defined by Engelder [44]. On the whole, the fracture density has higher value in limbs than that of hinge area and has higher value in southern limb than that of northern limb. However, near the fault zone, there has the highest fracture density value (9.36/m), which can be called the fracture damage zone (Fig 2e and 2f).
Based on the above statistics, the fractures present in the Bashijiqike Formation can be subdivided based on orientation into four distinct, mutually abutting fracture sets oriented NW-SE (set I), NNE-SSW (Set II), NE-SW (set III) and nearly EW (set IV), among which the former two sets are dominant (Fig 2f).

Core samples
Lithologically, K 1 bs member in the Keshen gas field were dominated by lithic feldspar sandstone, followed by feldspar lithic sandstone and little feldspathic sandstone. The sandstones are medium grained (0.5-0.25 mm) and fine grained (0.25-0.05 mm), medium cemented, medium to well sorted, containing much cement between the particles. The mineralogical composition of sandstone is 50% quartz, 12% feldspar, 5% clay, and 33% debris, respectively. This composition indicates that the K 1 bs belonged to medium hard to hard brittle sandstone with strong elastic-brittle characteristics. Forty-eight core plugs were drilled from K 1 bs along the direction of the bedding plane. All plugs were processed according to ISRM-suggested methods with a core diameter of 25 mm and length of 50 mm, ends ground flat to 0.01 mm with vertical deviation of less than 0.25˚.

Rock mechanical measurements
As listed in Tables 1 and 2, the plugs were subjected to a series of rock mechanical experiments that included uniaxial compression tests, triaxial compression tests, direct shear tests, and splitting tests. A Uniaxial/Triaxial Mechanics Instrument (TAWA-2000, USA) was used to simulate the failure process of plugs. This is a hybrid loading instrument with a maximum axial pressure of 2,000 kN and a maximum confining pressure of 140 MPa used for uniaxial compression and triaxial compression tests, with a loading rate of 0.01 mm/s. Over the course of loading, the confining pressure was manually set to seven stages at an interval of 5 MPa (0, 5, 10, 15, 20, 25 and 30 MPa) as the axial pressure increased up to peak stress. At all times, axial loads exceeded confining pressure by no more than one-tenth of the rock uniaxial compression strength (UCS). Before the CT scanning tests, uniaxial compression tests of several spare samples were conducted to make estimation of uniaxial compression strength or peak strength (σ c ). In this way, the complete stress-strain curves and mechanical parameters of sandstone were acquired, including the shear strength (τ s ), cohesive strength (C 0 ), and friction angle (φ) (Tables 1 and 2) (S1 Appendix).

CT Scanning experiments
The CT (Computer Tomography) scanning experiments were conducted at a Micro XCT-400 CT scanner at the Research Center of Oil & Gas Flow in Reservoir, China University of Petroleum. The mechanical deformation process of six core samples under classical uniaxial tests was monitored real time by the CT scanner (Fig 3). As shown in Fig 3a, the core compression device was made of beryllium that has low attenuation under X-ray. During the multi-stage loading test, the sample was loaded to various percentages of the peak stress (i.e., peak strength σ c ) acquired from the preliminary uniaxial tests, and the crucial load increment between two cycles was determined by each rapid generation of micro-fractures during uniaxial compression processes. Seven scans were performed in the loading process, i.e., at the initial loading, 25%, 50%, 65%, 85%, 100%, and 110% of stress peak, and post-peak, respectively, with a CT slice thickness of 2 mm, 25 slices per scan (S2 Appendix). For a grey-level CT image, the bright color denotes the low-density part such as fractured zone or pores, and the dark color denotes the high-density rock matrix, as shown in Fig 3b. The scanned 2-D images were reconstructed into 3-D geometry and then imported to Avizo (ThermoFisher Scientific, USA) for further image analysis (Fig 3c). In this study, it is important to distinguish induced different-sized fractures and natural vugs from CT images. As irregular morphology (e.g., nearly round, calabash-shaped or curviplanar) and various sizes after image processing, the micro-scale porethroats were calculated to be 5-90 μm in diameter and primarily distributed in isolated points, and locally in strips (Fig 3c). In contrast, the original micro-fractures developed in intact rock samples were measured to be 5-15 μm wide and 300 μm long, mainly took on three kinds of shapes: short intragranular fracture, straight rectangular fracture, and crooked grain-edged fracture with sharp boundaries between grains. Then, the fractures and vugs in processed CT images at initial loading stage could be easily classified with visual observation, and the original vug volume was necessarily subtracted from subsequent total damage or fracture volume during each deformation stage. Finally, the strain energy density and different-scale fracture density of the low-permeability tight sandstone samples corresponding to each crucial fracture growth stage were calculated, described in detail, and translated into a series of stress-strain curves for the purpose of identifying relationships between stress and fracture parameters (Table 3).

Evolution of fractures under uniaxial compression
According to the statistical results in Fig 4a, at the beginning the pre-existing micro-fractures in rocks began to close with compression, allowing rocks to become dense, fracture aperture to reduce significantly. In this stage, most of the original micro-fracture length was less than 0.3 cm. When the applied stress was 25% and 50% of the peak value, a few small new microfractures with random orientations began to initiate, among which the fractures with length ranging from 0.2 cm to 0.3 cm had relatively higher growth rate than other fractures of different scales. When the applied stress was increased to 65% of the peak value, growth rate difference gradually appeared among all different-scale fractures, such as the number of small-scale micro-fractures (length<0.3 cm) still steadily increased, whereas that of relatively middle-scale fractures (length ranging from 0.4-0.6 cm) rapidly increased at this stress stage, and that of fractures with length ranging from 0.3 to 0.4 cm decreased instead. Undoubtedly, new generating fractures in this stage tended to grow preferentially in the stress direction. When the applied stress was increased to 85% of the peak value even after this point, there was a rapid increase in the number of all scales of fractures (except the length<1 cm), which were of strong orientation and approximately equal in length. At the same time, after this critical stress value, large-scale fractures (length exceeding 1cm) began to appear, and their number gradually increased. This indicated that fractures tended to develop along anisotropy when the angle between pre-existing micro-fractures and stress direction was low (<0˚) (Fig 4b). However, during these stages, apertures of most of the fractures were still less than 100 μm, and they belonged to micro-scale fractures. When the applied stress was increased to the peak value or even after this point, it was interesting that only the number of large-scale fractures (length even reaching 5 cm and aperture even exceeding 100 μm) suddenly increased, in contrast, the other levels of fractures were almost invariable in number. It could be seen clearly for Figs 3 and 4 that these seemingly randomly distributed micro-fractures formed in previous stages began to coalesce with one another along certain predominant direction until the formation of macro-fractures, which further propagated along the stress direction approximatively (Fig 4b).
From the above evidence, it could be concluded that a large number of new fractures appeared at the compression stress levels of approximately 0.85 σ c (i.e., point B in Fig 4c and  4d), which could be defined as the threshold values for micro-fracture inception and development during stress loading. Zhao [45] conducted uniaxial compression experiments on Mechanics Instrument with sandstone samples under various stress conditions and find that micro-fractures increases as stress boosted, especially when stress reaches or exceeds 0.85 σ c or 0.88 σ c . At this stress level, the number of fractures exploded instantaneously and the number of middle-sized fractures (0.6 cm>length>0.2 cm) increased faster than that of small-scale and large-scale fractures. At the same time, the number of fractures with small angles with respect to principal stress (�30˚) increased faster than that of fractures with larger angles (>30˚) (Fig 4b). While the axial stress reached peak strength at Point C, the overall macro-rupture occurred along with a significant decrease in stress and a significant increase in strain (Figs 4c and 4d). When the stress exceeded the peak strength of rock sample, it began to fall until the residual strength remained, thus the structural fault was created, accompanied by a large amount of frictional energy.

Geomechanical modeling of fracture density under uniaxial stress state
In Earth's crust, when a rock is subjected to a 3-D stress condition, the strain energy density at a given point within the rock mass can be expressed as follows [46], [21]: where σ 1 , σ 2 and σ 3 were the maximum principal stress, intermediate principal stress and minimum principal stress, respectively (MPa); E was the elastic modulus and μ was the Poisson's ratio.
According to brittle fracture mechanics theory and maximum tensile stress theory, brittle rock will break when elastic strain energy accumulated in the brittle material equals to the energy demanded for generating fractures per unit volume of the element [47], [48]. Generally, brittle macro-fractures occur with strain energy releasing, especially when the surrounding three-dimensional stress state reaches the rock's strength [49]. At this time, part of the strain energy will be released as the surface energy of the new fractures while the rest will be released in form of elastic waves. Nonetheless, compared with the fracture surface energy, elastic wave energy was so weak that it can be neglected. We assume that all of the fractures in this study were caused by the releasing energy and therefore (Fig 4d), the law of conservation of energy must be satisfied: where D vf was the fracture volume density per unit volume (m 2 /m 3 ); $ f was the strain energy density for new fractures (J/m 3 ), also was regarded as the residual strain energy density by current strain energy density per unit volume minus the elastic strain energy density for new fractures ($ e ) (J/m 3 ); V was the unit cell volume (m 3 ); S was the surface area of the new fractures (m 2 ); J was the energy per unit area required for fractures, i.e., fracture surface energy (here, this energy was different from and had a far lower value than the theoretical value of molecular dissociation); coefficients a ¼ 1 J , and b ¼ À $ e J were relative coefficients, which could be derived by two methods: (1) Based on the values of J, $ e was calculated by physical concept and then a and b were obtained through conversion values. (2) Through experimental data, curve fitting and a statistical approach were employed to obtain the coefficients. In practice, the former has clearer physical meaning.
In most cases, there exist three types of stress environments in Earth's crust: Type Ia, Type II, and Type III and different mechanical behaviors appear under different stress environments [50]. As stated above, with a uniaxial compressive strength level of 0.85 σ c as the peak period for the initiation of fracture swarm, the corresponding strain density energy was close to $ e in theory. In this case, we assumed that $ e is strain energy density at σ = 0.85σ c . This value would be entered in a formula to calculate J and then the values of a and b could be confirmed. Based on the preceding mechanical experiments on tight sandstone in the Keshen area, the average uniaxial compressive strength of predominant mid-fine sandstone samples was identified, and we calculated $ e as 1.133 × 10 3 J/m 3 . Here, Eq (2) could be transformed into J ¼ $À $ e D vf . Calculated values of J were presented in Table 4.
Averaging J in Table 4 as � J ¼ 1055:27 J m 2 and substituting $ e ¼ 1:133 � 10 3 J m 3 into Eq (2), a quantitative relationship between strain density and fracture volume density was developed: To test and verify the reliability of $ e , a curve-fitting method was employed to calculate a and b. A comparison plot of the curve-fitting results versus data obtained by the mechanical experiments (Table 3 ) was shown in Fig 5, from which the following linear relationship was obtained: where the correlation coefficient reached 0.995, and the relative errors for coefficients a and b by two different methods were 12.6% and 17.3%, respectively, which are both within the permitted error range. It should be noted that regardless of whether the rock mass is intact or fractured, the fitting curve exhibited a positive linear relationship between parameters. Small  deviations were explained by anisotropic mechanical properties, i.e., a combination of preexisting weakness planes and unruptured rocks. Therefore, it was reliable that the necessarily overcoming elastic strain energy density for new fractures ($ e ) could be replaced by the corresponding strain energy density at 0.85σ c both in theory and practice.
Geomechanical modeling of fracture parameters under paleostress field Fracture geometric parameters. As Huang [51] and Song [33] demonstrated previously, the relationship between fracture volume density and strain energy density under triaxial compression experiments is no longer presented as a regular linear relationship. In other words, the energy accumulated in the brittle rocks must not only overcome the molecular internal cohesive forces, but also overcome the energy expenditure that naturally occurs to resist the confining pressure [51], [33]. Based on the above geomechanical model under uniaxial compression, the relationships between fracture volume density and strain energy density under triaxial stress state, tensile stress state were derived as follows: 1. Commonly, tectonic fractures are divided into tensile and shear fractures based on the formation mechanism, and they can be discriminated with the Griffith Criterion and Coulomb-Navier Criterion (also called the Mohr-Coulomb Criterion), respectively [52], [33], [53], [13]. For the fractures under confining pressure conditions, if there existed only compressive stresses (σ 3 � 0), the Mohr-Coulomb criterion would be selected, which was The Mohr-Coulomb Criterion suggested that a shear fracture only formed if the rock cohesion (C 0 ) was exceeded, which was depended on the magnitude of normal stress along a fracture plane. The relationships between fracture density, aperture and strain energy density, stress-strain were written as follows 8 > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > : Where ϕ was the internal friction angle (˚); θ is the angle between normal to the newly formed fracture plane and the maximum principal stress (˚) (Fig 6a); σ 1 was the maximum principal stress (MPa); σ 2 was the intermediate principal stress (MPa); σ 3 was the minimum principal stress (MPa); σ p was the rupture stress under action of σ 3 , different from the maximum principal stress; σ t was the tensile strength (MPa); J 0 was fracture surface energy with no confining pressure or under uniaxial compressive stress (J/m 2 ), ΔJ was the additional surface energy caused by confining pressure σ 3 (J/m 2 ); E was the elasticity modulus with no confining pressure (GPa); b was the fracture aperture (i.e. paleo-aperture) (m), referring to Fig 6b; D If was the fracture linear density (1/m); ε 3 was the tensile strain under current state of stress, dimensionless parameter; ε 0 was the maximum tensile strain, dimensionless parameter, corresponding to tensile strain when fracture beginning to form; E 0 was the proportionality coefficient related to lithology; and L 1 , L 2 , L 3 were side length of the selecting representing element volume (Fig 6a). 2. When there existed tensile stress, for brittle tight sandstone material the Griffith Criterion was used. When σ 3 � 0 and (σ 1 + 3σ 3 ) � 0, the applied failure criterion was given as When (σ 1 + 3σ 3 ) < 0, the failure criterion was simplified to σ 3 = −σ T , θ = 0.
If the failure criterion was reached, the relationships between fracture density, aperture and strain energy density, stress-strain were expressed as

> > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > :
According to the above description, when (σ 1 + 3σ 3 ) > 0, by Griffith criterion, there should be: cos 2y ¼ s 1 À s 3 2ðs 1 þs 3 Þ and y ¼ In a 3-D global coordinate system, the calculation of fracture strike and dip depends on a reasonable projection of coordinates, where the x-axis aligns with the east-west direction, the y-axis aligns with the north-south direction, and the z-axis aligns with vertical direction. Thus, the strike angle (α) could be characterized in the x-o-z plane using the normal vector of its fracture surface. Here � n ¼ fl m ng denoted the cosine of the normal vector projected in the x-o-z plane, whose inclination angle with respect to z-axis (α Z ) was obtained via the following: where α Z = arctan(−l/n) (Fig 7). Geologically speaking, in the a 3-D global coordinate system, the inclination angle between the lx + my + nz = 0 plane and y = 0 plane (α dip ) was defined as the inclination angle between the fracture surface and x-y plane, which was given by: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi l 2 þ m 2 þ n 2 p ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 0 2 þ 1 2 þ 0 2 p ¼ jmj ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Fracture seepage parameters. Generally, for multiple sets of fractures, the relationship between fracture porosity and fracture volume density and aperture was illustrated as: where m was the number of fractures; b i was the aperture of group i (m); D vfi was the volume density of group i (m 2 /m 3 ); φ ft was the total fracture volume, as decimal type. When the fracture permeability is calculated under the principal stress space, the principal stress value of the space node can be obtained by FE simulation, but the principal stress direction is inconsistent with the coordinate axis direction of global Cartesian coordinate system. Therefore, the mechanical model of fractures permeability cannot be directly applied. For arbitrarily distributed 3-D fractured media regions, according to Kronecker conversion formula [54] and principle of coordinate transformation [55], the permeability of three coordinate axes for a group of parallel fractures in the global coordinate system could be characterized as: Thus, the total fracture permeability was calculated using Eq (15), namely: where, m was the number of fractures; D lfi was the linear fracture density of group i (1/m); K fXi , K fYi a and K fZi were the permeability in the x-axis, y-axis and z-axis directions, respectively (m 2 ); α 11i , α 12i and α 13i were the angles between σ 1 and axis-x, y, z for group-i fractures (˚); α 21i , α 22i and α 23i were the angles between σ 2 and axis-x, y, z for group-i fractures (˚); α 31i , α 32i and α 33i were the angles between σ 3 and axis-x, y, z for group-i fractures (˚); θ was the angle between the short axis direction of fracture surface and σ 1 , namely the rupture angle of rock (˚); K T fX , K T fY and K T fZ were the total permeability in the x-axis, y-axis and z-axis directions, respectively (m 2 ); (n f1i , n f2i , n f3i ) were the three components of unit normal vector of group-i fracture surfaces, respectively, which could be calculated by Eq 16:

Fracture aperture under present-day stress
For our modeling purposes, we assumed that the paleostress field generated fractures, while the present-day stress field induces minor changes in fracture size but does not produce new fractures. We made further modification to the existing palaeomechanical models. Considering the influences of the normal stress and shear stress on the aperture of a fracture, Hicks et al. [56] and Jing et al. [57] derived an equation to calculate the aperture of a fracture due to the current in situ stress field as follows: where b 0 and b m were the original and current aperture of the fracture (m), respectively; σ' n was the effective normal stress (MPa), namely was the result of the normal stress acting perpendicular to fracture plane minus the fluid pressure acting inside fracture [8], [9]; b res was the residual aperture of the fracture (m); and σ nref was the corresponding effective normal stress (MPa) when the fracture aperture decreased by 90%. As for the parameter σ nref , many researchers [58], [59], [60] have proposed different determination ranges based on their experimental data. Through a series of permeability tests under variable confining pressure, Qin [60] showed that the reasonable value 30 MPa should be applied to low-permeability sandstones with uniaxial compressive strength ranging from 30 to 70 MPa. Because of minimum effects on total seepage capability, or difficult to be determined by conventional means, the parameter b res could be often ignored in calculations. Additionally, as a critical parameter, the effective stress s 0 n was the result of the normal stress acting perpendicular to fracture plane (compression assumed positive) minus the fluid pressure acting inside fracture [8], [9]: s 0 n ¼ s n À P. Therefore, the equation for calculating fracture porosity under present geostress field is modified as: Generally, the higher the filling degree in a fracture, the lower its aperture and porosity are expected to be. Here we introduced a novel method of mineral filling coefficient to characterize the filling degree in fractures form the perspective of natural hydraulic fracturing. In detail, under two limit conditions that when a fracture was most or full filled by minerals (calcite, quartz, gypsum or mud) the filling coefficient C was equal to 1, and when the fracture was almost unfilled, the C was zero. When the fracture was approximately half or more than filled, the C was assigned the value 0.5, but when the deposited minerals accounted for a quarter of the fracture volume, the C might be assigned 0.25 as the estimator. Therefore, this leads to two simpler expressions for estimating the effective fracture aperture and porosity relative to filling degree required for opening-mode fracturing: where b fe and ϕ ef were the effective fracture aperture and porosity, respectively; b m was the original fracture aperture before fillings occurred. In this way, the calculation formula of effective fracture permeability under the current stress field is modified to:

Geomechanical simulation of tectonic stress field
Geomechanical modeling with Finite Element Method (FEM) was conducted to simulate the paleostress and current stress distributions for prediction of fracture development and distribution. The general-purpose finite element code ANSYS (15.0, ANSYS, Inc., USA) was selected for this study because it is well-suited for analyzing geomechanical problems over a wide range of scales in one, two, and three dimensions [61], [23], [62], [27]. The initial threedimensional model geometries were constructed for the restored geologic area constrained by field measurements during Pliocene Kuqa Period [42] and incorporated the generalized mechanical stratigraphy of the Bashijiqike Formation (K 1 bs 1 , K 1 bs 2 and K 1 bs 3 ). The average thickness of the K 1 bs 1 , K 1 bs 2 and K 1 bs 3 member within the FE model was set to be 50 m, 180 m and 65 m, respectively. The target zone associated with the upper cover layers (i.e. K, E, N Formations) and basement (i.e. J Formation) were considered as an isolated body or boundary conditions for simulation (Fig 8a). All model layers were discretized using primarily threenode triangle plane strain elements along with some four-node quadrilateral elements (311,226 total elements in the model).
Based on regional analysis and acoustic emission of rocks, the Kuqa depression has primarily experienced three important thrusting movements during the Himalayan Period and Neotectonics Period, which can be defined as Early Initial Compressional Stage, Mid-term Strong Thrusting and Uplift Stage with greater slip displacement and Late Uplift and Recoil Stage with higher fold amplitude, respectively [23], [21], [43], [63]. It was important that the Kuqa tectonic movement after Pliocene sedimentary was the crucial time of fracture generation in the Keshen reservoir, and the direction of maximum tectonic stress in the period was identified as nearly NNW352˚ [21], [63], with a magnitude of about 370 MPa. Then, based on paleotectonic deformation, a methodology that included both inverse and forward methods was used to identify the magnitude of minimum paleostress field, which was about 70 MPa. Considering the Kuqa period after Pliocene sedimentary as the crucial time of fracture formation in the Bashijiqike reservoir, the cover layers are set to approximately 5500 m thick. Consequently, in the first step, we selected the southern edge (hinterland side) as the fixed boundary in the simulated stratigraphy so that the Keshen anticline could uplift if the stress state was appropriate (Fig 8a). The second step imposed a deformation of fault-related fold by forming a displacement boundary condition along the north (foreland) side of the model. The base of the model was pinned, a gravity load was applied to the entire model domain, and the system was allowed to reach equilibrium.
The state of in-situ stress is described by a stress tensor, which comprises three orthogonal principal stresses, each with an orientation and magnitude [44], [64], [65]. In this study, XMAC logging, acoustic emission test, hydraulic fracturing test and differential deformation analysis were used to determine the magnitude of maximum horizontal stress (σ H ), minimum horizontal stress (σ h ) and vertical stress (σ v ) of drilling wells in the Keshen area ( Table 5). The horizontal stress orientation in the study area was determined from paleomagnetic correction, wave velocity anisotropy test, and drilling-induced fractures (DIF). As shown in Table 6, the results indicated that the orientation of maximum horizontal stress in wells had obvious zoning characteristics, mainly ranging from NW322˚to NE 51˚. Interestingly, the direction of σ H in the central region was nearly north-south direction, distinctly different from the that of western and eastern regions. From Table 5 we could see that the maximum horizontal stress values at different depths of the study stratum were significantly higher than that of vertical stress, which indicated that a Type III stress environment (i.e., σ H >σ v >σ h ) was widely developed in Keshen area at present. The study area is primarily affected by the southern Indian plate and northern Kazakhstan-Junggar plate, with approximately twice the force as that from the North China plate [66]. After several design improvements and repeated tests, the reasonable scheme for mechanical boundary condition (including stress and displacement constraints) were finally confirmed. A magnitude of 122 MPa was applied to the E90˚and W270˚boundaries of the nested model, respectively, as horizontal σ h (Fig 8b). Simultaneously, A magnitude of 158 MPa was applied to the N360˚and S180˚boundaries of the nested model, respectively, as horizontal σ H . The vertical stress could be calculated and applied automatically in software based on gravitational acceleration. As for the variations of stress orientation with the slipping of boundary faults, we should take into account the extra shear stresses while applying loading conditions. In detail, we applied dextral and sinistral shearing stresses to the west and east boundaries, respectively, both with a magnitude of 130 MP.
To date, there have been hardly any published anisotropic models of the paleo-stress field of compressional basins, no far any other stress evolution models [67], [68]. Here our study was focused on the paleo-folding stress field, as opposed to building a detailed and complicated model for the present-day field. For simplicity, the mechanical parameters density of model, such as Young's modulus, and Poisson's ratio variables were set to the same value for each layer, respectively. Based on the above mechanical experiments, the material properties (i.e., density, Young's modulus, Poisson's ratio, friction angle, and cohesion) were assigned to the elements representing the various lithologies ( Table 6). The processing method used to determine the material properties of fault zones is important to the outcome of the geomechanical modeling because it directly influences the results of numerical simulation, for instance, the distributions of stresses and fractures [27]. Generally, the Young's modulus of fault zone is just 65-85% of the elastic modulus of a corresponding normal stratum [27], [28], which indicates that the strength parameters of fault zone is obviously lower than that of intact rocks in tight sandstone areas.

Prediction of fracture distribution
According to above deduced geomechanical models with respect to the composite rock failure criterions, the initial fracture parameters (density, aperture porosity and permeability) related to stress, strain and strain energy density can be calculated or directly extracted from the numerical simulations of paleostress field during the Kuqa tectonic movement of the Late Himalayan stage. Based on the simulated ancient fracture parameters and considering the impact of current stress field on fracture properties, we could predict the spatial distributions of fracture aperture and porosity for each element/layer, which were displayed on Ansys platform. As for the spatial filling coefficients of fracture in modeling process, one effective method was to use cokriging which considered the paleofluid migration direction to interpolate the filling degree of different layers in the whole region, which was primarily constrained by well-point data.
In this paper, compressional stress was positive and tensile stress was negative. In the K 1 bs, the maximum principal stress (σ 1 ), indicative of compression, ranged from 360 MPa to 418 MPa (Fig 9a). The minimum principal stress (σ 3 ), indicative of both compression and tension, mainly ranged from -5 MPa to 71.4 MPa (Fig 9b). Fig 9b showed that weak tensile stress (i.e., negative values) was only distributed in the top of Keshen anticline, and consistent with the long axis of fault-related fold, indicating the most probable development zones of tension fractures. The distributions of three principal stresses were similar to each other, all of them were primarily fault-controlled and secondly fold-controlled. For example, after rapid energy releasing in the rock, the fault zone showed lower stress values, and a transformation of stress field appeared near major faults, such as the higher stress values in the footwall and lower stress values in the hanging wall. Similarly, in the simulated current stress field, the distributions of maximum horizontal stress (σ H ), minimum horizontal stress (σ h ), and vertical stress (σ v ) were similar, and the gravity and boundary stresses played the major role within Keshen area (Fig 9c and 9d). The maximum horizontal stress was horizontal, with values varying between 90 MPa and 263 MPa, indicative of compression. Relative lower σ H values mainly were located in the middle part of the gas field, and gradually increased to both east and west directions. Whereas the σ h values decreased from north to south in the gas field, which indicated that the stress concentration was primarily come from the resistance of southern plate.
After superposition of paleostress field and current stress field in Keshen area, the current simulated linear fracture density ranged mainly from 1 to 9 m -1 in the plane (Fig 10a, 10b and  10c). Horizontally, areas with well-developed tectonic fractures were mainly located in fold limbs, fault zones (such as the Wells A2-7, B2-3, B2-8 and B2-9), locations of changes in the orientation of faults, zones on footwall (such as the Wells A1-1, A1-2), and southern regions around front limb. Vertically, the bottom layer had a relative higher density than the top layer, such as fracture linear density of K 1 bs 3 is higher than that of K 1 bs 2 and K 1 bs 1 , indicating that it was highly fractured and more brittle. Additionally, as a comparison, the high-value areas of current fracture aperture were roughly counter to the areas with higher fracture linear density, such as on top of the fold fracture aperture was high, but the fracture density was low, and in northern limb of fold the fracture density was high, but the fracture aperture was low (Fig 10d). Vertically, the current fracture aperture decreased with increasing strata depth, which was just opposite to the distribution of fracture density. However, the predicted well- developed and relatively well-developed areas of fracture density coincided with the areas with high fracture aperture, including the eastern fault zones near the Wells A2-6 and A2-7.
As an important factor in determining the quality and yield of tight gas reservoirs, the distribution of current fracture porosity was very similar with that of fracture aperture. As shown in Fig 10e, the most well-developed porosity zones were mainly located in (1) the northwest, southeast, and middle of the study area; (2) around anticlines; (3) around faults; and (4) the upper thin sandstones. In the anticline area, vertically, fracture porosity gradually decreased from the anticline core to the wings. Comparison of different-scale faults showed that the influence of faults at the high point on fracture porosity was more significant than that in lower position. In fractures associated with faults at the lower position, core observations and other evidence suggested that cement might be easily deposited simultaneously with fracture opening and propagation, partially filling and decreasing the large aperture parts [14]. As illustrated in Fig 10e, the fractures simulated in the K 1 bs can be divided into four sets, namely set I (NW-SE), Set II (NNE-SSW), set III (NE-SW) and set IV (nearly EW), where the former two sets were predominant, which basically coincided with orientations in drilling cores and FMI. Set I fracture strike 295˚, oblique to the fold axial trend and boundary faults with acute angle, and were interpreted as a regional fracture set that was present before the Keshen anticline and with the stress field of sinistral rotation and compressive shearing. Set II fractures striking 20å nd set III fractures striking 65˚near fault zones, were interpreted as composite conjugate fracture set associated with faulting related to the NNW-oriented Himalayan compression during initial anticline growth. Set IV fractures striking 96˚, nearly parallel to the fold trend, were found mainly within the crest and were interpreted to have formed in response to local tensile stresses during folding.
In general, the high permeability zone in Keshen area was concentrated in the high structure and fault zone, and the permeability value decreased with depth (Fig 11a, 11b and 11c). At the same time, the permeability of fractures in different directions was mainly related to the occurrence of fractures. The fractures in the study area were dominated by vertical joints, so the permeability value in the vertical direction was the highest. The main trend of fractures in the study area was nearly S-N direction, so the permeability in the north-south direction was obviously higher than the east-west direction. The cause of this change in fracture permeability could be briefly explained by stress state transitions (Fig 9). In the shallow part of Keshen anticline, the maximum principal stress was north-south direction, the minimum principal stress was vertical, and the intermediate principal stress was east-west. Under this condition, the permeability had the characteristics that the east-west direction was greater than the north-south direction, and the north-south direction was greater than the vertical direction. In the deep part, after reaching a certain depth, the stress field changed, the vertical stress increased to the intermediate principal stress, and the minimum principal stress was the east-west direction. It was observed that it was only necessary to rotate the shallow geostress model 90˚along the σ 1 axis to obtain the deep-ground stress model, and the permeability in the corresponding direction also changed. As a result, the permeability in the corresponding direction also changed, the maximum permeability appeared in the vertical direction, the minimum permeability became east-west, and the north-south permeability was centered.
As shown in Table 7, the overall geomechanical FE modeling results at the wells were in agreement with in-situ core observations and formation micro-imaging (FMI). The average fracture density at Wells A2-7 and B2-8 in K 1 bs was as high as 9.36/m and 5.62/m, respectively, according to the core and FMI in-situ measurements, while modeling results yielded 9.04/m and 6.18/m, for the same two wells, respectively. Comparison diagram of all 12 sample spots in the Bashijiqike Formation indicated little error between the computed and in-situ results. Only a few samples showed relative errors in average linear fracture density greater than 20% and these larger relative errors were distributed along fault ending, such as Wells B2-9, A2-5 and A2-8. Compared with the traditional core observation method, important in-situ physical parameters (i.e., aperture and porosity) of fractures are more easily obtained by FMI interpretation and CT scanning accurately and continuously.
It is worthwhile to note that in Table 7, the average relative error between the predicted fracture apertures and the measured fracture apertures in was not more than 20%, except for one well exceeding 35% (S3 Appendix). Similarly, based on data analysis, the majority of predicted fracture porosities had nice fit with the measured fracture porosities (average relative error<20%), indicating that the present prediction of tectonic fractures by geomechanical superposition method in the brittle tight sandstone regions was believable. As for the existing larger errors, it was inferred that probably the stress concentration phenomena induced by lithology variation or fault zone structure were still hard to be correctly reflected by normal FE mesh partition, which directly led to the uncertainty of simulation results including the development and distribution of tectonic fractures. As shown in Fig 11d, in the eastern part of the Keshen anticline, although there was large difference between the unobstructed flow before and after the reservoir reformation, the distribution trend of the latter was basically consistent with the fracture permeability distributions.

Conclusions
An adapted geomechanical method was used to calculate the multiple parameters of fracture governing the development of the tight sandstone reservoirs for unconventional resources in the Kuqa Depression within Tarim Basin, China. This method was selected for its superior applicability in developed fracture areas after key tectonic movement and its adaptability to fit the FE modeling. Its predictions also agreed well with present in-situ core observations and FMI interpretation.
From the rock mechanics experiments under continuous loading, fracture evolution of tight sandstone (sometimes containing a small number of micro-fractures) was summarized into three stages: (1) initial compaction stage, (2) propagation stage, and (3) coalescence stage extracted from three-dimensional CT scan. When reaching peak shear strength, tight sandstone rocks experienced volume expansion and micro-fractures rapidly connect to each other. This value tended to be τ = 0.85 σ c (σ c < 200 MPa), which was the critical sliding value mentioned in Byerlee's law [69]. From a microcosmic point of view, there existed only three failure modes in brittle rocks: shear fracturing, tensional fracturing, and composite fracturing, which depended mainly on confining pressure, principal stress intensity, and the pre-existing plane of weakness. After many empirical experiments, a relationship between the fracture volume density and stress-strain of tight sandstone reservoirs was finally established. To this end, considering the superimposed effect of paleostress and current stresses, the expressions or geomechanical models were obtained for characterizing the ancient fracture parameters after key tectonic movement and present-day fracture property parameters under different stress states. These expressions/models could be incorporated into a program running on a FE simulation platform, which had become a popular and effective modeling approach in the prediction of fracture generation and spatial distribution. The geomechanical model was, however, limited to the strength prediction for fractured anisotropic rocks and triaxial testing conditions. Further study was suggested to extend the composite modified criterion for anisotropic rock masses and polyaxial testing conditions, with emphasis on the effect of heterogeneous pre-existing fractures.