Hanyi Wang

Numerical Investigation of Fracture Spacing and Sequencing Effects on Multiple Hydraulic Fracture Interference and Coalescence in Brittle and Ductile Reservoir Rocks

Keywords: hydraulic fracture, fracture spacing, fracture sequencing, fracture interference, fracture coalescence, brittle, ductile, cohesive zone method (CZM), extended finite element method (XFEM)

For unconventional resources exploration and development, hydraulic fracture pattern and associated dimensions are critical in determining well stimulation efficiency and ultimate recovery. When creating arrays of hydraulic fractures along horizontal wells, stress field changes induced by hydraulic fractures themselves can lead to fracture interference and coalescence. The resulting complex fracture geometry may compromise or improve the effectiveness of the stimulation job, depending on the nature of the context. Currently, the prevailing approach for hydraulic fracture modeling also relies on Linear Elastic Fracture Mechanics (LEFM), which uses stress intensity factor at the fracture tip as fracture propagation criteria. Even though LEFM can predict hard rock hydraulic fracturing processes reasonably, but often fails to give accurate predictions of fracture geometry and propagation pressure in quasi-brittle and ductile rocks, such as poorly consolidated sands and clay-rich shales. In this study, a fully coupled hydraulic fracture propagation model based on the Extended Finite Element Method (XFEM), Cohesive Zone Method (CZM) and Mohr-Coulomb theory of plasticity is presented, to investigate the interference and coalescence of fluid-driven hydraulic fractures that initiated from horizontal wells. The results indicate that fracture spacing and the relative timing of fracture initiation control whether the fractures compete against each other to form a divergent pattern or coalesce into a single, primary fracture. Fracture growth can be arrested after fracture tips pass by when simultaneously fracturing adjacent horizontal wells. Even though the in-elastic rock deformation due to shear failure can strongly impact fracture geometry and fracturing pressure, it has limited influence on hydraulic fracture interaction patterns.

Introduction

The combination of horizontal drilling and multistage hydraulic fracturing technology has made possible the current flourishing oil and gas production from shale reservoirs in the United States, as well as the global fast growing investment in unconventional resources exploration and development. Accurate prediction of multiple fracture propagations and a comprehensive understanding of hydraulic fracture interference from a single well and between neighboring wells are vital in optimizing well placement and hydraulic fracturing treatment design (such as determine fracture spacing, fracturing sequence, fracture half-length, etc.). The hydraulic fracturing process is primarily aimed at well stimulation, but it provides other benefits as well. These include sand production control, connection of layered and laminated formations by vertical penetration, and afford the ability to spread depletion deep into the reservoir, thus significantly increasing economically recoverable reserves [START_REF] Economides | Reservoir Stimulation[END_REF].

When fracturing a horizontal wellbore, completions are distinguished by the methods of isolating fracturing section and diverting fracturing fluid to perforation tunnels. For open-hole completion, external casing packers are used as the isolation tool, and ports on sliding-sleeve are opened by dropped balls or coiled tubing. For cemented-casing completion, composite bridge plugs or internal casing packers are the isolation methods, and perforations or coiled tubing activated ports divert fracturing fluid. In these completions, starting from the toe of the horizontal well and working back to the vertical part of the well, a horizontal well is generally fractured multiple times, through a series of stages, to create multiple fractures to increase the reservoir contact area, as shown in Fig. 1. In each stage, the well is perforated in multiple perforation clusters with the aim of creating a separate fracture within each cluster. The number of clusters and the number of perforations within each cluster are selected based on the designed injection rate during fracturing. The process of dividing a horizontal well into different stages and propagating multiple fractures simultaneously from each cluster in a certain stage is often called "multistage fracturing".

Fig. 1 Illustration of multistage horizontal well fracturing

In most cases, the minimum stress is the minimum horizontal stress, resulting in vertical fractures that are perpendicular to that direction. Therefore, horizontal well trajectory presents two limiting cases (shown in Fig. 2): one drilled along the minimum stress, which would result in transverse fractures; and the other drilled at a right angle, which would result in longitudinal fractures. Both of these configurations may indicate multiple fracture treatments. Typically, any fracture that is more than approximately 10º away from the wellbore will behave as if it is transverse. Consequently, unless the wellbore azimuth is very tightly controlled, transverse fractures are the most likely scenario.

Fig. 2 Fracture configurations from a horizontal well: multiple transverse fractures and a longitudinal fracture Beside multistage hydraulic fracturing techniques along a single horizontal well, a number of different multi-well completion techniques have been proposed to improve stimulation efficacy, as shown in Fig. 3. The "simul-frac" technique is where companion wells are stimulated at the same location at the same time, whereas the "zipper-frac" technique employs companion wells that are stimulated in staggered locations at the same time. The intention with these techniques is to alter either or both the stress field and the pore pressure field to enhance the shear failure of weak planes and promote fracture complexity. and C, Modified zipper-frac [START_REF] Nagel | Quantitative Evaluation of Completion Techniques on Influencing Shale Fracture 'Complexity[END_REF] Traditionally, the hydraulic fracture is believed in a plane perpendicular to the least in-situ principal stress, and multiple fractures are parallel and separated from each other. With the assumptions of bi-wing fracture geometry and rock's tensile failure mechanism, analytical models are developed and widely used for hydraulic fracture design and analysis, such as the KGD model [START_REF] Khristianovich | Formation of vertical fractures by means of highly viscous fluids[END_REF][START_REF] Geertsma | A rapid method of predicting width and extent of hydraulic induced fractures[END_REF] and PKN model [START_REF] Perkins | Widths of hydraulic fractures[END_REF][START_REF] Nordgren | Propagation of vertical hydraulic fractures[END_REF]). However, under certain circumstances, the hydraulic fracture may evolve in a complex fashion due to various reasons, such as heterogeneous formation properties, intersection with natural fractures, and stress interference between multiple hydraulic fractures themselves. This, in turn, can significantly alter the relationship between the injection history and the fracture geometry that predicted by planar fracture models. So it is crucial to model non-planar hydraulic fracture propagation in order to understand its impact on the completion process and to ensure that undesirable situations do not arise.

