The Hanyi Wang

Matteo Marongiu-Porcu

Michael J Schlumberger

Economides

Poroelastic and Poroplastic Modeling of Hydraulic Fracturing in Brittle and Ductile Formations*

The prevailing approach for hydraulic fracture modeling relies on Linear Elastic Fracture Mechanics (LEFM). Generally, LEFM that uses stress intensity factor at the fracture tip, gives reasonable predictions for hard rock hydraulic fracturing processes, but often fails to give accurate predictions of fracture geometry and propagation pressure in formations that can undergo plastic failures, such as poorly consolidated/unconsolidated sands and ductile shales. The reasons are that the fracture process zone ahead of the crack tip, elasto-plastic material behavior and strong coupling between flow and stress cannot be neglected in these formations. Recent laboratory testing has revealed that in many cases fracture propagation conditions cannot be described by traditional LEFM models. Rather, fractures develop in cohesive zones. In this study, we developed a fully coupled poro-elasto-plastic hydraulic fracturing model using cohesive zone method, which is not only able to model fracture initiation and growth considering process zone effects, but can also capture the effects of plastic deformation in the bulk formation. The impact of the formation plastic properties on the fracture process is investigated, and the results are compared with existing models. In addition, the effects of different parameters on fracture propagation in ductile formations are also investigated through parametric study. The results indicate that plastic and highly deforming formations exhibit higher breakdown and propagation pressure. The more plastic the formation (lower cohesion strength), the higher is the net pressure required to propagate the fracture due to increased plastic deformation. Also, lower cohesion strength leads to shorter and wider fracture geometry. The effect of formation plasticity on a hydraulic fracture is mostly controlled by initial stress contrast, cohesion strength of formation rock and pore pressure. We also found that altering the effective fracture toughness can only partially mimic the consequences of increased toughness ahead of the fracture tip in ductile formations, but it fails to capture the effect of plastic deformation within the whole affected area, which can lead to underestimating the fracture width and overestimating the fracture length. For a more accurate modeling of fracturing in ductile formations, the whole plastic deformation region induced by the propagating fracture should be considered, especially when shear failure areas are large.

Introduction

Hydraulic fracturing has been widely used as a common practice to enhance the recovery of hydrocarbons from low to very low permeability reservoirs, as well as prevent sand production in high permeability reservoirs [START_REF] Economides | Reservoir Stimulation[END_REF]. Development of theoretical models of hydraulic fracturing started several decades ago, and provided a basis for hydraulic fracturing design, optimization and diagnostics. As a stimulation tool, 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. The mathematical formulation of the problem is represented by a set of nonlinear integro-differential equations. Also, the problem has a moving boundary where the governing equations degenerate and become singular. The complexity of the problem often restricts researchers to consider only simple fracture geometries. An assumption commonly made is that the rock is linearly elastic and the loss of fluid into the rock is evaluated using a one-dimensional diffusion equation. The mechanics of fracture opening and fluid loss are considered independent and their interaction is not considered, such as the KGD model [START_REF] Geertsma | A rapid method of predicting width and extent of hydraulic induced fractures[END_REF][START_REF] Khristianovich | Formation of vertical fractures by means of highly viscous fluids[END_REF] and PKN model [START_REF] Nordgren | Propagation of vertical hydraulic fractures[END_REF][START_REF] Perkins | Widths of hydraulic fractures[END_REF]. In order to include the effect of fluid diffusion in porous media, many numerical models have been proposed based on the theory of poroelasticity [START_REF] Adachi | Plain strain propagation of a hydraulic fracture in permeable rock[END_REF][START_REF] Carrier | Finite Element Modeling of Fluid-Driven Fracture in Permeable Medium[END_REF][START_REF] Detournay | The near-tip region of a fluid-driven fracture propagating in a permeable elastic solid[END_REF]. Although classic hydraulic fracturing simulators based on LEFM are convenient to use and can give reasonable predictions for brittle formations, they often fail to predict fracturing pressure and geometry with enough certainty in formations that can undergo plastic failures, such as poorly consolidated/unconsolidated sands and ductile shales.

By the proposition of [START_REF] Orowan | Energy criteria of fracture[END_REF] and [START_REF] Irwin | Analysis of stresses and strain near the end of crack traversing a plate[END_REF], the medium around a crack tip can be divided into two regions, defined as "Process Zone" and "Universal Singularity", as illustrated in Fig. 1. All microscopic details (complex processes of bond breaking and non-elastic behavior in front of the crack tip) are lumped together and hidden in the Process Zone. The assumptions behind the LEFM are that the material is linear elastic in the entire region right up to the point of crack tip, and the Process Zone is negligibly small compared to the macroscopic dimensions of the crack, so the stress and strain fields in the region of Universal Singularity (stress surrounding the crack tip develops a square root singularity) can have a singular form which is given by the so-called Stress Intensity Factors.

Fig. 1-Regions around the tip of a propagating crack

