A three-phase XFEM model for hydraulic fracturing with cohesive crack propagation

Computers and Geotechnics 69 (2015) 82–92

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

A three-phase XFEM model for hydraulic fracturing with cohesive crack propagation Saeed Salimzadeh a, Nasser Khalili b,⇑ a b

Department of Earth Science and Engineering, Imperial College, London, United Kingdom School of Civil and Environmental Engineering, The University of New South Wales, Sydney 2052, Australia

a r t i c l e

i n f o

Article history: Received 24 November 2014 Received in revised form 7 April 2015 Accepted 5 May 2015

Keywords: Coupled XFEM formulation Hydraulic fracturing Cohesive crack Three-phase

a b s t r a c t A three-phase hydro-mechanical model for hydraulic fracturing is proposed. Three phases include: porous solid, fracturing fluid and host fluid. Discontinuity is handled using extended finite element method (XFEM) while cohesive crack model is used as fracturing criterion. Flow through fracture is defined as one-dimensional laminar flow, and flow through porous medium (host reservoir) is defined as two-dimensional Darcy flow. Coupling between two fluids in each space, fracture and pore, is captured through capillary pressure–saturation relationship, while the identical fluids in fracture and pore are coupled through a so-called leak-off mass transfer term. Coupling between fluids and deformation is captured through compatibility of volumetric strain of fluids within fracture and pore, and volumetric strain of the matrix. Spatial and temporal discretisation is achieved using the standard Galerkin method and the finite difference technique, respectively. The model is verified against analytical solutions available from literature. The leaking of fracturing fluid into the medium and suction of porous fluid into the fracture around the tip, are investigated. Sensitivity analyses are carried out for cases with slow and fast injection rates. It is shown that the results by single-phase flow may underestimate the leak-off. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Hydraulic fracturing is a complex multi-physics problem. A fluid is injected into the host rock; the fracture initiates and propagates due to the induced hydraulic loading. The simulation of the process requires several components: (i) flow of the fracturing fluid within the fracture; (ii) the mechanical deformation of the porous medium due to the applied hydraulic load; (iii) the fracture propagation; (iv) the leak-off flow to the surrounding porous media (host rock); and (v) the fluid flow through the host medium. The first three components are widely addressed in the literature, while last two are considered only in few studies albeit with some restrictions. The flow through planar fracture is commonly modelled using lubrication theory [5] for an incompressible Newtonian fluid obeying cubic law. The lubrication equation is a nonlinear partial differential equation that relates the fracture width and the fracture pressure gradient. The fracturing flow model can also be extended to include power-law fluids [1,12] and slightly compressible fluids. However, in these models the fracture is assumed to be fully saturated, involving a single-phase fluid flow, which is not consistent

⇑ Corresponding author. http://dx.doi.org/10.1016/j.compgeo.2015.05.001 0266-352X/Ó 2015 Elsevier Ltd. All rights reserved.

with the actual field conditions where the fracturing fluid can have properties substantially different than those of the host fluid. Furthermore, depending on the in situ stresses the fracturing fluid front and the fracture front may not coincide, giving rise to the so-called fluid lag phenomenon [19,1,22]. In the lag region, the host fluid enters the fracture under the developed suction. Therefore, a second fluid is required to complete the fracture flow model. To incorporate discontinuity in a rock mass, numerous numerical models have been proposed in the literature using a variety of techniques such as discrete element [11]; boundary element method [3]; finite element method with remeshing [10,18]; and the extended finite element method-XFEM [6,33]. In XFEM, finite element space is enriched by additional functions inspired from the analytical or asymptotic solution of the problem. The enrichment is performed by activating extra degrees of freedom in nodes surrounding the discontinuity. XFEM has been applied to many engineering problems involving three-dimensional cracks, dynamic crack propagation, formation of shear bands, cohesive cracks, steady-state thermo-elastic fracture mechanics as well as temperature discontinuities [44,36,4,16,42]. However, only few attempts have been made to apply XFEM to coupled flow-deformation analysis. de Borst and co-workers [13,37,38] proposed a method for incorporating discontinuity in the pressure field; the flow field was enriched in a similar manner to the

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

displacement field, but with a different enrichment function. de Borst et al. [13] used the Heaviside step function as the enrichment function while Rethore et al. [37] proposed the ‘‘distance to the discontinuity’’ as the enrichment function for the pressure field. Rethore et al. [38] extended the model to unsaturated media with the passive gas phase, which was later expanded to active gas phase by Mohammadnejad and Khoei [35]. Nevertheless, in all these contributions the fluid flow through fracture was either ignored [13] or modelled using laminar flow [37,28] or Darcy’s relation with smeared fracture properties [38,35]. In addition, the flow gradient along the fracture was based on the enriched fracture pressures which may not represent the actual fluid pressure within the fracture as multiple enriched pressures are computed to represent the pressure at a single fracture node. The Biot’s coefficient was also applied for both porous media and fracture pressures which is not appropriate for the fracture pressures. In reality, fracture pressures are felt along the fracture faces in their entirety, without reduction from Biot’s coefficient. Dahi Taleghani [12] used XFEM to model the fracture discontinuity, but simulated fracture flow through a one-dimensional laminar flow. The coupling between fracture flow and mechanical model was achieved weakly through a successive procedure where the results of fracture flow were used to update the solution from the mechanical model and vice versa. Furthermore, the leak-off flow from the fracture to the porous media was ignored similar to the contributions of Lecampion [31], Weber et al. [45] and Gordeliy and Peirce [22]. Furthermore, these models were developed for single-phase fracture flow. In general, little attention is given in the literature to the leak of fracturing fluid into the host rock, particularly in the presence of two fluids within the fracture; i.e. the fracturing fluid and the host fluid in the lag region. The leak-off has a detrimental effect on the hydraulic fracturing process as it decreases the efficiency of the process by decreasing the fracture volume [17]. A substantial loss of fracturing fluid can occur in high permeability reservoirs due to leak-off [1]. The common leak-off model used is the one proposed by Carter [9] which implies that the leak-off is one-dimensional in the direction orthogonal to the fracture plane and into the infinite homogenous host medium. This model tends to underestimate the leak-off due to its underlying assumption of single phase 1D fluid flow normal to the fracture planes [24,32,30]. Other analytical methods have been proposed to calculate leak-off flow [23,32], however they are usually based on over-simplified and restrictive assumptions as to the shape of the fracture and prescribed boundary conditions. In case of two-phase flow, to authors’ best knowledge, no attempt has been made to investigate the leak-off in partially saturated medium. In this paper, a fully coupled three-phase model for hydraulic fracturing in porous media is proposed. Two fluids: fracturing fluid and host fluid are considered. The independent fluid flow for the porous host medium makes the present model capable of simulating 1D, 2D and 3D leak-off fluid flow through high/low permeable, homogeneous/inhomogeneous host medium with transient pressure boundary condition at the fracture-rock boundary. The fracture discontinuity in the deformation model is handled using XFEM technique, while cohesive crack model is used as fracture propagation criteria. Coupling between two fluids in each space, fracture and pore, is captured through capillary pressure–saturation relationship, while the identical fluids in fracture and pore are coupled through a leak-off mass transfer term. Particular attention is given to the flow simulation of two fluids: flow of fracturing fluid into the porous medium and its interaction with the non-wetting host fluid as well as the suction of the host fluid into the lag zone ahead of the fracturing fluid front. Although, under specific conditions i.e. significant in situ stresses, the size of lag