In recent years, the extended finite element method (XFEM) has emerged as a powerful numerical procedure for the analysis of fracture problems. This method was developed to ease difficulties in solving problems with localized features that are not efficiently resolved by mesh refinement and help alleviate shortcomings of the finite element method and has been used to model the propagation of various discontinuities, such as cracks and material interfaces. A key advantage of XFEM is that in such problems the finite element mesh does not need to be updated to track the crack path, and discontinuous enrichment functions are used to approximate the displacement discontinuity across the crack surface, which significantly reduced the computational costs and projection errors associated with conventional finite element methods that restricting the discontinuities to mesh edges. Since the introduction of this method, the XFEM has been used to model non-planar hydraulic fracture propagations by many authors [START_REF] Lecampion | An Extended Finite Element Method for Hydraulic Fracture Problems[END_REF][START_REF] Dahi-Taleghani | Numerical Modeling of Multistranded-Hydraulic-Fracture Propagation : Accounting for the Interaction between Induced and Natural Fractures[END_REF][START_REF] Leonhart | Extended finite element method for hydro-mechanical analysis of crack propagation in porous materials[END_REF][START_REF] Gordeliy | Implicit level set schemes for modeling hydraulic fractures using the XFEM[END_REF].

However, an assumption commonly made in all these proposed XFEM models are that the loss of fluid into the rock is independent of fracture opening, so the pressure diffusion and porous behavior of the rock deformation are not fully coupled.

In addition, all these models are under the assumptions of linear elastic fracture mechanics (LEFM), which uses stress intensity factor at the fracture tip as fracture propagation criteria. Although hydraulic fracturing simulators based on LEFM can give reasonable predictions for hard rock formations, they fail to predict fracture net pressure and geometry with enough certainty in ductile formations, even in the case of simple, planar fracture geometry. Numerous study and surveys have indicated that the net pressure observed in the field in ductile formations are much higher than that predicted by LEFM and this disparity is even larger in poorly consolidated formations. These observations have triggered a series of dedicated studies which looked into the importance of the plastic deformation in hydraulic fracturing [START_REF] Papanastasiou | The influence of plasticity in hydraulic fracturing[END_REF][START_REF] Germanovich | Modeling multi-segmented hydraulic fracture in two extreme cases: No leak-off and dominating leak-off[END_REF][START_REF] Papanastasiou | The effective fracture toughness in hydraulic fracturing[END_REF][START_REF] Van Dam | Impact of rock plasticity on hydraulic fracture propagation and closure[END_REF]. Laboratory studies have also shown that deformation of very weakly cemented sands and clay-rich shales can even occur through elastic-visco-plastic constitutive behavior, and the deformations of these types of rocks cannot be predicted by linear elasticity [START_REF] Sone | Visco-plastic Properties of Shale Gas Reservoir Rocks[END_REF]. In addition, even for brittle material, stress intensity factor alone can't fully characterize fracture propagation in all cases [START_REF] Chitsiriphanit | On Determination of Apparent Fracture Toughness and Fracture Process Zone[END_REF].

In order to overcome the above limitations, cohesive zone method (CZM) has been adopted by many authors to model fracture initiation and propagation. The conception of the cohesive zone was first introduced by [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture: general ideas and hypothesis, axially symmetric cracks[END_REF][START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF] to investigate fracture propagation in perfectly brittle materials. In order to investigate fracture damage behavior in ductile materials with small scale of plasticity, a fracture process zone was proposed by [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF], that adopt a critical opening condition as a fracture criterion. The physical meaningless tip singularity predicted by LEFM can be resolved with the idea of cohesive zone, which is a region ahead of the crack tip that is characterized by micro-cracking along the crack path, and the main fracture is formed by the interconnection of these micro-cracks due to damage evolution. [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] introduced the concept of fracture energy into the cohesive crack model and proposed a number of traction-displacement relationships. Based on this conception, [START_REF] Mokryakov | Analytical solution for propagation of hydraulic fracture with Barenblatt's cohesive tip zone[END_REF] proposed an analytical solution for hydraulic fracture with Barenblatt's cohesive tip zone, the results demonstrate that the derived solutions from cohesive tip zone model can fit the pressure log much more accurately than LEFM for the case of fracturing soft rock. [START_REF] Wang | Poroelastic and Poroplastic Modeling of Hydraulic Fracturing in Brittle and Ductile Formations[END_REF] developed a cohesive poro-elasto-plastic hydraulic fracture model for both brittle and ductile rocks. Their work indicates that plastic damage during fracturing execution can lead to higher propagation pressure, shorter and wider fracture geometry. Sensitivity analysis also found that the effects of shear failure on fracture propagation are mostly controlled by in-situ stress, rock cohesion strength and reservoir pore pressure. In addition, their work reveals that the effects of plastic damage in reservoir rocks cannot be well represented by imposing an artificially increased toughness at the fracture tip (i.e., increase critical fracture energy inside the cohesive zone), especially when the shear failure does not constrain itself to the very vicinity of fracture tip. However, despite the successful implementations of CZM to model hydraulic fracture initiation and propagation, all these proposed models require a pre-define path for fracture propagation, which severely limits their applicability in predicting complex hydraulic fracture geometry. [START_REF] Wang | Numerical Modeling of Non-Planar Hydraulic Fracture Propagation in Brittle and Ductile Rocks using XFEM with Cohesive Zone Method[END_REF]2016a) proposed a CZM based XFEM model, to investigate fracture re-orientation, multi-fracturing process and production well interference. His works indicate that hydraulic fracture prefers to propagate in formations with lower Young's modulus if these formations behave elastically, and the explanation for the abnormally high treatment pressure encountered in some fields (such as clay-rich Eagle Ford and the Haynesville Shales) reveals that shear failures are actually happening during the execution of a stimulation job, and these effects induced by the in-elastic rock behaviors cannot be neglected [START_REF] Wang | Numerical Modeling of Non-Planar Hydraulic Fracture Propagation in Brittle and Ductile Rocks using XFEM with Cohesive Zone Method[END_REF].

It is also concluded that hydraulic fracture tends to propagate towards the pressure depleted region (lower pore pressure due to production), even though the effective stress is higher in the depleted zone (Wang 2016a).

