Petroleum Hanyi Wang

Poro-Elasto-Plastic Modeling of Complex Hydraulic Fracture Propagation: Simultaneous Multi-Fracturing and Producing Well Interference*

Keywords: hydraulic fracturing, poro-elasto-plastic, extended finite element (XFEM), cohesive zone method (CZM), Mohr-Coulomb plasticity, brittle, ductile, multiple fracturing, well interference, complex fracture

With the increasing wide use of hydraulic fractures in the petroleum industry, it is essential to accurately predict the behavior of fractures based on the understanding of fundamental mechanisms governing the process. For effective reservoir exploration and development, hydraulic fracture pattern, geometry and associated dimensions are critical in determining well stimulation efficiency. In shale formations, non-planar, complex hydraulic fractures are often observed, due to the activation of preexisting natural fractures and the interaction between multiple, simultaneously propagating hydraulic fractures. The propagating of turning non-planar fractures due to the interference of nearby producing wells has also been reported. Current numerical simulation of hydraulic fracturing generally assumes planar crack geometry and weak coupling behavior, which severely limits the applicability of these methods in predicting fracture propagation under complex subsurface conditions. In addition, the prevailing approach for hydraulic fracture modeling also relies on Linear Elastic Fracture Mechanics (LEFM) by assuming the rock behaves pure elastically until complete failure. 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 ductile and quasi-brittle rocks, such as poorly consolidated/unconsolidated sands and ductile shales, even in the form of simple planar geometry. In this study, we present a fully coupled poro-elasto-plastic model for hydraulic fracture propagation based on the theories of extend finite element (XFEM), cohesive zone method (CZM) and Mohr-Coulomb plasticity, which is able to capture complex hydraulic fracture geometry and plastic deformations in reservoir rocks explicitly. To illustrate the capabilities of the model, example simulations are presented including ones involving simultaneously propagating multiple hydraulic fractures and producing well interference. The results indicate that both stress shadow effects and producing well interference can alter hydraulic fracture propagation behavior substantially, and plastic deformations in ductile reservoir rocks can indeed make a significant difference in fracturing pressure and final fracture geometry.

Introduction

Hydraulic fracturing has been widely used as a common practice to enhance the recovery of hydrocarbons from low permeability reservoirs and prevent sand production in high permeability reservoirs [START_REF] Economides | Reservoir Stimulation[END_REF]. Hydraulic fracturing, which creates extensive and conductive fracture path that extended from wellbore to the deep formation, combine with horizontal drilling, allow formerly uneconomical unconventional reservoirs to be commercially viable. Understanding hydraulic fracture initiation and propagation from wellbores is essential for performing efficient hydraulic fracture stimulation design and treatment.

As a reservoir stimulation tool to enhance production, the problem of hydraulic fracturing is in essence one of predicting the shape of the fracture as a function of time, given the fluid pressure at the wellbore or the flow rate into the fracture. Even today, modeling fluid-driven fracture propagation is still a challenging problem, due to its moving boundaries and strong coupling effects among different mechanisms. The mathematical formulation of the problem is represented by a set of nonlinear integro-differential equations, with moving boundaries where the governing equations degenerate and become singular. The complexity of the problem often restricts researchers to consider only simple fracture geometries, 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, non-planar fashion due to various reasons, such as heterogeneous formation properties, intersection with natural fractures, initiated at an unfavorable orientation and stress interference with other hydraulic fractures. This in turn can significantly alter the relationship between the injection history and the crack 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/or 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 [START_REF] Karihaloo | Modelling of Stationary and Growing Cracks in FE Framework without Remeshing: A State-of-the-Art Review[END_REF]. Since the introduction of this method, the XFEM has been used to model complex, nonplanar 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 and the mechanics of fracture opening are considered independent and their interaction is ignored, so the pressure diffusion and porous behavior of the rock deformation are not fully coupled. In addition, all these models are under the assumption 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 formation, 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 of 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] Papanastasiou | The effective fracture toughness in hydraulic fracturing[END_REF][START_REF] Germanovich | Modeling multisegmented hydraulic fracture in two extreme cases: No leak-off and dominating leak-off[END_REF][START_REF] Van Dam | Impact of rock plasticity on hydraulic fracture propagation and closure[END_REF]. Besides the commonly known poorly consolidated/unconsolidated sands, shale rocks can also exhibit in-elastic behavior. Laboratory studies have also shown that deformation of weakly cemented sands and 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, field study also indicates that it is more difficult to initiate and propagate a hydraulic fracture in ductile shales [START_REF] Parker | Haynesville Shale: Hydraulic Fracture Stimulation Approach[END_REF].