Numerous laboratory studies have shown that deformation of weakly cemented sands and ductile 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]. Ductile rocks often exhibit soft behavior (low Young's modulus) and according to the theory of LEFM, the lower the Yong's modulus, the lower the propagating pressure. However, on the contrary, many field study and surveys found that the net pressures observed in ductile formations are much higher than that predicted by LEFM, especially in unconsolidated sands, and these observations have triggered a series of dedicated studies which looked into the importance of the plastic deformation in hydraulic fracturing [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 influence of plasticity in hydraulic fracturing[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]. The disparities between field observations and the results predicted by LEFM indicate that the complex phenomena in the process zone beyond the fracture tip, elasto-plastic rock behavior and strong coupling between flow and stress cannot be neglected in these ductile formations.

Besides the commonly known poorly consolidated/unconsolidated sands, the fact that shale rocks can also exhibit in-elastic behavior has been largely ignored in both industry and academia. Not all shales are as brittle as the Barnett shale; depending on mineral compositions (such as clay, carbonates, quartz, etc.), many shales can exhibit a pronounced ductile behavior, such as the Eagle Ford and the Haynesville Shale. The ductile nature of these shale plays poses a 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], and the production from hybrid fracturing treatments (creating highly conductive primary fractures) wells outperforms the water-frac (creating large fracture surface area) wells in the Eagle Ford shale [START_REF] Stegent | Engineering a Successful Fracture-Stimulation Treatment in the Eagle Ford Shale[END_REF]; this demonstrates that creating a high-conductivity fracture may be more critical to long-term production than creating large surface area in ductile shales. Such trends can be explained by the fact that in clay-rich ductile shale, the level of complexity of the induced hydraulic fracture geometry is relatively low, dominant "bi-wing" fractures are more likely to occur and unpropped fractures will most likely not be conductive after initial production. Thus, being able to predict fracture geometry and fracturing pressure with enough confidence is not only important in poorly consolidated/unconsolidated sands, but also crucial in ductile shales, because of the long term production hinges on the portion of the primary fracture that can be propped with enough size and concentration of proppants. A completely different scenario (beyond the domain of interest of this study) opens up when a hydraulic fracture propagates in naturally fractured brittle formations, where the interaction with pre-existing natural fractures becomes the dominant phenomenon [START_REF] Weng | Modeling of Complex Hydraulic Fractures in Naturally Fractured Formation[END_REF].

In examining the plentiful literature concerning fracture mechanics, the cohesive zone modeling has been adopted by many authors to model fracture initiation and propagation by considering the effects of in-elastic behavior in front the fracture tip. The cohesive zone is a region ahead of the crack tip that is characterized by micro-cracking along the crack path. The main fracture is formed by interconnection of these micro cracks due to damage evolution. The physical meaningless tip singularity predicted by LEFM can be resolved with the idea of cohesive zone model [START_REF] Mokryakov | Analytical solution for propagation of hydraulic fracture with Barenblatt's cohesive tip zone[END_REF]. The conception of 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. Mathematically, the [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] model is similar to the Barenblatt model with the exception that different conditions (the condition of constant cohesive zone length and the condition of maximal cohesive zone opening) were applied to the cohesive zone parameters. Because cohesive properties in the cohesive zone are material parameters that do not depend on fracture dimensions and loading. The classical Griffith's brittle elastic fracture model ( 1921) is a limited case of Barenblatt's 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.

In recent works of adapting cohesive zone method for modeling of hydraulic fracturing, [START_REF] Dean | Hydraulic fracture predictions with a fully coupled geomechanical reservoir simulator[END_REF] proposed a fully coupled geo-mechanical reservoir simulator that incorporates two different fracture propagation criteria: a conventional LEFM criterion based on stress intensity factor and a cohesive strength criterion based on energy release rate, within the cohesive zone in front of propagating hydraulic fracture tip. [START_REF] Chen | Cohesive zone finite element based modeling of hydraulic fractures[END_REF] established a model based on the existing pore pressure cohesive finite elements to investigate the propagation of a penny-shaped hydraulic fracture in an infinite elastic medium. The finite element results match well with the analytical solutions where the fracture process is dominated by rock fracture toughness. [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 zone based on the KGD model by assuming impermeable elastic rock without considering leak-off. It shows that the derived solutions from cohesive zone model can fit the pressure log much more accurately than LEFM for the case of fracturing ductile formations. [START_REF] Sarris | Modeling of Hydraulic Fracturing in a Poroelastic Cohesive Formation[END_REF] developed a poroelastic cohesive zone model and investigated the effects of permeability, injection rate and formation compressibility on fracture geometry, it is concluded that higher pressures are needed to extend a fracture in a poroelastic medium than in an elastic medium. [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 include the process zone effect on fracturing in ductile rock. The results show that with the increase of effective fracture toughness, the cohesive zone model gives predictions that are more conservative on fracture length as compared with pseudo 3D and PKN models.

Even though all these models are able to account for in-elastic behavior beyond the fracture tip by lumping all the effects into the cohesive zone, they do not consider plastic deformation in the formation bulk. In reality, depending on the in-situ stresses and rock properties, the plastic deformation can occur at a certain distance from the fracture tip, without limiting itself to the very close vicinity of the fracture tip. In this study, we developed a hydraulic fracture model for both poroelastic and poroplastic formations based on the cohesive zone method, where the effects of plastic deformation in the whole simulated region are included. The physical process involves full coupling of the fluid leak-off from fracture surface and diffusion into the porous media, the rock deformation and the fracture propagation. The model is then solved by finite element method. We investigated how the plastic failures and fluid/rock interactions can impact the fracture geometry, breakdown and propagation pressure in ductile formations and the results are compared with that in brittle formations. Additional content and further discussions on fracture reorientation, multiple fracturing and producing well interference in both brittle and ductile reservoir rocks are presented in companionate articles (Wang 2015a;2015b)

Mathematical Model

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.

General Equations

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. Most hydraulic fracturing fluids for high permeability formations exhibit power law rheological behavior and temperature-related properties, while for low permeability unconventional formations, water is commonly used as hydraulic fracturing fluid. In order to avoid additional complexity added by fluid behavior, incompressible and Newtonian fluid is assumed in this study. The continuity equation which imposes the conservation of mass in one dimensional flow is [START_REF] Boone | A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media[END_REF] where 𝑣 𝐿 is leak-off velocity, 𝐶 𝐿 is leak-off coefficient, 𝑡 is the time elapsed since the start of the leakoff process, 𝑉 𝐿 is the cumulative fluid volume that passes though the surface 𝐴 𝐿 during the time period from time zero to time 𝑡, the integration constant 𝑆 𝑝 is called spurt loss, which represents a sudden fluid loss into the formation when new fracture surface area are created and a stable leakoff status has not been established yet within the newly create fracture surface area. The two coefficients 𝐶 𝐿 and 𝑆 𝑝 , can be determined from laboratory tests or, preferably, from the evaluation of Fracture Calibration Tests (FCT), which is also often referred to as MiniFrac, Mini-Fall-off (MFO) and Diagnostic Fracture Injection Test (DFIT) in literature. [START_REF] Notle | A General Analysis of Fracturing Pressure Decline With Application to Three Models[END_REF], [START_REF] Valkó | Fluid-Leakoff Delineation in High-Permeability Fracturing[END_REF] and [START_REF] Marongiu-Porcu | Comprehensive Fracture Calibration Test Design[END_REF] presented vast theoretical modeling and field example analyses on this specific subject. For high permeability formations, spurt loss can contribute to the majority of the fluid leakoff volume even with the addition of leakoff prevention additives. For low permeability formations, spurt loss is negligible unless the fissure system is considered, if present. Setting 𝐴 𝐿 to be the unit surface area, and discretizing Eq. 3 with respect to time, the difference of cumulative leakoff volume 𝑉 𝐿 between each time intervals can be calculated. When the discretized time intervals are small enough, the changes of cumulative leakoff volume 𝑉 𝐿 divided by time interval are approximate to the leakoff rate per unit surface area , 𝑞 𝑙, within that specific time interval. For a flow between parallel plates, the conservation equation of momentum balance can be used to relate 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]: where ∇𝑝 is the formation pressure gradient vector, 𝒒 𝒎 is the fluid flux velocity vector in the porous media, 𝒌 is formation permeability tensor. When the porous medium considered is isotropic, the permeability tensor will appear in its diagonal form, with equal valued elements on the diagonal. Biot's poroelastic theory ( 1941) is used to couple rock deformation and pore pressure. Together with fluid mass conservation, the whole system of fluid-rock deformation can be solved with appropriate initial and boundary conditions.

𝑞 𝑓 = 𝑢𝑤 =

Fracture Initiation and Propagation

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 is negligible 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 shear failure, which put the adequacy of such assumptions 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 a critical value. The cohesive method proposed 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] 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 [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] brittle elastic fracture model is a limited case of Barenblatt's 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. In essence, cohesive zone is way of representing the process in front of fracture tip mathematically.

The constitutive behavior of the cohesive zone is defined by the traction-separation relation, which includes initial loading, initiation of damage, and the propagation of damage that leading to eventual failure at the bonded surface. The behavior of the interface prior to initiation of damage is often described as linear elastic in terms of a penalty stiffness that degrades under tensile and/or shear loading but is unaffected 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. A cohesive potential function G is defined so that the traction is given by where T is traction force and 𝛿 is the displacement across a pair of cohesive surfaces. For Mode-I fracture (induced by forces that are perpendicular to the crack), the traction force T can be interpreted as a tensile stress in front of fracture tip. Based on different types of materials, different forms of traction-separation relations in cohesive zones can be adopted to simulate the dynamic fracture process. [START_REF] Tomar | Bounds for element size in a variable stiffness cohesive finite element model[END_REF] proposed a bilinear cohesive law, as shown in Fig. 2. It assumes that the material exhibit linear elastic behavior before the traction reaches the tensile strength T Max or the separation of cohesive surfaces exceeds the displacement of damage initiation, 𝛿 0 . Beyond 𝛿 0 , the traction reduces linearly to zero up to the displacement of complete failure 𝛿 𝑓 (linear softening), and any unloading takes place irreversibly. 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 IC by [START_REF] Kanninen | Advanced fracture mechanics[END_REF]: 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 may play an important role in ductile formations under certain loading conditions. The traction-separation law presented above for Mode-I based fracture criteria can be easily extended to mixed-mode fracture criteria. For an isotropic formation, the traction-separation responses in different modes are assumed to be the same in this study, as shown in Fig. 3. where t n , t s , t t refer to the nomal, the first, and the second shear stress components on the cohesive interface; 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 correspond to the displacements of initial damage in the normal, the first, and the second shear stress direction and δ n f , δ s f , δ t f are the displacements of complete failure in these three directions. 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 ( 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 in this model 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

T = 𝜕G 𝜕𝛿
G I c = K IC 2 E (1 -𝜈 2 )
G c = G I c + (G II c -G I c ) ( G shear G

Simulation Model

When fluid is pumped into the reservoir at a sufficiently high pressure, the reservoir rock will break down and open. Continued pumping of the fluid at this pressure will propagate the induced crack into the reservoir. In this model, a pre-defined surface made up of elements that support the cohesive zone traction-opening calculation is embedded in the rock and the hydraulic fracture grows along this surface. This typical cohesive zone model is illustrated in Fig. 4. The fracture process zone (unbroken cohesive zone) is defined within the separating surfaces where the surface tractions are nonzero (Traction-Separation law is still valid). The fracture is fully filled with fluid in the broken cohesive zone where no traction from rock fracture exists but fluid pressure is acting on the open fracture surfaces. When the energy release rate reaches the critical value G c in front of the propagating fracture tip (unbroken cohesive zone), the associated cohesive element will be completely damaged and fluid filled fracture will propagate. 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]. 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]. [START_REF] Nassir | Geomechanical Coupled Modeling of Shear Fracturing in Non-Conventional Reservoir[END_REF] The numerical solution procedure is schematically illustrated in Fig. 6. The coupled system of non-linear equations is solved numerically by Newton-Raphson method, which requires evaluating Jacobian matrix, start from initial guess values to approximate to the solution through successive iterations. 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 finite element package ABAQUS. In the following part, we'll compare the simulation results from our poro-elasto-plastic model with other existing models from literature to reveal how the fracture geometry, the net pressure and the fluid efficiency can be impacted by plastic deformations due to shear failure, then we carry out a set of parametric studies to determine which parameters mostly control the influence of plasticity on fracture propagation in ductile formations.

Results and Analysis

In this section we present our simulation results to demonstrate the importance of plasticity in modeling hydraulic fracturing in ductile formations. Our simulation results are compared with two other cohesive zone based hydraulic fracture models, which are the poroelastic model [START_REF] Sarris | Modeling of Hydraulic Fracturing in a Poroelastic Cohesive Formation[END_REF]) and the poroelastic model with effective toughness [START_REF] Yao | Linear Elastic and Cohesive Fracture Analysis to Model Hydraulic Fracture in Brittle and Ductile Rocks[END_REF].