In this study, by utilizing fully coupled CZM based XFEM hydraulic fracture propagation model, the initiation and propagation of multiple hydraulic fractures in horizontal wells are investigated. First, an introduction to the physical models and mathematical formulation will be presented. Then, synthetic simulation cases are constructed to examine what factors control multiple hydraulic fracture interaction patterns from a single horizontal well. And the following part will discuss the scenario where hydraulic fractures initiated from two adjacent wells. Finally, conclusion and discussion remarks are presented.

Mathematical Framework

Following fracture initiation, further fluid injection results in fracture propagation. the geometry of the created fracture can be approximated by taking into account the mechanical properties of the rock, the properties of the fracturing fluid, the conditions with which the fluid is injected (rate, pressure), and the stresses and stress distribution in the porous medium. In describing fracture propagation, which is a particularly complex phenomenon, two sets of laws are required: (i) Fundamental principles such as the laws of conservation of momentum, mass and energy, and (ii) criteria for propagation that include interactions of rock, fluid and energy distribution.

Fluid Flow

A wide variety of fluids have been used for fracturing including water, aqueous solutions of polymers with or without crosslinkers, gelled oils, viscoelastic surfactant solutions, foams, and emulsions. Many hydraulic fracturing fluids exhibit power law rheological behavior and temperature-related properties. In most low permeability reservoirs, slick water is often used a fracturing fluid. In order to avoid additional complexity added by fluid behavior, incompressible and Newtonian fluid is assumed in this study. Fig. 4 shows a sketch of fluid-driven hydraulic fracture with vary aperture.

𝒒 𝒇 = - 𝑤 3 12𝜇 ∇𝑝 𝑓 ( 1
)
where 𝑤 is the crack aperture, 𝑝 𝑓 is the fluid pressure inside the fracture and 𝜇 is the average fluid velocity over the crosssection of fracture. Pressure drop along the fracture can be determined by Eq. 1 with local flow rate and local fracture width.

The conservation of the fluid mass inside the fracture can be described by the Reynolds (lubrication) equation:

∇𝒒 𝒇 - ∂𝑤 ∂𝑡 + 𝑞 𝑙 = 0 (2)
where 𝑞 𝑙 is the local fluid loss in rock formation per unit fracture surface area. The local flow rate 𝑞 𝑓 can be determined by taking fluid leak-off into consideration. Pressure dependent leak-off model is used in this study to describe the normal flow from fracture into surrounding formations:

𝑞 𝑙 = 𝑐 𝑙 (𝑝 𝑓 -𝑝 𝑚 ) (3)
where 𝑝 𝑚 is pore pressure in the formation and 𝑐 𝑙 is fluid leak-off coefficient, which can be interpreted as the permeability of a finite layer of filtrate cake on the cracked surfaces. This leak-off model is suitable for low permeability, homogenous reservoirs. For practical applications, empirical Carter's Leak-off Model [START_REF] Valkó | Hydraulic Fracture Mechanics[END_REF] is normally used to represent the average leak-off behavior over the entire fracture surface, where the effects of variations in lithology among multiple layers and formation heterogeneity are included, by honoring field data from Fracture Calibration Tests (which is also often referred to as MiniFrac, Mini-Fall-off and Diagnostic Fracture Injection Test). [START_REF] Wang | Poroelastic and Poroplastic Modeling of Hydraulic Fracturing in Brittle and Ductile Formations[END_REF] have discussed how to implement Carter's Leak-off Model numerically.

Darcy's Law is used to describe fluid diffusion in the porous media:

𝒒 𝒎 = - 𝒌 𝜇 ∇𝑝 𝑚 (4)
where 𝒒 𝒎 is the fluid flux velocity vector in the porous media, 𝒌 is formation permeability tensor.

Coupled Deformation-Diffusion Phenomena

The basic theory of poroelasticity in which the fully coupled linear elastic rock deformation and pore pressure equations were initially introduced by the pioneering work of [START_REF] Biot | General theory of three dimensional consolidation[END_REF]. Since then, many researchers have contributed to its further development. In fluid-filled porous media, the total stresses σ 𝑖,𝑗 are related to the effective stresses 𝜎 𝑖,𝑗 ′ through:

σ 𝑖,𝑗 = 𝜎 𝑖,𝑗 ′ + 𝛼𝑝 𝑚 . ( 5
)
The effective stresses govern the deformation and failure of the rock, the poroelastic constant 𝛼 is a rock property that is independent of the fluid properties. [START_REF] Ghassemi | Theree Dimensional Poroelastic Hydraulic Fracture Simulation Using the Displacement Discontinuity Method[END_REF] demonstrated that variations in the value of poroelastic constant 𝛼 has negligible influence on fracture geometry. In this article, the poroelastic constant 𝛼 is assumed to be 1 and the equilibrium equation in the form of virtual work principle for the volume under its current configuration at time 𝑡 can be written as [START_REF] Zienkiewicz | The Finite Element Method[END_REF]:

∫ (𝜎 ′ + 𝑝 𝑚 𝐈) 𝑉 : 𝛿𝜀𝑑𝑉 = ∫ 𝒕𝛿𝑣𝑑𝑆 + 𝑆 ∫ 𝒇𝛿𝑣𝑑𝑉 𝑉 (6)
where 𝜎 ′ , 𝛿𝜀 are effective stress and virtual rate of deformation respectively. 𝒕 and 𝒇 are the surface traction per unit area and body force per unit volume, I is unit matrix. This equation is discretized using a Lagrangian formulation with displacements as the nodal variables. The porous medium is thus modeled by attaching the finite element mesh to the solid phase that allows liquid to flow through. A continuity equation required for the fluid, equating the rate of increase in fluid volume stored at a point to the rate of volume of fluid flowing into the point within the time increment:

𝑑 𝑑𝑡 (∫ 𝜌 𝑓 ∅𝑑𝑉 𝑉 ) + ∫ 𝜌 𝑓 ∅𝒏 𝑆 𝒒 𝒎 𝑑𝑆 = 0, (7)
where 𝜌 𝑓 and ∅ are the density of the fluid and porosity of the porous media respectively. 𝒏 is the outward normal to the surface S. The continuity equation is integrated in time using the backward Euler approximation and discretized with finite elements using pore pressure as the variable.

