Effect of Natural Fractures on Stress Evolution of Unconventional Reservoirs Using a Finite Element Method and a Fracture Continuum Method

Refracturing, temporary plugging, and infilling well design play an important role in the development of reservoirs. The prediction of stress distribution can provide the basic guiding theory for the design and implementation of these techniques. In this paper, a fully-coupled three-dimensional production model based on the finite element method (FEM) and fracture continuum method (FCM) for naturally fractured reservoirs is presented to study the effects of fluid consumption on the reservoir stress. Furthermore, the effects of natural fractures on the stress re-distribution and stress re-orientation are also studied. The model also considers the influence of natural fractures on the permeability, and the effect of the effective stress on natural fracture openings, pore-elastic deformation, and fluid consumption. An analytical solution model and Eclipse were used for the comparison, which verifies the accuracy of the model results. Based on two cases of one cluster of fractures and three clusters of non-planar fractures, the research results revealed that natural fractures have a significant influence on the surrounding drainage, stress distribution, and stress re-orientation during the development. Under the influence of natural fractures, the production of the fluid along the direction of natural fractures is significantly easier, and it is highly probably that the insufficient consumption area is perpendicular to the direction of natural fractures. Compared with the conventional model, the stress distribution in the proposed model is deflected to a certain extent under the flow mode dominated by natural fractures, which is significantly prominent in the non-planar fracture model. Due to the effect of natural cracks, the absolute values of the stress, displacement, and stress difference in this model are relatively larger than those in the conventional model. Moreover, the re-orientation angles of the maximum principal stress are significantly different. After considering the natural cracks, there was an increase in the change in re-orientation and the re-orientation range. The research findings reported in this paper can be used to predict the initiation, extension, and steering process of temporary plugging fracturing fractures and refracturing fractures in fractured reservoirs.