The idea of using effective toughness is to modify the poroelastic model, by lumping all plastic deformation effects into an increased toughness beyond the fracture tip [START_REF] Papanastasiou | The effective fracture toughness in hydraulic fracturing[END_REF]. In principle, the effective toughness can be determined either by investigating experimentally the formation rock softening behavior, or by matching the predicted net pressure to the actual net pressure data from a FCT. In this study, the effective toughness is determined by matching the poroelastic model calculated net pressure to the one predicted by our fully coupled poro-elasto-plastic model. The final adjusted effective toughness for the poroelastic model is 3 times larger (3 MPa√𝑚) than the original toughness (1 MPa√𝑚) declared in Table 1. Fig. 7 shows the net pressure predicted by three different models investigated in this study. The predicted net pressure from poroelastic model with increased effective toughness matches the net pressure from poro-elasto-plastic model significantly well, while the predicted net pressure from poroelastic model keeps itself at much lower values throughout the entire simulated injection. This means that higher net pressure levels are actually needed to initiate and propagate hydraulic fractures in plastic formations. It should be mentioned that not only ductile formation itself can lead to difficulties in initiating hydraulic fracture [START_REF] Parker | Haynesville Shale: Hydraulic Fracture Stimulation Approach[END_REF]), field study [START_REF] Yuan | High Pressure Stimulation -Impact of Hydraulic Fracture Geometry to Unconventional Gas Appraisal and Development in Compressional Settings[END_REF]) also indicates that perforations-associated damaged rock can behave like a plastic rubber buffering between the reservoir and fracture entry, which can lead to significant increase of formation breakdown pressure.

Fig. 7-Net pressure prediction with different models

However, as demonstrated in Fig. 8, the predicted fracture is wider in the poro-elasto-plastic model than in the poroelastic model with increased effective toughness, even though both models have the same matched net pressure. This is due to the fact that the effective toughness method considers uniquely the increased toughness within the very close vicinity of fracture tip due to rock plasticity, but it fails to capture the effects of plastic deformation that result from shear failure in the bulk formation induced by the instantaneous fracture propagation. So using CZM alone by increasing fracture toughness (i.e., increase critical fracture energy in the cohesive zone) is inadequate in modeling fracture propagation in ductile formations when the plastic damage is not constrained to the very close vicinity of fracture tip. Nevertheless, both models predict wider fracture widths than the poroelastic model, which does not consider formation plasticity at all.

Fig.8-Maximum fracture width prediction with different models

Table 2 shows the simulation results of 30 minutes of injection for the three different models, and we can observe that formation plasticity has a substantial impact on the fracture geometry after 30 minutes of injection. The predicted fracture length is reduced by nearly 38 % compared to the one predicted by the poroelastic model when plastic deformations are incorporated. In addition, the estimated fracture fluid efficiency is larger when plastic effects are considered because shorter fracture lengths lead to smaller leak-off surfaces. We can also observe that the use of the effective toughness can only partially mimic rock plasticity, even though it still provides wider and shorter fractures, but it fails to capture the effects of plastic deformations within the whole affected area. For a more accurate modeling of fracture propagation in ductile formations, the effects of the whole plastic deformation region should be included, especially when the shear failure area are large. Whenever this is not done, the fracture width is underestimated and the fracture length is overestimated, although good net pressure matches can be generated

Parametric Study

In this section we study the effects of different parameters on fracture propagation in ductile formations. The parameters include injection rate, leak-off coefficient, pore pressure, horizontal stress contrast and cohesion strength. In the parametric study only one parameter changes while all other input parameters are kept consistent with the base case.

Injection Rate

The injection rate is increased to 5 times the base case. Fig. 9 shows the impact of such injection rate on fracture geometry. It can be observed that wider fracture profiles can be obtained with faster propagation speed for larger injection rates, which in turn results in larger net pressures and wider plastic deformation. Because the fracture width does not increase significantly at a higher injection rate, the fracture propagation speed can be considered as approximately proportional to injection rate. At the end of the 2-minute injection, the fracture half-length with Qx5 injection rate is about 5 times longer than the base case.

Fig. 9-Impact of injection rate on fracture geometry

Leak-off Coefficient