XFEM Approximation

The extended finite element method was first introduced by [START_REF] Belytschko | Elastic Crack Growth in Finite Elements with Minimal Remeshing[END_REF] to help alleviate shortcomings of the finite element method for modeling the propagation of various discontinuities. By introducing the concept of partition of unity by [START_REF] Melenk | The Partition of Unity Finite Element Method: Basic Theory and Applications[END_REF] and extend local constructions to the whole space under the conventional finite element method context, the presence of discontinuities can be presented by special enriched functions in conjunction with additional degrees of freedom. Thereby, it enables the accurate approximation of fields that involve jumps, kinks, singularities, and other non-smooth features within elements.

In the context of fracture analysis, the enrichment functions generally consist two parts: asymptotic functions that capture the singularity around the crack tip and a discontinuous function that depicts the displacement jump across the fracture surfaces.

The displacement vector u can be approximated with the partition of unity enrichment [START_REF] Fries | Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description[END_REF]:

𝒖 = ∑ 𝑁 𝐼 (𝑥)[𝒖 𝐼 + 𝐻(𝑥)𝒂 𝐼 + ∑ 𝐹 𝛼 4 𝛼=1 (𝑥)𝒃 𝐼 𝛼 ] 𝑁 𝐼=1 ( 8
)
where 𝑁 𝐼 (𝑥) are the usual nodal shape functions, which have a value of 1 at the node whose number it bears and zero at all other nodes. 𝒖 𝐼 , is the usual nodal displacement vector associated with the continuous part of the finite element solution, it is applicable to all the nodes in the model. 𝒂 𝐼 , is the nodal enriched degree of freedom vector and 𝐻(𝑥) is the associated discontinuous jump function across the crack surfaces. The product of 𝒂 𝐼 and 𝐻(𝑥) is valid for nodes whose shape function support is cut by the fracture interior. 𝒃 𝐼 𝛼 , is also the nodal enriched degree of freedom vector, but only exist at the fracture tip, and 𝐹 𝛼 (𝑥), is the associated elastic asymptotic crack-tip functions. The product of 𝒃 𝐼 𝛼 and 𝐹 𝛼 (𝑥) is only used in the situation for the nodes whose shape function support is cut by the fracture tip. Fig. 5 illustrates the discontinuous jump function across the crack surfaces

𝐻(𝑥) = { 1 if (x -x * ) • 𝐧 ≥ 0 -1 otherwise (9)
where x is a sample (Gauss) point, x * is the point on the crack closest to x , and 𝐧 is the unit outward normal to the crack at x * . In Fig. 2, 𝑟 and θ denote local polar coordinates system with its origin at the crack tip and θ = 0 is tangent to the crack at the tip. The asymptotic crack tip function, 𝐹 𝛼 (𝑥), can be determined by [START_REF] Lecampion | An Extended Finite Element Method for Hydraulic Fracture Problems[END_REF]:

𝐹 𝛼 (𝑥) = [√𝑟sin θ 2 , √𝑟cos θ 2 , √𝑟sinθ sin θ 2 , √𝑟sinθ cos θ 2 ] ( 10
)
The above functions can reproduce the asymptotic mode I and mode II displacement fields in LEFM, which represent the neartip singular behavior in strains and stresses. The use of asymptotic crack-tip functions is not restricted to crack modeling in an isotropic elastic material. The same approach can be used to represent a crack along a bi-material interface, impinged on the bimaterial interface, or in an elastic-plastic power law hardening material. However, in each of these three cases, different forms of asymptotic crack-tip functions are required depending on the crack location and the extent of the inelastic material deformation. The different forms for the asymptotic crack-tip functions are discussed by [START_REF] Sukumar | Partition of Unity Enrichment for Bimaterial Interface Cracks[END_REF], [START_REF] Sukumar | Modeling Quasi-Static Crack Growth with the Extended Finite Element Method Part I: Computer Implementation[END_REF], and [START_REF] Elguedj | Appropriate Extended Functions for X-FEM Simulation of Plastic Fracture Mechanics[END_REF], respectively. Accurately modeling the crack-tip singularity requires constantly keeping track of where the crack propagates and is cumbersome because the degree of crack singularity depends on the location of the crack in a non-isotropic material. However, this problem can be avoided if the moving cracks are modeled with the cohesive zone method, which will be discussed in the following section.

Cohesive Zone Method

The propagation criterion introduced within the context of LEFM assumes that the fracture process zone, a region near the fracture tip where behavior of the material in not elastic (e.g. region of plastic deformation, microcracking, etc.), is small compared to the fracture size, and a fracture can propagate only if the stress intensity factor exceeds the material toughness. However, fracture propagation in quasi-brittle and ductile materials can induce a significant plastic deformation around the fracture tip due to share failure, which put the adequacy of such assumption into questioning. Moreover, even for brittle materials, where the fracture process zone can be lumped into a single point, the presence of an initial crack is needed for LEFM to be applicable. So bodies with blunt notches, instead of sharp cracks, cannot be analyzed using LFEM. The cohesive zone model is a suitable and simple model for process zone, and for materials that fail by crack growth and coalescence, the appearance of the process zone in a cross section normal to the crack edge is shown in Fig. 6 (left). The process zone can be decomposed into elements and the behavior of each single element is governed by relationships between boundary forces and displacements. At any given time, the location of the process zone is characterized by two tips: a mathematical tip where the crack opening is zero and a physical tip where the cohesive tractions are zero, as shown Fig. 6 (right). Before any loading is applied to the specimen, these two tips coincide. As the loading is imposed, the fracture process zone develops. The description of the cohesive zone model will be discussed in two parts. In the first part, the constitutive equations characterizing the displacement and stress fields for a given location of the crack path are presented. In the second part, the conditions for crack initiation and growth are expressed. Based on loading conditions, a crack can generate according to three different failure modes, which are often named as Mode I (tension failure), Mode II (shear failure under sliding) and Mode III (shear failure under tearing). For mode-I plane strain fracture, the critical fracture energy G I c equals the area under the traction-separation curve, which can be related to the rock fracture toughness K 𝐼𝐶 by [START_REF] Kanninen | Advanced fracture mechanics[END_REF]:

G I c = K 𝐼𝐶 2 𝐸 (1 -𝜈 2 ) ( 11
)
where E is Young's modulus of formation and 𝜈 is Poisson's ratio. Damage is assumed to initiate when one of the stress components reaches the value of maximum strength of rock material in that direction, which can be represented by a quadratic law

{ 〈𝑡 𝑛 〉 𝑡 𝑛 0 } 2 + { 𝑡 𝑠 𝑡 𝑠 0 } 2 + { 𝑡 𝑡 𝑡 𝑡 0 } 2 = 1 ( 12
)
where t n , t s , t t refer to the normal, the first, and the second shear stress components; and t n 0 , t s 0 , t t 0 represent the tensile strength of the rock material when the deformation is purely perpendicular to the interface and the shear strength of rock material in the first and the second shear direction. The symbol 〈 〉 used in the above equation represents the Macaulay bracket with the usual interpretation.The Macaulay brackets are used to signify that a pure compressive deformation or stress state does not initiate damage.

The mode mix of the deformation fields in the cohesive zone quantify the relative proportions of normal and shear deformation. The Benzeggagh-Kenane fracture criterion (1996) is implemented to determine the mixed-mode damage evolution during fracture propagation. This criterion is suitable for the situation when the critical fracture energy of rock material along the first and the second shear directions are similar. The combined energy dissipated due to failure G c , is defined as

G c = G I c + (G II c -G I c ) ( G shear G total ) η ( 13
)
where

G shear = G II c + G III c , G total = G shear + G I c . And G I c , G II c , G III c
are the work done by the tractions and their conjugate relative displacements in the normal, first, and second shear directions. For an isotropic formation, the traction-separation responses in different modes are assumed to be the same (G II c = G I c ), so the cohesive response is insensitive to parameter η. The fracture will propagate when the energy release rate reaches the value of G c . And the newly introduced crack is always orthogonal to the maximum local tensile stress (maximum hoop stress) direction when the fracture criterion is satisfied.

Mohr-Coulomb Theory of Plasticity

An important motivation behind the introduction of a cohesive zone in fracture mechanics is to avoid the non-physical, singular stress fields at the crack tip. Even though CZM method can be used to capture the effects of microscopic in-elastic deformation in front fracture tip by introducing higher critical fracture energy [START_REF] Yao | Linear Elastic and Cohesive Fracture Analysis to Model Hydraulic Fracture in Brittle and Ductile Rocks[END_REF]), previous study by [START_REF] Wang | Poroelastic and Poroplastic Modeling of Hydraulic Fracturing in Brittle and Ductile Formations[END_REF] has demonstrated that the impact of plastic deformation on hydraulic fracture propagation cannot be well captured by imposing an artificially increased toughness at the fracture tip (i.e., increase the critical fracture energy in the cohesive zone), so CZM method by itself has limitations to fully represent the effects of plastic deformations, especially when shear induced plastic failures are not constrained to the very close vicinity of the fracture tip. In order to account for any potential shear damage in the whole simulation domain, Mohr-Coulomb criterion is implemented. Associative behavior with constant friction and dilatation angle is considered. These assumptions are justified by the presence of high confining stresses prior to crack propagation and to a decrease in the initial in-situ mean pressure near the crack tip during propagation. The Mohr-Coulomb criterion assumes that yield occurs when the shear stress on any plane in a material reaches the same value as shear strength, which is defined as:

𝜏 𝑓 = 𝑐 + 𝜎 𝑛 tanφ (14)
where 𝜏 𝑓 is shear strength, 𝑐 is cohesion strength, 𝜎 𝑛 is the stress normal to a specific plane, φ is the friction angle. And the ratio of rate of plastic volumetric strain to the rate of plastic shear distortion is controlled by dilation angle [START_REF] Vermeer | Non-Associated Plasticity for Soils[END_REF]. In essence, CZM serves as a vehicle for modeling fracture propagation based on local energy release rate to avoid stress singularity at the fracture tip, and the impact of shear damage is captured by the Mohr-Coulomb flow theory of plasticity within the entire simulation domain. Ductile shales and sandstones can have cohesion strengths that ranges between 3MPa and 8MPa, while the cohesion strength of poorly consolidated sands and shales with planes of weaknesses is normally below 3MPa [START_REF] Afrouz | Practical Handbook of Rock Mass Classification Systems and Modes of Ground Failure[END_REF]). The effective rock cohesion in a naturally fractured shale formation should be less than 1 MPa to justify a reasonable size of stimulated reservoir volume (SRV) induced by shear failure observed in the field [START_REF] Nassir | Geomechanical Coupled Modeling of Shear Fracturing in Non-Conventional Reservoir[END_REF].

Simulation Model for Single Well Fracturing

In order to investigate the propagating behaviors of multiple hydraulic fractures, a synthetic model is constructed as shown in Fig. 8. The horizontal well is drilling in the direction of minimum horizontal stress and two set of perforations are aligned in the direction of maximum horizontal stress, so that transverse hydraulic fractures can be created. The domain length of the simulation is 100 meters, which is large enough to minimize the impact of boundary effects. All the outer boundaries have zero displacement along the direction that perpendicular to its surface, and constant pore pressure condition is also imposed on the outer boundaries. The whole simulated domain is saturated with reservoir fluid and incompressible Newtonian fluid is injected at a constant rate at the perforation entries.

Fig. 8 Graphical representation of simulation model for two perforations

Besides modeling fracture propagation in brittle rocks that undergoes purely elastic deformation, fracture propagation in ductile rocks will be also included. Table 1 shows all the input parameters for the simulation cases, unless otherwise specified.

The coupled system of non-linear equations that presented in the previous section is solved numerically by Newton-Raphson method. All the variables in the system are updated at the end of each time increment and input as initial values at the start of the next increment. The program for numerical calculations was developed using FORTRAN and commercial finite element package.

Results and Analysis