Introduction
In oil and gas exploitation, the pore pressure in reservoirs gradually decreases in accordance with the development. At this point, the pore pressure gradient due to the heterogeneity of the reservoir pore pressure distribution induces a change in the initial reservoir stress distribution [1]. In conventional reservoir research, the reservoir is generally considered as a porous medium elastomer. When evaluating the coupling law between fluid flow and geological stress in the porous elastic model, the fluid flow is generally considered to occur in three main directions. However, the fluid flow is due to the pressure gradient. For reservoirs such as shale, wherein natural fractures develop and dominate the reservoir fluid flow, the fluid flow is mainly along the direction of natural fractures. In this flow mode, the conventional permeability representation cannot accurately describe the flow of fluid. Therefore, the relationship between the natural fracture permeability and permeability tensor [2], the relationship between the permeability and fracture width [3], and the relationship between the effective stress and natural fracture width were established [4]. The permeability tensor can then be used to describe the flow of the fluid in other directions, and the variation in the width of the natural fracture can be used to describe the permeability evolution during production. Hence, a three-dimensional fluid-solid coupling mathematical model in porous media was developed as an equivalent simulator for the analysis and evaluation of the stress evolution and stress steering in fractured reservoirs during depletion.
Several studies have been conducted on the fluid-solid coupling in porous media; however, the effect of natural fractures on the reservoir stress was considered in a few of these studies. The study on the coupling between fluid and geomechanics in porous media can be divided into two parts. The first aspect is the coupling of fluid and geomechanics during injection [5][6][7][8][9][10]. The other aspect is the coupling of the fluid and geological stress in the depletion [11][12][13]. The change in the reservoir stress due to the fluid loss during hydraulic fracturing is the result of the coupling between the flow and geomechanics, which is a process that occurs within a short-time period. Ghassemi et al. established a three-dimensional fluidsolid coupling model for hydraulic fracturing and then solved the model using the finite element and displacement discontinuity methods. The distribution of the pore pressure and stress in the vicinity of the hydrofractures was obtained, and the failure of the rocks in the vicinity of the hydrofractures was analyzed using the Mohr-Coulomb failure criterion [14]. Salimzadeh et al. proposed a fully coupled threedimensional (3D) finite element model for the hydraulic fracturing of permeable rocks and then evaluated the applicable range of classical analytic solutions in the limit case using this model [15]. Gao and Ghassemi established a pore elastic model in the process of full 3D heterogeneous hydraulic fracturing using the coupling method of fluid flow and geomechanics and then solved it using the finite element method. The induced stress distribution and steering around hydraulic fractures in heterogeneous reservoirs were confirmed as functions of time and space [16].
In addition to changes in the reservoir stress during fluid injection, changes were observed when the reservoir was depleted. Oil and gas production is a significantly long-term reservoir consumption process [17][18][19]. Gupta et al. evaluated the influence of production consumption and completion migration on the reservoir stress, and the findings of the study indicated that the stress is more easily diverted under lowstress difference conditions. Moreover, the fracture characteristics, pressure exhaustion, and interwell interference have a significant influence on the infill well design and refracturing design in the latter stages [20]. Using a state-of-the-art coupled symmetric Galerkin boundary-element method (SGBEM) and the finite element method (FEM) model, Roussel and Sharma simulated the influence of prior production from the offset wells on the propagation direction of fractures initiated in an infill horizontal well. A stress reorientation zone was observed near the offset well after several years of production, which exhibited a significant influence on the fracture extension in the infill well [11]. Safari and Ghassemi proposed a numerical evaluation model for the changes in stress due to the depletion of the infill wells and then analyzed the bending of hydraulic fractures in the disturbed stress field and the connectivity of the hydraulic fractures between wells [21]. The failure effect and critical parameters that determine the fracture bending were systematically evaluated [21,22]. The results revealed that the depletion perturbed the stress tensor of the stratum around the fractured horizontal well. The perturbation stress field is related to stress/formation anisotropy, fluid fluidity, pore pressure, bottom hole working pressure, and Biot's constant. Sangnimnuan et al. evaluated the influence of the hydrofracture morphology on the stress evolution in the production process using the embedded discrete fracture method (EDFM) and the finite volume method (FVM) [1]. The study revealed that the fracture morphology has a significant influence on the stress distribution, surrounding drainage, and stress reorientation. Moreover, with a decrease in the stress difference, the steering is easier to achieve. The propagation path and radius of curvature of the hydrofractures can be predicted using the fracture propagation model combined with the varying stress field. In this paper, reservoir simulation was integrated with rock geomechanics to predict the well poststimulation productivities [23]. Several numerical and field studies have been presented to study the impact of the reservoir depletion on the pore pressure and the stress variation in the infill well zones [24]. Rezaei et al. presented refracturing analysis using a 2D poroelastic plain strain displacement discontinuity model built on the work of Zhang, Ghassemi and Zhang, and Chun [25][26][27][28]. Kumar and Ghassemi presented a geomechanical perspective to better understand the problem of "frac-hits" in horizontal well refracturing and to design solutions for it using geomechanics analysis and modeling [13]. The numerical analysis was based on a fully coupled 3D model "GeoFrac3D" with the capabilities to simulate multistage fracturing of multiple horizontal wells.
However, for reservoirs with natural fractures, the evolution of the reservoir-induced stress is influenced by the pressure consumption, hydrofracture morphology, and natural fractures, among others. Moreover, the influence of natural fractures is significant. Most researchers know that natural fractures have a significant impact on production during oil and gas production. However, when natural fractures and fluid-solid coupling are considered at the same time, the treatment of natural fractures becomes a difficult task. In summary, no research was found on the effect of the natural fracture heterogeneity on the induced stress in reservoir production. This paper summarizes and presents the results of the analysis conducted on the effect of natural fractures on the stress and stress diversion under the condition of plane fractures and nonplane fractures. This study was a combination and expansion of the fluid flow/geomechanics coupling theory and natural fracture heterogeneity in the pore elastic model. Compared with previous studies, our model is closer to the actual geological situation of the reservoir, so the results are more reliable and accurate. The research findings reported in this paper can be used to predict the initiation, extension, and steering process of temporary plugging fracturing fractures and refracturing fractures in fractured reservoirs.