Now we increase the leak-off coefficient 100 times and compare the results with the base case. From Fig. 10, we see that a larger leak-off rate leads to much shorter fracture length, because more fracturing fluid filtrates into the formation, but this appears to have a limited impact on the net pressure and the fracture width at the wellbore, as shown in Fig. 11. Slightly higher net pressure and smaller fracture width are obtained with higher leak-off rate due to poroelastic backstress (i.e., the induced additional confining stress around the fracture surface because of the increased local pore pressure). Even though the leak-off rate has been increased by 100 times, the reservoir permeability in this example is still large enough to allow the pore fluid to move freely and dissipate energy within the same timeframe as the fracture opening time, thus, the changes in pore pressure around the fracture and the poroelastic backstress induced by the extra leak off is relatively small. This shows that, ultimately, the leak-off coefficient has a larger impact on the fracturing fluid material balance (how much injection fluid is contributing to create fracture volume versus how much is leaking into the formation) than on the constitutive response of rock deformation, as long as the permeability is large enough that can prevent the buildup of local pore pressure.

Pore Pressure

The formation pore pressure is increased form 10 MPa (base case) to 20 MPa, while maintaining the same far field stresses. Fig. 12 shows the comparison between simulated net pressures with different pore pressures. The results clearly indicate that the lower the pore pressure, the lower the net pressure is needed to initiate and propagate the fracture. In turn Fig. 13 shows that with higher formation pore pressure, the shear failure induced plastic strain increases (PEEQ denotes equivalent plastic strain), and the plastic deformation area becomes larger, which leads to a shorter and wider fracture geometry. The results also indicate that even a small amount plastic deformation (with PEEQ on the order of 10 -3 ) can make a significant difference in fracturing pressure and propagation geometry.

Horizontal Stress Contrast

Trivially, the closer the initial reservoir stress conditions to the Mohr-Coulomb shear failure surface, more likely shear failure events will be triggered around the fracture during propagation. Here we investigate how the horizontal stress contrast impacts the plastic deformation by increasing such contrast from 1.1 (base case) to 2, and the results are shown in Fig. 14. Compared to the simulation results from base case as shown in Fig. 13 (lower illustration), it can be observed that a higher contrast between the initial minimum and maximum horizontal stresses leads to a shorter and wider fracture geometry with larger plastic deformation area. When the rock in-situ stress is far away from the shear failure envelop, additional stress induced by a hydraulic fracture may still not be large enough to cause shear failure at a certain distance from the fracture. However, when the rock is close to the shear failure envelop under initial condition, stress disturbance induced by the propagating fracture is able to push the rock into the failure envelop more easily, and a larger area around the propagating fracture can be impacted.

Fig. 14-Impact of higher horizontal stress contrast on plastic strain and fracture geometry

Cohesion Strength

The cohesion strength reflects how well the rock grains are cemented, as well as rock resistance to shear failure. Reservoir porous rocks under differential stress can exhibit significant permanent deformation due to intergranular slip and rotation that overcomes the elastic deformation of the grains and grain contacts if the cohesion strength is low. In order to examine the impact of rock cohesion strength on hydraulic fracture propagation, we set the cohesion strength to be 2 MPa, 4 MPa and 6 MPa, to represent formations with relative low (unconsolidated sands or shale with weak planes), moderate (ductile shale with micro-fractures) and high cohesion strength (soft sandstone) in this study.

Fig. 15 and Fig. 16 show the impact of the cohesion strength on the net pressure and the fracture width. It can be found that the more plastic (lower cohesion strength) the formation is, the higher the net pressure is needed to initiate and propagate the fracture, furthermore, lower cohesion strength also leads to wider fracture geometry.

Conclusions and Discussions

In this study, we first presented a brief literature review on hydraulic fracturing modeling based on the theory of Linear Elastic Fracture Mechanics and Cohesive Zone Method, and discussed the associated challenges when dealing with ductile rocks that can undergo plastic deformations during hydraulic fracture propagation. Then, a novel poro-elasto-plastic model that combines Cohesive Zone Method and Mohr-Coulomb flow theory of plasticity was presented; this model is capable of modeling hydraulic fracture propagation in both brittle and ductile formations. More importantly, the impact of formation plastic properties on the fracturing process is deeply investigated and the results are compared with established poroelastic and modified poroelastic models. From the results of this research, we can conclude that:

1. The effects of formation plasticity on a hydraulic fracture are mostly controlled by in-situ stress contrast, cohesion strength of formation rock and pore pressure. 2. Ignoring plastic deformation when modeling hydraulic fracture propagation in formations that can undergo plastic deformations, such as poorly consolidated/unconsolidated sands or ductile shales, can lead to inaccurate prediction of fracture geometry, net pressure and fluid efficiency. 3. The method of using effective toughness (e.g., increases critical fracture energy in the cohesive zone) to account for formation plasticity can only partly mimic the consequences of increased toughness within the very close vicinity ahead of fracture tip, but it fails capturing the effect of plastic deformation that induced by shear failure within the whole affected area, so it may underestimate the fracture width and overestimate fracture length. 4. In comparison to brittle rock, ductile and highly deforming formations exhibit higher breakdown and propagation pressure. The more plastic the formation (lower cohesion strength), the higher is the net pressure required to initiate and propagate the fracture due to increased energy that absorbed by plastic deformations. Also lower cohesion strength lead to shorter and wider fracture geometry. 5. For effective modeling of hydraulic fracturing propagation in ductile formations, the effects of plastic deformation in the whole stimulated reservoir volume should be included. So using Cohesive Zone Method alone, which only accounts for in-elastic behavior at the fracture tip, is not adequate for modeling the overall effects of plastic damage within the entire affected region.

Even though the plastic behavior of poorly consolidated/unconsolidated sands and its impact on hydraulic fracturing has already attracted much attention in both academia and industry, the plasticity of clay rich, ductile shales has been largely overlooked in hydraulic fracturing design and execution, despite many laboratory and field studies have shown that shale rocks can exhibit plastic or even time dependent elastic-visco-plastic behavior. From the results of our study, it seems that the impact of formation plastic deformations on hydraulic fracturing is not trivial. If these ductile formations are treated as pure elastic, then our optimized hydraulic fracturing design may render suboptimal well performance, because of the unsound fracture geometry prediction and the corresponding unappropriated choice of proppant size and concentration. In some cases, the unexpected higher breakdown pressure needed in job execution may not be achievable at all, due to surface equipment and operational limitations.

For heterogeneous formations, the rock properties can transit between brittle and ductile along the section of a horizontal well and the formation themselves can as well be laminated with both brittle and ductile rock layers. A better understanding of the properties of ductile rocks and their impact on hydraulic fracture propagation and well performance is crucial to optimize hydraulic fracturing design in these formations. How these ductile rock intervals impact hydraulic fracture interactions (stress shadow effects), fracture height containment and interactions with natural fractures requires future investigations.

The poro-elasto-plastic hydraulic fracturing propagation model proposed in this study can be used as a powerful tool to address these issues and the results from our simulation study agree well with some field observations. Like for any other model, field data that combines reservoir characterization, fracture diagnostics and post-production analysis will be needed to calibrate this predictive model for a specific reservoir.

Nomenclature

  = 0, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (1) where 𝑞 𝑓 is the local flow rate per unit height along the fracture in direction 𝑥 , 𝑞 𝑙 is the local fluid loss in rock formation per unit fracture surface area and 𝑤 is the crack opening. Eq. 1, the local flow rate 𝑞 𝑓 can be determined by taking fluid leak-off into consideration. Carter's Leak-off model is used in this study to describe the leak-off phenomena (Howard and Fast 1970): . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (2) 𝑉 𝐿 𝐴 𝐿 = 2𝐶 𝐿 √𝑡 + 𝑆 𝑝 , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (3)

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (4) where 𝑝 is the fluid pressure inside the fracture and 𝑢 is the average fluid velocity over the cross-section of fracture. Pressure drop along the fracture can be determined by Eq. 4 with local flow rate and local fracture width. Darcy's Law is used to describe fluid diffusion in the isotropic, porous media: 𝒒 𝒎 = -𝒌 𝜇 ∇𝑝, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (5)

  , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (6)

