Simulation of tectonic stress field and prediction of fracture distribution in shale reservoir*

In this paper, a finite element-based fracture prediction method for shale reservoirs was proposed using geostress field simulations, uniaxial and triaxial compression deformation tests, and acoustic emission geostress tests. Given the characteristics of tensile and shear fractures mainly developed in organic-rich shales, Griffith and Coulomb – Mohr criteria were used to calculate shale reservoirs' tensile and shear fracture rates. Furthermore, the total fracture rate of shale reservoirs was calculated based on the ratio of tension and shear fractures to the total number of fractures. This method has been effectively applied in predicting fracture distribution in the Lower Silurian Longmaxi Formation shale reservoir in southeastern Chongqing, China. This method provides a new way for shale gas sweet spot optimization. The simulation results have a significant reference value for the design of shale gas horizontal wells and fracturing reconstruction programs.


Introduction
For low-porosity and low-permeability shale reservoirs, the nano-scale pores in the matrix have basically no seepage capability. Therefore, fractures not only provide important space for hydrocarbon storage, but also provide efficient channels for hydrocarbon migration [1][2][3]. The great success of the marine organic-rich shale gas industry in North America shows that natural fractures can promote the large-scale accumulation of hydrocarbons in shale reservoirs [4][5][6].
Fractures are a key factor in obtaining high yields in shale reservoirs [7][8][9][10][11][12]. A large number of oilfield data around the world show that the degree of fracture development in tight reservoirs is closely related to productivity [13][14][15][16]. For example, the degree of fracture development of Paleozoic marine shale in North America is positively correlated with total gas content and free gas content. The success rate of shale gas exploration in the fractured zone is high. In addition, the natural gas productivity in the organic-rich marine shale reservoirs of the Lower Paleozoic in the Sichuan Basin of China is also positively correlated with the degree of fracture development [17][18][19].
In this paper, a finite element-based fracture prediction method for shale reservoirs was proposed using geostress field simulations, uniaxial and triaxial compression deformation tests, and acoustic emission geostress tests. This technology has achieved good application effects in the prediction of fracture distribution in the Lower Silurian Longmaxi Formation shale reservoir in southeastern Chongqing, China. Moreover, it provides a new way for shale gas sweet spot optimization, and the simulation results have important reference value for the design of shale gas horizontal wells and fracturing reconstruction programs.

Materials and methods
Experiments. In this paper, acoustic emissions and rock mechanics experimental tests were used to obtain the rock mechanical properties and paleostress of the target shale. The tests were completed in the Beijing SGS Rock Physics Laboratory. A GCTS petrophysical testing system was used for rock mechanics testing. The pressure sensor error of this instrument was less than 1 %, the displacement sensor ranges were between ±50 mm, and the strain accuracy was 0.0001 mm. In addition, the acoustic emission instrument was the SAMOS TM acoustic emission detection system. Its core component is the PCI-8 acoustic emission function card that processes the PCI bus in parallel. It has 8 channels of realtime acoustic emission feature extraction, waveform acquisition and processing capabilities on one board. Modern digital signal processing technology (DSP) is adopted, which is currently the most advanced acoustic emission processing system in the world.
Finite element model. This paper used the finite element method to simulate the tectonic stress field and then predicted the plane distribution of tectonic fractures based on rupture principles. The core technology of this method is to establish an accurate geological model, mechanical model, and calculation model of the simulated area. The measured rock mechanical property parameters and paleostress values were used to calibrate the fake stress field (Fig. 1). The organic-rich shale mainly develops tensile and shear fractures. Therefore, Griffith and Coulomb -Mohr failure criteria were used to calculate the tensile and shear failure rate, respectively. Finally, proposed the comprehensive rupture rate based on the coupling results of tensile and shear ruptures (Fig. 1). (1)

Fig. 2. Grid model of Longmaxi Formation shale in southeastern Chongqing area Рис. 2. Сеточная модель сланцев формации Лунмаси в юго-восточной части района Чунцин
Equation (1) can be simplified as: (2) In the formula, Ni, Nj and Nm are the morphological function or shape function of the element displacement, [N] is the shape function matrix, and [δ] e is the nodal displacement component matrix.
The strain of the element is a geometric equation: When equation (1) is substituted into equation (3), the strain matrix of the element can be obtained: (4) In the formula, the conversion matrix [B] is a geometric matrix.
For each element, the maximum principal stress is obtained through coordinate transformation: The maximum principal stress σ1 and the minimum principal stress σ3 can be obtained by solving the above formula.