Mathematic Model
2.1. Natural Fracture Continuum Method. In naturally fractured reservoirs such as shale and tight sandstone, the 2 Geofluids permeability of the reservoir matrix is significantly small and it is difficult to extract the fluid based only on the matrix [23,29]. In engineering, multicluster fracturing is generally employed in horizontal wells to activate reservoirs with natural fractures. The natural fractures then become the main channel of the oil and gas production and then dominate the flow mode of the reservoir fluid. Given that the permeability of natural fractures is significantly greater than that of the matrix, the contribution of the matrix permeability to fluid flow is neglected in this model and only the influence of natural fractures is considered [29]. The distribution angle of natural fractures developed in shale reservoirs is typically consistent, i.e., the inclination angle and approximation angle are generally distributed within a certain range of the main direction [30]. Therefore, it was assumed that the approximation angle and inclination angle of allnatural fractures in the model are consistent. In this model, the fracture continuum method proposed by McKenna and Reeves is used to describe the characteristics of natural fractures. The permeability of the element grid is characterized by the natural fracture inclination angle, approximation angle, opening angle, and fracture spacing [31]. The advantage of this method is that the natural fracture grid cannot be processed separately, and the results were obtained under the condition that the element is sufficiently small to meet the precision requirements [3,32]. Based on previous research [3], the natural fracture permeability was dispersed to the continuous unit. The permeability tensor of natural fractures defines the permeability for each grid cell in the model domain. For one fracture set, the permeability tensor can be expressed as follows: where k is the permeability tensor; n 1 , n 2 , and n 3 are the units normal to the fracture plane in the x, y, and z directions, respectively; and k nf is the permeability of the natural fracture. Moreover, n 1 , n 2 , and n 3 can be expressed as follows: sin ξ π 180 , n 2 = cos ς π 180 cos ξ π 180 , where ς = 90°− dip and ξ = strike − 90°. Using Equations (1) and (2), the permeability tensor can be used to represent the heterogeneity caused by natural fractures. This heterogeneity is determined by the opening, spacing, inclination, and approach angle of natural fractures. In the model, the fluidsolid coupling of the depletion is considered and the stress distribution of the reservoir changes with respect to time, thus leading to changes in the natural fracture permeability. The natural fracture permeability is mainly determined by the fracture opening as follows [3]: where δ is the reduction of the fracture aperture (fracture closure) and d is the fracture spacing.
To accurately describe the natural fracture permeability as a function of time, a coupling model of the natural fracture opening and effective stress was introduced. Bandis et al. evaluated the crack opening and effective normal stress and then proposed a semi-logarithmic relationship [33]: where σ n is the normal stress on the crack surface and ω and χ represent the constants in the semilog fracture closure/stress relationship.
The effective normal stress that acts on the crack surface can be expressed as follows: where θ is the angle between the normal to the fracture and the σ H direction and σ H ′ and σ h ′ , respectively, represent the maximum effective principal stress and the minimum effective principal stress, MPa. By combining Equations (4) and (5), the following was obtained [4]: where n is the effective stress ratio, which can be expressed as

Principle of Effective Stress.
Rock is a porous medium composed of a solid skeleton and pores. In addition, there are fluids (oil, gas, and water) stored in the pores of the skeleton. The fluid in the rock can bear or transfer pressure, which is defined as the pore pressure. Moreover, the stress transmitted through the contact surface between the rock particles is the effective stress. The Biot effective stress can be expressed as follows [34,35]: where σ ij ′ is the effective stress, σ ij is the total stress, δ ij is the Kronecker symbol, and α represents the Biot constant, which can be defined as follows [36]: where K v and K s represent the volume compression modulus of the rock and the compression modulus of the solid particles, respectively. Moreover, α is typically close to 1; thus, α = 1 was used.

Stress Balance Equation.
The theory of saturatedunsaturated seepage evaluates the movement of two or more fluids through pores. This study addressed the flow of a single-phase fluid in a porous medium. Based on the principle of virtual work, it can be known that at a certain moment, the virtual work of the rock mass of a unit is equal to the virtual work generated by the force that acts on the unit (physical force and surface force); thus, where F is the surface force, f is the body force, and δε and δu are the virtual strain and virtual displacement, respectively. The constitutive relationship expressed in an incremental form is where D ep is the elastoplastic matrix and dε l is the particle compression due to the pore fluid pressure, which only deforms in the positive direction without shear direction. The specific expression is as follows: where m = ½1, 1, 1, 0, 0, 0 T . Formula (7) can be expressed as follows: For α = 1, the combination of Formula (12) and the sum Formula (9) can be expressed as follows: Given that the seepage continuity equation contains a time-related term, to couple the stress with the seepage, the time derivative of the virtual work Equation (13) is required. The specific expression is as follows: where g is the acceleration vector of gravity, c f is the compressibility of the fluid, and ϕ is the matrix porosity.
In the process of oil and gas development, hydraulic fracturing is generally used to generate hydrofractures to increase the productive drainage area. As shown in the Figure 1, a hydrofracture was located at the y-z section. Hence, the continuity equation can be expressed as follows: where q Γ represents the interchange term between the fracture and matrix.

