Prediction of Fracture Opening Pressure in a Reservoir Based on Finite Element Numerical Simulation: A Case Study of the Second Member of the Lower Triassic Jialingjiang Formation in Puguang Area, Sichuan Basin, China

All stages of oil and gas exploration and development involve the study of in-situ stress. Since the traditional two-dimensional and three-dimensional homogeneous models can no longer fulfil the requirements of research and production, numerical simulation of the stress field has become an effective study method. In this study, we took the Jia 2 member in Puguang area as a case to establish a geological model and a mechanical model based on the tectonic framework and the distribution characteristics of the rock mechanical parameters, respectively, and loaded the model with the present-day in-situ stress state calculated from the logging data as the boundary conditions. The simulation results show that 1) the orientation of the maximum horizontal principal stress in the study area is near E-W, and the in-situ stress orientation is locally deflected due to the influence of faults; and 2) the magnitude of in-situ stress is predominantly affected by the burial depth and lithology, and the minimum horizontal principal stress, maximum horizontal principal stress, and differential stress are mainly concentrated in the ranges of 30–60, 50–80, and 10–40 MPa, respectively. We also analysed the opening sequence of the multiple fracture systems during development, using the present-day stress field model. The analysis revealed that the E-W fractures will open first, and the continuously increasing operating pressure will lead to formation breakdown, producing a fracture network.


INTRODUCTION
In-situ stress is a natural force that objectively occurs in a rock mass (Jing et al., 2011). It is narrowly defined as the internal stress confined in the Earth's crust (Xie et al., 2008), and broadly defined as the state of stress in the Earth in different geological periods. According to the time of occurrence, in-situ stress can be divided into paleo-stress and present-day in-situ stress. Paleo-stress typically refers to the stress before the Middle Pleistocene of the Quaternary, whereas the present-day in-situ stress refers to the stress state in the present formation (James, 1983;Zeng and Tian, 1998). Every process of oil and gas exploration, development, and production involves the study of in-situ stress. Paleo-stresses controlled the type of basin, as well as the generation, development, and combination, of geological structures. They also affected the migration and accumulation of oil and gas. Present-day in-situ stress influences the injectionproduction well arrangement, hydraulic fracturing plan, hydraulic fracture propagation, casing deformation, and wellbore stability in oil and gas development.
At present, the methods for evaluating in-situ stress states can be generally divided into three categories: field measurement and analysis, logging, and stress field numerical simulation. The field measurement and analysis methods are currently the most accurate means to obtain in-situ stress data, mainly including those based on casing stress relief (Cai et al., 1997;Qin et al., 2018), acoustic emission Li and Dong, 2009), hydraulic fracturing (Zhou et al., 2015), differential strain, seismic velocity anisotropy (Cheng et al., 2008) and micro-seismic wave detection (Liu et al., 2006). These methods have high accuracy but the measurement points are discrete and not able to represent low-stress states of the whole or intact rock mass. Moreover, they cannot be applied on a large scale because of the high cost. In the logging method, the maximum horizontal principal stress orientation is detected by imaging logging based on the information on wellbore collapse and drilling-induced fractures. The maximum and minimum horizontal principal stress magnitudes are then calculated with the conventional logging data and the corresponding stress field model. Thus, the continuous variation patterns of the maximum and minimum horizontal principal stresses with depth are obtained. This method is easy to implement but the calculation results generally deviate from the true values to an extent; they need to be corrected with experimental results to improve the calculation accuracy. In addition, this method can only be applied to calculate the in-situ stress state in a single well, not the inter-well stress state. Numerical simulation of the stress field is one of the most effective methods to analyse the tectonic stress field quantitatively (Hu and Li, 2015). It is the preferred method to reconstruct the structural evolution process and the spatial distribution of stress fields. The currently available numerical simulation techniques include the finite element method, finite difference method, boundary element method, and discrete element method, among which the finite element method is the optimum for stress field modelling. In recent years, owing to the substantial development of finite element theory and software technology, the performance of this method has been significantly improved from qualitative analysis to quantitative calculation, from two-dimensional to three-dimensional. Similarly, it is necessary to combine the numerical simulation approach with other methods when calculating the distribution characteristics of the stress field.
Since the three-dimensional numerical simulation of the stress field can more realistically reflect the subsurface morphology, the three-dimensional finite element numerical simulation method has been widely applied in China to evaluate the differences in vertical stress caused by the weight of the overlying strata. However, when applying this method, researchers generally simulate the entire formation as a homogeneous body (Tan et al., 1997;Shen et al., 2004;Li et al., 2007;Su and Li, 2009;Sang et al., 2010;Tian et al., 2011;Zhan et al., 2011;Dai et al., 2013;Zeng et al., 2013;Meng et al., 2014;Zu et al., 2014;Hu and Li, 2015;Wang et al., 2015;Gong et al., 2021a), ignoring the heterogeneity of the layer due to the variations in sedimentary facies and lithology. Yet, strong heterogeneity can lead to redistribution of the in-situ stress magnitude and abrupt changes in the in-situ stress direction across the layer (Zeng et al., 2010;Zu et al., 2014;Gong et al., 2019;Gong et al., 2021b), thereby directly affecting the exploitation effect of the reservoir (Gao et al., 2000;Su and Li, 2009;Zeng et al., 2013;Li et al., 2019;Gong et al., 2021a;Gong et al., 2021b). This study attempted to establish a three-dimensional, heterogeneous present-day in-situ stress model for the Jia 2 member in Puguang area based on the tectonic framework, distribution of the rock mechanical parameters, and present-day in-situ stress state calculated from logging data, which can provide an essential basis for well pattern deployment and wellbore-instability prevention during drilling. As a next step, we applied this model to predict the fracture opening sequence to provide a foundation for designing a rational hydraulic fracturing plan.

