Study on induced stress of hydraulic fracturing in fractured‐porous elastic media

The study on induced stress of hydraulic fracturing is an important part of hydraulic fracturing design, temporary plugging and steering fracturing design, and completion design in oil and gas development. Based on fracture continuum method (FCM), cohesive unit method, and finite volume method (FVM), the propagation behavior of multiple hydraulic fractures (HFs) in fractured‐porous media reservoir is simulated from the perspective of fluid and geological stress coupling. The influence of natural fractures (NFs) on stress distribution and stress inversion sensitive factors are studied. NFs in the model are equivalent to physical models controlled by tension and shear bonds. The proposed model can deal with any number of NFs, which greatly improves the computational efficiency of the model. Then PKN analytic solution and Abaqus software results were used to verify the correctness of our model. The results show that the NF will make the reservoir stress show strong heterogeneity after fracturing, and the prominent feature is the regional stress mutation. The sudden change of stress shows the sudden increase or decrease in stress; with the increase in the initial stress difference, the stress inversion becomes more difficult. When the stress difference exceeds 3 MPa, the stress inversion area is almost 0. The stress interference increases with the decrease in cluster spacing. The stress inversion region between HFs decreases, while the external stress inversion region increases; Biot coefficient has great influence on stress and stress inversion. With the decrease in Biot coefficient, stress inversion region will decrease obviously. When the Biot coefficient is 0.4, there is no stress inversion area in the whole reservoir. Poisson's ratio has little influence on stress inversion. When Poisson's ratio changes between 0.1 and 0.3, the change range of stress is within 1 MPa, and the stress inversion area is almost unchanged. Our results can be helpful for understanding the stress distribution after hydraulic fracturing in fractured‐porous media.