Model Calculation and Validation
3.1. Finite Element Dispersion. The shape function can be defined as follows: where B and N u represent the shape function of the displacement, N p is the shape function of pore pressure, u is the nodal displacement, and p w is the pore pressure of the element node.

Geofluids
By substituting Equation (17) into Equation (15), the solid phase finite element formula can be obtained after simplification, as follows: where In the analysis of the seepage field, there are two types of boundary conditions. The first is the flow boundary condition, and the second is the pore pressure boundary condition. The flow boundary conditions can be expressed as follows: where n c is the normal unit of the flow boundary. This boundary condition can be used as the outer boundary condition. When there are production fractures in the reservoir after fracturing, the flow boundary condition can also be considered as the internal boundary condition or fracture boundary condition. The pore pressure boundary conditions can be expressed as follows: where p fb is the pore pressure value at the known boundary. The Galerkin method can be used, which can be expressed as follows: where a and b are arbitrary functions, A is the governing equation, and B is the continuity equation through the boundary.
By substituting Equation (15) as A, Equation (20) as B, and the formal function expression (17) into Equation (22) and by setting a = −b, after the simplification, the following can be obtained: where

Geofluids
By combining Equation (18) with Equation (23), the geomechanic-seepage coupling equation can be expressed as follows:  Figure 1) and the same model was established by Eclipse software for a simulation comparison. The permeability, porosity, and reservoir thickness of the model were 5 md, 0.1, and 50 m, respectively. The boundary of the model was closed; thus, the reserves in the model were constant. The other main parameters were the same as those shown in Table 1. In the proposed model, the flow model was established without considering natural fractures. Fluid-solid coupling is not considered in Eclipse software. The daily output and cumulative oil production after one year of production were obtained by adopting the method of constant pressure production, as shown in Figure 2. The results of the model were approximately consistent with those of Eclipse. Due to the different factors considered by the models, the change process was slightly different. Therefore, the model in this paper demonstrated a high accuracy.

Case Studies
4.1. Basic Parameters. As discussed in this section, a threedimensional reservoir physical model was established to evaluate the influence of anisotropy due to natural fractures on the reservoir stress during the oil and gas development. To fully demonstrate the significant influence of natural fractures on the generation of induced stress, a model was established that does not consider the influence of k xz , k yz , and k xy . However, the influence of the permeability in the main direction was considered for comparison. As shown in Figure 3, the length, width, and height of the model were 400 m, 400 m, and 100 m, respectively, and the hydraulic fracture was located at the cross section of the y-z-axis. The domain was discretized to 50 four-node cells in the x-and y-directions and 20 four-node cells in the z-direction. For the treatment of the reservoir boundary, the nodes on the left and right surfaces of the model were free to move along the x-axis and not the y-axis. All the nodes at the top and bottom outer surfaces were free to move along the y-axis and not the x-axis. Only the vertical displacement of all the nodes on the outer surface of the top and bottom could be achieved. Moreover, the six boundaries were nonflow boundaries, i.e., the pore pressure in the reservoir was continuously reduced in the process of production. In all the calculations below, the downhole flow velocity was calculated using the Peaceman equation [37]. The other parameters in the calculation are shown in Table 1.