83

zone can be negligible and the fracture may be assumed to be fully filled by the fracturing fluid, however, for the sake of completeness, the proposed formulation includes two fluids within the fracture. Coupling between fluids and mechanics is satisfied through the compatibility of volumetric strain of fluids within fractures and pores, and volumetric strain of the matrix. Spatial and temporal approximations of the governing equations are achieved using the standard Galerkin method and finite difference technique, respectively. The model is verified against several benchmark solutions from the literature. The salient features of the model are demonstrated through simulation of hydraulic fracture process in porous medium. 2. Governing equations 2.1. Conceptual model The model proposed comprises three interacting models: deformation model, flow model for porous medium (host flow model) and fracture flow model. The deformation model is based on the theory of linear elasticity for the stress range of interest, satisfying the conditions of equilibrium. The fracture discontinuity is handled using XFEM technique in conjunction with the cohesive fracture mechanics. The cohesive crack model is preferred in this study due to the ease of implementation and the ability to capture the analytical solution of KGD model [21]. However, the proposed model can also be used in conjunction with the Linear Elastic Fracture Mechanics (LEFM). The host flow model is based on the Darcy’s law which is combined with the mass balance equation of wetting, w, and non-wetting fluids, nw i.e. fracturing fluid and host fluid [39]. The fracture flow model is based on the lubrication theory which is extended to slightly compressible fluids. The two flow models are coupled through leak-off mass transfer term. The host flow is considered continuous across the fracture discontinuity with discontinuous flow velocity at the fracture. The discontinuity in the flow velocity is imposed through leak-off term. The leak-off is defined as a function of the permeability of the porous medium with respect to each fluid, dynamic viscosity of the fluid constituent and fluid pressure gradient between fracture and surrounding porous medium. The deformation model and host flow model are coupled through effective stress equation, while fracture flow model and deformation model are coupled through compatibility of volumetric strain of fluids within fracture and the fracture width, as well as through the hydraulic loading on the fracture planes. 2.2. Deformation model The deformation model is expressed satisfying the condition of equilibrium on a representative elementary volume (REV) of the porous medium. For quasi-static conditions, the linear momentum balance equation for this elemental volume may be written as

div r þ F ¼ 0

ð1Þ

where r is the total external stress, and F is the body force per unit volume. The stress and strain relationship of the element is in turn expressed as

dr0 ¼ Dde

ð2Þ 0

in which r is the effective stress, D is the tangent drained stiffness matrix, e is the strain of the fractured porous medium. The effective stress for a porous medium saturated by two fluids, i.e. wetting ðwÞ and non-wetting ðnwÞ fluids, is expressed as

  dr0 ¼ dr þ bw dppw þ bnw dppnw I

ð3Þ

84

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

in which, b is the incremental effective stress parameter, p is the fluid pressure and I is the second-order identity tensor. The incremental effective stress parameters will be quantified as [40]

 cs  bw ¼ w 1  c  cs  bnw ¼ ð1  wÞ 1  c

@ ðsvÞ @s

dt ¼

ð5Þ

in which T represents the tangential modulus matrix of the discontinuity. It is noted that only cohesive tractions normal to the discontinuity ðtn Þ is considered in this study and the shear cohesive traction and shear relative displacement along the fracture discontinuity are considered zero. In addition to cohesive traction, in fluid saturated fracture discontinuities, the hydraulic loading is applied to the fracture faces. The hydraulic loadings can be either positive for pressurised fluid within the fracture or negative due to suction in the lag region. Therefore, the natural boundary conditions for loading at the fracture discontinuity is expressed as

ð6Þ

where, s is the suction (capillary pressure). Assuming infinitesimal deformations, strain is related to displacement by

1 2

e ¼ ð r u þ ur Þ

ð7Þ

where u donates the displacement vector of the porous medium. Combining Eqs. (1)–(3) and (7) yields the differential equation describing the deformation field for a continuum porous medium saturated by wetting and non-wetting fluids

  1 div Dðrdu þ durÞ  bw dppw I  bnw dppnw I þ dF ¼ 0 2

ð8Þ

Following the cohesive fracture mechanics, the nonlinear behaviour of the fracturing medium in the cohesive zone is governed by a traction-separation law as

t ¼ tðwÞ

ð9Þ

where t is the cohesive traction and w is the fracture width, the difference in the displacement vector between two faces of the fracture discontinuity. The fracture width is in turn a function of discontinuous part of the displacement field (also known as enriched part within extended finite element method-XFEM)

w ¼ f ðuE Þ ¼ uþE þ uE

@tn @tn @w duE ¼ TduE dw ¼ @w @w @uE

ð4Þ

in which, w is the unsaturated incremental effective stress parameter; c and cs are compressibility coefficients of medium and solid grains respectively. Incremental effective stress parameters are subjected to the constraint, bpw þ bpnw ¼ 1  ccs , known as Biot’s coefficient. Khalili et al. [26] showed that the unsaturated incremental effective stress parameter, w, may be approximated from the total form of the effective stress parameter for unsaturated soils, v, as

w¼

transferred across the cohesive zone a decaying function of the fracture width. Assuming mode I failure, the linearised traction-separation law can then be written as

f nc dr:nc ¼ dt  dp

−

+

ð12Þ

f is the average fluid pressure within the fracture, and nc is where p the outward unit vector normal to the discontinuity (Fig. 1). Noticed that in the above equation Biot’s coefficient is not applied to fracture loading, and that the fluid pressure in the fracture is treated identical to an external pressure applied to fracture planes. Such a pressure does not require scaling according to Biot’s coefficient. 2.3. Fracture flow model Independent fracture flow model is considered for fracture discontinuity in this research. Fracture pressures are calculated directly, which may be linked to multiple pore pressures within the surrounding porous medium. The objective is to obtain a more realistic representation of fracture flow compared to current models presented in the literature in which gradient of fracture flow is calculated entirely based on the enriched component of the fluid pressures [13,37,38,28]. Assuming the fracture length is much larger than the fracture width, the average velocity of fluid n along the fracture is calculated using the cubic law as [47]

v f n ¼ kfrn

ð10Þ

in which, uE is the discontinuous (enriched) part of displacement at the fracture discontinuity and superscription signs + and  show two sides of the discontinuity (Fig. 1). When the failure limit of a quasi-brittle material is exceeded, the material exhibits a softening behaviour. This is captured by rendering the cohesive traction

ð11Þ

w2 @pf n 12ln @l

ð13Þ

where w is the fracture width, l is the fluid dynamic viscosity, kfr is the relative permeability, pf is the fracture fluid pressure, and l is the dimension along the fracture. Notice that the conductivity coef  2 ficient kfrn 12wl may be substituted by other theoretically or expern