First, the simulation results are compared with analytical solutions from KGD model [START_REF] Khristianovich | Formation of vertical fractures by means of highly viscous fluids[END_REF][START_REF] Geertsma | A rapid method of predicting width and extent of hydraulic induced fractures[END_REF], by configuring the model such that: the dimensions of the domain of analysis are much bigger than the fracture aperture and length, the permeability is defined to minimize the influence of poroelastic effects ahead of the fracture tip, cohesive properties are selected to ensure a small cohesive zone relative to the size of the fracture, and the rock behaves linear elastically. The results are shown in Fig. 9. it can be noticed that that the increase of fracture width at the wellbore agrees well with the asymptotic solution, and the pressure decline after breakdown is slightly higher than that of analytical solution. This is because the porous behavior of the medium surrounding the fracture and changes in the pore fluid pressure cannot be eliminated completely, which in turn affects the effective stress distribution. The oscillation of the simulation curves correspond to halts and sequels in the fracture propagation process step by step for the series of time increments. These fluctuations can be minimized by reducing the size of discretized mesh.

Fig. 9 Simulation results comparison between CZM based XFEM and KDG model

Next, the impact of fracturing spacing, and sequence of fracture propagation on fracture interference and coalescence are presented and discussed. Fig. 10 shows displacement magnitude and fracture geometry after 5 and 20 minutes of injection. The injection rates at both perforations are set to be the same, so both hydraulic fractures initiated and propagated simultaneously.

The results show that the hydraulic fractures tend to propagate away from each other, due to stress interference. When both fractures propagating at the same time and the fracture width are sustained by high pressure inside fracture, local distribution of stresses is significantly altered. The fractures orientated themselves to propagate in the direction of least resistance. It should be mentioned that the displace field presented in all the following simulation results are enlarged during post process, so a better visual observation of the fracture geometry can be obtained. Fig. 11 shows the formation pore pressure distribution, along with fracture geometry. After 5 minutes of injection, there's a high-pressure zone that concentrates around the two perforation sets. This is because the fluid leaks off from the fracture into the surrounding formations, and the induced deformations of pore space compress the pore fluid, where the permeability of rock matrix is not high enough to allow fluid to move freely at the same time scale. After 20 minutes of injection, it can be noticed that the high-pressure zone occurs at the fracture tips when fractures are propagating. This can be explained by the fact that when a complete new fracture surface is generated within an element at the fracture tip, there is a sudden fluid leak off into the adjacent formation, due to a large pressure difference between the fluid inside the fracture and the surrounding formation pore fluid. The leak off rate will gradually stabilize itself in a declining trend as the pressure difference becomes smaller over time. It is obvious that the simultaneously propagating fractures can strongly interference with each other because of the local stress disturbance; otherwise, the fractures would propagate along a straight line in the direction that perpendicular to the minimum horizontal stress. Fig. 12 shows the impact of fracture spacing on fracture interference. When the fracture spacing is 5 meters, other than 19 meters, the fracture turns more shapely. Compared to Fig. 10, it can be observed that the closer the fracture spacing, the most severe the fracture interference. Closer inspection also reveals that fracture width is narrower when fracture spacing is only 5 meters, the reduction in fracture width can adversely affect fracturing execution process and lead to suboptimal fracture geometry, because of proppant bridging and premature screen-out may occur from narrower than expected/designed fracture width. Besides the commonly known poorly consolidated/unconsolidated sands, shale rocks can also exhibit pronounced ductile behavior. Depending on mineral compositions (such as clay, carbonates, quartz, etc.), the ductility of shales poses a serious challenge in hydraulic fracturing design and execution. Field studies indicate that it is more difficult to initiate a hydraulic fracture in ductile Haynesville shales [START_REF] Parker | Haynesville Shale: Hydraulic Fracture Stimulation Approach[END_REF]), so it is crucial to incorporate in-elastic deformations into hydraulic fracturing models in these contexts. Next, the impact of shear damage on hydraulic fracture propagation behavior are investigated, by allowing formation rocks undergo plastic deformation when local stress state reaches the Mohr-Coulomb failure envelope. Fig. 13 shows the plastic deformation area and fracture geometry. It can be observed that the plastic damage occurs during fracture propagation and it does not prevent fracture interference of simultaneously propagating fractures. Fig. 14 shows the evolution of fracture width and pressure at the wellbore in brittle (with poro-elastic effects) and ductile (with poro-elastic-plastic effects) rocks, the results indicate that higher pressure is needed to fracturing ductile rocks and larger fracture width can be created in ductile rocks. This implies that longer injection time and larger proppant size may be more suitable for clay-rich ductile shales when optimizing hydraulic fracture design and execution.

Fig. 13 Plastic damage and fracture geometry with fracture spacing=9m

Fig. 14 Fracture width and pressure at wellbore with and without shear damage