Geofluids
comparison. Compared with the boundary conditions in the model, the calculated parameters were consistent with the proposed model. During the production, the pore pressure decreases continuously due to the continuous fluid extraction and the pressure is the lowest near the fracture. Figure 5 compares the distribution of the reservoir pore pressure after 5 years of production with and without considering the influence of natural fractures. As can be seen from the figure, due to the influence of natural fractures, there was a change in the fluid flow pattern. Compared with the conventional model, the pore pressure dissipation area was greater and it was easier to reduce the pore pressure along the direction of natural fractures for the formation of a clockwise deflecting pressure drop funnel. Therefore, under the natural fracture distribution condition, the fluid in the upper right and lower left corners could be easily extracted, whereas the fluid in the upper left and lower right corners could not be easily extracted. This can serve as an essential guide to the design of infill wells.
According to the principle of effective stress, the reservoir stress changes under the influence of the constant dissipation of the pore pressure. For the conventional model, the distribution of the stress σ xx , σ yy , and σ zz in the reservoir after 5 years of production is shown in Figures 6(a), 6(b), and 6(c). Due to fluid consumption, the reservoir stress distribution exhibited a high complexity. Moreover, there was a decrease in σ xx near the crack as a result of the pressure dissipation, whereas σ xx increased on the left and right parts of the domain to support the pressure depletion in the x-direction. In addition, there was a decrease in σ yy near the crack as a result of the pressure dissipation, whereas σ yy increased at the top and bottom parts of the domain to support the pressure depletion in the y-direction. Moreover, σ zz mainly exhibited a decrease in stress in the central region of the domain. Considering the influence of natural fractures on the stress evolution, the results are shown in Figures 7(a), 7(b), and 7(c). Compared with the conventional model, the stress distribution of the proposed model is similar to that of the conventional model. However, considering the influence of tensors, the pressure dissipation of the proposed model was more rapid and the changes in stress were greater. The most significant difference is that the distributions of σ xx , σ yy , and σ zz exhibited a degree of deflection. The stress distribution deflection was found to be significantly different from the stress distribution near the fractures in the conventional model, and this difference has a significant influence on the fracture initiation direction and fracture extension of refracturing. As can be seen in Figures 8(a) and 8(b), the distribution of the shear stress on the x-y section exhibited a clockwise deflection when compared with the conventional   The main difference was that there was a more significant change in the absolute value of the displacement after the influence of natural fractures was considered. Figures 11(a), 11(b), and 11(c) present the σ xx , σ yy , and σ zx stress distributions of the proposed model and conventional model after 5 years of production on the z-x section, respectively. As can be seen from the figure, the stress distribution of the conventional model exhibited a symmetrical distribution. Compared with the conventional model, the overall trend of the stress distribution in the proposed model was in good agreement. However, the stress distribution was symmetrically distributed under the influence of natural fractures, and there was a relative increase in the variation range of the stress value.
Significant attention has been directed toward the reservoir stress evolution and surrounding drainage in the process of fluid recovery, given their significance for the latter refrac-turing and completion of infill wells. The change in the direction of the horizontal principal stress direction has a significant influence on the design of refracturing. Due to the heterogeneity of the permeability and natural fractures, the heterogeneity consumption of the pore pressure may result in stress reorientation and changes in the direction of the horizontal maximum principal stress. As can be seen from Figures 12(a) and 12(b), a region with a stress difference less than 0 was observed near the fracture after one year of production. In addition, the stress in this region was reversed. It can also be found that the time when the stress difference changes the most is not the 20 years with the longest production life but the 10th year. This is because with an increase in the production life, most of the fluid in the easily flowing region was produced, whereas the fluid in the previously relatively noneasily flowing region was gradually produced. Therefore, there was a decrease in the difference between the pore pressures in the regions, which led to a decrease in the stress differences and an increase in the influence range. Figures 12(c) and 12(d) present the shear stress profiles on Paths A and B. Compared with the conventional model, the shear stress absolute value of the proposed model is influenced more significantly by natural cracks.
Refracturing is one of the most effective methods for latter reservoir development. The distribution of the reservoir stress after the production and the direction of the horizontal 8 Geofluids maximum principal stress are two critical factors that determine the fracture pressure, fracture initiation direction, and fracture extension trajectory of refracturing. The stress evolution in the production process was previously analyzed. For the solution of the change in reorientation of the maximum horizontal principal stress (σ H max ) on the x-y section, the action of the vertical principal stress (σ zz ) was neglected; thus, the original direction of the maximum horizontal principal stress is shown in Figure 13(a). After 5 and 10 years of production, there was a significant change in the direction of σ H max . In the conventional model, the change area of the principal stress direction was in the shape of a vertical ellipse, whereas the change area of the σ H max direction in the proposed model was in the shape of an inclined ellipse. The degree of inclination is related to the location of natural fractures. The area with the largest reorientation change of σ H max was near the hydrofracture. Compared with the conventional model, the change in reorientation of σ H max in the proposed model was influenced more significantly by natural fractures (θ 1 ′ > θ 1 , θ 2 ′ > θ 2 ). In the conventional x-y section (d) Figure 9: (a) and (b) present u x and u y on the x-y section without considering natural fractures, respectively, and (c) and (d) present u x and u y on the x-y section considering natural fractures, respectively. 9 Geofluids model, the change in reorientation of σ H max was 90° (  Figure 14(a)); however, the change in reorientation of σ H max in the proposed model was 110°under the influence of natural cracks (Figure 14(b)). Similar to the change in the stress difference, the change in the reorientation was the largest at the 10th year, not the 20th year. Moreover, the range of the reorientation areas is proportional to production time.

