Numerical simulation of the stress field and fault sealing of complex fault combinations in Changning area, Southern Sichuan Basin, China

In order to quantitatively study the impact of faults on shale gas productivity, taking Changning area in the south of Sichuan Basin as an example, this study analyzes the impact of fault scale and strike on shale gas productivity under complex structural environment. The typical well profile is selected to establish a three‐dimensional numerical model, and the structural stress distribution characteristics of common structural combinations in this area are analyzed by finite‐element numerical simulation. If, IFT, and IC parameters are selected to quantitatively calculate the sealing of faults. It is found that NE‐trending faults and faults with fault drop grade of 1 in the study area have the greatest damage to shale preservation. The combination of wide and gentle syncline and small fault, low‐angle anticline, and large fault is a typical structural style in this area. The structural combination styles of well N216 and well N217 are similar, but well N216 is subjected to large compressive stress. The calculation shows that the fault plugging performance of well N216 is better than that of well N217.


| INTRODUCTION
The exploration and development of shale gas are currently key aspects of oil and gas exploration in the world. China has developed rapidly in this field and has successfully realized commercial exploitation in the Sichuan Basin. By 2020, the annual output of shale gas exceeded 200 × 10 8 m 3 . However, this output still represents only 60% of the expected production capacity because southern China is affected by superposed multistage tectonic movements, the geological conditions are complex, and faults are widely developed. [1][2][3][4][5] After early hydrocarbon generation, the shale gas in a reservoir escapes to varying degrees during later evolution. Different from the relatively simple burial conditions in the United States, the shale reservoir of the Wufeng Longmaxi Formation in the Sichuan Basin usually has a complex structural type involving folds and various types of fault combinations, as shown by drilling. [6][7][8][9] The sealing property of faults developed in shale reservoirs is one of the key factors determining whether the formed reservoirs can be preserved in an area with complex geological conditions such as the Sichuan Basin, which shows the effects of superimposition of multistage tectonic movements. [10][11][12] In situ stress is a natural stress existing in the crust where no engineering disturbance has occurred and is mainly composed of a self-weight stress field and a tectonic stress field. All geological structures such as joints, faults, and folds indicate crustal rock deformation and fractures caused by tectonic stress. [13][14][15] At the same time, various geological structures also have different effects on the distribution pattern of the crustal stress field. Studying the influence of different types of geological structures on the distribution pattern of the crustal stress field is very relevant for understanding the state of the original stress field of shale reservoirs against the background of complex structures in South China. The study of fold deformation and fault sealing is one of the important topics in the field of oil and gas exploration and development. [16][17][18][19][20] There is a very close relationship between in situ stress and fault sealing. Previous research has led to a deep understanding of the qualitative relationship between in situ stress and fault sealing. In an environment with strong compressive stress, the pores in the rock mass become smaller under pressure, and the sealing performance of a section is good, which creates a caprock for oil and gas. In contrast, when a section is in a tensile stress environment, the sealing performance of the section is poor, and a fault opens, becoming a channel for oil and gas escape. [21][22][23][24][25][26][27][28][29] In the quantitative evaluation of fault sealing, obtaining the stress value on the fault plane is the key to solving the problem. At present, the most intuitive and accurate method to obtain the stress value is to measure the crustal stress by field hydraulic fracturing tests, but this method is restricted by the location and number of well points, and the measured in situ stress is not the normal stress on the fault plane. The numerical simulation of tectonic stress field based on finite-element analysis can better solve this problem, and can inverse the stress field of different tectonic elements in different periods. [30][31][32][33][34][35][36][37][38][39][40] By using the results of structural interpretation, establishing the corresponding geological model, and carrying out numerical simulation calculations of the tectonic stress field, the tectonic stress field of specific structures and the normal and shear stresses on the fault plane can be obtained, and the characteristics of stress on a section in different periods can be reproduced. [41][42][43][44] On this basis, through the comparison of section normal and shear stresses and section rock compressive and shear strengths, the quantitative study of fault sealing has the advantages of simplicity, effectiveness, and easy visualization. This method can not only investigate the sealing properties of faults in different geological historical periods but also quantitatively evaluate the sealing properties of faults in different intervals and horizons. [45][46][47][48][49][50] Based on the analysis of fault characteristics and existing data on drilling productivity in the study area, this study analyzes the impact of different grades of faults on shale gas productivity, determines the typical structural combinations of high-yield wells in the area, and constructs a corresponding geological model by using the interpreted results of seismic profiles that pass through wells; according to the finite-element method of numerical analysis, the distribution characteristics of the structural stress field for a complex structural combination are analyzed. According to the calculated values of the normal stress and shear stress on a section, the sealing strengths of different faults in different layers are quantitatively evaluated.