imentally adopted relationships. If the size of the lag zone is negligible, the fracture can be assumed to be fully saturated by the fracturing fluid, so the relative permeability in the above equation is equal to 1. The fracture width is equal to the discontinuous displacement across the fracture. The mass balance equation for a slightly compressible fracture fluid n may hence be written as

 @   @ q Sf n v f n w þ qn Sf n w þ Ln ¼ 0 @l n @t

ð14Þ

in which, q is the fluid density, Sf n is the degree of saturation of the fluid n and L is the source/sink term due to leak-off. The leak-off leads to the mass transfer coupling between the fracture and host (porous) media flow. Assuming that the fracture fluid has a Newtonian behaviour, using Darcy’s law the leak-off flow per unit area of the fracture wall can be written as

Ln ¼ qn Fig. 1. Schematic representation of the problem with the body X, boundary C, fracture boundary Cc , fracture pressure, pf , and pore pressure pp .

kpn

ln

rpn

ð15Þ

where kpn and l are the permeability of porous media with respect to fluid n and fluid dynamic viscosity, respectively; rp is the

85

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

pressure gradient between fracture and porous medium, and is defined as

rpn ¼



pf n  ppn y



ð16Þ

in which, y is the distance and pp is the fluid pressure in the porous medium. Substituting fluids velocity into mass balance equation and after some manipulation we have

!   @pf n @Sf n kpn pf n  ppn @ w @pf n @w Sf n kfrn þ Sf n þw þ ¼ Sn wcn @l @t 12ln @l @t @t ln y 3

ð17Þ where, cn is the fluid compressibility. Assuming that the fluid constituent is barotropic, we have

@pf n @ qn ¼ qn c n @t @t

ð18Þ

For the case of single incompressible fluid with no leak-off, i.e. Sf n ¼ 1; kfrn ¼ 1; cn ¼ 0 and Ln ¼ 0, the Eq. (17) reduces to lubrication equation [5]

  @pf 1 @ @w ¼ w3 12l @l @t @l

@Sf n @Sf n @sf ¼ @t @sf @t

ð20Þ

in which, sf ¼ pfnw  pfw is suction within the fracture. The final governing equations for the flow of wetting and non-wetting fluids within fracture discontinuity may be written as

  @pfw @pfnw @ w3 @pfw @w ¼ Sfw Sfw kfrw þ kfw;fw þ kfw;fnw @l @t 12lw @l @t @t   kpw pfw  ppw þ lw y   @pfw @ w3 @pfnw @w ¼ Sfnw Sfnw kfrnw þ kfnw;fw þ kfnw;fnw @l @t 12lnw @l @t   @pfnw kpnw pfnw  ppnw  þ @t lnw y

kfnw;fnw

@Sfw @sf

@Sfw ¼ Sfnw wcnw  w @sf

kfw;fnw ¼ kfnw;fw ¼ w

@Sfw @sf

ð21Þ

ð22Þ



@u @t

vs ¼

ð28Þ

where un and u are the displacement vectors of the fluid n and solid skeleton, respectively. The mass balance equation of the fluid n within porous medium may be written as

  @   div qn nn v n þ qn nn  Ln ¼ 0 @t

ð29Þ

Combining the equation of linear momentum balance (Darcy’s law) with the mass balance equation of the fluid within porous medium, and with some manipulation, the flow model for unsaturated porous media can be written as [27,41]

div



kpw 

lw

rppw þ qw g

 kpnw 

lnw





@ppw @ðdiv uÞ þ kpw;pw @t @t   @ppnw kpw ppw  pfw ð30Þ þ þ kpw;pnw @t lw y

¼ bw



rppnw þ qnw g



@ppw @ ðdiv uÞ þ kpnw;pw @t @t   @ppnw kpnw ppnw  pfnw þ þ kpnw;pnw ð31Þ @t lnw y

¼ bnw

  @Spw kpw;pw ¼ nSpw cw þ bw  nSpw cs  n @sp

ð32Þ

  @Spw kpnw;pnw ¼ nSpnw cnw þ bnw  nSpnw cs  n @sp

ð33Þ

@Spw @sp

ð34Þ

where, n is the porosity and Spn is the degree of saturation of fluid n in the porous medium. The pressure field is assumed continuous across the fracture discontinuity, with the discontinuity in the Darcy flow velocity normal to the fracture taken into account through leak-off loading defined in (15). 3. Numerical model 3.1. Discretisation in space

ð24Þ ð25Þ

rppn þ qn g ln where v r is the relative velocity vector, kpn

and

ð23Þ

A mathematical description of fluid flow through host porous medium can be obtained by combining linear momentum balance equation (Darcy’s law) with the mass conservation equation of the fluid. Neglecting the inertial and viscous effects, the equation of linear momentum balance (Darcy’s law) for the fluid n in porous media may be written as

kpn 

@un @t

kpw;pnw ¼ kpnw;pw ¼ n

2.4. Host flow model

v rn ¼ 

vn ¼

div

The term on the RHS of (17) may be written in terms of suction (capillary pressure) as

ð27Þ

where nn and v n are the volume fraction and the absolute velocity of the fluid n, and v s is the solid skeleton velocity. The absolute velocities v n and v s are defined as

ð19Þ

@Sf n @t

kfw;fw ¼ Sfw wcw  w

vrn ¼ nn ðv n  v s Þ

ð26Þ

is the porous medium permeability tensor with respect to the fluid n and g is the vector of gravitational acceleration. The relative velocity of the fluid n with respect to moving solid is given by

Spatial discretisation has been performed using the standard Galerkin method with displacements and fluids pressures defined as primary variables. Extended finite element method (XFEM) is used to model the fracture discontinuity within displacement field at element level. Arbitrary discontinuities can be handled using XFEM without the need for remeshing. In XFEM, the displacement field consists of two parts, one continuous ðuC Þ and the other discontinuous ðuE Þ

u ¼ uC þ uE

ð35Þ

The continuous part is the standard displacement field corresponding to the situation without any cracks. The discontinuous displacement field (also known as enriched part) is based on local partitions of unity and allows the element to include a discontinuity. Subscriptions C and E stand for continuous and enriched part of the displacement throughout the paper. Heaviside step function HðxÞ, also known as ‘jump’ function, is typically used to enrich the nodes in fully cracked elements if their support was cut by the crack into two disjoint pieces [33]. HðxÞ is given as [42]

 HðxÞ ¼ signððx  xc Þ:nc ÞsignðxÞ ¼

þ1 x > 0 1 x < 0

ð36Þ

86

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

where, xc is the nearest projection of point x onto the fracture line, and nc is the outward unit normal to the fracture at point xc . Using standard Galerkin method, the displacement and pressures within an element may be approximated from nodal values as

bC uC ¼ N u

ð37Þ

bE uE ¼ N E H ð x Þ u b ppn ¼ N p pn

ð38Þ

bfn pf n ¼ Nc p

ð40Þ

ð39Þ

where, N; NE and Nc are the standard shape function vector, the shape function for enriched nodes, and shape function for the fracbC; u b E; p b p and p b f are the nodal values of continuous and disconture, u tinuous displacement, porous media pressure and fracture pressure, respectively. The set of discretised equations can be written as



