Numerical Analysis of the Stress Shadow Effects in Multistage Hydrofracturing Considering Natural Fracture and Leak-Off Effect

: As a critical technological approach, multistage fracturing is frequently used to boost gas recovery in compact hydrocarbon reservoirs. Determining an ideal cluster distance that effectively integrates pre-existing natural fractures in the deposit creates a fracture network conducive to gas movement. Fracturing fluid leak-off also impacts water resources. In our study, we use a versatile finite element–discrete element method that improves the auto-refinement of the grid and the detection of multiple fracture movements to model staged fracturing in naturally fractured reservoirs. This computational model illustrates the interaction between hydraulic fractures and pre-existing fractures and employs the nonlinear Carter leak-off criterion to portray fluid leakage and the impacts of hydromechanical coupling during multistage fracturing. Numerical results show that sequential fracturing exhibits the maximum length in unfractured and naturally fractured models, and the leak-off volume of parallel fracturing is the smallest. Our study proposes an innovative technique for identifying and optimizing the spacing of fracturing clusters in unconventional reservoirs.


Introduction
The critical energy scenario in China necessitates the urgent augmentation of unconventional hydrocarbon reserves, including gas hydrates [1].In response to the rational and efficient extraction of tight reservoirs, the advent of horizontal well fracturing technology was facilitated by hydraulic fracturing [2].In 2010, M. J. Mayerhofer proposed the Stimulated Reservoir Volume (SRV) concept while studying the microseismic technology and fracture changes in Barnett shale [3].As the injection of fracturing fluid continues, the net pressure of the primary fracture continually increases.When this pressure surpasses the difference in the principal horizontal stress in the in situ stress and the tensile strength of the rock, secondary fractures will form [3,4].During the segmented fracturing process of horizontal wells, the already-formed hydraulic fractures introduce a stress shadow effect, also known as induced stress.This stress shadow effect increases the subsequent perforation fracturing pressure, altering the original geo-stress state and the fracture expansion direction.The fracture formation process exhibits the "3S" effect of stress accumulation, stress transfer, and stress shadow [5,6].In the 1980s, Warpinski and Teufel began studying the interaction between fractures, providing an analytical solution for the stress field of a hydraulic fracture with constant height and finite length [7].In 2004, Fisher confirmed this phenomenon through microseismic mine research and explained the stress shadow effect using microseismic detection methods.He pointed out that when hydraulic fractures interact with natural fractures, various forms such as crossing, termination, and offset occur, with the impact range of the stress shadow effect being 1.5 times the height of the fracture [8][9][10][11][12].Numerous studies have been conducted on the topic commonly referred to as the stress shadow area.A literature review was also undertaken on this subject, detailing the stress disturbances generated in the zone between fractures, confirmed by applying microseismic technology [12][13][14].Hence, the stress shadow phenomenon is the principal element of forming an extensive, intricate fracture network, thereby increasing complexity within the shale reservoir [15,16].Roussel and Sharma pointed out that a preceding hydraulic fracture would alter the stress field distribution of subsequent fractures within a specific fracture spacing, causing the fracture to turn and become tangent to the primary fracture [14,15,17,18].Negal simulated the stress shadow effect of multiple fractures under varying reservoir conditions, spacing, and in situ stress ratios using continuous and discrete elements [19].Therefore, deciphering the intrinsic mechanism of the stress shadow effect and utilizing it to enhance the fractures' complexity and connectivity are paramount.
The primary objective of volumetric hydraulic fracturing in large-scale shale reservoirs is to facilitate the communication between artificial and natural fractures, thereby establishing a three-dimensional fracture network.This has been confirmed by means of micro-seismic measurements and surface tilt-meter measurements [7,[9][10][11]20].Numerous studies have scrutinized the interplay between hydraulic fractures and pre-existing fractures.Blanton conducted a series of experiments exploring the impact of the intersection angle between hydraulic and natural fractures and the differential stress on the state of fracture crossing [21].Gu and others proposed a simple standard for the intersection extension of hydraulic fractures and frictional natural fractures based on linear elastic fracture mechanics theory, verified through triaxial hydraulic fracturing experiments [22,23].Olson and Bahorich replaced rocks with gypsum blocks for hydraulic fracturing experiments, studying the crossing process of hydraulic fractures and bonded natural fractures [24,25].Should there exist an abundance of natural fractures, they would be prioritized as the principal determinant in altering the trajectory of the fractures in contrast to other influences.
Various models and simulators that consider the coupling of multiple physical fields have been suggested to emulate the process of hydraulic fracturing [26].Although some new methods have been developed, such as the phase field method [27] and the peridynamic method [28], methods including the finite difference (FD) method [29], finite element (FE) method [30], the extended finite element method (XFEM) [31,32], the displacement discontinuity method (DDM) [33], and the discrete element method (DEM) [34], are still standard and effective methods for hydraulic fracturing simulation.Fu [35] suggested a linked hydro-mechanical model utilizing the finite element and finite volume method to address the geomechanics and hydrodynamics issue.Riahi [36] proposed a model using DEM to anticipate the progression of intricate fractures on diverse natural fracture patterns.Coupled thermal-hydro-mechanical (THM) FEM techniques were implemented to simulate the flow of fluid and heat reaction in fracture problems [37,38].Li [39] constructed a thoroughly integrated three-dimensional thermo-hydro-mechanical (THM) hydraulic fracturing (HF) model to optimize multi-well completion design, utilizing a hybridized finite volume/finite element approach, contemplating both HF propagation and natural reactivation.Owen [40] presented a unified finite element-discrete element (FE-DE) methodology to address large-scale pragmatic issues, encompassing fracture solids and particulate media.This synergized finite/discrete element approach offers a more accurate resolution in the context of hydraulic fracturing challenges.Traditional methods, such as the continuum FEM models, are dependent only on the combination of damage and the tensile model to capture porous flow; the integration of the proppant transport model presents challenges when executing a pure continuum FEM [41].The discrete element approach often assumes that the edge of an element pre-determines the direction of fracture propagation, thereby precluding the formation of genuinely complex fractures that are frequently observed in highly heterogeneous reservoirs [42,43].In addressing the difficulties associated with the traditional FE and DE methods, this study employs an adaptive FE-DE method that dynamically updates geometries and meshes to enhance fracture propagation modeling, providing more accurate stress analyses and reliable fracture paths [44].
The conventional metric utilized for assessing the efficacy of multistage fracturing is typically the volume of the hydraulic fracture [39].Indeed, the most credible gauge to appraise the fracturing outcome is the volume of gas generated post-hydro-fracturing, flow-back, and a series of procedures.However, few techniques can accurately portray the dynamics inherent in the exhaustive procedure of multistage hydro-fracturing and gas generation.There remain substantial challenges in evaluating the fracturing impact based on gas production or recovery rate.In this investigation, we employ the adaptive FE-DE technique incorporated in the ELFEN TGR software, utilizing local re-meshing technology to enhance precision.To juxtapose the impact of the stress shadow effect with the natural fracture, we simulated two scenarios employing two fracturing sequences, namely, sequential and simultaneous fracturing, adjusting the perforation cluster space from 12.5 m to 100 m.