| GEOLOGICAL BACKGROUND
The study area is located on the southern margin of the Sichuan Basin, extending to southern Gaoxian County in the west and reaching the Xuyong area in the east. This area is structurally located in the transitional region between the fold belt in southern Sichuan and the fault-fold belt in Loushan. The Changning anticline is a secondary structural unit in this area, with a WNW-ESE anticline axis, steep northeastern limb, and gentle southwestern wing. The anticline is south of the Loushan fold belt, northeast of the Sichuan high and steep fold belt, and adjacent to the Shuanglong-Luochang syncline to the west. The structural development in the study area is complex, with many faults and nose structures. The faults in the northeastern part of the study area have long extents and large fault throws, while small faults are mainly developed in the south and west; the main fault strikes are NE and NNW, which are quite different from the fold strike. 10,51,52 The faults in the region are superimposed under the influence of tectonic stresses with multiple stages and directions; this situation has led to multistage and multiform fault and fold combinations, which have great importance for the enrichment of shale gas and improvements in reservoir performance and productivity ( Figure 1).

| SAMPLES AND METHODS
This study is mainly carried out by the combination of fine seismic interpretation and three-dimensional finiteelement numerical simulation. According to the interpretation results of crosswell seismic profile in the study area, the effects of different structural combinations on shale gas enrichment and preservation are analyzed. Based on the seismic interpretation results, a three-dimensional geological model is established, and the influence of complex structural combinations on the distribution of in situ stress field is analyzed. Calculate the normal stress and shear stress on the fault plane near N216 and N217 wells, and select three parameters: fault sealing coefficient (I f ), fault tightness index (I FT ), and fault shear index (I C ) to quantitatively characterize the fault sealing, and analyze the control effect of fault sealing on productivity. 53

| Fault sealing coefficient
The fault sealing coefficient (I f ) refers to the ratio of normal stress and fluid pressure on the fault plane, namely, where σ is the normal stress on the section, f is the anomalous pressure coefficient of the formation, ρ w is the density of water, which is 1.0 g/cm 3 , g is the gravity acceleration of 9.8 m/s 2 , and Z is the stratum depth.
When I f > 1, the fault is closed, and the greater the I f value is, the higher the degree of sealing. When I f ≤ 1, a fault opens and becomes a passage for oil and gas migration, and the smaller the I f value is, the higher the degree of opening.

| Fault tightness coefficient
Predecessors have made qualitative research on the relationship between the stress state of the section and the compressive strength of the section material. The research shows that when the compressive stress of the section is greater than the compressive strength of the rock mass, the cracks in the rock mass are compressed, the pores are reduced, and the sealing is enhanced [54][55][56] ; and others have quantitatively characterized this process and proposed the concept of the fault closure index. The fault closure index (I FT ) can be expressed as the ratio of the normal stress σ on the cross section to the compressive strength σ P of the rock in the fault zone, namely, Surface geological map of the study area. 10

| Fault shear index
The fault shear index (I C ) can be defined as the ratio of the shear stress τ on the fault plane to the shear strength σ C of the fault zone material 55 : Here, σ C is the sum of the intrinsic shear strength of the fault zone material and the friction force on the fault plane, that is, where c is the inherent shear strength in MPa and ω is the internal friction coefficient.
The significance of the fault shear index is that when I C > 1, shear slip will occur on both sides of the fault, which will cause a change in the physical properties of the fault zone. If I C > 1 and the normal stress on the fault plane is compressive, the two plates of the fault plane will undergo compressive shear dislocation, which will compress the material in the fault zone. The original oil and gas seepage channels such as pores and fractures are closed under shear, which strengthens the sealing degree of faults and improves the ability to block oil and gas. In contrast, when the normal stress is tensile, tensile shear dislocation occurs on both sides of the fault, and the scale of original oil and gas seepage channels such as pores and fractures expands under tensile action. At the same time, due to the shear dislocation of both sides, new fractures can be created in the fault zone, the fault opening degree can be improved, and faults can become favorable channels for oil and gas migration. However, when I C ≤ 1, the fault undergoes only shear dislocation, the two sides of the fault will not produce relative displacement, and the physical properties of the fault zone will not change obviously; thus, the influence on fault sealing can be ignored.

| Distribution characteristics of faults
After the Indosinian movement, the study area was mainly affected by three tectonic movements accompanied by compressive stress in three directions, and a fold and fracture system finally formed that is oriented mainly in the NE and NW directions, with a few oriented nearly E-W. NE-trending faults mainly formed from the late Yanshanian to early Himalayan movements, and NW-trending planar conjugate shear fractures and NEtrending sectional shear fractures formed by NW-trending compression in the study area. Among these fractures, the shear fractures in the NE-trending section further expanded and connected, forming many NE-trending faults. NE-trending faults are extremely well developed in well block N216, with dip angles greater than 60° and long extents. NW and nearly E-W faults are less developed, and the lengths of the faults are small. From well block N201 to well block N209, there are few small-scale faults, and the shale reservoir is little affected by structural activities. The formations in this area are relatively stable (Figure 2).