KCC

KCE

"

KEC KEE þ TEE

#

X b_ C u  _u b E



n¼w;nw

CC CE



!



0 b_ pn þ b_ f n ¼ P_ p p C E fn pn

ð41Þ

" # X  u   b_ C b pn  b pn  Lpn;f n p b f n ¼ Q pn b_ pf  Lpn;pn p  CTC CTE pn Mpn;pf p  Hpn p _u bE f¼w;nw ð42Þ " # X  u   b_ C bfn  b f n  Lf n;pn p b pn ¼ Q f n b_ f f  Lf n;f n p  0 CTE f n Mf n;f f p  Hf n p _u bE f¼w;nw ð43Þ

in which,

Z

KCC ¼

Z

KCE ¼

ZX

KEC ¼

Z

KEE ¼ TEE ¼

X

Z

X

ð44Þ

BT1C DB1E dX

ð45Þ

BT1E DB1C dX

ð46Þ

BT1E DB1E dX

ð47Þ

NTE TNE dC

ð48Þ

BT2C bn NdX

ð49Þ

X

Cc

CCpn ¼ CEpn ¼ CEf n ¼ Hpn ¼

Z Z

X

BT2E bn Nd

ZX Z Z

X

ð50Þ

NTE Sf n Nc d

ð51Þ

BT3

kpn

ln

B3 dX

@Nc T w3 @Nc Sf n kfrn dl @l 12 ln @l Cc Z Mpn;pf ¼ NT kpn;pf NdX ZX NTc kf n;f f Nc dC Mf n;f f ¼ Cc Z kpn NT NdX Lpn;pn ¼ ln y Cc Z kpn NT N dX Lpn;f n ¼ ln y c Cc Z kpn Lf n;f n ¼ NTc Nc d C l Cc ny Z kpn Lf n;pn ¼ NTc NdC l Cc ny Hf n ¼

@N @x

3

0

6 0 B1C ¼ 6 4

@N 7 7 @y 5

@N @y

ð60Þ

@N @x

2 @NE HðxÞ 6 B1E ¼ 6 4

ð52Þ ð53Þ

@NE HðxÞ 7 7 @y 5

0 @NE HðxÞ @y

ð56Þ ð57Þ

ð61Þ

@NE HðxÞ @x

in which, x and y are horizontal and vertical coordinates, respectively. 3.2. Discretisation in time To obtain numerical solution of the governing equations, the rate form of the discretised equations is integrated over the time domain. The time integration is performed using a generalised trapezoidal method such that for an arbitrary function f over a time interval Dt one can write tþDt

    t tþDt t f ðtÞdt ¼ ð1  hÞf þ hf Dt ¼ h Df þ f Dt

ð62Þ

t tþDt

where the f and f are values of function f in time t and t þ Dt . The time integration parameter h, ranging between 0 and 1, defines the type of integration scheme. h ¼ 0 gives forward, h ¼ 1 backward, h ¼ 0:5 Crank–Nicolson and h ¼ 2=3 Galerkin scheme. The Galerkin scheme provides unconditional stability and first order of accuracy. It also provides better stability than the Crank–Nicolson scheme in spite of the lower accuracy. With Crank–Nicolson scheme the solution can display oscillatory response, particularly in the vicinity of the perturbation [25,29]. Applying Eq. (62) to Eqs. (41)–(43) yields



KCC

"

KCE

KEC KEE þ TEE

# !

X CC

bC 0 Du b pn þ b f n ¼ DP  Dp Dp bE CE f n CE pn Du n¼w;nw

ð63Þ   CTC

CTE

"

pn

bC Du bE Du

# b pn   Hpn hDt D p

X

b pf Mpn;pf D p

f¼w;nw

  b pn  Lpn;f n hDt D p bfn  Lpn;pn hDt D p   b tpn Dt þ Lpn;pn p b tpn  Lpn;f n p b tf n Dt þ Q pn Dt ¼ Hpn p

ð54Þ ð55Þ

3

0

@x

t

X

C

Cc

2

Z

BT1C DB1C d

X

K is the element stiffness matrix modified for softening cohesive tractions through TEE ; C is the flow-deformation coupling matrix, H is the fluid flow matrix, M is the coupling mass matrix, L is the leak-off mass matrix, P is the vector of nodal forces, Q is the vector of nodal fluid flux; B1 ; B2 ¼ dT B1 , and B3 ¼ rN, are derivatives of the shape functions. d is the identity vector, r is the gradient operator and l is the dimension along the fracture discontinuity. The dot (_) represents the rate form. The components of matrices B1C and B1E in two-dimensional plane-strain condition will be [7]

"  ½ 0 CE f n

bC Du bE Du

# bfn   Hf n hDt D p

X

ð64Þ

bff Mf n;f f D p

f¼w;nw

  b f n  Lf n;pn hDt D p b pn  Lf n;f n hDt D p   b tf n Dt þ Lf n;f n p b tf n  Lf n;pn p b tpn Dt þ Q f n Dt ¼ Hf n p

ð65Þ

The above equations can be arranged in a matrix form of SX ¼ F in which S is the general stiffness matrix defined as

ð58Þ

2

K

Cp

Cf

ð59Þ

6 CT 4 p

Hp

Lp 7 5

S¼

CTf

LTp

Hf

3 ð66Þ

87

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

in which,

K¼

KCC

KCE

KEC

KEE þ TEE

"

CC pw

CCpnw

CEpw

CEpnw

0

0

CEfw

CEfnw

Cp ¼ " Cf ¼ Hp ¼

ð67Þ # ð68Þ

4. Verification and application

# ð69Þ

Hpw hDt  Mpw;pw  Lpw;pw hDt Mpw;pnw Mpnw;pw Hpnw hDt  Mpnw;pnw  Lpnw;pnw hDt



ð70Þ

Hf ¼

The four-node quadrilateral elements together with two-point Gaussian quadrature rule are used for numerical approximation and numerical integration. Numerical integration of the elements bisected by the fracture discontinuity is done using common method of element partitioning in which numerical integration is done separately within each sub-domain.



Hfn hDt  Mfw;fw  Lfw;fw hDt

Mfw;fnw

Mfnw;fw

Hfnw hDt  Mfnw;fnw  Lfnw;fnw hDt

A three-point bending test is selected to validate the fracture propagation with softening cohesive tractions, followed by hydraulic fracturing examples to test the hydro-mechanical coupling aspects of the model with leak-off feature. For all examples only one fracture is allowed to propagate, however, the model can be extended to simulate multiple fractures and crack branching. 4.1. Three-point bending test

ð71Þ

Lp ¼

LTp ¼



Lpw;fw hDt

0

0

Lpnw;fnw hDt

Lfw;pw hDt

0

0

Lfnw;pnw hDt

ð72Þ

ð73Þ

X is the generalised vector of unknown nodal values

8 9 bC > > Du > > > > > > > bE > > > D u > > > > > > < Dp b pw = X¼ > b pnw > > Dp > > > > > > > > b fw > D > p > > > > : b > ; D p fnw