| INTRODUCTION
With the rapid consumption of conventional oil and gas resources, unconventional oil and gas resources are more and more concerned by oilfield developers. However, because unconventional reservoirs are generally characterized by low porosity, low permeability, and NF development, the recovery rate of such reservoirs is generally not high. In order to improve unconventional oil and gas recovery, horizontal wellstaged multicluster fracturing is often used in oil fields. 1,2 Reservoir reconstruction is a complicated process, which is different from the conventional hydraulic fracturing process. For reservoirs where NFs are well developed, such as shale, not only do HFs expand after fracturing, but NFs also break under fluid pressure. This, in turn, leads to a significant increase in tight reservoir permeability, which promotes fracturing fluids loss. For the simulation of hydraulic fracturing behavior in unconventional reservoirs, a variety of factors, such as NFs, fluid filtration, reservoir deformation, temperature change, and damage, need to be considered. 3 In order to accurately understand the mechanics, seepage and coupling mechanism in the fracturing process, these factors should be considered as much as possible in the modeling process.
At present, most scholars are keen on the simulation of dynamic fracture extension behavior, which mainly includes two aspects: general nonplanar fracture extension and complex nonplanar extension of fractures. Generally, nonplanar fracture propagation refers to a computational model that determines the direction of fracture propagation according to the physical state of the fracture tip at each extension step and the preselected criteria of rock fracture and fracture propagation. Under this computational model, cracks are not confined to a standard plane, so there is no analytical solution at present. For the problem of nonplanar crack propagation, many scholars have developed a variety of numerical methods to match it, mainly including finite element method, extended finite element method (EFEM), displacement discontinuity method (DDM) derived from boundary element method, discrete element method (DEM), FVM, and phase field method (PFM). Belytschko et al, 4,5 Olson, 6 Wu et al, 7,8 Shi et al 9,10 and Zhou et al 11,12 have more prominent studies in these fields. Guo et al 13 established a 3D model of multicrack extension based on FEM. This model takes into account the stress interference effect between HFs, which leads to the limitation of fracture length opening, but does not consider its nonplanar extension behavior. Ren et al 14 established a multicluster fracture extension model based on the pseudo-three-dimensional model and used the discontinuous displacement method in the boundary element method to solve the coupling problem. The model takes into account the nonuniform extension of each hydrofracture. Shi et al 9,10 and Wang et al 15 established a multicrack synchronous expansion model based on the XFEM and calculated the fracture extension step length with implicit tip iterative progressive solution. The influence of formation stress conditions on the turning extension of HF is also studied. Based on discrete element model and discrete element software, some scholars have established two-dimensional and three-dimensional HF extension models. [16][17][18] By simulating the elastic and cementing action between rocksolid particles and the action of coupling pore fluid on solid particles, Sarmadivaleh 16 proposed a two-dimensional discrete element model to calculate the process of microscopic failure between cemented particles to the final formation of macroscopic cracks. Based on the microscopic characteristics of rock particles, the models simulate the macroscopic fracture extension behavior from the physical mechanism, so it can simulate the complex behavior of multiple three-dimensional fracture turning extension. 19,20 The complex nonplanar fracture model is further developed based on nonplanar fracture model. This model is no longer limited to the study of the nonplanar extension of HFs. Its core is to deal with the intersection and extension behavior of many HFs and hundreds of NFs. 21 As the number of natural cracks increases exponentially, the division and processing of model grids and the optimization and upgrading of numerical calculation methods are all major challenges for complex nonplanar fracture models. As the frontier of fracture extension, complex nonplanar fracture model has become the focus of research in recent years. Based on the particle discrete element method, Damjanac et al 17,18 regard the NF as an ellipse of different areas composed of particles embedded in the reservoir. This model can well consider the coupling between the fluid inside and outside the HF and the deformation of the reservoir, but the main disadvantage is that the extension direction of the HF is preset, and the extension mode after the intersection is characterized by crossing, while the discrete NF does not have the expansion performance. Based on the finite element method, Zou et al 22 and Li and Zhang 23 set some matrix grids to represent the NFs and distinguished them by defining that the NF grids and matrix grids have different failure conditions. This method can accurately represent the flow of fluid in the main HF and NF, but it cannot realize the coupling of fluid flow and reservoir deformation in all fractures, and does not consider the impact of fluid leakage, and can only approximate the arbitrary direction of fracture propagation. Wu et al 8 and Olson and Wu 24 established a model of intersection and expansion of multiple clusters of HFs and NFs by using the simplified three-dimensional displacement discontinuity method, which can truly simulate the formation of fracture network under the interference of interwell fracture extension. In the model, the stress interference during the extension of multiple fractures is considered, which can optimize well spacing and accurately realize the extension of fractures in any direction under the influence of stress shadow. However, the coupling between fracture leakage fluid and reservoir deformation is not considered in the model. Based on the XFEM method, Shi et al 10 simulated the intersection and extension law of HFs and NFs, as well as the fracture network expansion behavior in fractured reservoirs. Compared with conventional FEM, XFEM has higher computational efficiency and stronger grid flexibility, so it is more suitable for the simulation of fracture network expansion behavior. Based on DEM method and PFC software, Yoon et al 25 and Zhang et al 26 characterized the NF structure with "Smooth Joint," constructed a digitalized NF rock mass, and simulated the planar expansion behavior of HFs in it. The influence of NF angle, mechanical properties and other parameters on the intersection, and propagation of HF and NF are studied and analyzed. Due to the huge amount of calculation involved in DEM method, this model is only used for the expansion simulation of small-size fracture network, such as the simulation of indoor core fracturing experiment.
Through the analysis of previous studies, it is found that most of the current studies focus on the extension behavior of HFs under the influence of NFs and joints in nonintersecting or intersecting state. However, due to the limitation of the solution method, all the characteristics of HF, NF, and matrix cannot be fully considered. For example, when the number of NF lines increases sharply, the calculation time for solving the process of complex fracture extension under the condition of coupled fluid and stress is prone to nonconvergence and exponential increase in the calculation time. 27,28 Therefore, most of the current studies on complex fracture network can only be applied to mechanical research, and it is difficult to apply to engineering. In addition, when fluid flows in HF and NF, the fluid filtration is generally treated as the filtration loss is zero or the filtration loss is calculated through the filtration coefficient. This treatment is inconsistent with the actual situation. In most hydraulic fracturing simulations, the reservoir is considered as an elastomer rather than a porous elastic medium. The stress obtained under this condition ignores the induced stress caused by pore pressure. In fact, there is a pressure gradient between the fluid in the fracture and the fluid in the matrix. Increased fluid pressure in the matrix creates a "back stress" that tends to close the fracture. 29 This is also ignored in most studies. In addition, in many studies, most stress nebulae of complex fracture networks are not given, and the authenticity of coupling is questioned. Salimzadeh and Khalili 30,31 did a better job of this, and they proposed an XFEM model that included two independent flow models in the fracture and the rock matrix, with a leak-off mass transfer between fracture and rock matrix to link the two. The leak-off depends on the pressure gradient in the matrix adjacent to the fracture, as well as on the fluid viscosity and matrix permeability. Norbeck et al, 32 using an embedded fracture model, also considered two flow domains for matrix and fracture in two dimensions, and linked them through a similar mass transfer term. Gao and Ghassemi 33 developed a fully three-dimensional coupled model of pore elasticity during hydraulic fracturing. The model takes into account both the induced stress caused by HF opening and the induced stress caused by fluid filtration. The model is used to analyze the stress difference between homogeneous reservoir and heterogeneous reservoir. Salimzadeh et al 34 did the similar work with Gao and Ghassemi. However, the effect of NFs on reservoir stress redistribution was not considered in their study.
Based on the above analysis, we propose a relatively simple and efficient hydraulic fracturing model for fractured-porous elastic reservoirs based on the FCM, cohesive unit method, fluid-solid coupling theory, and FVM. The coupling effect of fluid flow and geological stress under random distribution of NFs can be considered in the model. In this paper, the NF is equivalent to a physical model controlled by tension bond and shear bond. After considering the influence of NFs, the reservoir stress showed local stress mutation and strong heterogeneity. Finally, we analyze the sensitivity of stress inversion. The results are of great significance for understanding stress redistribution, fracture propagation, and coupling between HFs, matrix, and NFs after fracturing.