In previous simulation cases, the hydraulic fractures are initiated and propagated at the same time with the same injection rate in both perforation sets. In reality, the hydraulic fractures may not initiate and propagate simultaneously, due to formation heterogeneity, perforations flaws or the existence of weak planes and natural fractures that intersect perforation tunnels. In general, fracture initiation pressure from different perforations depends on their location around the wellbore, and the shape of the perforation itself. Perforation breakdown will be sequential and related to wellbore pressure, moving from the easiest to break (least pressure) to more difficult. In actual operations, each stage of the horizontal well contains multiple perforation clusters. These clusters compete with each other for the injected slurry. And the variable mode of fracture initiation and extension from individual clusters makes the process somewhat unstable. In order to examine the interactions of sequential propagating fractures, the injection rate at each perforation set are modified, as shown in Fig. 15. For the right perforation set, the injection rate increases gradual to the target rate in 50 seconds, and for the left perforation set, the injection rate reaches the target level after 240 seconds. Because the right perforation set receives more injection fluid initially, the associated fracture will initiate and propagate earlier. Fig. 16 shows the fracture geometry at two different injection time frames when fracturing spacing is 19 meters. As expected, the hydraulic fracture initiated and propagated earlier from the right perforation set, which leads to longer fracture geometry at the same injection time. The effects of fracture interference still exist, and the fractures propagated away from each other to minimize the fracturing resistance. However, if the fracturing spacing is reduced, the results can be completely different, as shown in Fig. 17 and Fig. 18. When fracture spacing is 9 meters, the hydraulic fracture initiated and propagated first form the right perforation tunnel, but both fractures have the propensity to propagate parallel to each other. In addition, the fracture width of the right hydraulic at the section that closes to the wellbore is significantly reduced, due to the impact of the left hydraulic fracture, which initiated and propagated later. When fracture spacing is only 5 meters, the left hydraulic fracture tends to coalesce into the right hydraulic fracture. These results agree well with field observations that when hydraulic fractures initiate and propagate from a perforation cluster, where multiple perforation tunnels are created alongside each other within a small interval (around 2 ft distance), only a one dominating hydraulic fracture will be formed away from the wellbore, because even small perturbations of the fluid received in each perforation tunnels will lead to fracture coalescence. So in closely spaced perforations, a fracture initiating from one perforation may make a minor deviation from its path and connect with another perforation. In contrast, when hydraulic fractures initiate and propagate from multiple clusters, where the space between clusters are normally around 250 to 500 ft, hydraulic fractures always push against each other during propagation. The simulation results also agree with laboratory experiment. Fig. 20 shows the fracture fluid pressure and final geometry for four closely spaced hydraulic fractures that propagated sequentially inside a cubic block from experiment [START_REF] Kear | Three Dimensional Forms of Closely-Spaced Hydraulic Fractures[END_REF]. It indicates that the fracturing pressure follows a KGD trend after breakdown and fractures tend to merge into each other at a distance away from the injection entries. So from the above analysis, it is clear that fracture spacing and fracture initiation sequence have a profound impact on the hydraulic fracture interactions. Depending on the formation properties, the numbers of perforation clusters per stage, the spacing between perforation clusters and the fracturing operational procedures, the prorogating hydraulic fractures can strongly interference with each other, to possibly undermine or strength our stimulation efficiency.

Fig. 20

Fracture pressure records and cross section of experiment block [START_REF] Kear | Three Dimensional Forms of Closely-Spaced Hydraulic Fractures[END_REF]

Simultaneously Fracturing Multiple Wells

In previous sections, we only investigated hydraulic fracture initiation and propagation from a single horizontal well. In the following part, simultaneous fracturing two horizontal wells are examined. Simultaneous multiple fracturing (Simul-frac) is a hydraulic fracturing technology by fracturing two or more parallel horizontal wells simultaneously to improve oil and gas production. Simul-frac has been successfully used in shale developments and both initial production rates and estimated ultimate recovery are significantly improved over a single well for most simul-frac wells [START_REF] Mutalik | Case History of Sequential and Simultaneous Fracturing of the Barnett Shale in Parker Country[END_REF]. Although field practices have yielded a significantly improvement from simul-frac over stand-along wells, and microseismic data also show a more complex fracturing network [START_REF] Waters | Simultaneous Hydraulic Fracturing of Adjacent Wells in the Woodford Shale[END_REF], the reasons for simul-frac success are still not clear. In order to investigate fracture interactions of simultaneous fracturing from two horizontal wells, a model scenario is constructed, as shown in Fig. 21. Two hydraulic fractures initiate and propagate from their respective perforation set at the top and bottom horizontal well. The fracture offset between these two perforations is 9 m. All the input parameters and boundary conditions are the same as previous sections. Fig. 22 shows the fracture geometry at two different injection time frames. It can be observed that after 5 minutes of injection, both fractures continue propagating along a straight line that perpendicular to the direction of minimum horizontal stress. However, when fracture tip passes by each other after 30 minutes of injection, these two fractures intertwined with each other.

With closer inspection, we can notice that the fracture tips re-oriented themselves to be perpendicular to the main fracture path, which arrested the fracture growth and prevents it from penetrating deeper into the formation, which results in wider fracture geometry as a consequence of continuous injection of fluid.

Fig. 22

Displacement magnitude and fracture geometry with fracture offset=9m

Fig. 23 shows the direction of minimum horizontal stress map (the red line segment denotes the direction) after 5 and 15 minutes of injection. Initially, the direction of minimum horizontal stress is along the x-direction, however, as fracture propagating, it alters the local in-situ stresses around it, and the direction of minimum horizontal stress can be even reversed when close to the fracture surface. When hydraulic fractures propagate into the "stress alteration zone" that induced by the other fracture, fracture reorients itself to propagate perpendicular to the new direction of local minimum horizontal stress. Fig. 24 shows the plastic strain and fracture geometry when shear damage is considered. The fracture growth follows the same trend as in above elastic case. However, after 10 minutes of injection, the shear damage zone start to overlap and connect to each other and finally, the fracture tips reoriented and ceased to propagate deeper into the opposite wells. Unlike hydraulic fractures initiate and propagate from a single horizontal well, hydraulic fractures in simultaneous fracturing wells will deflect from the direction of far-field maximum principal stress and attract each other, which can prevent main fractures pass through the formation into opposite wells. It is possible that the plastic damage induced by shear failure can alter the local permeability in the hydraulic fracture stimulated volume, if the shear damage mostly results from natural fractures or pre-existing weak planes. Even though this presented model is able to capture the main hydraulic fracture evolution explicitly and the relate local permeability variations to shear failure mechanisms by using the concept of effective cohesion strength [START_REF] Nassir | Geomechanical Coupled Modeling of Shear Fracturing in Non-Conventional Reservoir[END_REF]) in stimulated reservoir volume (SRV), but how to quantify the rock effective permeability altered by the propagating hydraulic fracture in the stimulated regime is still poorly understood, and currently, the only way to measure the overall rock effective permeability after hydraulic fracturing is through Fracture Calibration Tests (FCT). In addition, both the rock matrix permeability and induced fracture conductivity can be pressure and stress dependent and vary with time and space [START_REF] Wang | Impact of Shale-Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption, and Geomechanics[END_REF][START_REF] Wang | What Factors Control Shale Gas Production Decline Trend: A Comprehensive Analysis and Investigation[END_REF], so investigate the shear failure induced permeability change inside the associated SRV opens up a completely new issue and hence, will not be discussed further. But if shear damage does increase local permeability by opening up natural fracture, then fracture interactions between adjacent horizontal wells can indeed improve fracture complexity, enhances reservoir contact.