ð74Þ

and F is the generalised force vector defined as

8 > > > > > > > > > > > > > <

9 > > > > > > > > > > > > > =

D PC D PE

  b tpw Dt þ Lpw;pw p b tpw  Lpw;fw p b tfw Dt þ Q pw Dt Hpw p   F ¼ Hpnw p b tpnw Dt þ Lpnw;pnw p b tpnw  Lpnw;fnw p b tfnw Dt þ Q pnw Dt > > > > > >   > > > > > > b tfw Dt þ Lfw;fw p b tfw  Lfw;pw p b tpw Dt þ Q fw Dt > > H p fw > > > > > > > >   > > > > : Hfnw p b t Dt þ Lfnw;fnw p b t  Lfnw;pnw p bt Dt þ Q fnw Dt ; fnw

fnw

This example is adopted from Wells and Sluys [46]. A simply supported beam is loaded symmetrically at the centre of the beam on the top edge. Young’s modulus and Poisson’s ratio are assumed as E ¼ 100 MPa and t ¼ 0, respectively. Tensile strength is f t ¼ 1:0 MPa and cohesive fracture energy is Gf ¼ 0:1 kN/m. The model dimensions and the mesh are shown in Fig. 2. The normal traction force tn transmitted across a discontinuity is made an exponentially decaying function of the fracture width, w, as

  f tn ¼ f t exp  t w Gf

The traction-separation diagram for this example is shown in Fig. 3. As the stresses vary intensely around the fracture tip, the smoothing of stresses is crucial in order to accurately determine the principal stress direction. The stresses in the vicinity of the fracture tip are averaged using the Gaussian averaging method

R rxdX r ¼ RX X xdX

x¼

  r2 exp  2a2 ð2pÞ3=2 a3 1

ð81Þ

pnw

   Hinþh þ Linþh DtXn þ DPnþh þ Q nþh Dt

2.5

2

1.5

1

ð76Þ 0.5

where

Xinþh ¼ ð1  hÞXin þ hXinþ1 The iterations are repeated until consecutive values of within a specified tolerance e

P -3 3 x 10

Height (m)

The non-diagonal components of stiffness matrix S are populated with the coupling matrices Cf for mechanical-fracture flow coupling, Cp for mechanical-formation flow coupling and Lp for fracture-formation flow coupling. The set of equations in (63)– (65) are highly nonlinear. Nonlinearities arising from fluid pressure dependency of constitutive coefficients are handled using Picard method such that the coefficient matrices appearing in the stiffness matrix S as nonlinear functions of unknown variables X, are updated in every iteration i as

iþ1

i

X

nþ1  Xnþ1 < e

ð80Þ

in which x is the Gaussian weight function that smoothens the stress field at the tip and is defined as

ð75Þ

iþ1 Xnþ1 ¼ S1 nþh

ð79Þ

ð77Þ Xinþ1

agree to

ð78Þ

0

0

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Length (m) Fig. 2. Three-point bending beam. Red line shows the propagated fracture and blue circles show the enriched nodes.

88

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

1

tn (MPa)

0.8 0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

1

Fracture width, w (mm) Fig. 3. Traction-separation law used for three-point bending test.

Vertical Dimension (m)

506 1.20 1.00

Wells and Sluys (2001) Present Study

P (N)

0.80 0.60

504 502 500 498 496 494

0.40

0

4

6

8

10

12

Horizontal Dimension (m)

0.20 0.00

2

Fig. 5. Model used for hydraulic fracturing examples (top) and the mesh configuration (bottom). Red line shows the propagated hydraulic fracture.

0

0.5

1

1.5

u (mm)

The fracturing fluid with dynamic viscosity lw ¼ 1  104 kPa.s

Fig. 4. Load–displacement diagram for three-point bending test.

and compressibility coefficient cw ¼ 1  106 kPa1 is injected into the fracture at the rate Q w ¼ 1 m3/min (slow-injection). The

where r is the distance to the fracture tip and parameter a determines how quickly the weight function decays away from fracture tip. The parameter a is commonly taken equal to three times the element size [46] or approximately 1% of Hillerborg’s characteristic length [15] defined as

lc ¼

E Gf 1  t2 f 2t

ð82Þ

which also identifies the size of the cohesive zone. To properly capture the distribution of tractions in the cohesive zone, a sufficient number of elements should be used to discretise the cohesive zone. Moes and Belyteschko [34] suggested using a minimum of about two elements in the cohesive zone. For this example the size of the cohesive zone is lc ¼ 10 mm. Fig. 4 shows the load–displacement response for this example. Also included in this figure is the result by Wells and Sluys [46]. Very good agreement is found between both results. The area under the load–displacement curve is equal to 0.3026 Nmm which agrees well with the total cohesive fracture energy (0.30 Nmm). 4.2. Hydraulic fracturing examples In these examples, a viscous fluid i.e. fracturing fluid (wetting fluid) is injected into a fracture with initial length of 1 m. The host medium is assumed to be saturated with gas as non-wetting fluid at atmospheric pressure. The model dimensions and mesh configuration is shown in Fig. 5. The model parameters used are: Young’s modulus E ¼ 30 GPa, Poisson ratio t ¼ 0:25, tensile strength f t ¼ 500 kPa, and cohesive fracture energy Gf ¼ 0:1 KN/m. The fracture is assumed to have plane-strain behaviour with fracture height of 100 m.

dynamic viscosity of gas is assumed as its compressibility is calculated from

cnw ¼

1 patm þ pnw

lnw ¼ 1  108 kPa.s and ð83Þ

where patm ¼ 100 kPa is the reference atmospheric pressure and pnw is the gas pressure. The saturation-suction relationship as well as relative permeabilities of fracturing fluid and host gas for the porous medium is calculated from the van Genuchten–Mualem (VGM) model as

 Sw ¼ Srw þ ð1  Srw Þ 1 þ

s

sref m 2

  krw ¼ Se1=2 1  1  S1=m e  2m krnw ¼ ð1  Se Þ1=2 1  S1=m e

1=ð1mÞ !m ð84Þ ð85Þ ð86Þ