| Fault scale and productivity
Faults of different scales have various influences on shale gas productivity. Given the lengths and throws of faults, the faults in the Sichuan Basin and its periphery can be divided into four grades according to scale. The first grade is the basin-controlling fault, which has a large scale, long extent, great depth, and early development, and plays an integral role in controlling the tectonic evolution and sedimentation of the basin. The first grade is represented by the Huayingshan fault and Qiyueshan fault. The second grade is the zone-controlling fault, which governs the distribution of structural belts, has a large scale, extends a long distance, and cuts some strata. Typical second-grade faults are the Taihechang structural belt and Lintanchang structural belt. The third grade is local structure-controlling fractures, and the fault throw is generally 100-500 m. The fourth grade is a small fault in a local structure, and the fault throw is usually <100 m. There is currently no uniform standard for fault classification, and the classification schemes vary according to the characteristics of different research areas and diverse research purposes.
According (4) C = C + is large. This area has a gentle syncline structure, with good overall preservation conditions and high productivity. Through statistical analysis of the data, when the gas well location is close to the large-scale NE-trending fault, the productivity is greatly affected (as shown in the table below), and the productivity of adjacent gas wells is low, or no gas is produced ( Figure 4A). Among them, only wells N216 and N217 have a production capacity of more than 100,000 m 3 /day ( Figure 4B). The reasons may be that the sealing performance of the fault near the well locations is good and that the shale gas around the well is well preserved.

| Typical fault-fold combination model
Good preservation conditions are the key to "reservoir formation and production control" of shale gas, and shale gas will accumulate in stable areas with weak tectonic movement and poorly developed faults. The southern Sichuan Basin has experienced transformations involving multistage and multicycle tectonic movements. Different stages and intensities of tectonic movements have resulted in different strata fold deformation, fracture degree, and denudation degree. A variety of fault combination patterns have formed, which have greatly damaged the preservation conditions of reservoirs. However, in areas with strong tectonic movement and well-developed faults, faults may also become closed zones, and shale gas in The syncline structure located in the middle of the study area is a stable deeply buried feature with relatively wide and gentle strata, little fault development, and great burial depth; the shale reservoir is only slightly damaged. Shale gas has the characteristics of retention in belts and enrichment in the center of the syncline.
Taking well N201 (the tested output is 150,000 m 3 /day) as an example, the structural style is gentle and stable, the well is far away from the denudation area, and the degree of fault development is low ( Figure 5). The fault throws of F46 and F54 near the N201 well are less than 200 m, and the faults do not penetrate the upper overburden; these faults are far from the well site at a distance of more than 4 km; thus, the degree of damage to shale preservation conditions is low. Wide and gentle syncline gas reservoirs easily form A-shaped natural fractures with favorable physical properties.

| Low-angle anticline + large fault
The western part of the study area is characterized by strong structural deformation and large-scale and dense fault development. Taking well N216 (test production 20.29 × 10 4 m 3 /day) and well N217 (test production 11.12 × 10 4 m 3 /day) as examples, the Longmaxi Formation at the core of the anticline is relatively flat with a low dip angle. Large-angle torsion occurs only near the turning ends of both limbs where faults and cracks are likely to develop.

F I G U R E 7
Interpretation results of C-C' seismic profile combinations encountered are similar, but the production gap is large. According to the field hydraulic fracturing test data, the current maximum principal stress in well N216 is 64.8 MPa, which is greater than that in well N217 at 47.6 MPa ( Table 1). The normal stress and shear stress on the fault planes near the two wells are different, the fault sealing capacities also differ, and the preservation conditions are diverse.

| Finite-element numerical simulation
Obtaining the normal stress of the fault plane is the key to solving the problem. The typical sections of wells N216 and N217 are selected for three-dimensional modeling analysis, and the normal and shear stresses on the fault plane are obtained by the numerical simulation method.

| Geological model
According to the results of seismic interpretation on the section and the planar distribution characteristics of faults, the main large-scale faults that have developed near the well location are selected to reasonably simplify the geological model. The model is a multilayer structure, and the rock mechanics parameters of Longmaxi formation are obtained by experiment. The mechanical parameters of other strata and fault zones are determined according to the manual of rock mechanics and relevant literature, combined with the rock mass characteristics of the study area. [57][58][59] The rock mechanical parameters of each stratum are in the Table 2. [60][61][62] The main body in this study is the Longmaxi Formation strata and faults, which have highly accurate grid division. The remaining formation model grids are rougher. The model of the 216 profile is divided into 269,105 units and 50,736 nodes ( Figure 8A). The geological model of profile 217 is divided into 236,776 units and 44,240 nodes ( Figure 8B).