www.nznj.ru
Griffith and Coulomb -Mohr criterions. Under the action of regional tectonic stress, there are two main types of ruptures inside the rock: tensile and shear fractures [8][9].
The expression of the plane rupture criterion of Griffith theory is: When σ1 + 3σ3 ≥ 0, the rupture criterion is: (6) When σ1 + 3σ3 < 0, the rupture criterion is: In the formula, σ1 is the maximum principal stress, MPa; σ3 is the minimum principal stress, MPa, and σT is the tensile stress of the rock, MPa.
The Coulomb -Mohr criterion believes that the shear failure on a plane is related to the combination of the normal stress σ and the shear stress τ. The Coulomb -Mohr shear rupture criterion can be expressed as: In the formula, |τ| is the shear strength of the rock, MPa; σ is the normal stress, MPa; C is the cohesive force, MPa; φ is the internal friction angle, ˚; tanφ is the internal friction coefficient.
Comprehensive rupture rate. In this paper, the tensile rupture rate It and the shear rupture rate In were introduced to characterize different types of fractures: = / . (9) In the formula, σT is the effective tensile stress, MPa; σt is the tensile strength of the rock, MPa. When It ≥ 1, tensile ruptures will occur. = /| |. (10) In the formula, τn is the effective shear stress, MPa, and |τ| is the shear strength of the rock, MPa. When In ≥ 1, shear ruptures will occur.
The rupture mode of shale is a comprehensive reflection of tensile and shear stresses [9][10]. Therefore, in order to better quantitatively characterize the development degree of structural fractures in shale reservoirs, a comprehensive fracture coefficient was proposed. = ( + )/2. (11) In the formula, a and b are the ratios of tensile and shear fractures respectively. In this paper, a : b = 3 : 2. Similarly, when Iz ≥ 1, the rock reaches a fractured state, and the higher the comprehensive fracture rate value of shale, the greater the fractured degree.

Results
Palaeo-stress based on acoustic emission. In the simulation of in-situ stress, the assignment of reasonable rock mechanics parameter attributes of the geological model of the target layer is essential. Furthermore, the assigned geological model is converted to a mechanical model. According to regional tectonic movement and acoustic emission tests, it is believed that during the Yanshan period, the tectonic activity in southeastern Chongqing was the strongest (148.8 MPa maximum tectonic stress); followed by the Himalayan movement (122.5 MPa maximum tectonic stress) ( Table 1).
Rock mechanics parameters. The faults in the study area were divided into first-order, secondorder, and third-order faults. At the same time, the fold areas were split into slot folds, battlement folds, and barrier folds. The rock mechanics test results of different types of shales in the Longmaxi Formation in the study area are shown in Table 2.
Mechanical properties of fault and fold zones. The fault zone was defined as a "weak zone" whose elastic modulus was 50-70 % of the ordinary sedimentary strata. At the same time, the Poisson's ratio was more extensive than that of the ordinary sedimentary strata, and the differences between them were between 0.02 and 0.1. The folding zone was identified as a "tough zone", and its elastic modulus was 1.5 to 3 times that of the normal sedimentary formation. At the same time, the Poisson's ratio was smaller than the normal sedimentary formation, and the differences between them were between 0.01 and 0.15 (Table 3).  Discussion Tectonic stress field distribution. It can be seen from the simulation results of the tectonic stress field (Fig. 3) that the maximum principal stresses of the Longmaxi Formation shale reservoir in southeastern Chongqing were concentrated between -217.404 and -4.109 MPa. Positive values were defined as tensile stress, and negative values were defined as compressive stress. The maximum principal stress inside the fault zone was lower than that of ordinary sedimentary strata, and the stress intensity values were mainly distributed between -46.768 and -4.19 MPa. For areas with underdeveloped faults, the maximum principal stress value distribution ranged from -103.647 to -46.768 MPa.
The rocks inside the fold zone are severely deformed, especially the rocks at the shaft and turning ends of the folds are more severely deformed. The stress in these structural parts will