in which the residual saturation Srw ¼ 0, the empirical curve-fitting parameter m ¼ 0:4396, the reference suctions sref ¼ 10 and sref ¼ 100 kPa for the fracture and porous media, respectively, and w Srw . The VGM model was used by the effective saturation Se ¼ S1S rw other researchers for similar problems [38,35]. The simulated examples for slow-injection scenario include: (i) the case without leak-off (kp ¼ 0Þ, (ii) the case with leak-off and low medium permeability (kp ¼ 1  1013 m2) and (iii) the case with leak-off and high medium permeability (kp ¼ 1  1012 m2). The medium porosity in all examples is n ¼ 0:01. Geertsma and de Klerk [21] obtained an approximate analytical solution for the case without leak-off and with negligible toughness (KGD model).

89

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

 3 16 lq 1 wwell ¼ A t3 E0  13 1 Dpwell ¼ B lE02 t3 L=2 ¼ C

 0 3 16 Eq

l

2

t3

ð87Þ ð88Þ

5000 KGD Spence and Sharp (1985) kp = 0 kp = 1E-13 kp = 1E-12

4000

Well Pressure (kPa)

Later, Spence and Sharp [43] extended the solution to include fracture toughness. In the analytical solution, the fracture width at the wellbore wwell , net well pressure Dpwell and the fracture half-length L=2 are given as

3000 2000 1000 0

ð89Þ

0

1

2

3

4

5

Time (min)

12

2

kp ¼ 1  10 m , the equilibrium is reached within five minutes, consequently the fracture width and fracture half-length is stabilised and no fracture propagation is observed. This implies that

Fracture Width at the Welbore (m)

0.003 KGD Spence and Sharp (1985) kp = 0

0.002

kp = 1E-13 kp = 1E-12

0.001

0 0

1

2

3

4

5

Time (min) Fig. 7. Fracture width at the wellbore versus elapsed time in slow-injection examples.

30

Fracture Half-Length (m)

the most permeable case with kp ¼ 1  1012 m2 which represents the most suitable case for this solution. For the case without leak-off (kp ¼ 0Þ very good agreement is found between present model results and KGD analytical solutions for net well pressure and fracture width at the wellbore. KGD solution is for viscosity dominated regime (zero toughness), and the parameters used in the present simulation also falls into the viscosity dominated regime as the dimensionless viscosity and toughness [14] are M ¼ 308 and K ¼ 0:239, respectively. Viscosity dominated regime is in accordance with hydraulic fracturing treatment, where the energy dissipation associated with the flow of viscous fracturing fluid often dominates over the energy dissipation associated with the rock damage [30]. Spence and Sharp [43] calculated higher net well pressure due to including finite toughness. The fracture half-length (distance from fracture origin to the fracture tip) is larger than the analytical solution due to the pressure lag ahead of the fracturing fluid front as shown in Fig. 9. The size of the lag region can be affected by the in situ stresses [19,20]. The in situ stresses are considered zero in these simulations to allow the formation of lag zone to include two phases within the fracture. However, as shown in the next section, if the lag zone is negligible, a single fluid can be considered for the fracture. Negative pressure (suction) is induced in the lag region. If the position of the fracturing fluid front is considered instead of the fracture tip, as shown in Fig. 8, a perfect match is observed between present model and KGD results. Again, Spence and Sharp [43] calculated slightly lower fracture length, while Carbonell et al. [8] calculated higher fracture length for the limiting case of viscosity-storage regime. When the leak-off is allowed by considering permeable host medium, a significant portion of the injected fracturing fluid is dissipated into the surrounding media. This causes smaller fracture width and lower fracture propagation. The net well pressure on the other hand remains higher and maintains a steady leak-off flow towards host medium which becomes equal to the injected flow rate as time elapses. The higher the medium permeability, the sooner the equilibrium. For the case of medium permeability of

Fig. 6. Well pressure versus elapsed time in slow-injection examples.

20

KGD

Spence and Sharp (1985)

Carbonell et al. (1999)

Adachi and Detournay(2008)

kp = 0

kp = 0 (Fluid front)

kp = 1E-13

kp = 1E-12

10

0 0

1

2

3

4

5

Time (min) Fig. 8. Fracture half-length versus elapsed time in slow-injection examples.

0.001 Physical Fracture

0.0008

Fracture Width (m)

where q ¼ Q =h is the fracturing fluid injection rate per unit height of the fracture, E0 ¼ 1Em2 is the plane-strain elastic modulus, l is the dynamic viscosity of fracturing fluid and t is the injection time. The coefficients A ¼ 2:36; B ¼ 1:09 and C ¼ 0:539 for Geertsma and de Klerk [21] model and A ¼ 2:70; B ¼ 1:57 and C ¼ 0:515 for Spence and Sharp [43] model. The results from present study for the net well pressure, the fracture width and the fracture half-length are shown in Figs. 6–8, respectively. Included in these figures are the analytical solutions by Geertsma and de Klerk [21] and Spence and Sharp [43]. In Fig. 8 the solution by Carbonell et al. [8] for the limiting case of storage-viscosity regime and the solution by Adachi and Detournay [2] for leak-off dominated regime are also included. The leak-off dominated solution is calculated for

Fluid Front

0.0006

Fluid lag region

0.0004 0.0002

Fluid saturated region

0 -0.0002 -0.0004

Cohesive tractions zone

-0.0006 -0.0008 -0.001 0

5

10

15

20

25

Fracture Half-Length (m) Fig. 9. Physical fracture shape and fluid front in hydraulic fracturing example with no leak-off (t = 5 min).

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

Vertical Dimension (m)

0.9 506

0.8

504

0.7 0.6

502

0.5

500

0.4 498 0.3 496 0.2 494

0.1 0

1

2

3

4

5

Fracture Width at the Welbore (m)

90

0.005 0.004 0.003

kp = 1E-15

kp = 1E-13

kp = 5E-13

kp = 8E-13

kp = 1E-12

0.002 0.001 0

0

20

40

60

80

100

120

Time (sec) Fig. 12. Fracture width at the wellbore versus elapsed time in fast-injection examples.

6

Horizontal Dimension (m)

Vertical Dimension (m)

0.9 506

0.8

504

0.7 0.6

502

0.5 500 0.4 498 0.3 496

0.2

494

0.1

0

1

2

3

4

5

6

Horizontal Dimension (m) Fig. 10. The invaded zone showed by contours of degree of saturation of hydraulic fracturing fluid within host medium, top: kp ¼ 1  1012 m2, bottom: kp ¼ 1  1013 m2.

Well Pressure (kPa)

5000 4000 3000

kp = 1E-15

kp = 1E-13

kp = 5E-13

kp = 8E-13

kp = 1E-12

2000 1000 0

0

20

40

60

80

100

120

Time (sec) Fig. 11. Well pressure versus elapsed time in fast-injection examples.

all the injected fracturing fluid will be leaked into the host medium. The model results for the case with kp ¼ 1  1012 m2 is in good agreement with the analytical solution by Adachi and Detournay [2] for leak-off dominated regime in early time response; however, in late time response the two models diverge as the current model approaches an equilibrium in which there is no further fracture growth and the entire injected fracturing fluid is dissipated into the host medium, while the analytical solution predicts fracture growth proportional to square root of time. In addition, the analytical solution is for single-phase fluid flow, whereas the numerical solution is for two phase flow. As is shown

later in this paper, single-phase models underestimate the leak-off in partially saturated media, and a two-phase fluid flow model is required to correctly evaluate the leak-off. The shape of the invaded zone, i.e. the zone within host medium saturated by fracturing fluid, also depends on the permeability of the host medium. In Fig. 10, the contours of the degree of saturation of fracturing fluid within host medium are shown for the two examples with permeable host media. For the case with higher permeability, lower pressure gradient is required to push fracturing fluid into the porous medium so the invaded zone shape is short and dispersed, whereas for the case with the lower permeability it is elongated and thin. In the latter case, higher pressure gradient is required to push fracturing fluid into the porous medium, so the preferred path for the fluid is along the fracture. This causes longer fracture propagation and lower injection pressure. The injection rate in the previous examples is relatively low (Q w ¼ 1 m3/min), which produces higher leak-off. In the next set of examples, the injection rate is increased by a factor of 5 (fast-injection) to investigate the sensitivity of the results including the leak-off and its effect on the fracture length and width. In these examples the following model parameters are used: Young’s modulus E ¼ 50 GPa, Poisson ratio t ¼ 0:25, tensile and cohesive fracture energy strength f t ¼ 1000 kPa, Gf ¼ 0:1 KN/m. The lag zone is assumed to be negligible due to in situ stresses, so the fracture is fully saturated by the fracturing fluid and a single phase fluid is assumed within the fracture. This assumption increases the efficiency of the model by reducing the pressure oscillation in the fracture. The number of the iterations required for convergence is reduced to an average number of three iterations per time step compared to an average number of five iterations in the previous examples with two-phase fracture flow. Five cases with porous medium permeability of kp ¼ 1015 ; 1013 ; 5  1013 ; 8  1013 and 1012 m2 are considered and the results for injection pressure, fracture width and fracture half-length are shown in Figs. 11–13. Injection pressure is generally increased compared with the previous slow-injection examples due to the increase in the medium stiffness. However, the increase in the leak-off due to increase in the porous medium permeability affects the pressure response differently in early and late time. As shown in Fig. 11, the increase of leak-off decreases the pressure faster in the early time, however, it approaches to higher value in the late time which is consistent with the previous slow-injection examples, where the higher permeability developed higher injection pressure. The fracture width again is clearly affected by the leak-off where the higher leak-off created lower fracture width as shown in Fig. 12. Similar trend is observed for fracture half-length, where the leak-off reduces the induced fracture length as shown in Fig. 13. However, for both

91

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

50

40

kp = 1E-15

kp = 1E-13

kp = 5E-13

kp = 8E-13

kp = 1E-12

Adachi and Detournay (2008)

Fracture Half-Length (m)

Fracture Half-Length (m)

50

30 20 10 0

0

20

40

60

80

100

kp = 1E-13

40

kp = 1E-12 Adachi and Detournay (2008)

30 20 10 0

120

kp = 5E-13

0

20

40

Time (sec) Fig. 13. Fracture half-length versus elapsed time in fast-injection examples.

Fracture Opening (m)

Well Pressure (kPa)

3000 kp = 1E-13

2000

kp = 5E-13 kp = 1E-12

1000

120

kp = 1E-13

0.001 kp = 5E-13 kp = 1E-12

0 -0.001 -0.002 0

10

20

30

40

50

Fracture Half-Length (m) 0

20

40

60

80

100

120 Fig. 17. Fracture shape at the end of simulation in single-phase examples.

Time (s) Fig. 14. Well pressure versus elapsed time in single-phase examples.

Fracture Width at the Welbore (m)

100

0.002

4000

presented in Fig. 16 is the analytical solution by Adachi and

0.005 kp = 1E-13

0.004 kp = 5E-13

0.003

kp = 1E-12

0.002 0.001 0

80

Fig. 16. Fracture half-length versus elapsed time in single-phase examples.

5000

0

60

Time (min)

0

20

40

60

80

100

120

Time (min) Fig. 15. Fracture width at the wellbore versus elapsed time in single-phase examples.

fracture width and length, the deviation due to the leak-off starts after a certain amount of time elapsed (lag time) which the lag time depends on the permeability of the porous medium. The lower the permeability, the higher the lag time is. Included in Fig. 13 is the analytical solution by Adachi and Detournay [2] for 12

2

the case with kp ¼ 10 m . Again, the analytical solution which is for single-phase fluid underestimates the leak-off and predicts higher fracture length. Three more cases are analysed in which the medium as well as the fracture is assumed to be fully saturated with a single fluid. The model parameters and injection rate are as per the fast-injection examples. The permeability of the medium is assumed as kp ¼ 1013 ; 5  1013 and 1012 m2. Constant pressure boundary condition is assumed for porous medium pressure on the far boundaries. The results for injection pressure, fracture width and half-length are presented in Figs. 14–16, respectively. Also

Detournay [2] for the most permeable case (kp ¼ 1012 m2). In general, good agreement is found between the analytical solution and present model which further confirms the validity of the present model with the leak-off. The single-phase model generally predicts higher well pressure compared with the results of two-phase model as shown in Fig. 14. Higher fracture width and length are also predicted by the single-phase model. Fracture width varies with change in the medium permeability in a similar trend to the one predicted by the two-phase model in which higher permeability (higher leak-off) results in the lower fracture width (Fig. 15). However, the observed trend for the well pressure and fracture length is different from that predicted by the two-phase model. As shown in Figs. 14 and 16, higher leak-off is associated by lower pressure and a little bit higher fracture length. It may seem a bit unusual to have lower induced fracture length when the leak-off is lower, so a close look at the fracture shape may explain why. The fracture shape at the end of simulation is shown in Fig. 17. Although, the lower leak-off cases have lower fracture length, they have higher volume which is consistent with the amount of wasted fracturing fluid due to the leak-off. The ratio of the fracture volume over the injected fluid volume for the three cases with kp ¼ 1013 ; 5  1013 and 1012 m2 are 0.9403, 0.8455 and 0.7984, respectively. Furthermore, the single-phase model predicts lower sensitivity of the results to change in medium permeability which is in contrast with the two-phase model results. This is due to the fact that single-phase model underestimates the leak-off. This demonstrates the importance of the two-phase model with leak-off. 5. Conclusions A three-phase hydro-mechanical model for hydraulic fracturing is proposed. Three phases include: porous solid, fracturing fluid and host fluid. The fracture discontinuity is captured using

92

S. Salimzadeh, N. Khalili / Computers and Geotechnics 69 (2015) 82–92

extended finite element method (XFEM) with cohesive fracture propagation. Two fluids are considered including fracturing fluid and host fluid. Flow through fracture is defined as one-dimensional laminar flow, and flow through porous medium (host reservoir) is defined as two-dimensional Darcy flow. Coupling amongst the phases within the model is achieved using: (i) effective stress concept for linking deformation and host flow models, (ii) leak-off mass transfer for coupling of fracture flow and host flow models, and (iii) hydraulic loading and the fracture volumetric compatibility for coupling of fracture flow and deformation models. The proposed model is discretised spatially using standard Galerkin method and temporally using finite difference technique. The discretised model is implemented into well-known software and verified against several examples available in the literature. The features of the model in capturing leak-off flow are demonstrated through modelling several hydraulic fracturing examples. The impact of the host medium permeability and associated leak-off flow on the net well pressure, fracture propagation and fracture width is studied and results are presented. Results for both slow and fast injection cases are presented. It is shown that the model is able to reduce to fully saturated fracture as well as fully saturated porous medium. It is demonstrated that models based on single-phase flow underestimate the leak-off. References [1] Adachi J, Siebrits E, Peirce A, Desroches J. Computer simulation of hydraulic fractures. Int J Rock Mech Min Sci 2007;44:739–57. [2] Adachi J, Detournay E. Plane strain propagation of a hydraulic fracture in a permeable rock. Eng Fract Mech 2008;75:4666–94. [3] Aliabadi MH, Brebbia CA. Advanced formulations in boundary element methods. Elsevier Applied Science; 1993. [4] Asferg JL, Poulsen PN, Nielsen LO. A consistent partly cracked XFEM element for cohesive crack growth. Int J Numer Meth Eng 2007;72(4):464–85. [5] Batchelor GK. An introduction to fluid dynamics. Cambridge, UK: Cambridge University Press; 1967. [6] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601–20. [7] Budyn E, Zi G, Moës N, Belytschko T. A method for multiple crack growth in brittle materials without remeshing. Int J Numer Meth Eng 2004;61:1741–70. [8] Carbonell R, Desroches J, Detournay E. A comparison between a semianalytical and a numerical solution of a two-dimensional hydraulic fracture. Int J Solids Struct 1999;36(31–32):4869–88. [9] Carter E. Optimum fluid characteristics for fracture extension. Drill Prod Pract 1957:261–70. [10] Carter BJ, Wawrzynek PA, Ingraffea AR. Automated 3-d crack growth simulation. Int J Numer Meth Eng 2000;47(1–3):229–53. [11] Curtin WA, Scher H. Mechanics modeling using a spring network. J Mater Res 1990;5(3):554–62. [12] Dahi-Taleghani A. Analysis of hydraulic fracture propagation in fractured reservoirs: an improved model for the interaction between induced and natural fractures, PhD dissertation, Austin,Texas: The University of Texas at Austin; 2009. [13] de Borst R, Rethore J, Abellan M. A numerical approach for arbitrary cracks in a fluid-saturated medium. Arch Appl Mech 2006;75:595–606. [14] Detournay E. Propagation regimes of fluid-driven fractures in impermeable rocks. Int J Geomech 2004;4:1–11. [15] Dias-da-Costa D, Alfaiate J, Sluys LJ, Júlio E. A comparative study on the modelling of discontinuous fracture by means of enriched nodal and element techniques and interface elements. Int J Fract 2010;161:97–119. [16] Duflot M. The extended finite element method in thermoelastic fracture mechanics. Int J Numer Methods Eng 2008;74:827–47. [17] Economides MJ, Nolte KG, Ahmed U. Reservoir stimulation. Englewood Cliffs, New Jersey: Prentice Hall; 1989.

[18] Fu P, Johnson SM, Carrigan CR. An explicitly coupled hydro-geomechanical model for simulating hydraulic fracturing in arbitrary discrete fracture networks. Int J Numer Anal Meth Geomech 2013;37:2278–300. [19] Garagash DI, Detournay E. Near tip processes of a fluid-driven fracture. ASME J Appl Mech 2000;67:183–92. [20] Garagash D, Detournay E, Adachi J. Multiscale tip asymptotics in hydraulic fracture with leak-off. J Fluid Mech 2011;669:260–97. [21] Geertsma J, de Klerk F. A rapid method of predicting width and extent of hydraulically induced fractures. J Petrol Technol 1969;21:1571–81. [22] Gordeliy E, Peirce A. Coupling schemes for modeling hydraulic fracture propagation using the XFEM. Comput Meth Appl Mech Eng. 2013;253:305–22. [23] Gordeyev YN, Entov VM. The pressure distribution around a growing crack. J Appl Math Mech 1997;61(6):1025–9. [24] Hagoort J, Weatherill BD, Settari A. Modeling the propagation of waterfloodinduced hydraulic fractures. Soc Pet Eng J 1980:293–301. [25] Hughes TJR. The finite element method. Linear static and dynamic finite element analysis. Englewood Cliffs, NJ: Prentice-Hall Inc; 1987. [26] Khalili N, Geiser F, Blight GE. Effective stress in unsaturated soils: review with new evidence. Int J Geomech 2004;4(2):115–26. [27] Khalili N. Two phase fluid flow through fractured porous media with deformable matrix. Water Resour Res 2008;44. Art. No. W00C04. [28] Khoei AR, Vahab M, Haghighat E, Moallemi S. A mesh-independent finite element formulation for modelling crack growth in saturated porous media based on an enriched-FEM technique. Int J Fract 2014;188:79–108. [29] Khoshghalb A, Khalili N, Selvadurai APS. A three-point discretisation technique for parabolic partial differential equations. Int J Numer Anal Meth Geomech 2011;35:406–18. [30] Kovalyshen Y. Fluid-driven fracture in poroelastic medium, Ph.D. thesis; University of Minnesota; 2010. [31] Lecampion B. An extended finite element method for hydraulic fracture problems. Commun Numer Meth Eng 2009;25:121–33. [32] Mathias S, van Reeuwijk M. Hydraulic fracture propagation with 3-D leak-off. Transp Porous Med 2009;80:499–518. http://dx.doi.org/10.1007/s11242-0099375. [33] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46:131–50. [34] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69:813–33. [35] Mohammadnejad T, Khoei AR. Hydro-mechanical modelling of cohesive crack propagation in multiphase porous media using the extended finite element method. Int J Numer Anal Meth Geomech 2013;37:1247–79. [36] Réthoré J, Gravouil A, Combescure A. A combined space time extended finite element method. Int J Numer Meth Eng 2005;64:260–84. [37] Rethore J, de Borst R, Abellan MA. A two-scale approach for fluid flow in fractured porous media. Int J Numer Meth Eng 2007;71:780–800. [38] Rethore J, de Borst R, Abellan MA. A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks. Comput Mech 2008;42:227–38. [39] Salimzadeh S, Khalili N. Coupling reservoir simulation in naturally fractured reservoir: implicit versus explicit formulation. In: IACMAG 2011, Melbourne, Australia; 2011. [40] Salimzadeh S, Khalili N. Consolidation of unsaturated lumpy clays. J Geo-Eng Sci 2014;2:67–82. http://dx.doi.org/10.3233/JGS-141317. [41] Salimzadeh S, Khalili N. A Three-dimensional numerical model for double porosity media with two miscible fluids including geomechanical response. Int J Geomech 2015. http://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000494. [42] Shao Q, Bouhala L, Younes A, Nunez P, Makardi A, Belouettar S. An XFEM model for cracked porous media: effects of fluid flow and heat transfer. Int J Fract 2014;185:155–69. [43] Spence DA, Sharp P. Self-similar solutions for elastohydrodynamic cavity flow. Proc R Soc Lond A 1985;400:289–313. [44] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modeling. Int J Numer Meth Eng 2000;48:1549–70. [45] Weber N, Siebert P, Willbrand K, Feinendegen M, Clauser C, Fries TP. The XFEM with an explicit-implicit crack description for hydraulic fracture problems. Eff Sustain Hydraul Fract 2013. http://dx.doi.org/10.5772/56383. Chapter 35. [46] Wells GN, Sluys LJ. A new method for modelling of cohesive cracks using finite elements. Int J Numer Meth Eng 2001;50(12):2667–82. [47] Witherspoon PA, Wang JSY, Iwai K, Gale JE. Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res 1980;16:1016–24.