| Model boundary conditions and load conditions
When simulating the evolution of the stress field of the composite structure, the approach is divided into two steps that are carried out successively. After several trial calculations, when the simulated calculation value at the well location is close to the measured value, the T A B L E 1  [62][63][64][65] Step I simulates the self-weight stress field in the stage of pressure from the overlying strata. A vertical stress σ V is applied to the top of the model as a load, and the vertical displacement of the bottom surface of the model is fixed. The normal horizontal displacements of the four sides are fixed.
Step II simulates the tectonic stress field imposed by the horizontal tectonic stress. At this stage, the fixed right side of the model is released, and the horizontal tectonic stress σ2 is applied as the load. The vertical displacements of the bottom surface and top surface of the model are fixed, and the normal horizontal displacements of the other three side surfaces are fixed. The measured maximum horizontal principal stress consists of two parts, as shown in Equations (5) and (6).
In the formula, σ H is the maximum measured horizontal principal stress, σ 1 is the horizontal stress component generated by the self-weight stress field; σ 2 is the horizontal stress component generated by the tectonic stress field, μ is Poisson's ratio, and σ V is the vertical stress.  Figures 9A and 10A). The reliability of the model is high. The maximum principal stress on the fault plane is horizontal, and the minimum principal stress is vertical, which is in the stress state of a reverse fault and is consistent with the actual geological exploration results ( Figures 9B and 10B).
The maximum principal stress generally increases with increasing depth, forming a stress concentration area on one side of the fault end and a stress reduction area on the other side. The stress direction at the end is variable. The normal stresses on the F11 and F13 fault planes near well N216 are 60.2 and 73.4 MPa, respectively. The average normal stresses on the F21, F20, and F54 fault planes near well N217 are 43, 41, and 37 MPa, respectively.
In the composite structure in which the fault truncates one limb of the anticline, the fault (excluding the end) replaces the influence of the anticline axis so that the stress concentration originally seen in the anticline axis is transferred to the fault and its vicinity, and the anticline axis no longer presents an obvious stress concentration area.

| Distribution of the shear stress
The shear stresses in wells N216 and N217 are affected by the structural shapes; the shear stress on the whole fault plane is lower than those in the strata on both sides, and there is a shear stress concentration area at the end of the fault. In the strata fault block near the inner side of the fault end, the shear stress is low. The average shear stresses on the F11 and F13 faults near well N216 in the Longmaxi Formation are 19.83 and 16.97 MPa, respectively ( Figure 11). The average shear stresses on the F21, F54, and F20 fault planes near well N217 are 9.7, 9.6, and 10.2 MPa, respectively ( Figure 12).

| Calculation of section closure
The normal and shear stresses on the section can be analyzed by finite-element numerical simulation based on the actual structural style. According to the normal and shear stresses on the section and the comparison of rock compressive and shear strength in the section, a study of fault sealing can be carried out, which has the advantages of simplicity, effectiveness, and easy visualization. 55,56,[66][67][68][69][70] According to the normal and shear stresses on the section calculated by numerical simulation, the section sealing capacity is comprehensively analyzed by using three quantitative evaluation parameters: I f , I FT , and I C ( the gas wells are more than 2 km away from the grade 1 fault, and the impact of the fault on gas well productivity is not obvious. 2. In the study area, the gas wells are mainly located in the N201 and N209 well areas. The main structural style is the combination of small faults and a wide and gentle syncline, and the overall preservation conditions are good. Faults of different scales are densely developed in the N216 well area. The structural style of typical wells N216 and N217 with good productivity is the combination of large-scale reverse faults and low-angle anticlines, forming a horst structural style. The nearby fault is a reverse fault, which is in a state of compressive stress with good fault sealing, forming a trapped gas reservoir. In the complex structural area where faults are relatively well developed, the horst structure combination that formed under structural compression is beneficial to the preservation of shale. 3. For the same structural style, the greater the compressive in situ stress is, the better the fault sealing performance, the less the amount of shale gas that escapes, and the higher the productivity. The strata of the Longmaxi Formation bear a large tectonic compressive stress, and the maximum principal stress is mainly horizontal. The normal stresses on the F11 and F13 fault planes in the N216 well are 60 and 71 MPa, and the shear stresses are 18.3 and 21.7 Mpa, respectively. The normal stresses on the F21, F20, and F54 fault planes in well N217 are 43, 41, and 37 MPa, respectively, and the shear stresses are 6.7 and 9.2 MPa. According to the quantitative results of calculations using the selected section sealing parameters, the F11 and F13 faults sealing near the N216 well location are strong, and the degree of damage to the shale reservoir is low. This is the key factor for the high productivity of well N216.