GEOLOGICAL SETTINGS
The Puguang gas field is located in the north-eastern part of the Sichuan Basin, with an exploration area of 1,116 km 2 (Xu et al., 2004;Qian and Zhong, 2009). To its north and east are the Micangshan thrust belt and Dabashan arcuate thrust belt, respectively. This gas field is not only controlled by the tectonic deformation of the basement of the Sichuan Basin, but also influenced by the deformation of the Micangshan and Dabashan thrust belts ( Figure 1). The entire study area is developed with NE-trending structures. Several boundary faults running through this area divide the region into multiple secondary units, including two positive structures, the Maoba-Leiyinpu and Shuangshixi-Qingxichang structural belts, and two negative structures, the Puguang main body and Xuanhan syncline, showing a pattern of uplift and depression. Three major sets of décollements are developed, namely, from top to bottom, the gypsum salt rock from the lower part of the Leikoupo Formation to the fourth member of the Jialingjiang Formation, the Silurian shale, and the shale of the Cambrian  (Tang et al., 2008;Zu et al., 2017).
The Puguang area is developed with gas-bearing systems of terrestrial and marine facies. Currently, the marine facies, Changxing and Feixianguan Formations, are the main production layers in the Puguang gas field. The Jialingjiang Formation is an important regional gasproducing formation in the Sichuan Basin with a good exploration prospect; the Jia 2 member has the best gas logging show. This layer has considerable undulations, with a burial depth between 2,500 and 5,500 m. The lithology of the Jia 2 member is complex, including dolostone, limestone, gypsum rock, and mudstone. According to the core analysis and testing data, the porosity of the dolostone is generally >2%, with an average of 5.8%; the maximum porosity of the limestone is 2.18%, and the average is 1.72%; and the permeability of the matrix is generally <0.25 × 10 −3 μm 2 . Thus, this layer belongs to an ultralow-porosity and low-permeability reservoir.

METHOD OF NUMERICAL SIMULATION OF TECTONIC STRESS FIELD
In this study, we used the finite element method to perform numerical simulation of the tectonic stress field (Zeng et al., 2010;Gong et al., 2019;Gong et al., 2021b). This method is based on the variational principle and approximate interpolation discretization, breaking through the limitations of traditional numerical simulation that requires tectonic reconstruction from the aspects of geometry and kinematics. It can simulate heterogeneous geological features using the mechanical properties of different materials in the calculation. Moreover, the elements of the subdivided geological body are interconnected, and the number of elements can be set as required. The more divided the elements are, the closer the model is to the actual geological body. The whole body moves under the boundary stress, which significantly improves the accuracy of model operation. The finite element numerical simulation of the stress field includes three major steps: establishing a geological model, constructing a mechanical model, and analysing boundary conditions.