Fig. 2 -

 2 Fig. 2-Traction-Separation law and fluid inside cohesive surface where 𝐸 is Young's modulus of formation and 𝜈 is Poisson's ratio. Alternative, the critical fracture energy can be determined by laboratory experiment. For bilinear traction-separation law, displacement of complete failure 𝛿 𝑓 can be determined by: 𝛿 𝑓 = 2G I c T Max = 2K IC 2 (1 -𝜈 2 ) ET Max . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (8) The stress-displacement relation in the initial part of the loading curve corresponds to linear elastic deformation as follows: σ = δT Max 𝛿 0 = 𝛿K n , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (9) where K 𝑛 is the stiffness of the stress-displacement relation in the loading regime and is assumed to be equivalent to the stiffness of formation according to the following equation: σ = 𝛿K 𝑛 = Eδ 𝑑 , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (10) where d is the initial thickness of cohesive surfaces and the post-peak linear softening regime is given by the following: σ = T Max [1 -(𝛿 -𝛿 0 ) (𝛿 𝑓 -𝛿 0 ) ] , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (11)

Fig. 3 -

 3 Fig. 3-Traction-Separation law for mix-modesDamage is assumed to initiate when one of the stress components reaches the value of maximum strength of rock material, which can be represented by a quadratic law{ 〈𝑡 𝑛 〉 𝑡 𝑛 0 } 2 + { 𝑡 𝑠 𝑡 𝑠 0 } 2 + { 𝑡 𝑡 𝑡 𝑡 0 } 2

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (16) 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 isotropic failure assumed in this study where G II c = G I c , the cohesive response is insensitive to parameter η, and the fracture will propagate when the energy release rate reaches G I c .

Fig. 4 -Fig. 5 -

 45 Fig. 4-Cohesive zone embedded along fracture path (Modified from Chen et al., 2009)In this study, a two dimensional plane strain fluid-driven fracture in porous rock is simulated and a close view of this model is shown in Fig.5. (The actual simulation domain is much larger to reduce the influence of boundary effects). The initial separation along the fracture path is set zero that represents intact formation and the cohesive elements at the perforation points are defined as initially open to allow entry of the fluid at the perforation tunnel(s), so that the initial flow and fracture growth are possible. The whole area is saturated with reservoir fluid and incompressible Newtonian fluid is injected at a constant rate at the perforation. Grid size along the fracture path is 0.2 m and all the displacements of outer boundaries are fixed. It should be emphasized that the CZM used in this article is implemented by an artificially pre-defined path of cohesive element that embedded along element edges. The cohesive element glues 2 conjugated cohesive surfaces of the reservoir pressuredisplacement element. Under different local stress conditions, the behavior of these cohesive surfaces are governed by the traction-separation laws in normal or shear directions, as depicted by Fig2 (left side) and Fig3. When the cohesive element is completely damaged, the bounded two cohesive surfaces detached and hydraulic fracture propagates.

Fig. 6 -

 6 Fig. 6-Solution flow diagram for fracture propagation procedure.

FigFig. 11 -

 11 Fig. 10-Impact of leak-off coefficient on fracture length

Fig. 12 -

 12 Fig. 12-Impact of pore pressure on net pressure

Fig. 15 -

 15 Fig. 15-Impact of cohesion strength on net pressure

Fig. 17 -

 17 Fig. 17-Plastic strain and fracture geometry with 6MPa and 4MPa cohesion strength

𝑐=

  𝐴 𝐿 = Leakoff surface, 𝐿 2 , 𝑚 2 𝑐 = Cohesion strength, m/L𝑡 2 , Pa 𝐶 𝐿 = Carter's leakoff coefficient, 𝐿/√𝑡, m/√𝑠 𝑑 = Initial thickness of cohesive surfaces, L, m 𝐸 = Young's modulus, m/Lt 2 , Pa 𝐺 = Potential function 𝐺 I Fracture critical energy in Mode I, mL 2 /t 2 , J 𝐺 II 𝑐 = Fracture critical energy in Mode II, mL 2 /t 2 , J 𝐺 III 𝑐 = Fracture critical energy in Mode III, mL 2 /t 2 , J 𝒌 = Formation permeability tensor in formation, L 2 , md 𝐾 𝑛 = Cohesive stiffness, m/L 2 t 2 , Pa/m 𝑝 = Pore pressure, m/L𝑡 2 , Pa 𝑞 𝑓 = Flow rate along fracture per unit height, 𝐿 2 /𝑡, 𝑚 3 /(𝑠 • 𝑚) 𝑞 𝑙 = Leakoff rate per unit fracture surface area, 𝐿 2 /𝑡, 𝑚 3 /(𝑠 • 𝑚 2 ) 𝒒 𝒎 = Fluid flux velocity in formation matrix, L/t, m/s 𝑆 𝑝 = Spurt loss, L, m 𝒕 = Surface traction on cohesive interface, m/Lt 2 , Pa 𝑇 = Traction force, m/L𝑡 2 , Pa 𝑇 𝑀𝑎𝑥 = Cohesive tensile strength, m/Lt 2 , Pa 𝑉 𝐿 = Cumulative leakoff volume, 𝐿 3 toughness in Mode I, m/√Lt 2 , Pa√m 𝜇 = Fluid viscosity, m/Lt, cp 𝜈 = Poisson's ratio 𝑣 𝐿 = Leakoff rate, L/t, m/s 𝑉 𝐿 = Cumulative leak-off volume, 𝐿 3 , 𝑚 3 σ = Stress, m/Lt 2 , Pa 𝜏 𝑓 = Shear strength, m/Lt 2 , Pa φ = Friction angle

Table 1 -Input parameters for base case simulation

 1 incorporated elasto-plastic behavior in modeling shear fracturing in shale reservoirs; this work shows that the effective rock cohesion in naturally fractured shale formations should be less than 1 MPa to justify a reasonable size of stimulated reservoir volume (SRV) induced by shear failure observed in the field. Table1shows all the input parameters for the synthetic base case simulation in this study, which would fit in the context of a ductile shale formation that contains microfractures.

	Input Parameters	Value
	Young's modulus	20 GPa
	Poisson's ratio	0.25
	Fluid viscosity	1 cp
	Tensile strength	2 MPa
	Formation effective permeability	1 md
	Injection rate per unit thickness	0.0005 (m 3 /s/m)
	Specific weight of fluid	9.8 kN/m 3
	Initial pore pressure	20 MPa
	Maximum horizontal stress	42 MPa
	Minimum horizontal stress	37 MPa
	Vertical stress	65 MPa
	Fracture toughness	1 MPa√m
	Carter's leakoff coefficient	2.00E-06 m/√s
	Spurt loss coefficient	2.5E-04 m
	Porosity	0.1
	Friction angle	27°
	Dilation angle	8°
	Cohesion strength	4 MPa

Table 2 -Simulation results of 30-minute injection with different models

 2

		Poroelastic (Sarris and Papanastasiou 2012)	Poroelastic Using Effective Toughness (Yao 2012)	Poro-Elasto-Plastic
	Fracture length (m)	80.6	58.8	50.2
	Fracture width at wellbore (mm)	11.52	16.95	20.67
	Fluid efficiency (%)	75.20	80.72	84.06