Fig. 3. Distribution of maximum principal stress of Longmaxi Formation shale in southeastern Chongqing area Рис. 3. Распределение максимумов основных нормальных напряжений в сланцах формации Лунмаси в юго-восточной части района Чунцин
www.nznj.ru be highly concentrated under the premise that there is no fault damage and cannot release the stress. Suppose the fold is clamped by reverse faults, such as the southwestern Huayuan and Pengshui west-trending fault fold belt, or the faultrelated folds adjacent to the fault, such as Longshan and Xiushan areas. In that case, the maximum principal stress value will be higher. In addition, the closer the fold is to the fault, the more obvious the stress gradient changes. In addition to the fold mentioned above belts, the fault belt's end and the fault's turning end are also the transition areas from the broken rocks inside the fault belt to continuous strata. The rocks in these areas are at the edge of ruptures. Therefore, the stress value is higher. The maximum principal stress distribution in these areas ranges from -217.404 to -103.647 MPa. The shear stresses of the target shale reservoir in the study area ranged from -4.707 to 49.222 MPa (Fig. 4). Among them, positive values were defined as left-handed and defined negative values as right-handed. The structures of the study area showed obvious strike-slip characteristics, and the NNE-trending "S"-shaped faults and folds had the attributes of counterclockwise rotation and twisting. It can see from Figure 4 that the shear stress value in the study area is mainly positive, reflecting the counterclockwise left-handed shear stress field in the Himalayan period in southeast Chongqing. The simulation results are consistent with the compression-torsional strike-slip structural deformation characteristics shown in the Himalayan period in the study area.
Prediction of fracture distribution. In this paper, the degrees of fracture development in shale reservoirs were divided into five levels ( Table 4).
According to Figure 5, there are widely distributed fractures of grade I-IV in the eastern part of the study area, and the fracture development coefficients are mainly 1.4-4. Among them, the sizeable trough-shaped fold axis in the northern part of Huayuan is an area with highly developed fractures (level IV), and the fracture development coefficient has even reached above 4 (the formation was severely broken). The southern (Xiushan, south of Huayuan) and northwestern (Lianhu area) areas are favorable areas for level II and level I fractures, respectively. The western part (Wulong, Pengshui) mainly develops firstlevel fractures, and the fracture development coefficients are between 1-1.4.

Fig. 5. Distribution of fracture development coefficient in Longmaxi Formation shale in southeastern Chongqing area Рис. 5. Распределение коэффициента развития трещин в сланцах формации Лунмаси в юго-восточной части района Чунцин
The areas with high TOC content and brittle mineral content in Longmaxi Formation shale reservoirs are mainly located in the deposition center of black shale, namely Lianhu-Qianjiang and South Longshan areas. With the same comprehensive fracture coefficient, shale fractures with high TOC content and brittle mineral content are more developed. The eastern and southern parts of the study area have the most developed fractures, especially in the areas adjacent to the faults and the relatively strong-deformed grooved fold shafts, where some normal tensile faults have appeared. The southern part of the study area is the development zone of sandy shelf facies shale. The rock elastic modulus is high, and Poisson's ratio is low. It is prone to develop fractured under the action of external forces.

Conclusions
(1) The core of the numerical simulation method of the tectonic stress field lies in the establishment of the accurate geological model, mechanical model, and mathematical model. Given the particularity of shale reservoirs, the geological model must be used as the basis during the simulation process, and shale types and rock mechanical properties must be classified according to rock facies.
(2) For the interior of the fold belt in the study area, especially the shale reservoir near the axis of the fold and the turning end, has suffered severe structural deformation, which is a highly concentrated area of stress. The end of the fault zone and the turning end are the continuous transition area from the broken shale inside the fault zone to the ordinary sedimentary strata, which is at the edge of rupture and has high-stress values. The black carbonaceous, siliceous, and calcareous shale of shallow sea shelf facies with stable distribution and weak structural deformation in the deposition center have high elastic modulus and low Poisson's ratio. These brittle shales are prone to develop structural fractures.
(3) The quantitative prediction of shale fracture distribution cannot be based on a single factor as the criterion. Otherwise, it will cause onesided and limited results. The total fracture rate that affects the development of fractures in shale reservoirs should be considered as much as possible. This study found that fracture development areas are mostly concentrated in high-stress areas with severe structural deformation.