In order to model hydraulic fracture growth in quasi-brittle and ductile formations, cohesive zone method (CZM) has been adopted by many authors to model fracture initiation and propagation. The conception of cohesive zone was first introduced by Barenblatt [17,[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 Dugdale [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 interconnection of these micro cracks due to damage evolution. These proposed cohesive zone models assumed that the cohesive zone length and the cohesive stress distribution are material parameters that do not depend on fracture dimensions and loading. The classical Griffith's brittle elastic fracture model [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF][START_REF] Griffith | The theory of rupture[END_REF] is a limited case of cohesive zone model: if the cohesive zone length tends to zero (with respective increase of the cohesive stress) then, in the limit, the stress-strain state corresponding to the classic square root asymptotes will be obtained.

Based on this conception, Mokryakov [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 zone model can fit the pressure log much more accurately than LEFM for the case of fracturing soft rock. Yao [START_REF] Yao | Linear Elastic and Cohesive Fracture Analysis to Model Hydraulic Fracture in Brittle and Ductile Rocks[END_REF] developed a 3D cohesive zone model to predict fracture propagation in brittle and ductile rocks and the effective fracture toughness method was proposed to consider the fracture process zone effect on the ductile rock fracture. The results show that in ductile formations, the cohesive zone model gives predictions that are more conservative on fracture length as compared with pseudo 3D and PKN models. Wang et al. [START_REF] Wang | Poroelastic versus Poroplastic Modeling of Hydraulic Fracturing[END_REF] developed a hydraulic fracture model for both brittle and ductile rock fracturing that combines the cohesive zone method and Mohr-Coulomb flow theory of plasticity. Unlike previous studies, their model not only considers in-elastic behavior at the fracture tip inside the cohesive zone, but also captures plastic deformations in the bulk formation. Their work indicates that plastic deformations during fracturing execution can lead to higher initiation, propagation pressure and shorter, wider fracture geometry. It also found that effects of formation plasticity on fracture propagation are mostly controlled by in-situ stress, rock cohesion strength and pore pressure, and the effects plastic deformations in the bulk formation cannot be fully represented by the mothed of using effective fracture toughness. However, all these CZM models proposed in literature require a pre-define path for fracture propagation, which severely limits the applicability of these methods in predicting complex fracture geometry.

Besides CZM, Continue Damage Model (CDM) has been proposed to capture damage and plastic deformations in reservoir rocks during hydraulic fracture propagation. Busetti et al. [START_REF] Busetti | Damage and plastic deformation of reservoir rocks: Part 1. Damage fracturing[END_REF][START_REF] Busetti | Damage and plastic deformation of reservoir rocks: Part 2. Propagation of a hydraulic fracture[END_REF] proposed an elastic-plastic continue damage model for hydraulic fracture propagation, branching and segmentation. A local damage variable that ranges from 0 (intact rock) to 1(complete failure) is used as an indicator of fracture propagation path. Even though CDM has the capacity to model shear and tensile damages induced by complex fracture propagation, it fails to capture the fracture width and length explicitly, and the simulation results are highly mesh size dependent, because the fracture morphology is only represented by damage variables that smeared across the discretized simulation domain. Without knowing the fracture width and length explicitly, it is not possible to determine the appropriate proppant size, proppant concentration, fluid efficiency and slurry pumping schedule, which renders this method unattractive in practical hydraulic fracturing design and post-fracture production prediction.

In this study, a fully coupled poro-elasto-plastic complex hydraulic fracturing model that based on XFEM and CZM is presented, which is not only be able to capture non-planar fracture geometry evolution explicitly, but also be able to capture plastic deformations in the reservoir rocks. Finite element package ABAQUS is used as a platform to develop the numerical simulation program. The phenomenon of stress shadow effects between simultaneously propagating fractures and producing well interference are investigated with this presented model in both brittle and ductile formations.

Mathematical Frameworks

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 order to avoid additional complexity added by fluid behavior, incompressible and Newtonian fluid is assumed in this study. Fig. 1 shows a sketch of fluid-driven hydraulic fracture with vary 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 μ [START_REF] Boone | A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media[END_REF]:

𝒒 𝒇 = - 𝑤 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 adjacent 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 one-dimensional that defines the leak-off rate perpendicular to the fracture surface, so the leak-off direction changes with the orientation of the fracture surface.

In addition, the fluid diffusion process and rock deformations are fully coupled. 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 vector.

Coupled Deformation-Diffusion Phenomena

The basic theory of poroelasticity in which the fully coupled linear elastic rock deformation and pore pressure equations was initially introduced by the pioneering work of Biot [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. Ghassemi [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 study, 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:

∫ (𝜎 ′ + 𝑝 𝑚 𝐈) 𝑉 : 𝛿𝜀𝑑𝑉 = ∫ 𝒕𝛿𝑣𝑑𝑆 + 𝑆 ∫ 𝒇𝛿𝑣𝑑𝑉 𝑉 ( 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.

The inelastic rock material behavior follows the Mohr-Coulomb flow theory of plasticity for a cohesive frictional dilatant material. Associative behavior with constant 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φ (8)
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, Concrete and Rock[END_REF].

XFEM Approximation

The extended finite element method was first introduced by Belytschko and Black [START_REF] Belytschko | Elastic Crack Growth in Finite Elements with Minimal Remeshing[END_REF] to help alleviate shortcomings of the finite element method in modeling the propagation of various discontinuities. By introducing the concept of partition of unity by Melenk and Babuska [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 nonsmooth 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 ( 9
)
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. 2 illustrates the discontinuous jump function across the crack surfaces Fig. 2 Illustration of normal and tangential coordinates for a smooth crack And the discontinuous jump function has the following form:

𝐻(𝑥) = { 1 if (x -x * ) • 𝐧 ≥ 0 -1 otherwise ( 10
)
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 ] ( 11
)
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 bimaterial interface, impinged on the bimaterial interface, or in an elastic-plastic material. However, in each of these three cases different forms of asymptotic cracktip 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 have been discussed by many authors [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][START_REF] Elguedj | Appropriate Extended Functions for X-FEM Simulation of Plastic Fracture Mechanics[END_REF]. However, 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 next section.

Cohesive Zone Method

There is a number of fracture propagation criteria have been proposed by previous studies, which is usually described by either a stress condition or an energy condition. The propagation criterion introduced within the context of LEFM assumes that the 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 ductile formations can induce a significant plastic deformation around the fracture due to share failure, which put the adequacy of such assumption into questioning. The cohesive zone model that takes the process zone into consideration is able to capture non-linear fracture mechanics behavior based on energy condition, and a fracture will propagate when the energy release rate in the process zone reaches the critical fracture energy, and the critical fracture energy can be estimated by correlations or laboratory experiment.

The constitutive behavior of the cohesive zone is defined by the traction-separation relation, which includes the stage of initial loading, initial damage, and the damage evolution that leading to the final failure at the bounded interface. The behavior of the interface prior to initiation of damage can be described as linear elastic but with a penalty stiffness that degrades under tensile and shear loading if the corresponding stress reaches a critical value, but the stiffness is not affected by pure compression. Laboratory experiments can be used to derive these relations by investigating post-peak behavior with principle/shear stress and axial/shear strain data.

In this study, bilinear cohesive law [START_REF] Tomar | Bounds for element size in a variable stiffness cohesive finite element model[END_REF] is used in our model, as shown in Fig. 3. where t n , t s , t t refer to the nomal, 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; δ n 0 , δ s 0 , δ t 0 corespond to the displacement of initial damage in the normal, the first, and the second shear stress direction and δ n f , δ s f , δ t f are the displacement of complete failure in these three directions. It assumes that the material exhibit linear elastic behavior before the traction reaches the tensile strength/shear strength or the separation of cohesive surfaces exceeds the displacement of damage initiation. Beyond that, the traction reduces linearly to zero up to the displacement of complete failure, and any unloading takes place irreversibly.

Fig. 3 Linear Traction-Separation law for different modes

Based on loading conditions, a crack can generate according 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 ) ( 12
)
where 𝐸 is Young's modulus of formation and 𝜈 is Poisson's ratio. Alternative, the critical fracture energy can be determined by laboratory experiment. Typically, conventional hydraulic fracture models based on LEFM employ Mode-I based fracture criteria that only consider tensile failure mechanism in fracture propagation. However, shear failure can play an important role in ductile formations under certain loading conditions. In this study, the combined effects of different modes will be used to define the damage initial and the propagation criteria in the following discussions.

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 ( 13
)
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 stress components of the traction-separation model are affected by the damage according to

𝐭 = { (1 -𝐷)𝒕 ̅ damage initated 𝒕 ̅ no damage occurs ( 14
)
where 𝐭 are stress components, 𝒕 ̅ are the stress components predicted by the elastic traction-separation behavior for the current strains without damage. D is a scalar damage variable, which has an initial value of 0 and monotonically increases to 1 as damage developing, represents the overall damage that comes from the combined effects of different traction-separation modes in the rock material. For linear softening as shown in Fig. 3, the evolution of the damage variable, D, reduces to [START_REF] Turon | A damage model for the simulation of delamination in advanced composites under variable-model loading[END_REF]:

D = 𝛿 𝑚 𝑓 (𝛿 𝑚 𝑚𝑎𝑥 -𝛿 𝑚 0 ) 𝛿 𝑚 𝑚𝑎𝑥 (𝛿 𝑚 𝑓 -𝛿 𝑚 0 ) ( 15
)
where 𝛿 𝑚 is effective displacement, defined as

𝛿 𝑚 = √〈𝛿 𝑛 〉 2 + 𝛿 𝑠 2 + 𝛿 𝑡 2 ( 16
)
The mode mix of the deformation fields in the cohesive zone quantify the relative proportions of normal and shear deformation. The the Benzeggagh-Kenane fracture criterion [START_REF] Benzeggagh | Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus[END_REF] is implemented to determine the mixed-mode damage evolution during fracture propagation. This criterion is suitable for 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 ) η ( 17
)
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 in this study, where 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 direction when the fracture criterion is satisfied. Unlike these methods in previous studies [START_REF] Mokryakov | Analytical solution for propagation of hydraulic fracture with Barenblatt's cohesive tip zone[END_REF][START_REF] Yao | Linear Elastic and Cohesive Fracture Analysis to Model Hydraulic Fracture in Brittle and Ductile Rocks[END_REF][START_REF] Wang | Poroelastic versus Poroplastic Modeling of Hydraulic Fracturing[END_REF], which require that the cohesive zone surfaces align with element boundaries and the cracks propagate along a set of predefined paths, The combination of XFEM and cohesive zone method can model fracture initiation and propagation in any arbitrary path in both brittle and ductile material, by taking the advantages of both XFEM and cohesive zone methods. Under such combination, the near tip asymptotic singularity is replaced by cohesive zone and only the displacement jump across an entire fractured element is considered. Therefore, the crack has to propagate across an entire element at a time to avoid the need to model the stress singularity.

Phantom nodes, which are superposed on the original real nodes, are introduced to represent the discontinuity of the cracked elements [START_REF] Rabczuk | A new crack tip element for the phantom-node method with arbitrary cohesive cracks[END_REF], as illustrated in Fig. 4. When the element is intact, each phantom node is completely constrained to its corresponding real node. When the element is cut through by a crack, the cracked element splits into two parts. Each part is formed by a combination of some real and phantom nodes depending on the orientation of the crack. Each phantom node and its corresponding real node are no longer tied together and can move apart. Fig. 4 The principle of the phantom node method [START_REF] Rabczuk | A new crack tip element for the phantom-node method with arbitrary cohesive cracks[END_REF] The magnitude of the separation is governed by the cohesive law until the cohesive strength of the cracked element is zero, after which the phantom and the real nodes move independently. To have a set of full interpolation bases, the part of the cracked element that belongs in the real domain, Ω 0 , is extended to the phantom domain, Ω 𝑝 ,. Then the displacement in the real domain can be interpolated by using the degrees of freedom for the nodes in the phantom domain. The jump in the displacement field is realized by simply integrating only over the area from the side of the real nodes up to the crack; i.e., Ω 0 + and Ω 0 -. This method provides an effective and attractive engineering approach and has been used for simulation of the initiation and growth of multiple cracks in solids [START_REF] Song | A Method for Dynamic Crack and Shear Band Propagation with Phantom Nodes[END_REF][START_REF] Remmers | The Simulation of Dynamic Crack Propagation using the Cohesive Segments Method[END_REF]. It has been proven to exhibit almost no mesh dependence if the mesh is sufficiently refined.

3 Numerical Simulation (the section of near wellbore scale simulation in the original manuscript is replaced with this new section that contains two field scale ones) In this section, two field scale cases are investigated and discussed to demonstrate the capability of the presented poro-elastoplastic hydraulic fracture model in both brittle and ductile formations. The first simulation scenario involves multi-fracturing in horizontal wells, where simultaneous propagating hydraulic fractures interfere with each other. The second simulation scenario examines the hydraulic fracture propagation behavior when a nearby producing well has depleted the surrounding regions. It is important to understand how hydraulic fracture behaves in these complex subsurface conditions that commonly lead to the development of complex hydraulic fracture, which poses significant challenges in optimizing hydraulic fracture design and execution.

Since the impact of formation heterogeneity properties is beyond the scope of this study, the reservoir rocks are assumed to be homogenous and isotropic in the following simulation cases. CZM is implemented to describe the damage initiation and evolution in both brittle and ductile formations, by using predefined cohesive law and energy based fracture propagation criteria. For brittle formations, rock behaves linear elastically until fracture initiated. However, for ductile formations, additional constitutive equation is needed. Even though CZM method can be used to capture the effects of microscopic inelastic 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 Wang et al. [START_REF] Wang | Poroelastic versus Poroplastic Modeling of Hydraulic Fracturing[END_REF] has demonstrated that the impact of plastic deformation on hydraulic fracture propagation cannot be well captured by imposing an artificial 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 plastic deformations are not constrained to the very close vicinity of the fracture tip. In order to account for any potential plastic damage in the whole simulation domain, Mohr-Coulomb criterion is implemented. When the rock in-situ stress state reaches the Mohr-Coulomb failure envelope, plastic deformation occurs in the region. 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 rock plastic deformations is captured by the Mohr-Coulomb flow theory of plasticity within the entire simulation domain.

It is possible that the plastic deformations induced by shear failure may alter the local permeability in the hydraulic fracture stimulated regime, if micro-fractures or weak planes are pre-existing. 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, 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), which is also often referred to as MiniFrac, Mini-Fall-off (MFO) and Diagnostic Fracture Injection Test (DFIT) in literature. Investigate the shear damage induced permeability change inside the associated stimulated reservoir volume (SRV) opens up a complete new issue and hence, will not be discussed further. In the following simulations, the rock matrix permeability is considered as a constant property.

Simultaneous Hydraulic Fracture Propagation

Horizontal wells with multiple fractures are now commonly used in low-permeability, unconventional reservoirs. The spacing between perforations, the number and orientation of transverse fractures all have a major impact on well production. The hydraulic fracture propagation behavior is controlled by the reservoir rock properties, formation in-situ stresses and operation conditions, however, the propagating or propped hydraulic fracture itself can also change the direction and magnitude of insitu stresses within the stimulated area. The stress interference (also known as stress shadow effects) between multiple hydraulic fractures can significantly impact the optimal fracture spacing [START_REF] Roussel | Optimizing Fracture Spacing and Sequencing in Horizontal Well Fracturing[END_REF] and treatment design [START_REF] Wu | Simultaneous Multifracture Treatments: Fully Coupled Fluid Flow and Fracture Mechanics for Horizontal Wells[END_REF]. In this section, the stress shadow effects on the behavior of simultaneous propagating fractures, as well as the impact of plastic damages in ductile reservoir rocks are investigated. Fig. 5. shows the simulation model for simultaneous hydraulic fracture propagation at the beginning of injection. The fracture path at the perforation points are defined as initially open to allow entry of the fluid at the perforation tunnels, so that the initial flow and fracture growth are possible. The whole simulated domain is saturated with reservoir fluid and constant total injection rate and equal pressure conditions are applied at the perforations. All the outer boundaries have zero displacement along the direction that perpendicular to its surface, and constant pore pressure condition is imposed on the outer boundaries.

Fig. 5 Simulation model for simultaneous hydraulic fracture propagation

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]. Detailed discussions and sensitivity analysis on plastic properties on planar hydraulic fracture propagation have been presented in previous study [START_REF] Wang | Poroelastic versus Poroplastic Modeling of Hydraulic Fracturing[END_REF]. Table 1 shows all the input parameters for the presented synthetic simulation in this section. In order to examine the impact of elastic and plastic rock deformations on the behaviors of simultaneous propagating hydraulic fractures, three cases with different rock properties are examined, which are hard rock (E=20 GPa), soft rock (E=10 GPa), ductile rock (E=10 GPa and plastic deformation is allowed when failure criteria is satisfied)

Cohesion strength

First, the results of the presented model is compared with KGD asymptotic solutions [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. Fig. 6 demonstrates that the simulation results match well with analytical solutions when restrictive conditions are applied. It should be mentioned that 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. Dynamic meshing or post process techniques can be used to smooth the results of simulation data.

Fig. 6 Model verification with KGD analytical solutions, E=10 GPa Fig. 7 shows the fracture propagation path and the induced shear stress distribution for soft rock with 10 GPa Young's modulus. It can be noted that the fractures first initiated along the direction of perforation that perpendicular to the initial minimal horizontal stress, and then they gradually change the propagation direction and starts to propagate away from each other, due to local stress changes that result from stress shadow effects. It can be also observed that the shear stress is intensified in front of fracture tip, where multiple shear zones are developed due to local stress disturbance by propagating fractures. Because of the symmetric nature of the simulated domain, the shear stress distribution is also symmetric. It should be mentioned that the displace field presented in all the following simulation results are enlarged by a scale factor of 50, so a better visual observation of the fracture geometry and formation deformations can be obtained. Fig. 8 shows the reservoir pore pressure distribution at two different stages. It can be observed that a high pressure zone can be developed at the tip of both hydraulic fractures. This phenomena can be reasonably explained by the fact that that when a complete new fracture surface is generated within a cell at the fracture tip, there is a sudden fluid leak off into the adjacent formation, which increase the local pore pressure. Because fracturing fluid continues leak-off into the surrounding formations through exposed fracture surface, the overall pore pressure is higher around the propagating hydraulic fractures. It is also interesting to notice that the pore pressure in the area that between the simultaneous propagating hydraulic fractures is relative high, because unlike the leak-off fluid closer to the left and right boundary that can be dissipated into far field, the trapped leak-off fluid in this region can only diffuse through a narrow window that between the two hydraulic fracture tips. So besides the stress interference that transmitted by rock matrix itself, the elevated pore pressure in the region between simultaneous propagating hydraulic fractures can also push the fractures propagating away from each other. Fig. 9 shows the fracture propagation path and the induced shear stress distribution for hard rock with 20 GPa Young's modulus. Compare to Fig. 7, It can be noted that even though the shear stress follows the same pattern, but the fractures are significantly longer at the end of injection (t=1800 seconds). This is due to the fact that the more stiff the formation rock, the thinner the fracture width. And with the same amount of injection volume, longer fracture geometry is expected for hard reservoir rocks.

Fig. 9 Shear stress intensified during hydraulic fracture propagation for E=20 GPa (left side t=180 seconds; right side t=1800 seconds) Fig. 10 shows the plastic strain and plastic deformation area when in-elastic behavior is considered in ductile rocks. It can be observed that the closer to the propagating hydraulic fractures, the more severe the plastic deformation. It can be also noticed that, when compared to Fig. 7, the fractures have shorter and wider geometry. Even though plastic failure happens when local in-situ stress state reaches the Mohr-Coulomb failure envelope, stress shadow effects still push the fractures propagating away from each other. The ductility of reservoir rocks does not prevent the occurrence of stress shadow effects.

Fig. 10 Plastic damage during hydraulic fracture propagation for E=10 GPa when plastic deformation is allowed (left side t=180 seconds; right side t=1800 seconds)

Fig. 11 shows the pressure and fracture width at the wellbore during simulation. Because of the symmetry of this presented model, the pressure and fracture width is represented by only one of the fractures. In addition, in order to compare the behavior of single propagating fracture and simultaneous fractures, the simulation results from a single propagating fracture in the same simulation domain are included as a reference. First, consider the cases of brittle rocks when the formation only behaves linear elastically. It can be observed that fracture propagation pressure is lower and fracture geometry is wider when only one fracture is propagating in the simulation domain. The stress shadow effects that induced by two simultaneous propagating fractures increase the fracture propagation pressure and reduce the fracture width. On top of that, if double the reservoir rock stiffness, fracture with can be further reduced. So in brittle formations, the stress shadow effects are more severe when the reservoir rock's Young's modulus is high. However, If the ductile nature of the formation rock is considered (inelastic deformations can happen due to shear failure), the fracture propagation behavior and fracture geometry can be significantly different from that predicted by models that only capture the effects of poro-elasticity. It can be clearly observed that plastic deformation increases the breakdown and propagation pressure, and a more interesting phenomenon is that: even under stress shadow effects, fracture width of the two simultaneous propagating fractures in ductile rocks is wider than that of a single propagating fracture in brittle rocks. So it is obvious that the occurrence of plastic deformations can offset the effects of stress shadow on the reduction of fracture width if shear induced plastic failure happens.

Fig. 11 Pressure and fracture width at the wellbore for different simulation cases

Producing Well Interference

Not all wells are hydraulic fractured at initial reservoir conditions, and in many cases, hydraulic fracturing can be severely impacted by producing wells, because these producing wells may have depleted the adjacent regimes and altered the in-situ stress and pore pressure conditions. This can also happen when fracturing an infill well where existing nearby horizontal wells have depleted the surrounding reservoir drainage area. It has been reported that poorly hydraulic fracturing job due to producing well interference can lead to significant reduction in stimulation efficiency and gas production [START_REF] Peza | How Fracture Interference Impacts Woodford Shale Gas Production[END_REF]. In this section, a synthetic case of hydraulic fracturing that affected by a nearby producing well is investigated, as shown in Fig. 12. The producing well with a hydraulic fracture of 40 meters length is located that the right corner of the simulation domain (the other hydraulic fracture of the producing well that outside areas outside the simulation domain are not included). All the outer boundaries have zero displacement along the direction that perpendicular to its surface, the same as in previous section. All the input parameters are provided in Table 2. In order to present the effects of pore pressure depletion by the producing well, a constant production pressure of 2 MPa is imposed on the hydraulic fracture of the producing well, while all the other boundaries have symmetric conditions. A pressure transient period is simulated until most of the simulation domain is impacted by the producing well; the result is shown in Fig. 13. It can be observed that the pressure around the hydraulic fracture of the producing well is around 2 MPa while the pressure at the left boundary of simulation domain is 15 MPa, which is the original pore pressure. This pressure distribution is assumed to be the in-situ pressure at the time of hydraulic fracturing execution and will be mapped to the following hydraulic fracture propagation simulations as the initial reservoir pressure condition Fig. 14 shows the fracture propagation path and the induced shear stress distribution. It can be noted that the fractures first initiated along the direction of perforation and propagated along a straight line for around 20 meters, and then it gradually change its propagation direction towards the producing well. It can be also observed that the distribution of shear stress is not symmetric on both sides of the fracture due to unevenly distributed reservoir pore pressure and the turning fracture geometry. Fig. 15 shows the pore pressure distribution after 60 minutes of injection. It can be observed that the pressure is significantly lower around the producing well than the far field region on the left side of the simulation domain. On one hand, the low pressure zone has lower local stress; on the other hand, low pressure zone can attract more fracturing fluid to divert its flow path. These two mechanisms draw the hydraulic fracture to propagate toward the producing well. Unlike the case of simultaneous propagating fracture, where stress shadow effects lead to non-planar hydraulic fracture, the unevenly distributed pore pressure due to the interference of nearby producing well can also lead to the development of non-planar hydraulic fracture.

Fig. 15 Pore pressure distribution within simulation domain for brittle rock, t=3600 seconds Fig. 16 shows the plastic strain and plastic deformation area if formation in-elastic behavior is considered. Even though Fig. 14 indicates the distribution of shear stress is unsymmetrical on both sides of the propagating fracture, but the plastic deformation seems occurs evenly on both sides of the fracture. This is because only the intensified shear stress region near the propagating fracture tip that reaches the failure envelope matters as regard to the shear failure and the magnitude of shear stress is almost the same at both sides near the fracture tip.

Fig. 16 Plastic strain and plastic deformation area within simulation domain for ductile rock, t=3600 seconds Fig. 17 shows the pressure and fracture width at the injection wellbore during simulation. Once again, it demonstrates that plastic deformations induced by shear failure increase fracture propagating pressure and fracture width.

Fig. 17 Pressure and fracture width at the wellbore for brittle and ductile reservoir rocks

Conclusions and Discussions

In this study, a fully coupled, poro-elasto-plastic hydraulic fracture numerical model is presented based on the combination of XFEM based CZM and Mohr-Coulomb theory of plasticity. The physical process involves fully coupling of the fluid leak-off from fracture surface and diffusion into the porous media, the rock deformation and the fracture propagation. XFEM is implemented to determine the arbitrary solution dependent fracture path and CZM is used to model fracture initiation and damage evolution without generating stress singularity at the fracture tip. In addition, Mohr-Coulomb criterion is implemented to account for any potential plastic damage in the whole simulation domain. Two representative synthetic field case scenarios are investigated with the presented model, which provides some fundamental insight into the behavior of complex fracture propagation in both brittle and ductile formations. The results indicate that the stress interference of multiple simultaneous propagating fractures not any increases fracture propagation pressure and reduces fracture width, but also lead to the development of complex fracture geometry. Even though the ductile nature of reservoir rocks does not prevent the occurrence of stress interference by multiple propagating fractures, plastic deformation during fracture propagation can offset the effects of stress shadow on the reduction of fracture width. When hydraulic fracture initiated or propagating in a region that the pore pressure is influenced by nearby producing well, the fracture tend to propagate towards the producing well where the surrounding pore pressure is relative low. In all cases, plastic deformations induced by the propagating fracture lead to higher propagation pressure and shorter, wider fracture geometry. The model presented in this article can be used as a powerful tool for complex hydraulic propagation modeling to investigate fracture propagation behaviors in various reservoir stress and geological conditions, but due to the nature of the involved complexity in hydraulic fracturing modeling and subsurface uncertainties, both laboratory experiment and field data (eg., treatment pressure, fracture diagnostic data, micro seismic imaging and post-fracture production rate etc.) are needed to calibrate this fully coupled model before it is applicable to guide hydraulic fracture design and field operations for a specific reservoir.

Nomenclature

Fig. 1

 1 Fig. 1 Sketch of a plane-strain hydraulic fracture with varying aperture

Fig. 7

 7 Fig. 7 Shear stress intensified during hydraulic fracture propagation for E=10 GPa (left side t=180 seconds; right side t=1800 seconds)

Fig. 8

 8 Fig. 8 Pore pressure distribution during hydraulic fracture propagation for E=10 GPa (left side t=180 seconds; right side t=1800 seconds)

Fig. 12

 12 Fig. 12 Simulation model for producing well interference

Fig. 14

 14 Fig. 14 Shear stress distribution during hydraulic fracture propagating (left side t=180 seconds; right side t=3600 seconds)

𝑐=

  𝑐 𝑙 = Leak-off coefficient, m/Pa•s 𝑑 = Initial thickness of cohesive surfaces, m 𝐸 = Young's modulus, Pa 𝒇 = Body force per unit volume, N/m 3 𝐺 = Potential function 𝐺 I Fracture critical energy in Mode I, J 𝐺 II 𝑐 = Fracture critical energy in Mode II, J 𝐺 III 𝑐 = Fracture critical energy in Mode III, J 𝐈 = Unit Matrix 𝒌 = Formation permeability vector in formation, md 𝐾 𝑛 = Cohesive stiffness, Pa/m 𝑝 𝑓 = Fluid pressure inside fracture, Pa 𝑝 𝑚 = Pore pressure in formation matrix, Pa 𝒒 𝒇 = Flow rate along fracture per unit height, 𝑚 3 /𝑠 𝑞 𝑙 = Leakoff rate per unit fracture surface area, 𝑚 3 /(𝑠 • 𝑚 2 ) 𝒒 𝒎 = Fluid flux velocity in formation matrix, m/s 𝒕 = Surface traction per unit area, pa/m 2 𝑇 = Traction force, Pa damage initiation, m 𝛿 𝑓 = Displacement at failure, m 𝛿 𝑚 = Effective displacement, m 𝜀 = Normal strain 𝐾 𝐼𝐶 = Rock fracture toughness in Mode I, Pa√m 𝜇 = Fluid viscosity, cp 𝜈 = Poisson's ratio 𝑣 𝐿 = Leakoff rate, m/s 𝜌 𝑓 = Fluid mass density, kg/m 3 σ = Stress, Pa 𝜎 ′ = Effective stress, Pa 𝜏 𝑓 = Shear strength, Pa ∅ =Porosity of porous media

Table 1

 1 Input Parameters for model simultaneous hydraulic fracture propagation

	Input Parameters	Value
	Young's modulus	10 GPa , 20 GPa
	Poisson's ratio	0.25
	Fluid viscosity	1 cp
	Tensile strength	1 MPa
	Formation permeability	1 md
	Injection rate per thickness	0.0005 (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
	Leak-off coefficient	1E-14 m 3 /s/Pa
	Porosity	0.1
	Friction angle	27 o
	Dilation angle	8 o

Table 2

 2 Input Parameters for model of producing well interference

	Input Parameters	Value
	Young's modulus	15 GPa
	Poisson's ratio	0.25
	Fluid viscosity	1 cp
	Tensile strength	1 MPa
	Formation permeability	10 md
	Injection rate per thickness	0.0005 (m 3 /s/m)
	Specific weight of fluid	9.8 kN/m 3
	Initial pore pressure	15 MPa
	Producing well pressure	2MPa
	Maximum horizontal stress	30 MPa
	Minimum horizontal stress	25 MPa
	Vertical stress	35 MPa
	Critical Fracture Energy	100 J/m 2
	Leak-off coefficient	1E-13 m 3 /s/Pa
	Porosity	0.1
	Friction angle	27 o
	Dilation angle	8 o
	Cohesion strength	3 MPa