| NF continuum model
In naturally fractured reservoirs, such as shale and tight sandstone, the permeability of the reservoir matrix is very small and it is difficult to extract the fluid relying only on the matrix. 35,36 In engineering, multicluster fracturing in horizontal wells is generally adopted to activate reservoirs with NFs. At this time, NFs become the main channel of oil and gas production and dominate the flow mode of reservoir fluid. 37 In this model, the FCM proposed by McKenna and Reeves 38 is used to characterize the characteristics of NFs. The permeability of element grid is characterized by the NF inclination angle, approach angle, opening angle, and fracture spacing. The advantage of this method is that the NF grid can no longer be processed separately, and the results obtained under the condition that the element is small enough can also meet the precision requirements. [39][40][41] Based on Chen's research, the NF permeability is dispersed to the continuous unit. The permeability tensor of NFs defines permeability for each grid cell in the model domain. For one fracture set, the permeability tensor is shown as follows 42 : where k ij is the permeability tensor of grid center node and n 1, 2, 3 is the unit normal to the fracture plane in the x, y, and z direction, respectively. k nf is the permeability of NF。n 1 ,n 2 ,n 3 can be expressed as follows: where = 90 • −dip and = strike − 90 • . Through Equations (1) and (2), the permeability tensor can be used to represent the heterogeneity caused by NFs. This heterogeneity is determined by the opening, spacing, inclination, and approach angle of NFs. NF permeability is mainly determined by fracture opening 43 : where w NF is the reduction of fracture aperture (fracture closure). d is the fracture spacing. When there are multiple cracks in the same grid area, as shown in Figure 1A,B, the permeability tensor can be expressed as the sum of the permeability tensor of each NF, namely: where N is the number of fracture sets.
The k m xx , k m yy , and k m zz values are calculated for each grid block based on the parameters of fracture aperture, spacing, strike, and dip probability distributions that are defined for each fracture set n.