Hanyi Wang

Numerical Modeling of Non-Planar Hydraulic Fracture Propagation in Porous Medium using XFEM with Cohesive Zone Method*

Keywords: hydraulic fracturing, XFEM, cohesive zone method, Non-planar, complex fracture, brittle, ductile

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 unconventional resources exploration and development, hydraulic fracture pattern, geometry and associated dimensions are critical in determining well stimulation efficiency. Multi-stranded non-planar hydraulic fractures are often observed in shale gas stimulation sites. Nonplanar fractures propagating from wellbores inclined from the direction of maximum horizontal stress have also been reported. Current computational methods for the simulation of hydraulic fractures generally assume single, symmetric and planar crack geometries, which severely limits the applicability of these methods in predicting complex fracture geometry. In addition, 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 ductile/unconsolidated formations, even in the form of simple planar geometry. In this study, we present a fully coupled nonplanar hydraulic fracture propagation model in permeable medium based on the Extended Finite Element Method (XFEM), and the arbitrary solution-dependent fracture path cam be determined by solving a set of discontinuity equations. The Cohesive Zone Method (CZM), which is able to model crack initiation and growth by considering process zone effects, is implemented at the crack tip to enable the modeling of fracture propagation in both brittle and ductile formations. To illustrate the capabilities of the model, example simulations are presented including ones involving non-planar hydraulic fracture growth that initiated from an unfavorable perforation direction. The results indicate that the model is able to capture the development of non-planar hydraulic fracturing propagation in both brittle and ductile formations. The method represents a useful step towards the prediction of non-planar, complex hydraulic fractures and can provide us a better guidance of design wells and hydraulic fractures that will better drain reservoir volume in formation with complex stress conditions and heterogeneous properties.

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]. Recently, a significant increase in shale gas production has resulted from hydraulic fracturing, which is used to create extensive artificial fractures around (generally cased and cemented) wellbores. When combined with horizontal drilling, the hydraulic fracturing may allow formerly unpractical shale layers to be commercially viable. And understanding fracture initiation and propagation from wellbores is essential for performing efficient hydraulic fracture stimulation treatment.

As a stimulation tool, 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. The mathematical formulation of the problem is represented by a set of nonlinear integro-differential equations. Also, the problem has a moving boundary 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 by other hydraulic fractures. This in turn can significantly alter the relationship between the injection history and the crack geometry that predicted by plannar 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 greatly facilitates the discretization of complex fractures since the finite element mesh is not required to fit crack surfaces. The displacement discontinuity across the crack surfaces approximated through discontinuous enrichment functions. Likewise, the singularity of the linear elastic solution along the crack front is approximated by properly selected enrichment functions [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, 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]Gordeliy and Peirce, 2013).

However, an assumption commonly made in all these proposed XFEM models are that the loss of fluid into the rock is evaluated using a one-dimensional diffusion equation and the mechanics of fracture opening and fluid loss are considered independent and their interaction is not considered, 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 formation, they fail to predict fracture net pressure and geometry with enough certainty in soft/ductile formation, even in the case of simple, planar fracture geometry. Numerous study and surveys have indicated that the net pressure observed in the field in soft/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 weak non-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 order to model hydraulic fracture propagation in soft/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 [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 interconnection of these micro cracks due to damage evolution. 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, it demonstrates 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] 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 pore pressure cohesive zone model gives predictions that are more conservative on fracture length as compared with pseudo 3D and PKN models. [START_REF] Wang | Poroelastic versus Poroplastic Modeling of Hydraulic Fracturing[END_REF] developed a hydraulic fracture model for both poroelastic and poroplastic formations with the cohesive zone method. Their work indicates that plastic deformation during fracturing execution can lead to higher propagation pressure, shorter and wider fracture geometry. It also found that effects of formation plasticity on fracture propagation are mostly controlled by in-situ stress, plastic property and pore pressure, and the effects plastic deformation in the bulk formation cannot be fully represented by imposing an artificial increased toughness at the fracture tip. 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.

In this study, a fully coupled, non-planar hydraulic fracture model was developed and for the first time, XFEM combined with CZM was introduced together to model hydraulic fracture propagation in brittle and ductile formations. 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 determined solution dependent fracture path and CZM is used to model fracture initiation and damage evolution. Example simulations are presented involving non-planar fracture growth and propagation that initiated from a non-favorable perforation angle to illustrate the capabilities of the model.

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. Most 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. 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: 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 material on the cracked surfaces. Darcy's Law is used to describe fluid diffusion in the porous media: 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 was 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: 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 no 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 [START_REF] Zienkiewicz | The Finite Element Method[END_REF]: 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 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 | Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus[END_REF]. It is an extension of the conventional finite element method based on the concept of partition of unity by [START_REF] Melenk | The Partition of Unity Finite Element Method: Basic Theory and Applications[END_REF], which allows local enrichment functions to be easily incorporated into a finite element approximation. The presence of discontinuities is ensured by the 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.

For the purpose of fracture analysis, the enrichment functions typically consist of the near-tip asymptotic functions that capture the singularity around the crack tip and a discontinuous function that represents the jump in displacement across the crack surfaces. The approximation for a displacement vector function u with the partition of unity enrichment is [START_REF] Fries | Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description[END_REF]: where 𝑁 𝐼 (𝑥) are the usual nodal shape functions; the first term on the right-hand side of the above equation, 𝒖 𝐼 , is the usual nodal displacement vector associated with the continuous part of the finite element solution; the second term is the product of the nodal enriched degree of freedom vector, 𝒂 𝐼 , and the associated discontinuous jump function 𝐻(𝑥) across the crack surfaces; and the third term is the product of the nodal enriched degree of freedom vector, 𝒃 𝐼 𝛼 , and the associated elastic asymptotic crack-tip functions, 𝐹 𝛼 (𝑥). The first term on the right-hand side is applicable to all the nodes in the model; the second term is valid for nodes whose shape function support is cut by the crack interior; and the third term is used only for nodes whose shape function support is cut by the crack tip. 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]: 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 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. 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. Therefore, moving cracks are modeled with the cohesive zone method and phantom nodes in this study.