Conclusions and discussions

In-situ stress field has dominant influence on the creation of a hydraulic fracture, and the fracture follows the local path of least resistance, not the global path. The propagating fracture can alter the stress field in the region surrounding the fracture, and under the circumstances of multiple propagating fractures, undesirable fracture interactions may lead to ineffective well stimulation treatments, which increases well cost, reduces production and hampers well economics. So an accurate understanding of the interaction effects between multiple hydraulic fractures is crucial to optimize completion designs and reservoir development. In this study, a fully coupled hydraulic fracture propagation model based on the Extended Finite Element Method (XFEM), Cohesive Zone Method (CZM) and Mohr-Coulomb theory of plasticity is presented, which is able to capture the development of complex hydraulic fracturing propagation in both brittle and ductile rocks. The phenomena of fracture interference and coalescence from multistage fracturing and simultaneous fracturing are investigated. The simulation results agree well with field observations and laboratory experiments. The main conclusions from the study include:

1. When multiple fractures propagation from a single horizontal well, fracture spacing and the sequence of fracture initiation control the interactions between multiple propagating fractures. When fractures initiate and propagate simultaneously, fractures always deviate from each other regardless of fracture spacing. When fracture initiate and propagate sequentially, large fracture spacing lead to competing fractures while small fracture spacing (can be on the order of meters) result in fracture coalescence.

2. When multiple fractures propagate simultaneously from parallel horizontal wells, the fracture tips tend to re-orient themselves to be perpendicular to the main fracture after passing by each other. This lead to the development of complex fracture geometry and prevent fractures from penetrating deeper into the direction of opposite wells.

3. Even though plastic damage due to shear failure can strongly impact fracture length, width and fracturing pressure, it has limited influence on fracture interaction patterns.

Compared to brittle rocks, ductile rocks are more difficult to breakdown and higher pressure is needed for fracture propagation. The locations of ductile rocks can significant impact the evolution of hydraulic fracture geometry. In a typical sandstone or carbonate reservoirs, shale layers are normally served as barriers to contain the growth of fracture height, due to their larger Poisson's ratio and higher in-situ horizontal stresses. In shale reservoir, fracture height growth can also be contained due to its laminated properties. For example, shales that contain more content of clay tend to be more ductile and it can constrain fracture to propagate only within adjacent brittle shale layers. If ductile and brittle rocks are distributed alternatively along the direction of horizontal wellbore, then well stimulation efficiency can be significantly impacted. For example, in one stage of fracturing, if some of perforation clusters located in the ductile section while others located in brittle section, then it is often the case that fractures only initiate and propagate from clusters in the brittle section, where less pressure is needed. At the end, fractures only grow inside brittle rocks and the reservoir is not evenly stimulated as designed.

In some worse cases, the fracture may even interfere other nearby wells because of much longer penetration length, due to the less than expected number of simultaneous propagating fractures (while fluid injection volume remains the same). So a better understanding of formation rock properties and its heterogeneities on appropriate scales is crucial for hydraulic fracture modeling and the optimization design of well completion in unconventional reservoirs.

Fig. 3

 3 Fig. 3 Common shale completion schemes-A, Simultaneous hydraulic fracturing (simul-frac); B, Sequential fracturing (zipper-frac);and C, Modified zipper-frac[START_REF] Nagel | Quantitative Evaluation of Completion Techniques on Influencing Shale Fracture 'Complexity[END_REF]

Fig. 4

 4 Fig. 4 Sketch of a plane-view hydraulic fracture with varying aperture For a flow between parallel plates, local tangential flow rate 𝒒 𝒇 can be determined by the pressure gradient to the fracture width for a Newtonian fluid of viscosity μ (Boone and Ingraffea 1990):

Fig. 5

 5 Fig. 5 Illustration of normal and tangential coordinates for a smooth crack And the discontinuous jump function has the following form:

Fig. 6 Fig. 7

 67 Fig. 6 Illustration of cohesive zone model for fracture process zone The constitutive behavior of the cohesive zone is defined by the traction-separation relation, as shown in Fig.7. It assumes that the material damage imitated when traction reaches the tensile strength/shear strength in an enriched element. Beyond that, the traction reduces linearly to zero up to the displacement of complete failure, and any unloading takes place irreversibly.

Fig. 10

 10 Fig. 10 Displacement magnitude and fracture geometry with fracture spacing=9m

Fig. 11

 11 Fig. 11 Pore pressure and fracture geometry with fracture spacing=9m

Fig. 12

 12 Fig. 12 Displacement magnitude and fracture geometry after 20 minutes of injection with fracture spacing=5m and 19m

Fig. 15

 15 Fig. 15Fracturing fluid injection rate at perforations

Fig. 16 Fig. 18

 1618 Fig. 16 Displacement magnitude and fracture geometry with fracture spacing=19m

Fig. 19

 19 Fig. 19 Plastic damage and fracture geometry with fracture spacing=5m

Fig. 21

 21 Fig. 21 Graphical representation of simulation model for simul-frac, fracture offset=9m

Fig. 23

 23 Fig. 23 Direction of minimum horizontal stress with fracture offset=9m

Fig. 24

 24 Fig. 24 Shear damage and fracture geometry after at different time frames with fracture offset=9m

Table 1

 1 Input parameters

	Input Parameters	Value
	Fracture spacing	9 m
	Perforation length	0.5 m
	Elastic modulus	15 GPa
	Poisson's ratio	0.25
	Fluid viscosity	1 cp
	Tensile strength	1 MPa
	Formation effective permeability	1 md
	Injection rate per unit reservoir thickness	1.0E-3 (m 3 /s/m)
	Specific weight of fluid	9.8 kN/m 3
	Initial pore pressure	10 MPa
	Maximum horizontal stress	25 MPa
	Minimum horizontal stress	20 MPa
	Vertical stress	30 MPa
	Critical fracture energy	100 J/m 2
	Pressure dependent leak-off coefficient	4E-14 m 3 /s/Pa
	Porosity	0.1
	Friction angle	30 o
	Dilation angle	10 o
	Cohesion strength	2 MPa