| Rock failure criterion
During the process of hydraulic fracturing, a large amount of fracturing fluid enters the reservoir.
These fracturing fluids are filtered into NFs and matrices, changing pore, and fluid pressure in the reservoir. Tensile or shear failure occurs when pore pressure or fluid pressure in NFs meets the set failure criteria, as shown in Figure 1C-F. This kind of damage will cause the change in NF opening and then greatly increase the permeability of NF, thus playing the role of reservoir reconstruction. According to Warpinski criterion, 44 shear failure occurs when the shear stress acting on the fracture plane is greater than the shear strength of NFs.
Here is as follows: where is the shear stress on the NF wall surface, MPa; 0 is the cohesion of NFs, MPa; is the friction coefficient of NF, dimensionless; basic is the basic friction angle, usually between 30° and 40°; n is the normal stress on the NF surface, MPa; p NF is the fluid pressure in the NF, MPa. When the matrix grid is small enough, the pore pressure of the grid can be accelerated to equal the fluid pressure in the NF. 39,45 As fracturing time increases, fluid pressure in the reservoir and in NFs continues to rise. When its value is greater than the normal stress of the NF surface, tensile failure occurs to the NF, that is as follows: (2) n 1 = cos 180 sin 180 n 2 = cos 180 cos 180 where K t is the tensile strength of natural cracks, MPa. Based on Warpinski criterion, the normal stress and shear stress on the NF plane under the influence of NF inclination angle and approach angle can be expressed as follows 46,47 : where f is the force exerted on the natural crack wall surface, MPa; n is unit normal vector of NF surface, dimensionless.

| The opening of a NF
The opening of NF in shear failure can be expressed as follows 48 : where K s is shear directional stiffness of NF, MPa/m; φ dil represents shear expansion angle of NF, °; Δ represents the effective shear stress of NF, MPa; the effective shear stress can be expressed as follows: where eff represents the effective stress, MPa.
When tensile failure occurs, the NF opening of tensile failure can be expressed as follows: where K n is normal stiffness of natural crack, MPa/m.
For the tension failure zone of NF, when shear stress exists in the NF surface, the zone is an overlapping zone of tension and shear failure in real time. According to linear elasticity theory, the total fracture opening includes the normal opening of NF opening, shear expansion opening, and original opening, which can be expressed as follows: For the NF zone with only shear slip, its opening is the shear expansion opening and the original opening, which can be expressed as follows: By substituting Equation (11) and Equation (12) into Equation (3), the permeability after activation of NFs can be obtained.

| HF propagation model
In our model, the matrix cell grid is an orthogonal grid, and the HF unit is embedded in the grid layer with a 0-thickness cohesive unit, as shown in Figure 2. Therefore, the crack extension trajectory in this paper is consistent with the set cohesion unit. Under the influence of heterogeneous distribution of NFs, the stress distribution of reservoir and HF surface will show strong heterogeneity. The calculation of HF parameters will be the result of coupling fluid pressure within the fracture, reservoir permeability, pore pressure, and normal stress on the fracture wall.

| Material balance equation
According to the principle of material balance, the material balance equation in a single HF is as follows 49 : where q is the flow rate within the HF, m 3 /s; h f represents HF height, m; w HF is HF opening, m; s represents the coordinate of fracture length direction, m; t is the time, s; q L is the filtration rate of fracturing fluid, m/s。 At the same time, according to the principle of material balance, the amount of fracturing fluid injected should be equal to the increment of fracture volume plus the loss of fracturing fluid, so the global material balance equation is as follows: The total injection amount of fracturing fluid should be equal to the sum of flows of each HF: where q T is the total pump injection flow, m 3 /s; q i represents the flow rate of HF i, m 3 /s; L HF,i is the length of HF i, m; m is the number of HF.

| Flow equation
Hydraulic fracturing measures are carried out with slippery water. The fluid has Newtonian fluid characteristics and is incompressible, so the pressure drop equation in HF can be expressed as follows: where p HF is the fluid pressure in the HF, MPa; μ is liquid viscosity, mPa·s.

| Fracture width equation
In conventional fracturing simulation, the fracture width in the above equation can be obtained by the following equation 34 : where ν is Poisson's ratio, dimensionless; z represents the direction coordinate of fracture height, m; c is the closure stress of the fracture wall surface, MPa; G is shear modulus, MPa.
However, when there are multiple HFs extending at the same time and the interference of NFs, the stress interference effect between each fracture should be considered. Therefore, in the model in this paper, when calculating the stress change in the elastic medium through the FVM, it is necessary to obtain the displacement of each HF wall surface (cohesive element surface) and finally obtain the dynamic fracture width through displacement. The slit width of this model can be expressed as follows 50,51 : where n represents the unit normal vector on the HF surface, dimensionless; u + , u − represents the displacement on the left and right sides of the HF, as shown in Figure 2. The HF parameters can be determined by coupling Equation (19) with the intrafracture flow Equation (17) and the material balance Equation (15). When the stress on the HF wall is strongly heterogeneous, the displacement of the left and right walls of the HF may be different, resulting in uneven distribution of fracture width.