𝒖 = ∑
𝐹 𝛼 (𝑥) = [√𝑟sin θ 2 ,
A key development that facilitates treatment of cracks in an XFEM analysis is the description of crack geometry, because the mesh is not required to conform to the crack geometry. The level set method (Osher and Fronts, 1988), which is a powerful numerical technique for tracking interfaces and shapes, fits naturally with the XFEM and makes it possible to model arbitrary crack growth without remeshing. To fully characterize a crack, two orthogonal signed distance functions are needed. The first,𝜙, describes the crack surface, while the second,𝜓, is used to construct an orthogonal surface so that the intersection of the two surfaces gives the crack front, as illustrated in Fig. 3. Crack growth is modelled by appropriately updating the 𝜙 and 𝜓 functions at each iteration. In the context of XFEM, no explicit representation of the boundaries or interfaces is needed for the crack geometry because they are entirely described by the nodal data.

Fig. 3-Representation of a nonplanar crack by two signed distance functions (Stolarska et al., 2001)

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 soft/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 [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].

The constitutive behavior of the cohesive zone is defined by the traction-separation relation, which include initial loading, initiation of damage, and the propagation of damage leading to eventual failure at the bonded surface. The behavior of the interface prior to initiation of damage is often described as linear elastic in terms of a penalty stiffness that degrades under tensile and/or shear loading but is unaffected 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 proposed by [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. 4. 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 there 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. 4-Linear Traction-Separation law for different modes

Based on loading type, there are three basic crack propagation modes in fracture process. They are named as Mode I (tension, opening), Mode II (shear, sliding), Mode III (shear, tearing), respectively [START_REF] Simpson | Failure Mechanics Of Anisotropic Shale During Brazilian Tests[END_REF]. In generally, a crack can generate according to any of these modes or combination of them. 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 (Kanninen and Popelar 1985): where E is Young's modulus of formation and ν is Poisson's ratio. 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 soft/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.

G I c = K 𝐼𝐶 2 𝐸 ( 1
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 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. 4, the evolution of the damage variable, D, reduces to ( 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] Belytschko | Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus[END_REF] is implemented in this model 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 where

{ 〈𝑡 𝑛 〉 𝑡 𝑛 0 } 2 + { 𝑡 𝑠 𝑡 𝑠 0 } 2 + { 𝑡 𝑡 𝑡 𝑡 0 } 2 =
G c = G I c + (G II c -G I c ) ( G shear G
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 XFEM-based cohesive zone method can be used to simulate crack initiation and propagation along an arbitrary, solution-dependent path in brittle and ductile materials, since the crack propagation is not tied to the element boundaries in a mesh. In this case the near-tip asymptotic singularity is not needed, and only the displacement jump across a cracked 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 (Van de Meer and Sluys, 2009), as illustrated in Fig. 5. 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. 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 by [START_REF] Song | A Method for Dynamic Crack and Shear Band Propagation with Phantom Nodes[END_REF] and [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.

Simulation Model

Because of well stability problems, isolating well from unwanted zones and several other operational considerations, usually wellbores are cased and then perforated; therefore, majority of hydraulic fracturing treatments are performed through perforations. Hydraulic fractures initiate and propagate along a preferred fracture plane (PFP), which is the path of least resistance. In most cases, stress is largest in the vertical direction, so the PFP is vertical and is perpendicular to the minimum horizontal stress. Because perforations control the initial onset of fracture, if the perforations are misaligned with the direction of PFP, the fracture will reorient itself to propagate in the direction of least resistance. This can result in a tortuous path of smaller width, which can result in early screenout from proppant bridging and, eventually, not optimal stimulation treatments. Fig. 5 shows the reoriented fractures due to unfavorable perforations.

Fig. 5-Reoriented fractures due to unfavorable perforations (Almaguer et al., 2002).

In order to demonstrate our XFEM based CZM model's capability of capturing non-planar hydraulic fracturing propagation, a near wellbore region with fractures initiated from two different perforation angles are simulated, as shown in Fig. 6. The fracture path at the perforation points are defined as initially open to allow entry of the fluid at the perforation tunnel(s), so that the initial flow and fracture growth are possible. A dynamic pressure that equals to the fluid pressure at the perforation entry is imposed on the inner surface of wellbore, which will be implicitly calculated by the model during analysis. 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. The whole simulated regime is saturated with reservoir fluid and incompressible Newtonian fluid is injected at a constant rate at the perforation entries.

Fig. 6-Graphical representation of wellbore, perforations and stresses

Besides modeling fracture propagation in brittle formation that undergoes purely elastic deformation, fracture propagation in soft/ductile formation will be also included in this study. 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 here. 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 solution procedure is schematically illustrated in Fig. 6. The coupled system of non-linear equations 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 finite element package ABAQUS. Because the fracture has to propagate across an entire element at a time during the simulation, mesh size has to sufficient small enough to minimize this effect. In this study, the mesh size monotonically decreases from the outer boundary (0.05 m) to the wellbore (0.004 m).

Results and Analysis

In this section, simulation results are presented to demonstrate the capability of XFEM based CZM to model non-planar hydraulic fracture propagation near the wellbore region, when the direction of perforation is not aligned with the PFP. In addition, fracture propagation in both brittle (elastic rock properties) and ductile (elasto-plastic rock properties) formations are also examined. Fig. 7 shows the fracture first initiated along the direction of perforation, and then it gradually changes its propagation direction to align itself with the direction of PFP until it hits the simulation boundary. It can be also observed that the shear insitu stress is intensified in front of fracture tip, where a two wing pattern, shear zone is developed due to local stress disturbance by propagating fracture. It should be mentioned that the displace field presented in all the following results are enlarged by a scale factor of 50, so the fracture geometry and formation deformations can be seen more clearly in the figures. Fig. 8 shows the formation pore pressure distribution at two different stages. The left figure depicts a moment before hydraulic fracture initiates, when pressure inside the wellbore and perforation channel gradually builds up due to the constant injection of fracturing fluid. Then the fluid begins diffuse into the formation, which leads to higher pore pressure around the wellbore. The right figure depicts a moment when fracture is propagating. It can be observed that a high pressure zone is developed at the fracture tip and a low pressure zone is developed in in front of the fracture tip. The higher pressure zone can be explained by the fact 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, due to a large pressure difference between the fluid inside the fracture and the surrounding formation pore fluid. The low pressure zone is developed because of the newly generated fracture volumes by the microcracks (not completely damaged cells) in front of the fracture tip, where it cannot be filled by the fracture fluid. 9 shows the leak-off rate and cumulative leak-off volume through the surfaces of fractured cells. From the left figure, it can be noticed that the leak-off rate is higher in the newly fractured cells and lower when it is closer to the wellbore (except for those cells that affected by the simulation boundary conditions). This is mainly because the pressure difference between the fluid inside the fracture and the surrounding pores is higher within the newly fractured region. However, as time goes on, leakoff rate will gradually decline because of the increase of surrounding pore pressure due to fluid diffusion. From the right figure, it can be observed that the total cumulative leak-off volume is higher within the fractured cells that closer to the wellbore, this is due to the fact that these cells are fractured earlier, so the exposure time the fracture surface to the fracture fluid is longer and more fracture fluid can leak-off though these surfaces. 11 shows the wellbore pressure and fracture width at the wellbore during hydraulic fracture propagation in both elastic and elasto-plastic formations. It can be observed that the wellbore pressure is almost the same in both formations at initial stage, but pressure in plastic formation gradually surpass the pressure in elastic formation. This is because the plastic deformation region is small initially, so its impact on fracture breakdown pressure is limited, however, as fracture propagates deeper into the formation, more area deformations plastically and in turn, the fracture geometry and wellbore pressure are impacted, this is also the reason for the wider fracture width, that is observed in plastic formation. It can also be noticed that the simulation stops around 16 seconds in the elastic formation, but it runs 1 more second in plastic formation, before the fracture tip hits the simulation boundary. This can be explained by the fact that with wider fracture geometry in the plastic formation, the fracture length is shorter, with the same fluid injection rate in both formations. For more details on how the plastic deformations can impact fracture geometry and fracture propagation pressure beyond the near wellbore region, readers can refer to [START_REF] Wang | Poroelastic versus Poroplastic Modeling of Hydraulic Fracturing[END_REF]. 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. Even though smaller mesh size and simulation time interval can smooth the curves, it can also exert too much burden on computational efforts. So in this paper, the discretization of time and space is a compromise between solution accuracy, numerical convergence and computation time, and the simulation results presented here are still valid to draw general conclusions. Even though the model developed in this study can be applied to a broader range of hydraulic fracture applications, such as to investigate fracture initiated from deviated wells, re-fracture process, multiple fractures propagation and interactions under stress shadow, fracture propagation in heterogeneous formations, etc. The investigation of this paper is focusing on the fracture initialization and development in the near wellbore region, to demonstrate the capability of this presented model in capturing the non-planar fracture propagation in both brittle and ductile formations. More applications and case studies of this fully coupled numerical model on hydraulic fracture treatment design and analysis will be investigated in future research efforts.

Conclusions

In this study, a fully coupled, non-planar hydraulic fracture numerical model was developed and the concept of XFEM combined with CZM was introduced together to model hydraulic fracture process in permeable medium. 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. Example simulations indicated that the presented model is able to capture the initialization and development of non-planar hydraulic fracturing propagation in both brittle and ductile formations. The model developed in this study represents a useful step towards the prediction of non-planar, complex hydraulic fracture evolution and can provide us a better guidance of design wells and hydraulic fractures that will better drain reservoir volume in formations with complex stress conditions and heterogeneous properties.

Fig. 1 -

 1 Fig. 1-Sketch of a plane-strain 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): 𝒒 𝒇 = -𝑤 3 12𝜇 ∇𝑝 𝑓 , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (1)

  ∇𝑝 𝑚 , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (4)

σ

  𝑖,𝑗 = 𝜎 𝑖,𝑗 ′ + 𝛼𝑝 𝑚 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (5)

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (6)

  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)

  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 , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (9)

  -𝜈 2 ), . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (11)

  1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (12) 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (13)

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (16)

Fig. 4 -

 4 Fig. 4-The principle of the phantom node method (Dassault Systè mes, 2014)

Fig. 6 -

 6 Fig. 6-Solution flow diagram for fracture propagation procedure.

Fig. 7 -Fig. 8 -

 78 Fig. 7-Shear stress developed during hydraulic fracture propagation.

Fig. 9 -

 9 Fig. 9-Leak-off rate and cumulative leak-off volume through fractured cells Fig. 10 shows the plastic strain and plastic deformation area during hydraulic fracture propagation when formation plasticity is considered in soft/ductile formations. It can be observed that the closer to the fracture surface, the more server the plastic deformation, and plastic deformation area becomes larger when the fracture aligned itself with the direction of PFP. The reason lies behind the plastic deformation is the induced shear stress in front the fracture tip, as depicted by Fig 7.

Fig. 10 -

 10 Fig. 10-Plastic strain and plastic deformation during hydraulic fracture propagationFig.11shows the wellbore pressure and fracture width at the wellbore during hydraulic fracture propagation in both elastic and elasto-plastic formations. It can be observed that the wellbore pressure is almost the same in both formations at initial stage, but pressure in plastic formation gradually surpass the pressure in elastic formation. This is because the plastic deformation region is small initially, so its impact on fracture breakdown pressure is limited, however, as fracture propagates deeper into the formation, more area deformations plastically and in turn, the fracture geometry and wellbore pressure are impacted, this is also the reason for the wider fracture width, that is observed in plastic formation. It can also be noticed that the simulation stops around 16 seconds in the elastic formation, but it runs 1 more second in plastic formation, before the fracture tip hits the simulation boundary. This can be explained by the fact that with wider fracture geometry in the plastic formation, the fracture length is shorter, with the same fluid injection rate in both formations. For more details on how the plastic deformations can impact fracture geometry and fracture propagation pressure beyond the near wellbore region, readers can refer to[START_REF] Wang | Poroelastic versus Poroplastic Modeling of Hydraulic Fracturing[END_REF]. 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. Even though smaller mesh size and simulation time interval can smooth the curves, it can also exert too much burden on computational efforts. So in this paper, the discretization of time and space is a compromise between solution accuracy, numerical convergence and computation time, and the simulation results presented here are still valid to draw general conclusions.

Fig. 11 -

 11 Fig. 11-Fracture width and wellbore pressure for elastic and plastic formations

  off coefficient, 𝐿 2 𝑡/𝑚, m/Pa•s 𝑑 = Initial thickness of cohesive surfaces, L, m 𝐸 = Young's modulus, m/L𝑡 2 , Pa 𝒇 = Body force per unit volume, m/L 2 /t 2 , N/m 3 𝐺 = Potential function 𝐺 I 𝑐 = Fracture critical energy in Mode I, mL 2 /t 2 , J 𝐺 II 𝑐 = Fracture critical energy in Mode II, mL 2 /t 2 , J 𝐺 III 𝑐 = Fracture critical energy in Mode III, mL 2 /t 2 , J 𝐈 = Unit Matrix 𝒌 = Formation permeability tensor in formation, L 2 , md 𝐾 𝑛 = Cohesive stiffness, m/L 2 𝑡 2 , Pa/m 𝑝 𝑓 = Fluid pressure inside fracture, m/L𝑡 2 , Pa 𝑝 𝑚 = Pore pressure in formation matrix, m/L𝑡 2 , Pa 𝒒 𝒇 = Flow rate along fracture per unit height, 𝐿 3 /𝑡, 𝑚 3 /𝑠 𝑞 𝑙 = Leakoff rate per unit fracture surface area, 𝐿/𝑡, 𝑚 3 /(𝑠 • 𝑚 2 ) 𝒒 𝒎 = Fluid flux velocity in formation matrix, L/t, m/s 𝒕 = Surface traction per unit area, m/L 3 t 2 , pa/m 2 𝑇 = Traction force, m/L𝑡 2 , Pa 𝑇 𝑀𝑎𝑥 = Cohesive tensile strength, m/L𝑡 2 , Pa 𝑉 𝐿 = Cumulative leakoff volume, 𝐿 3 , 𝑚 3 𝑤 = Fracture width, L, m 𝛼 = Poroelastic constant 𝛾 = toughness in Mode I, m/√Lt 2 , Pa√m 𝜇 = Fluid viscosity, m/Lt, cp 𝜈 = Poisson's ratio 𝑣 𝐿 = Leakoff rate, L/t, m/s 𝜌 𝑓 = Fluid mass density, m/L 3 , kg/m 3 σ = Stress, m/L𝑡 2 , Pa 𝜎 ′ = Effective stress, m/L𝑡 2 , Pa 𝜏 𝑓 = Shear strength, m/L𝑡 2 , Pa ∅ =Porosity of porous media

  𝒇 -∂𝑤 ∂𝑡 + 𝑞 𝑙 = 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)

  𝑁 𝐼 (𝑥)[𝒖 𝐼 + 𝐻(𝑥)𝒂 𝐼 + ∑ 𝐹 𝛼

	4
	(𝑥)𝒃 𝐼 𝛼 ]
	𝛼=1

𝑁 𝐼=1 , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (8)

Table 1

 1 shows all the input parameters for the simulation.

	Input Parameters	Value
	Wellbore diameter	0.2 m
	Boundary Length	2 m
	Perforation Length	0.1 m
	Perforation Angle	60 o
	Elastic modulus	20 GPa
	Poisson's ratio	0.25
	Fluid viscosity	1 cp
	Tensile strength	2 MPa
	Formation permeability	50 md
	Flow rate	0.00005 (m 3 /s/m)
	Specific weight of fluid	9.8 kN/𝑚 3
	Initial pore pressure	10 MPa
	Maximum horizontal stress	23 MPa
	Minimum horizontal stress	20 MPa
	Vertical stress	30 MPa
	Critical Fracture Energy	100 J/m 2
	Leak-off coefficient	5E-8 m 3 /s/Pa
	Porosity	0.2
	Friction angle	27 o
	Dilation angle	8 o
	Cohesion strength	3 MPa
	Table 1-Input parameters for simulation