Adaptive FE-DE Modeling for Multistage Fracturing
We built an integrated model capable of dealing with hydraulic fracturing problems in tight reservoirs, considering pre-existing fractures, hydro-mechanical coupling, and leak-off effects [44].The embraced numerical method is a hybrid FEM/DEM coupled hydro-mechanical approach for simulating hydraulic fracture.The governing equations of the methodology, geo-mechanical equations, proppant transport, fracture fluid flow, and the effect of leak-off are briefly introduced as follows.

Numerical Methods for TH Coupling
An ELFEN TGR analysis is split into five critical phases according to the hydraulic fracture field practice [45]: 1.
Initiation of model effective stress, pore-forces, and fracturing fluid pressures; 2.
Flow-back and cleanup operations; 5.
Gas production.
Considering the interplay between fracturing fluid and reservoir gas, the simulation requires a sophisticated multi-phase flow process to mimic the complete physical procedure.The simulation accomplishes the integration of solid stress, fluid, and fracture fluid flow.

Basic Equations
Assuming that the acceleration of fluid relative to solid and the convective term can be neglected, the momentum balance equation for the solid field is as follows [46]: The equation for the liquid permeation of porous flow in the rock structure, integrating mass conservation with Darcy's principle, can be articulated through Equation (2) [47].Equation (3) controls the fracture fluid flow within the fracture network [46].
Assuming that we follow the parallel plate theory, k f r can be calculated through the following formula [48]: Water 2024, 16, 1308 4 of 20 The relationship of storage coefficient S f r , fracture stiffness K f r n , and the bulk modulus of the fracturing fluid K f r f [46] can be written as: The equation for mass conveyance can be expressed as: For the evaluation of gas extraction and manufacturing, the gas permeation formula and the gas grid formula, which amalgamate mass preservation with Darcy's principle, are provided as follows [47]:

Model of Leak-Off and Proppant Movement
In our model, the matrix of tight reservoir rock was assumed to be isotropic.Figure 1 illustrates a nearly brittle substance's joint uniaxial tension stress-strain behavior.The stress-strain curve gradient H and Young's modulus E during the damage phase are defined as follows; more information can be found in reference [49].
Assuming that we follow the parallel plate theory,   can be calculated through the following formula [48]: The relationship of storage coefficient   , fracture stiffness    , and the bulk modulus of the fracturing fluid    [46] can be written as: The equation for mass conveyance can be expressed as: For the evaluation of gas extraction and manufacturing, the gas permeation formula and the gas grid formula, which amalgamate mass preservation with Darcy's principle, are provided as follows [ In our model, the matrix of tight reservoir rock was assumed to be isotropic.Figure 1 illustrates a nearly brittle substance's joint uniaxial tension stress-strain behavior.The stress-strain curve gradient H and Young's modulus E during the damage phase are defined as follows; more information can be found in reference [49].
The liquid cannot be compacted or condensed.

2.
The flow is locally analogous to the movement between two smooth, parallel surfaces.