| Fracture propagation criterion
According to the large number of rock mechanic experimental results, the crack initiation mode can be divided into the open-, slip-, and tearing-type fractures according to their causes. In the model, the cohesive element method is used to predetermine the propagation path of HF, which is perpendicular to the direction of minimum horizontal principal stress. Under the condition of single-point injection, the initiation of the HFs is considered as an open type. The stress intensity factor of the open fracture is expressed as follows: where K I represents the stress intensity factor of the open type at the crack end, MPa ⋅ √ m; L i is the length of the HF i , m; E (k) is the integral coefficient of the open-type fracture, dimensionless.
When the stress intensity factor satisfies the following conditions, the fracture can be extended.
where K Ic is the fracture toughness of the rock, MPa ⋅ √ m.

| Fluid filtration
In order to establish a complete fluid-solid coupled model, Darcy's law is used to establish a coupling fluid filtration model with permeability, fluid pressure, and pore pressure as the variable terms 52,53 : where S represents the exchange area for fluid filtration, m 2 ; f represents the permeability coefficient, 1/m; p p is the pore pressure, MPa.

| Coupling model of fluid and geomechanics
As shown in Figure 3, fluid filtration into the reservoir will lead to increased pore pressure in the reservoir matrix and fluid pressure in NFs, which will have two effects. On the one hand, with the increase in pore pressure, when the failure criteria (maximum tensile stress criterion and/or Coulomb criterion) are satisfied, the NF will be damaged, and the NF permeability will change. This will result in a significant increase in the overall permeability of the reservoir. On the other hand, nonuniform pore pressure increases, resulting in nonuniform change in induced stress. The influence of these two aspects will eventually affect the calculation of HF width and fracture expansion. In our model, the reservoir matrix has the characteristics of porous media. The activated NF will make the permeability of porous media show strong heterogeneity. In order to study the change of reservoir stress from a macroscopic perspective, it is assumed that it satisfies the coupling model of fluid and geological stress deduced by Biot. 54,55 The governing equation of this coupled system is based on mass conservation and linear-momentum balance. The mechanical deformation equation can be expressed as follows: where and b represent the total stress tensor and the single-phase fluid bulk density respectively, kg∕m 3 ; g is gravitational acceleration, m∕s 2 . Combined with Biot's theory, considering the influence of pore elasticity, the relationship between stress and strain and between strain and pore pressure can be expressed as follows: 56 where the subscript i represents initial state; C dr is the Rank 4 elastic tensor; b is the Biot coefficient; I is the Rank 2 identity tensor; v = tr( ) is the volumetric strain; and = 1 2 ∇u +∇ T u is the strain tensor. u is the displacement vector containing three components, m; V and q represent fluid-flow rate (m/s) and a source/sink term (q = q f ∕S), 1/s; M is the Biot modulus, MPa. The relationship between the Biot modulus and the Biot coefficient can be shown as follows 57,58 : where c f is fluid compressibility, MPa −1 ; m represents porosity, dimensionless; K v is the bulk modulus of solid grain, MPa; K dr represents the drained bulk modulus (MPa), which can be show as Similarly, Lame constants and can be shown as Substitute the volumetric mean total stress v , Equation (11) can be shown as where v = 1 3 tr represents the trace of the stress tensor. Based on Darcy's Law, the fluid flow velocity can be expressed as follows:

F I G U R E 3 Coupling diagram of fracture parameters and reservoir deformation
If ignore the gravitational term, substituting Equation (14) in Equation (12), the relationship between strain and pore pressure can be obtained: Equations (24) and (33) are called fixed-strain split. During the calculation, the iteration stops when the convergence criteria are reached for both equations. As demonstrated by Kim et al, 59 (29), (30), (33), and (24) into Equation (23), ignore the gravitational term, and the relationship between displacement and pore pressure can be obtained.

The Equation (34) can be expressed in a different form:
where n and n − 1 are the current timestep and the previous timestep.