Establishment of Geological Model
Establishing a geological model is the prerequisite and basis for numerical simulation. On the one hand, the geological model determines the framework of the mechanical model. On the other, the deviation between the geological model and the actual geological conditions determines the simulation accuracy, which is critical to the correct analysis of the in-situ stress distribution pattern.
The establishment of a geological model needs to follow certain principles: 1) When selecting the target horizon of the work area, the range of the work area shall be extended outward as far as possible, and the real tectonic relief in the work area should be reflected as much as possible based on the tectonic framework. At the same time, it is also necessary to refer to the distribution of sedimentary facies in the work area and divide the geological elements under the joint constraints of facies boundaries and faults. 2) Based on the regional stress state, the exterior of the geological body in the work area is wrapped with a regular cube, whose surface planes facing three directions are perpendicular to the three principal stress directions in the area, respectively. 3) Since numerous faults of different sizes are developed in the work area, displaying these faults one by one would generate excessive workload and also cause errors in software processing. Therefore, the faults need to be screened. The ones whose length is larger than 1 km that influence and reflect the tectonic pattern of the study area are retained, and the minor ones whose length is less than 1 km are ignored.
The geological model in this study was established based on the three-dimensional present-day tectonic framework of the Jia 2 member in the Puguang area. According to the single-well data of the Jia 2 member, the thickness of this layer ranges 146-293 m, with an average of 220 m; hence, the model thickness was set to 220 m.