Nonplanar Multifracture Study.
In the actual segmented multicluster fracturing process, there is stress interference between the fractures when multiple fractures are simultaneously extended, which leads to the case wherein the middle fractures are restrained and the fractures on both sides are extended in a nonplanar manner. Hence, the displacement discontinuity method was used to simulate the extended trajectories of three clusters of hydrofractures with spacings of  10 Geofluids 50 m. The three hydraulic fractures are considered as production fractures. It was assumed that all three clusters of fractures pass through the reservoir in the z-direction. The other calculation parameters are the same as those in Table 1. Under this condition, the reservoir stress distribution after 1, 5, and 10 years of production was simulated. Figure 15 presents the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production using the conventional model. The pore pressure distribution in the domain was axially symmetric. However, under the influence of natural fractures, there was an initial decrease in the pore pressure along the strike direction of the natural  Due to the influence of the nonplanar fracture geometry, the stress distribution after production exhibited a high complexity. Without considering the effect of natural fractures, the stress distribution of the conventional model was axially symmetric. Compared with a single fracture, there was an increase in the rates of the pore pressure dissipation and stress decrease and there was an increase in the size of the influence range, as shown in Figures 17(a), 17(b), and 17(c). The stress (σ xx ) on the left and right of the domain and the stress (σ yy ) on the top and bottom of the domain increased to support the pressure depletion of the middle area. When the influence of natural fractures was considered, the stress distribution in the proposed model exhibited a significant degree of deflection, as shown in Figures 18(a), 18(b), and 18(c). The deflection of the stress distribution confirms that natural fractures have a significant influence on the reservoir fluid flow and production.
Compared with a single crack, the distribution of the principal stress was more complex. The change in reorientation of σ H max was θ 1 ′ > θ 1 , θ 2 ′ > θ 2 , as shown in Figure 19. Due to the influence of the nonplanar crack geometry, the rotation angle of the maximum principal stress direction in the conventional model reached 135°. In the proposed model, the rotation angle of the principal stress direction near the crack was concentrated at 140°. The reorientation region affected by the natural fractures was larger than those of the conventional models. Therefore, under the action of the nonplanar fracture geometry and natural fractures, the evolution of the reservoir stress tends to be complex in the production process. From a comparison of Figures 20(a) and Figure 18(a), the profile of the change in reorientation on Path B was found to be significantly different due to the influence of the hydraulic fracture geometry. During the first year of production, the central region did not turn; however, it turned in accordance with an increase in the production time. From a comparison of Figures 20 and Figure 18, there was a similarity in that the greatest change in the reorientation of σ H max occurred in the 10th year and the change in reorientation was correlated with the reservoir depletion time.

Conclusion
A fully coupled 3D pore elastic model based on the FEM was developed to analyze the stress evolution and stress reorientation in heterogeneous porous media in the reservoir depletion process. In the model, the continuous natural fracture model is used to characterize the effect of natural fractures on the reservoir fluid flow. The relationship between the natural fracture width and effective stress was used to establish a fluid cluster single fracture model and a three-cluster nonplane fracture model for an example analysis. There was a good agreement between the analytical solutions and x-y section y (c) Figure 15: (a), (b), and (c) present the pore pressure distribution after 1, 5, and 10 years of nonplanar fracture production without considering the influence of natural fractures. 12 Geofluids    14 Geofluids numerical results. By comparing the conventional models, it was confirmed that natural cracks have a significant influence on the stress evolution, surrounding drainage, and change in reorientation in the depletion process. Under the influence of natural fractures, the fluids in reservoirs along the direction of the natural fractures are easier to recover, and areas that are perpendicular to the natural fracture surface are prone to underconsumption. In addition, under the flow mode dominated by natural fractures, the stress distribution in the proposed model is deflected to a certain extent, especially in the nonplanar fracture model. Due to the action of natural cracks, the absolute values of the stress, displacement, and stress difference in the proposed model are relatively greater and the change in the reorientation of the maximum principal stress steering is significantly different. After considering the natural cracks, the steering angle and steering range increased. Moreover, this directly results in significant differences in the design of the late-stage fracturing temporary plugging steering and refracturing technology for conventional reservoirs. The simulation results revealed that it is very important to consider the effects of natural fractures on the stress redistribution, surrounding drainage, and reorientation in production. The findings of this study can serve as a theoretical guide for the design of refracturing, temporary plugging steering, and infilling wells and help predict the direction of new fractures and the trajectory of fracture extensions.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.