| Workflow of calculation
The goal of this work is to build a numerical model of hydraulic fracturing in fractured reservoirs, in which NFs and the reservoir matrix are considered as porous media with pore elasticity through the FCM. Then, the effects of NFs on reservoir stress, reservoir displacement, and wall morphology of HFs were studied by coupling pore elasticity effect with fluid flow in HFs. In the calculation, the length of each extension of HF is determined artificially, so the time required is determined by the material balance equation. It is worth noting that in the fracture expansion module, we assume the initial fluid pressure and time within the HF, and then, through the nonlinear iteration of pressure and time, the fluid pressure, time, and fracture width meeting, the accuracy requirements can be obtained. During the iteration, the crack width is obtained by coupling fluid and geomechanics model, because we use the fully coupled nonlinear iterative method to calculate. Therefore, the fracture extension, fluid and geomechanics coupling, and NF failure should be iterated in the same time step until the overall error is satisfied. The calculation flow of the model is shown in Figure 4. In the figure, i represents the number of iterations. Solid and dotted lines represent the workflow path and parameter coupling path.

| VALIDATION
In this section, the fluid and geomechanics coupling numerical model is analyzed by using the presence volume method. Then, using PKN analytic solution, Abaqus software was compared with our crack propagation model.

| Numerical model
The numerical model is discretized in time and space by FVM. The discrete time is an implicit method with first-order precision, while the discrete space contains implicit and explicit methods, in which the majority is dependent on the Gaussian linearization method. 60 Tang et al 61 expressed the volume of each control unit into an integral form and gave the discrete details. Geomechanics Equation (22) can be rewritten as In Equation (37), the left-hand side is the implicit surface diffusion term. The first term on the right of the equality is explicit surface-diffusion term, and the second term is the explicit pressure-coupling term. The explicit constant term represents the initial state. Substituting the filter term into the source and sink term, and the fluid flow Equation (36) can be rewritten as In Equation (38), the first term to the right of the equal sign is the pressure term from a previous timestep and the second term is the explicit displacement coupling term. The last term represents source/sink term.

| HF propagation
Based on the numerical model of hydraulic fracturing of fractured reservoirs in this paper, we established a single HF expansion physical model. The effects of NFs on HF propagation and reservoir stress are considered in the model. Then, the results of fracture parameters calculated by PKN model and Abaqus finite element software were compared with those calculated by our model. The comparison results are shown in Figure 5, which shows the change in fracture width and HF length with time at the injection point. As can be seen from the figure, the variation trend of fracture length and width is basically consistent with that obtained by the analytic solution of PKN model and Abaqus software. The difference between them is caused by the PKN model and Abaqus software's difference in filtration treatment. In addition, the fracture width and length of our model fluctuated with time, which was caused by reservoir filtration and strong heterogeneity of stress distribution. The main parameters in the validation are shown in Table 1.

| Case study and model parameters
The study on the distribution of reservoir stress during hydraulic fracturing is of great significance for understanding the stress interference and the design of temporary plugging fracturing diversion. To study the effect of hydraulic fracturing on stress redistribution in fractured-porous media reservoirs, a multifracture extension model was established ( Figure 6). The length, width, and height of the physical model are 400 m, 60 m, and 40 m, respectively. In the model, the length of grid cells along the X-axis is 0.5 m, and the length of grid cells along the Y-axis is 1 m. The pore pressure at the boundary is constant, and the displacement perpendicular to the boundary is zero. NFs are randomly distributed in fractured reservoirs, and their positions are controlled by the approach angle and inclination angle. In this paper, the inclination angle is set as 90°, and the distribution of the approach angle is shown in Figure 8A. In addition, NFs are set as virtual fracture surfaces with tensile and shear resistance, and their destruction is determined by Warpinski criterion. The permeability of the reservoir matrix varies from 0.01 mD to 0.1 mD, as shown in Figure 7B. The main geological parameters and construction parameters in the simulation are shown in Table 2.