Establishment of Mechanical Model
Based on the geological model, the mechanical elements are divided and assigned with different rock mechanical parameters according to the distribution of the parameters (mainly the elastic modulus and Poisson's ratio) that affect the state of the stress field across the layer. The rock mechanical parameters can be acquired using the static and dynamic methods. Since the parameters obtained by the static method are discrete and costly, the core data of the Jia 2 member in the Puguang area are limited, whereas the logging data are relatively abundant. Therefore, we mainly used the logging data to calculate the rock mechanical parameters in this study (Barree et al., 2009).
where E is the elastic modulus, GPa; μ is the Poisson's ratio, dimensionless; ρ is the rock density, g/cm 3 ; t s is the transverse wave time difference, μs/ft; and t p is the longitudinal wave time difference, μs/ft. After fitting the single-well transverse waves with the multiple regression method, we calculated the rock mechanical parameters for 27 wells and established an equivalent mechanical model for a single well based on the litho-stratification data in the vertical direction. Among the rock mechanical parameters, the response of the rock's elastic modulus to the in-situ stress distribution is more pronounced than the Poisson's ratio [24] . Therefore, considering the calculation accuracy and the computing power of the workstation, the geological model was divided into different mechanical elements based on the criteria of elastic modulus: E ≥ 54 GPa, 54GPa > E ≥ 50 GPa, 50GPa > E ≥ 46 GPa, and 46GPa > E. To eliminate the influence of boundary effects, we confined the exterior of the modelled region with a rectangular frame as a loading and constraining object. The rectangular frame is a separate mechanical element whose rock mechanical parameters are determined based on the mean values of the internal geological body, and the corresponding elastic modulus and Poisson's ratio are assigned to each mechanical element.
The treatment of faults is a critical issue. As the samples from fault zones are extremely loose and fragmented, their mechanical parameters cannot be obtained from rock mechanics experiments. According to the theory of tectonic deformation, the fault fracture zone, during its formation, is the dominant stress release zone of the entire geological body. After formation, it is the weakest part of the geological body. In previous studies, 40-70% Young's modulus of the standard rock mass and a Poisson's ratio with an increase of 0.02 were typically taken as the mechanical parameters of the fault zones. The more complex the tectonics of the study area, the more developed the faults, and the smaller the Young's modulus and the larger the Poisson's ratio of the fault zone (Chen and Tian, 1998;Wang et al., 1999;Zu et al., 2018). On this basis, the mechanical model is subdivided into a grid. The geological body is represented by a finite number of elements interconnected by nodes and lines and transformed into a mathematical model for calculation. The element subdivision needs to be performed according to the geological characteristics; therefore, the same element cannot cross two different geological bodies. Because of if the acute angles is less than 5°, the software will not able to calculate the stress of the element around the corner. So, particularly small acute angles (<5°) in the modelled geological bodies and elements should be avoided. In this study, the strata of the Jia 2 member in the Puguang area were subdivided into 2.78 million elements.

Boundary Condition Analysis
For in-situ stress simulation, the boundary conditions mainly include stress boundary conditions and displacement boundary conditions. Since it is impossible to obtain accurate and reliable data on strain rate induced by the present-day in-situ stress state, we selected the stress boundary conditions as the boundary constraints of the finite element modelling in this study, i.e., the boundary conditions were calculated by loading the stress.
The present-day in-situ stress direction was obtained by analysing the orientation of drilling-induced fractures and the wellbore collapse direction from cross-dipole acoustic logs and dip meter logs of 30 wells in the Puguang area. The results show that the present-day maximum horizontal principal stress orientation is near E-W in the Jia 2 member. The present-day in-situ stress magnitude is calculated based on conventional logging data. According to the structural and geological characteristics of the Puguang area, we employed Huang's model and the combined spring model to calculate the principal horizontal stress in 22 wells and obtained the average gradient of principal stress magnitude varying with depth in the Jia 2 member. The maximum horizontal principal stress σ 1 had an average gradient of 0.0166 MPa/m, with an E-W orientation; and the minimum horizontal principal stress σ 2 had an average gradient of 0.0115 MPa/m, with a N-S orientation. In addition, the gravity applied in the vertical direction and the constrained bottom were also boundary conditions for the loading calculations (Figure 2).

RESULTS OF STRESS FIELD NUMERICAL SIMULATION Error Analysis
The accuracy of model calculation results is evaluated by comparative analysis with actual data. The simulated in-situ stress direction was tested with the maximum horizontal principal stress direction calculated based on the wellbore collapse and drilling-induced fracture information from logging data interpretation ( Table 1). The average error in the overall in-situ stress direction is 6.25°, and the maximum error in a single well does not exceed 15°. The simulated maximum and minimum principal stress magnitudes were tested with the in-situ stress magnitude calculated from the single-well data ( Table 2). The overall maximum and minimum principal stress magnitudes have an average error of 7.29% and 8%, respectively, and the single-well error does not exceed 20%, indicating that the simulation results are reliable.

Stress Field Distribution Pattern
The simulation results reveal that the present-day in-situ stress state in the Jia 2 member of the Puguang area is dominated by E-W compressive stress (Figure 3). The minimum horizontal principal stress in the Jia 2 member is mainly concentrated between 30 and 60 MPa, the maximum horizontal principal stress is generally between 50 and 80 MPa (Figure 4), and the differential stress predominantly lies between 10 and 40 MPa ( Figure 5). The in-situ stress distribution is remarkably controlled by the burial depth and lithology, and the overall distribution pattern of the maximum and minimum horizontal principal stresses is generally consistent. The high-value areas of principal stress magnitude are concentrated in the north-eastern and western Puguang area, where the tectonic deformation is less intense and the tectonic relief is gentle. The maximum and minimum horizontal principal stresses are small in the structurally high part of the central area.

DISCUSSION
The study and evaluation of in-situ stress are carried out throughout the exploration and development of oil and gas fields. Especially during oil and gas field development, the insitu stress has a significant impact on the well pattern deployment, water injection development design, perforation scheme planning, and hydraulic fracturing treatment, as well as the later-stage casing deformation analysis and wellbore stability evaluation (Wan et al., 2009). In water injection Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 839084 5 development, it is critical to determine the fracture opening pressure. In particular, the fracture opening pressure and formation fracturing pressure in the areas with multiple fracture systems should be comprehensively considered ; their magnitudes are closely related to the in-situ stress. In this study, the opening pressure of the fractures with different orientations and the formation fracturing pressure were calculated based on the simulation results of the present-day stress field to determine the appropriate fracture opening sequence.

The Angles Between the Maximum Principle
In-Situ Stress and the Azimuth of Nature Fracture As Figure 3 shows, the in-situ stress orientation is locally deflected due to the influence of faults, and the effects of faults on the in-situ stress state are particularly prominent in the case of similar rock mechanical parameters (Zu et al., 2014). The fault zones are generally shattered zones manifested as weak geological bodies. Under the same regional compressive stress, the fault zone absorbs the stress and causes the deflection of the in-situ stress direction in surrounding areas. The trend of deflection gradually converges to the fault strike, and the deflection of in-situ stress direction is more significant at the fault tips than in the middle section of the fault, with a maximum local deflection of 20-30°. Because the state of the in-situ stress indeed varies between the elements separated by the fault, as a result, the defection pattern varies. When looking from the centre to any of the tip of the fault, on the left side around the fault tips, the orientation of the principle of the in-situ stress rotates clockwise, and on the opposite side, i.e., the right side, the orientation of the principle of the in-situ stress rotates anti-clockwise. In areas away from the fault, the degree of deflection gradually decreases. In addition, the scale of the fault and the angle between the fault strike and the regional maximum horizontal principal stress also affect the degree of present-day in-situ stress field deflection in the vicinity of the fault zone.
The strike, dip, and dipping direction of the tectonic fractures in the Jia 2 member of the Puguang area were counted by imaging logging, and the statistics ( Figure 6) show that there are mainly three sets of tectonic fractures in the Jia 2 member, trending E-W, NW-SE, and N-W, respectively ( Figure 6A). Among them, the near E-W fractures are high-angle fractures dipping towards due north or due south. Their dip is concentrated between 69.81°and 89.08°, with an average of 80°. The fractures with near NW-SE orientations dip towards NE or SW, with a value of 25-48°and an average of 40°. The near N-S trending fractures dip towards W, with a value of 20-40°and an average of 30°( Figure 6B).

The Magnitude of Opening Pressure and the Formation Fracturing Pressure
The opening pressure of each fracture system and the formation fracturing pressure can be calculated using the three-dimensional present-day stress field model, and the opening pressure of the fractures with different orientations can be expressed as (Zeng, 2004): Hρ s g sin θ + Hρ s gcosθ − Hρ o g + σ H sin θ sin β where P i is the fracture opening pressure, MPa; μ is the Poisson's ratio of the rock, dimensionless; H is the burial depth of the  fracture, m; ρ s is the rock density, kg/m 3 ; g is the gravitational acceleration, N/kg; θ is the dip of the fracture (°); β is the angle between the present-day maximum horizontal principal stress and the fracture strike (°); ρ o is the pore fluid density, kg/m 3 ; σ H is the maximum horizontal principal stress, MPa; and σ h is the minimum horizontal principal stress, MPa. The formation fracturing pressure is : P f is the formation fracturing pressure, P o is the pore fluid pressure, and σ T is the tensile strength of the rock.
The calculation results show that the opening pressure of the fractures with near E-W, near N-S, and NW-SE orientations is concentrated in the ranges of 33-49, 78-107, and 84-117 MPa, respectively. In comparison, the formation fracturing pressure lies between 41 and 64 MPa, higher than the opening pressure of the near E-W fractures but lower than that of the other two sets of fractures. Therefore, during the hydraulic fracturing operation, the E-W fractures will open first and serve as the main seepage channels. As the fracturing pressure gradually increases, the formation will begin to rupture, producing a fracture network. A contour map of the fracture opening pressure in the Puguang area can be plotted based on the calculated opening pressures of different fracture systems (Figure 7). Designing the appropriate water injection pressure based on the fracture opening pressure in different areas is an effective measure to enhance oil and gas recovery.

CONCLUSION
1) The numerical simulation results show that the present-day in-situ stress state in the Jia 2 member of the Puguang area is dominated by E-W compressive stress, and the in-situ stress orientation is locally deflected due to the influence of faults. 2) The minimum horizontal principal stress, maximum horizontal principal stress, and differential stress in the study area are mainly concentrated in the ranges of 30-60, 50-80, and 10-40 MPa, respectively, and the in-situ stress distribution is remarkably controlled by the burial depth and lithology. 3) Among the opening pressures of the fractures with different orientations in the Jia 2 member, calculated using the threedimensional heterogeneous present-day in-situ stress model, the opening pressure of the near E-W fractures is the smallest, concentrated between 33 and 49 MPa, which is lower than the formation fracturing pressure. During hydraulic fracturing, the E-W fractures will open first and become the main seepage channels. The continuously increasing operating pressure will lead to formation breakdown, producing a fracture network.  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 839084 8