3.
The flow is streamlined, characterized by a low Reynolds number.
Water 2024, 16, 1308 5 of 20 The fluid behavior can be written using the non-Newtonian Power Law model: During the hydraulic fracturing process, the injected fluid not only serves for fracturing purposes but also permeates into the minute pore throat radii of the shale grains.Incorporating the effect of leak-off constitutes another step towards enhancing the reliability of the numerical method.Over time, the exposed fracture surface exhibits this behavior through a one-dimensional Carter leak-off model [50].The relationship of initial volume loss V sp per unit and spurt time t sp is: where the leak-off coefficient, C, is provided by: Additional leak-off parameters can be discovered in the cited materials [46].Proppant is another important reason for the stress shadow phenomenon.The 1D proppant equation is [51]: The concluding component of this proppant transportation formula pertains solely to the fluid, as there is no proppant seepage into the matrix.The flux for the fluid phase within the fracture could be substituted by a defined mixed flux as follows; the other related methods and algorithms can be found in reference [44].

Coupling Strategy
The method employs a staged coupling scheme, explicitly [52] solving the mechanical governing equation and implicitly [53][54][55][56] solving the fracture fluid flow equation, enhancing both simulation efficiency and stability.The fluid pressures are not equal in the structure and network fields, (p e n c1 ̸ = p i n c1 .The implicit analysis of the network field progresses sequentially in time from coupling times t c1 and t c2 , yielding fluid pressures p i n c1 and p i n c2 .Figure 2 illustrates the communication between the implicit and explicit solvers, which occurs across successive coupling times t c1 and t c2 .

Adaptive Re-Meshing and Coarsening Strategy
As a fracture tip advances, automatic re-meshing sustains a refined mesh around the tip.After the fracture tip has moved on, the necessity for a refined mesh diminishes, allowing for a coarsening process to simplify the grid along the formed break segments.This practice reduces the total number of elements within the domain, thereby augmenting efficiency [45].The parameters consist of two components: fracture insertion options and coarsening data.The fracture insertion options dictate the re-meshing at the advancing fracture tip, encompassing mesh density, density factor, and bubble size.On the other hand, the coarsening data control the coarsening of the mesh trailing the advancing fracture tip, comprising coarsening frequency, density, the threshold factor of coarsening density,

Adaptive Re-Meshing and Coarsening Strategy
As a fracture tip advances, automatic re-meshing sustains a refined mesh around the tip.After the fracture tip has moved on, the necessity for a refined mesh diminishes, allowing for a coarsening process to simplify the grid along the formed break segments.This practice reduces the total number of elements within the domain, thereby augmenting efficiency [45].The parameters consist of two components: fracture insertion options and coarsening data.The fracture insertion options dictate the re-meshing at the advancing fracture tip, encompassing mesh density, density factor, and bubble size.On the other hand, the coarsening data control the coarsening of the mesh trailing the advancing fracture tip, comprising coarsening frequency, density, the threshold factor of coarsening density, non-coarsening zone size, and maximum coarsening zone size.The definition of each parameter can be found in reference [57].

Numerical Models and Procedure
As demonstrated in Figure 3a, a 2D model of multistage fracturing was utilized, extending across a 1000 m length and 600 m height and encompassing five clusters.The initial perforation fracture lengths are approximately 2 m near the horizontal well.The inter-cluster spacing 's' was systematically varied from 12.5 m to 100 m to appraise the stress interference derived from neighboring perforation clusters.The target reservoir lies at a perpendicular depth of 3000 m, situated within an environment featuring normal crustal stress and a moderate temperature of 45 °C.The horizontal in situ stress  ℎ in the x direction and vertical in situ stress   in the y direction are 40 MPa and 44 MPa, respectively.To authenticate the influence of natural fractures on the propagation of fractures within the reservoir, we incorporated two types of natural fractures near the perforation clusters in the model.This was done to simulate the stress variations typically seen among pre-existing fractures.Such a simulation aims to understand better the behavior of fracture propagation under realistic conditions, where natural fractures within the reservoir can significantly impact the direction and extent of new fractures created during the fracturing process.As evidenced by Figure 3b, the realm of the pre-existing natural fracture is centrally situated within the reservoir, measuring 500 m in length and 200 m in height.The size and location of this pre-existing natural fracture zone are strategically chosen to align with the refined mesh constraint, which is meant to enhance the precision and efficiency of the computation.Doing so ensures that the computational resources are focused

Numerical Models and Procedure
As demonstrated in Figure 3a, a 2D model of multistage fracturing was utilized, extending across a 1000 m length and 600 m height and encompassing five clusters.The initial perforation fracture lengths are approximately 2 m near the horizontal well.The inter-cluster spacing 's' was systematically varied from 12.5 m to 100 m to appraise the stress interference derived from neighboring perforation clusters.The target reservoir lies at a perpendicular depth of 3000 m, situated within an environment featuring normal crustal stress and a moderate temperature of 45 • C. The horizontal in situ stress S h in the x direction and vertical in situ stress S v in the y direction are 40 MPa and 44 MPa, respectively.To authenticate the influence of natural fractures on the propagation of fractures within the reservoir, we incorporated two types of natural fractures near the perforation clusters in the model.This was done to simulate the stress variations typically seen among pre-existing fractures.Such a simulation aims to understand better the behavior of fracture propagation under realistic conditions, where natural fractures within the reservoir can significantly impact the direction and extent of new fractures created during the fracturing process.As evidenced by Figure 3b, the realm of the pre-existing natural fracture is centrally situated within the reservoir, measuring 500 m in length and 200 m in height.The size and location of this pre-existing natural fracture zone are strategically chosen to align with the refined mesh constraint, which is meant to enhance the precision and efficiency of the computation.Doing so ensures that the computational resources are focused on the most critical areas of the model, providing a higher resolution of results where it matters most.Maintaining consistency in the simulations, the dimensions of the refined mesh area within the reservoir remained constant across models for various scenarios.This allows for a meaningful comparison of results across different scenarios, as the same level of precision is applied in each case.Adaptive refinement at the fracture tips in the original finite element meshes is a strategic approach to improve the accuracy and reliability of the resulting fracture network.This adaptive refinement focuses computational resources on the areas of most significant importance, which are the fracture tips in this case.mesh area within the reservoir remained constant across models for various scenarios.This allows for a meaningful comparison of results across different scenarios, as the same level of precision is applied in each case.Adaptive refinement at the fracture tips in the original finite element meshes is a strategic approach to improve the accuracy and reliability of the resulting fracture network.This adaptive refinement focuses computational resources on the areas of most significant importance, which are the fracture tips in this case.

Un-fractured model
Naturally fractured model Regarding the stress shadow effect, fracture sequence is another significant triggering element.In our model, we hypothesized two types of fracturing sequences: sequential and simultaneous, as illustrated in Figure 4.For the sequential sequence, fracture propagation is stage-by-stage, where fracturing advances from the cluster on the right (the first) to the one on the left (the fifth).Conversely, the simultaneous sequence sees all clusters commencing fracturing concurrently.Here, the volume of the injected fluid equals that of sequential fracturing.Still, the volume for each cluster is allocated by initial conditions and stress interaction within the fracture zone-the ultimate fracture network results from stress re-orientation and fluid volume rivalry between the perforation clusters.To negate the influence of fluid competition, we introduced a parallel fracturing case, depicted in Figure 4c, as a comparator within the simultaneous sequence.In this sequence, the total volume of injected fluid remains constant.However, each cluster receives an equal fluid volume allocation.Consequently, fractures at different clusters propagate independently, and propagation deflection is induced solely by the stress shadow effect.By incorporating this scenario into simultaneous fracturing, we can delve deeper into the internal mechanisms of the stress-shadow effect.Regarding the stress shadow effect, fracture sequence is another significant triggering element.In our model, we hypothesized two types of fracturing sequences: sequential and simultaneous, as illustrated in Figure 4.For the sequential sequence, fracture propagation is stage-by-stage, where fracturing advances from the cluster on the right (the first) to the one on the left (the fifth).Conversely, the simultaneous sequence sees all clusters commencing fracturing concurrently.Here, the volume of the injected fluid equals that of sequential fracturing.Still, the volume for each cluster is allocated by initial conditions and stress interaction within the fracture zone-the ultimate fracture network results from stress re-orientation and fluid volume rivalry between the perforation clusters.To negate the influence of fluid competition, we introduced a parallel fracturing case, depicted in Figure 4c, as a comparator within the simultaneous sequence.In this sequence, the total volume of injected fluid remains constant.However, each cluster receives an equal fluid volume allocation.Consequently, fractures at different clusters propagate independently, and propagation deflection is induced solely by the stress shadow effect.By incorporating this scenario into simultaneous fracturing, we can delve deeper into the internal mechanisms of the stress-shadow effect.

Basic Parameters
We simulated five fracture stages for each sequential fracturing case and changed the space between the perforation cluster from 12.5 m to 100 m with a two-fold increase.The duration of the fracturing, flow-back, and gas production stages for the sequential and simultaneous sequences are listed in Table 1 and Table 2. Typical formations and fracture

Basic Parameters
We simulated five fracture stages for each sequential fracturing case and changed the space between the perforation cluster from 12.5 m to 100 m with a two-fold increase.The duration of the fracturing, flow-back, and gas production stages for the sequential and simultaneous sequences are listed in Tables 1 and 2. Typical formations and fracture properties were taken from the literature and are shown in Table 3. Table 4 lists the parameters for the local re-meshing and coarsening strategy.Every instance was performed utilizing ELFEN TGR on a desktop PC with an Intel Xeon Bronze 3.40 GHz CPU.Several numerical results could be calculated using the method and model described in the previous section, in which the morphology of cracks and first principal stress can reflect the inner mechanism during the fracturing process.When the stress shadow effect takes effect, the shear stress among the fractures will change significantly so that the orientation of the fracture propagating would deflect.As a result, we decided to select these typical results to analyze the influence of multi-factors (cluster-spacing, fracturing sequence, and natural fractures) on the stress shadow effect in multistage hydraulic fracturing.As mentioned before, we set the space to 100 m to avoid the stress shadow effect caused by the previous fractured stages; in this simulation, the space between the adjoined perforation clusters ranges from 12.5 m to 100 m; the aim was to observe the variety in stress field across the fracture area when we change the spacing.

Effect of Fracture Sequence and Cluster Spacing
As depicted in Figure 5, the induced fracture network displays varying inter-fracture spacings in sequential and simultaneous fracturing.Furthermore, the image also delineates the distribution of the first principal stress following the completion of the fracturing operation in the initially unfractured model.In the context of the sequential fracturing scheme, the fracture propagation sequence initiates from the first-stage cluster.It concludes at the fifth, with each stage being subject to a slick water injection volume of 150 m 3 .The cumulative injection volume across all five clusters amounts to 750 m 3 .Figure 5A(a) illustrates that the initial three perforation clusters exhibit predominantly vertical fractures.An initial deflection is observed at the upper segment of the fracture from the fourth cluster.The fifth fracture converges with its precedent due to the influence of stress shadowing, thereby forming an elongated, unified fracture.
When the spacing is increased to 25 m, the initial fracture follows the direction of the maximum stress, and the fracture's morphology closely resembles a vertical straight line.As the fracturing operation unfolds, a deflection from the initial fracture is observable in the second and subsequent fractures, with the deflection angle incrementally increasing as fracturing progresses.During the fourth and fifth stages, the fractures evolve into a unique, single-wing morphology, propagating laterally away from the horizontal well.The deflection angle culminates at its maximum value at the culmination of the final fracturing stage.
Figure 5A(c-e) illustrates the conditions across spacing intervals of 50 m, 75 m, and 100 m.Within these intervals, the deflection angle demonstrates a diminution trend, and the fracture propagation orientation increasingly aligns with the vertical direction.Notably, Water 2024, 16, 1308 10 of 20 fractures initially favor propagation along the trajectory of principal stress, as depicted vertically in the figures.Nevertheless, once the pre-existing fractures are imbued with proppant, the influence of stress shadowing on the directional orientation of subsequentstage fracturing is comparatively marginal, particularly in juxtaposition with simultaneous fracturing.This suggests nuanced interplays of stress fields and fracture propagation directions during staged fracturing operations.
as fracturing progresses.During the fourth and fifth stages, the fractures evolve into a unique, single-wing morphology, propagating laterally away from the horizontal well.The deflection angle culminates at its maximum value at the culmination of the final fracturing stage.
Figure 5A (c-e) illustrates the conditions across spacing intervals of 50 m, 75 m, and 100 m.Within these intervals, the deflection angle demonstrates a diminution trend, and the fracture propagation orientation increasingly aligns with the vertical direction.Notably, fractures initially favor propagation along the trajectory of principal stress, as depicted vertically in the figures.Nevertheless, once the pre-existing fractures are imbued with proppant, the influence of stress shadowing on the directional orientation of subsequent-stage fracturing is comparatively marginal, particularly in juxtaposition with simultaneous fracturing.This suggests nuanced interplays of stress fields and fracture propagation directions during staged fracturing operations.Under simultaneous fracturing, an identical total fluid volume of 750 m 3 , as applied in sequential fracturing, is infused into the horizontal well.This results in concurrent propagation across all fractures, thereby inducing a significant stress shadow effect.Figure 5B provides a comprehensive depiction of the primary principal stress cloud and fractures under the conditions of simultaneous fracturing.This comparison underlines the substantial impact of fracturing strategy on stress distribution and fracture propagation patterns.When the spacing is 25 m, the fracture originating from the final cluster extends well into the far-field region, surpassing the lengths of the other fractures.It is observed that fractures at the boundary areas (first and fifth) propagate significantly further than those located centrally (second to fourth), leading to a dramatic transformation of the stress field surrounding the central zone.The stress shadow effect in this central area significantly inhibits fracture propagation, to the extent that further growth of the aperture is precluded.In contrast, this effect facilitates the propagation of fractures in the boundary areas, enabling them to extend into the far-field region.The stress shadow effect intensifies with decreasing fracture spacing, indicating its critical role in fracture propagation dynamics.
Figure 6 offers an insightful representation of the critical outcomes relating to shear stress during fracture propagation and deflection.As depicted in Figure 6a, fractures begin propagation with the shear stress almost equally distributed on either side of the fractures.After the injection of the fracturing fluid, the stress field initiates a redistribution process, Water 2024, 16, 1308 11 of 20 causing the deflection of the fractures, notably in the first and fifth stages, as shown in Figure 6b.As the fracturing process continues into the subsequent stages represented in Figure 6c-f, the propagation of fractures in the fifth stage is facilitated by the stress shadow effect and the redistribution of slick water volume, making it more straightforward for the fractures to expand in the disturbed area as compared to the unbroken matrix.This is attributed to the stress shadow effect, which constrains the growth of the central fractures and encourages the propagation of fractures distanced from the wellbore.Consequently, under the conditions of a 25 m spacing, the fracture originating from the fifth stage can extend to the furthest extent within the reservoir.
stress shadow effect and the redistribution of slick water volume, making it more straightforward for the fractures to expand in the disturbed area as compared to the unbroken matrix.This is attributed to the stress shadow effect, which constrains the growth of the central fractures and encourages the propagation of fractures distanced from the wellbore.Consequently, under the conditions of a 25 m spacing, the fracture originating from the fifth stage can extend to the furthest extent within the reservoir.
The stress shadow effect between neighboring clusters diminishes upon increasing the fracture spacing.This weakening effect enhances the propagation of central fractures while curtailing the length of edge fractures, stabilizing the total fracture length.Figure 5B (c) illustrates that, at a 50 m spacing, the reduced stress shadow effect triggers the propagation of central fractures, with their lengths incrementally increasing as the spacing widens.Remarkably, once the fracture spacing extends to 100 m, the length of the fracture generated at the third stage surpasses even that of the fifth stage at the edge.This results in the most extended fracture length when compared to other spacings.The stress shadow effect between neighboring clusters diminishes upon increasing the fracture spacing.This weakening effect enhances the propagation of central fractures while curtailing the length of edge fractures, stabilizing the total fracture length.Figure 5B(c) illustrates that, at a 50 m spacing, the reduced stress shadow effect triggers the propagation of central fractures, with their lengths incrementally increasing as the spacing widens.Remarkably, once the fracture spacing extends to 100 m, the length of the fracture generated at the third stage surpasses even that of the fifth stage at the edge.This results in the most extended fracture length when compared to other spacings.
In parallel fracturing, each cluster is allocated an identical fluid volume of 150 m 3 , as in the case of sequential fracturing, ensuring the uniform propagation of fractures throughout the fracturing process.In tandem, the stress shadow effect intensifies, markedly modifying the direction of fracture propagation.Owing to the independent fluid supply designated for each cluster, as depicted in Figure 7a, the stress interaction between adjacent perforation clusters gives rise to a patchwork-like morphology, contrasting with that observed in Figure 5B(a) at the equivalent spacing.This phenomenon is attributed to the stress superposition effects occurring around the fracture tips, which not only halt the propagation of cracks at one end but also significantly deflect the propagation direction.
Once the spacing extends to 50 m, as shown in Figure 7b, the deflection angle of the central fracture approximates the vertical direction.In contrast, the fractures in the first and fifth stages deviate from the vertical direction, angling towards the reservoir due to the absence of the stress shadow effect at both ends.This trend continues in the subsequent intervals (Figure 7c,d), with the deflection angle decreasing in direct proportion to the increase in spacing.
Once the spacing extends to 50 m, as shown in Figure 7b, the deflection angle of the central fracture approximates the vertical direction.In contrast, the fractures in the first and fifth stages deviate from the vertical direction, angling towards the reservoir due to the absence of the stress shadow effect at both ends.This trend continues in the subsequent intervals (Figure 7c,d), with the deflection angle decreasing in direct proportion to the increase in spacing.

Effect of Pre-Existing Natural Fractures on Fracture Propagation
The existence of pre-existing natural fractures plays a pivotal role in the development of complex fracture networks within naturally fractured reservoirs.Table 5 presents the input parameters for the natural fractures as featured in the Discrete Fracture Network (DFN) model.Except for degree, parameters such as spacing, length, and persistence remain identical for both sets.
Sequential fracturing reveals significant differences in morphologies and stress fields compared to the un-fractured model.Figure 8A (a) showcases the final fracture network along with the first principal stress cloud maps at a spacing of 12.5 m.Notably, fractures initially propagate in the direction of maximum in situ stress until they intersect with natural fractures, a stark difference from what is observed in Figure 5 (a).

Effect of Pre-Existing Natural Fractures on Fracture Propagation
The existence of pre-existing natural fractures plays a pivotal role in the development of complex fracture networks within naturally fractured reservoirs.Table 5 presents the input parameters for the natural fractures as featured in the Discrete Fracture Network (DFN) model.Except for degree, parameters such as spacing, length, and persistence remain identical for both sets.Sequential fracturing reveals significant differences in morphologies and stress fields compared to the un-fractured model.Figure 8A(a) showcases the final fracture network along with the first principal stress cloud maps at a spacing of 12.5 m.Notably, fractures initially propagate in the direction of maximum in situ stress until they intersect with natural fractures, a stark difference from what is observed in Figure 5(a).
A noteworthy observation is that the convergence of fractures from adjacent clusters consolidates into a single extension, as demonstrated in the central position of the figure.The fracturing fluid migrates from the newly formed cracks to the previous stage, thereby reactivating the fracture in the upper stage and enabling its continued propagation.As the spacing between clusters increases and the distribution of pre-existing natural fractures is factored in, fractures are more likely to consolidate within smaller spacings than larger ones.
When changing the spacing from 12.5 m to 100 m within the naturally fractured model, it becomes clear that the stress shadow effect in fracture propagation diminishes compared to the un-fractured model.Concurrently, pre-existing natural fractures emerge as the predominant factor.As the spacing narrows, fractures from adjacent clusters are more likely to intersect with other fractures by propagating along the natural fractures.However, as the spacing widens, fractures tend to propagate vertically to reactivate distant natural fractures after negotiating the natural fractures.
Persistence (m) 15 15 When changing the spacing from 12.5 m to 100 m within the naturally fractured model, it becomes clear that the stress shadow effect in fracture propagation diminishes compared to the un-fractured model.Concurrently, pre-existing natural fractures emerge as the predominant factor.As the spacing narrows, fractures from adjacent clusters are more likely to intersect with other fractures by propagating along the natural fractures.However, as the spacing widens, fractures tend to propagate vertically to reactivate distant natural fractures after negotiating the natural fractures.Under simultaneous fracturing in the naturally fractured model, a pronounced stress shadow effect amongst the fractures hinders extensive propagation into far-field areas, confining the fractures primarily to areas proximal to the horizontal well.As depicted in Under simultaneous fracturing in the naturally fractured model, a pronounced stress shadow effect amongst the fractures hinders extensive propagation into far-field areas, confining the fractures primarily to areas proximal to the horizontal well.As depicted in Figure 8B, the path of fracture fluid propagation is initially dictated by the initial perforations, owing to the distribution of natural fractures, before it transitions into the adjacent natural fractures.The length of the fractures is predominantly determined by the distance between the perforation cluster and the natural fractures' location.Due to the proximity of the perforations, fractures intersect with those adjacent to them.However, as the spacing widens, this intersection phenomenon weakens but is not accompanied by long-distance propagation.The newly formed fractures are distributed on both sides of the wellbore.This suggests a complex relationship between perforation spacing, natural fracture distribution, and fracture propagation, which is crucial for effective fracturing strategies.

Discussion
The total fracture lengths under different fracturing sequences and perforation spacings are depicted in Figure 9a.Sequential fracturing exhibits the maximum total length across various perforation spacings among the three fracturing sequences.This can be attributed to the comparatively weaker stress shadow effect during each stage of fracture expansion in sequential fracturing.In subsequent stages, the stress shadow effect only alters the direction of expansion, with minimal impact on the overall length.For simultaneous fracturing, the total length of hydraulic fractures is the shortest and remains relatively consistent across different perforation spacings.This is due to the intense stress shadow effect between fractures during simultaneous fracturing, which restricts the overall expansion of fractures, resulting in a total length of approximately half that of sequential and parallel fracturing.The fracturing effect of parallel fracturing lies between simultaneous and sequential fracturing.Owing to the individual supply of fracturing fluid at each stage, the fractures expand fully in this stage; the stress shadow effect generated by adjacent fractures alters the expansion direction.The average fracture length per stage in sequential fracturing is 10.8% longer than in parallel fracturing.The fracture volume and leak-off volume of three fracture sequences is shown in Figure 10; we observe a consistent increase in leak-off volume as the process progresses.This reveals that the high leak-off effect significantly influences the dynamics of fractures during hydraulic fracturing.In sequential fracturing, upon completion of fracture extension in the current stage, the in situ stress in the formation exerts considerable compressive force on the formed fracture surfaces, leading to fracture closure and a reduction in fracture volume.The decline phase of the fracture volume occurs between each stage of sequential fracturing.Meanwhile, although minor fluctuations in fracture volume are observed in simultaneous and parallel fracturing, the fracture volume ultimately reaches its maximum value at the end of the fracturing process.Owing to the high leak-off effect, the volume of fluid loss constitutes a substantial proportion of the total fracturing fluid volume across different fracturing sequences, accounting for approximately 40% to 80% of the total fracturing.
Due to the intense stress shadow effect, simultaneous fracturing inhibits rapid growth in fracture length.Under the same leak-off parameters, more fracturing fluid is retained within the fracture.Thus, even though partial closure occurs under in situ stress, the overall fracture volume in simultaneous fracturing still surpasses that in sequential fracturing.On the other hand, parallel fracturing, which involves individual fluid supply at each perforation location, results in the most enormous fracture volume among the three fracturing sequences; consequently, it has the smallest fluid loss volume.
The leak-off volume generally decreases with increasing perforation spacing, as shown in Figure 11.Parallel fracturing exhibits the smallest fluid leak-off, amounting to just 43% of the total fracturing fluid volume when the spacing is 50 m.Compared to simultaneous fracturing and sequential fracturing, the former records the most significant fluid leak-off volume, peaking at 79% of the total fracturing fluid volume.The fluid leakoff volume of sequential fracturing lies between simultaneous and parallel fracturing.However, as the perforation spacing increases, the total fluid leak-off volume shows a declining trend.This is attributable to the diminishing stress shadow effect, allowing more fracturing fluid to be retained within the hydraulic fractures.The total fracture lengths in the naturally fractured model under sequential and simultaneous fracturing are illustrated in Figure 9b.It can be observed that when the perforation spacing is less than 25 m, the total length of hydraulic fractures in the naturally fractured model is not significantly different between sequential and simultaneous fracturing, amounting to 127.1 m and 86.91 m, respectively.When the spacing is small, the fracture networks formed by sequential and simultaneous fracturing interconnect after initial expansion to create a common network, resulting in negligible differences between the two.However, when the perforation spacing increases to over 50 m, sequential fracturing, compared to simultaneous fracturing, can connect to natural fractures to the greatest extent during each expansion stage, thereby forming a more complex fracture network.The fracture lengths are 3.1 times, 3.3 times, and 3.8 times that of simultaneous fracturing, respectively.This observation reflects the importance of perforation spacing and fracturing methods in exploiting the presence of natural fractures to maximize hydraulic fracture network complexity and length.Similar results can be derived when comparing the numerical model results with laboratory experiments [44,58,59].
The fracture volume and leak-off volume of three fracture sequences is shown in Figure 10; we observe a consistent increase in leak-off volume as the process progresses.This reveals that the high leak-off effect significantly influences the dynamics of fractures during hydraulic fracturing.In sequential fracturing, upon completion of fracture extension in the current stage, the in situ stress in the formation exerts considerable compressive force on the formed fracture surfaces, leading to fracture closure and a reduction in fracture volume.The decline phase of the fracture volume occurs between each stage of sequential fracturing.Meanwhile, although minor fluctuations in fracture volume are observed in simultaneous and parallel fracturing, the fracture volume ultimately reaches its maximum value at the end of the fracturing process.Owing to the high leak-off effect, the volume of fluid loss constitutes a substantial proportion of the total fracturing fluid volume across different fracturing sequences, accounting for approximately 40% to 80% of the total fracturing.Due to the intense stress shadow effect, simultaneous fracturing inhibits rapid growth in fracture length.Under the same leak-off parameters, more fracturing fluid is retained within the fracture.Thus, even though partial closure occurs under in situ stress, the overall fracture volume in simultaneous fracturing still surpasses that in sequential fracturing.On the other hand, parallel fracturing, which involves individual fluid supply at each perforation location, results in the most enormous fracture volume among the three fracturing sequences; consequently, it has the smallest fluid loss volume.
The leak-off volume generally decreases with increasing perforation spacing, as shown in Figure 11.Parallel fracturing exhibits the smallest fluid leak-off, amounting to just 43% of the total fracturing fluid volume when the spacing is 50 m.Compared to simultaneous fracturing and sequential fracturing, the former records the most significant fluid leak-off volume, peaking at 79% of the total fracturing fluid volume.The fluid leak-off volume of sequential fracturing lies between simultaneous and parallel fracturing.However, as the perforation spacing increases, the total fluid leak-off volume shows a declining trend.This is attributable to the diminishing stress shadow effect, allowing more fracturing fluid to be retained within the hydraulic fractures.

Conclusions
This study used a recently developed adaptive FE-DE model considering HM coupling and natural fractures [39] to analyze the influence of multi-factors (cluster-spacing, fracturing sequence, and natural fractures) on the stress shadow effect in multistage hydraulic fracturing.Our analysis encompassed an evaluation of the efficacy of various fracturing approaches: sequential, simultaneous, and parallel fracturing.This was done considering the influence of differential distances between the perforation cluster and the natural fractures.The conclusive total fracture length and morphology were taken as numerical outcomes for this assessment.Based on the numerical results presented, the suggestions and concluding observations are as follows: (1) Within the range of spacings evaluated in our model, the stress shadow effect modulates the direction of fracture propagation to some extent without significantly altering the total length in sequential fracturing scenarios.The impact of the stress shadow effect on fracturing is more subdued in this context when compared with other scenarios scrutinized in this study.(2) The prominence of the stress shadow effect increases in simultaneous fracturing scenarios when implemented with closer spacing, precipitating hydraulic interconnections between neighboring perforation clusters.This facilitates the extension of edge fractures into farther areas and amplifies the leak-off effect.On the contrary, larger spacing encourages the development of central fractures without contributing to an increase in total length.In simultaneous settings, parallel fracturing results in a more extended fracture length.(3) The distance between adjacent clusters primarily dictates the stress perturbation effect surrounding the natural fractures (NFs).In scenarios with closer spacing, adjacent fractures converge and propagate in a specific direction.In contrast, propagation is primarily along the maximum in situ stress direction in more extensive spacing settings.(4) The role of the leak-off effect in the hydraulic fracturing process is critical.It impacts the fracture length and volume and significantly influences the stress shadow effect among various fractures.The leak-off volume generally decreases with increasing perforation spacing attributable to the diminishing stress shadow effect among different fracturing sequences.
The findings from this research augment our understanding of the fundamental mechanisms underlying the stress shadow effect in multistage hydraulic fracturing.This process necessitates intricate computational coupling with various factors.Nevertheless, subsequent challenges persist in optimizing cluster spacing and well spacing in multiwell designs.These challenges need to be addressed from the microseismicity and gas production standpoint.We aim to explore these aspects in our future research endeavors.

21 Figure 2 .
Figure 2. Communication of network and structure fields during successive coupling time.

Figure 2 .
Figure 2. Communication of network and structure fields during successive coupling time.

Water 2024 ,
16, x FOR PEER REVIEW 8 of 21

Figure 8 .
Figure 8. Dynamic propagation of cracks and first principal stress (Pa) of sequential and simultaneous fracturing in the naturally fractured model for different spacings: (a) 12.5 m, (b) 25 m, (c) 50 m, (d) 75 m, and (e) 100 m.

Figure 8 .
Figure 8. Dynamic propagation of cracks and first principal stress (Pa) of sequential and simultaneous fracturing in the naturally fractured model for different spacings: (a) 12.5 m, (b) 25 m, (c) 50 m, (d) 75 m, and (e) 100 m.

Figure 9 .
Figure 9. Fracture length of multistage fracturing in the un-fractured and the naturally fractured model at different spacings: (a) un-fractured model and (b) naturally fractured model.

Figure 9 .
Figure 9. Fracture length of multistage fracturing in the un-fractured and the naturally fractured model at different spacings: (a) un-fractured model and (b) naturally fractured model.

Figure 10 .
Figure 10.Evolution of total fracture volume and leak-off volume of different fracturing sequences.

Figure 10 .
Figure 10.Evolution of total fracture volume and leak-off volume of different fracturing sequences.

Figure 10 .
Figure 10.Evolution of total fracture volume and leak-off volume of different fracturing sequences.

Figure 11 .
Figure 11.Evolution of total leak-off volume of different fracturing sequences.
Dynamic viscosity coefficient of the pore fluid µ n Dynamic viscosity coefficient of the fracturing fluid ρ g Liquid density of the pore fluid ρ fn Liquid density of the fracturing fluid K g Bulk modulus of the pore fluid K f r f Bulk modulus of the fracturing fluid )

Table 1 .
Duration of sequential fracturing and gas production stages.

Table 2 .
Duration of simultaneous fracturing and gas production stages.

Table 4 .
Computational parameters for mesh refinement and coarsening.

Table 5 .
Parameters of the naturally fractured model.