| Effects of NFs on stress redistribution
Based on the model established, the stress contour map after hydraulic fracturing in fractured and porous media reservoirs was obtained, as shown in Figure 8. The pore pressure distribution in Figure 8A is different from that of homogeneous reservoir after fracturing. In fractured and porous media reservoirs, the contour lines of pore pressure are highly heterogeneous and extremely uneven. At the same time, the isolines in Figure 8B,C,D,E also show strong heterogeneity. Stress in fractured-porous media reservoir shows local abrupt change behavior after fracturing. For better identification, Figure 8F is a local magnification diagram of the stress contour map along the x-axis ( Figure 8E). The stress distribution in the Figure 8F shows strong stress fluctuation, and the stress contour is staggered. The mutant contour area is distributed in blocks. In order to reflect the sudden change of stress more vividly, the horizontal principal stress along the x-axis and shear stress are represented as three-dimensional cloud maps, as shown in Figure 9A,B. For the conventional homogeneous model, the locations with stress mutations occur at the tip of the HF, while other regions generally do not. As can be seen from Figure 9, in our model, in addition to the phenomenon of stress mutation at the tip of HF, many stress mutations also occur in other local areas. The sudden stress in these local areas is shown as the sudden increase and the sudden decrease in the stress value ( Figure 9A,B). With the development of hydraulic fracturing, NFs are gradually destroyed under the action of pore pressure. This damage changes the permeability distribution of the porous media reservoir, resulting in a strong heterogeneity of fluid flow in the reservoir. It is this highly heterogeneous fluid flow that causes the reservoir to exhibit sudden stress changes. This sudden change in stress is a characteristic phenomenon of reservoir stress in fractured-porous media during hydraulic fracturing.

| Effect of initial stress difference on stress inversion
According to elastic mechanics, since there is no shear stress, the maximum horizontal principal stress under the initial condition is perpendicular to the horizontal wellbore; that is, it forms a 90° included angle with the x-axis, and the minimum horizontal principal stress forms a 0° included angle with the x-axis. The inversion of horizontal principal stress means that the angle between the maximum horizontal principal stress and the x-axis becomes approximately 0°, while the angle between the direction of the minimum horizontal principal stress and the x-axis becomes approximately 90°. The maximum horizontal principal stress ( max ), minimum horizontal principal stress ( min ), and included angle ( ) in the Figure 10 can be obtained from the following equation: The study of stress inversion is of great significance for understanding the mechanical mechanism in porous elastic media, predicting the trajectory of HFs, and designing repeated fracturing. As can be seen from Figure 10, the angle between the maximum horizontal principal stress in the green area and the x-axis approaches 0°, which means that the original minimum horizontal principal stress in this area has changed to the maximum horizontal principal stress due to hydraulic fracturing and stress interference between multiple fractures. At the same time, stress inversion did not occur in the area adjacent to the HF wall, but occurred in the area a certain distance away from the HF wall. This is quite different from the conventional crack propagation model. Compared with Figure 10A,B,C, with the increase in stress difference, the stress reversal region decreases sharply. Figure 11 shows the distribution curve of the angle ( ) along the x-axis. When the stress difference is 1 MPa, there are more points staying at 0°, which indicates that the stress reversal occurs in these areas. With the increase in stress difference, when the stress difference is 3 MPa, only three points have stress reversal behavior. However, when the stress difference is 4 MPa, there is no point at the position of 0°, which indicates that there is no stress inversion in the whole area. Although there is no stress inversion, the direction of stress will be offset due to shear stress. As the shear stress is small ( Figure 8C), the deviation angle is also small.

| Effect of cluster spacing on stress inversion
Cluster spacing has a great influence on the stress distribution during hydraulic fracturing. Reasonable control cluster spacing can be prepared for secondary fracturing design. Figure 12 shows the effect of cluster spacing on stress inversion. Within a certain range of cluster spacing, with the decrease in cluster spacing, the stress inversion region between HFs decreases gradually, while the external stress inversion region increases slowly. When the cluster spacing is 5 m, the stress inversion region between HFs completely disappears. Figure 13 also shows the distribution curve of the included angle along the X-axis. As can be seen from the figure, when the cluster spacing is 5 m, no point between HFs is at the position of 0°. When the cluster spacing is 7-11 m, the point of 0° begins to appear in the middle. Therefore, with the decrease in cluster spacing and the enhancement of stress interference, the stress inversion between HFs is weakened, while the stress inversion in the external region is enhanced.

| Effect of Biot coefficient on stress inversion
For reservoirs with different characteristics, their Biot coefficients (b) are different. In general, the Biot coefficient can be determined by Equation (27). The effect of Biot coefficient on stress is shown in Figure 14. With the increase in b, the horizontal principal stress along the X-axis increases gradually, and its stress variation range is within 2 MPa. In addition, with the increase in b, the horizontal principal stress along the Y-axis increases gradually in region A, while decreases gradually in region b. The variation range of stress in region A is within 4.5 Pa, while the variation range of stress in region B is within 1.5 mPa. Figure 15 shows the effects of F I G U R E 8 (A) represents pore pressure distribution; (B) represents the horizontal principal stress along the y-axis; (C) represents the shear stress induced by hydraulic fracturing; (D) represents the stress difference distribution after fracturing. (E) represents the contour map of horizontal principal stress along the x-axis; (F) represents an enlarged sketch of a local area different Biot coefficients on the stress inversion region. As can be seen, with the decrease in the effective coefficient, the area of the stress inversion area gradually decreases. When the Biot coefficient is 0.4, there is no stress inversion area in the whole reservoir. In general, Biot coefficient has a significant effect on stress reorientation.

| Effect of Poisson's ratio on stress inversion
Based on the established model, the stress distribution when Poisson's ratio is 0.1-0.3 is calculated under the same condition of other parameters. The result is shown in Figure 16. Figure 16A shows the relationship between the horizontal principal stress and Poisson's ratio along the direction of the horizontal wellbore. With the increase in Poisson's ratio, the stress interference is weakened, so the F I G U R E 1 1 The influence of stress difference on the angle between the horizontal maximum principal stress and the x-axis F I G U R E 1 2 The angle distribution between the maximum principal stress and the x-axis after fracturing at different cluster spacing F I G U R E 1 3 The influence of cluster spacing on the angle between horizontal maximum principal stress and x-axis stress gradually decreases. The relation curve between horizontal principal stress along the y direction and Poisson's ratio is shown in Figure 16B. The influence of Poisson's ratio on stress shows two trends. As shown in Figure 16C, with the increase in Poisson's ratio, the stress in this region gradually decreases. In the other region, with the increase in Poisson's ratio, the stress gradually increases, as shown in Figure 16D. Overall, Poisson's ratio has little influence on the stress. When Poisson's ratio varies between 0.1 and 0.3, the corresponding stress varies within 1 MPa. Therefore, when Poisson's ratio changes between 0.1 and 0.3, the area of stress inversion changes very little, as shown in Figure 17.

| CONCLUSIONS
In this paper, the extension behavior of three clusters of HFs in fractured and porous media reservoirs is simulated based on cohesive elements, FCM and FVM. In our model, the reservoir matrix is considered as a porous medium, and NFs are virtually present, with the ability to resist tension and shear. Activation of NFs was determined by Warpinski criterion. Then, the HF extension, fluid flow, and stress calculation are obtained by fully coupled nonlinear iteration. This paper focuses on the study of stress effects in naturally fracturedporous media reservoirs and the effects of multiple factors on stress inversion. The following conclusions are obtained. F I G U R E 1 4 The effect of Biot coefficient on horizontal principal stress along the x (a) and y (b) axes F I G U R E 1 5 The angle distribution between the maximum principal stress and the x-axis after fracturing under different Biot coefficients 1. The NF will make the reservoir stress show strong heterogeneity after fracturing, and the prominent feature is the regional stress mutation. The sudden change in stress shows the sudden increase or decrease in stress. 2. With the increase in the initial horizontal stress difference, the stress inversion becomes more difficult. When the stress difference exceeds 3 MPa, the stress inversion area is almost 0. The stress interference increases with the decrease in cluster spacing. The stress inversion region between HFs decreases, while the external stress inversion region increases. 3. With the decrease in Biot coefficient, the area of stress inversion region will decrease obviously; Poisson's ratio has little influence on stress inversion. When Poisson's ratio changes between 0.1 and 0.3, the change range of stress is within 1 MPa, and the stress inversion area is almost unchanged.
The results of our study are of great significance for understanding the stress distribution after hydraulic fracturing in fractured and porous media reservoirs, analyzing the fracture extension trajectory and the design of secondary fracturing.